1、電力系統(tǒng)穩(wěn)態(tài)分析課程設(shè)計題目名稱兩機五節(jié)點網(wǎng)絡(luò)潮流計算方法牛拉法和pq法學生姓名學 號系、專 業(yè)電氣工程系電氣工程及其自動化指導(dǎo)教師2012年12月28日 TOC o 1-5 h z 摘要第一章原理簡介31.1對潮流分析的簡介31.1.1潮流計算方法分析比較31.2 MATLAB 簡介41.2.1矩陣的運算51.3牛頓拉夫遜法計算潮流分布6第二章程序及結(jié)果102.1設(shè)計資料及參數(shù)102.1.1牛頓拉夫遜法的程序框圖132.2用Matlab設(shè)計程序142.2.1程序的編寫142.2.2程序運行結(jié)果19 HYPERLINK l bookmark112 o Current Document 2.2.

2、3p_q法程序編寫22 HYPERLINK l bookmark148 o Current Document 總結(jié)32參考文獻32電力系統(tǒng)穩(wěn)態(tài)分析課程設(shè)計1.1對潮流分析的簡介潮流分析是研究電力系統(tǒng)的一種最基本和最重要的計算。最初，電力 系統(tǒng)潮流計算是通過人工手算的，后來為了適應(yīng)電力系統(tǒng)日益發(fā)展的需要，采 用了錯誤!未指定書簽。交流計算臺。隨著電子數(shù)字計算機的出現(xiàn)，1956年Ward 等人編制了實際可行的計算機潮流計算程序。這樣，就為日趨復(fù)雜的大規(guī)模電 力系統(tǒng)提供了極其有力的計算手段。經(jīng)過幾十年的時間，電力系統(tǒng)潮流計算已 經(jīng)發(fā)展得十分成熟。潮流計算是研究電力系統(tǒng)穩(wěn)態(tài)運行情況的一種計算，是根 據(jù)

3、給定的運行條件及系統(tǒng)接線情況確定整個電力系統(tǒng)各個部分的運行狀態(tài)，如 各母線的電壓、各元件中流過的功率、系統(tǒng)的功率損耗等等。電力系統(tǒng)潮流計 算是計算系統(tǒng)動態(tài)穩(wěn)定和靜態(tài)穩(wěn)定的基礎(chǔ)。在電力系統(tǒng)規(guī)劃設(shè)計和現(xiàn)有電力系 統(tǒng)運行方式的研究中，都需要利用電力系統(tǒng)潮流計算來定量的比較供電方案或 運行方式的合理性、可靠性和經(jīng)濟性。1.1.1潮流計算方法分析比較高斯-賽德爾潮流計算法原理簡單，編程實現(xiàn)容易，特別是對于配網(wǎng) 潮流有其獨特優(yōu)勢。但是高斯-賽德爾潮流計算法在牛頓法以及各種解耦法出 現(xiàn)以后似乎成了一種邊緣性的方法。牛頓-拉夫遜法的優(yōu)點是收斂速度快，若 初值選擇較好，算法將具有平方收斂特性，一般迭代45次便可

4、以收斂到一個 非常精確的解，而且其迭代次數(shù)與所計算的網(wǎng)絡(luò)規(guī)模基本無關(guān)。牛頓-拉夫遜 法也具有良好的收斂可靠性對于呈病態(tài)的系統(tǒng)，牛頓-拉夫遜法均能可靠地收斂。牛頓法的缺點是每次迭代的計算量和所需的內(nèi)存量較大。這是因為雅可比 陣元素的數(shù)目約為2(n- 1)X2(n- 1 )個(直角坐標)，且其數(shù)值在迭代過程中不斷 變化。不過，內(nèi)存占用量及每次迭代所需的時間與程序設(shè)計技巧密切相關(guān)。牛 頓-拉夫遜法的可靠收斂取決于一個良好的啟動初值，如果初值選擇不當，算法 有可能根本不收斂或收斂到一個無法運行的解點上。對于正常運行的系統(tǒng)，各 節(jié)電電壓一般均在額定值附近，偏移不會太大，并且各節(jié)電的相角差也不大， 所以對

5、各節(jié)電可以采用統(tǒng)一的電壓初值。P- Q分解法是為了改進牛頓-拉夫遜法 在內(nèi)存占用量及計算速度方面的不足，P- Q分解法根據(jù)電力系統(tǒng)實際運行狀態(tài) 的物理特點，對極坐標形式的牛頓-拉夫遜法修正方程式進行了合理的簡化。 它無論在內(nèi)存占用量還是計算速度方面都比牛頓-拉夫遜法有較大的改進，是目 前計算速度最快的潮流算法。1.2 Matlab的簡介MATLAB是一種交互式、面向?qū)ο蟮某绦蛟O(shè)計語言，廣泛應(yīng)用于工業(yè) 界與學術(shù)界，主要用于矩陣運算，同時在數(shù)值分析、自動控制模擬、數(shù)字信號 處理、動態(tài)分析、繪圖等方面也具有強大的功能。MATLAB程序設(shè)計語言結(jié)構(gòu) 完整，且具有優(yōu)良的移植性，它的基本數(shù)據(jù)元素是不需要定

6、義的數(shù)組。它可以 高效率地解決工業(yè)計算問題，特別是關(guān)于矩陣和矢量的計算。MATLAB與C語 言和FORTRAN語言相比更容易被掌握。通過M語言，可以用類似數(shù)學公式的 方式來編寫算法，大大降低了程序所需的難度并節(jié)省了時間，從而可把主要的 精力集中在算法的構(gòu)思而不是編程上。另外，MATLAB提供了一種特殊的工具：工具箱(TOOLBOXES ).這些工具箱主要包括：信號處理(SIGNALPROCESSING)、控制系統(tǒng)(CONTROL SYSTEMS )、神經(jīng)網(wǎng)絡(luò)(NEURALNETWORKS)、模糊邏輯(FUZZY LOGIC)、小波(WAVELETS)和模擬(SIMULATION)等等。不同領(lǐng)域

7、、不同層次的用戶通過相應(yīng)工具的學習和應(yīng) 用，可以方便地進行計算、分析及設(shè)計工作。MATLAB設(shè)計中，原始數(shù)據(jù)的填 寫格式是很關(guān)鍵的一個環(huán)節(jié)，它與程序使用的方便性和靈活性有著直接的關(guān)系。原始數(shù)據(jù)輸入格式的設(shè)計，主要應(yīng)從使用的角度出發(fā)，原則是簡單明了，便于 修改。1.2.1矩陣的運算常數(shù)與矩陣的運算即是同該矩陣的每一元素進行運算。但需注意進行數(shù)除 時，常數(shù)通常只能做除數(shù)?；竞瘮?shù)運算中，矩陣的函數(shù)運算是矩陣運算中最實用的部分，常用的主 要有以下幾個:det(a)eig(a)求矩陣a的行列式求矩陣a的特征值inv(a) 或 a 八(-1)求矩陣a的逆矩陣rank(a)求矩陣a的秩trace(a)求矩

8、陣a的跡(對角線元素之和)我們在進行工程計算時常常遇到矩陣對應(yīng)元素之間的運算。這種運算不同 于前面講的數(shù)學運算，為有所區(qū)別，我們稱之為數(shù)組運算。數(shù)組的加、減與矩 陣的加、減運算完全相同。而乘除法運算有相當大的區(qū)別，數(shù)組的乘除法是指 兩同維數(shù)組對應(yīng)元素之間的乘除法，它們的運算符為“.*”和“/”或“.”。 前面講過常數(shù)與矩陣的除法運算中常數(shù)只能做除數(shù)。在數(shù)組運算中有了“對應(yīng) 關(guān)系”的規(guī)定，數(shù)組與常數(shù)之間的除法運算沒有任何限制。另外，矩陣的數(shù)組 運算中還有幕運算（運算符為）、指數(shù)運算（exp）、對數(shù)運算（log）、和開 方運算（sqrt）等。有了 “對應(yīng)元素”的規(guī)定，數(shù)組的運算實質(zhì)上就是針對數(shù) 組

9、內(nèi)部的每個元素進行的。矩陣的幕運算與數(shù)組的幕運算有很大的區(qū)別。1.3牛頓拉夫遜法計算潮流分布一、牛頓拉夫遜法求解過程大致可以分為以下步驟：（1）形成節(jié)點導(dǎo)納矩陣（2）將各節(jié)點電壓設(shè)初值U（3）將節(jié)點初值代入相關(guān)求式，求出修正方程式的常數(shù)項向量（4）將節(jié)點電壓初值代入求式，求出雅可比矩陣元素（5）求解修正方程，求修正向量（6）求取節(jié)點電壓的新值（7）檢查是否收斂，如不收斂，則以各節(jié)點電壓的新值作為初值自第3步 重新開始進行狹義次迭代，否則轉(zhuǎn)入下一步（8）計算支路功率分布，PV節(jié)點無功功率和平衡節(jié)點柱入功率。二、直角坐標系計算（1）牛頓-拉夫遜法潮流計算的公式。把牛頓法用于潮流計算，采用直角坐標形

10、式。其中電壓和支路導(dǎo)納可表示為:v =e + jfi i iYij=氣+電(1-1)PQ節(jié)點的輸出有功功率和無功功率是給定的，則第i節(jié)點的給定功率 設(shè)為匕和Qs (稱為注入功率)。假定系統(tǒng)中的第1、2、m節(jié)點為PQ節(jié)點， 對其中每一個節(jié)點的N-R法表達式：KP = P - P = P - e (Ge - B f ) - f (G f + B e ) = 0 TOC o 1-5 h z i is i is i ij j ij j i ij j ij jj=1j=1AQ = Q - Q = Q - f (Ge - Bf ) + e (G f + Be ) = 0i is i is iij j ij

11、 j iij j ij jj =1j =1i =(1、2、m)(1-2)PV節(jié)點的有功功率和節(jié)點電壓幅值是給定的。假定系統(tǒng)中的第m+1、m+2、n-1節(jié)點為PV節(jié)點，則對其中每一 PV節(jié)點可以列寫方程：AP = P - P = P - e (G e - B f ) - f (G f + B e ) = 0j i ij j ij jj=1(1-3)i is i is i ij j ij jTAU2 = U2 U2 = U2 -(e2 + f 2) is i is i ii =(m+1、m+2、n-1)形成雅可比矩陣。當j=i時，對角元素為l = 2 (G c B f) GeB f = N dei

12、j j ij j ii i ii i iii j=l=一乙(G f -B e ) + B e -G f =HOJIJ / J ii I II I IIi j=ldAQ Y=2L(G f -B e ) + B e -G f =Lde歹 j ij j H i H i iii j=l= e -B f) + G e +B f =Jdfij j ij j ii i a ii iii j=l5AI/2j- = lede ，(1-4)當時,矩陣非對角元素為:SAPdAQ = = -(Ge+Bf) = N =-J dedfu 1 u i ij ijSAP dNQ=B e -G f =H =Ldf de ij

13、i ij i ij ij TOC o 1-5 h z jJdAU 2 dAU 2i- =i_ = 0dedfjj(1-5)二、極坐標計算對于潮流計算中待于求出功率的節(jié)點功率方程組，在某個近似解附近用泰 勒級數(shù)展開略去二階及以上的高階得到已矩陣表示的修正方程:APAQ式中節(jié)點n為節(jié)點數(shù)，m為PV節(jié)點，雅閣比矩陣是（2n-2m-2）階非奇異矩陣，雅可比矩陣各元素表示如下:膏=播=陽。呦峋=警=屏國一牛頓拉夫遜極坐標潮流計算的修正方程的迭代方程為:2.1設(shè)計資料及參數(shù):課程名稱電力系統(tǒng)穩(wěn)態(tài)分析課程設(shè)計設(shè)計題目兩機五節(jié)點網(wǎng)絡(luò)潮流計算牛拉法指導(dǎo)教師劉景霞時間1周一、教學要求電力系統(tǒng)穩(wěn)態(tài)分析課程設(shè)計以設(shè)計

14、和優(yōu)化電力系統(tǒng)的潮流分析為重點，提高學生綜合能力為目標，盡可能結(jié)合實際工程進行。設(shè)計內(nèi)容的安排要充分考慮學?，F(xiàn)有的設(shè)備， 設(shè)計時間及工程實際需要，并使學生初步學會運用所學知識解決工程中的實際問題。二、設(shè)計資料及參數(shù)（一）設(shè)計原始資料1、待設(shè)計電氣設(shè)備系統(tǒng)圖2、電力系統(tǒng)網(wǎng)絡(luò)各元件參數(shù)3、電力系統(tǒng)電氣元件的使用規(guī)范4、電力工程電氣設(shè)計手冊（二）設(shè)計參考資料1、電力系統(tǒng)穩(wěn)態(tài)分析，陳珩，中國電力出版社，2007，第三版2、電力系統(tǒng)分析，韓禎祥，浙江大學出版社，2005,第三版3、電力系統(tǒng)分析課程實際設(shè)計與綜合實驗，祝書萍，中國電力出版社，2007，第一版三、設(shè)計要求及成果根據(jù)給定的參數(shù)或工程具體要求，

15、收集和查閱資料；學習相關(guān)軟件（軟件自選）。在給定的電力網(wǎng)絡(luò)上畫出等值電路圖。運用計算機進行潮流計算。編寫設(shè)計說明書?；疽螅壕帉懗绷饔嬎愠绦颍辉谟嬎銠C上調(diào)試通過（？）；運行程序并計算出正確結(jié)果（？）；寫出課程設(shè)計報告（包括以下內(nèi)容）（1份）（1）程序框圖；（2）源程序；（3）符號說明表；（4）算例及計算結(jié)果編寫計算說明書（1份）。四、進度安排根據(jù)給定的參數(shù)或工程具體要求，收集和查閱資料（半天）學習軟件（MATLAB或C語言等）（一天半）編程計算復(fù)雜系統(tǒng)潮流計算（三天）編寫計算設(shè)計書（一天）五、評分標準課程設(shè)計成績評定依據(jù)包括以下幾點：1）工作態(tài)度（占10%）；2）基本技能的掌握程度（占20%

16、）；3）程序編寫是否合理是否有運行結(jié)果（40%）；4）課程設(shè)計說明書編寫水平（占30%）。5）分為優(yōu)、良、中、合格、不合格五個等級。4O.()l+j0.03寸 zorso0.6-jO. i系統(tǒng)接線圖其中節(jié)點1為平衡節(jié)點，節(jié)點2、3、4、5為PQ節(jié)點。2.1.1牛頓拉夫遜法程序框圖電力系統(tǒng)穩(wěn)態(tài)分析課程設(shè)計2.2.1程序的編寫%說明：為了使節(jié)點按照先PQ,再PV節(jié)點，最后平衡節(jié)點的次序編號，以便與公式對照， 節(jié)點1與節(jié)點5對調(diào)。%節(jié)點阻抗矩陣Z=0,0.04+0.12i,0,0.08+0.24i,0;0.04+0.12i,0,0.06+0.18i,0.06+0.18i ,0.02+0.06i;0,

17、0.06+0.18i,0,0.01+0.03i,0.08+0.24i;0.08+0.24i,0.06+ 0.18i,0.01+0.03i,0,0;0,0.02+0.06i,0.08+0.24i,0,0;%求互導(dǎo)納for m=1:5for n=1:5if Z(m,n)=0y(m,n)=0;elsey(m,n)=1/Z(m,n);endendend%求導(dǎo)納for m=1:5for n=1:5if m=nY(m,n)=-y(m,n);elseY(m,n)=sum(y(m,:);endendend%導(dǎo)納矩陣G=real(Y);B=imag(Y);%計算各節(jié)點功率的不平衡量設(shè)U=E+jF ； Y=G+B

18、j；YU=1,1,1,1,1.06;E=real(U);F=imag(U);%設(shè) S=P+Qj；S=-0.60-0.10i;0.20+0.20i;-0.45-0.15i;-0.40-0.05i;0 P=real(S);Q=imag(S);k=0;C=1;while C 0.00001E(5)=1.06;F(5)=0;for m=1:4for n=1:5%計算Pi,Qi,設(shè)Pi=Pt;Qi=Qt,按照書上的公式：Pt(n) = (E(m)*(G(m,n)*E(n)-B(m,n)*F(n)+F(m)*(G(m,n)*F(n)+B(m,n)*E( n);Qt(n) = (F(m)*(G(m,n)*E

19、(n)-B(m,n)*F(n)-E(m)*(G(m,n)*F(n)+B(m,n)*E( n);end%設(shè)P,Q的改變量為dP,dQdP(m)=P(m)-sum(Pt);dQ(m)=Q(m)-sum(Qt);end%計算Hii,Nii,Jii,Lii,由公式4-41b 左側(cè)公式實現(xiàn)，sum(Ai), sum(Bi)用于實現(xiàn) 公式中的sigerma從j到n的求和； for m=1:4for n=1:5Bi(n)=G(m,n)*F(n)+B(m,n)*E(n);Ai(n)=G(m,n)*E(n)-B(m,n)*F(n); endH(m,m)=sum(Bi)-(B(m,m)*E(m)+G(m,m)*F

20、(m)+2*G(m,m)*F(m);N(m,m)=sum(Ai)-(G(m,m)*E(m)-B(m,m)*F(m)+2*G(m,m)*E(m);J(m,m)=-2*B(m,m)*F(m)+sum(Ai)-(G(m,m)*E(m)-B(m,m)*F(m);L(m,m)=-2*B(m,m)*E(m)-(sum(Bi)-(B(m,m)*E(m)+G(m,m)*F(m); end%設(shè)雅可比矩陣為JJ，以下語句用來實現(xiàn)雅可比矩陣中對角線上元素日N J L的排列 for m=1:4JJ(2*m-1,2*m-1)=H(m,m);JJ(2*m-1,2*m)=N(m,m);JJ(2*m,2*m-1)=J(m,m

21、);JJ(2*m,2*m)=L(m,m);end%以下語句用于實現(xiàn)雅可比矩陣非對角線上元素的排列for m=1:4for n=1:4if m=nelseH(m,n)=-B(m,n)*E(m)+G(m,n)*F(m);N(m,n)=G(m,n)*E(m)+B(m,n)*F(m);J(m,n)=-B(m,n)*F(m)-G(m,n)*E(m);L(m,n)=G(m,n)*F(m)-B(m,n)*E(m);JJ(2*m-1,2*n-1)=H(m,n);JJ(2*m-1,2*n)=N(m,n);JJ(2*m,2*n-1)=J(m,n);JJ(2*m,2*n)=L(m,n);endendendJJ%設(shè)由

22、P,Q的改變量組成的8x1矩陣為PQ，由E,F的改變量組成的8x1矩陣為dU for m=1:4PQ(2*m-1)=dP(m); PQ(2*m)=dQ(m);enddU=inv(JJ)*PQ;% 求逆矩陣C=max(abs(dU);%abs絕對值或復(fù)數(shù)絕對值for n=1:4F(n)=F(n)+dU(2*n-1);E(n)=E(n)+dU(2*n);endfor n=1:5U(n)=E(n)+(F(n)*i;endk=k+1;k, dU=dU,PQ,Uend%計算S(5),也就是題目中的S1,即平衡節(jié)點功率for m=1:5I(m)=Y(5,m)*U(m);endS(5)=U(5)*sum(c

23、onj(I)Y =3.7500 -11.2500i-2.5000 + 7.5000i0-1.2500 +3.7500i0-2.5000 + 7.5000i10.8333 -32.5000i-1.6667 + 5.0000i-1.6667 +5.0000i-5.0000 +15.0000i0-1.6667+5.0000i12.9167-38.7500i-10.0000+30.0000i-1.2500 +3.7500i-1.2500+ 3.7500i-1.6667+5.0000i-10.0000+30.0000i12.9167-38.7500i00-5.0000+15.0000i-1.2500+

24、3.7500i06.2500 -18.7500iS =-0.6000 - 0.1000i0.2000 + 0.2000i-0.4500 - 0.1500i-0.4000 - 0.0500i0JJ =11.25003.7500-7.5000-2.500000-3.7500-1.2500-3.750011.25002.5000-7.5000001.2500-3.7500-7.5000-2.500033.400010.5333-5.0000-1.6667-5.0000-1.66672.5000-7.5000-11.133331.60001.6667-5.00001.6667-5.000000-5.0

25、000-1.666738.975012.8417-30.0000-10.0000001.6667-5.0000-12.991738.525010.0000-30.0000-3.7500-1.2500-5.0000-1.6667-30.0000-10.000038.750012.91671.2500-3.75001.6667-5.000010.0000-30.0000-12.916738.7500dU =-0.1076-0.09220.00930.0141-0.04730.0430-0.08630.0154PQ =-0.6000-0.4000-0.1000-0.05000.50001.1000-

26、0.37500.07501.0093-0.1076i1.0430 - 0.0473i1.0154 - 0.0863i1.01410.0922i1.0600JJ =00-3.650511.05164.3956-7.3011-3.3304-1.6652-5.595610.85163.3304-7.3011001.6652-3.6505-7.7040-2.962133.159413.0913-5.1360-1.9747-5.1360-12.9621.9747-7.7040-12.580433.60831.9747-5.13601.9747-5.013600-4.9331-2.123838.38481

27、6.0299-29.5988-12.7427002.1238-4.9331-16.888838.078812.7427-29.5988-3.6876-1.6135-4.9168-2.1513-29.5009-12.907938.1553161.6135.2727-3.68762.1513-4.916812.9079-29.5009-17.072738.0553k =2dU =0.0032-0.0131-0.0004-0.00750.0017-0.01010.0021-0.PQ =0.01640108-0.0636-0.0770-0.0220-0.0008-0.03170.0103-0.0362

28、U =0.9963 - 0.1044i1.0355 - 0.0477i1.0053 - 0.0846i1.0033 -0.0902i1.0600JJ =10.97744.3253-7.2110-3.273800-3.6055-1.6369-5.496010.65563.2738-7.2110001.6369-3.6055-7.6466-2.946632.933412.9530-5.0977-1.9644-5.0977-1.96442.9466-7.6466-12.584333.33711.9644-5.09771.9644-5.097700-4.8857-2.098338.049415.829

29、9-29.3142-12.5898002.0983-4.8857-16.693837.679012.5898-29.3142-3.6499-1.5922-4.8665-2.1230-29.1988-12.738037.799716.06211.5922-3.64992.1230-4.866512.7380-29.1988-16.844337.6306dU =1.0e-003 *0.0241-0.1731-0.0025-0.08920.0098-0.12310.0128-0.1327PQ=1.0e-003*-0.0034-0.8423-0.5236-0.2008-0.0815-0.31980.0

30、241-0.39060.9961 - 0.1044i1.0354 - 0.0477i1.0052 - 0.0845i1.0032 -0.0901i1.0600JJ =10.97644.3243-7.2098-3.273100-3.6049-1.6366-5.495110.65293.2731-7.2098001.6366-3.6049-7.6459-2.946432.930712.9517-5.0973-1.9643-5.0973-1.2.94649643-7.6459-12.583933.33401.9643-5.09731.9643-5.097300-4.8851-2.098038.045

31、115.8278-29.3106-12.5883002.0980-4.8851-16.691937.674012.5883-29.3106-3.6494-1.5920-4.8658-2.1227-29.1950-12.736337.795216.1.59200599-3.64942.1227-4.865812.7363-29.1950-16.842137k =4dU =1.0e-007 *0.0232.6252-0.2889-0.0013-0.13520.0054-0.18760.0087-0.2038PQ =1.0e-006 *-0.0240-0.1486-0.0424-0.0189-0.0

32、165-0.0415-0.0042-0.0525U =0.9961 - 0.1044i1.0354 - 0.0477i1.0052 - 0.0845i1.0032 -0.0901i1.0600S =-0.6000 - 0.1000i0.2000 + 0.2000i-0.4500 - 0.1500i-0.4000 - 0.0500i1.2982 + 0.2445i2.2.3pq法程序編寫基于MATLAB的PQ分解法潮流計算程序MATLAB已廣泛應(yīng)用于自動控制、數(shù)學運算、 信號分析、計算機技術(shù)、圖像信號處理、財務(wù)分析、航天工業(yè)和生物醫(yī)學工程等領(lǐng)域。由 于MATLAB語言功能強大、人際界面友好、編程

33、效率高、強大而智能化的作圖功能，且具 有編程語句簡潔、靈活、表達和運算能力強等顯著特點。程序清單%本程序的功能是用PQ分解法進行潮流計算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 = ）；X=input（ 請輸入由節(jié)點號及其對地阻抗形成的矩陣:X=）；na=input（ 請輸入 PQ 節(jié)點號:na=）;Y=zeros（n）;YI= zeros（

34、n）;e= zeros（1,n）;f= zeros（1,n）;V= zeros（1,n）;O= zeros（1,n）;for i=1:nif X (i、2 ) ?H0、PHX (i、l)、Y (p、p) /X (i、2)、endendfor iHlmlif Bl (i、7 Hophbi(i、l)、qHBl(i、2)、else PH Bl(i、2)、qHBl(i、l)、endY(p、q)HY(p、q)lBIHimag (Yrfor iHlmS (i) HB2 (i、1)IB2 (i、2)、 BI (iH-) H BI (iH-) + B2 (i、 end PHreal(S)、QHimag (S)

35、、 for iHlm e (i) Hreal (B2 (iL)、 f (i) Himag (B2 (iL)、 v(i)HB2 (i、4)、 end for iHlm if B2 (i、6)h H2V (i) Hsqrt (e (i) ？2 + f (i) 72 。(i) Hatan (f(i) /e(i) ) * end end for iH2m if iH Hn B (iH-) Hl /B (iH-)、 else IclHi+l.for j1=IC1:nB(i,j1)= B(i,j1)./B(i,i); end B(i,i)=1./B(i,i); for k=i+1:n for j1=i+1

36、:nB(k,j1)=B(k,j1)-B(k,i)*B(i,j1); end endendendp=0;q=0;for i=1:nif B2(i,6)= =2 p=p+1;k=0;for j1=1:nif B2(i1,6)= =2 k=k+1;A(p,k)=BI(i,j1);end endendendfor i=1:naif i= =naA(i,i)=1./A(i,i);else k=i+1;for j1=k:naA(i,j1)= A(i,j1)./A(i,i); end A(i,i)=1. /A(i,i);for k=i+1:nafor j1=i+1:naA(k,j1)= A(k,j1)-A(

37、k,i)*A(i,j1); end end end end ICT2=1;ICT1=0;kp=1,kq=1;K=1;DET=0;ICT3=1;while ICT2=0|ICT3=0ICT2=0;ICT3=0;for i=1:n if i=isbC(i)=0;for k=1:nC(i)=C(i)+V(k)*(G(i,k)*cos(O(i)-O(k)+BI(i,k)*sin(O(i)-O(k);endDP1(i)=P(i)-V(i)*C(i);DP(i)=DP1(i)./V(i);DET=abs(DP1(i);if DET=prICT2=ICT2+1;endendendNp(K)=ICT2;if

38、ICT2=0for i=2:nDP(i)=B(i,i)*DP(i);if i =nIC1=i+1;for k=ic1:nDP(k)=DP(k)-B(k,i)*DP(i);endelsefor LZ=3:iL=i+3-LZ;IC4=L-1;for MZ=2:IC4I=IC4+2-MZ;DP(I)=DP(I)-B(I,L)*DP(L);endendendendfor i=2:nO(i)=O(i)-DP(i);endkq=1;L=0;for i=1:nif B2(i,6)= =2 C(i)=0;L=L+1; for k=1:n C(i) = C(i) +V(k)*(G(i,k)*sin(O(i)-O

39、(k)+BI(i,k)*cos(O(i)-O(k); end DQ1(i)=Q(i)-V(i)*C(i); DQ(L)=DQ1(i)./V(i); DET=abs(DQ1(i); if DET=pr ICT3=ICT3+1; end end end end else kp=0; if kq=0; L=0; for i=1:nif B2(i,6)= =2 C(i)=0;L=L+1; for k=1:nC(i)=C(i)+V(k)*(G(i,k)*sin(O(i)-O(k)-BI(i,k)*cos(O(i)-O(k); end DQ1(i)=Q(i)-V(i)*C(i); DQ(L)= DQ1(i

40、)./V(i); DET=abs(DQ1(i); end end end end Nq(K)=ICT3; if ICT3=0 L=0; for i=1:na DQ(i)=A(i,i)*DQ(i); if i= =na for LZ=2:i L=i+2-LZ;IC4=L-1;for MZ=1:IC4I=IC4+1-MZ;DQ(I)=DQ(I)-A(I,L)*DQ(L);endendelseIC1=i+1;for k=IC1:naDQ(k)=DQ(k)-A(k,i)*DQ(i);endendendL=0;for i=1:nif B2(i,6)= =2L=L+1;V(i)=V(i)-DQ(L);en

41、dendkp=1;K=K+1;elsekq=0;if kp=0K=K+1;endendfor i=1:nDy(K-1,i)=V(i);endenddisp(迭代次數(shù))；disp (K)；disp (每次沒有達到精度要求的有功功率個數(shù)為)；disp (NP)；disp (每次沒有達到精度要求的無功功率個數(shù)為) disp ( Nq);for k=1:nE(k)=V(k)*cos(O(K)+V(k)*sin(O(k)*j;O(k)=O(k)*180./pi;end disp(各節(jié)點的電壓標么值E為(節(jié)點號從小到大排)：)；disp(E);disp(各節(jié)點的電壓U大小(節(jié)點號從小到大排)：。為；dis

42、p(U);disp(各節(jié)點的電壓相角O (節(jié)點號從小到大排)：。為；disp(O);for p=1:nC(p)=0;for q=1:nC(p)= C(p)+conj(Y(p,q)*conj(E(q);endS(p)=E(p)*C(p);enddisp(各節(jié)點的功率S (節(jié)點號從小到大排)：。為；disp(S)；disp( 各條支路的首端功率Si為(順序同您輸入B1時的一樣)：)為；for i=1:nlif B1(i,6)= =0p=B1(i,1);q=B1(i,2);else p=B1(i,2);q=B1(i,1);endSi(p,q)=E(p)*(conj(E(p)*conj(B1(i,4)

43、./2)+(conj(E(p)*B1(i,5)-conj(E(q)*conj(1./(B1(i,3)*B1(i,5);disp(Si(p,q);enddisp( 各條支路的末端功率Sj (順序同您輸入B1時的一樣)：)為；for i=1:nlif B1(i,6)= =0p=B1(i,1);q=B1(i,2);else p=B1(i,2);q=B1(i,1);endSj(q,p)=E(q)*(conj(E(q)*conj(B1(i,4)./2)+(conj(E(q)*B1(i,5)-conj(E(p)*conj(1./(B1(i,3)*B1(i,5);disp(Sj(q,p);enddisp(

44、各條支路的功率損耗DS (順序同您輸入B1時的一樣)：)為；for i=1:nlif B1(i,6)= =0p=B1(i,1);q=B1(i,2);else p=B1(i,2);q=B1(i,1);endDS(i)=Si(p,q)+Sj(q,p);disp(DS(i);endfor i=1:KCs(i)=i;for j=1:nDy(K,j)=Dy(K-1,j);endenddisp( 以下是每次迭代后各節(jié)點的電壓值(如圖所示)；plot(Cs,Dy),xlable( 迭代次數(shù))，ylable( 電壓)，title( 電壓迭代次數(shù)曲線)；設(shè)計驗證現(xiàn)以一個例子來說明該軟件得使用方法：電力系統(tǒng)接線，

45、如4.1所示，試求潮流分布n為節(jié)點數(shù)、nl為支路數(shù)、isb為平衡母線節(jié)點號(固定為1)、pr為誤差精度。輸入由支路參數(shù)形成的矩陣B1矩陣B1的每行是由下列參數(shù)構(gòu)成的：某支路的首端號P；末端號Q，且P Q；支路的阻抗(R +jX)；支路的對地容抗；支路的變比K；折算到哪一側(cè)的標志(如果支路的首端P處于高壓側(cè)則請輸入“1”，否則請輸入“0”)。Figure 4.1 schematic請輸入各節(jié)點參數(shù)形成的矩陣B2 矩陣B2的每行是由下列參數(shù)構(gòu)成的：節(jié)點所接發(fā)電機的功率SG；節(jié)點負荷的功率SL；節(jié)點電壓的初始值；PU節(jié)點電壓U的給定值；節(jié)點所接的無功補償設(shè)備的容量；節(jié)點分類標號igl。輸入數(shù)據(jù)：請輸入節(jié)點數(shù)：n=5請輸入支路數(shù)：nl=5請輸入平衡母線節(jié)點號：isb=1請輸入誤差精度：pr=0.00001請輸入由支路參數(shù)形成的矩陣：B1=1 2 0.03i 0 1.050;2 3 0.08 + 0.3i 0.5i 1 0;2 4 0.1+0.35i 0 1 0;3 4 0.04+0.25i 0.5i 1 0;3 5 0.015i 0 1.05 1 請輸入由節(jié)點參數(shù)形成的矩陣：B2=0 0 1.05 1.05 0 1;0 3.7 + 1.3i 1.05 0 0 2;0 2+1i 1.05 0 0 2;0 1.6+0.8i 1.05 0