newton_cankao.m

自己写的牛顿—拉夫逊潮流计算源码,紧供参考

clear
clc
%本程序的功能是用牛顿-拉夫逊法进行潮流计算
% n=input('请输入节点数:n=');
% nl=input('请输入支路数:nl=');
% ph=input('请输入平衡母线节点号:ph=');
% pre=input('请输入误差精度pre=');
% Z=input('请输入由支路参数形成的矩阵Z=[节点号 节点号 支路阻抗    系数       对地导纳]');
% A=input('请输入各节点参数形成的矩阵:A=[节点注入功率  节点初始电压   节点类型      ]');



%具体例子
n=5;
nl=7;
ph=1;
pre=1e-5;
Z=[1  2  0.02+0.06i  1   0; 
    1 3  0.08+0.24i 1    0;
    2  3 0.06+0.18i 1    0;
    2 4 0.06+0.18i 1    0;
    2 5 0.04+0.12i 1   0;
    3 4 0.01+0.03i 1  0;
    4 5 0.08+0.24i 1  0];
Y=zeros(n);
for i=1:nl
     p=Z(i,1);q=Z(i,2);
    Y(p,q)=Y(p,q)-1./(Z(i,3)*Z(i,4));
    Y(q,p)=Y(p,q);
    Y(q,q)=Y(q,q)+1./(Z(i,3)*Z(i,4)^2);
    Y(p,p)=Y(p,p)+1./Z(i,3);
end %形成节点导纳矩阵   
disp('节点导纳矩阵Y为:');
disp(Y);
G=real(Y);B=imag(Y);

A=[0           1.06      1  ;   %1号平衡节点,
    0.2+0.2i      1      2   ;
    -0.45-0.15i   1      2  ;
    -0.4-0.05i    1      2  ;
    -0.6-0.1i    1       2   
   ];   
for i=1:n
    e(i)=real(A(i,2));
    f(i)=imag(A(i,2));
   s(i)=A(i,1);
end
P=real(s);Q=imag(s);
count=0;flag=1;NO=2*n;N=NO-1;
%*******************************************开始求雅可比矩阵*********************
while   flag~=0              %  循环标志flag
        flag=0;
        fprintf('\n*************************************第  %d  次  迭  代******************************\n\n',count);
        for i=1:n
             if i~=ph
             C(i)=0;
             D(i)=0;
                for j=1:n
                C(i)=C(i)+G(i,j)*e(j)-B(i,j)*f(j);
                D(i)=D(i)+G(i,j)*f(j)+B(i,j)*e(j);
                end
            P1=C(i)*e(i)+f(i)*D(i);              %disp('P1=');disp(P1);
            Q1=f(i)*C(i)-D(i)*e(i);           %disp('Q1=');disp(Q1);       
            V1=e(i)^2+f(i)^2;
 %******************************若节点不是PV节点,即为PQ节点时 ***************           
            if A(i,3)~=3               
                DP=P(i)-P1;             %求DPi,DQi
                DQ=Q(i)-Q1;       
           
            for j=1:n     
                 if j~=ph & j~=i       %节点不为平衡节点时非对角线元素
                         H1=-B(i,j)*e(i)+G(i,j)*f(i); %disp('H1');
                         N1=G(i,j)*e(i)+B(i,j)*f(i);
                         J1=-N1;
                         L1=H1;                       
                         p=2*i-3;q=2*j-3;H(p,q)=H1;H(p,N)=DP;  m=p+1;                     %行成矩阵块
                         H(m,q)=J1;H(m,N)=DQ;q=q+1;H(p,q)=N1;H(m,q)=L1; 
                        
                  elseif j==i & j~=ph  %节点不为平衡节点时对角线元素
                         HH1=D(i)-B(i,i)*e(i)+G(i,i)*f(i);
                         NN1=C(i)+G(i,i)*e(i)+B(i,i)*f(i);
                         JJ1=C(i)-G(i,i)*e(i)-B(i,i)*f(i);
                         LL1=-D(i)-B(i,i)*e(i)+G(i,i)*f(i);                    
                       
                         p=2*i-3;q=2*j-3;H(p,q)=HH1;H(p,N)=DP;m=p+1;                 %行成矩阵块
                         H(m,q)=JJ1;H(m,N)=DQ;q=q+1;H(p,q)=NN1;H(m,q)=LL1;
 
                   end
               end
   %**********************若节点是PV节点*****************************              
           else                              
                 DP=P(i)-P1;
                 DV=V(i)^2-V1;
                 for j=1:n
                    if j~=ph & j~=i               %节点不为平衡节点时非对角线元素

                         HV1=-B(i,j)*e(i)+G(i,j)*f(i); 
                         NV1=G(i,j)*e(i)+B(i,j)*f(i);
                         R1=0;
                         R2=0;
                         p=2*i-3;q=2*j-3;H(p,q)= HV1;H(p,N)=DP;m=p+1;
                         H(m,q)=R1;H(m,N)=DV;q=q+1;H(p,q)= NV1;H(m,q)=R2;
                      
                     elseif j==i & j~=ph               %节点不为平衡节点时对角线元素
                        HHV1=D(i)-B(i,i)*e(i)+G(i,i)*f(i);
                        NNV1=C(i)+G(i,i)*e(i)+B(i,i)*f(i);
                        RR1=-2*e(i);
                        RR2=-2*f(i);
                        p=2*i-3;q=2*j-3;H(p,q)=HHV1;H(p,N)=DP;m=p+1;
                        H(m,q)=RR1;H(m,N)=DV;q=q+1;H(p,q)= NNV1;H(m,q)=RR2;
                     end
                end
            end
         end
      end   
DJ=H(:,[1:NO-2]);   %雅可比矩阵矩阵DJ
disp('雅克比矩阵DJ:');
disp(DJ);
DW=H([1:NO-2],N); %修正方程常数项矩阵DW
disp('DW=');
disp(DW)
DY=DJ\DW;%修正量DY
 %输出电压修正量
for i=1:2:NO-3
    fprintf('电压修正量Df=  %.6f      De=  %.6f\n',DY(i),DY(i+1));   
  
end
for i=1:NO-2
          DET=abs(DW(i));
          if DET>=pre   
           flag= flag+1; disp('flag=');disp(flag);
           end
end
 count=count+1;
F=zeros(n,1);E=zeros(n,1) ;
 for i=2:n
        f(i)=f(i)+DY(2*i-3);
        e(i)=e(i)+DY(2*i-2);       
end

E=complex(e,f);
  fprintf('节点电压第%d的近似值:',count);
   disp(E);  
 %各节点的功率  
 for p=1:n
    C(p)=0;
    for q=1:n
        C(p)=C(p)+conj(Y(p,q))*conj(E(q));
    end
    S(p)=E(p)*C(p);
  end    
end
%************************while循环结束*********************
disp('各节点的功率S为(节点号从小到大排列):');
disp(S);
disp('平衡节点功率为:');
disp(S(ph));

%以下计算线路功率
disp('线路始端功率Sij');
disp('E='); disp(E);
for i=1:nl
    p=Z(i,1);q=Z(i,2);
    Sij(i)=E(p)*[conj(E(p))*conj(Z(p,5))+[conj(E(p))- conj(E(q))]*conj(Y(p,q))];
    disp(Sij(i));       %计算支路首端功率Si
end

disp('线路末端功率Sji');
for i=1:nl
        p=Z(i,1);q=Z(i,2);
   
       Sji(i)=E(q)*[conj(E(q))*conj(Z(q,5))+[conj(E(q))- conj(E(p))]*conj(Y(q,p))];
disp(Sji(i));
end                      %计算支路末端功率Sj

disp('各条支路的功率损耗为Sij+Sji=:');             
for i=1:nl

   DS(i)=Sij(i)+Sji(i);
   disp(DS(i));
end