关于用matlab计算配电网三相潮流计算，请大家帮忙检验下！
我根据网上的一篇关于单相潮流计算写出下面的程序：
clear
n=33;%节点数
b=32;%支路数
Sb=10;%基准功率
Ub=12.66;%基准电压
Zb=Ub^2/Sb;
%支路自阻抗
Z1=[1,0,1,0.0935+0.0477*i,0.0933+0.0475*i,0.0931+0.0474*i
2,1,2,0.5003+0.2548*i,0.4989+0.2541*i,0.4979+0.2536*i
3,2,3,0.3714+0.1891*i,0.3704+0.1886*i,0.3696+0.1882*i
4,3,4,0.3868+0.1970*i,0.3856+0.1964*i,0.3849+0.1960*i
5,4,5,0.8312+0.7176*i,0.8288+0.7154*i,0.8271+0.7140*i
6,5,6,0.1900+0.6380*i,0.1894+0.6262*i,0.1890+0.6249*i
7,6,7,0.7220+0.2386*i,0.7199+0.2379*i,0.7185+0.2374*i
8,7,8,1.0454+0.7510*i,1.0423+0.7488*i,1.0403+0.7473*i
9,8,9,1.0596+0.7510*i,1.0565+0.7488*i,1.0544+0.7473*i
10,9,10,0.1995+0.0659*i,0.1989+0.0657*i,0.1985+0.0656*i
11,10,11,0.3800+0.1256*i,0.3788+0.1252*i,0.3781+0.1250*i
12,11,12,1.4900+1.1723*i,1.4856+1.1688*i,1.4826+1.1665*i
13,12,13,0.5497+0.7235*i,0.5480+0.7214*i,0.5470+0.7200*i
14,13,14,0.5998+0.5338*i,0.5980+0.5323*i,0.5969+0.5312*i
15,14,15,0.7514+0.5531*i,0.7491+0.5515*i,0.7477+0.5504*i
16,15,16,1.3083+1.7468*i,1.3044+1.7416*i,1.3018+1.7382*i
17,16,17,0.7429+0.5826*i,0.7407+0.5808*i,0.7393+0.5797*i
18,1,18,0.1664+0.1588*i,0.1659+0.1583*i,0.1656+0.1580*i
19,18,19,1.5267+1.3757*i,1.5222+1.3716*i,1.5192+1.3689*i
20,19,20,0.4156+0.4855*i,0.4144+0.4841*i,0.4135+0.4831*i
21,20,21,0.7195+0.9513*i,0.7174+0.9485*i,0.7159+0.9466*i
22,2,22,0.4579+0.3129*i,0.4566+0.3119*i,0.4557+0.3113*i
23,22,23,0.9114+0.7197*i,0.9087+0.7176*i,0.9069+0.7161*i
24,23,24,0.9094+0.7116*i,0.9067+0.7095*i,0.9049+0.7081*i
25,5,25,0.2060+0.1049*i,0.2054+0.1046*i,0.2050+0.1044*i
26,25,26,0.2884+0.1468*i,0.2876+0.1464*i,0.2870+0.1461*i
27,26,27,1.0748+0.9477*i,1.0717+0.9449*i,1.0695+0.9430*i
28,27,28,0.8162+0.7111*i,0.8138+0.7090*i,0.8122+0.7066*i
29,28,29,0.5151+0.2623*i,0.5135+0.2616*i,0.5125+0.2610*i
30,29,30,0.9890+0.9774*i,0.9860+0.9745*i,0.9841+0.9726*i
31,30,31,0.3151+0.3637*i,0.3142+0.3662*i,0.3136+0.3655*i
32,31,32,0.3461+0.5381*i,0.3450+0.5365*i,0.3444+0.5355*i];
Z1_1=Z1(:,[4,5,6])/Zb;
%支路互阻抗
Z2=[1,0,1,0.0009+0.0004*i,0.0013+0.0007*i,0.0011+0.0005*i
2,1,2,0.0049+0.0025*i,0.0073+0.0037*i,0.0059+0.0030*i
3,2,3,0.0036+0.0018*i,0.0054+0.0027*i,0.0043+0.0022*i
4,3,4,0.0038+0.0019*i,0.0057+0.0029*i,0.0045+0.0023*i
5,4,5,0.0081+0.0070*i,0.0122+0.0106*i,0.0098+0.0084*i
6,5,6,0.0018+0.0061*i,0.0028+0.0092*i,0.0022+0.0074*i
7,6,7,0.0071+0.0023*i,0.0106+0.0035*i,0.0085+0.0028*i
8,7,8,0.0103+0.0074*i,0.0154+0.0110*i,0.0123+0.0088*i
9,8,9,0.0104+0.0074*i,0.0156+0.0110*i,0.0125+0.0088*i
10,9,10,0.0019+0.0006*i,0.0029+0.0009*i,0.0023+0.0007*i
11,10,11,0.0037+0.0012*i,0.0056+0.0018*i,0.0044+0.0014*i
12,11,12,0.0146+0.0115*i,0.0220+0.0173*i,0.0176+0.0138*i
13,12,13,0.0054+0.0071*i,0.0081+0.0106*i,0.0064+0.0085*i
14,13,14,0.0059+0.0052*i,0.0088+0.0078*i,0.0070+0.0063*i
15,14,15,0.0074+0.0054*i,0.0111+0.0081*i,0.0088+0.0065*i
16,15,16,0.0128+0.0172*i,0.0193+0.0258*i,0.0154+0.0206*i
17,16,17,0.0073+0.0057*i,0.0109+0.0086*i,0.0087+0.0068*i
18,1,18,0.0016+0.0015*i,0.0024+0.0023*i,0.0019+0.0018*i
19,18,19,0.0150+0.0135*i,0225+0.0203*i,0.0180+0.0162*i
20,19,20,0.0040+0.0047*i,0.0061+0.0071*i,0.0049+0.0057*i
21,20,21,0.0070+0.0093*i,0.0106+0.0140*i,0.0085+0.0112*i
22,2,22,0.0045+0.0030*i,0.0067+0.0046*i,0.0054+0.0036*i
23,22,23,0.0089+0.0070*i,0.0134+0.0106*i,0.0107+0.0085*i
24,23,24,0.0089+0.0070*i,0.0134+0.0105*i,0.0107+0.0084*i
25,5,25,0.0020+0.0010*i,0.0030+0.0015*i,0.0024+0.1044*i
26,25,26,0.0028+0.0014*i,0.0042+0.0021*i,0.0034+0.0017*i
27,26,27,0.0105+0.0093*i,0.0158+0.0140*i,0.0127+0.0112*i
28,27,28,0.0080+0.0070*i,0.0120+0.0105*i,0.0096+0.0084*i
29,28,29,0.0050+0.0025*i,0.0076+0.0038*i,0.0060+0.0031*i
30,29,30,0.0097+0.0096*i,0.0146+0.0144*i,0.0116+0.0115*i
31,30,31,0.0031+0.0036*i,0.0046+0.0054*i,0.0037+0.0043*i
32,31,32,0.0034+0.0053*i,0.0051+0.0079*i,0.0040+0.0063*i];
Z2_2=Z2(:,[4,5,6])/Zb;
%节点功率
Sl=[1,0,1,32+19*i,33+20*i,35+21*i
2,1,2,30+13*i,31+15*i,29+13*i
3,2,3,45+30*i,0+0*i,35+24*i
4,3,4,20+10*i,20+10*i,20+10*i
5,4,5,20+6*i,20+7*i,20+7*i
6,5,6,65+33*i,70+34*i,65+33*i
7,6,7,70+34*i,65+33*i,65+33*i
8,7,8,20+7*i,18+6*i,22+7*i
9,8,9,21+7*i,20+7*i,0+0*i
10,9,10,14+9*i,16+11*i,15+10*i
11,10,11,20+11*i,20+12*i,20+12*i
12,11,12,21+12*i,19+11*i,20+12*i
13,12,13,40+28*i,38+27*i,42+25*i
14,13,14,0+0*i,19+3*i,20+3*i
15,14,15,19+6*i,20+7*i,21+7*i
16,15,16,19+6*i,21+7*i,20+7*i
17,16,17,30+14*i,30+13*i,30+13*i
18,1,18,33+15*i,29+13*i,28+12*i
19,18,19,29+13*i,28+12*i,33+15*i
20,19,20,29+12*i,30+13*i,31+15*i
21,20,21,28+12*i,33+15*i,29+13*i
22,2,22,30+16*i,31+17*i,29+17*i
23,22,23,130+60*i,140+70*i,150+70*i
24,23,24,150+70*i,130+70*i,140+60*i
25,5,25,20+8*i,20+8*i,20+9*i
26,25,26,18+7*i,22+9*i,20+9*i
27,26,27,19+6*i,22+8*i,19+6*i
28,27,28,38+23*i,42+25*i,40+22*i
29,28,29,60+180*i,70+210*i,70+210*i
30,29,30,45+20*i,51+23*i,54+27*i
31,30,31,70+33*i,72+35*i,68+32*i
32,31,32,20+13*i,20+14*i,20+13*i];
Sl_1=Sl(:,[4,5,6])*0.001/Sb;
Z=[Z1_1,Z2_2];
k=0;
V=ones(n,3);
t=0;
%迭代开始
while t<b&k<10
%计算节点注入电流
for l=1:b
j=Z1(l,3);
ua=V(j+1,:);
I(j,:)=conj(Sl_1(j,:)./ua);
end
%回推计算支路电流
J=zeros(b,3);
l=b;
J(l)=J(l)+I(l);
for jj=1:b-1
l=l-1;
for m=l+1:b
if Z1(m,2)==Z1(l,3)
J(l)=J(l)+J(m);
end
end
J(l)=J(l)+I(l);
end
%前推计算节点电压
for l=1:b
j=Z1(l,3)+1;
i=Z1(l,2)+1;
V(j,:)=V(i,:)-([Z(l,1),Z(l,4),Z(l,5);Z(l,4),Z(l,2),Z(l,6);Z(l,5),Z(l,6),Z(l,3)]*J(l,:).').';
end
%收敛判定
t=0;
for j=2:n
SS=V(j,:).*conj(I(j-1,:));
r=max(SS-Sl_1(j-1,:));
dp=abs(real(r));
dq=abs(imag(r));
L1=(dp<0.000001)&(dq<0.000001);
F(j-1,1)=L1;
if L1==1
t=t+1;
end
end
k=k+1;
end
%输出结果
disp('输出各节点电压')
disp(V);
disp('迭代次数')
disp(k);
disp('显示收敛节点情况：‘1’—收敛，‘0’—不收敛');
disp(F);
运行结果为：
输出各节点电压
1.0000 1.0000 1.0000
0.9991 + 0.0001i 1.0000 + 0.0000i 1.0000 + 0.0000i
0.9946 + 0.0005i 0.9999 + 0.0000i 0.9999 + 0.0000i
0.9923 + 0.0008i 0.9999 + 0.0000i 0.9999 + 0.0000i
0.9901 + 0.0011i 0.9999 + 0.0000i 0.9999 + 0.0000i
0.9845 + 0.0006i 0.9998 + 0.0000i 0.9998 + 0.0000i
0.9834 - 0.0006i 0.9998 - 0.0000i 0.9998 - 0.0000i
0.9820 - 0.0004i 0.9998 - 0.0000i 0.9997 - 0.0000i
0.9801 - 0.0007i 0.9998 - 0.0000i 0.9997 - 0.0000i
0.9784 - 0.0009i 0.9998 - 0.0000i 0.9997 - 0.0000i
0.9782 - 0.0009i 0.9998 - 0.0000i 0.9997 - 0.0000i
0.9778 - 0.0008i 0.9998 - 0.0000i 0.9997 - 0.0000i
0.9761 - 0.0011i 0.9998 - 0.0000i 0.9996 - 0.0000i
0.9754 - 0.0015i 0.9998 - 0.0000i 0.9996 - 0.0000i
0.9751 - 0.0016i 0.9998 - 0.0000i 0.9996 - 0.0000i
0.9747 - 0.0017i 0.9998 - 0.0000i 0.9996 - 0.0000i
0.9740 - 0.0021i 0.9997 - 0.0000i 0.9996 - 0.0000i
0.9738 - 0.0021i 0.9997 - 0.0000i 0.9996 - 0.0000i
0.9989 + 0.0000i 1.0000 + 0.0000i 1.0000 - 0.0000i
0.9977 - 0.0004i 1.0000 - 0.0000i 0.8790 + 0.0521i
0.9975 - 0.0005i 1.0000 - 0.0000i 0.8790 + 0.0521i
0.9973 - 0.0006i 1.0000 - 0.0000i 0.8790 + 0.0521i
0.9935 + 0.0003i 0.9999 + 0.0000i 0.9999 + 0.0000i
0.9913 - 0.0002i 0.9999 - 0.0000i 0.9999 - 0.0000i
0.9901 - 0.0005i 0.9999 - 0.0000i 0.9999 - 0.0000i
0.9840 + 0.0008i 0.9998 + 0.0000i 0.9998 + 0.0000i
0.9832 + 0.0011i 0.9998 + 0.0000i 0.9998 + 0.0000i
0.9798 + 0.0014i 0.9998 + 0.0000i 0.9997 + 0.0000i
0.9774 + 0.0018i 0.9998 + 0.0000i 0.9997 + 0.0000i
0.9763 + 0.0022i 0.9998 + 0.0000i 0.9997 + 0.0000i
0.9750 + 0.0018i 0.9998 + 0.0000i 0.9996 + 0.0000i
0.9748 + 0.0017i 0.9998 + 0.0000i 0.9996 + 0.0000i
0.9747 + 0.0016i 0.9998 + 0.0000i 0.9996 + 0.0000i
迭代次数
3
显示收敛节点情况：‘1’—收敛，‘0’—不收敛
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
大家帮我看看这程序对不对,谢谢了
