case39连续潮流

基于IEEE case39编写的连续潮流程序，用于计算电力系统的稳态电压裕度。

% function lianxuchaoliu
clear;
clc;
n=9;%节点数；
nl=9;%支路数；
isb=1;%平衡节点号；
pr=0.00001;%误差精度；
%依次是支路首端；末端，支路阻抗；对地电纳；支路变比；折算到哪一侧标志（高压侧为1；低压侧为0）；
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];
%依次是节点的发电机功率；负荷功率；节点电压初值；PV节点电压V给定值；节点无功补偿设备容量；节点分类标号（平衡1；PQ2；PV3）；
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];
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）;
end
for i=1:n
    S（i）=b2（i1）-b2（i2）;
    B（ii）=B（ii）+b2（i5）;
end
P=real（S）;Q=imag（S）;
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）;

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件      14448  2011-03-07 21:02  lianxuchaoliu.m

----------- ---------  ---------- -----  ----

                14448                    1