pqsjcx.m

一个潮流计算软件,可供电力行业有关人员参考使用.

%P-Q法解潮流，对于多接点系统，工作出色！
T=[]      %变压器线路数据
Le=[1 3 0.08 0.24 0
    2 5 0.04 0.12 0
    2 4 0.06 0.18 0
    3 4 0.01 0.03 0
    2 3 0.06 0.18 0
    4 5 0.08 0.24 0
    1 2 0.02 0.06 0]     %线路支路数据
Fd=[3 0.45 0.15
    4 0.4 0.05
    5 0.6 0.1]       %符荷功率
Gr=[2 0.2 0.2
    1 0 0]
Gcv=[]       %PV节点
THSLACK=[1 1.06 0]        %定义平衡节点
S=[]
if max(S) ==0
    S=[]
end
if max(T)==0
    T=[]
end
if max(Gcv)==0
    Gcv=[]
end
%=======================数据读入,形成导纳阵==================================
if size(T)~=0
    Jn1=max(T(:,1))
    Jn3=max(T(:,2))
else
    Jn1=0
    Jn3=0
end
Jn2=max(Le(:,1))
Jn4=max(Le(:,2))
Jn5=max(Fd(:,1))
JN=[Jn1 Jn2 Jn3 Jn4 Jn5]
Jn=max(JN)                     %计算出系统的节点数
Ps=zeros(Jn,1)
Qs=zeros(Jn,1)
if size(Gcv)~=0
    PVn=Gcv(:,1)                   %PV节点号
else
    PVn=0
end
Bs=THSLACK(1,1)                %平衡节点号
V=ones(1,Jn)                   %预设各个节点的电压幅值为1
AN=zeros(1,Jn)                 %预设各个节点的相角为0
if PVn~=0
    for m=1:length(PVn)
        V(1,PVn(m))=Gcv(m,4)
    end
end
V(1,Bs)=THSLACK(1,2)                   %平衡节点的幅值
AN(1,Bs)=THSLACK(1,3)                  %平衡节点的角度


for x=1:size(Fd(:,1))                  %读入各个节点的功率输出为负值，输入为正值
    Ps(Fd(x,1),1)=Ps(Fd(x,1),1)-Fd(x,2)
    Qs(Fd(x,1),1)=Qs(Fd(x,1),1)-Fd(x,3)
end
for x=1:size(Gr(:,1))
    Ps(Gr(x,1),1)=Ps(Gr(x,1),1)+Gr(x,2)
    Qs(Gr(x,1),1)=Qs(Gr(x,1),1)+Gr(x,3)
end
if size(Gcv)~=0
    for x=1:size(Gcv(:,1))
        Ps(Gcv(x,1),1)=Ps(Gcv(x,1),1)+Gcv(x,2)
        Qs(Gcv(x,1),1)=Qs(Gcv(x,1),1)+Gcv(x,3)
    end
end
if size(S)~=0
    for x=1:size(S(:,1))
        Ps(S(x,1),1)=Ps(S(x,1),1)-S(x,2)
        Qs(S(x,1),1)=Qs(S(x,1),1)-S(x,3)
    end
end
[k1,k2]=size(Le)
for x=1:Jn
    Y(x,x)=0
end
for x= 1:k1
    nf=Le(x,1)
    nt=Le(x,2) 
    y(nf,nt)=1/(Le(x,3)+Le(x,4)*i)
    Y(nf,nt)=-y(nf,nt)
    Y(nt,nf)=-y(nf,nt)
    B1(nf,nt)=-imag(y(nf,nt))
    B1(nt,nf)=-imag(y(nf,nt))
    B1(nf,nf)=B1(nf,nf)-B1(nf,nt)
    B1(nt,nt)=B1(nt,nt)-B1(nf,nt)
    B2(nf,nt)=1/Le(x,4)
    B2(nt,nf)=1/Le(x,4)
    B2(nf,nf)=B2(nf,nf)-1/Le(x,4)+Le(x,5)
    B2(nt,nt)=B2(nt,nt)-1/Le(x,4)+Le(x,5)
    Y(nf,nf)=Y(nf,nf)+y(nf,nt)+Le(x,5)*i
    Y(nt,nt)=Y(nt,nt)+y(nf,nt)+Le(x,5)*i
end
[k1,k2]=size(T)
if k1~=0
    for x= 1:k1
        nf=T(x,1)
        nt=T(x,2)
        y(nf,nt)=1/(T(x,3)+T(x,4)*i)
        Y(nf,nt)=-y(nf,nt)/T(x,5)
        Y(nt,nf)=-y(nf,nt)/T(x,5)
        B1(nf,nt)=-imag(y(nf,nt)/T(x,5))
        B1(nt,nf)=-imag(y(nf,nt)/T(x,5))
        B1(nf,nf)=B1(nf,nf)+imag(y(nf,nt)/T(x,5)^2)
        B1(nt,nt)=B1(nt,nt)+imag(y(nf,nt))
        B2(nf,nt)=-1/(T(x,4)*T(x,5))
        B2(nt,nf)=B2(nf,nt)
        B2(nf,nf)=B2(nf,nf)-1/(T(x,4)*T(x,5)^2)
        B2(nt,nt)=B2(nt,nt)-1/T(x,4)
        Y(nf,nf)=Y(nf,nf)+y(nf,nt)/T(x,5)^2
        Y(nt,nt)=Y(nt,nt)+y(nf,nt)
    end
end
    if Bs==Jn
        B1=B1(1:(Jn-1),1:(Jn-1))
    elseif Bs==1
        B1=B1(2:Jn,2:Jn)
    else
        J1=B1(1:(Bs-1),1:(Bs-1))
        J2=B1(Bs+1:Jn,Bs+1:Jn)
        J3=B1(Bs+1:Jn,1:(Bs-1))
        J4=B1(1:(Bs-1),Bs+1:Jn)
        B1=[J1,J4;J3,J2]
    end
    A=ones(1,Jn)
    B21=[]
    for k=1:Jn
        if k==Bs
            A(1,k)=0
        end
        if size(Gcv)~=0
            for x=1:size(Gcv(:,1))
                if k==Gcv(x,1)
                    A(1,k)=0
                end
            end
        end
        if A(1,k)==1
            B21=[B21;B2(k,:)]
        end
    end
    B22=[]
    for k=1:Jn
        if A(1,k)==1
            B22=[B22,B21(:,k)]                  %形成B22矩阵
        end
    end
    
%===================================================================
Z=inv(Y)
G=real(Y)
B=imag(Y)
for C=1:30
    for n=1:Jn
        w=0
        for k=1 : Jn
            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)                  
        DP(Bs,C)=0
    end
    for n=Jn:-1:1
        if n==Bs
            DPV(n,:)=[]
        end
    end
   %=======================修正方程=============================
    VAN=[]
    VAN=-inv(B1)*DPV(:,C)                      %根据修正公式修正Dangle
    if Bs==Jn
        VAN=[VAN;0]
    elseif Bs==1
        VAN=[0;VAN]
    else
        VAN=[VAN(1:(Bs-1));0;VAN((Bs+1):end)]
    end
    for k=1:Jn
       DAN(k)=VAN(k)/V(C,k)
       AN(C+1,k)=AN(C,k)+DAN(k)
    end
    for n=1:Jn
        w=0
        for k=1 : Jn
            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)
        if A(1,n)==0,n==Bs
            DQ(n,:)=0
        end
    end
    
    for n=Jn:-1:1
         if A(1,n)==0,n==Bs
            DQV(n,:)=[]
        end
    end
    DV=-inv(B22)*DQV(:,C)
    
%     if Bs==Jn
%         DV=[DV;0]
%     elseif Bs==1
%         DV=[0;DV]
%     else
%         DV=[DV(1:(Bs-1));0;DV((Bs+1):end)]
%     end
    m=1
    for k=1:Jn
        if A(1,k)==1
            V(C+1,k)=V(C,k)+DV(m)
            m=m+1
        elseif A(1,k)==0
            V(C+1,k)=V(C,k)
        end
    end
    if max(abs(DQ(:,C)))<0.00001,max(abs(DP(:,C)))<0.00001            %判断误差
        break
    end
end
     %V(k)=e(k)+f(k)*i     %电压向量
 for k=1:Jn 
    RE(k,1)=k
    RE(k,2)=V(C+1,k)                        %电压幅值
    RE(k,3)=AN(C+1,k)*180/pi%电压相角
    RE(k,4)=P(C,k)                             %节点有功，正值为输入，负值为输出
    RE(k,5)=Q(C,k)
    RE(k,6)=V(C+1,k)*cos(AN(C+1,k))
    RE(k,7)=V(C+1,k)*sin(AN(C+1,k))
    RE(k,8)=RE(k,6)+i*RE(k,7)          %电压降落的复数表示形式
end
for k1=1:Jn                                   %各个支路电流
    for k2=1:Jn
        if k1~=k2 &Y(k1,k2)~=0
          DUz(k1,k2)=RE(k1,8)-RE(k2,8)       %各个支路的电压降落
          DUzGE(k1,k2)=conj(DUz(k1,k2))       %电压降落的共轭复数
          YzGE(k1,k2)=conj(Y(k1,k2))
          yzGE(k1,k2)=YzGE(k1,k2)*(-1)
          Sz(k1,k2)=RE(k1,8)*DUzGE(k1,k2)*yzGE(k1,k2) %求取各线路功率
          Yz(k1,k2)=Y(k1,k2)
          Iz(k1,k2)=DUz(k1,k2)*Yz(k1,k2)     % 求取个支路电流
          IzGE(k1,k2)=conj(Iz(k1,k2))
          Pz=real(Sz)       %  各个线路上的有功功率
          Qz=imag(Sz)         % 各个线路上的无功功率
        end
    end  
end
 
for k1=1:Jn
    for k2=1:Jn
         
         DSz1(k1,k2)=Sz(k1,k2)+Sz(k2,k1)     %求各线路上的功率损耗
         DSz2(k1,k2)=DUz(k1,k2)*IzGE(k1,k2)    %用另一种方法求取各线路上的功率损耗
         DPz=real(DSz1)      %各个线路上的有功损耗
         DQz=imag(DSz1)       %各个线路上的无功损耗
     end
 end   
 
k3=1
for k1=1:Jn                                  
    for k2=1:Jn
        if Iz(k1,k2)~=0
            Id(k3,1)=k1
            Id(k3,2)=k2
            Id(k3,3)=real(Iz(k1,k2))
            Id(k3,4)=imag(Iz(k1,k2))
            k3=k3+1
        end
    end
end
clc
Y
B
RE
DUz
Iz
Pz
Qz
DPz
DQz