隱式改進歐拉法在暫態穩定分析中的應用

蔣德瓏 尹淑萍 王克文 代旭光

(鄭州大學電氣工程學院，河南 鄭州 450001)

0 引言

經濟和技術的進步使電力系統得到了迅速發展，但也對電力系統穩定性提出了更高的要求［1］。多年來，大量文獻研究了電力系統的專業化和自動化，在電力系統暫態穩定自動化分析中更是取得了長足的進步［2-11］。如采用計算出擾動后各發電機轉子間相對角度隨時間的變化曲線(即搖擺曲線)來判斷系統的穩定性［6-9］;對暫態穩定程序進行模塊化處理，使程序具有較好的擴展性和靈活性［10-11］。

相對顯式解法，隱式解法有其明顯的優越性。在總結前人研究的基礎上，本文將經典的改進歐拉法與隱式解法相結合，用于電力系統的暫態穩定計算，以有效地改善計算的精度和速度，提高暫態穩定自動化分析的水平，為電力系統的安全穩定運行提供更為有力的保障。

1 牛頓法潮流計算

牛頓法是數學中解決非線性方程式的典型方法，它通過對非線性方程逐次線性化求解方程的根，并能保持較好的收斂性。潮流計算是研究電力系統暫態穩定的基礎，它根據給定的發電運行方式和系統接線方式來確定整個電力系統各部分的運行狀態，包括各母線的電壓、各元件流過的功率和系統的功率損耗等，為暫態穩定分析提供初始運行狀態。當前，牛頓法被廣泛應用于電力系統潮流計算，其計算的核心問題是修正方程式的建立和求解［12］。本文采用經典的牛頓法潮流計算方法，其求解過程大致分為以下幾步:首先給定節點導納矩陣YB和節點電壓初值e(0)、f(0)，代入功率偏差量的表達式中，求解 ΔP(0)、ΔQ(0)、(ΔV2)(0)，再求出修正方程式的系數矩陣(雅可比矩陣)各元素;接著通過解修正方程式得到修正量Δe(0)、Δf(0)，然后對節點電壓 e(1)=e(0)-Δe(0)、f(1)=f(0)-Δf(0)進行修正，進而利用 e(1)、f(1)求得 ΔP(1)、ΔQ(1)、(ΔV2)(1);最后判斷結果是否收斂，如收斂，則求出各支路潮流并輸出結果，否則以e(1)、f(1)為初值重新迭代。

2 穩定分析數學模型的確定

整個電力系統的模型，在數學上可以統一描述成如下一般形式的微分-代數方程組:式中:x為微分方程組中描述系統動態特性的狀態變量;y為代數方程組中系統的運行參量。本文重點關注發電機微分方程的求解。

2.1 發電機數學模型

考慮到暫態穩定的計算精度和系統的復雜性，本文采用忽略定子電磁暫態但考慮轉子阻尼繞組作用的五階發電機模型(E'q、E″q、E'd、δ、ω)，即考慮 f繞組、D繞組、Q繞組的電磁暫態以及轉子運動的機電暫態。

轉子電壓平衡方程為:

轉子運動方程為:

式中:T'd為勵磁繞組本身的時間常數;T″d為勵磁繞組短路、定子繞組開路時d軸阻尼繞組的時間常數;T″q為定子繞組開路時q軸阻尼繞組的時間常數;TJ為發電機組的慣性時間常數;E'q為發電機q軸的暫態電動勢;E″q為發電機q軸的次暫態電動勢;E″d為發電機d軸的次暫態電動勢;δ為功角;ω為轉子角頻率。將式(2)與式(3)聯立，即可得到同步發電機的五階數學模型。

2.2 隱式改進歐拉法

對于一般的非線性微分方程式:

確定時間步長 Δt=h，由(tn，xn)點推算(tn+1，xn+1)點時，改進歐拉法的遞推公式即為:

改進歐拉法的核心思想是取某時段始點導數值與終點導數值的平均值作為該時刻折線段的斜率，進而得到更加精確的計算結果。而用f(xn，tn)和f(xn+1，tn+1)代替積分曲線［tn，tn+1］區間始點和終點的斜率，即可將式(5)進一步變換為:

式(6)中只有xn+1為未知數，因而利用求解代數方程式的方法就可以求得xn+1，這就得到了隱式改進歐拉法。

傳統的改進歐拉法是利用已知量，根據遞推公式逐次遞推計算出相應時段終點的函數值，因此是一種顯式解法。隱式改進歐拉法不是給出遞推公式，而是首先把微分方程化為差分方程，然后利用求解差分方程的方法確定函數值，其特點就是把微分方程的求解問題轉換成一系列代數方程的求解過程。因此，隱式解法有兩個明顯的優點:一是當微分方程和代數方程需要聯立求解時，利用隱式解法便于消除交接誤差;二是隱式解法可以采取較大的步長。

2.3 差分方程

對于電力系統暫態計算涉及到的發電機、電動機和控制器等數學模型中的微分方程，皆可以將其化為差分方程［13-14］。下面以同步發電機的數學模型為例，利用隱式改進歐拉法，將五階發電機模型中的微分方程化為差分方程。

由此得到的用于描述各發電機、電動機、控制器等差分方程組，可與用于描述電力網絡的代數方程組聯立求解，發揮了牛頓法所具有的數值穩定、無交接誤差和收斂速度快等特點。

3 初值計算及運行參數的求解

根據上述潮流計算法，可以計算出電力系統遭受擾動以前的正常運行狀態，得到電力網絡的運行參數，并由此求出電力系統中有關元件的狀態參數初始值。這些都是進行電力系統暫態穩定分析的先決條件。

3.1 初值的計算

在進行暫態穩定分析之前，必須求出暫態過程計算所需要的初值條件，即首先要確定求解微分方程所需要的初值。

在簡化的暫態穩定計算程序中，初值的計算包括求系統擾動瞬間發電機的暫態電勢、轉子角度、原動機的機械功率和等值導納等。這些參數在電力系統擾動的瞬間是不會突變的［10］。因此，可以由擾動前的正常運行狀態計算得到。

3.2 與網絡方程聯立求解

在電力系統中，電力網絡將系統中看起來相互獨立的所有暫態元件聯系在一起，在暫態過程中的任一時刻，各暫態元件注入網絡的電流不但由其本身的特性決定，且整個電力網絡必須滿足基爾霍夫定律。其中，前者由各暫態元件自身的代數方程描述，后者反映在電力網絡方程中。因此，為了求解網絡方程，需要列出各暫態元件自身的代數方程，并對其進行處理，從而與網絡方程聯立求解。將發電機作為電流源，其節點注入電流的表達式為:

對于負荷模型，只需對導納矩陣的相應對角塊進行變化，即可得到新的網絡方程。

4 算例分析

利用上述方法編寫暫態穩定計算程序，對典型的9節點三機系統進行仿真計算，三機系統示意圖如圖1所示。

圖1 三機系統示意圖Fig.1 Three-machine system

圖1中，系統有3臺發電機、3個負荷和9條支路，頻率為60 Hz，支路數據和發電機參數見文獻［14］。

假設系統在0 s受到擾動，在線路5-7靠近母線7處發生三相接地短路，則故障在5個周波(約0.08333 s)后由斷開線路5-7而被消除。

本文采用牛頓法潮流計算得出正常運行情況下的系統潮流，如表1所示。以其作為初始運行條件，由隱式改進歐拉法暫態穩定分析程序可計算各發電機功角在故障過程中的變化情況。

表1 潮流計算結果Tab.1 Results of load flow calculation

在暫態穩定計算程序中，各發電機用五階模型模擬，各負荷用恒定阻抗模擬，電力網絡用導納矩陣描述，微分方程用隱式改進歐拉法求解，網絡方程用直接法求解。程序計算了從短路故障開始后，系統在2 s內的暫態過程，得出計及發電機凸極效應時各發電機的功角和相對搖擺角，其中數值積分采用0.001 s的步長。

系統故障過程中各發電機的功角數據如表2所示，計算時間步長取 0.01 s。表 2 中，1-1、2-2、3-3分別表示發電機①、②、③的功角。限于篇幅，表2中只截取了以0.2 s為間隔的數據。

表2 功角數據Tab.2 Data of power angle

系統故障過程中各發電機功角的相對角度如表3所示，即發電機相對搖擺角，計算時間步長取0.01 s，其中，下標2-1表示發電機②相對于發電機①的角度，3-1表示發電機③相對于發電機①的角度。限于篇幅，表3中只截取了以0.4 s為間隔的數據，“*”表示極值時刻。

表3 相對搖擺角數據Tab.3 Data of relative swing angle

由表2和表3可知，雖然發電機的功角不斷增大，但其相對搖擺角的變化并不大，第一擺時的最大搖擺角度 為 δ2-1=126.47259°(t=0.59s)、δ3-1=94.32671°(t=0.64 s)，第二擺時的最大搖擺角度為δ2-1=126.74058°(t=1.88 s)、δ3-1=95.87463°(t=1.92 s)，比第一擺時的角度稍大，2 s內各發電機間的最大相對搖擺角均小于180°。

為了便于分析，利用Matlab7.1對2 s內電力系統發生三相短路故障的功角暫態過程進行仿真，發電機功角和相對搖擺角仿真結果如圖2和圖3所示。圖2表示系統故障過程中各發電機的功角變化曲線，反映出3臺發電機的功角逐漸增大，說明有功出力逐漸增大。圖3表示系統故障過程中各發電機的功角的相對變化，反映出在計及發電機凸極效應時，各發電機間的最大相對搖擺角均小于180°。通過與圖2的比較可以發現，雖然3臺發電機的功角逐漸增大，但其相對搖擺角的變化并沒有超過180°，說明系統在發生該故障時仍能保持穩定，即系統是暫態穩定的，與表2和表3中數據反映的情況一致。

算例結果證明，該方法能夠較好地應用于電力系統暫態穩定的判定過程，同時用隱式解法對微分方程的處理，也更有利于在暫態計算過程中對各變量進行修正，使發電機功角的計算結果更加接近真實值。

5 結束語

本文采用隱式改進歐拉法對電力系統暫態計算中各元件的數學模型進行處理，通過聯立求解的方法進行暫態穩定計算，最后對電力系統多機運行方式下發電機功角的暫態過程進行了分析。

與傳統的改進歐拉法相比，該方法具有積分精度高和無交接誤差的優點，可以滿足長期仿真計算對數值精度的要求。同時，本文介紹的方法保持了電力系統內在的模塊性，在程序設計上具有良好的可擴展性與靈活性。由于選擇了較為理想的系統，忽略了勵磁系統和調相機對電力系統暫態的影響，采用線性化模型得到的功角搖擺曲線還存在一定的誤差，因此，下一步可以采用保留高階項的模型，以求進一步改善計算精度。

［1］倪以信，陳壽孫，張寶霖.動態電力系統的理論和分析［M］.北京:清華大學出版社，2002.

［2］孫暉，趙菁.建立專業電力自動化監控系統的探討［J］.自動化儀表，2003，24(6):15-18.

［3］唐雯，羅安，唐杰，等.電力監控系統圖形繪制及數據處理一體化平臺［J］.自動化儀表，2007，28(8):56-58.

［4］張海生，范征宇.S7-300 PLC和VC++.NET在電力監控系統中的應用［J］.自動化儀表，2006，27(8):50-52.

［5］程若發，馮士芬，江曉舟.基于ARM的同步發電機故障錄波系統的設計［J］.自動化儀表，2007，28(9):11-14.

［6］西安交通大學，清華大學.電力系統計算［M］.北京:水利電力出版社，1978.

［7］Al-Taee M A，Al-Azzawi F J，Al-Taee A A，et al.Real-time assessment of power system transient stability using rate of change of kinetic energy method［C］//IEEE Proceedings of Generation Transmission and Distribution，2001，148(6):505-510.

［8］吳穎，趙巖，蔣傳文，等.風電場穿透功率極限計算方法及發展［J］.自動化儀表，2008，29(11):7-11.

［9］趙艷雷，童建忠.改進歐拉法在電力系統暫態分析中的應用和軟件設計［J］.微計算機信息，2004，20(5):98-99.

［10］Kasztenny B，Kezunovic M.A method for linking different modeling techniques for accurate and efficient simulation［J］.IEEE Transactions on Power Systems，2000，15(1):65-72.

［11］Suyono H，Nor K M，Yusof S.Component based development for transient stability power system simulation software［C］//Proceedings of Power Engineering Conference，PECon，2003:16-21.

［12］陳珩.電力系統穩態分析［M］.2版.北京:中國電力出版社，1995.

［13］武東亞，王克文，馬潔.基于概率潮流的多運行方式下發電機功角的暫態過程分析［J］.電氣應用，2009，28(6):64-67.

［14］王錫凡，方萬良，杜正春.現代電力系統分析［M］.北京:科學出版社，2003.