clear?all?
clc
N=16;%輸入系統(tǒng)節(jié)點(diǎn)數(shù)
%輸入系統(tǒng)數(shù)據(jù)（支路、節(jié)點(diǎn)）矩陣A
%第一列是支路號
%第二列是支路開關(guān)狀態(tài)；用K=1或0表示分段開關(guān)閉合或聯(lián)絡(luò)開關(guān)斷開
%第三列支路首端節(jié)點(diǎn)號
%第四列支路末端節(jié)點(diǎn)號
%第五列支路電阻（標(biāo)幺值）
%第六列支路電抗（標(biāo)幺值）
%第七列支路末端節(jié)點(diǎn)負(fù)荷有功（標(biāo)幺值）
%第八列支路末端節(jié)點(diǎn)負(fù)荷無功（標(biāo)幺值）
%第九列支路末端節(jié)點(diǎn)連接的電容器容量（標(biāo)幺值）
A=[?1???1???1???2???0.075???0.10???0.02???0.016?????0
????2???1???1???3???0.11????0.11???0.04???0.027?????0
????3???1???1???4???0.11????0.11???0.01???0.009?????0
????4???1???2???5???0.09????0.18???0.02???0.008?????0.0025
????5???1???2???6???0.08????0.11???0.03???0.015?????0.002
????6???1???5???7???0.04????0.04???0.015??0.012?????0
????7???1???3???8???0.08????0.11???0.05???0.03??????0.012
????8???1???8???9???0.08????0.11???0.045??0.02??????0.037
????9???1???8??10???0.11????0.11???0.006??0.001?????0.006
???10???1???3??11???0.11????0.11???0.01???0.009?????0.018
???11???1???4??12???0.09????0.12???0.01???0.007?????0
???12???1???4??13???0.08????0.11???0.01???0.009?????0?
???13???1??13??14???0.04????0.04???0.021??0.01??????0.018
???14???0???6??10???0.04????0.04???0.006??0.001?????0.006
???15???0??11??12???0.04????0.04???0.001??0.007?????0
???16???0???7??14???0.09????0.12???0.021??0.01??????0.018?????];
[mAnA]=size（A）;
a=0;b=0;
for?i=1:mA
???K=A（i2）;
????if?K==1
????????a=a+1;
????????C（a1）=A（i1）;
????else?if?K==0
????????????b=b+1;
????????????D（b1）=A（i1）;
????????end?
????end
end??????
C%C為閉合開關(guān)投入的支路
D%D為聯(lián)絡(luò)開關(guān)打開的支路
[mCnC]=size（C）
B=zeros（mCnA）;
for?i=1:mC
????B（i:）=A（i:）;
end
B
Vbus=ones（N1）;%給定電壓初值
Vbus（11）=1.0;
Vbus2=Vbus;%迭代后的電壓矩陣
[mBnB]=size（B）
k=0;%迭代次數(shù)
Ploss=zeros（mB1）;%存支路的有功損耗
Qloss=zeros（mB1）;%支路無功損耗
I=zeros（mB1）;%支路電流
F=zeros（N1）;
%%本段程序?qū)⒅分匦屡判???s1為排好序的支路矩陣
TempB=B;
n=1;
s2=[];????
while?~isempty（TempB）
????[snB]=size（TempB）;%s為支路數(shù)??
????m=1;
????while?s>0
???????i=find（TempB（:3）==TempB（s4））;%末端節(jié)點(diǎn)是否為其他支路首端節(jié)點(diǎn)
????????if?isempty（i）
????????????s1（n:）=TempB（s:）;%如果i是空集則該節(jié)點(diǎn)為葉節(jié)點(diǎn)
????????????n=n+1;
????????else?s2（m:）=TempB（s:）;%如果i不是空集則該節(jié)點(diǎn)為非葉節(jié)點(diǎn)
????????????m=m+1;
????????end
????????s=s-1;
???end
????TempB=s2;?????????????%將s2賦值給TempB.重復(fù)上述判斷，直到TempB為空集為止
????s2=[];
end
s1
?%前推進(jìn)行支路功率計(jì)算
?????while?k<10??%（一般取3~5次迭代即可收斂）
????Pij1=zeros（N1）;?%該支路首端節(jié)點(diǎn)及與其相連的其他后續(xù)支路的功率情況
????Qij1=zeros（N1）;
??????for?s=1:mB????????
????????ii=s1（s3）;??????%取排好序的支路首節(jié)點(diǎn)???
????????jj=s1（s4）;??????%取排好序的支路尾節(jié)點(diǎn)
????????Pload=s1（s7）;
????????Qload=s1（s8）-s1（s9）;?%節(jié)點(diǎn)有功和無功負(fù)荷
????????R