【正文】 amp。%PQ節(jié)點的Jia參數(shù)=Nia2real(I(i)) Ykb(2*(i1),2*(a1))=Ykb(2*(i1)1,2*(a1)1)+2*imag(I(i))。%PQ節(jié)點Hia參數(shù) Ykb(2*(i1)1,2*(a1))=real(Y(i,a))*real(U(i))imag(Y(i,a))*imag(U(i))real(I(i))。amp。%PQ節(jié)點的Jia參數(shù) Ykb(2*(i1),2*a)=Ykb(2*(i1)1,2*a1)+2*imag(I(i))。%PQ節(jié)點Hia參數(shù) Ykb(2*(i1)1,2*a)=real(Y(i,a))*real(U(i))imag(Y(i,a))*imag(U(i))real(I(i))。Ybk中的ii1。%PQ節(jié)點的Jia參數(shù)=Nia Ykb(2*(i1)1,2*(a1))=Ykb(2*(i1)1,2*(a1)1)。%PQ節(jié)點Hia參數(shù) Ykb(2*(i1)1,2*(a1))=real(Y(i,a))*real(U(i))imag(Y(i,a))*imag(U(i))。amp。%PQ節(jié)點的Jia參數(shù)=Nia Ykb(2*(i1),2*(a1))=Ykb(2*(i1)1,2*(a1)1)。%PQ節(jié)點Hia參數(shù) Ykb(2*(i1)1,2*(a1))=real(Y(i,a))*real(U(i))imag(Y(i,a))*imag(U(i))。aa1。 a=m)%如果a號節(jié)點是平衡節(jié)點之后的PQ節(jié)點。%PQ節(jié)點的Lia參數(shù)Hia else if(a=s+1 amp。%PQ節(jié)點Nia參數(shù) Ykb(2*(i1),2*a1)=Ykb(2*(i1)1,2*a)。a不變 Ykb(2*(i1)1,2*a1)=imag(Y(i,a))*real(U(i))real(Y(i,a))*imag(U(i))。 for a=1:n if(a~=i) if(a=s1)%如果a號節(jié)點是平衡節(jié)點之前的PQ節(jié)點。 i=m)%如果i號節(jié)點是平衡節(jié)點之后的PQ節(jié)點。 else if(i=s+1 amp。%PQ節(jié)點的Lia參數(shù)=Hia+2imag(I(i))。%PQ節(jié)點Nia參數(shù) Ykb(2*i,2*(a1)1)=Ykb(2*i1,2*(a1))2*real(I(i))。 a=n) Ykb(2*i1,2*(a1)1)=imag(Y(i,a))*real(U(i))real(Y(i,a))*imag(U(i))imag(I(i))。 elseif(a=m+1 amp。%PQ節(jié)點的Jia參 數(shù)=Nia2real（I(i)）； Ykb(2*i,2*(a1))=Ykb(2*i1,2*(a1)1)+2*imag(I(i))。%PQ節(jié)點Hia參數(shù) Ykb(2*i1,2*(a1))=real(Y(i,a))*real(U(i))imag(Y(i,a))*imag(U(i))real(I(i))。amp。%PQ節(jié)點的Jia參數(shù)Ykb(2*i,2*a)=Ykb(2*i1,2*a1)+2*imag(I(i))。%PQ節(jié)點Hia參數(shù) Ykb(2*i1,2*a)=real(Y(i,a))*real(U(i))imag(Y(i,a))*imag(U(i))real(I(i))。%PV節(jié)點Jia參數(shù)=NiaYkb(2*i1,2*(a1))=Ykb(2*i1,2*(a1)1)。%PQ節(jié)點Hia參數(shù) Ykb(2*i1,2*(a1))=real(Y(i,a))*real(U(i))imag(Y(i,a))*imag(U(i))。 a=n)%Ykb中的aa1。%PQ節(jié)點的Lia參數(shù)=Hia elseif(a=m+1 amp。%PQ節(jié)點Nia參數(shù) Ykb(2*i,2*(a1)1)=Ykb(2*i1,2*(a1))。i不變 Ykb(2*i1,2*(a1)1)=imag(Y(i,a))*real(U(i))real(Y(i,a))*imag(U(i))。amp。%PQ節(jié)點的Jia參數(shù)=NiaYkb(2*i,2*a)=Ykb(2*i1,2*a1)。 %PQ節(jié)點Hia參數(shù) Ykb(2*i1,2*a)=real(Y(i,a))*real(U(i))imag(Y(i,a))*imag(U(i))。 end end i=1。 elseif(i=s+1) for a=1:n I(i)=I(i)+Y(i,a)*U(a)。 %節(jié)點的不平衡量的dPQV2計算結(jié)束 for i=1:n%計算節(jié)點注入電流I； if(i=s1) for a=1:n I(i)=I(i)+Y(i,a)*U(a)。%PV節(jié)點有功功率不平衡量 dPQV2(2*(i1))=U_PV(im)^2abs(U(i))^2。amp。%PQ節(jié)點有功功率不平衡量 dPQV2(2*(i1))=imag(S_PQ(i))imag(A)。amp。%PQ節(jié)點有功功率不平衡量 dPQV2(2*i)=imag(S_PQ(i))imag(A)。 A=U(i)*A。)。 for a=1:n A=A+(Y(i,a)39。flag=zeros(1,nm)。%PV節(jié)點的無功功率；Q_PV=Q_PV.*。%平衡節(jié)點的視在功率S=zeros(1,n)。k=1000。dU=zeros(2*(n1),1)。%不平衡量向量；Ykb=zeros(2*(n1),2*(n1))。牛頓——拉夫遜潮流計算開始\n39。%利用高斯塞爾德法得到比較準確的節(jié)點電壓的初值。% Qmax；PV節(jié)點無功功率最大值% Qmin：PV節(jié)點無功功率最小值% epxl: 最大電壓變量dumax的最大值，即規(guī)定的潮流計算的精度，實數(shù)。 end end end s_node u_node=U2 牛頓——拉夫遜潮流計算程序%牛頓——拉夫遜潮流計算。 else if(i=m) s_node(i)=S_PQ(i)。endSs=Ss*U(s)。)endfor i=1:n %計算平衡節(jié)點的視在功率 Ss=Ss+Y(s,i)39。endendif(k==0) fprintf(39。else k=k1。迭代次數(shù)%g\n39。 endendi=1。 end end end du=abs(U(i)u)。 flag(im)==1) U(i)=((P_PV(im)Q_PV(im)*j)/U(i)39。 elseif(i=m+1 amp。/U(i)39。A)/Y(i,i)。,i)。 fprintf(39。 elseif(Q_PV(im)Qmin) Q_PV(im)=Qmin。PV節(jié)點%gQ大于Qmax轉(zhuǎn)化為PQ節(jié)點\n39。 flag(im)=1。 U(i)=abs(U(i))*(cos(ang)+sin(ang)*j)。 Q_PV(im)=Qmin) Upv=((P_PV(im)Q_PV(im)*j)/U(i)A)/Y(i,i)。 if(Q_PV(im)=Qmax amp。*Y(i,i)*U(i)+U(i)39。amp。 end end a=1。 else %第i節(jié)點不是平衡節(jié)點，進行以下處理 A=0。%保存U(i)的初值，為求du和dumax做準備； if(i==s)%第i個節(jié)點是否為平衡節(jié)點？ % i=i+1。%平衡節(jié)點的視在功率for a=m+1:n U(a)=U_PV(am)。%PV節(jié)點的無功功率；Q_PV=Q_PV.*。%最大電壓變化量，實數(shù)flag=zeros(1,nm)。%將電壓初始值的平衡節(jié)點的初始電壓設為給定值k=1000。)U=ones(1,n)。function [u_node]=DYgaosai(Y,m,n,s,S_PQ,P_PV,U_PV,us,Qmax,Qmin,epxl)%**************************************************************************%以下為高斯塞爾德潮流計算程序fprintf(39。參考文獻參考文獻[1] 陳珩. 電力系統(tǒng)穩(wěn)態(tài)分析 [M]. 南京：南京工學院，1985051.[2] 王錫凡. 現(xiàn)代電力系統(tǒng)分析 [M]. 北京：科學出版社，20033.[3] 西安交通大學. 電子數(shù)字計算機的應用：電力系統(tǒng)計算 [M]. 西安：西安交通大學，1987101.[4] 南京工學院. 電力系統(tǒng) [M]. 南京：南京工學院，1980081.[5] 吳天明，彭彬. MATLAB電力系統(tǒng)設計與分析 [M]. 國防工業(yè)出版社，200401.[6] 林暉. VSCMTDC和交流混合的電力系統(tǒng)潮流計算 [D]. 湖南大學電氣與信息工程學院，.[7] 梁海平. 智能潮流分析系統(tǒng)的研究與開發(fā) [D]. 華北電力大學電力系統(tǒng)及其自動化專業(yè)，.[8] 許建明. 經(jīng)典法和遺傳算法在電力系統(tǒng)有功優(yōu)化中的研究 [D]. 南昌大學信息工程學院，.[9] 趙芳 陳士方 張圣集. 保留非線性的電力系統(tǒng)概率潮流計算 [J]. 電力情報，2000年03期.[10] 成海濱 沈茂亞. 電力系統(tǒng)最優(yōu)潮流算法研究綜述 [J]. 電力電氣，2006年第25卷第11期.[11] 羅瑋. 電力市場環(huán)境下基于最優(yōu)潮流的輸電容量充裕度研究 [D]. 上海交通大學電力系統(tǒng)及其自動化專業(yè)，[12] 郭金偉. 電力市場環(huán)境下基于最優(yōu)潮流的節(jié)點實時電價和購電份額研究 [D]. 北京交通大學，2008.[13] 劉方. 電力系統(tǒng)動態(tài)最優(yōu)潮流的模型與算法研究[D]. 重慶大學電氣工程學院，[14] 郗忠梅 李有安 趙法起 張博. 基于Matlab的電力系統(tǒng)潮流計算[J]. 山東農(nóng)業(yè)大學學報，2010年41卷第2期.[15] 樊宇璐 李世作 張志斌. 基于擬牛頓法的電力系統(tǒng)潮流計算[J]. 電氣開關，2009年第2期.[16] Yan Sun, NaiShan Flow Calculation Method for Islanded Power Network [J]. Department of Electrical Engineering Guangxi University ,Nanning, China,.[17] Andrey Pazderin, Sergey Flow Calculation by Combination ofNewtonRaphson Method and Newton’s Method in Optimization.[J].URAL STATE TECHNICAL UNIVERSITY – UPI,[18]D. A. Arzamastsev, P. I. Bartolomej, A. M. Holjan, Automatic control system and optimization of modes of electric power systems, Moscow:Vysshaya Shkola, 1983.[19]V. I. Tarasov, Nonlinear minimization methods for