電力系統(tǒng)靜態(tài)安全分析

1、南京理工大學(xué)電力系統(tǒng)穩(wěn)態(tài)分析學(xué) 院：專 業(yè)：學(xué) 生 姓 名：題 目：日 期:指 導(dǎo) 教 師:自動化學(xué)院 電力系統(tǒng)及其自動化 陸未 學(xué) 號： 111101120 基于MATLAB的牛頓-拉夫遜 潮流計算研究 2011 年 12 月8 日 楊偉摘 要電力系統(tǒng)穩(wěn)態(tài)分析是研究電力系統(tǒng)運行和規(guī)劃方案最重要和最基本的手段，其任務(wù)是根據(jù)給定的發(fā)電運行方式及系統(tǒng)接線方式求解電力系統(tǒng)的穩(wěn)態(tài)運行情況，包括各母線定的電壓、各元件中通過的功率等等。在電力系統(tǒng)運行方式和規(guī)劃方案研究中，都需要通過潮流計算來分析比較運行方式或規(guī)劃方案的可行性、可靠性和經(jīng)濟(jì)性。牛頓-拉夫遜法是電力系統(tǒng)潮流計算的常用算法之一，它收斂性好，迭代

2、次數(shù)少，在電力系統(tǒng)潮流計算中得到了廣泛的應(yīng)用。本文介紹了電力系統(tǒng)潮流計算方法中的牛頓-拉夫遜法的相關(guān)知識及其基本原理，并用MATLAB編寫程序，最后通過算例來驗證該程序的正確性。關(guān)鍵詞：電力系統(tǒng)潮流計算；牛頓-拉夫遜法；MATLABI1 緒論電力系統(tǒng)穩(wěn)態(tài)分析是研究電力系統(tǒng)運行和規(guī)劃方案最重要和最基本的手段，其任務(wù)是根據(jù)給定的發(fā)電運行方式及系統(tǒng)接線方式求解電力系統(tǒng)的穩(wěn)態(tài)運行情況，包括各母線定的電壓、各元件中通過的功率等等。在電力系統(tǒng)運行方式和規(guī)劃方案研究中，都需要通過潮流計算來分析比較運行方式或規(guī)劃方案的可行性、可靠性和經(jīng)濟(jì)性。1.1潮流計算概述電力系統(tǒng)潮流計算也分為離線計算和在線計算兩種，前者

3、主要用于系統(tǒng)規(guī)劃設(shè)計和安排系統(tǒng)的運行方式，后者則用于正在運行系統(tǒng)的經(jīng)常監(jiān)視及實時控制。利用電子數(shù)字計算機(jī)進(jìn)行電力系統(tǒng)潮流計算從50年代中期就已經(jīng)開始。在這20年內(nèi)，潮流計算曾采用了各種不同的方法，這些方法的發(fā)展主要圍繞著對潮流計算的一些基本要求進(jìn)行的。對潮流計算的要求可以歸納為下面幾點：（1）計算方法的可靠性或收斂性；（2）對計算機(jī)內(nèi)存量的要求；（3）計算速度；（4）計算的方便性和靈活性。電力系統(tǒng)潮流計算問題在數(shù)學(xué)上是一組多元非線性方程式求解問題，其解法都離不開迭代。因此，對潮流計算方法，首先要求它能可靠地收斂，并給出正確答案。由于電力系統(tǒng)結(jié)構(gòu)及參數(shù)的一些特點，并且隨著電力系統(tǒng)不斷擴(kuò)大，潮流計

4、算的方程式階數(shù)也越來越高，對這樣的方程式并不是任何數(shù)學(xué)方法都能保證給出正確答案的。這種情況成為促使電力系統(tǒng)計算人員不斷尋求新的更可靠方法的重要因素。1.2 潮流計算的意義及其發(fā)展電力系統(tǒng)潮流計算是電力系統(tǒng)分析中的一種最基本的計算，是對復(fù)雜電力系統(tǒng)正常和故障條件下穩(wěn)態(tài)運行狀態(tài)的計算。潮流計算的目標(biāo)是求取電力系統(tǒng)在給定運行狀態(tài)的計算。即節(jié)點電壓和功率分布，用以檢查系統(tǒng)各元件是否過負(fù)荷。各點電壓是否滿足要求，功率的分布和分配是否合理以及功率損耗等。對現(xiàn)有電力系統(tǒng)的運行和擴(kuò)建，對新的電力系統(tǒng)進(jìn)行規(guī)劃設(shè)計以及對電力系統(tǒng)進(jìn)行靜態(tài)和 1暫態(tài)穩(wěn)定分析都是以潮流計算為基礎(chǔ)。潮流計算結(jié)果可用如電力系統(tǒng)穩(wěn)態(tài)研究，安

5、全估計或最優(yōu)潮流等對潮流計算的模型和方法有直接影響。實際電力系統(tǒng)的潮流技術(shù)那主要采用牛頓-拉夫遜法。在運行方式管理中，潮流是確定電網(wǎng)運行方式的基本出發(fā)點；在規(guī)劃領(lǐng)域，需要進(jìn)行潮流分析驗證規(guī)劃方案的合理性；在實時運行環(huán)境，調(diào)度員潮流提供了多個在預(yù)想操作情況下電網(wǎng)的潮流分布以校驗運行可靠性。在電力系統(tǒng)調(diào)度運行的多個領(lǐng)域都涉及到電網(wǎng)潮流計算。潮流是確定電力網(wǎng)絡(luò)運行狀態(tài)的基本因素，潮流問題是研究電力系統(tǒng)穩(wěn)態(tài)問題的基礎(chǔ)和前提。在用數(shù)字計算機(jī)解電力系統(tǒng)潮流問題的開始階段，普遍采取以節(jié)點導(dǎo)納矩陣為基礎(chǔ)的逐次代入法。這個方法的原理比較簡單，要求的數(shù)字計算機(jī)內(nèi)存量比較下，適應(yīng)50年代電子計算機(jī)制造水平和當(dāng)時電力

6、系統(tǒng)理論水平。但它的收斂性較差，當(dāng)系統(tǒng)規(guī)模變大時，迭代次數(shù)急劇上升，在計算中往往出現(xiàn)迭代不收斂的情況。這就迫使電力系統(tǒng)計算人員轉(zhuǎn)向以阻抗矩陣為基礎(chǔ)的逐次代入法。阻抗法改善了系統(tǒng)潮流計算問題的收斂性，解決了導(dǎo)納法無法求解的一些系統(tǒng)的潮流計算，在60年代獲得了廣泛的應(yīng)用。阻抗法的主要缺點是占用計算機(jī)內(nèi)存大，每次迭代的計算量大。當(dāng)系統(tǒng)不斷擴(kuò)大時，這些缺點就更加突出。為了克服阻抗法在內(nèi)存和速度方面的缺點，60年代中期發(fā)展了以阻抗矩陣為基礎(chǔ)的分塊阻抗法。這個方法把一個大系統(tǒng)分割為幾個小的地區(qū)系統(tǒng)，在計算機(jī)內(nèi)只需要存儲各個地區(qū)系統(tǒng)的阻抗矩陣及它們之間聯(lián)絡(luò)線的阻抗，這樣不僅大幅度地節(jié)省了內(nèi)存容量，同時也提高

7、了計算速度?？朔杩狗ㄈ秉c的另一途徑是采用牛頓-拉夫遜法。這是數(shù)學(xué)中解決非線性方程式的典型方法，有較好的收斂性。在解決電力系統(tǒng)潮流計算問題時，是以導(dǎo)納矩陣為基礎(chǔ)的，因此，只要我們能在迭代過程中盡可能保持方程式系數(shù)矩陣的稀疏性，就可以大大提高牛頓法潮流程序的效率。自從60年代中期，在牛頓法中利用了最佳順序消去法以后，牛頓法在收斂性。內(nèi)存要求。速度方面都超過了阻抗法，成為60年代末期以后廣泛采用的優(yōu)秀方法。1.3 本文的主要工作本文介紹了電力系統(tǒng)潮流計算方法中的牛頓-拉夫遜法的相關(guān)知識及其基本原理，并用MATLAB編寫程序，最后通過算例來驗證該程序的正確性。 22 牛頓-拉夫遜法潮流計算基本理論牛

8、頓法（又稱牛頓-拉夫遜法）是求解非線性方程式的典型方法。該方法有較好的收斂性，迭代次數(shù)少，在電力系統(tǒng)潮流計算中也得到應(yīng)用。目前，牛頓法潮流計算是最為廣泛、效果最好的一種潮流計算方法。2.1 牛頓法基本原理牛頓法是把非線性方程式的求解過程變成反復(fù)對相應(yīng)的線性方程式的求解過程，即非線性問題通過線性化逐步近似，這就是牛頓法的核心。下面以非線性方程式的求解過程來進(jìn)行說明。設(shè)有非線性方程式f(x)=0 （1.1）設(shè)x(0)為該方程式的初值， Dx(0)為初值x(0)的修正量。如果求得Dx(0)，則就可以得到真解x=x(0)-Dx(0) （1.2）為此，將式f(x(0)-Dx(0)=0 （1.3）按泰勒級

9、數(shù)展開f(x(0)-Dx)=f(x)-f(x)Dx(0)(0)'(0)(0)f''(x(0)+(Dx)2-L （1.4） 2!如果x(0)接近真值，則Dx(0)相對來說是足夠小，所以可以略去所有Dx(0)的高次項。因此f(x(0)-f'(x(0)Dx(0)=0 （1.5）可得 Dx(0)f(x(0)='(0) （1.6） f(x)由于式（1.5）是式（1.4）的簡化結(jié)果，所以由式（1.6）得到Dx(0)后，還不能得到方程式（1.1）的真解。實際上，用Dx(0)對x(0)修正后得到的x(1)：x(1)=x(0)-Dx(0) （1.7）只是向真正解更逼近一些

10、?，F(xiàn)在如果再以作為初值x(1)，求解3f(x(1)-f'(x(1)Dx(1)=0 （1.8） 就能得到更趨近真正解的x(2)x(2)=x(1)-Dx(1) （1.9）這樣反復(fù)下去，就構(gòu)成了不斷求解非線性方程式的逐次線性化過程。第t次迭代時的參數(shù)方程為f(x(t)-f'(x(t)Dx(t)=0 （1.10）或 f(x(t)=f'(x(t)Dx(t) （1.11）上式左端可以看成是近似解x(t)引起的誤差，當(dāng)f(x(t)®0時，就滿足了原方程式（1.1），因而x(t)就成為該方程的解。式中f'(x(t)是函數(shù)f(x)=0 在x(t)點的一次導(dǎo)數(shù)，也就是曲線

11、在x(t)點的斜率，如圖1.1所示，tga(t)=f'(x(t) （1.12）修正量Dx(t)則是由x(t)點的切線與橫軸的交點來確定，由圖1.1可以直觀的看出牛頓法的求解過程。圖1.1 牛頓法的幾何解釋現(xiàn)在把牛頓法推廣到多變量非線性方程組的情況。設(shè)有變量x1,x2,L,xn的非線性聯(lián)立方程組：4f1(x1,x2,L,xn)=0üf2(x1,x2,L,xn)=0ïïý (1.13)Mïfn(x1,x2,L,xn)=0ïþ給定各變量初值x1,x2,L,xn并使其滿足(0)(0)f1(x1(0)-Dx1(0),x2-Dx

12、2,L,xn(0)f2(x1(0)-Dx1(0),x2(0)-Dxn)=0üï(0)(0)(0)-Dx2,L,xn-Dxn)=0ïý （1.14） ï(0)(0)(0)-Dx2,L,xn-Dxn)=0ïþ(0)(0)(0)(0)(0)(0)，假設(shè)Dx1,Dx2,L,Dxn(0)(0)(0)為其修正量，(0)M(0)fn(x1(0)-Dx1(0),x2對以上n個方程式分別按泰勒級數(shù)展開，當(dāng)忽略Dx1,Dx2,L,Dxn的二次項和高次項時，可以得到所組成üé¶f1ù¶f1

13、2;f1(0)(0)(0)f1(x,x,L,x)-êDx1+Dx2+LDxnú=0ï¶x¶x¶xïêú20n0ë10ûïéù¶f¶f¶fï(0)(0)(0)(0)f2(x1(0),x2,L,xn)-ê2Dx1(0)+2Dx2+L2Dxnú=0ï¶x20¶xn0ý （1.15）êúë¶x10ûïM&#

14、239;ïé¶fnù¶fn¶fn(0)(0)(0)(0)(0)(0)fn(x1,x2,L,xn)-êDx1+Dx2+LDxnú=0ïï¶x20¶xn0êúë¶x10ûþ(0)1(0)2(0)n式中：(0)¶fi¶xi(0)為函數(shù)fi(x1,x2,L,xn)對自變量xj的偏導(dǎo)數(shù)在點（Dx1,Dx2,L,Dxn(0)）處的值。把上式寫成矩陣形式：é¶f1êê&#

15、182;x1(0)(0)(0)éf1(x1,x2,L,xn)ùêêúê¶f2(0)(0)(0)êf2(x1,x2,L,xn)ú=ê¶xêMúê1êúêM(0)(0)(0)êëfn(x1,x2,L,xn)úûêê¶fnêë¶x1¶f1¶fùL1ú¶x20¶xn0ú

16、;0(0)éùDx1ú¶f2¶fê(0)úL2úêDx2ú （1.16） ¶x20¶xn0úê0úúMúê(0)úëDxnúûúê¶fn¶fnúL¶x¶xn0ú200û5這是變量Dx1,Dx2,L,Dxn過它可以解出Dx1,Dx2,L,Dxn(0)(0)(0)(0)(0)的線性方程組，稱為

17、牛頓法的修正方程，通，并可以進(jìn)一步求得(0)x1(1)=x1(0)-Dx1(0)üï(1)(0)(0)x2=x2-Dx2ïý (1.17)Mï(1)(0)(0)ïxn=xn-Dxnþ式中x1,x2,L,xn向真正解逼近了一步，如果再以它們作為初值重復(fù)解式（1.16）修正方程式，等到更接近真解的x1,x2,L,xn，如此迭代下去，并按式（1.17）進(jìn)行修正，直到滿足收斂要求為止并停止迭代計算，這就構(gòu)成了牛頓法的迭代過程。一般第t次迭代式的修正方程為é¶f1êê¶x1(t)(t)

18、(t)éf1(x1,x2,L,xn)ùêêúê¶f2(t)(t)(t)f(x,x,L,x)ê21ú=ê¶x2nêMúê1êúêM(t)(t)(t)êëfn(x1,x2,L,xn)úûêê¶fnêë¶x1¶f1¶f1ùLú¶x¶x2ntttú(t)ú&

19、#233;Dx1ù¶f2¶fê(t)úL2úêDx2ú (1.18) ú¶x¶x2ntttêúúMúê(t)úëDxnúûúê¶fn¶fnúL¶x2t¶xntútû(2)(2)(2)(1)(1)(1)上式可以簡寫為F(X(t)=J(t)DX(t) (1.19)式中(t)(t)éf1(x1(t),x2

20、,L,xn)ùêú(t)(t)(t)f(x,x,L,x)ú (1.20) 212nF(X(t)=êêMúêú(t)(t)(t)êëfn(x1,x2,L,xn)úû為第t次迭代時函數(shù)的誤差相量；6J(t)é¶f1êê¶x1êê¶f2=ê¶x1êêMêê¶fnêë¶x1¶f1

21、2;f1ùLú¶x¶x2tntútú¶f2¶f2úL¶x2t¶xntú (1.21) túúú¶fn¶fnúL¶x¶xntú2ttû稱為的第t次迭代時的雅克比矩陣；éDx1(t)ùê(t)úDx2ú (1.22) =êêMúêú(t)êëDxnú&#

22、251;DX(t)為第t次迭代時的修正量相量。同樣，也可以寫出類似（1.17）的算式X(t+1)=X(t)-DX(t) (1.23)這樣反復(fù)交替的解式（1.19）及式（1.23）就可以使X(t+1)逐步趨近方程式的真正解。為了判斷收斂情況，可采用一下兩個不等式中的一個：DX(t)<e1 (1.24)F(X(t)<e2 (1.25)式中，e1,e2為預(yù)先給定的很小正數(shù)。2.2 牛頓法潮流求解在潮流計算中，節(jié)點極坐標(biāo)功率方程式Pi=ViåVj(Gijcosqij+Bijsinqij) (1.26)jÎiQi=ViåVj(Gijcosqij-Bijsinqi

23、j) (1.27)jÎi式中：jÎi表示å號后的節(jié)點j都直接與i節(jié)點相連，并且包括j=i的情況。節(jié)點功率誤差DPi=Pis-ViåVj(Gijcosqij+Bijsinqij) (1.28)jÎi7DQi=Qis-ViåVj(Gijcosqij-Bijsinqij) (1.29)jÎi式中：Pis，Qis為節(jié)點i給定的有功功率及無功功率。當(dāng)采用的是極坐標(biāo)模型時，待求量是各節(jié)點電壓的幅值Vi和角度qi。對于PV節(jié)點來說，節(jié)點i電壓幅值Vi是給定的，不再作為變量。同時，該點不能預(yù)先給定無功功率Qis，這樣，方程式中DQi也就失去

24、了約束作用。因此，在迭代過程中應(yīng)該取消與PV節(jié)點有關(guān)的無功功率方程式。只有當(dāng)?shù)Y(jié)束后，即各節(jié)點電壓向量求得之后，才利用這些方程式來求各PV節(jié)點應(yīng)維持的無功功率。同理，由于平衡節(jié)點電壓幅值和相角都是給定量，因此與平衡節(jié)點有關(guān)的方程式也不參與迭代過程。迭代結(jié)束后，利用平衡節(jié)點的功率方程式來確定其有功功率和無功功率。設(shè)系統(tǒng)節(jié)點總數(shù)為n，PV節(jié)點共r個，PQ節(jié)點共 n-r-1個，平衡節(jié)點1個（在潮流計算中，該類節(jié)點一般在系統(tǒng)中只設(shè)一個）。為了求解方便，把平衡節(jié)點排在最后，即設(shè)為第n個節(jié)點，則潮流要解的方程式應(yīng)包括 DPi=Pis-ViåVj(Gijcosqij+Bijsinqij)， i=

25、1,2,L,n-1，共n-1個方程式。jÎi Qi=ViåVj(Gijcosqij-Bijsinqij)，i=1,2,L,n-r-1，共n-r-1個方程式。jÎi以上方程式的待求量為各節(jié)點電壓的角度qi及電壓幅值Vi，其中qi共有n-1個。由于Vi中不包括PV節(jié)點的電壓幅值，所以共有n-r-1個。這樣未知量共有2n-r-2個，恰好可由以上2n-r-2個方程式求出。將式(1.28)、式(1.29)按泰勒級數(shù)展開，略去高次項，可得éDPùéH111êDPúêH212êúê

26、4;MúêMêúêHn-1,1DPn-1úê=êêDQ1úêK11êúêDQêK21ê2úêMúêMêúêDQên-r-1ûúêëëKn-r-1,1H12H22MHn-1,2K12K21MLLMLLLMH1,n-1H2,n-1MHn-1,n-1K1,n-1K2,n-1MN11N21MNn-1,1L11L21M

27、Ln-r-1,1N12N22MNn-1,2L12L22MLLMLLLMN1,n-1ùéDq1ùêúN2,n-1úDq2úêúúêúMMúêúNn-1,n-1úêDqn-1úL1,n-1úêDV1/V1úúêúL2,n-1úêDV2/V2úúêúMMúêúLn-r-1,n

28、-1úDV/Vên-r-1n-r-1ûúûëKn-r-1,2LKn-r-1,n-1Ln-r-1,2L(1.30)8éDPùéH或 êú=êëDQûëKNùéDqúêDVL/ûëù(1.31) úVû式中，電壓幅值的修正量采用DV1/V1,DV2/V2,LDVn-r-1/Vn-r-1的形式?jīng)]有特殊意義，只是為了使雅克比矩陣中各元素具有比較相似的表達(dá)式。利用簡單

29、的微分運算對式（1.28）、式（1.29）取偏導(dǎo)，可以得到雅克比矩陣中各元素的表達(dá)式 非對角元素（j¹i）Hij=¶DPi=-VVij(Gijsinqij-Bijcosqij) (1.32) ¶qi¶DPiVj=-VVij(Gijcosqij+Bijsinqij) (1.33) ¶Vj¶DPiVj=VVij(Gijcosqij+Bijsinqij) (1.34) ¶Vj¶DQiVj=-VVij(Gijsinqij-Bijcosqij) (1.35) ¶VjNij=Kij=Lij= 對角元素Hii=

30、2;DPi=ViåVj(Gijsinqij-Bijcosqij)=Vi2Bii+Qi (1.36) ¶qijÎij¹iNii=¶DPiVi=-ViåVj(Gijcosqij+Bijsinqij)-2Vi2Gii=-Vi2Gii-Pi (1.37) ¶VijÎij¹iKii=¶DPiVi=-ViåVj(Gijcosqij+Bijsinqij)=Vi2Gii-Pi (1.38) ¶VijÎij¹iLii=¶DQiVi=-ViåVj(Gijs

31、inqij-Bijcosqij)+2Vi2Bii=Vi2Bii-Qi (1.39) ¶VijÎij¹i2.3 牛頓法潮流計算步驟下面討論的是用極坐標(biāo)形式的牛頓法求解過程，大致分為以下幾個步驟： 形成節(jié)點導(dǎo)納矩陣； 給各節(jié)點電壓設(shè)初值（Vi(0),qi(0)）； 將節(jié)點電壓初值代入式(1.28)、式(1.29)，求出修正方程式的常數(shù)項向量9DPi,DQi； 將節(jié)點電壓初值代入式(1.32)式(1.39)，求出雅可比矩陣元素； 求解修正方程式(1.30)，求出變量的修正向量DV,Dq； 根據(jù)式(1.17)或式(1.23)，求取節(jié)點電壓的新值；圖1.2 牛頓法計算潮流的

32、程序框圖 檢查是否收斂，由式(1.24)和式(1.25)可知，若電壓趨近于真解時，功率偏移量將趨于零.如不收斂，則以各節(jié)點電壓的新值作為初值自第步重新開始下一次迭代，否則轉(zhuǎn)入下一步。 計算支路功率分布，PV節(jié)點無功功率和平衡節(jié)點注入功率。 輸出結(jié)果，并結(jié)束。10圖1.2為牛頓法計算潮流的程序框圖。圖中Kmax為事先給定的最大迭代次數(shù)，當(dāng)實際迭代次數(shù)K>Kmax時，即認(rèn)為計算不收斂。下面就據(jù)此流程圖編寫程序。113 牛頓-拉夫遜法潮流計算編程本文討論用極坐標(biāo)形式的牛頓法潮流的求解。下面就幾個子函數(shù)作簡要說明。3.1 潮流計算數(shù)據(jù)文件建立在一般的潮流計算中，需輸入的原始數(shù)據(jù)包括節(jié)點（母線）數(shù)

33、據(jù)；支路（包括接地支路，變壓器支路等）數(shù)據(jù)。3.1.1 節(jié)點數(shù)據(jù)節(jié)點數(shù)據(jù)按一下的格式進(jìn)行輸入格式：節(jié)點編號，節(jié)點電壓，節(jié)點相角，節(jié)點注入有功，節(jié)點注入無功，節(jié)點類型其中，節(jié)點類型取值為1、2、3，分別對應(yīng)PQ節(jié)點、PV節(jié)點和平衡節(jié)點；節(jié)點注入功率如果是負(fù)荷，取負(fù)值；否則取正值；節(jié)點編號的原則是先PQ節(jié)點，然后PV節(jié)點，最后平衡節(jié)點。3.1.2 支路數(shù)據(jù)對于支路(包括對地支路和變壓器），原始數(shù)據(jù)格式定義如下格式：節(jié)點I，節(jié)點J,，線路電阻R，線路電抗X，線路電導(dǎo)G，線路電納B，變壓器變比K其中，K=0表示普通線路；K>0，表示變壓器非標(biāo)準(zhǔn)變比在j側(cè)；K<0，表示變壓器非標(biāo)準(zhǔn)變比在i側(cè)

34、。于是，可以通過以上的格式進(jìn)行節(jié)點和線路數(shù)據(jù)輸入。在實際中，我們通過把節(jié)點和線路數(shù)據(jù)存在文件中，通過子程序1將數(shù)據(jù)文件的打開進(jìn)行數(shù)據(jù)輸入，程序代碼詳見附錄A。3.2 節(jié)點重新排序在前面講到的電力系統(tǒng)的潮流計算中，總是假設(shè)節(jié)點編號依據(jù)PQ節(jié)點、PV節(jié)點、平衡節(jié)點的次序進(jìn)行編號的。但是，實際電力系統(tǒng)的節(jié)點編號具有一定的隨意性，并不總是按照這個原則形成原始數(shù)據(jù)的。因此，當(dāng)原始節(jié)點數(shù)據(jù)不按照PQ節(jié)點、PV節(jié)點、平衡節(jié)點的次序進(jìn)行編號時，需要對節(jié)點重新編號。該程序見子程序2，程序代碼詳見附錄A。 12在潮流計算完成后還需要重新將計算結(jié)果重新還原為原來的節(jié)點編號。該程序見子程序9，程序代碼詳見附錄A。3.

35、3 計算節(jié)點導(dǎo)納矩陣由電力系統(tǒng)分析的知識可知，節(jié)點導(dǎo)納矩陣的計算可以歸納如下：1) 節(jié)點導(dǎo)納矩陣是方陣，其階數(shù)就等于網(wǎng)絡(luò)中除去參考節(jié)點外的節(jié)點數(shù)。參考節(jié)點一般取大地，編號為零。2) 節(jié)點導(dǎo)納矩陣是稀疏矩陣，其各行非對角元素中非零元素的個數(shù)等于對應(yīng)節(jié)點所連接的不接地支路數(shù)。3) 節(jié)點導(dǎo)納矩陣的對角元素，即各該節(jié)點的自導(dǎo)納等于相應(yīng)節(jié)點所連支路的導(dǎo)納之和。Yii=åyijjÎi (3.1)式中，yij為節(jié)點i與節(jié)點j間支路阻抗的倒數(shù)；當(dāng)節(jié)點i有接地支路時，還應(yīng)包括j=0的情況。4) 節(jié)點導(dǎo)納矩陣的非對角元素Yij等于節(jié)點i與節(jié)點j之間的導(dǎo)納的負(fù)值。Yij=-1=-yij zij

36、(3.2)5) 節(jié)點導(dǎo)納矩陣一般是對稱矩陣，這是網(wǎng)絡(luò)的互易特性所決定的。從而，一般只要求求取這個矩陣的上三角或下三角部分。按照以上的算式，對于實際網(wǎng)絡(luò)可以根據(jù)給定的支路參數(shù)和連接情況，直觀而簡單地求出導(dǎo)納矩陣。對于電力系統(tǒng)中的支路參數(shù)，可以用一下的算式來表示Zt=R+jX或Yt=1/ZtYm=G+jB/2 (3.3)設(shè)導(dǎo)納矩陣為Y，開始時設(shè)Y=0。3.3.1 普通線路對于普通的線路來說，K=0,j¹0，等值電路如圖3.1所示。由圖可知 Y(I,I)=Y(I,I)+Yt+Ym；Y(J,J)=Y(J,J)+Yt+Ym；13Y(I,J)=Y(I,J)-Yt； Y(J,I)=Y(I,J)；y

37、t圖3.1 普通線路等值電路3.3.2對地支路對于對地支路來說，K=0,j=0，等值電路如圖3.2所示。由圖可知 Y(I,I)=Y(I,I)+Ym；圖2.2 對地支路等值電路3.3.3變壓器支路對于對地支路來說，因非標(biāo)準(zhǔn)變比究竟在哪一側(cè)的不同，有下面兩種情況： 1）變壓器支路（k>0，非標(biāo)準(zhǔn)變比在j側(cè)），等值電路如圖3.3所示。由圖可知Y(I,I)=Y(I,I)+Yt+Ym； Y(J,J)=Y(J,J)+K*K*Yt； Y(I,J)=Y(I,J)-K*Yt； Y(J,I)=Y(I,J)；yt/k*圖2.3 變壓器等值電路（k>0，非標(biāo)準(zhǔn)變比在j側(cè)）2）變壓器支路（k<0，非標(biāo)

38、準(zhǔn)變比在i側(cè)），等值電路如圖3.4所示。由圖可知14Y(I,I)=Y(I,I)+Yt+Ym;Y(J,J)=Y(J,J)+K*K*Yt;Y(I,J)=Y(I,J)-K*Yt;Y(J,I)=Y(I,J);k*ytk*)yt圖2.4 變壓器等值電路（k<0，非標(biāo)準(zhǔn)變比在i側(cè)）由上述的基本原理和計算公式來編寫程序，可以得到形成節(jié)點導(dǎo)納矩陣的子程序3，程序代碼詳見附錄A。3.4 功率偏差的計算按式(1.28)和式(1.29)進(jìn)行計算，就可很容易地得到功率偏差。該程序見子程序3，程序代碼詳見附錄A。3.5 計算雅克比矩陣按式(1.32)和式(1.39)進(jìn)行計算，就可得到雅克比矩陣。該程序見子程序4，

39、程序代碼詳見附錄A。3.6 列主元消去法求解方程組在求解式(1.30)或式(1.31)時，運用列主元消去法進(jìn)行求解。對于求解線性方程組Ax=b等價于求解PA=LUìLy=Pb íîUx=y (3.4)由上述的基本原理和計算公式來編寫程序，可以得到求解方程組的子程序5，程序代碼詳見附錄A。3.7 計算節(jié)點功率注入按式(1.26)和式(1.27)進(jìn)行計算，就可很容易地得到PV節(jié)點的無功和平衡節(jié)點的功率。該程序見子程序7，程序代碼詳見附錄A。153.8 計算支路潮流潮流計算的目的不僅僅是為了得到節(jié)點信息，更重要的是希望得到系統(tǒng)中每條線路的功率流動情況。計算每條線路流過的

40、功率的計算公式如下： 3.8.1 普通線路由普通線路的等值電路圖3.1，可得ggggægöIij=çUi-Uj÷YT+UiYm=Ui(YT+Ym)-UjYTèø(3.5) (3.6) (3.7)Sij=UiIij=U(YT+Y2i*gg*2j*gg*m)-UUigg*jYT*Sji=UjIji=U(YT+Y3.8.2 變壓器支路*m)-UgjUiYT*g*對于變壓器支路等值電路如圖3.3和圖3.4所示，它的計算同普通線路。 3.8.3 對地支路由對地支路的等值電路圖3.2，可得Ii0=UiYmgg*g(3.8) (3.9)Si0=Ui

41、Ii0=Ui2Ym*由上述的基本原理和計算公式來編寫程序，可以得到計算支路潮流的子程序8，程序代碼詳見附錄A。3.9 牛頓-拉夫遜法潮流計算程序根據(jù)圖1.2 所示的牛頓法計算潮流的程序框圖編寫程序，可以得到牛頓法計算潮流的主程序，程序代碼詳見附錄A。此外，潮流結(jié)果的輸出。該程序見子程序10，程序代碼詳見附錄A。164 算例4.1系統(tǒng)模型電力系統(tǒng)模型如圖4.1所示，系統(tǒng)的參數(shù)見圖。利用該模型進(jìn)行牛頓-拉夫遜潮流計算。圖4.1 電力系統(tǒng)模型在程序求解潮流時，數(shù)據(jù)的輸入采用M文件的輸入方式。節(jié)點和支路數(shù)據(jù)保存在lw_system.m中，如圖4.2所示。圖4.2 原始節(jié)點和支路數(shù)據(jù)4.2 計算結(jié)果17在MATLAB的Command Window中輸