function?lianxuchaoliu
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）;
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）=P（i）-e（i）*C-f（i）*D;
????????????w（2*i-1）=Q（i）-f（i）*C+e（i）*D;
????????else?if?b2（i6）==3%P，V節點；
????????????????w（2*i）=P（i）-e（i）*C-f（i）*D;
????????????????w（2*i-1）=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）=-1*（G（ij）*e（i）+B（ij）*f（i））;
????????????????????Jac（2*i-12*j）=（G（ij）*e（i）+B（ij）*f（i））;
????????????????????Jac（2*i2*j）=B（ij）*e（i）-G（ij）*f（i）;
????????????????????Jac（2*i-12*j-1）=B（ij）*e（i）-G（ij）*f（i）;
????????????????else?if?j==i
????????????????????????m=0;h=0;
????????????????????????for?r=1:n
????????????????????????????m=m+G（ir）*e（r）-B（ir）*f（r）;
????????????????????????????h=h+G（ir）*f（r）+B（ir）*e（r）;
????????????????????????end
????????????????????????Jac（2*i2*j-1）=-1*m-G（ii）*e（i）-B（ii）*f（i）;
????????????????????????Jac（2*i-12*j）=-1*m+G（ii）*e（i）+B（ii）*f（i）;
????????????????????????Jac（2*i2*j）=-1*h+B（ii）*e（i）-G（ii）*f（i）;
????????????????????????Jac（2*i-12*j-1）=h+B（ii）*e（i）-G（ii）*f（i）;
????????????????????end
????????????????end
????????????else?if?b2（i6）==3%P，V節點；
????????????????????if?j~=i
????????????????????????Jac（2*i2*j-1）=-1*（G（ij）*e（i）+B（ij）*f（i））;
????????????????????????Jac（2*i-12*j）=0;
????????????????????????Jac（2*i2*j）=B（ij）*e（i）-G（ij）*f（i）;