Lab Code

clear all
clc
%Initialization
V=[];
V(1,:)=[1,1,1];
w=[];
w(1)=1;
Ploss(1) = 0;
%Grid Parameters
D_w = 100;
D_v = 20;
w0=1; %nominal value
V0 = 1;
Vbase= 381;
Sbase= 10000;
fbase = 50;
Ibase = Sbase/(sqrt(3)*Vbase);
Zbase = (sqrt(3)*Vbase^2)/Sbase;
P_4 = 1;
P_5 = 1;
l_1=[0.01:0.01:0.05] ;
% r_1=[0.01:0.01:0.05];
%Algorithm
for a =1:5
%     r1=r_1(a);
    r1=0.0103 ;
    l1=l_1(a);
for i=2:10000000000
%Y-Matrix
    [Y_BUS Z_12 Z_23] = Y_bus(w(i-1),r1,l1);


 %All Buses are drop Buses, Calculate PGk,QGk for 3 bus system:
 for n =1:3%n stands for bus index; i denotes number of iterartion
        Pg(i,n) = D_w*(w0 - w(i-1)); % P at Generator[1 2 3]
        Qg(i,n) = D_v*(V0 - abs(V(i-1,n))); % Q at Generator[1 2 3] %行数代表循环的index列数代表busindex
 end
 %Bus 1&3 contains load:add the load

    Pb(i,:) = -[P_4 0 P_5] + Pg(i,:); % P injected at Bus[1 2 3]
    Qb(i,:) = Qg(i,:); % Q injected at Bus[1 2 3]

%  %Calculate Vk(i+1)using Pb,Qb, directly come from book fomular

    V(i,1)= (1/Y_BUS(1,1))*((Pb(i,1)-(Qb(i,1)*j))/conj(V(i-1,1))-Y_BUS(1,2)*V(i-1,2)-Y_BUS(1,3)*V(i-1,3));   %Voltage at Bus 1
    V(i,2)= (1/Y_BUS(2,2))*((Pb(i,2)-(Qb(i,2)*j))/conj(V(i-1,2))-Y_BUS(2,1)*V(i,1)-Y_BUS(2,3)*V(i-1,3));   %Voltage at Bus 2
    V(i,3)= (1/Y_BUS(3,3))*((Pb(i,3)-(Qb(i,3)*j))/conj(V(i-1,3))-Y_BUS(3,1)*V(i,1)-Y_BUS(3,2)*V(i,2)); %Voltage at Bus 3
%
%   %Calculate power loss in the system
%
    I12(i) = (V(i,1)-V(i,2))*(-Y_BUS(1,2)); %Current from Bus1 to Bus2
    I23(i) = (V(i,2)-V(i,3))*(-Y_BUS(2,3)); %Current from Bus2 to Bus3
%
    Ploss(i) = ((abs(I12(i))^2)*(real(Z_12)))+((abs(I23(i))^2)*(real(Z_23)));
%
    w(i) = w0-((P_4 + P_5 + Ploss(i))/(3*D_w));
   %convert pu
    Vbus= [abs(V(i,1)*Vbase), abs(V(i,2)*Vbase), abs(V(i,3)*Vbase)];
    Vbusangle = [angle(V(i,1))*(180/pi), angle(V(i,2))*(180/pi), angle(V(i,3))*(180/pi)];
    I = [abs(I12(i))*Ibase, abs(I23(i))*Ibase];
    Iangle = [angle(I12(i))*(180/pi), angle(I23(i))*(180/pi)];
    PG = [Pg(i,1)*Sbase, Pg(i,2)*Sbase, Pg(i,3)*Sbase];
    QG= [Qg(i,1)*Sbase, Qg(i,2)*Sbase, Qg(i,3)*Sbase];
    f= (w(i)*fbase);

    %stopping condition.
      if norm(V(i,:)-V(i-1,:))<1e-9 && norm(w(i)-w(i-1))<1e-9

        break;

      end

end

 VBus(a,:)= Vbus;
 PG_chnageload(a,:)= PG;
 QG_chnageload(a,:)=QG;
 f_chnageload(a,:)=f;
end


figure(1)
plot(l_1,VBus(:,1),'-*')
hold on
plot(l_1,VBus(:,2),'-*')
hold on
plot(l_1,VBus(:,3),'-*')
xlabel('Ll1(p.u.)') % x-axis label
ylabel('Voltage (V)')
title('The voltage at each bus')
legend('Bus1','Bus2','Bus3','location','northwest')

figure(2)
plot(l_1, QG_chnageload(:,1),'-*')
hold on
plot(l_1, QG_chnageload(:,2),'-*')
hold on
plot(l_1, QG_chnageload(:,3),'-*')
xlabel('Ll1(p.u.)')
ylabel('Reactive Power (Var)')
title('The reactive power of each genrator')
legend('G1','G2','G3','location','northwest')
%
figure(3)
plot(l_1, PG_chnageload(:,1),'-*')
hold on
plot(l_1, PG_chnageload(:,2),'-*')
hold on
plot(l_1, PG_chnageload(:,3),'-*')
xlabel('Ll1(p.u.)')
ylabel('Active Power (W)')
title('The active power of each generator')
legend('G1','G2','G3','location','northwest')
%
% %
figure(4)
plot(l_1,  f_chnageload,'-*')

xlabel('Ll1(p.u.)') % x-axis label
ylabel('frequency (Hz)')
title('The frequency of the microgrid')
legend('frequency f')

function [Y_BUS,Z_12,Z_23] = Y_bus(w,r1,l1)

Z_12=r1+1j*w*l1;
Z_23=0.0358+1j*w*0.0262;
Y11=inv(Z_12);
Y12=-inv(Z_12);
Y13=0;
Y21=-inv(Z_12);
Y22=inv(Z_12)+inv(Z_23);
Y23=-inv(Z_23);
Y31=0;
Y32=-inv(Z_23);
Y33=inv(Z_23);
Y_BUS=[Y11 Y12 Y13;Y21 Y22 Y23;Y31 Y32 Y33];

end