ljf_chaoliujisuan_090326.m

正确的潮流计算程序,经多个算例验证了正确性,值得一看

clear
clc
n=5;%input('请输入节点数:n=');
nl=7;%input('请输入支路数:nl=');
isb=1;%input('请输入平衡母线节点号:isb=');
pr=1e-5;%input('请输入误差精度:pr=');
B1=[1,2,0.02+j*0.06,0,1,0;1,3,0.08+j*0.24,0,1,0;2,3,0.06+j*0.18,0,1,0;2,4,0.06+j*0.18,0,1,0;2,5,0.04+j*0.12,0,1,0;3,4,0.01+j*0.03,0,1,0;4,5,0.08+j*0.24,0,1,0];%input('请输入由支路参数形成的矩阵:B1=');变压器侧为1，否则为0
B2=[0,0,1.06,1,0,1;0,0.2+j*0.2,1,1,0,2;0.45+j*0.15,0,1,1,0,2;0.4+j*0.05,0,1,1,0,2;0.6+j*0.1,0,1,1,0,2];%input('请输入各节点参数形成的矩阵:B2=');
X=[1,0;2,0;3,0;4,0;5,0];%input('请输入由节点号及其对地阻抗形成的矩阵:X=');
Y=zeros(n);U=zeros(1,n);cta=zeros(1,n);V=zeros(1,n);O=zeros(1,n);S1=zeros(nl);
for i=1:n
    if X(i,2)~=0;
        p=X(i,1);
        Y(p,p)=X(i,2);
    end
end
for i=1:nl
    if B1(i,6)==0
        p=B1(i,1);q=B1(i,2);
    else p=B1(i,2);q=B1(i,1);
    end
    Y(p,q)=Y(p,q)-1./(B1(i,3)*B1(i,5));
    Y(q,p)=Y(p,q);
    Y(q,q)=Y(q,q)+1./(B1(i,3)*B1(i,5))+(B1(i,4)*B1(i,5))./2;
    Y(p,p)=Y(p,p)+1./(B1(i,3)*B1(i,5))+(B1(i,4)*B1(i,5))./2;
end %求导纳矩阵
G=real(Y);B=imag(Y);
for i=1:n
    cta(i)=angle(B2(i,3));%电压角度
    U(i)=abs(B2(i,3));
    %V(i)=B2(i,4);
end
for i=1:n
    S(i)=B2(i,2)-B2(i,1);
    B(i,i)=B(i,i)+B2(i,5);
end
P=real(S);Q=imag(S);
ICT1=0;IT2=1;
while IT2~=0
    IT2=0;t1=1;t2=1;
    for i=1:n
        if i~=isb
            C(i)=0;
            D(i)=0;
            for j1=1:n
                C(i)=C(i)+U(i)*U(j1)*(G(i,j1)*cos(cta(i)-cta(j1))+B(i,j1)*sin(cta(i)-cta(j1)));
                D(i)=D(i)+U(i)*U(j1)*(G(i,j1)*sin(cta(i)-cta(j1))-B(i,j1)*cos(cta(i)-cta(j1)));
            end
            DP(t1)=P(i)-C(i);
            t1=t1+1;
            if B2(i,6)==2
                DQ(t2)=Q(i)-D(i);
                t2=t2+1;
            end 
        end
    end
    t1=t1-1;t2=t2-1;
    DPQ=[DP';DQ']; %求DP,DQ
    for i=1:t1+t2
        if abs(DPQ(i))>pr
            IT2=IT2+1;
        end
    end
    H=zeros(t1,t1);N=zeros(t1,t2);K=zeros(t2,t1);L=zeros(t2,t2);
    for i=1:t1+1
        for j1=1:t1+1
            if j1~=isb&j1~=i
               H(i,j1)=U(i)*U(j1)*(G(i,j1)*sin(cta(i)-cta(j1))-B(i,j1)*cos(cta(i)-cta(j1)));
           elseif j1~=isb&j1==i
               H(i,j1)=-U(i)^2*B(i,j1)-D(i);
           end
        end
    end
    for i=1:t1+1
        for j1=1:t2+1
            if j1~=isb&j1~=i
                N(i,j1)=U(i)*U(j1)*(G(i,j1)*cos(cta(i)-cta(j1))+B(i,j1)*sin(cta(i)-cta(j1)));
            elseif j1~=isb&j1==i
                N(i,j1)=U(i)^2*G(i,j1)+C(i);
            end
        end
    end
    for  i=1:t2+1
        for j1=1:t1+1
            if j1~=isb&j1~=i
                K(i,j1)= -U(i)*U(j1)*(G(i,j1)*cos(cta(i)-cta(j1))+B(i,j1)*sin(cta(i)-cta(j1)));
            elseif j1~=isb&j1==i
                K(i,j1)=-U(i)^2*G(i,j1)+C(i);
            end
        end
    end
    for i=1:t2+1
        for j1=1:t2+1
            if j1~=isb&j1~=i
                L(i,j1)=U(i)*U(j1)*(G(i,j1)*sin(cta(i)-cta(j1))-B(i,j1)*cos(cta(i)-cta(j1)));
            elseif j1~=isb&j1==i
                L(i,j1)=-U(i)^2*B(i,j1)+D(i);
            end
        end
    end
    H(isb,:)=[];
    H(:,isb)=[];
    N(isb,:)=[];
    N(:,isb)=[];
    K(isb,:)=[];
    K(:,isb)=[];
    L(isb,:)=[];
    L(:,isb)=[];
    J=[H,N;K,L];%求雅可比矩阵
    modify=J\DPQ;
    Dcta=modify(1:t1);
    U1=U;
    U1(isb)=[];
    DU=diag(U1,0)*modify([t1+1:t1+t2],:);
    t4=1;
    for i=1:n
        if i~=isb
            cta(i)=cta(i)+Dcta(t4);
            t4=t4+1;
        end
    end
    t5=1;
    for i=1:n
        if B2(i,6)==2
            U(i)=U(i)+DU(t5);
            t5=t5+1;
        end
    end
    ICT1=ICT1+1;
end   %修正原值
cta1=cta*180/pi;
for i=1:n
    UU(i)=U(i)*cos(cta(i))+j*U(i)*sin(cta(i));
end

%求平衡节点功率
Ip=0;
for q=1:n
    Ip=Ip+conj(Y(isb,q))*conj(UU(q));
end
Sp=UU(isb)*Ip;
%disp('--------------------------------------------------------------------------------');
%disp('按公式计算全部线路功率，结果如下:');
for i=1:nl
    if B1(i,6)==0
        p=B1(i,1);q=B1(i,2);
    else p=B1(i,2);q=B1(i,1);
    end
    Si(p,q)=UU(p)*(conj(UU(p))*conj(B1(i,4)./2)+(conj(UU(p))-conj(UU(q)))*conj(1./(B1(i,3)*B1(i,5))));%各条支路首端功率Si
    Sj(q,p)=UU(q)*(conj(UU(q))*conj(B1(i,4)./2)+(conj(UU(q))-conj(UU(p)))*conj(1./(B1(i,3)*B1(i,5))));%各条支路末端功率Si
    %Si(q,p)=UU(q)*(conj(UU(q))*conj(B1(i,4)./2)+(conj(UU(q))-conj(UU(p)))*conj(1./(B1(i,3)*B1(i,5))));
    f=[p,q,Si(p,q)];
    %disp(f);
end
dSz=Sp+sum(S);%网络总损耗
for i=1:nl
    if B1(i,6)==0
        p=B1(i,1);q=B1(i,2);
    else p=B1(i,2);q=B1(i,1);
    end
    DS(i)=Si(p,q)+Sj(q,p);%各条支路功率损耗DS
    %(DS(i));
end
for p=1:n                                      %各节点功率，注入为正，吸收为负
    c(p)=0;
    for q=1:n
    c(p)=c(p)+conj(Y(p,q))*conj(UU(q));
    end
    s(p)=UU(p)*c(p);
end