lianxuchaoliu.m

电力系统连续潮流程序

function lianxuchaoliu
clear;
clc;
n=9;%节点数；
nl=9;%支路数；
isb=1;%平衡节点号；
pr=0.00001;%误差精度；
b1=[1 4 0.0576i 0 1.05 1; 4 5 0.017+0.092i 0.158i 1 0; 5 6 0.039+0.17i 0.358i 1 0; 3 6 0.0586i 0 1.05 1; 6 7 0.0119+0.1008i 0.209i 1 0; 7 8 0.0085+0.072i 0.149i 1 0; 8 2 0.0625i 0 1.05 1; 8 9 0.032+0.161i 0.306i 1 0; 9 4 0.01+0.085i 0.176i 1 0];
%依次是支路首端；末端，支路阻抗；对地电纳；支路变比；折算到哪一侧标志（高压侧为1；低压侧为0）；
b2=[0 0 1.05 1.05 0 1; 1.63 0 1.05 1.05 0 3; 0.85 0 1.05 1.05 0 3; 0 0 1 0 0 2; 0 0.9+0.3i 1 0 0 2; 0 0 1 0 0 2; 0 1+0.35i 1 0 0 2; 0 0 1 0 0 2; 0 1.25+0.5i 1 0 0 2];
%依次是节点的发电机功率；负荷功率；节点电压初值；PV节点电压V给定值；节点无功补偿设备容量；节点分类标号（平衡1；PQ2；PV3）；
Y=zeros(n);%求导纳阵；
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)^2)+b1(i,4)./2;
    Y(p,p)=Y(p,p)+1./b1(i,3)+b1(i,4)./2;
end
%disp('系统的导纳阵为：');
%disp(Y);
G=real(Y);B=imag(Y);
for i=1:n
    e(i)=real(b2(i,3));
    f(i)=imag(b2(i,3));
    v(i)=b2(i,4);
end
for i=1:n
    S(i)=b2(i,1)-b2(i,2);
    B(i,i)=B(i,i)+b2(i,5);
end
P=real(S);Q=imag(S);
w=zeros(2*n,1);Jac=zeros(2*n);
t=0;
while t==0
    for i=1:n
    if b2(i,6)~=isb
        C=0;D=0;
        for j=1:n
            C=C+G(i,j)*e(j)-B(i,j)*f(j);
            D=D+G(i,j)*f(j)+B(i,j)*e(j);
        end
        if b2(i,6)==2%P，Q节点；
            w(2*i)=P(i)-e(i)*C-f(i)*D;
            w(2*i-1)=Q(i)-f(i)*C+e(i)*D;
        else if b2(i,6)==3%P，V节点；
                w(2*i)=P(i)-e(i)*C-f(i)*D;
                w(2*i-1)=v(i)^2-(e(i)^2+f(i)^2);
            end
        end
    else
        w(2*i-1)=0;
        w(2*i)=0;
    end
end
%disp(w);
w1=w(3:2*n);
for i=1:n
    for j=1:n
        if b2(i,6)~=isb
            if b2(i,6)==2%P，Q节点；
                if j~=i
                    Jac(2*i,2*j-1)=-1*(G(i,j)*e(i)+B(i,j)*f(i));
                    Jac(2*i-1,2*j)=(G(i,j)*e(i)+B(i,j)*f(i));
                    Jac(2*i,2*j)=B(i,j)*e(i)-G(i,j)*f(i);
                    Jac(2*i-1,2*j-1)=B(i,j)*e(i)-G(i,j)*f(i);
                else if j==i
                        m=0;h=0;
                        for r=1:n
                            m=m+G(i,r)*e(r)-B(i,r)*f(r);
                            h=h+G(i,r)*f(r)+B(i,r)*e(r);
                        end
                        Jac(2*i,2*j-1)=-1*m-G(i,i)*e(i)-B(i,i)*f(i);
                        Jac(2*i-1,2*j)=-1*m+G(i,i)*e(i)+B(i,i)*f(i);
                        Jac(2*i,2*j)=-1*h+B(i,i)*e(i)-G(i,i)*f(i);
                        Jac(2*i-1,2*j-1)=h+B(i,i)*e(i)-G(i,i)*f(i);
                    end
                end
            else if b2(i,6)==3%P，V节点；
                    if j~=i
                        Jac(2*i,2*j-1)=-1*(G(i,j)*e(i)+B(i,j)*f(i));
                        Jac(2*i-1,2*j)=0;
                        Jac(2*i,2*j)=B(i,j)*e(i)-G(i,j)*f(i);
                        Jac(2*i-1,2*j-1)=0;
                    else if j==i
                            m=0;h=0;
                            for r=1:n
                                m=m+G(i,r)*e(r)-B(i,r)*f(r);
                                h=h+G(i,r)*f(r)+B(i,r)*e(r);
                            end
                            Jac(2*i,2*j-1)=-1*m-G(i,i)*e(i)-B(i,i)*f(i);
                            Jac(2*i-1,2*j)=-2*f(i);
                            Jac(2*i,2*j)=-1*h+B(i,i)*e(i)-G(i,i)*f(i);
                            Jac(2*i-1,2*j-1)=-2*e(i);
                        end
                    end
                end
            end
        else 
            Jac(2*i-1,2*j-1)=0;
            Jac(2*i,2*j)=0;
            Jac(2*i-1,2*j)=0;
            Jac(2*i,2*j-1)=0;
        end
    end
end
%disp(Jac);
Jac1=Jac(3:2*n,3:2*n);
for k=1:2*n-2
    m=0;
    for i=k+1:2*n-2
        m=Jac1(i,k)./Jac1(k,k);
        w1(i)=w1(i)-m*w1(k);
        for j=k+1:2*n-2
            Jac1(i,j)=Jac1(i,j)-m*Jac1(k,j);
        end
    end
end
E(2*n-2)=w1(2*n-2)./Jac1(2*n-2,2*n-2);
for i=2*n-3:-1:1
    c=0;
    for k=i+1:2*n-2
        c=c+Jac1(i,k)*E(k);
        E(i)=(w1(i)-c)./Jac1(i,i);
    end
end
%disp(E);
for i=1:n-1
    e(i+1)=e(i+1)-E(2*i-1);
    f(i+1)=f(i+1)-E(2*i);
end
%disp(e);
%disp(f);
for i=1:n-1
    b(2*i-1)=abs(E(2*i-1));
    b(2*i)=abs(E(2*i));
end
KB=max(b);
%disp(KB);
if KB<pr
    t=1;
else
    t=0;
end
end
%disp(e);
%disp(f);
for i=1:n
    fz(i)=sqrt(e(i)^2+f(i)^2);
end
disp(fz);
JJ1=zeros(2*n-1,2*n-1);Jac1=zeros(2*n-2,2*n-2);ttt=0;ccc=1;%迭代次数；
L=0.05;T=0;%L是步长，T是参数；
while ttt==0
%确定扩展矩阵
for i=1:n-1
    if b2(i+1,6)==3
        JJ1(2*i,2*n-1)=1.5*1.05;
        JJ1(2*i-1,2*n-1)=0;
    else if b2(i+1,6)==2
            if b2(i+1,2)==0
                JJ1(2*i-1,2*n-1)=0;
                JJ1(2*i,2*n-1)=0;
            else
                JJ1(2*i,2*n-1)=-real(b2(i+1,2));
                JJ1(2*i-1,2*n-1)=-imag(b2(i+1,2));
            end
        end
    end
end
for i=1:n
    for j=1:n
        if b2(i,6)~=isb
            if b2(i,6)==2% P,Q节点；
                if j~=i
                    Jac(2*i,2*j-1)=-1*(G(i,j)*e(i)+B(i,j)*f(i));
                    Jac(2*i-1,2*j)=(G(i,j)*e(i)+B(i,j)*f(i));
                    Jac(2*i,2*j)=B(i,j)*e(i)-G(i,j)*f(i);
                    Jac(2*i-1,2*j-1)=B(i,j)*e(i)-G(i,j)*f(i);
                else if j==i
                        m=0;h=0;
                        for r=1:n
                            m=m+G(i,r)*e(r)-B(i,r)*f(r);
                            h=h+G(i,r)*f(r)+B(i,r)*e(r);
                        end
                        Jac(2*i,2*j-1)=-1*m-G(i,i)*e(i)-B(i,i)*f(i);
                        Jac(2*i-1,2*j)=-1*m+G(i,i)*e(i)+B(i,i)*f(i);
                        Jac(2*i,2*j)=-1*h+B(i,i)*e(i)-G(i,i)*f(i);
                        Jac(2*i-1,2*j-1)=h+B(i,i)*e(i)-G(i,i)*f(i);
                    end
                end 
            else if b2(i,6)==3%P,V节点，
                    if j~=i
                        Jac(2*i,2*j-1)=-1*(G(i,j)*e(i)+B(i,j)*f(i));
                        Jac(2*i-1,2*j)=0;
                        Jac(2*i,2*j)=B(i,j)*e(i)-G(i,j)*f(i);
                        Jac(2*i-1,2*j-1)=0;
                    else if j==i
                            m=0;h=0;
                            for r=1:n
                                m=m+G(i,r)*e(r)-B(i,r)*f(r);
                                h=h+G(i,r)*f(r)+B(i,r)*e(r);
                            end
                            Jac(2*i,2*j-1)=-1*m-G(i,i)*e(i)-B(i,i)*f(i);
                            Jac(2*i-1,2*j)=-2*f(i);
                            Jac(2*i,2*j)=-1*h+B(i,i)*e(i)-G(i,i)*f(i);
                            Jac(2*i-1,2*j-1)=-2*e(i);
                        end
                    end
                end
            end
        else 
            Jac(2*i-1,2*j-1)=0;
            Jac(2*i,2*j)=0;
            Jac(2*i-1,2*j)=0;
            Jac(2*i,2*j-1)=0;
        end
    end
end
Jac1=Jac(3:2*n,3:2*n);
for i=1:2*n-2
    for j=1:2*n-2
        JJ1(i,j)=Jac1(i,j);
    end
end
if ccc==1%以负荷为连续参数；
    for j=1:2*n-2
        JJ1(2*n-1,j)=0;
    end
    JJ1(2*n-1,2*n-1)=1;
    w2=zeros(2*n-1,1);
    for i=1:2*n-2
        w2(i,1)=0;
    end
    w2(2*n-1,1)=1;
end
%disp(JJ1);
if ccc>1%以切向量中分量最大（绝对值最大）的状态变量选定为连续参数；
    for i=1:2*n-2
        Jd(i)=abs(d(i));
    end
    for i=1:2*n-3
        if Jd(i)>=Jd(i+1)
            zd=Jd(i);ek=i;
        else if Jd(i)<=Jd(i+1)
                zd=Jd(i+1);ek=i+1;
            end
        end
    end
    for j=1:2*n-1
        JJ1(2*n-1,j)=0;
    end
    JJ1(2*n-1,ek)=1;
    w2=zeros(2*n-1,1);
    for i=1:2*n-1
        w2(i,1)=0;
    end
    if d(ek)>0
        w2(ek,1)=1;
    else if d(ek)<0
            w2(ek,1)=-1;
        end
    end
end
for k=1:2*n-1
    m=0;
    for i=k+1:2*n-1
        m=JJ1(i,k)./JJ1(k,k);
        w2(i)=w2(i)-m*w2(k);
        for j=k+1:2*n-1
            JJ1(i,j)=JJ1(i,j)-m*JJ1(k,j);
        end
    end
end
d(2*n-1)=w2(2*n-1)./JJ1(2*n-1,2*n-1);
for i=2*n-2:-1:1
    c=0;
    for k=i+1:2*n-1
        c=c+JJ1(i,k)*d(k);
        d(i)=(w2(i)-c)./JJ1(i,i);
    end
end
%disp(d);
for i=1:n-1
    e(i+1)=e(i+1)+L*d(2*i-1);
    f(i+1)=f(i+1)+L*d(2*i);
end
T=T+L*d(2*n-1);
disp(d(2*n-1));
%disp(e);
%disp(f);
%disp(T);
%对预估的近似解进行校正；
tt=0;
while tt==0
    for i=1:n
        if b2(i,6)~=isb
            C=0;D=0;
            for j=1:n
                C=C+G(i,j)*e(j)-B(i,j)*f(j);
                D=D+G(i,j)*f(j)+B(i,j)*e(j);
            end
            if b2(i,6)==2%  P,Q节点；
                if b2(i,2)~=0
                    wj(2*i)=P(i)-T*real(b2(i,2))-e(i)*C-f(i)*D;
                    wj(2*i-1)=Q(i)-T*imag(b2(i,2))-f(i)*C+e(i)*D;
                else
                    wj(2*i)=P(i)-e(i)*C-f(i)*D;
                    wj(2*i-1)=Q(i)-f(i)*C+e(i)*D;
                end
            else if b2(i,6)==3%P,V节点；
                    wj(2*i)=P(i)+T*1.5*1.05-e(i)*C-f(i)*D;
                    wj(2*i-1)=0;
                end
            end
        else
            wj(2*i)=0;
            wj(2*i-1)=0;
        end
    end
    wj1=wj(3:2*n);
    for i=1:n
        for j=1:n
            if b2(i,6)~=isb
                if b2(i,6)==2%P,Q节点；
                    if j~=i
                        JJac(2*i,2*j-1)=-1*(G(i,j)*e(i)+B(i,j)*f(i));
                        JJac(2*i-1,2*j)=(G(i,j)*e(i)+B(i,j)*f(i));
                        JJac(2*i,2*j)=B(i,j)*e(i)-G(i,j)*f(i);
                        JJac(2*i-1,2*j-1)=B(i,j)*e(i)-G(i,j)*f(i);
                    else if j==i
                            m=0;h=0;
                            for r=1:n
                                m=m+G(i,r)*e(r)-B(i,r)*f(r);
                                h=h+G(i,r)*f(r)+B(i,r)*e(r);
                            end
                            JJac(2*i,2*j-1)=-1*m-G(i,i)*e(i)-B(i,i)*f(i);
                            JJac(2*i-1,2*j)=-1*m+G(i,i)*e(i)+B(i,i)*f(i);
                            JJac(2*i,2*j)=-1*h+B(i,i)*e(i)-G(i,i)*f(i);
                            JJac(2*i-1,2*j-1)=h+B(i,i)*e(i)-G(i,i)*f(i);
                        end
                    end
                else if b2(i,6)==3%P,V节点；
                        if j~=i
                            JJac(2*i,2*j-1)=-1*(G(i,j)*e(i)+B(i,j)*f(i));
                            JJac(2*i-1,2*j)=0;
                            JJac(2*i,2*j)=B(i,j)*e(i)-G(i,j)*f(i);
                            JJac(2*i-1,2*j-1)=0;
                        else if j==i
                                m=0;h=0;
                                for r=1:n
                                    m=m+G(i,r)*e(r)-B(i,r)*f(r);
                                    h=h+G(i,r)*f(r)+B(i,r)*e(r);
                                end
                                JJac(2*i,2*j-1)=-1*m-G(i,i)*e(i)-B(i,i)*f(i);
                                JJac(2*i-1,2*j)=-2*f(i);
                                JJac(2*i,2*j)=-1*h+B(i,i)*e(i)-G(i,i)*f(i);
                                JJac(2*i-1,2*j-1)=-2*e(i);
                            end
                        end
                    end
                end
            else
                JJac(2*i-1,2*j-1)=0;
                JJac(2*i,2*j)=0;
                JJac(2*i-1,2*j)=0;
                JJac(2*i,2*j-1)=0;
            end
        end
    end
    JJac1=JJac(3:2*n,3:2*n);
    for k=1:2*n-2
        m=0;
        for i=k+1:2*n-2
            m=JJac1(i,k)./JJac1(k,k);
            wj1(i)=wj1(i)-m*wj1(k);
            for j=k+1:2*n-2
                JJac1(i,j)=JJac1(i,j)-m*JJac1(k,j);
            end
        end
    end
    E1(2*n-2)=wj1(2*n-2)./JJac1(2*n-2,2*n-2);
    for i=2*n-3:-1:1
        c=0;
        for k=i+1:2*n-2
            c=c+JJac1(i,k)*E1(k);
            E1(i)=(wj1(i)-c)./JJac1(i,i);
        end
    end
    %disp(E1);
    for i=1:n-1
        e(i+1)=e(i+1)-E1(2*i-1);
        f(i+1)=f(i+1)-E1(2*i);
    end
    for i=1:n-1
        bx(2*i-1)=abs(E1(2*i-1));
        bx(2*i)=abs(E1(2*i));
    end
    KB1=max(bx);
    if KB1<pr
        tt=1;
    else
        tt=0;
    end
end
%disp(e);
%disp(f);
if d(2*n-1)>0
    ttt=0;ccc=ccc+1;
else
    ttt=1;
end
end
disp(T);
disp(ccc);
%disp(e);
%disp(f);
for i=1:n
    fz1(i)=sqrt(e(i)^2+f(i)^2);
end
disp(fz1);