clear;
clc;
n=9;%節點數；
nl=9;%支路數；
isb=1;%平衡節點號；
pr=0.00001;%誤差精度；
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];
%依次是支路首端；末端，支路阻抗；對地電納；支路變比；折算到哪一側標志（高壓側為1；低壓側為0）；
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];
%依次是節點的發電機功率；負荷功率；節點電壓初值；PV節點電壓V給定值；節點無功補償設備容量；節點分類標號（平衡1；PQ2；PV3）；
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）;??????????%G：導納矩陣實部?B:導納矩陣虛部
for?i=1:n
????e（i）=real（b2（i3））;???????%節點電壓實部
????f（i）=imag（b2（i3））;???????%節點電壓虛部
????v（i）=b2（i4）;?????????????%PV節點電壓V給定值
end
for?i=1:n
????S（i）=b2（i1）-b2（i2）;?????%節點的發電機功率-負荷功率
????B（ii）=B（ii）+b2（i5）;????%節點無功補償設備容量?
end
P=real（S）;Q=imag（S）;??????????%ΔP有功?ΔQ無功?
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-1）=P（i）-e（i）*C-f（i）*D;
????????????w（2*i）=Q（i）-f（i）*C+e（i）*D;
????????else?if?b2（i6）==3%P，V節點；
????????????????w（2*i-1）=P（i）-e（i）*C-f（i）*D;
????????????????w（2*i）=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）=-G（ij）*e（i）-B（ij）*f（i）;??%-M（ij）
????????????????????Jac（2*i-12*j）=G（ij）*e（i）+B（ij）*f（i）;???????%-N（ij）
????????????????????Jac（2*i2*j）=G（ij）*f（i）-B（ij）*e（i）;?????????%-L（ij）
????????????????????Jac（2*i-12*j-1）=G（ij）*f（i）-B（ij）*e（i）;?????%-H（ij）
????????????????else?if?j==i
????????????????????????m=0;h=0;
????????????????????????for?r=1:n
????????????????????????????m=m+G（ir）*e（r）-B（ir）*f（r）;?????????%a（i）
????????????????????????????h=h+G（ir）*f（r）+B（ir）*e（r）;?????????%b（i）
????????????????????????end
????????????????????????Jac（2*i2*j-1）=m-G（ii）*e（i）-B（ii）*f（i）;???%-M（ij）
????????????????????????Jac（2*i-12*j）=m+G（ii）*e（i）+B（ii）*f（i）;???%-N（ij）
????????????????????????Jac（2*i2*j）=-1*h-B（ii）*e（i）+G（ii）*f（i）;??%-L（ij）
????????????????????????Jac（2*i-12*j-1）=h-B（ii）*e（i）+G（ii）*f（i）;?%-H（ij）
????????????????????end
????????????????end
????????????else?if?b2（i6）==3%P，V節點；
????????????????????if?j~=i
????????????????????