%用于一個平衡節點，其他為PQ節點的多節點系統。
clear
%=======================請輸入數據==================================
JDS=14    %輸入系統的節點數
 %線路支路數據(包括普通的線路阻抗和變壓器阻抗，線路與變壓器的模型均采用Pi型)  
%線路兩端的節點號，線路（或變壓器）阻抗，線路（或變壓器）對地導納1/2*B
SN=100
ZL=[1 2 0.01938 0.05917 0.02640
    2 3 0.04699 0.19797 0.02190
    2 4 0.05811 0.17632 0.01870
    1 5 0.05403 0.22304 0.03460
    2 5 0.05695 0.17388 0.01700
    3 4 0.06701 0.17103 0.01730
    4 5 0.01135 0.04211 0.00640
    5 6 0.00000 0.25202 0.00000
    4 7 0.00000 0.20912 0.00000
    7 8 0.00000 0.17615 0.00000
    4 9 0.00000 0.55618 0.00000
    7 9 0.00000 0.11001 0.00000
    9 10 0.03181 0.08450 0.00000
    6 11 0.09498 0.19890 0.00000
    6 12 0.12291 0.15581 0.00000
    6 13 0.06615 0.13027 0.00000
    9 14 0.12711 0.27038 0.00000
    10 11 0.08205 0.19207 0.00000
    12 13 0.22092 0.19988 0.00000
    13 14 0.17093 0.34802 0.00000]                            
 %符荷功率流出-節點號，輸出有功，輸出無功,注意流出系統為正
PQF=[2 21.70 12.07
    3 94.20 19.00
    4 47.80 -3.90
    5 7.60 1.60
    6 11.20 7.50
    7 0.00 0.00
    8 0.00 0.00
    9 29.50 16.60
    10 9.00 5.80
    11 3.50 1.80
    12 6.10 1.60
    13 13.50 5.80
    14 14.90 5.90]
PQG=[2 40.00 42.40
    3 0.00 23.39
    4 0.00 0.00
    5 0.00 0.00
    6 0.00 12.24
    7 0.00 0.00
    8 0.00 17.36
    9 0.00 0.00
    10 0.00 0.00
    11 0.00 0.00
    12 0.00 0.00
    13 0.00 0.00
    14 0.00 0.00]
 S=[1 1.06 0]           %定義平衡節點-節點號，電壓，角度
%=======================數據讀入完成==================================
V=ones(1,JDS)                   %預設各個節點的電壓幅值為1,各個節點的相角為0,輸入有功無功為0
AN=zeros(1,JDS)                
Ps=zeros(JDS,1)
Qs=zeros(JDS,1)
V(1,1)=S(1,2)                   %平衡節點的幅值
AN(1,1)=S(1,3)                  %平衡節點的角度
for n=1:size(PQF(:,1))
    PQF(n,2)=PQF(n,2)/SN
    PQF(n,3)=PQF(n,3)/SN
end
for n=1:size(PQG(:,1))
    PQG(n,2)=PQG(n,2)/SN
    PQG(n,3)=PQG(n,3)/SN
end

if size(PQF)~=0
 for x=1:size(PQF(:,1))                  %讀入各個節點的功率輸出為負值，輸入為正值
    Ps(PQF(x,1),1)=Ps(PQF(x,1),1)-PQF(x,2)
    Qs(PQF(x,1),1)=Qs(PQF(x,1),1)-PQF(x,3)
 end
end
if size(PQG)~=0
 for x=1:size(PQG(:,1))                  
    Ps(PQG(x,1),1)=Ps(PQG(x,1),1)+PQG(x,2)
    Qs(PQG(x,1),1)=Qs(PQG(x,1),1)+PQG(x,3)
 end
end
for x=1:JDS
    Y(x,x)=0
end
for x= 1:size(ZL(:,1))                   %取K1K2為每線路的兩端點節點號
    k1=ZL(x,1)
    k2=ZL(x,2) 
    
    YA(k1,k2)=1/(ZL(x,3)+ZL(x,4)*i)             %計算導納絕對值，然后推算出導納陣
    Y(k1,k2)=-YA(k1,k2)
    Y(k2,k1)=-YA(k1,k2)
    Y(k1,k1)=Y(k1,k1)+YA(k1,k2)+ZL(x,5)*i
    Y(k2,k2)=Y(k2,k2)+YA(k1,k2)+ZL(x,5)*i
    
    B1(k1,k2)=imag(Y(k1,k2))                 %取B1時考慮線路電阻，計算出導納Y，然后取G+jB的B部分
    B1(k2,k1)=imag(Y(k1,k2))
    B1(k1,k1)=B1(k1,k1)-B1(k1,k2)
    B1(k2,k2)=B1(k2,k2)-B1(k1,k2)
        
    B21(k1,k2)=1/ZL(x,4)              %B2是忽略電阻，只算電抗的倒數
    B21(k2,k1)=1/ZL(x,4)
    B21(k1,k1)=B21(k1,k1)-1/ZL(x,4)+ZL(x,5)
    B21(k2,k2)=B21(k2,k2)-1/ZL(x,4)+ZL(x,5)
end
    B11=B1(2:JDS,2:JDS)%B1去掉平衡節點的行和列     
    B22=B21(2:JDS,2:JDS)
    G=real(Y)
    B=imag(Y)      
    %到此，算出了初步的各節點電壓相角輸入有功無功，算出了Y，B1和B2
%===================================================================             
for C=1:50
    for n=1:JDS
        w=0
        for k=1 : JDS   %Pi=Ui*SUM(j=1to n)Uj*(GijcosDij+BijsinDij)
            w=w+V(C,k)*(G(n,k)*cos(AN(C,n)-AN(C,k))+B(n,k)*sin(AN(C,n)-AN(C,k))) 
        end
        P(C,n)=V(C,n)*w
        DP(n,C)=Ps(n,1)-P(C,n)                             %P的誤差
        DPV(n,C)=DP(n,C)/V(C,n)                  
    end
    DP(1,C)=0                 %令第一行為0，因為以上計算不算平衡節點在內
    DPV(1,:)=[]  
   %=======================修正方程=============================
    VAN=[]
    VAN=-inv(B11)*DPV(:,C)                      %DPV為DELTA P /U =-B‘ U DELTA D
    VAN=[0;VAN]
    for k=1:JDS
       DAN(k)=VAN(k)/V(C,k)          %算出角度變化值
       AN(C+1,k)=AN(C,k)+DAN(k)
    end
    
    for n=1:JDS
        w=0
        for k=1 : JDS         %Qi=Ui*SUM(j=1to n)Uj*(GijsinDij-BijcosDij)
            w=w+V(C,k)*(G(n,k)*sin(AN(C+1,n)-AN(C+1,k))-B(n,k)*cos(AN(C+1,n)-AN(C+1,k)))
        end
        Q(C,n)=V(C,n)*w
        DQ(n,C)=Qs(n,1)-Q(C,n)                %Q的誤差
        DQV(n,C)=DQ(n,C)/V(C,n)
    end
        DQ(1,C)=0   %把非PQ點的元素全部去掉
        DQV(1,:)=[]
    DV=-inv(B22)*DQV(:,C)            %算出電壓差值
    DV=[0;DV]
    for k=1:JDS        
        n=k
        V(C+1,k)=V(C,k)+DV(n)                   %算出電壓新值，進行下一次迭代
    end 
    if max(abs(DQ(:,C)))<0.00001,max(abs(DP(:,C)))<0.00001            %判斷誤差
        break
    end
end
  %=======================開始計算結果============================
 for k=1:JDS 
    ANS(k,1)=k
    ANS(k,2)=V(C+1,k)                        %電壓幅值
    ANS(k,3)=AN(C+1,k)*180/pi                %電壓相角
    ANS(k,4)=P(C,k)                          %節點有功，正值為輸入，負值為輸出
    ANS(k,5)=Q(C,k)  
 end
%輸出計算結果
clc
Y
B11
B22
ANS