function?reactivepower
clear
baseMVA=100;????????????????????????%功率基值
%%讀Data1中數據?
load?Data1.txt
Bus=Data1（:1）;????????????????????????%節點號
Vtype=Data1（:2）;????????????????????%節點類型
Pload=Data1（:3）;????????????????????%負載有功
Qload=Data1（:4）;????????????????????%負載無功
Vset=Data1（:8）;??????????????????????%母線電壓幅值
Pgen=Data1（:14）;??????????????????%發電機發出有功
Qgen=Data1（:15）;??????????????????%發電機發出無功
Qsh=Data1（:16）;??????????????????%并聯電容電納標幺值
%%讀Data2中數據
load?Data2.txt
II=Data2（:1）;??????????????????????
JJ=Data2（:2）;????????????????????????%支路端點號
R=Data2（:3）;?????????????????????????%兩點間電阻
X=Data2（:4）;?????????????????????????%兩點間電抗
B=Data2（:5）/2;?????????????????????%線路對地電納
K=Data2（:9）;?????????????????????????%變壓器非標準變壓比
Ltype=Data2（:14）;????????????????%線路類型
%%求導納矩陣Y
y1=zeros（14）;
y2=zeros（14）;
y3=zeros（14）;
lin=length（II）;?????????????????????%支路數
for?x=1:lin
????switch?Ltype（x）
????????case?1
????????????y1（II（x）JJ（x））=1/（R（x）+i*X（x））;
????????????y1（JJ（x）II（x））=y1（II（x）JJ（x））;
????????????y3（II（x）JJ（x））=i*B（x）;
????????????y3（JJ（x）II（x））=i*B（x）;
????????case?2
????????????y1（II（x）JJ（x））=1/（（R（x）+i*X（x））*K（x））;
????????????y1（JJ（x）II（x））=y1（II（x）JJ（x））;
????????????y2（II（x）JJ（x））=（1-K（x））/（（R（x）+i*X（x））*K（x）^2）;
????????????y2（JJ（x）II（x））=（K（x）-1）/（（R（x）+i*X（x））*K（x））;
????end??
end
clear?x
Y=zeros（14）;
for?x=1:14
????Y（xx）=sum（y1（x:））+sum（y2（x:））+sum（y3（x:））+i*Qsh（x）;
end
clear?x;
Y=Y-y1;
G=real（Y）;
B=imag（Y）;
%%設電壓初值
U=Vset;
e=real（U）;
f=imag（U）;
%%
Ps=zeros（114）;
Qs=zeros（114）;
D=ones（261）;
for?x=1:14
????Ps（x）=（Pgen（x）-Pload（x））/baseMVA;
????Qs（x）=（Qgen（x）-Qload（x））/baseMVA;
end
clear?x;
N=0;
Jacbi=zeros（26）;
while?max（abs（D））>0.000001
???????for?x=2:14????????????????????????????????????????????%節點功率及電壓不平衡量
????????switch?Vtype（x）
????????????case?1???????????????????????????????????????????%PQ節點
????????????????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???????????????????????????????????????????%PV節點
????????????????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:14???????????????????????????????????????????????????%求雅克比矩陣
??????for?J=2:14
??????????if?I~=J??????????????????????????????????????????????%非對角元素
??????????????Jacbi（（2*I-3）（2*J-3））=B