兩機五節(jié)點網(wǎng)絡(luò)潮流計算方法牛拉法和pq法電力系統(tǒng)穩(wěn)態(tài)分析課程設(shè)計1

1、電力系統(tǒng)穩(wěn)態(tài)分析課程設(shè)計題 目 名 稱 兩機五節(jié)點網(wǎng)絡(luò)潮流計算方法牛拉法和pq法 學(xué) 生 姓 名 學(xué) 號 系 、專 業(yè) 電氣工程系電氣工程及其自動化 指 導(dǎo) 教 師 2012年12月28日目 錄摘要.第一章 原理簡介.3 1.1對潮流分析的簡介.3 1.1.1 潮流計算方法分析比較.31.2 matlab簡介.4 1.2.1 矩陣的運算51.3牛頓拉夫遜法計算潮流分布.6第二章 程序及結(jié)果.102.1 設(shè)計資料及參數(shù).10 2.1.1 牛頓拉夫遜法的程序框圖132.2 用matlab設(shè)計程序.14 2.2.1 程序的編寫142.2.2程序運行結(jié)果.192.23p_q法程序編寫22 總結(jié).32

2、參考文獻.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ù)給定的運行條件及系統(tǒng)接線情況確定整個電力系統(tǒng)各個部分的運行狀態(tài)，如各母線的電壓、各元件中流過的功率、系統(tǒng)的功率損耗等等。電力系統(tǒng)潮流計算

3、是計算系統(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 次便可以收斂到一個非常精確的解，而且其迭代次數(shù)與所計算的網(wǎng)絡(luò)規(guī)模基本無關(guān)。牛頓 - 拉夫遜法也具有良好的收斂可靠性 對于呈病態(tài)的系統(tǒng)，牛頓

4、-拉夫遜法均能可靠地收斂。牛頓法的缺點是每次迭代的計算量和所需的內(nèi)存量較大。這是因為雅可比陣元素的數(shù)目約為 2(n- 1)2(n- 1)個（直角坐標），且其數(shù)值在迭代過程中不斷變化。不過，內(nèi)存占用量及每次迭代所需的時間與程序設(shè)計技巧密切相關(guān)。牛頓-拉夫遜法的可靠收斂取決于一個良好的啟動初值，如果初值選擇不當，算法有可能根本不收斂或收斂到一個無法運行的解點上。對于正常運行的系統(tǒng)，各節(jié)電電壓一般均在額定值附近，偏移不會太大，并且各節(jié)電的相角差也不大，所以對各節(jié)電可以采用統(tǒng)一的電壓初值。p- q分解法是為了改進牛頓-拉夫遜法在內(nèi)存占用量及計算速度方面的不足，p- q分解法根據(jù)電力系統(tǒng)實際運行狀態(tài)的物

5、理特點，對極坐標形式的牛頓 - 拉夫遜法修正方程式進行了合理的簡化。它無論在內(nèi)存占用量還是計算速度方面都比牛頓-拉夫遜法有較大的改進，是目前計算速度最快的潮流算法。1.2 matlab 的簡介 matlab是一種交互式、面向?qū)ο蟮某绦蛟O(shè)計語言，廣泛應(yīng)用于工業(yè)界與學(xué)術(shù)界，主要用于矩陣運算，同時在數(shù)值分析、自動控制模擬、數(shù)字信號處理、動態(tài)分析、繪圖等方面也具有強大的功能。matlab程序設(shè)計語言結(jié)構(gòu)完整，且具有優(yōu)良的移植性，它的基本數(shù)據(jù)元素是不需要定義的數(shù)組。它可以高效率地解決工業(yè)計算問題，特別是關(guān)于矩陣和矢量的計算。matlab與c語言和fortran語言相比更容易被掌握。通過m語言，可以用類似

6、數(shù)學(xué)公式的方式來編寫算法，大大降低了程序所需的難度并節(jié)省了時間，從而可把主要的精力集中在算法的構(gòu)思而不是編程上。另外，matlab提供了一種特殊的工具：工具箱（toolboxes）.這些工具箱主要包括：信號處理（signal processing）、控制系統(tǒng)（control systems）、神經(jīng)網(wǎng)絡(luò)（neural networks）、模糊邏輯(fuzzy logic)、小波(wavelets)和模擬（simulation）等等。不同領(lǐng)域、不同層次的用戶通過相應(yīng)工具的學(xué)習(xí)和應(yīng)用，可以方便地進行計算、分析及設(shè)計工作。matlab設(shè)計中，原始數(shù)據(jù)的填寫格式是很關(guān)鍵的一個環(huán)節(jié)，它與程序使用的方便性和

7、靈活性有著直接的關(guān)系。原始數(shù)據(jù)輸入格式的設(shè)計，主要應(yīng)從使用的角度出發(fā)，原則是簡單明了，便于修改。 1.2.1 矩陣的運算常數(shù)與矩陣的運算即是同該矩陣的每一元素進行運算。但需注意進行數(shù)除時，常數(shù)通常只能做除數(shù)?；竞瘮?shù)運算中，矩陣的函數(shù)運算是矩陣運算中最實用的部分，常用的主要有以下幾個：det(a) 求矩陣a的行列式eig(a) 求矩陣a的特征值inv(a)或a (-1) 求矩陣a的逆矩陣rank(a) 求矩陣a的秩trace(a) 求矩陣a的跡（對角線元素之和）我們在進行工程計算時常常遇到矩陣對應(yīng)元素之間的運算。這種運算不同于前面講的數(shù)學(xué)運算，為有所區(qū)別，我們稱之為數(shù)組運算。數(shù)組的加、減與矩陣

8、的加、減運算完全相同。而乘除法運算有相當大的區(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ù)組內(nèi)部的每個元素進行的。矩陣的冪運算與數(shù)組的冪運算有很大的區(qū)別。1.3牛頓拉夫遜法計算潮流分布一、牛頓拉夫遜法求解過程大致可以分為以下步驟：（1）形成節(jié)點導(dǎo)納矩陣（2）將各

9、節(jié)點電壓設(shè)初值u（3）將節(jié)點初值代入相關(guān)求式，求出修正方程式的常數(shù)項向量（4）將節(jié)點電壓初值代入求式，求出雅可比矩陣元素（5）求解修正方程，求修正向量（6）求取節(jié)點電壓的新值（7）檢查是否收斂，如不收斂，則以各節(jié)點電壓的新值作為初值自第3步 重新開始進行狹義次迭代，否則轉(zhuǎn)入下一步（8）計算支路功率分布，pv節(jié)點無功功率和平衡節(jié)點柱入功率。二、直角坐標系計算 (1)牛頓-拉夫遜法潮流計算的公式。把牛頓法用于潮流計算，采用直角坐標形式。其中電壓和支路導(dǎo)納可表示為： （1-1） (2)pq節(jié)點的輸出有功功率和無功功率是給定的，則第i節(jié)點的給定功率設(shè)為和（稱為注入功率）。假定系統(tǒng)中的第1、2、m節(jié)點為

10、pq節(jié)點，對其中每一個節(jié)點的n-r法表達式: =（1、2、m） （1-2） (3)pv節(jié)點的有功功率和節(jié)點電壓幅值是給定的。假定系統(tǒng)中的第m+1、m+2、n-1節(jié)點為pv節(jié)點，則對其中每一pv節(jié)點可以列寫方程： （1-3） =（m+1、m+2、n-1）(4)形成雅可比矩陣。 當j=i時，對角元素為 （1-4） 當時,矩陣非對角元素為： （1-5）三、極坐標計算 對于潮流計算中待于求出功率的節(jié)點功率方程組，在某個近似解附近用泰勒級數(shù)展開略去二階及以上的高階得到已矩陣表示的修正方程： 式中節(jié)點n為節(jié)點數(shù)，m為pv節(jié)點，雅閣比矩陣是（2n-2m-2）階非奇異矩陣，雅可比矩陣各元素表示如下：牛頓拉夫遜

11、極坐標潮流計算的修正方程的迭代方程為：2.1設(shè)計資料及參數(shù)：課程名稱電力系統(tǒng)穩(wěn)態(tài)分析課程設(shè)計設(shè)計題目兩機五節(jié)點網(wǎng)絡(luò)潮流計算牛拉法指導(dǎo)教師劉景霞時間1周一、教學(xué)要求電力系統(tǒng)穩(wěn)態(tài)分析課程設(shè)計以設(shè)計和優(yōu)化電力系統(tǒng)的潮流分析為重點，提高學(xué)生綜合能力為目標，盡可能結(jié)合實際工程進行。設(shè)計內(nèi)容的安排要充分考慮學(xué)?，F(xiàn)有的設(shè)備，設(shè)計時間及工程實際需要，并使學(xué)生初步學(xué)會運用所學(xué)知識解決工程中的實際問題。二、設(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，第

12、三版2、電力系統(tǒng)分析，韓禎祥，浙江大學(xué)出版社，2005，第三版3、電力系統(tǒng)分析課程實際設(shè)計與綜合實驗，祝書萍，中國電力出版社，2007，第一版三、設(shè)計要求及成果1.根據(jù)給定的參數(shù)或工程具體要求，收集和查閱資料；學(xué)習(xí)相關(guān)軟件（軟件自選）。2.在給定的電力網(wǎng)絡(luò)上畫出等值電路圖。3.運用計算機進行潮流計算。4.編寫設(shè)計說明書?；疽螅?.編寫潮流計算程序；2.在計算機上調(diào)試通過（？）；3.運行程序并計算出正確結(jié)果（？）；4.寫出課程設(shè)計報告（包括以下內(nèi)容）（1份）（1）程序框圖；（2）源程序；（3）符號說明表；（4）算例及計算結(jié)果5.編寫計算說明書（1份）。四、進度安排根據(jù)給定的參數(shù)或工程具體要求

13、，收集和查閱資料（半天）學(xué)習(xí)軟件（matlab或c語言等）（一天半）編程計算復(fù)雜系統(tǒng)潮流計算（三天）編寫計算設(shè)計書（一天）五、評分標準課程設(shè)計成績評定依據(jù)包括以下幾點：1） 工作態(tài)度（占10%）；2） 基本技能的掌握程度（占20%）；3） 程序編寫是否合理是否有運行結(jié)果(40%)；4） 課程設(shè)計說明書編寫水平(占30%)。5） 分為優(yōu)、良、中、合格、不合格五個等級。系統(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é)點

14、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,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:5 for n=1:5 if z(m,n)=0 y(m,n)=0; else y(m,n)=1/z(m,n); end endend%求導(dǎo)納for m=1:5 for n=1:5 if m

15、=n y(m,n)=-y(m,n); else y(m,n)=sum(y(m,:); end end end%導(dǎo)納矩陣g=real(y);b=imag(y);%計算各節(jié)點功率的不平衡量設(shè)u=e+jf ；y=g+bj；y u=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:4 for n=1:5%計算pi,qi,設(shè)

16、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(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:4 for

17、n=1:5 bi(n)=g(m,n)*f(n)+b(m,n)*e(n); ai(n)=g(m,n)*e(n)-b(m,n)*f(n); end h(m,m)=sum(bi)-(b(m,m)*e(m)+g(m,m)*f(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 %

18、設(shè)雅可比矩陣為jj，以下語句用來實現(xiàn)雅可比矩陣中對角線上元素h 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);jj(2*m,2*m)=l(m,m);end%以下語句用于實現(xiàn)雅可比矩陣非對角線上元素的排列for m=1:4 for n=1:4 if m=n else h(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(

19、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); end endendjj%設(shè)由p,q的改變量組成的81矩陣為pq，由e,f的改變量組成的81矩陣為dufor m=1:4 pq(2*m-1)=dp(m); pq(2*m)=dq(m);enddu=inv(jj)*pq;%求逆矩陣c=max(abs(du);%abs絕對值或復(fù)數(shù)絕對值for n=1:4 f(n)=f(n)+du(2*n-1); e(n)=e(n)+du(2*n);en

20、d for n=1:5 u(n)=e(n)+(f(n)*i; end k=k+1;k, du=du,pq,uend%計算s(5),也就是題目中的s1，即平衡節(jié)點功率for m=1:5 i(m)=y(5,m)*u(m); end s(5)=u(5)*sum(conj(i)y = 3.7500 -11.2500i -2.5000 + 7.5000i 0 -1.2500 + 3.7500i 0 -2.5000 + 7.5000i 10.8333 -32.5000i -1.6667 + 5.0000i -1.6667 + 5.0000i -5.0000 +15.0000i 0 -1.6667 + 5.

21、0000i 12.9167 -38.7500i -10.0000 +30.0000i -1.2500 + 3.7500i -1.2500 + 3.7500i -1.6667 + 5.0000i -10.0000 +30.0000i 12.9167 -38.7500i 0 0 -5.0000 +15.0000i -1.2500 + 3.7500i 0 6.2500 -18.7500is = -0.6000 - 0.1000i 0.2000 + 0.2000i -0.4500 - 0.1500i -0.4000 - 0.0500i 0 jj = 11.2500 3.7500 -7.5000 -2.

22、5000 0 0 -3.7500 -1.2500 -3.7500 11.2500 2.5000 -7.5000 0 0 1.2500 -3.7500 -7.5000 -2.5000 33.4000 10.5333 -5.0000 -1.6667 -5.0000 -1.6667 2.5000 -7.5000 -11.1333 31.6000 1.6667 -5.0000 1.6667 -5.0000 0 0 -5.0000 -1.6667 38.9750 12.8417 -30.0000 -10.0000 0 0 1.6667 -5.0000 -12.9917 38.5250 10.0000 -

23、30.0000 -3.7500 -1.2500 -5.0000 -1.6667 -30.0000 -10.0000 38.7500 12.9167 1.2500 -3.7500 1.6667 -5.0000 10.0000 -30.0000 -12.9167 38.7500k = 1du = -0.1076 0.0093 -0.0473 0.0430 -0.0863 0.0154 -0.0922 0.0141pq = -0.6000 -0.1000 0.5000 1.1000 -0.3750 0.0750 -0.4000 -0.0500u = 1.0093 - 0.1076i 1.0430 -

24、 0.0473i 1.0154 - 0.0863i 1.0141 - 0.0922i 1.0600 jj = 11.0516 4.3956 -7.3011 -3.3304 0 0 -3.6505 -1.6652 -5.5956 10.8516 3.3304 -7.3011 0 0 1.6652 -3.6505 -7.7040 -2.9621 33.1594 13.0913 -5.1360 -1.9747 -5.1360 -1.9747 2.9621 -7.7040 -12.5804 33.6083 1.9747 -5.1360 1.9747 -5.1360 0 0 -4.9331 -2.123

25、8 38.3848 16.0299 -29.5988 -12.7427 0 0 2.1238 -4.9331 -16.8888 38.0788 12.7427 -29.5988 -3.6876 -1.6135 -4.9168 -2.1513 -29.5009 -12.9079 38.1553 16.2727 1.6135 -3.6876 2.1513 -4.9168 12.9079 -29.5009 -17.0727 38.0553k = 2du = 0.0032 -0.0131 -0.0004 -0.0075 0.0017 -0.0101 0.0021 -0.0108pq = 0.0164

26、-0.0636 -0.0770 -0.0220 -0.0008 -0.0317 0.0103 -0.0362u = 0.9963 - 0.1044i 1.0355 - 0.0477i 1.0053 - 0.0846i 1.0033 - 0.0902i 1.0600 jj = 10.9774 4.3253 -7.2110 -3.2738 0 0 -3.6055 -1.6369 -5.4960 10.6556 3.2738 -7.2110 0 0 1.6369 -3.6055 -7.6466 -2.9466 32.9334 12.9530 -5.0977 -1.9644 -5.0977 -1.96

27、44 2.9466 -7.6466 -12.5843 33.3371 1.9644 -5.0977 1.9644 -5.0977 0 0 -4.8857 -2.0983 38.0494 15.8299 -29.3142 -12.5898 0 0 2.0983 -4.8857 -16.6938 37.6790 12.5898 -29.3142 -3.6499 -1.5922 -4.8665 -2.1230 -29.1988 -12.7380 37.7997 16.0621 1.5922 -3.6499 2.1230 -4.8665 12.7380 -29.1988 -16.8443 37.630

28、6k = 3du = 1.0e-003 * 0.0241 -0.1731 -0.0025 -0.0892 0.0098 -0.1231 0.0128 -0.1327pq = 1.0e-003 * -0.0034 -0.8423 -0.5236 -0.2008 -0.0815 -0.3198 0.0241 -0.3906u = 0.9961 - 0.1044i 1.0354 - 0.0477i 1.0052 - 0.0845i 1.0032 - 0.0901i 1.0600 jj = 10.9764 4.3243 -7.2098 -3.2731 0 0 -3.6049 -1.6366 -5.49

29、51 10.6529 3.2731 -7.2098 0 0 1.6366 -3.6049 -7.6459 -2.9464 32.9307 12.9517 -5.0973 -1.9643 -5.0973 -1.9643 2.9464 -7.6459 -12.5839 33.3340 1.9643 -5.0973 1.9643 -5.0973 0 0 -4.8851 -2.0980 38.0451 15.8278 -29.3106 -12.5883 0 0 2.0980 -4.8851 -16.6919 37.6740 12.5883 -29.3106 -3.6494 -1.5920 -4.865

30、8 -2.1227 -29.1950 -12.7363 37.7952 16.0599 1.5920 -3.6494 2.1227 -4.8658 12.7363 -29.1950 -16.8421 37.6252k = 4du = 1.0e-007 * 0.0232 -0.2889 -0.0013 -0.1352 0.0054 -0.1876 0.0087 -0.2038pq = 1.0e-006 * -0.0240 -0.1486 -0.0424 -0.0189 -0.0165 -0.0415 -0.0042 -0.0525u = 0.9961 - 0.1044i 1.0354 - 0.0

31、477i 1.0052 - 0.0845i 1.0032 - 0.0901i 1.0600 s = -0.6000 - 0.1000i 0.2000 + 0.2000i -0.4500 - 0.1500i -0.4000 - 0.0500i 1.2982 + 0.2445i2.2.3pq法程序編寫基于matlab 的pq 分解法潮流計算程序matlab 已廣泛應(yīng)用于自動控制、數(shù)學(xué)運算、信號分析、計算機技術(shù)、圖像信號處理、財務(wù)分析、航天工業(yè)和生物醫(yī)學(xué)工程等領(lǐng)域。由于matlab語言功能強大、人際界面友好、編程效率高、強大而智能化的作圖功能，且具有編程語句簡潔、靈活、表達和運算能力強等顯著特點。程

32、序清單%本程序的功能是用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(n);e= zeros(1,n);f= zeros(1,n);v= zeros(1,n);o= zeros(1,n

33、);for i=1:nif x(i,2)=0; p=x(i,1); y(p,p)=1./x(i,2); endendfor i=1:nlif b1(i,6)= =0 p=b1(i,1);q=b1(i,2);else p= b1(i,2);q=b1(i,1);endy(p,q)=y(p,q)-1./(b1(i,3)*b1(i,5);yi(p,q)= yi(p,q) -1./b1(i,3);y(q,p) =y(p,q);yi(q,p) = yi(p,q);y(q,q)= y(q,q)+1./(b1(i,3)* b1(i,5)2)+b1(i,4)./2;yi(q,q)= yi(q,q) +1./b1

34、(i,3);y(p,p)= y(p,p) +1./b1(i,3)+ b1(i,4)./2;yi(p,p)= yi(p,p) +1./b1(i,3);end %求導(dǎo)納矩陣g=real(y);b=imag(yi);bi=imag(y);for i=1:n s(i)=b2(i,1)-b2(i,2); bi(i,i)= bi(i,i)+ b2(i,5);endp=real(s);q=imag(s);for i=1:ne(i)=real(b2(i,3);f(i)=imag(b2(i,3);v(i)=b2(i,4);endfor i=1:nif b2(i,6)= =2 v(i)=sqrt(e(i)2+f(

35、i)2); o(i)=atan(f(i)./e(i);endendfor i=2:nif i= =n b(i,i)=1./b(i,i);else ic1=i+1; for j1=ic1:n b(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:n b(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:n if b2(i1,6)= =2 k=k+1; a(p,k)=bi(

36、i,j1); end endendendfor i=1:naif i= =na a(i,i)=1./a(i,i);else k=i+1; for j1=k:na a(i,j1)= a(i,j1)./a(i,i); end a(i,i)=1. /a(i,i); for k=i+1:na for j1=i+1:na a(k,j1)= a(k,j1)-a(k,i)*a(i,j1); endendendendict2=1;ict1=0;kp=1,kq=1;k=1;det=0;ict3=1;while ict2=0|ict3=0ict2=0;ict3=0; for i=1:nif i=isb c(i)=

37、0; for k=1:n c(i)=c(i)+v(k)*(g(i,k)*cos(o(i)-o(k)+bi(i,k)*sin(o(i)-o(k); end dp1(i)=p(i)-v(i)*c(i); dp(i)=dp1(i)./v(i); det=abs(dp1(i); if det=pr ict2=ict2+1; endendendnp(k)=ict2; if ict2=0 for i=2:n dp(i)=b(i,i)*dp(i); if i =n ic1=i+1; for k=ic1:n dp(k)=dp(k)-b(k,i)*dp(i); end else for lz=3:i l=i+3

38、-lz; ic4=l-1; for mz=2:ic4 i=ic4+2-mz; dp(i)=dp(i)-b(i,l)*dp(l); end end endendfor i=2:n o(i)=o(i)-dp(i);endkq=1;l=0;for i=1:n if b2(i,6)= =2c(i)=0;l=l+1;for k=1:n c(i)= c(i) +v(k)*(g(i,k)*sin(o(i)-o(k)+bi(i,k)*cos(o(i)-o(k);enddq1(i)=q(i)-v(i)*c(i);dq(l)=dq1(i)./v(i);det=abs(dq1(i);if det=pr ict3=i

39、ct3+1;endendendendelse kp=0; if kq=0;l=0;for i=1:n if 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(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);endendendendnq(k)=ict3;if ict3=0 l=0; for i=1:nadq(i)=a(i,i)*dq(i);if i= =na for lz=2

40、:i l=i+2-lz; ic4=l-1; for mz=1:ic4 i=ic4+1-mz; dq(i)=dq(i)-a(i,l)*dq(l); endendelse ic1=i+1; for k=ic1:na dq(k)=dq(k)-a(k,i)*dq(i);endendendl=0;for i=1:nif b2(i,6)= =2 l=l+1; v(i)=v(i)-dq(l);endend kp=1;k=k+1;else kq=0; if kp=0k=k+1; endendfor i=1:n dy(k-1,i)=v(i);endenddisp(迭代次數(shù))；disp（k）；disp（每次沒有達

41、到精度要求的有功功率個數(shù)為）；disp（np）；disp（每次沒有達到精度要求的無功功率個數(shù)為）；disp（nq）;for k=1:n e(k)=v(k)*cos(o(k)+v(k)*sin(o(k)*j; o(k)=o(k)*180./pi;enddisp(各節(jié)點的電壓標么值e為（節(jié)點號從小到大排）：);disp(e);disp(各節(jié)點的電壓u大?。ü?jié)點號從小到大排）：)為；disp(u);disp(各節(jié)點的電壓相角o（節(jié)點號從小到大排）：)為；disp(o);for p=1:n c(p)=0; for q=1:nc(p)= c(p)+conj(y(p,q)*conj(e(q); end s

42、(p)=e(p)*c(p);enddisp(各節(jié)點的功率s（節(jié)點號從小到大排）：)為；disp（s）；disp(各條支路的首端功率si為（順序同您輸入b1時的一樣）：）為；for i=1:nl if 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)./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時的一

43、樣）：）為；for i=1:nl if 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(各條支路的功率損耗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); end

44、ds(i)=si(p,q)+sj(q,p); disp(ds(i);endfor i=1:k cs(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)接線，如4.1所示,試求潮流分布（1）n為節(jié)點數(shù)、nl為支路數(shù)、isb為平衡母線節(jié)點號（固定為 1）、pr為誤差精度。（2）輸入由支路參數(shù)形成的矩陣 b1矩陣 b1的每行是由下列參數(shù)構(gòu)成的： 某支路的首端號

45、p； 末端號 q，且 p q； 支路的阻抗（r +jx）； 支路的對地容抗； 支路的變比 k； 折算到哪一側(cè)的標志（如果支路的首端 p 處于高壓側(cè)則請輸入“1”，否則請輸入“0”）。figure 4.1 schematic（3）請輸入各節(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 0 2;5 0 1.05 1.05 0 3請輸入由