clc
clear all;
format long;
Z=[0 1 0.0922+j*0.047 100+j*60
1 2 0.0493+j*0.2511 90+j*40
2 3 0.366+j*0.1864 120+j*80
3 4 0.3811+j*0.1941 60+j*30
4 5 0.8190+j*0.707 60+j*20
5 6 0.1872+j*0.6188 200+j*100
6 7 0.7114+j*0.2351 200+j*100
7 8 1.03+j*0.74 60+j*20
8 9 1.044+j*0.74 60+j*20
9 10 0.1966+j*0.065 45+j*30
10 11 0.3744+j*0.1238 60+j*35
11 12 1.468+j*1.155 60+j*35
12 13 0.5416+j*0.7129 120+j*80
13 14 0.5910+j*0.526 60+j*10
14 15 0.7463+j*0.5450 60+j*20
15 16 1.289+j*1.7210 60+j*20
16 17 0.7320+j*0.574 90+j*40
1 18 0.164+j*0.1565 90+j*40
18 19 1.5042+j*1.3554 90+j*40
19 20 0.4095+j*0.4784 90+j*40
20 21 0.7089+j*0.9373 90+j*40
2 22 0.4512+j*0.3083 90+j*50
22 23 0.8980+j*0.7091 420+j*200
23 24 0.8960+j*0.7011 420+j*200
5 25 0.2030+j*0.1034 60+j*25
25 26 0.2842+j*0.1447 60+j*25
26 27 10.59+j*0.9337 60+j*20
27 28 0.8042+j*0.7006 120+j*70
28 29 0.5075+j*0.2585 200+j*600
29 30 0.9744+j*0.9630 150+j*70
30 31 0.3105+j*0.3619 210+j*100
31 32 0.3410+j*0.5302 60+j*40];%33nodes datas
Sb=1;
Vb=10.5;
Zb=Vb^2/Sb;
Z(:,3)=Z(:,3)/Zb;
Z(:,4)=Z(:,4)/Sb/1000;%handle datas
V=ones(32,1);%set V
for i=1:32
I=(Z(:,4).')*(diag(1./V));
I=I';
I1=I;
for k=32:-1:1 %后推电流
A=(find(Z(:,1)==k));
if size(A)==0
I(k,1);
else
I(k,1)=sum(I(A,1))+I(k,1); %
end
end
V(1,1)=1.1-Z(1,3)*I(1,1);%前推电压
for L=2:32
B=find(Z(:,2)==L);%
V(L,1)=V(Z(B,1),1)-Z(B,3)*I(B,1);
end
dS=Z(:,4)-((diag(V(:,1)))*conj(I1(:,1)));
dP=real(dS);
dQ=imag(dS);
dPQ=[dP;dQ];
if max(abs(dPQ))<0.000001
break;
else
i=i+1;
end
end
Z(:,3)
I
V
V=abs(V)
i
