clear
%手動導入節點數據后保存為.mat文件
%所有數據均為標幺值
load（‘matlab.mat‘）
%Bus：節點號
%Vtype：節點類型
%Pload：負荷有功
%Qload：負荷無功
%Pgen：發電機有功
%Qgen：發電機無功
%Vset：節點電壓
%Qsh：電納標幺值
%II、JJ：支路端點號
%R：支路電阻
%X：支路電抗
%B：對地電納/2
%K：變壓器非標準變壓比
%Ltype：是否有變比支路
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%求導納矩陣Y
y1=zeros（30）;
y2=zeros（30）;
y3=zeros（30）;
lin=length（II）;%求支路數
for?x=1:lin
????switch?Ltype（x）
????????case?1
????????????y1（II（x）JJ（x））=1/（R（x）+1i*X（x））;
????????????y1（JJ（x）II（x））=y1（II（x）JJ（x））;
????????????y3（II（x）JJ（x））=1i*B（x）;
????????????y3（JJ（x）II（x））=1i*B（x）;
????????case?2
????????????y1（II（x）JJ（x））=1/（（R（x）+1i*X（x））*K（x））;
????????????y1（JJ（x）II（x））=y1（II（x）JJ（x））;
????????????y2（II（x）JJ（x））=（1-K（x））/（（R（x）+1i*X（x））*K（x）^2）;%IIII=y1+y2???????????
????????????y2（JJ（x）II（x））=（K（x）-1）/（（R（x）+1i*X（x））*K（x））;%JJJJ=y1+y2
????end
end
clear?x;
Y=zeros（30）;%定義導納矩陣
for?x=1:30%求對角線元素
????Y（xx）=sum（y1（x:））+sum（y2（x:））+sum（y3（x:））+1i*Qsh（x）;
end
clear?x;
Y=Y-y1;%非對角線元素
G=real（Y）;
B=imag（Y）;
%%%%%%%%%%%%%%%%%%
%設電壓初值
U=Vset;
e=real（U）;%實部
f=imag（U）;
%
Ps=zeros（301）;%注入有功
Qs=zeros（301）;%注入無功
D=ones（301）;
for?x=1:30
????Ps（x）=Pgen（x）-Pload（x）;
????Qs（x）=Qgen（x）-Qload（x）;
end
clear?x;
N=0;
%24*2+5+2=55
Jacbi=zeros（55）;
while?max（abs（D））>0.0000000001
for?x=2:30%節點功率及電壓不平衡量
????switch?Vtype（x）
????????case?1
????????????D（2*x-3）=Ps（x）-e（x）*（G（x:）*e-B（x:）*f）-f（x）*（G（x:）*f+B（x:）*e）;
????????????D（2*x-2）=Qs（x）-f（x）*（G（x:）*e-B（x:）*f）+e（x）*（G（x:）*f+B（x:）*e）;
????????case?2
????????????D（2*x-3）=Ps（x）-e（x）*（G（x:）*e-B（x:）*f）-f（x）*（G（x:）*f+B（x:）*e）;
????????????D（2*x-2）=Vset（x）.*Vset（x）-（e（x）.^2+f（x）.^2）;
????end
end
clear?m
for?I=2:30
????for?J=2:30
????????if?I~=J%非對角元素
????????????Jacbi（（2*I-3）（2*J-3））=B（IJ）*e（I）-G（IJ）*f（I）;
????????????Jacbi（（2*I-3）（2*J-2））=-（G（IJ）*e（I）+B（IJ）*f（I））;
??????????s