clear
clc
%極坐標形式的牛頓法潮流計算程序
%B1矩陣為支路參數矩陣，其中1、2：分別為支路首端、末端節點編號；3：支路阻抗；4：接地導納；5：支路變比；6：支路首端處于K側為11側為0；7：單側對地導納
%B2矩陣為節點參數矩陣，其中1：該節點的發電機功率；2：該節點負荷功率；3：節點電壓初始幅值；4：節點電壓初始相位；5：PV節點的電壓給定值；6：節點所接的無功補償設備容量；
%7：節點分類編號：1為平衡節點，2為PQ節點，3為PV節點
n=input（‘請輸入節點數:n=‘）;
nl=input（‘請輸入支路數:nl=‘）;
isb=input（‘請輸入平衡節點編號:isb=‘）;
pr=input（‘請輸入誤差精度:pr=‘）;
B1=input（‘請輸入支路參數矩陣:B1=‘）;
B2=input（‘請輸入節點參數矩陣:B2=‘）;
%創建矩陣
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）;%導納矩陣的非對角元
????????Y（qp）=Y（pq）;
????????Y（pp）=Y（pp）+1./B1（i3）+0.5*B1（i4）+B1（i7）;%導納矩陣的對角元
????????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
%分解導納的實部虛部
G=real（Y）;B=imag（Y）;
for?i=1:n
????theta（i1）=B2（i4）;
????V（i1）=B2（i3）;
end
%求出兩節點的相位差
Times=0;
%求出PQ節點的個數
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=zeros（n-11）;
dQ=zeros（m1）;
for?i=1:n-1
????dP（i1）=real（B2（i1）-B2（i2））-P（i）;
end
for?i=1:m
????dQ（i1）=imag（B2（i1）-B2（i2））-Q（i）;
end
H=zeros（n-1n-1）;
N=zeros（n-1m）;
K=zeros（mn-1）;
L=zeros（mm）;
%求出雅克比矩陣元素的表達式
for?i=1:n-1
????for?j=1:n-1
????????if?i~=j
????????????H（ij）=-V（i）*V（j）*（G（ij）*sin（theta（i）-theta（j））-B（ij）*cos（theta（i）-theta（j）））;?
????????else
????????????H（ii）=V（i）^2*B（ii）+Q（i）;
????????end
????end
end
for?i=1:n-1
???for?j=1:m
???????if?i~=j
???????????N（ij）=-V（i）*V（j）*（G（ij）*cos（theta（i）-theta（j））+B（ij）*sin（theta（i）-theta（j）））;
???????else
???????????N（ii）=-V（i）^2*G（ii）-P（i）;?
???????end
???end
end
for?i=1:m
????for?j=1:n-1
????????if?i~=j
????????????K（ij）=V（i）*V（j）*（G（ij）*cos（theta（i）-theta（j））+B（ij）*sin（theta（i）-theta（j）））;
????????else?
????????????K（ii）=V（i）^2*G（ii）-P（i）;
????????end
????end
end
for?i=1:m
????for?j=1:m
????????if?i~=j
????????????L（ij）=-V（i）*V（j）*（G（ij）*sin（theta（i）-theta（j））-B（ij）*cos（theta（i）-theta（j）））;
????????else
????????????L（ii）=V（i）^2*B（ii）-Q（i）;
????????end
????end
end
Jacbi=-[HN;KL];
Vd=diag（V）;
A=[dP;dQ];
M=Jacbi\A;
for?i=1:n-1
????dtheta（i1）=M（i）;
end
for?i=1:n-1
????theta（i1）=theta（i）+dtheta（i）;
end
for?i=1:m
????dV（i1）=V（i）*M（i+n-1）;
end
for?i=1:m
????V（i）=V（i）+dV（i）;
end
for?i=1:n-1
????w（i1）=theta（i）*180/pi;
end
Times=Times+1;
end
Y
dP
dQ
H
N
K
L
Jacbi
Vd
dtheta
dV
Times
w
V

?屬性????????????大小?????日期????時間???名稱
-----------?---------??----------?-----??----

?????文件???????3239??2015-04-18?19:29??PolarCoordinates?-?副本.m

-----------?---------??----------?-----??----

?????????????????3239????????????????????1