%一.讀取參數(shù)
clear??
format?short
f=load（‘F:\電力系統(tǒng)分析作業(yè)\123.txt‘）;%打開路徑下存放的IEEE14節(jié)點的參數(shù)數(shù)據(jù)
zljz=zeros（208）;
zljz=f（1:201:8）;
jdjz=zeros（148）;
jdjz=f（21:341:8）;
gnd=zeros（13）;?
gnd=f（351:3）;
%對地支路數(shù)
ngnd=f（361）;
%節(jié)點個數(shù)?
npoint=f（371）;??
%線路支路數(shù)?
nzl=f（381）;
error=1e-12;%容許誤差
Nmax=600;%最大迭代次數(shù)

%二.形成節(jié)點導(dǎo)納矩陣
??%PV節(jié)點個數(shù)
pvjd=0;
for?i=1:npoint??????????
???????if（jdjz（i2）==2）
????????pvjd=pvjd+1;
???????end
end
???%PQ節(jié)點個數(shù)?
?pqjd=npoint-pvjd-1;
?%矩陣維數(shù)
?ndim=npoint-1+pqjd;?
%找出平衡節(jié)點序號
?for?i=1:npoint????????????
???????if（jdjz（i2）==1）
???????????nphjd=i;
???????end
?end??
?%將平衡節(jié)點調(diào)整至最后一行
???jdjz1=zeros（18）;????
???jdjz1（1:）=jdjz（nphjd:）;
???jdjz（nphjd:）=jdjz（npoint:）;
???jdjz（npoint:）=jdjz1（1:）;
???for?k=1:nzl?????????????%更新線路
??????i=1;
??????for?i=1:2
??????????if（zljz（ki）==nphjd）
??????????zljz（ki）=npoint;
??????????elseif（zljz（ki）==npoint）
??????????????zljz（ki）=nphjd;
????????????end
????????end
???end
%初始化各導(dǎo)納矩陣，Y，YR（實部），YI（虛部）??
YR=zeros（npointnpoint）;
YI=zeros（npointnpoint）;
Y=zeros（npointnpoint）;
%1.不含變壓器的線路導(dǎo)納矩陣
k=1;
for?k=1:nzl????????????????
???i=zljz（k1）;
???j=zljz（k2）;
???mo=zljz（k3）^2+zljz（k4）^2;
????if（zljz（k7）==0）
???YR（ii）=YR（ii）+zljz（k3）/mo+0.5*zljz（k5）;
???YI（ii）=YI（ii）-zljz（k4）/mo+0.5*zljz（k6）;
???YR（ij）=YR（ij）-zljz（k3）/mo;
???YI（ij）=YI（ij）+zljz（k4）/mo;
???YR（ji）=YR（ji）-zljz（k3）/mo;
???YI（ji）=YI（ji）+zljz（k4）/mo;
???YR（jj）=YR（jj）+zljz（k3）/mo+0.5*zljz（k5）;
???YI（jj）=YI（jj）-zljz（k4）/mo+0.5*zljz（k6）;
???end
end
Y=YR+YI*1i;
%2.含變壓支路的線路導(dǎo)納矩陣
i=1;
for?i=1:nzl??
???if（zljz（i7）==1）
???????mo=zljz（i3）^2+zljz（i4）^2;
???????YR（zljz（i1）zljz（i1））=（zljz（i8）^2）*zljz（i3）/mo+YR（zljz（i1）zljz（i1））;
???????YI（zljz（i1）zljz（i1））=-（zljz（i8）^2）*zljz（i4）/mo+YI（zljz（i1）zljz（i1））;
???????YR（zljz（i2）zljz（i2））=zljz（i3）/mo+YR（zljz（i2）zljz（i2））;
???????YI（zljz（i2）zljz（i2））=-zljz（i4）/mo+YI（zljz（i2）zljz（i2））;
???????YR（zljz（i1）zljz（i2））=-zljz（i8）*zljz（i3）/mo+YR（zljz（i1）zljz（i2））;
???????YI（zljz（i1）zljz（i2））=-zljz（i8）*-zljz（i4）/mo+YI（zljz（i1）zljz（i2））;
???????YR（zljz（i2）zljz（i1））=-zljz（i8）*zljz（i3）/mo+YR（zljz（i2）zljz（i1））;
???????YI（zljz（i2）zljz（i1））=-zljz（i8）*-zljz（i4）/mo+YI（zljz（i2）zljz（i1））;
???end
end
??Y=YR+YI*1i;?
%3.對地并聯(lián)支路的導(dǎo)納矩陣修正
i=1;
??for?i=1:ngnd???????????????
???c=gnd（i1）;
????YR（cc）=YR（cc）+gnd（i2）;
????YI（cc）=YI（cc）-gnd（i3）;
end
???Y=YR+YI*1i;
disp（‘節(jié)點導(dǎo)納矩陣為：‘）?
Y
%三.B’B’’矩陣的生成
%B’矩陣生成
B=imag（Y）;?
ph=find（jdjz（:2）==1）;?
BB=B;?
BB（:ph）=[];?BB（ph:）=[];?
B1=BB;
%B’’矩陣的生成
phpv=find（jdjz（:2）~=3）;?
???????BB=B;??
?BB（:phpv）=[];?
?BB（phpv:）=[];??????
??B2=BB;??
%四.計算不平衡量dPdQ
f=0;%迭代標(biāo)志位
N=1;%迭代次數(shù)初始化
%計算各節(jié)點凈有功
Pjd=zeros（1npoint）;
i=1;
for?i=1:npoint
????Pjd（i）=jdjz（i3）-jdjz（i5）;
????end
%節(jié)點凈無功功率
Qjd=zeros（1npoint）;
i=1;
for?i=1:npoint
????Qjd（i）=jdjz（i4）-jdjz（i6）;
????end
%計算各支路功率有功和無功:統(tǒng)一潮流方程計算每個節(jié)點的功率大小
while（f==