牛頓拉夫遜法潮流計算

1、目錄摘要11.設計意義與要求21.1設計意義21.2設計要求22.牛頓拉夫遜算法32.1牛頓算法數(shù)學原理：32.2 直角坐標系下牛頓法潮流計算的原理43 詳細設計過程93.1節(jié)點類型93.2待求量93.3導納矩陣93.4潮流方程103.5修正方程114.程序設計144.1 節(jié)點導納矩陣的形成144.2 計算各節(jié)點不平衡量154.3 雅克比矩陣計算- 17 -4.4 LU分解法求修正方程- 19 -4.5 計算網絡中功率分布- 22 -5.結果分析- 22 -6.小結- 25 -參考文獻- 26 -附錄：- 27 -摘要潮流計算是電力網絡設計及運行中最基本的計算，對電力網絡的各種設計方案及各種運

2、行方式進行潮流計算，可以得到各種電網各節(jié)點的電壓，并求得網絡的潮流及網絡中各元件的電力損耗，進而求得電能損耗。在數(shù)學上是多元非線性方程組的求解問題，求解的方法有很多種。牛頓拉夫遜法是數(shù)學上解非線性方程式的有效方法，有較好的收斂性。將牛頓法用于潮流計算是以導納矩陣為基礎的，由于利用了導納矩陣的對稱性、稀疏性及節(jié)點編號順序優(yōu)化等技巧，使牛頓法在收斂性、占用內存、計算速度等方面都達到了一定的要求。本文以一個具體例子分析潮流計算的具體方法，并運用牛頓拉夫遜算法求解線性方程關鍵詞： 電力系統(tǒng) 潮流計算 牛頓拉夫遜算法1.設計意義與要求1.1設計意義潮流計算是電力系統(tǒng)分析中的一種最基本的計算，他的任務是對

3、給定運行條件確定系統(tǒng)運行狀態(tài)，如各母線上的電壓（幅值及相角）、網絡中的功率分布及功率損耗等。潮流計算的結果是電力系統(tǒng)穩(wěn)定計算和故障分析的基礎。具體表現(xiàn)在以下方面： (1)在電網規(guī)劃階段,通過潮流計算,合理規(guī)劃電源容量及接入點,合理規(guī)劃網架,選擇無功補償 方案,滿足規(guī)劃水平的大、小方式下潮流交換控制、調峰、調相、調壓的要求。 (2)在編制年運行方式時,在預計負荷增長及新設備投運基礎上,選擇典型方式進行潮流計算,發(fā) 現(xiàn)電網中薄弱環(huán)節(jié),供調度員日常調度控制參考,并對規(guī)劃、基建部門提出改進網架結構,加快基建進度的建議。 (3)正常檢修及特殊運行方式下的潮流計算,用于日運行方式的編制,指導發(fā)電廠開機方式

4、,有 功、無功調整方案及負荷調整方案,滿足線路、變壓器熱穩(wěn)定要求及電壓質量要求。 (4)預想事故、設備退出運行對靜態(tài)安全的影響分析及作出預想的運行方式調整方案。 總結為在電力系統(tǒng)運行方式和規(guī)劃方案的研究中，都需要進行潮流計算以比較運行方式或規(guī)劃供電方 案的可行性、可靠性和經濟性。同時，為了實時監(jiān)控電力系統(tǒng)的運行狀態(tài)，也需要進行大量而快速的潮流計算。因此，潮流計算是電力系統(tǒng)中應用最廣泛、最基本和 最重要的一種電氣運算。在系統(tǒng)規(guī)劃設計和安排系統(tǒng)的運行方式時，采用離線潮流計算；在電力系統(tǒng)運行狀態(tài)的實時監(jiān)控中，則采用在線潮流計算。1.2設計要求1）根據(jù)給定的運行條件，確定圖中電力系統(tǒng)潮流計算時各節(jié)點的

5、類型、待求量；2）求節(jié)點導納矩陣；3）給出潮流方程或功率方程的表達式；4）當用牛頓拉夫遜法計算潮流時，給出修正方程和迭代收斂條件；2.牛頓拉夫遜算法2.1牛頓算法數(shù)學原理：牛頓法 （Newton Method）：解非線性方程f(x)=0的牛頓(Newton) 法，就是將非線性方程線性化的一種方法。它是解代數(shù)方程和超越方程的有效方法之一。設有單變量非線性方程，給出解的近似值，它與真解的誤差為，則將滿足，即 將上式左邊的函數(shù)在附近展成泰勒級數(shù)，如果差值很小，二次及以上階次的各項均可略去得：這是對于變量的修正量的線性方程式，成為修正方程，解此方程可得修正量用所求得的去修正近似解，便得修正后的近似解同

6、真解仍然有誤差。為了進一步逼近真解，可以反復進行迭代計算，迭代計算通式是迭代過程的收斂判據(jù)為式中，和為預先給定的小正數(shù)。牛頓-拉夫遜法實質上就是切線法，是一種逐步線性化的方法，此法不僅用于求單變量方程，也適用于多變量非線性代數(shù)方程的有效方法。牛頓法至少是二階收斂的，即牛頓法在單根附近至少是二階收斂的，在重根附近是線性收斂的。 牛頓法收斂很快，而且可求復根，缺點是對重根收斂較慢，要求函數(shù)的一階導數(shù)存在。2.2 直角坐標系下牛頓法潮流計算的原理采用直角坐標時，節(jié)點電壓可表示為導納矩陣元素則表示為將上述表示式代入的右端，展開并分出實部和虛部，便得 假定系統(tǒng)中的第1，2，3，··&

7、#183;，m號節(jié)點為PQ節(jié)點，第i個節(jié)點的給定功率設為和，對該節(jié)點可列寫方程 （i=1,2，···，m） 假定系統(tǒng)中的第m+1，m+2，···，n-1號節(jié)點為PV節(jié)點，則對其中每一個節(jié)點可以列寫方程（i=m+1，m+2，···，n-1）第n號節(jié)點為平衡點，其電壓是給定的，故不參加迭代。以上兩個方程組總共包含了2（n-1）個方程，待求的變量有也是2（n-1）個。我們還可看到，上面兩個方程式已經具備了方程組的形式。因此，不難寫出如下的修正方程式 式中上述方程中雅克比矩陣的各元素，可以對上式求偏導數(shù)獲得。當時

8、當時修正方程式還可以寫成分塊矩陣的形式式中，和都是二維列向量；是介方陣。對于PQ節(jié)點對于PV節(jié)點從以上表達式可以看到，雅克比矩陣有以下特點：（1） 雅克比矩陣各元素都是節(jié)點電壓的函數(shù)，它們的數(shù)值將在迭代過程中不斷的改變。（2） 雅克比矩陣的子塊中的元素的表達式只用到導納矩陣中的對應元素。若，則必有。因此，式中分塊形式的雅克比矩陣同節(jié)點導納矩陣一樣稀疏，修正方程的求解同樣可以用稀疏矩陣的求解技巧。（3） 雅克比矩陣的元素或子塊都不具有對稱性。用牛頓-拉夫遜法計算潮流的流程：首先要輸入網絡的原始數(shù)據(jù)以及各節(jié)點的給定值并形成節(jié)點導納矩陣。輸入節(jié)點電壓初值和，置迭代計數(shù)k=0。然后開始進入牛頓法的迭代

9、過程。在進行第k+1次迭代時，其計算步驟如下：（1） 按上一次迭代計算出的節(jié)點電壓值和，計算各類節(jié)點的不平衡量、和。（2） 按條件校驗收斂，即< 如果收斂，迭代到此結束，轉入計算各線路潮流和平衡節(jié)點的功率，并打印輸出計算結果。不收斂則繼續(xù)計算。（3）計算雅克比矩陣的各元素。（4）解修正方程式，求節(jié)點電壓的修正量和。（5）修正各節(jié)點的電壓（6）迭代計數(shù)加1，返回第一步繼續(xù)迭代過程。輸入原始數(shù)據(jù)形成節(jié)點導納矩陣給定節(jié)點電壓初值k=0計算<是否計算雅克比矩陣各元素解修正方程式，求，計算平衡節(jié)點功率輸出圖1牛頓-拉夫遜法潮流計算程序框圖3 詳細設計過程3.1節(jié)點類型電力系統(tǒng)潮流計算中，節(jié)點

10、一般分為如下幾種類型：PQ節(jié)點：節(jié)點注入的有功功率無功功率是已知的PV節(jié)點：節(jié)點注入的有功功率已知，節(jié)點電壓幅值恒定，一般由無功儲備比較充足的電廠和電站充當；平衡節(jié)點：節(jié)點的電壓為1*exp(0°)，其注入的有功無功功率可以任意調節(jié)，一般由具有調頻發(fā)電廠充當。更復雜的潮流計算，還有其他節(jié)點，或者是這三種節(jié)點的組合，在一定條件下可以相互轉換。對于本題目，節(jié)點分析如下：節(jié)點1給出有功功率為2，無功功率為1， PQ節(jié)點。節(jié)點2給出有功功率為0.5，電壓幅值為1.0，PV節(jié)點。節(jié)點3電壓相位是0，電壓幅值為1，平衡節(jié)點。3.2待求量節(jié)點1待求量是V，； 節(jié)點2待求量是Q，；節(jié)點3待求量是P，

11、Q。3.3導納矩陣導納矩陣分為節(jié)點導納矩陣、結點導納矩陣、支路導納矩陣、二端口導納矩陣。結點導納矩陣：對于一個給定的電路(網絡)，由其關聯(lián)矩陣A與支路導納矩陣Y所確定的矩陣。支路導納矩陣：表示一個電路中各支路導納參數(shù)的矩陣。其行數(shù)和列數(shù)均為電路的支路總數(shù)。二端口導納矩陣：對應y于二端口網絡方程，由二端口參數(shù)組成節(jié)點導納矩陣：以導納的形式描述電力網絡節(jié)點注入電流和節(jié)點電壓關系的矩陣。它給出了電力網絡連接關系和元件特性的全部信息，是潮流計算的基礎方程式。 本例應用結點導納矩陣具體計算時，根據(jù)如下公式：由題給出的導納可求的節(jié)點導納矩陣如下：進而節(jié)點導納矩陣為： 3.4潮流方程網絡方程是潮流計算的基礎

12、，如果給出電壓源或電流源，便可解得電流電壓分布。然而，潮流計算中，這些值都是無法準確給定的，這樣，就需要列出潮流方程。對n個節(jié)點的網絡，電力系統(tǒng)的潮流方程一般形式是 （i=1,2，n）其中,即PQ分別為節(jié)點的有功功率無功功率。3.5修正方程計算節(jié)點1的不平衡量計算節(jié)點2的不平衡量節(jié)點3是平衡節(jié)點，其電壓是給定的，故不參加迭代。根據(jù)給定的容許誤差，按收斂判據(jù)進行校驗，以上節(jié)點1、2的不平衡量都未滿足收斂條件，于是繼續(xù)以下計算。修正方程式為： (n=3)以上雅可比矩陣J中的各元素值是通過求偏導數(shù)獲得的，對PQ節(jié)點來說，是給定的，因而可以寫出 對PV節(jié)點來說，給定量是，因此可以列出當時, 雅可比矩陣

13、中非對角元素為當時,雅可比矩陣中對角元素為:代入數(shù)值后的修正方程為：求解修正方程得：3.6收斂條件一輪迭代結束，根據(jù)收斂條件收斂判據(jù)，若等式成立，結果收斂，迭代結束，計算平衡節(jié)點的功率和線路潮流計算，否則繼續(xù)計算雅可比矩陣，解修正方程，直到滿足收斂判據(jù)。4.程序設計4.1 節(jié)點導納矩陣的形成導納矩陣元素則表示為/*計算導納矩陣*G11=1.25;B11=-5.5;G22=1.3;B22=-7;G33=1.55;B33=-6.5;G12=G21=-0.5;B12=B21=3;G13=G31=-0.75;B13=B31=2.5;G23=G32=-0.8;B23=B32=4;for(i=1;i<

14、;4;i+)for(j=1;j<4;j+)printf("%f+(%f)j",Gij,Bij);printf(" ");printf("n");/形成節(jié)點導納矩陣/*printf("n");4.2 計算各節(jié)點不平衡量假定系統(tǒng)中的第1，2，3···，m號節(jié)點為PQ節(jié)點，第i個節(jié)點的給定功率設為和，對該節(jié)點可列寫方程 （i=1,2，···，m） 假定系統(tǒng)中的第m+1，m+2，···，n-1號節(jié)點為PV節(jié)點，則對其中每一個節(jié)點

15、可以列寫方程 （i=m+1，m+2，···，n-1） 第n號節(jié)點為平衡點，其電壓是給定的，故不參加迭代，其計算程序如下：/計算各節(jié)點不平衡量loop1:printf("迭代次數(shù)k1=%dn",k1);for (i=1;i<3;i+)float a=0,b=0;for(j=1;j<4;j+)a+=Gij*ej-Bij*fj; b+=Gij*fj+Bij*ej;Pi=Psi-(ei*a+fi*b);/計算有功功率的增量Qi=Qsi-(fi*a-ei*b);/計算無功功率的增量V22=V2s*V2s-e2*e2;printf("

16、有功功率增量P1=%f",P1); printf(" ,");printf("n");printf("有功功率增量P2=%f",P2); printf(" ,");printf("n");printf("無功功率增量Q1=%f",Q1); printf(" ,");printf("n");printf("電壓增量V22=%f",V22);printf("n")4.3 雅克比矩陣計算上述方

17、程中雅克比矩陣的各元素，可以對計算各點不平衡量得公式中求偏導數(shù)獲得。當時當時以下為程序：/*形成雅克比矩陣*for(j=1;j<3;j+)if(1=j)float c=0,d=0;int m;for(m=1;m<4;m+)c+=G1m*em-B1m*fm; d+=G1m*fm+B1m*em;J1*N-1j*N-1=-c-G1j*e1-B1j*f1;J1*N-1j*N=-d+B1j*e1-G1j*f1;J1*Nj*N-1=d+B1j*e1-G1j*f1;J1*Nj*N=-c+G1j*e1+B1j*f1;elseJ1*N-1j*N-1=-G1j*e1-B1j*f1; J1*Nj*N=G

18、1j*e1+B1j*f1; J1*N-1j*N=B1j*e1-G1j*f1; J1*Nj*N-1=B1j*e1-G1j*f1;for(j=1;j<3;j+)if(2=j)float c=0,d=0;int m;for(m=1;m<4;m+)c+=G2m*em-B2m*fm; d+=G2m*fm+B2m*em;J2*N-1j*N-1=-c-G2j*e2-B2j*f2;J2*N-1j*N=-d+B2j*e2-G2j*f2;J2*Nj*N-1=-2*e2;J2*Nj*N=-2*f2;elseJ2*N-1j*N-1=-G2j*e2-B2j*f2;J2*Nj*N=0;J2*N-1j*N=B2

19、j*e2-G2j*f2;J2*Nj*N-1=0;printf("雅克比矩陣是:n");for(i=1;i<5;i+)for(j=1;j<5;j+)printf("%f",Jij);printf(" ");printf("n");4.4 LU分解法求修正方程LU分解,又稱Gauss消去法,可把任意方陣分解成下三角矩陣的基本變換形式(行交換)和上三角矩陣的乘積。其數(shù)學表達式為:A=LU。其中L為下三角矩陣的基本變換形式,U為上三角矩陣。若有矩陣Ax=b 把矩陣LU分解，求AX=b的問題就等價于求出A=LU后

20、：因為Ly=b可求y，再因為Ux=y，可求出x。原始的求法x=A（-1）*b,某些情況下，如果矩陣A中的數(shù)非常小，我認為不是因為大數(shù)除以小數(shù)誤差大么，1/A算出的誤差會很大。但LU可以把A分解成兩個都比A大的矩陣的乘積，1/L的誤差比1/A小的多。求修正方程的程序如下/*計算修正方程*for(i=1;i<M;i+)Lii=1;for(i=1;i<M;i+)U1i=J1i;Li1=Ji1/U11;for(n=2;n<M;n+)for(j=n;j<M;j+)sigma1=0;for(s=0;s<=n-1;s+)sigma1+=Lns*Usj;Unj=Jnj-sigma

21、1;for(i=n;i<M;i+)sigma2=0;for(s=0;s<=n-1;s+)sigma2+=Lis*Usn;Lin=(Jin-sigma2)/Unn;b1=P1;b2=Q1;b3=P2;b4=V22;for(i=1;i<M;i+)sigma1=0;for(n=1;n<=i-1;n+)sigma1+=Lin*yn;yi=bi-sigma1;for(i=M-1;i>=1;i-)sigma2=0;for(n=i+1;n<M;n+)sigma2+=Uin*xn;xi=(yi-sigma2)/Uii;xe1=-x1;xe2=-x3; xf1=-x2;xf2

22、=-x4; printf("節(jié)點電壓:n");for(i=1;i<3;i+)ei+=xei; fi+=xfi;for(i=1;i<3;i+)printf("e%d=",i); printf("%f",ei); printf(" ,");for(i=1;i<3;i+)printf("f%d=",i); printf("%f",fi); printf(" ,");printf("n")4.5 計算網絡中功率分布最后要計算出

23、平衡節(jié)點的功率和網絡中的功率分布。5.結果分析給定節(jié)點電壓初值,經過五次迭代過程后，得到程序的顯示結果如下（?。赫砜傻霉?jié)點電壓和不平衡功率的變化情況，分別于表1和表2所示：迭代計數(shù)k節(jié)點電壓1121314151表1 迭代過程中節(jié)點電壓變化情況迭代計數(shù)k節(jié)點不平衡量0-2.000000-1.0000000.5000000.0000001-0.148189-0.976948-0.0726420.0000002-0.086588-0.595284-0.0495680.0096833-0.653416-4.467996-0.365523-0.0392754-0.191594-1.288891-0.

24、1011960.0022475-0.084434-0.585319-0.0496760.008383表2 迭代過程中節(jié)點不平衡量變化情況進行了五次迭代，結果仍然沒有收斂。 經過查找相關的資料得到：“多年的實踐證明，牛頓法具有很好的二次收斂性，是求解多元非線性方程的經典算法，至今仍是電力系統(tǒng)潮流計算的主流。因此，一般認為算法不是導致不收斂的原因，潮流不收斂產生的主要原因是計算的初始條件給得不合理，導致潮流方程無解?！?中國自動化網. 改善調度員潮流計算收斂性的措施6.小結通過本次電力系統(tǒng)分析課程設計，使我了解了自己在哪些方面有缺陷。首先，在拿到本次課設的題目時，就看不懂題目！這就給了自己一個不小

25、的打擊。于是我認真地看電力系統(tǒng)分析下冊有關潮流計算的牛頓-拉夫遜法！了解了何謂PQ節(jié)點，PV節(jié)點，平衡節(jié)點等。大體上知道了運用牛頓-拉夫遜法的各個步驟！其次，讀題的過程中也遇到了一些麻煩:1.不知道圖中節(jié)點二的“0.5”和“1”分別指代的是什么？2.各個節(jié)點對應的是什么節(jié)點？通過自己上網查資料，在網上看到了一些人有關的說法以及相似的題目的解答，這給我不小的啟發(fā)！我認識到了查找資料的重要性！經過了自己的努力，我知道了以上的答案：1.節(jié)點二中的“0.5”是有功功率，“1”代表的是電壓幅值。2.節(jié)點1是PQ節(jié)點，節(jié)點2是PV節(jié)點，節(jié)點3是平衡節(jié)點。最后，遇到了最難難的地方：用程序來實現(xiàn)牛頓-拉夫遜算

26、法！起初是想用matlab程序來進行程序的編寫并進行仿真，但是查找了一些資料以后，由于自己對于matlab的了解并不深而且自己學習的知識也不夠扎實，給自己帶來了不小的麻煩最終無法完成牛頓-拉夫遜的算法。后來在網上找到了有關C語言的相關程序！經過了自己的仔細研讀，了解了各個函數(shù)的作用。經過了自己的改造將程序完整的編輯了出來，并實現(xiàn)的預期的功能!通過本次課程設計，我知道了自己學習的知識還不夠扎實！很多方面只是應付考試，到了讓你做東西的時候確實還是相當困難的!尤其是在編程方面的缺陷！現(xiàn)如今是一個軟硬件相結合的時代，其中軟件更具有競爭性。因此在今后的學習過程中要端正學習態(tài)度。做好每一個細節(jié)，不斷完善自

27、我，提高自身的學習的水平。為將來的學習和工作打下良好的基礎！參考文獻1 何仰贊等.電力系統(tǒng)分析上冊M武漢：華中理工大學出版社.2 何仰贊等.電力系統(tǒng)分析下冊M武漢：華中理工大學出版社.3 諸俊偉等.電力系統(tǒng)分析M.北京：中國電力出版社，1995.4 周全仁等.電網計算與程序設計M.長沙：湖南科學技術出版社,1983.5丁化成.單片機應用技術A.北京:北京航空航天大學出版社,2000.附錄：#include <stdio.h>#include <math.h>#include <stdlib.h>#define N 2#define M 5main()doub

28、le G44,B44,J55;float e4=0,1,1,1,f4=0,P4,Q4,Ps4=0,-2,0.5,xe3,xf3;float Qs4=0,-1,V2s=1,float V22,max,P3,Q3;float a1=0,b1=0;int i,j,n,s,k1=0;float LMM=0,UMM=0,sigma1,sigma2,bM,yM,xM;/*計算導納矩陣*G11=1.25;B11=-5.5;G22=1.3;B22=-7;G33=1.55;B33=-6.5;G12=G21=-0.5;B12=B21=3;G13=G31=-0.75;B13=B31=2.5;G23=G32=-0.8

29、;B23=B32=4;for(i=1;i<4;i+)for(j=1;j<4;j+)printf("%f+(%f)j",Gij,Bij);printf(" ");printf("n");/形成節(jié)點導納矩陣/*printf("n");/*/計算各節(jié)點不平衡量loop1:printf("迭代次數(shù)k1=%dn",k1);for (i=1;i<3;i+)float a=0,b=0;for(j=1;j<4;j+)a+=Gij*ej-Bij*fj; b+=Gij*fj+Bij*ej;P

30、i=Psi-(ei*a+fi*b);/計算有功功率的增量Qi=Qsi-(fi*a-ei*b);/計算無功功率的增量V22=V2s*V2s-e2*e2;printf("有功功率增量P1=%f",P1); printf(" ,");printf("n");printf("有功功率增量P2=%f",P2); printf(" ,");printf("n");printf("無功功率增量Q1=%f",Q1); printf(" ,");prin

31、tf("n");printf("電壓增量V22=%f",V22);printf("n");/*篩選出最大值*max=fabs(P1)>fabs(P2)?fabs(P1):fabs(P2);max=max>fabs(Q1)?max:fabs(Q1);max=max>fabs(V22)?max:fabs(V22);printf("max=%fn",max);/*while (k1<=4)/*形成雅克比矩陣*for(j=1;j<3;j+)if(1=j)float c=0,d=0;int m;

32、for(m=1;m<4;m+)c+=G1m*em-B1m*fm; d+=G1m*fm+B1m*em;J1*N-1j*N-1=-c-G1j*e1-B1j*f1;J1*N-1j*N=-d+B1j*e1-G1j*f1;J1*Nj*N-1=d+B1j*e1-G1j*f1;J1*Nj*N=-c+G1j*e1+B1j*f1;elseJ1*N-1j*N-1=-G1j*e1-B1j*f1; J1*Nj*N=G1j*e1+B1j*f1; J1*N-1j*N=B1j*e1-G1j*f1; J1*Nj*N-1=B1j*e1-G1j*f1;for(j=1;j<3;j+)if(2=j)float c=0,d=0;int m;for(m=1;m<4;m+)c+=G2m*em-B2m*fm; d+=G2m*fm+B2m*em;J2*N-1j*N-1=-c-G2j*e2-B2j*f2;J2*N-1j*N=-d+B2j*e2-G2j*f2;J2*Nj*N-1=-2*e2;J2*Nj*N=-2*f2;elseJ2*N-1j*N-1=-G2j*e2-B2j*f2;J2*Nj*N=0