连续潮流MATLAB程序

本程序为9节点连续潮流计算。 （百度百科）连续潮流需要注意的问题： 1 连续潮流中两个重要的问题是：连续参数的选择和步长的确定。参数选取不正确或是步长不合理都会导致潮流不收敛。 2 连续潮流所求出的负荷裕度是针对特定网络状态和特定的负荷增长方式而言的。当增长方式改变或是网络结构或是参数发生改变后，连续潮流求出的负荷裕度是会改变的。 3 至今为止，出现了基于负荷增长的连续潮流、基于网络断线故障的连续潮流等。

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（i6）==0
        p=b1（i1）;q=b1（i2）;
    else
        p=b1（i2）;q=b1（i1）;
    end
    Y（pq）=Y（pq）-1./（b1（i3）*b1（i5））;
    Y（qp）=Y（pq）;
    Y（qq）=Y（qq）+1./（b1（i3）*b1（i5）^2）+b1（i4）./2;
    Y（pp）=Y（pp）+1./b1（i3）+b1（i4）./2;
end
%disp（‘系统的导纳阵为：‘）;
%disp（Y）;
G=real（Y）;B=imag（Y）;
for i=1:n
    e（i）=real（b2（i3））;%节点电压初值实部
    f（i）=imag（b2（i3））;%节点电压初值虚部
    v（i）=b2（i4）;%PV节点电压V给定值
end
for i=1:n
    S（i）=b2（i1）-b2（i2）;%得到节点传输的功率本节点内发电机发出的的功率减去本节点上的负荷功率
    B（ii）=B（ii）+b2（i5）;%加上节点无功补偿设备容量
end
P=real（S）;Q=imag（S）;%p为传输的有功功率，q为传输的无功功率
w=zeros（2*n1）;Jac=zeros（2*n）;
t=0;
while t==0
    for i=1:n
    if b2（i6）~=isb
        C=0;D=0;
        for j=1:n
            C=C+G（ij）*e（j）-B（ij）*f（j）;
            D=D+G（ij）*f（j）+B（ij）*e（j）;
        end
        if b2（i6）==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（i6）==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（i6）~=isb
            if b2（i6）==2%P，Q节点；
                if j~=i
                    Jac（2*i2*j-1）=-1*（G（ij）*e（i）+B（ij）*f（i））;
                    Jac（2*i-12*j）=（G（ij）*e（i）+B（ij）*f（i））;
                    Jac（2*i2*j）=B（ij）*e（i）-G（ij）*f（i）;
                    Jac（2*i-12*j-1）=B（ij）*e（i）-G（ij）*f（i）;
                else if j==i
                        m=0;h=0;
                        for r=1:n
                            m=m+G（ir）*e（r）-B（ir）*f（r）;
                            h=h+G（ir）*f（r）+B（ir）*e（r）;
                        end
                        Jac（2*i2*j-1）=-1*m-G（ii）*e（i）-B（ii）*f（i）;
                        Jac（2*i-12*j）=-1*m+G（ii）*e（i）+B（ii）*f（i）;
                        Jac（2*i2*j）=-1*h+B（ii）*e（i）-G（ii）*f（i）;
                        Jac（2*i-12*j-1）=h+B（ii）*e（i）-G（ii）*f（i）;
                    end
                end
            else if b2（i6）==3%P，V节点；
                    if j~=i
                        Jac（2*i2*j-1）=-1*（G（ij）*e（i）+B（ij）*f（i））;
                        Jac（2*i-12*j）=0;
                        Jac（2*i2*j）=B（ij）*e（i）-G（ij）*f（i）;