clc;
clear all;
n=40;%n=input('请输入节点数');
b=39;%b=input('请输入支路数');
%disp('请输入支路阻抗和节点功率矩阵，第一列存支路号，第二列存首节点号，第三列存尾节点号，第四列存支路自阻抗,第五列存尾节点给定功率');
Sb=15;
Ub=23;
Zb=Ub^2/Sb;
Z=[1,0,1,0.169/Zb+i*0.655/Zb,0
2,1,2,0.864/Zb+i*0.751/Zb,0
3,1,3,0.196/Zb+i*0.655/Zb,0
4,2,4,1.374/Zb+i*0.774/Zb,0.522/Sb+i*0.174/Sb
5,2,5,0.864/Zb+i*0.751/Zb,0
6,3,6,0.444/Zb+i*0.439/Zb,0
7,3,7,0.196/Zb+i*0.655/Zb,0
8,5,8,0.864/Zb+i*0.751/Zb,0
9,5,9,0.864/Zb+i*0.751/Zb,0
10,6,10,1.374/Zb+i*0.774/Zb,1.08/Sb+i*0.36/Sb
11,6,11,0.864/Zb+i*0.751/Zb,0
12,7,12,0.279/Zb+i*0.051/Zb,0
13,7,13,0.279/Zb+i*0.051/Zb,0
14,8,14,1.374/Zb+i*0.774/Zb,0.36/Sb+i*0.12/Sb
15,8,15,1.374/Zb+i*0.774/Zb,0.63/Sb+i*0.21/Sb
16,9,16,1.374/Zb+i*0.774/Zb,0.45/Sb+i*0.15/Sb
17,9,17,0.864/Zb+i*0.751/Zb,0
18,11,18,1.374/Zb+i*0.774/Zb,0.54/Sb+i*0.18/Sb
19,11,19,1.374/Zb+i*0.774/Zb,0.675/Sb+i*0.225/Sb
20,12,20,1.374/Zb+i*0.774/Zb,0.45/Sb+i*0.15/Sb
21,12,21,0.444/Zb+i*0.439/Zb,0
22,13,22,0.444/Zb+i*0.439/Zb,1.8/Sb+i*0.6/Sb
23,13,23,0.444/Zb+i*0.439/Zb,0
24,17,24,1.374/Zb+i*0.774/Zb,0.675/Sb+i*0.225/Sb
25,17,25,1.374/Zb+i*0.774/Zb,0.63/Sb+i*0.21/Sb
26,21,26,1.374/Zb+i*0.774/Zb,0.54/Sb+i*0.18/Sb
27,21,27,0.444/Zb+i*0.439/Zb,0
28,23,28,1.374/Zb+i*0.774/Zb,0.855/Sb+i*0.285/Sb
29,23,29,0.864/Zb+i*0.751/Zb,0
30,27,30,1.374/Zb+i*0.774/Zb,0.63/Sb+i*0.21/Sb
31,27,31,0.864/Zb+i*0.751/Zb,0
32,29,32,1.374/Zb+i*0.774/Zb,0.54/Sb+i*0.18/Sb
33,27,33,1.374/Zb+i*0.774/Zb,0.72/Sb+i*0.24/Sb
34,31,34,1.374/Zb+i*0.774/Zb,0.765/Sb+i*0.225/Sb
35,31,35,0.864/Zb+i*0.751/Zb,0
36,35,36,1.374/Zb+i*0.774/Zb,0.54/Sb+i*0.18/Sb
37,35,37,0.864/Zb+i*0.751/Zb,0
38,37,38,1.374/Zb+i*0.774/Zb,0.675/Sb+i*0.225/Sb
39,37,39,1.374/Zb+i*0.774/Zb,0.522/Sb+i*0.174/Sb];
k=0;
V=ones(n,1);
t=0;
%迭代开始处
while t<b & k<10
%算节点注入电流
x1=Z(b,3);x=x1-n;
for l=1:b
j=Z(l,3);
ua=V(j+1,1);
I(j,1)=conj(Z(j,5)/ua);
end
%回推算支路电流
J=zeros(b,1);
l=b;
J(l)=J(l)+I(l);
for jj=1:b-1
l=l-1;
for m=l+1:b
if Z(m,2)==Z(l,3)
J(l)=J(l)+J(m);
end
end
J(l)=J(l)+I(l);
end
%前推算节点电压
for l=1:b
j=Z(l,3)+1;
i=Z(l,2)+1;
V(j,1)=V(i,1)-Z(l,4)*J(l,1);
end
%收敛判定
t=0;
for j=2:n
SS=V(j,1)*conj(I(j-1,1));
dp=real(SS-Z(j-1,5));
dq=imag(SS-Z(j-1,5));
S(j-1,1)=SS;
ddp=abs(dp);
ddq=abs(dq);
L1=(ddp<0.000001)&(ddq<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);
for j=1:b
if F(j,1)==0
disp('显示不收敛节点号、计算功率');
disp(j);disp(S(j,1));
end
end
for j=1:n
Vm(j,1)=abs(V(j,1));Va(j,1)=angle(V(j,1));
end
disp('输出各节点电压幅值');
disp(Vm);
disp('输出各节点电压相角');
disp(Va)
%返回开始处
