電力系統電磁暫態仿真算法研究綜述

楊明，張永明，張子騫，顧禹軒

(1. 上海電機學院 電氣學院， 上海 201306； 2. 格拉茨技術大學 電力系統研究所，奧地利格拉茨)

0 引 言

目前大規模新能源的并網以及新型電力電子器件的應用，電力系統數字仿真復雜性不斷提高[1-2]，電磁暫態特性也變得日趨復雜。在實時仿真系統中，使用數字計算的仿真必須無限接近物理實際系統，要求仿真時間與被仿真系統時間相對應并能在一個仿真步長內快速計算所有差分方程。實時數字仿真器(Real Time Digital Simulator，RTDS)是能將電力系統電磁暫態仿真實時進行分析的一種并行系統[3-4]。電磁暫態仿真分析根據各類模型的不同特性對電力系統中相關元器件進行精準建模，準確地獲得電力系統的動態特征，在分析電力系統中不同暫態過程時，這就要求電磁暫態仿真具備高性能仿真能力，目前基于Dommel算法的電磁暫態程序EMTP(Electromagnetic Transients Program，EMTP)通常應用在電力系統電磁暫態仿真中，我國中國電力科學研究院研發的EMTPE和加拿大的PSCAD/EMTDC都是基于EMTP程序而開發。為克服未來電網電磁暫態實時仿真大規模的挑戰，基于HYPERSIM的電磁暫態實時仿真技術的研究也有了新突破[5]。電力系統電磁暫態仿真程序軟件一般采用定步長的數值積分算法用于描述電力系統中的基本元器件，同時對各元器件進行離散化處理，考慮到仿真計算過程中的穩定性與精度問題，梯形法成為各種數值分析方法的首選方法，但梯形法最大的問題是在數值求解過程中產生數值振蕩現象。另外在電磁暫態仿真分析過程中也會因其它各種原因出現同樣現象，使仿真過程不準確。數值積分代換技術是解決電磁暫態問題的普遍方法，主要特征在于其簡單性、全面應用性以及計算高效性[6]。對于電力系統的剛性特性，算法會在求解穩定性上受到影響，非線性特性也會影響求解過程的計算效率。由于數值積分算法特性和性能的差異性，因此選擇合適的算法是電磁暫態仿真分析的基礎和核心。

近年來，一些新算法和新技術的研究使電磁暫態仿真在其穩定性、計算精度和仿真效率上得到提高，應用領域也更為廣泛。文章討論了目前應用在電磁暫態仿真中的經典積分算法和研究者提出的新型積分算法，并通過討論算法的性能指標對比分析了所應用方法的優缺點。

1 仿真理論基礎

1.1 理論基礎

電力系統中各類數學模型和與之對應的數值算法的求解構成了電力系統電磁暫態仿真過程[7]，例如圖1所示的電感支路。

圖1 電感支路及諾頓等效電路

由圖1可得其微分方程如式(1)所示：

(1)

對式(1)方程進行差分化，由此可得到差分方程[8-10]。例如對式(1)采用梯形積分法后可得：

(2)

式(2)差分方程可表示成如圖1所示諾頓等效電路(伴隨電路)，作用于當前時步電壓項部分等效為電導，前一時步電氣量并作用于當前時步電流項部分等效為電流源歷史項，通過節點方程聯立式(2)可得式(3)所示的電磁暫態仿真基本方程：

Gu=i

(3)

節點分析法在大規模電力系統中相較于狀態方程法在計算速度上更占優勢，例如傳統電磁暫態仿真程序即采用了節點分析法。

1.2 問題與措施

隨著電網規模的不斷擴大，電力系統大量應用電力電子裝置的趨勢使其復雜性日趨增加，傳統電磁暫態仿真性能已經無法滿足需求[11]，對電磁暫態仿真分析的精準建模能力要求也不斷提高。在仿真過程中，如開關動作過程中存在電感電容元件的換路過程，將會引起非狀態變量的突變，若繼續采用梯形積分法將會引起非原型的數值振蕩[12-14]；同樣，因開關動作不在整步長時刻點，定步長算法隱開關的延遲會引入非特征諧波，導致仿真失真[15]；此外控制系統與主系統之間產生的外部時延問題也將導致數值波動現象。由于控制系統中非線性環節產生的內部時延，也會導致在控制系統求解時出現數值不穩定問題[16]。

針對數值振蕩產生的原因分析，研究者提出了許多解決方法：(1)求出突變后的非狀態變量。但對于實 際復雜電力系統不切實際；(2)避免在網絡結構突變時引入非狀態變量。例如在采用后退歐拉法時，由于后退歐拉法在求解方程中不會出現t-Δt時刻的非狀態變量，因此在突變過程中，不會產生不合理的等值注入電流源，從而避免了數值振蕩問題，但這種解決方式精度不高[17]；(3)減小時間步長。由于數值振蕩的幅值隨1/Δt變化，因此較小的時間步長反而會造成數值振蕩問題，并不普遍適用；(4)添加阻尼電阻。該方法中阻尼電阻Rd與步長Δt、阻尼因子α選取有關，但選取值不能保證最優，影響計算精度，并且當電阻值選取太小時，反而會加劇數值振蕩[18]；(5)文獻[19]通過適當的控制模型來最小化或補償解決外部時延。另外文獻[20]提出了一種基于濾波的抑制方法解決外部時延。為消除內部時延，文獻[21]采用牛頓法來獲得控制系統的聯立求解。

圍繞電力系統電磁暫態仿真需求，數值積分算法是開展電力系統非線性剛性系統的基礎與核心，由于數值積分算法的不同特性，應考慮到算法的穩定性、仿真精度以及仿真效率。下一章則將著重分析所提及到的算法性能指標。

2 數值積分算法的性能指標

2.1 計算穩定性

(4)

對于線性微分動力方程，具備L穩定性的數值積分算法也能夠消除數值振蕩，基于對剛性系統的適應性以及仿真結果的準確性要求，選取電力系統仿真積分算法時，考慮計算穩定性是準確模擬電力系統電磁暫態仿真的前提。

2.2 計算精度

不同數值積分算法，相應的精度和誤差也不盡相同。在采用數值方法求解微分方程時會存在截斷誤差，例如當yn=y(xn)時，用數值積分算法計算yn+1的誤差：

Rn=y(xn+1)-y(xn)

(5)

式(5)中Rn即為局部截斷誤差。通過泰勒公式展開可得數值方法的局部截斷誤差為O(hP+1)，其中P為階數。當步長h越小時，P值就越高，而局部截斷誤差與P值成反比，Rn將不斷減小，因此計算精度就越高。顯然若在同等步長h的數值積分算法中，階數越高，計算精度則越高。

2.3 計算效率

數值積分算法的計算效率是滿足電磁暫態仿真實時性要求的關鍵。顯示積分法仿真效率較高，但穩定性差，很少應用于電磁暫態仿真中。文章主要選擇隱式積分法進行計算效率的對比。可通過采用高斯消元法求解離散后的方程[23]，以求解后的浮點運算數作為計算效率的指標。

文章對數值積分算法的討論中，對圖2所示RL電路進行求解計算，以式(6)所示微分方程一般形式為基礎，對運用在電力系統電磁暫態仿真的歐拉法、梯形法、均值法、CDA法、THTA法以及新型積分算法(塊廣義向后差分法、精細積分法、根匹配技術)進行討論。其中x(t)為狀態變量；f為時間t與x的函數。

(6)

圖2 進行計算的電路

3 經典積分算法研究

3.1 歐拉法

3.1.1 前向歐拉法

前向歐拉法的積分格式如式(7)所示:

yn+1=yn+hf(tn,yn)

(7)

前向歐拉法采用固定步長，計算簡單，但由于是顯示積分，等效電阻為負，因此穩定性差。在選取為固定步長后，當仿真步數越多，則計算誤差就越大，當步長較大時，前向歐拉法的精度并不高，因此在實際數值求解中不經常使用。前向歐拉法的穩定域如圖3所示。

圖3 前向歐拉法的穩定域

3.1.2 后向歐拉法

后向歐拉法采用固定步長，計算簡單，計算效率高，積分格式為：

yn+1=yn+hf(tn,yn+1)

(8)

圖4 后向歐拉法的穩定域

該方法在求解方程中不會出現t-Δt時刻的非狀態變量，因此在突變過程中，不會產生不合理的等值注入電流源，從而不會引起數值振蕩問題。后向歐拉法的等效電阻為Re=2L/Δt，具有強阻尼特性，可抑制諧振電路中的高頻振蕩，但對低頻振蕩的抑制作用不佳[24]，并且也會將原來仿真中本應該顯現出來的振蕩消除而不能夠準確的模擬整個仿真過程。因相位畸變的存在，在仿真的過程中也不能完全使用后向歐拉法，另外因阻抗的存在也會使得離散系統產生附加損耗。

對圖2中RL電路采用后向歐拉法時仿真波形如圖5所示。

圖5 后向歐拉法

后向歐拉法通常結合其他方法一同仿真計算發揮其穩定性優勢，基于MATLAB平臺的PSAT(Power System Analysis Toolbox)仿真工具即采用了此方法[25]。

3.2 梯形積分法

3.2.1 隱式梯形積分法

隱式梯形積分法的積分格式為：

(9)

局部截斷誤差采用泰勒級數展開方式分析如式(10)所示:

(10)

與泰勒級數相對比可知，隱式梯形積分法的局部截斷誤差為O(h3)，則P值為2，因此具有2階精度。

對圖2中RL電路采用隱式梯形積分法時仿真波形如圖6所示。

圖6 隱式梯形積分法

圖7 隱式梯形積分法的穩定域

對任意仿真步長均可以保證截斷誤差衰減，隱式梯形積分法是A穩定方法中具有最小誤差常數的一種數值方法[26]，在仿真中可利用其穩定性優勢，聯立同步發電機和網絡的差分與代數方程計算后，可有效消除接口誤差[24]。但由于梯形法不是L穩定性的，因此在仿真中若出現因網絡拓撲結構變化而導致非狀態變量突變時，將產生不衰減的數值振蕩現象[27-28]。

梯形積分法的研究已經比較成熟，但許多研究者將其與其他方法相組合形成新方法，為目前實時仿真分析提供了更多的擴展空間。針對梯形積分法產生數值振蕩問題，文獻[29-30]提出了基于梯形積分法的改進算法，即阻尼梯形法，此方法將梯形法和后向歐拉法加權混合，相當于在電感上并聯小電導，在電容上串聯小電阻。阻尼梯形法引入了阻尼因子α，精度可介于梯形法與后向歐拉法之間，該方法的精度與阻尼因子α的大小選取有關[30]，然而阻尼因子的選取并不是能夠完全確定的，都是通過估計得來，因此文獻[31]在分析阻尼梯形法誤差的基礎上，提出了阻尼梯形法的補償修正法，提高了仿真精度。文獻[32]提出的龍-庫-梯法將龍格-庫塔法與梯形法相結合，相互彌補不足，使仿真更加精確可靠，對數值振蕩也有良好的衰減作用。

隱式梯形積分法在求解過程中，其浮點數運算量與后向歐拉法相當，計算效率高。由于隱式梯形積分法對于電力系統這樣的剛性系統具備良好的適應性，因此在實際仿真中被廣泛應用。傳統電磁暫態EMTP仿真程序因其精度高與穩定性好等優勢而廣泛應用。另外電力系統商業計算程序BPA、PSASP也釆用了隱式梯形積分法。

3.2.2 均值法

為消除非原型數值振蕩，EMTP曾采用算后處理法，即在仿真計算過程中依舊采用梯形積分法，不對數值振蕩做特別處理，而在輸出時將振蕩的平均值作為仿真結果。如式(11)所示，由梯形積分法在t-Δt時刻步長與t時刻步長得到的值取平均后輸出[33]。與隱式梯形積分法對比表明具有便捷性和有效性，計算簡單，對數值振蕩有一定抑制效果，但穩定性較差。

(11)

對圖2中RL電路采用均值法時仿真波形如圖8所示。

圖8 均值法

3.3 CDA法

臨界阻尼調整法(Critical Damping Adjustment, CDA)是由文獻[34]提出。采用變步長，計算復雜，實際上是一種積分算法相結合的方式，仍以梯形法作為方程求解計算的基礎，主要計算步驟是僅在網絡突變發生時，采用后向歐拉法在兩個半步長中進行計算，此方法將數值振蕩從一開始就得以消除，因此對數值振蕩有一定抑制作用[35-36]。CDA法利用其臨界阻尼特性，避免了電導矩陣的變化，為電路中的所有元件帶來相同的等效電導，節省了計算工作量[37]，但算法的相互轉換卻導致系統編程更為繁瑣，切換時刻的確定造成了計算負擔，影響仿真速度，計算效率低。受組合方式的啟發，文獻[38]將一種新型4步線性多步法與梯形積分法結合切換計算，精度更高且能有效避免數值振蕩。

采用CDA法消除數值振蕩的前提條件是需檢測出突變現象的發生時刻，而在實際情況下突變類型多種多樣，如控制系統中電壓源和電流源因限幅環節的影響，難以檢測到突變現象情況[39]，此時CDA法無法準確抑制數值振蕩。另外因暫態仿真中延遲出現的突變時刻，會造成CDA法無法準確檢測與抑制振蕩，因此文獻[40]提出基于2級3階Radau II A法的并行求解方法，無需檢測突變時刻即可消除振蕩，且并行求解方式也提高了計算效率，該方法也適用于線性與非線性常微分初值問題。

對圖2中RL電路采用CDA法時仿真波形如圖9所示。

圖9 臨界阻尼調整法

另外CDA法中包含后向歐拉法，在長期模擬結果中會引入相移問題，計算時會導致精度降低并產生誤差[41]。在含有大量電力電子開關器件開關操作時，CDA法只有在正常計算時可保持2階精度，在切換算法消除數值振蕩過程中，算法精度會降至1階。因此文獻[42]提出了3S-DIRK變階變步長方法，可在正常計算時具備3階精度，其余計算可切換至具有L穩定性的算法抑制數值振蕩并可保證精度不低于2階。

實際應用中，軟件Microtran實現了CDA技術[43]，EMTP的DCG-EPRI版本(即EMTP版本3.0)中也采用了臨界阻尼調整法，雖然消除了數值振蕩，但仿真波形中會出現一些處理痕跡，結果并不理想[44]。

3.4 THTA法

THTA法(Trapezoidal History Term Averaging)采用變步長，具備A穩定性，計算較復雜。此方法仍以梯形積分法為基礎，在Δt/2的兩個步長中確定歷史源值的平均值。如式(12)所示，先采用梯形積分法在時間t中計算歷史電流源值，再計算平均值ITHTA(t)，并將此值作為新的歷史電流源值計算VLTHTA(t)，第二個半步長重復以上步驟計算，在數值穩定后，繼續采用隱式梯形積分法求解。

(12)

采用THTA法對各步長中的參數分析列舉，如表1所示。

表1 THTA法對各參數分析

對圖2中RL電路采用THTA法時仿真波形如圖10所示。

該方法在仿真過程中保持了梯形積分法精度高的優勢，消除了誤差，僅在Δt/2步長中出現一次抖動，之后保持穩定，對數值振蕩有一定抑制作用。THTA法與CDA法不同的是，它只需上一步驟的歷史源值產生一個新的值就能免受數值振蕩，計算效率更高。

圖10 THTA法

4 新型積分算法研究

目前大規模新能源發電接入電力系統后，呈現出電力電子化趨勢[45]，大量具有寬頻帶響應特性的電力電子裝置不斷接入至電力系統中，極大改變了電力系統的動態特性[46]，使電磁暫態仿真的分析復雜性不斷增加。電力系統電磁暫態仿真在經典積分算法的基礎上不斷發展改進，目前越來越來多的新型積分算法和技術不斷提出并獲得應用。

4.1 塊廣義向后差分法

塊廣義向后差分方法(Generalized Backward Differentiation Formula,GBDF)是塊邊界值方法中的一類，由文獻[47]提出。具備A穩定性和L穩定性[48]。精度可在1階(也稱為后向歐拉法)～5階之間變化，采用2或3階的GBDF方法時可求解初值問題，因具有較強的阻尼特性，可以得到在時間上很平滑的解，因此在電磁暫態仿真中可消除數值振蕩現象。當采用s級s階的GBDF方法求解式(6)一階微分方程的初值問題，計算格式可描述為式(13)～式(15)：

(13)

式中N為時間離散網格點數；h為時間積分步長；αi,i∈(0,s)的取值可參考文獻[49]。

(14)

(15)

同時采用Brugnano等提出的附加方法選擇策略與邊界值方法聯合求解待求點數N相匹配的N個方程，將不會對數值穩定性與計算精度產生影響。采用塊廣義向后差分法可允許較大的時間積分步長求解離散后的一階模型得到離散點的時域解，相較于經典積分算法中的臨界阻尼調整法，無需檢測網絡結構的突變現象，減輕了計算負擔，且能有效避免數值振蕩，適用于輸電線路電磁暫態仿真中[50]。

電力系統規模與復雜性的增加，對電力系統電磁暫態實時仿真要求也更加嚴格，為達到實時仿真分析計算要求，文獻[51]將塊廣義向后差分法與擴展的隱式梯形積分法結合應用在非線性常微分初值問題的數值計算中，解決了高維非線性初值問題面臨的維數災問題。文獻[52]利用塊廣義向后差分法，采用矩陣分解，不會因整體雅可比矩陣或多個分塊子矩陣被三角分解而降低效率，數值計算效率比經典微分求積法更高。

4.2 精細積分法

電力系統電磁暫態仿真過程需要解決線性或非線性系統，但最主要的還是先將線性時不變的系統精細的分析好[51]。文獻[53]以Duhamel積分為基礎，提出了矩陣指數計算的精細積分方法，在此基礎上又提出解決瞬態初始問題時保持高精度、高穩定性的求解方法。對于線性常微分方程，其解基于Duhamel積分可表示為式(16)所示：

(16)

通過等時間步長η將時間劃分，即可得到式(16)的遞推形式：

(17)

其中，eHη=I+Ht+(Hη)2/2+(Hη)3/3!+…；H為常數矩陣，非齊次項f為時間函數。選取不同的Duhamel數值積分方法可衍生出精細積分法的其他形式，雖然本質上為顯示積分，但例如高精度直接積分法利用拉格朗日高階插值多項式，在處理非線性問題上，精度和穩定性依然有優勢[54]。

精細積分法將微分變量的導數部分細分為線性部分和非線性部分，是用于計算矩陣指數的一種方法，對于非齊次方程，通過增維方法，將其擴容為齊次方程，最后采用精細積分法求解時，可避免矩陣的求逆計算，提高了計算速度，但精度不高。當在求解非線性剛柔耦合問題時，精細積分法處理剛性問題也有獨特優勢[55]。文獻[56]通過將精細積分法進行增維的方式應用到電磁暫態計算中，一定程度上解決了數值振蕩問題，但對于大規模電力系統，龐大的計算量增加了誤差，計算精度降低。文獻[57]提出一種基于泰勒級數的精細計算方法，整合的過程劃分為更小的部分，運用泰勒級數展開計算指數矩陣，但存在嚴重的舍入誤差，因此文獻[58]提出改進的泰勒級數精細積分法，并將此運用到電磁暫態數值計算中，相較于經典積分算法，既保留了算法的優點，又能夠避免數值振蕩問題，計算精度也得到提高。為處理采用隱式梯形積分法產生的問題，文獻[59]在精細積分算法的基礎上，提出了隱式精細積分算法，比經典積分算法精度更高，更易于實現。

精細積分法有較好的穩定性，但在時間上劃分的許多區間會使計算效率降低，在實際應用中影響很大，因此文獻[60]采用組合稀疏矩陣技術，提出了精細積分級數解的并行算法。為解決電力系統輸電線中過電壓數值模擬的計算效率問題，在精細積分法的基礎上，文獻[61]提出高精度直接積分法，相較于CDA法小步長仿真，大步長計算效率更高，具有良好的數值穩定性。為避免狀態矩陣求逆影響仿真速度，文獻[62]提出基于精細積分法的多步法，具備良好的穩定性和高精度優勢，也極大地提高了仿真效率。另外文獻[63]將非線性部分作線性化假設，使增維的精細積分法拓展到了非線性系統研究領域，應用領域進一步拓寬。

4.3 根匹配技術

基于截斷泰勒級數的積分方法在仿真階躍響應時，會產生數值振蕩問題，相較于經典積分算法，差分方程的指數形式(Root-matching Techniques，根匹配技術)是一種替代方法。

采用指數形式能夠在循環計算之前，預先計算好指數項并將其存儲，計算效率并比隱式梯形積分法在數值計算上更高效[64]。該方法通過確定電路結構計算得到其S域的傳遞函數H(s)，由傳遞函數可得到其差分方程的指數形式，與Dommel方法得到的差分方程類似，可將差分方程的指數形式看作一個諾頓等效形式。文章以RL電路為例，其傳遞函數如式(18)所示:

(18)

差分方程指數形式的諾頓表達式如式(19)所示：

(19)

傳統電磁暫態仿真直接在時域內進行差分求解，而根匹配法從頻域離散系統出發，通過反變換進行差分化，具有良好的穩定性，而且保持了原連續模型的動態特性，精度高于傳統積分算法。然而根匹配法并非適合任何元件，針對RL并聯與RC串聯回路并不能有效避免非原型振蕩以及根匹配法無法對獨立電感、電容元件進行差分處理的問題，文獻[65]提出了基于離散相似的電磁暫態仿真方法，首先將組合形式拆分為獨立元件，再添加正負虛擬電阻轉換為RL串聯與RC并聯的組合元件，經過離散相似法差分后，歷史項中不含非狀態變量，因此有效避免了數值振蕩。

根匹配技術總結有以下優點[66]：(1)能夠消除截斷誤差，并因此消除了數值振蕩。與采用的步長無關，能夠同時應用于電網絡方程和控制模塊；(2)能夠等效成諾頓等效電路，可與數值積分代換方法完全兼容，并且矩陣求解技術保持不變；(3)能夠提供高效和準確的時域仿真解，目前可應用于PSCAD和EMTDC仿真軟件中。

5 結束語

文章根據電磁暫態仿真分析的高性能需求討論了經典積分算法和新型積分算法，主要結論如下：

(1)后向歐拉法具備A穩定性和L穩定性，可避免數值振蕩，但精度低，阻抗的存在也會產額外損耗。隱式梯形積分法具備A穩定性和2階精度，計算效率高，在電力系統仿真中應用最為廣泛，但不是L穩定性，因此不具備消除數值振蕩的能力。臨界阻尼調整法能夠從源頭上消除數值振蕩，但需檢測突變發生時刻，算法切換增加了計算負擔，在計算精度和計算效率上不占優勢。以梯形法為基礎的THTA法僅一個時間步長即可消除數值振蕩并保持穩定，比CDA法計算更簡便有效;

(2)經典積分算法的研究主要集中在不同特性算法的組合運用以發揮各自優勢，為電力系統電磁暫態仿真提供了更多發展空間。如受臨界阻尼調整法的啟發，將隱式梯形積分法與一類新的低階、L穩定的線性多步法相結合，比較適用于含有大量電力電子裝置的電力系統電磁暫態仿真中;

(3)新型積分算法中精細積分法具備高精度和良好的穩定性。為滿足實時分析計算需求，在求解大規模問題中可采用并行計算技術提高精細積分法的計算效率。精細積分法本質上適用于線性定常部分的求解，利用其處理線性剛性問題的優勢可構造出更高性能數值方法。如何利用精細積分法優點求解非線性問題，在高性能電力系統電磁暫態仿真中尚有研究價值。根匹配技術不僅可以形成與Dommel算法結構相似的諾頓等效電路，可較為容易的應用在現有電磁暫態程序中并有效避免了非原型數值振蕩，也在精度和穩定性方面大大提高，但根匹配法為恒穩定算法，并非適用于任何元件，針對基于根匹配法的非線性以及存在耦合關系的元件處理還需進一步研究。