powerflow.txt

潮流计算

%==========================================================================
%==========================================================================
%==========================================================================
%潮流计算MATLAB 粗略程序                 [e-mail]   2009.3.27
%==========================================================================
%==========================================================================
%==========================================================================
%creat a new_data
t=0;
s=0;
r=0;
w=0;
number=input('How many node are there=');
% Convert Pq to a new array
for ii=1:number
    if data(ii,4)==1
        t=t+1;
            for jj=1:16
                new_data1(t,jj)=data(ii,jj);
            end;
            a(1,t)=ii;
            s=s+1;                                        %record the number of the PQ node
    end;
end;

%Convert pv to a new array
for ii=1:number
    if data(ii,4)==2
        t=t+1;
            for jj=1:14
                new_data1(t,jj)=data(ii,jj);
            end;
            a(1,t)=ii;
            r=r+1;                                        %record the number of the PV node
    end;
end;
%Convert set_v to a new array
for ii=1:number
    if data(ii,4)==3
        t=t+1;
            for jj=1:14
                new_data1(t,jj)=data(ii,jj);
            end;
            a(1,t)=ii;
            w=w+1;
    end;
end;


%creat a new_data2
[x,y]=size(data2)              
for ii=1:x
    for jj=1:2
        for mm=1:number
            if data2(ii,jj)==a(1,mm)
                new_data2(ii,jj)=mm;
            end;
        end;
    end;
end;

for ii=1:x
    for jj=3:14
        new_data2(ii,jj)=data2(ii,jj);
    end;
end;

%creat a Y 
Y=zeros(number,number);
YY=zeros(number,number);
yy=zeros(number,number);

for ii=1:x
   % for jj=1:14
        iii=new_data2(ii,1);
        jjj=new_data2(ii,2);
        if new_data2(ii,5)==2
            sub=new_data2(ii,6)./(new_data2(ii,7).*new_data2(ii,7)+new_data2(ii,6).*new_data2(ii,6))-new_data2(ii,7)./(new_data2(ii,7).*new_data2(ii,7)+new_data2(ii,6).*new_data2(ii,6))*i;
            Y(iii,jjj)=-sub./new_data2(ii,14);
            YY(iii,jjj)=sub./new_data2(ii,14);
            Y(jjj,iii)=-sub/new_data2(ii,14);
             YY(jjj,iii)=sub./new_data2(ii,14);
            yy(iii,jjj)=(1.-new_data2(ii,14))./(new_data2(ii,14).*new_data2(ii,14)).*sub;
            yy(jjj,iii)=(new_data2(ii,14)-1)./(new_data2(ii,14)).*sub;
        else
            Y(iii,jjj)=-new_data2(ii,6)./(new_data2(ii,7).*new_data2(ii,7)+new_data2(ii,6).*new_data2(ii,6))+new_data2(ii,7)./(new_data2(ii,7).*new_data2(ii,7)+new_data2(ii,6).*new_data2(ii,6))*i;
             YY(iii,jjj)=new_data2(ii,6)./(new_data2(ii,7).*new_data2(ii,7)+new_data2(ii,6).*new_data2(ii,6))-new_data2(ii,7)./(new_data2(ii,7).*new_data2(ii,7)+new_data2(ii,6).*new_data2(ii,6))*i;
            Y(jjj,iii)=-new_data2(ii,6)./(new_data2(ii,7).*new_data2(ii,7)+new_data2(ii,6).*new_data2(ii,6))+new_data2(ii,7)./(new_data2(ii,7).*new_data2(ii,7)+new_data2(ii,6).*new_data2(ii,6))*i;
            YY(jjj,iii)=new_data2(ii,6)./(new_data2(ii,7).*new_data2(ii,7)+new_data2(ii,6).*new_data2(ii,6))-new_data2(ii,7)./(new_data2(ii,7).*new_data2(ii,7)+new_data2(ii,6).*new_data2(ii,6))*i;
            yy(iii,jjj)=new_data2(ii,8)./2.*i;
            yy(jjj,iii)=new_data2(ii,8)./2.*i;
        end;
    %end;
end;

for iii=1:number
    Y(iii,iii)=0;
end;

%for ii=1:x
  %  for jj=1:14
    for     iii=1:number
        for jj=1:number
   % if iii~=jj
    Y(iii,iii)=Y(iii,iii)+YY(iii,jj)+yy(iii,jj);
      % end;
   end;
   Y(iii,iii)=Y(iii,iii)+new_data1(iii,16).*i;
end;


%creat B, G

for ii=1:number
        for jj=1:number
            G(ii,jj)= real(Y(ii,jj));      
        B(ii,jj)= imag(Y(ii,jj));
    end;
end;


%creat Initial_P Initial_Q Initial_V
for ii=1:(s+r)
    set_P(ii,1)=(new_data1(ii,9)-new_data1(ii,7))./100;
end;
for ii=1:s;
    set_Q(ii,1)=(new_data1(ii,10)-new_data1(ii,8))./100;
end;

for ii=1:r
    set_V(ii,1)=new_data1(ii+s,12).*new_data1(ii+s,12);%try to modify for sike of correcting
end;

Initial_p_q_v=[set_P;set_Q;set_V];
disp(Initial_p_q_v);
%creat Initial_e,Initial_f
for ii=1:number-1
    e(ii,1)=1;
    f(ii,1)=0.0;%change f to test used to be 1.0
end;

  e(number,1)=new_data1(number,12);
  f(number,1)=0;
 % e(64,1)=0.88;%test 118ieee 
 % f(64,1)=0.39395826829394;
   % f(14,1)=0;
   % e(10,1)=1.045;
   %e(11,1)=1.01;
    %e(12,1)=1.07;
    %e(13,1)=1.09;
%//////////////////////////////////////////////////////////////////////////
%/////////////////////////////////////////////////////////////////////////
%//////////////////////////////////////////////////////////////////////////
%//////////////////////////////////////////////////////////////////////////
% Start NEWTOWN CALULATION
for try_time=1:25
    
%Creat every node consume P Q and U
n=s;
m=r;
for ii=1:(n+m)
    sum1=0;
    for jj=1:(n+m+1)
        sum1=sum1+e(ii,1).*(G(ii,jj).*e(jj,1)-B(ii,jj).*f(jj,1))+f(ii,1).*(G(ii,jj).*f(jj,1)+B(ii,jj).*e(jj,1));
    end;
    p(ii,1)=sum1;
end;

for ii=1:n
    sum2=0;
    for jj=1:(n+m+1)
        sum2=sum2+f(ii,1).*(G(ii,jj).*e(jj,1)-B(ii,jj).*f(jj,1))-e(ii,1).*(G(ii,jj).*f(jj,1)+B(ii,jj).*e(jj,1));
    end;
    q(ii,1)=sum2;
end;

disp('q=');
disp(q);

u=zeros((n+m),1);
for ii=(n+1):(n+m)
    u(ii,1)=e(ii,1).*e(ii,1)+f(ii,1).*f(ii,1);
end;

for ii=n+1:(n+m)  
    extra_u((ii-n),1)=u(ii,1);
end;
disp('extra_u=');
disp(extra_u);

sum=[p;q;extra_u];
disp(sum)
disp(s);
disp(p);
%creat Jacobian 
disp(n);
disp(m);
for ii=1:(n+m)
    for jj=1:(n+m)
        if (ii~=jj)
            PF(ii,jj)=B(ii,jj).*e(ii,1)-G(ii,jj).*f(ii,1);
           
            PE(ii,jj)=-G(ii,jj).*e(ii,1)-B(ii,jj).*f(ii,1);
           
            
        else
            ss=0;
            qq=0;
            
            for num=1:(n+m+1)
                ss=ss+G(ii,num).*f(num,1)+B(ii,num).*e(num,1);
                qq=qq+G(ii,num).*e(num,1)-B(ii,num).*f(num,1);
            end;
            
            PF(ii,jj)=-ss+B(ii,jj).*e(ii,1)-G(ii,jj).*f(ii,1);%TEST+1
            PE(ii,jj)=-qq-G(ii,jj).*e(ii,1)-B(ii,jj).*f(ii,1);%TEST+1
            
        end;
    end;
end;


for ii=1:n
    for jj=1:m+n
        if (ii~=jj)
           
           QE(ii,jj)=B(ii,jj).*e(ii,1)-G(ii,jj).*f(ii,1);%TEST+1
            
            QF(ii,jj)=G(ii,jj).*e(ii,1)+B(ii,jj).*f(ii,1);%TEST+1
            
        else
            ss=0;
            qq=0;
            
           for num=1:(n+m+1)
                ss=ss+G(ii,num).*f(num,1)+B(ii,num).*e(num,1);
                qq=qq+G(ii,num).*e(num,1)-B(ii,num).*f(num,1);
            end;
            
            
            QF(ii,jj)=-qq+G(ii,jj).*e(ii,1)+B(ii,jj).*f(ii,1);%TEST+1
            QE(ii,jj)=ss+B(ii,jj).*e(ii,1)-G(ii,jj).*f(ii,1);%TEST+1
           
        end;
    end;
end;
%disp('QF');
%disp(QF);
%disp('QE');
%disp(QE);
UE=zeros((n+m),(n+m));
UF=zeros((n+m),(n+m));

for ii=n+1:n+m
    for jj=1:(n+m)
        if (ii~=jj)
           
            UE(ii,jj)=0;
            UF(ii,jj)=0;
        else
            ss=0;
            qq=0;
            
            for num=1:(n+m+1)
                ss=ss+G(ii,num).*f(num,1)+B(ii,num).*e(num,1);
                qq=qq+G(ii,num).*e(num,1)-B(ii,num).*f(num,1);
            end;
            
            
            UF(ii,jj)=-2.*f(ii,1);
            UE(ii,jj)=-2.*e(ii,1);
        end;
    end;
end;

for ii=(n+1):(n+m)
    for jj=1:(n+m)
    extra_UE((ii-n),jj)=UE(ii,jj);
    extra_UF((ii-n),jj)=UF(ii,jj);
    end;
end;
%disp('extra_UE');
%disp(extra_UE);
%disp('extra_Uf');
%disp(extra_UF);


Jacobian=[PF,PE;QF,QE;extra_UF,extra_UE];


%disp('Jacobian=');
%disp(Jacobian);

%creat substract result
substract_result=Initial_p_q_v-sum;
%disp('substract_result');
%disp(substract_result);
%calculate delta_f_e
delta_f_e=-inv(Jacobian)*substract_result;
%disp(delta_f_e);
for ii=1:number-1;
    f(ii,1)=f(ii,1)+delta_f_e(ii,1);
    e(ii,1)=e(ii,1)+delta_f_e(ii+number-1,1);
end;

if max(substract_result)<1e-4
    break;
end ;
end;

%disp('substract_result');
%disp(substract_result);
%disp('e=');

%disp(e);
%disp('f=');
%disp(f);

for ii=1:number
   uuu(ii,1)= e(ii,1).*e(ii,1)+f(ii,1).*f(ii,1);
U_RESULT(ii,1)=sqrt(uuu(ii,1));
end;

for ii=1:number
  for  jj=1:number
    if ii==a(1,jj)
        Old_Uresult(ii,1)=U_RESULT(jj,1)
    end;
end;
end;

for ii=1:number
   
 
        Old_Uresult(ii,2)=ii;
end;
%disp('U_result');

%disp(U_RESULT);
disp('=====================================');
disp('The last result is :')

disp('===========U===================BUS-NO.');
disp('U=')
disp(Old_Uresult);


%calculate the angle
PI=3.141592
for ii=1:number
    
    Angle(ii,1)=atan(f(ii,1)./e(ii,1))./PI*180;
    
end;

for ii=1:number
  for  jj=1:number
    if ii==a(1,jj)
        Old_Angle(ii,1)=Angle(jj,1);
        Old_Angle(ii,2)=ii;
    end;
end;
end;
disp('=======Angle===================BUS-NO.');
disp('Angle=');
disp(Old_Angle);
disp('=====Try-times=======================')
disp('Try-times=')
disp(try_time);

%disp('p====================');
%disp(p);

 % for  jj=1:number
 % if a(1,jj)==118
 %       man=jj
  %  end;

%end;
%disp('man=========');
%disp(man)

sum4=0;
for jj=1:number
    Y_conj(number,jj)=conj(Y(number,jj));
    sum4=sum4+Y_conj(number,jj).*(e(jj,1)-f(jj,1)*i);
end;
sum4=sum4*e(number,1);
disp('=======P Q Of The Swing node======BUS-NO');
%disp(sum4);

Blance_Q(1,1)=imag(sum4)*100;
Blance_Q(1,2)=a(1,number);
Blance_P(1,1)=real(sum4)*100;
Blance_P(1,2)=a(1,number);
disp('Q Of the Swing node= ');

disp(Blance_Q);
disp('P  Of the Swing node= ')
disp(Blance_P);
disp('=================================BUS-NO');
%calculate the Q of the P-V node

Q_PV_node=zeros(number,2);
Y_conj=conj(Y);
for ii=(s+1):(s+r)
    for jj=1:number
        
    Q_PV_node(ii,1)=Q_PV_node(ii,1)+(e(ii,1)+f(ii,1)*i).*(Y_conj(ii,jj).*(e(jj,1)-f(jj,1)*i));
    end;
end;
for ii=(s+1):(s+r);
Q_PV_node(ii,1)=Q_PV_node(ii,1).*100+new_data1(ii,8)*i;
end;

for ii=1:number
    old_number=a(1,ii);
    Q_PV_node_old(old_number,1)=Q_PV_node(ii,1);
end;
for ii=1:number
Q_PV_node_old(ii,1)=imag(Q_PV_node_old(ii,1));
end;
for ii=1:number
    Q_PV_node_old(ii,2)=ii;
end;

disp('Q gen=');

disp(Q_PV_node_old);