load flow decoupled

by

 

this program will help to find the load flow analysis of decoupled

decop.m
                                                                     
                                                                     
                                                                     
                                             
%%%%%%%%%%%%%%%% ==== DECOUPLED ESWAR's ==== %%%%%%%%%%%%%
clc
A=[ 1       2.0         1.0         0           0           1.04
    2       0.0         0.0         0.5         1.0         1   
    3       1.5         0.6         0.0         0           1.04];

seimp=0.02+0.08i; %% series impedence
shimp=1/seimp;    %% series admittance
shadm=0.02i;      %% shunt admittance
n=length(A(:,1)); %% length of matrix of given data
n1=input('enter the number of pv bus');
n2=n-n1-1;
Y=zeros(n,n);       %% initialisation of Y matrix
G=zeros(n,n);       %% initialisation of G matrix
B=zeros(n,n);       %% initialisation of B matrix
cP=zeros(n,1);      %% initialisation of caluculated P matrix
cQ=zeros(n,1);      %% initialisation of Caluculated Q matrix
sP=zeros(n,1);      %% initialisation of Specified P matrix
sQ=zeros(n,1);      %% initialisation of Specified Q matrix
dP=zeros(n,1);      %% initialisation of delta  P matrix
dQ=zeros(n,1);      %% initialisation of delta  Q matrix
th=[0;0;0];         %% initialisation of voltage Angle matrix
J=zeros(n,n);       %% initialisation of Jacobian matrix
J11=zeros(n1+n2,n1+n2);
J12=zeros(n1+n2,n2);
J21=zeros(n2,n1+n2);
J22=zeros(n2,n2);
%=================================calc of Y and Y angle ====
for i=1:n
    for j=1:n
        if(i==j)
             Y(i,j)=2*shimp+shadm;
              B(i,j)=imag(Y(i,j));
              G(i,j)=real(Y(i,j));
         else
            Y(i,j)=-shimp;
            B(i,j)=imag(Y(i,j));
            G(i,j)=real(Y(i,j));
        end
    end
end
disp(' Y matrix ')
disp(Y)
%============================= CP   &   CQ==
for i=1:n     
    for j=1:n
      cP(i,1)=cP(i,1)+A(i,6)*A(j,6)*(G(i,j)*cos(th(i,1)-th(j,1))+B(i,j)*sin(th(i,1)-th(j,1)));
      cQ(i,1)=cQ(i,1)+A(i,6)*A(j,6)*(G(i,j)*sin(th(i,1)-th(j,1))-B(i,j)*cos(th(i,1)-th(j,1)));
    end
end
disp(' Cp,,,,,,CQ')
disp(cP)
disp(cQ)
%==================== specified Q & specified P=====
    for i=2:n
         sP(i,1)=A(i,4)-A(i,2);
         sQ(i,1)=A(i,5)-A(i,3);
    end
    disp(' Sp,,,,,,SQ')
    disp(sP)
    disp(sQ)
%====================== delta P & delta Q  ===Power mismatches
for i=2:n
   dP(i,1)=sP(i,1)-cP(i,1);
   dQ(i,1)=sQ(i,1)-cQ(i,1);
end
disp(' Dp,,,,,,DQ')
disp(dP)
disp(dQ)  
%==================== Jacobian matrix ====
for i=1:n1+n2  
    ii=i+1;
    for j=1:n1+n2
        jj=j+1;
        if(i==j)
        J11(i,i)=-cQ(ii,1)-B(ii,ii)*A(ii,6)^2;
        else
        J11(i,j)=A(ii,6)*A(jj,6)*(G(ii,jj)*sin(th(ii,1)-th(jj,1))-B(ii,jj)*cos(th(ii,1)-th(jj,1)));
        end
    end
end
for i=1:n2  
    ii=i+1;
    for j=1:n2
        jj=i+1;
        if(i==j)
        J22(i,i)=cQ(ii,1)-B(ii,ii)*A(ii,6)^2;
        else
        J22(i,j)=J11(ii,jj);
        end
    end
end
disp(' Total jacobian')
J=[J11 J12;J21 J22];
disp(J)
%%================calculation of  delta( voltage angle) && delta V
ddP=zeros(n-1,1);
dth=zeros(n-1,1);
for i=1:n-1
    ddP(i,1)=dP(i+1,1);
end
ddQ=zeros(n2,1);
dV=zeros(n,1);
for i=1:n2
    ddQ(i,1)=dQ(i+1,1);
 end
M=[ddP;ddQ];
N=inv(J)*M;
for i=1:n1+n2
    dth(i+1,1)=N(i,1);
end
for i=1:n-n1-n2
    dV(i+1,1)=N(n1+n2+i,1);
end
%%=================== voltages && delta after interation
 th=th+dth;
 disp(' Voltage Angles ')
 disp(th);
 A(:,6)=A(:,6)+dV;
 disp(' voltages of interation 1')
 disp(A(:,6));

  


Contact us