程序.txt

牛顿拉夫逊算法进行潮流计算,希望能对大家的学习有所帮助

%牛顿法潮流计算程序
clear,clear all

l(1,3)=(0.0154+0.158i)*1; % 切除一回线路，阻抗乘2
l0(1,3)=0.388i;
l(1,29)=0.0001i;

t(2,9)=0.0217i;
tk(2,9)=1.075;

l(3,31)=0.0001i;

t(4,19)=0.064i;
tk(4,19)=1.025;

t(5,18)=0.0375i;
tk(5,18)=1.075;

t(6,15)=0.0001i;
tk(6,15)=1.00;

t(7,10)=0.001i;
tk(7,10)=1.00;

l(8,18)=0.001i;

l(9,22)=0.0559+0.218i;
l0(9,22)=0.1954i;
l(9,24)=0.0147+0.104i;

l(9,23)=0.0034+0.0131i;

t(10,9)=-0.002i;
tk(10,9)=1;
t(10,11)=0.018i;
tk(10,11)=1.025;


l(11,25)=0.0001i;

l(12,26)=0.0001i;
l(12,27)=0.0001i;

l(13,28)=0.0001i;

l(14,19)=0.0034+0.02i;

t(15,12)=0.018i;
tk(15,12)=1.025;
t(15,14)=-0.002i;
tk(15,14)=1;

l(16,18)=0.0033+0.0333i;
l(16,19)=0.0578+0.218i;
l0(16,19)=0.1887i;
l(16,20)=0.0165+0.0662i;
l0(16,20)=0.2353i;
l(16,21)=0.0374+0.178i;
l0(16,21)=0.164i;
l(16,29)=0.0001i;
l(16,16)=-1i;

t(17,13)=0.01i;
tk(17,13)=1.025;
t(17,16)=0.001i;
tk(17,16)=1.027;
t(17,33)=0.0337i;
tk(17,33)=1;


l(19,21)=0.0114+0.037i;
l(19,30)=(0.0196+0.0854i)*1;   % 切除线路，阻抗乘0
l0(19,30)=0.081i;

l(20,22)=0.0214+0.0859i;
l0(20,22)=0.3008i;

l(21,22)=0.015+0.0607i;
l0(21,22)=0.2198i;

l(22,23)=0.0537+0.19i;
l0(22,23)=0.1653i;

l(23,24)=0.0106+0.074i;

l(25,26)=0.0033+0.0343i;
l0(25,26)=1.8797i;
l(25,25)=0.7318i;

l(26,26)=0.7318i;

l(27,28)=0.00245+0.0255i;
l0(27,28)=1.395i;
l(27,27)=0.7318i;

l(30,31)=0.0001i;

t(32,22)=0.0124i;
tk(32,22)=1.075;

t(33,17)=0.0337i;
tk(33,17)=1;

t(34,30)=0.0438i;
tk(34,30)=1.025;

t(35,31)=0.0328i;
tk(35,31)=1.025;

t(36,24)=0.015i;
tk(36,24)=1.075;

t(36,36)=0;
tk(36,36)=0;
l(36,36)=0;
l0(36,36)=0;
ld(36,36)=0;
yt(36,36)=0;

for i=1: 36
    for j=1: 36
        if t(i,j)~=0
            %t(i,j)=t(i,j)/(tk(i,j)*tk(i,j));
           % tk(i,j)=1;
        end
    end
end
            


for i=1: 36
    for j=1: 36
        if l(i,j)~=0
        l(j,i)=l(i,j);
        ld(i,j)=1/l(i,j);
        ld(j,i)=1/l(i,j);
        l0(j,i)=l0(i,j);
        end
    end
end
%                    ;线路参数对称化
for i=1: 36
    for j=1: 36
        if t(i,j)~=0
        yt(i,j)=1/(t(i,j)*tk(i,j));
        end
    end
end     %  形成变压器支路导纳，未对称化。

for i=1: 36
    for j=1: 36
        if yt(i,j)~=0
        yt(j,i)=yt(i,j);
        end
    end
end           %变压器支路互导纳对称化。


for i=1: 36
    yt0(i)=0;
    y(i)=0;
    y0(i)=0;
    ty(i)=0;
end

for i=1: 36
       for j=1: 36
        if t(i,j)~=0 ;   %如果i,j之间是变压器支路，节点i和节点j分别计算对地导纳。
        yt0(i)=yt0(i)+(tk(i,j)-1)/(tk(i,j)*t(i,j));
        yt0(j)=yt0(j)+(1-tk(i,j))/(tk(i,j)*tk(i,j)*t(i,j));
         
        end
    end
end                                %  计算节点i变压器支路对地导纳

for i=1: 36
    for j=1: 36
        y(i)=y(i)+ld(i,j);
        y0(i)=y0(i)+l0(i,j);
        ty(i)=ty(i)+yt(i,j);
    end
end                        %;计算节点i线路互导纳、对地导纳、变压器支路导纳之和。

 

for i=1: 36
    Y(i,i)=y(i)+y0(i)+ty(i)+yt0(i);
end                                % 计算对角线元素

for i=1: 36
    for j=1: 36
        if i~=j
            Y(i,j)=-yt(i,j)-ld(i,j);           
        end
    end
end

B=imag(Y);
G=real(Y);

 

 

 

 


n=36;
m=31;
qs(9)=-2.21;    ps(9)=-3.76;%原为-3.76

 qs(16)=-2.3;    ps(16)=-5; %%
% qs(16)=-2.3 ;    ps(16)=-5 ;

qs(18)=-2.2;    ps(18)=-4.3;%%

qs(19)=-0.662; ps(19)=-0.864;
qs(20)=-0.474; ps(20)=-0.719;
qs(21)=-0.5;   ps(21)=-0.7;
qs(22)=-1.69;  ps(22)=-2.265;
qs(23)=-1.44;  ps(23)=-2.87;
qs(29)=-0.1*1.;  ps(29)=-5.2;%原为-5.2
  
qs(4)=0.7; %原为0.7
              ps(4)=1.6;
qs(5)=3.34;   ps(5)=4.3;
qs(2)=3.6;    %ps(2)=6*1.8;
              ps(2)=6;
              ps(32)=3.1; % 原为3.1   change to 6.1
              ps(33)=-0.01;
              ps(34)=2.25;
              ps(35)=3.06; %原为3.06  change to 8.5
              
              
              qs(37)=0; ps(37)=0;
              q(37)=0;  p(37)=0;


  


for i=1: n
    V(i)=1;
    delta(i)=0;
end   % 所有节点电压初值
%delta(36)=-0.;
%V(3)=1.1;
%V(4)=1.05;
 
for i=1: m+n-1
    xiuz(i)=0;
end   %修正量初值

%dpq1(1)=1;
%while dpq1(1)> 0.0001

for kk=1:16
    for i=1: m
    V(i)=V(i)+xiuz(n-1+i)*V(i);
    end   %修正节点电压幅值
    for i=1: n-1
    delta(i)=delta(i)+xiuz(i);
    end   %修正节点电压相角

    for i=1: n
    tem1(i)=0;
    end

    for i=1: n
        for j=1: n
        tem1(i)=tem1(i)+V(j)*(G(i,j)*cos(delta(i)-delta(j))+B(i,j)*sin(delta(i)-delta(j)));
        end
        p(i)=V(i)*tem1(i);      %  计算各节点有功注入，形成雅可比矩阵用
    end    
    
    for i=1: n-1
         dpq(i)=ps(i)-p(i);
     end     % 计算有功不平衡量 

     for i=1: n
    tem2(i)=0;
    end
   
    for i=1: n
        for j=1: n
        tem2(i)=tem2(i)+V(j)*(G(i,j)*sin(delta(i)-delta(j))-B(i,j)*cos(delta(i)-delta(j)));
        end
        q(i)=V(i)*tem2(i);         %  计算各节点无功注入，形成雅可比矩阵用
    end   
    
    for i=1: m
         dpq(n-1+i)=qs(i)-q(i);
     end        % 计算无功不平衡量


    for i=1: n-1
       for j=1: m
          if i==j
              J(i,j+n-1)=-V(i)*V(i)*G(i,j)-p(i);
          else
            J(i,j+n-1)=-V(i)*V(j)*(G(i,j)*cos(delta(i)-delta(j))+B(i,j)*sin(delta(i)-delta(j)));
           end
       end
     end       %计算雅可比矩阵 N；

    for i=1: m
        for j=1: n-1
           if i==j
            J(i+n-1,j)=V(i)*V(i)*G(i,j)-p(i);
           else
               J(i+n-1,j)=V(i)*V(j)*(G(i,j)*cos(delta(i)-delta(j))+B(i,j)*sin(delta(i)-delta(j)));
           end
        end
     end       %计算雅可比矩阵 K；

    for i=1: m
         for j=1: m
             if i==j
                J(n-1+i,n-1+j)=V(i)*V(i)*B(i,i)-q(i);
             else
                  J(n-1+i,n-1+j)=-V(i)*V(j)*(G(i,j)*sin(delta(i)-delta(j))-B(i,j)*cos(delta(i)-delta(j)));
             end
         end
      end       %计算雅可比矩阵 L；

    for i=1: n-1
         for j=1: n-1
            if i==j
               J(i,j)=V(i)*V(i)*B(i,j)+q(i);        %计算雅可比矩阵 H；
            else
                J(i,j)=-V(i)*V(j)*(G(i,j)*sin(delta(i)-delta(j))-B(i,j)*cos(delta(i)-delta(j)));
            end
         end
     end

   for i=1: m+n-1
       dpq1(i)=dpq(i);
   end
    
   xiuz=(-inv(J)*dpq1');
   delta1=delta*180/3.14159265;
   
   
   %
   
end


z=inv(Y);%阻抗矩阵备用



36节点系统算例潮流计算结果(编号与原接线不同)

电压

0.9914    0.9903    1.0195    1.0094    1.0333    0.9887    0.9895    0.9948    0.9899    0.9895    1.0196    1.0117   1.0010    0.9885    0.9887    0.9915    0.9672    0.9948    0.9945    1.0026    0.9988    1.0392    0.9777    1.0464    1.0196    1.0118    1.0117    1.0011    0.9914    1.0195    1.0195    1.0000    1.0000    1.0000    1.0000    1.0000
相位

-41.5130  -11.9921  -16.4671  -22.8737  -31.8242  -30.0071  -19.7508      -41.5325  -20.2014  -19.7508  -23.7891  -31.1677 -38.2770  -30.1359      -30.0071  -41.5135  -41.2169  -41.5325  -28.8754  -37.9659  -31.5197     -30.0790  -20.1593   -5.5066  -23.8103  -31.1464  -31.1951  -38.2493      -41.5284  -16.4535  -16.4519  -27.8000  -41.2368  -10.7673               -10.6603                  0