function?reactivepower
clear
baseMVA=100;????????????????????????%功率基值
%%讀Data1中數(shù)據(jù)?
load?Data1.txt
Bus=Data1（:1）;????????????????????????%節(jié)點號
Vtype=Data1（:2）;????????????????????%節(jié)點類型
Pload=Data1（:3）;????????????????????%負(fù)載有功
Qload=Data1（:4）;????????????????????%負(fù)載無功
Vset=Data1（:8）;??????????????????????%母線電壓幅值
Pgen=Data1（:14）;??????????????????%發(fā)電機發(fā)出有功
Qgen=Data1（:15）;??????????????????%發(fā)電機發(fā)出無功
Qsh=Data1（:16）;??????????????????%并聯(lián)電容電納標(biāo)幺值
%%讀Data2中數(shù)據(jù)
load?Data2.txt
II=Data2（:1）;??????????????????????
JJ=Data2（:2）;????????????????????????%支路端點號
R=Data2（:3）;?????????????????????????%兩點間電阻
X=Data2（:4）;?????????????????????????%兩點間電抗
B=Data2（:5）/2;?????????????????????%線路對地電納
K=Data2（:9）;?????????????????????????%變壓器非標(biāo)準(zhǔn)變壓比
Ltype=Data2（:14）;????????????????%線路類型
%%求導(dǎo)納矩陣Y
y1=zeros（14）;
y2=zeros（14）;
y3=zeros（14）;
lin=length（II）;?????????????????????%支路數(shù)
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）;
%%設(shè)電壓初值
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????????????????????????????????????????????%節(jié)點功率及電壓不平衡量
????????switch?Vtype（x）
????????????case?1???????????????????????????????????????????%PQ節(jié)點
????????????????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節(jié)點
????????????????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