5节点MATLAB潮流计算

%形成节点导纳矩阵
B1=[1 3 0.08 0.24;
    1 2 0.02 0.06;
    2 3 0.06 0.18;
    2 4 0.06 0.18;
    2 5 0.04 0.12;
    3 4 0.01 0.03;
    4 5 0.08 0.24;];
y=zeros（5）;Y=zeros（5）;
for i=1:7
    p=B1（i1）;
    q=B1（i2）;
    y（pq）=1.00/（B1（i3）+B1（i4）*sqrt（-1））;
end
for b=1:5
    for d=（b+1）:5
        y（db）=y（bd）;
    end
end
m=5;n=5;
for i=1:n
    for j=1:n
        if（i==j）
            for m=1:n
                Y（ii）=Y（ii）+y（im）;
            end
        else
            Y（ij）=-y（ij）;
        end
    end
end
%赋初始值
e（2）=1.00;e（3）=1.00;e（4）=1.00;e（5）=1.00;
f（2）=0;f（3）=0;f（4）=0;f（5）=0;
G=real（Y）;B=imag（Y）;
S（2）=0.20+0.20i;S（3）=-0.45-0.15i;S（4）=-0.40-0.05i;S（5）=-0.60-0.1i;
P=real（S）;Q=imag（S）;k=0;p=1;N1=5;
%设置精度和最大循环次数
while p>0.00001&k<100
    e（1）=1.06;
    f（1）=0;   
    for m=1:N1
        for n=1:N1           Pa（n）=（e（m）*（G（mn）*e（n）-B（mn）*f（n））+f（m）*（G（mn）*f（n）+B（mn）*e（n）））;            Qa（n）=（f（m）*（G（mn）*e（n）-B（mn）*f（n））-e（m）*（G（mn）*f（n）+B（mn）*e（n）））;
        end
            Pt（m）=sum（Pa）;Qt（m）=sum（Qa）;
            dP（m）=P（m）-sum（Pa）;
            dQ（m）=Q（m）-sum（Qa）;
    end  
    for m=1:N1
        for n=1:N1
            Bi（n）=G（mn）*f（n）+B（mn）*e（n）;
            Ai（n）=G（mn）*e（n）-B（mn）*f（n）;
        end
        H（mm）=sum（Bi）-（B（mm）*e（m）+G（mm）*f（m））+2*G（mm）*f（m）;
        N（mm）=sum（Ai）-（G（mm）*e（m）-B（mm）*f（m））+2*G（mm）*e（m）;
        J（mm）=-2*B（mm）*f（m）+sum（Ai）-（G（mm）*e（m）-B（mm）*f（m））;
        L（mm）=-2*B（mm）*e（m）-（sum（Bi）-（B（mm）*e（m）+G（mm）*f（m）））;
    end
    for