matlab pq.txt

基于matlab和C++解多平衡节点+普通PQ分解潮流计算问题

%用于一个平衡节点，其他为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