clear
clc
%極坐標(biāo)形式的牛頓法潮流計算程序
%B1矩陣為支路參數(shù)矩陣，其中1、2：分別為支路首端、末端節(jié)點編號；3：支路阻抗；4：接地導(dǎo)納；5：支路變比；6：支路首端處于K側(cè)為11側(cè)為0；7：單側(cè)對地導(dǎo)納
%B2矩陣為節(jié)點參數(shù)矩陣，其中1：該節(jié)點的發(fā)電機功率；2：該節(jié)點負(fù)荷功率；3：節(jié)點電壓初始幅值；4：節(jié)點電壓初始相位；5：PV節(jié)點的電壓給定值；6：節(jié)點所接的無功補償設(shè)備容量；
%7：節(jié)點分類編號：1為平衡節(jié)點，2為PQ節(jié)點，3為PV節(jié)點
n=input（‘請輸入節(jié)點數(shù):n=‘）;
nl=input（‘請輸入支路數(shù):nl=‘）;
isb=input（‘請輸入平衡節(jié)點編號:isb=‘）;
pr=input（‘請輸入誤差精度:pr=‘）;
B1=input（‘請輸入支路參數(shù)矩陣:B1=‘）;
B2=input（‘請輸入節(jié)點參數(shù)矩陣:B2=‘）;
%創(chuàng)建矩陣
Y=zeros（n）;
V=zeros（n1）;
for?i=1:nl
????if?B1（i6）==0?%支路不含變壓器
????????p=B1（i1）;q=B1（i2）;
????????Y（pq）=Y（pq）-1./B1（i3）;%導(dǎo)納矩陣的非對角元
????????Y（qp）=Y（pq）;
????????Y（pp）=Y（pp）+1./B1（i3）+0.5*B1（i4）+B1（i7）;%導(dǎo)納矩陣的對角元
????????Y（qq）=Y（qq）+1./B1（i3）+0.5*B1（i4）+B1（i7）;
????else?
????????p=B1（i2）;q=B1（i1）;%支路含有變壓器
????????Y（pq）=Y（pq）-1./（B1（i3）*B1（i5））;
????????Y（qp）=Y（pq）;
????????Y（pp）=Y（pp）+1./B1（i3）+B1（i7）;?%pi型等效電路
????????Y（qq）=Y（qq）+1./（（B1（i5）^2*B1（i3）））+B1（i7）;
????end
end
%分解導(dǎo)納的實部虛部
G=real（Y）;B=imag（Y）;
for?i=1:n
????theta（i1）=B2（i4）;
????V（i1）=B2（i3）;
end
%求出兩節(jié)點的相位差
Times=0;
%求出PQ節(jié)點的個數(shù)
m=0;
for?i=1:n
????if?B2（i7）==2
????????m=m+1;
????end
end
M=ones（n-1+m1）;
while?max（abs（M））>?pr
%形成潮流方程
P=zeros（n1）;
Q=zeros（n1）;
for?i=1:n
????for?j=1:n
????????P（i1）=V（i）*sum（（G（i:）.*cos（theta（i）-theta（:））‘+B（i:）.*sin（theta（i）-theta（:））‘）*V（:））;
????????Q（i1）=V（i）*sum（（G（i:）.*sin（theta（i）-theta（:））‘-B（i:）.*cos（theta（i）-theta（:））‘）*V（:））;
????end
end
%求出功率的不平衡量
dP=zero