%高斯——賽爾德潮流計算
%[s_nodeu_node]=gaosai（YmnsS_PQP_PVU_PVusQmaxQminexpl）
%返回值：s_node:??節點功率
%???????u_node：??節點電壓
%輸入值：Y：???導納矩陣
%????????m：???PQ節點數包含平衡節點
%????????n：???總結點數
%????????s：???平衡節點編號，一般定為1號編號
%????????S_PQ:?PQ節點的視在功率[1:m]其中平衡節點的視在功率暫定為0；
%????????P_PV：PV節點的有功功率[m+1:n]
%????????U_PV:?PV節點的電壓給定值只有幅值，沒有相位角，是實數[m+1:n]
%????????us:???平衡節點的電壓給定值，復數。
%????????Qmax；PV節點無功功率最大值
%????????Qmin：PV節點無功功率最小值
%????????epxl:?最大電壓變量dumax的最大值，即規定的潮流計算的精度，實數。
function?[u_node]=DYgaosai（YmnsS_PQP_PVU_PVusQmaxQminepxl）
%**************************************************************************
%以下為高斯塞爾德潮流計算程序
?fprintf（‘高斯——塞爾德潮流計算開始\n‘）
U=ones（1n）;%節點電壓的初始值；
U（s）=us;%將電壓初始值的平衡節點的初始電壓設為給定值
k=1000;%設置迭代次數
dumax=0;%最大電壓變化量，實數
flag=zeros（1n-m）;%PV節點在計算過程中轉化為PQ節點的標志，PV轉化為PQ節點時，相應的標志位置1
Q_PV=ones（1n-m）;%PV節點的無功功率；
Q_PV=Q_PV.*0.2;
Ss=0;%平衡節點的視在功率
for?a=m+1:n
????U（a）=U_PV（a-m）;%第a個節點是PV節點，暫時將PV節點電壓給定值賦給節點電壓初始值其相角暫時為0
end
while（k>=1）
for?i=1:n
????u=U（i）;%保存U（i）的初值，為求du和dumax做準備；
????if（i==s）%第i個節點是否為平衡節點？
???????%?i=i+1;%d第i個節點是平衡節點，跳過。
????else????%第i節點不是平衡節點，進行以下處理
????????A=0;
????????for?a=1:n
????????????if（a~=i）
?????????????A=A+Y（ia）*U（a）;
????????????end
????????end
????????a=1;
????????if（i>=m+1?&&?flag（i-m）==0）?%第i各節點是否為PV節點？
????????????Q_PV（i-m）=-imag（U（i）‘*Y（ii）*U（i）+U（i）‘*A）;
????????????if（Q_PV（i-m）<=Qmax?&&?Q_PV（i-m）>=Qmin）
????????????????Upv=（（P_PV（i-m）-Q_PV（i-m）*j）/U（i）-A）/Y（ii）;
????????????????ang=angle（Upv）;
????????????????U（i）=abs（U（i））*（cos（ang）+sin（ang）*j）;
%???????????????Q_PV（i-m）=-imag（U（i）‘*Y（ii）*U（i）+U（i）‘*A）;
????????????else
????????????????if（Q_PV（i-m）>Qmax）
????????????????????Q_PV（i-m）=Qmax;
????????????????????flag（i-m）=1;
????????????????????fprintf（‘PV節點%gQ大于Qmax轉化為PQ節點\n‘i）;
????????????????elseif（Q_PV（i-m）????????????????????Q_PV（i-m）=Qmin;
????????????????????flag（i-m）=1;
????????????????????fprintf（‘PV節點%gQ小于Qmin轉化為PQ節點\n‘i）;
????????????????end
????????????????U（i）=（（P_PV（i-m）-Q_PV（i-m）*j）/U（i）‘-A）/Y（ii）;
????????????end
????????else
????????????if（i<=m）%第i個節點是PQ節點
????????????????U（i）=（S_PQ（i）‘/U（i）‘-A）/Y（ii）;
????????????elseif（i>=m+1?&&?flag（i-m）==1）
????????????????U（i）=（（P_PV（i-m）-Q_PV（i-m）*j）/U（i）‘-A）/Y（ii）;
????????????end
????????end
????end
?du=abs（U（i）-u）;
??if（du>dumax）
??dumax=du;
??end
end
i=1;
if（dumax????fprintf（‘迭代次數%g\n‘1000-k+1）
????break;
else
????k=k-1;
????dumax=0;
end
end
if（k==0）
????fprintf（‘高斯——塞爾德潮流計算不能收斂！！！\n‘）
end
for?i=1:n?%計算平衡節點的視在功率
????????Ss=Ss+Y（si）‘*U（i）‘;
end
Ss=Ss*U（s）;%平衡節點的視在功率
????for?i=1:n%計算各個節點的視在功率
????????if（i==s）
????????????s_node（i）=Ss;
????????else
????????????if（i<=m）
????????????????s_node（i）=S_PQ（i）;
????????????elseif（i>=m+1）
????????????????s_node（i）=P_PV（i-m）+Q_PV（i-m）*j;
????????????end
????????end
????end
????s_node
????u_node=U
????fprintf（‘高斯——塞爾德潮流計算結束\n‘）

????
????????????
????????????

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

?????文件???????3531??2011-05-07?20:29??newton\DYgaosai.m

?????文件??????21808??2011-06-06?09:25??newton\newton.m

?????文件???????2597??2011-05-12?08:59??newton\write.m

?????文件??????16384??2011-06-06?09:25??newton\輸出數據.xls

?????目錄??????????0??2011-10-18?08:58??newton

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

????????????????44320????????????????????5