连续潮流matlab代码

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

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）;          %G：导纳矩阵实部 B:导纳矩阵虚部
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-1）=P（i）-e（i）*C-f（i）*D;
            w（2*i）=Q（i）-f（i）*C+e（i）*D;
        else if b2（i6）==3%P，V节点；
                w（2*i-1）=P（i）-e（i）*C-f（i）*D;
                w（2*i）=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）=-G（ij）*e（i）-B（ij）*f（i）;  %-M（ij）
                    Jac（2*i-12*j）=G（ij）*e（i）+B（ij）*f（i）;       %-N（ij）
                    Jac（2*i2*j）=G（ij）*f（i）-B（ij）*e（i）;         %-L（ij）
                    Jac（2*i-12*j-1）=G（ij）*f（i）-B（ij）*e（i）;     %-H（ij）
                else if j==i
                        m=0;h=0;
                        for r=1:n
                            m=m+G（ir）*e（r）-B（ir）*f（r）;         %a（i）
                            h=h+G（ir）*f（r）+B（ir）*e（r）;         %b（i）
                        end
                        Jac（2*i2*j-1）=m-G（ii）*e（i）-B（ii）*f（i）;   %-M（ij）
                        Jac（2*i-12*j）=m+G（ii）*e（i）+B（ii）*f（i）;   %-N（ij）
                        Jac（2*i2*j）=-1*h-B（ii）*e（i）+G（ii）*f（i）;  %-L（ij）
                        Jac（2*i-12*j-1）=h-B（ii）*e（i）+G（ii）*f（i）; %-H（ij）
                    end
                end
            else if b2（i6）==3%P，V节点；
                    if j~=i