%牛頓法潮流計算程序
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