計及多重開關的電磁暫態仿真插值算法

劉 濤，閆貽鵬，金 娜，王旭東，林濟鏗

（1．天津市電力公司電力科學研究院，天津 300384；2．上海市電力公司青浦供電公司，上海 201700；3．上海市電力公司金山供電公司，上海 201500；4.同濟大學電子與信息工程學院電力系，上海 201804）

計及多重開關的電磁暫態仿真插值算法

劉 濤1，閆貽鵬2，金 娜3，王旭東1，林濟鏗4

（1．天津市電力公司電力科學研究院，天津 300384；2．上海市電力公司青浦供電公司，上海 201700；3．上海市電力公司金山供電公司，上海 201500；4.同濟大學電子與信息工程學院電力系，上海 201804）

為解決快速電磁暫態仿真程序中因多重開關動作易引起迭代震蕩及仿真精度不高問題，提出一種能避免迭代震蕩及具有較高精度的快速電磁暫態仿真插值新算法。該算法首先通過內插技術確定開關動作的準確時刻；然后從最近的開關動作時刻開始對本步長內該開關動作點之后的區間進行后向歐拉法試探積分；進而利用外插回推技術得到開關動作后的變量值作為下一步積分的初始值；最后進行一步梯形法積分并插值回到積分整步點，完成一個完整時步積分計算。算例證明本算法能避免迭代震蕩并提高仿真精度，擁有應用前景。

插值技術；電磁暫態仿真；數值振蕩；多重開關

目前，電力電子元件及設備在電力系統中的應用非常廣泛。由于電力電子元件的快速多次反復動作特性造成了系統的網絡拓撲結構也是隨時間不斷變化，且很可能在一個很短的時間內出現多個開關動作。如何準確地反映及捕捉各個開關動作給網絡及系統帶來的影響及現象，是所有電磁暫態仿真研究及程序開發所面臨的挑戰。

基于定步長積分技術而開發的EMTP程序是目前應用最為廣泛的電磁暫態仿真程序[1，2]。在開關動作時，若實際開關動作時間出現在非整步點上，而積分過程卻只能在整點上反映開關動作情況，往往造成開關動作點延遲，并且會出現非特征諧波[3-5]和數值振蕩[6-10]。對于該問題，目前常采用的是插值算法[11-15]，即利用梯形積分中的分段線性化，在開關電流過零點時利用插值技術尋找到精確的開關動作時間。而對于在一個積分時步中包含有多個開關動作，即出現所謂的多重開關[16]時，當前的插值算法大都只是尋找一次開關動作或將多個動作近似到一個時間點上，從而不能準確反映實際情況，且也容易引起積分震蕩。如何更準確地模擬多重開關動作問題一直是備受關注的課題。

本文首先對包含有多重開關的電磁暫態積分過程中的數值振蕩和非特征諧波的產生機理進行理論分析；在此基礎上，針對多重開關動作問題，提出了計及多重開關動作的插值新算法用于電磁暫態仿真積分計算，以消除仿真計算時出現的數值振蕩和非特征諧波。

1 非特征諧波和數值振蕩

EMTP等電磁暫態仿真程序均采用定步長的梯形積分法，從而只能計算整步點變量值。然而電力電子開關的動作點不一定出現在整步點上，以圖1為例，當開關斷開后，由于電感電流不能突變，二級管導通，在電感電流衰減至零時二極管關斷，假設在td時刻二極管關斷，如圖2所示，則在第一個整仿真步長時刻（即圖2中t-Δt）程序計算得到的二極管電流為正，二極管繼續導通，在第二個整仿真步長時刻（即圖2中時刻t）程序計算得到二極管電流為負，這個信息直到下一個整仿真步長時刻（即圖2中t+Δt）才被處理，將二極管關斷。而二極管電流實際上是在第一個仿真步長和第二個仿真步長之間的時刻（即圖2中td）過零點。

圖1 簡單電力電子開關電路Fig.1 Simple power electronic switch circuit

圖2 固定步長仿真時開關動作示意Fig.2 Switch action schemes at fixed-step simulation

程序中開關動作時間上的延遲將造成電壓電流波形出現不真實的“尖峰”，即非特征諧波。因此，只有找到真正的開關動作點，讓仿真中電子開關的動作回歸到正確時刻才能消除非特征諧波。

對于圖1中電感支路，其梯形法積分公式為

如果開關在td時刻動作，電壓u將發生突變，即UL（td-）≠UL（td+），在仿真中應以td+時刻的變量求解td+Δt時刻的值，但實際計算時常用td-時刻非狀態變量值計算等值電流源，而由于UL（td-）≠UL（td+），這就造成了數值震蕩。

目前常采用向后歐拉法消除數值震蕩，向后歐拉法的公式為

從式（2）可看出，該方法避免了非狀態突變量UL（td-）的影響，因此能夠徹底地消除數值震蕩。

2 新的插值算法

以上分析可知，要仿真電力電子元件的動作過程，首先要尋找到正確的開關動作點，然后采用向后歐拉法，從而既能消除震蕩產生，又能正確反應開關動作。基于此思想，同時考慮到快速仿真中往往需采用大步長而導致一個積分時步中可能包括多重開關動作，本文提出一種新的插值算法，既能對多重開關的開斷準確反映，又能避免震蕩。

算法的核心思想是，在從t積分到t+Δt時刻時檢測到開關動作td，以t和t+Δt時刻的變量插值到td，得到td（0-）變量，接著檢查在td到t+Δt直接是否有開關動作，若有則以td和t+Δt為基準插值到新的開關動作點td′，再從新的開關動作點用一步1%的向后歐拉法以消除振蕩，再通過外插回推得到td（0+）時刻的變量，最后回到整步點，開始下一時步的仿真積分。具體過程如下：

（1）由t-Δt時刻使用一步梯形法積分法到t；

（2）由t時刻使用一步梯形法積分法到t+Δt，此時檢測到td時刻開關動作；

（3）由t和t+Δt時刻內插到td時刻值；

（4）改變開關狀態，重新形成導納矩陣，從td時刻以步長為t+Δt-td的向后歐拉法積分到t+Δt時刻，檢測td到t+Δt之間是否還有開關動作，若有則由td和t+Δt時刻內插到新的開關動作點，重復步驟（4），若沒有，則轉入步驟（5）；

（5）在td時刻，用向后歐拉法以1%步長積分一步到td+Δt/100；

（6）由td+Δt/100和td時刻值外插到td-Δt/100時刻值；

（7）在td+Δt/100時刻以1%步長的向后歐拉法積分一步到td；

（8）在td時刻以1%步長的向后歐拉法積分一步到td+Δt/100；

（9）在td+Δt/100時刻用一整步梯形法積分至td+Δt/100+Δt；

（10）由td+Δt/100和td+Δt/100+Δt時刻值內插到t+Δt整步時刻，恢復正常梯形法積分直到檢測到下一個開關動作點。插值變化過程如圖3所示。

圖3 新插值算法的插值變化過程Fig.3 Interpolating process for the new interpolationextrapolation algorithm

上述插值算法的解釋和說明：

（1）第（3）步通過內插技術得到開關動作前的變量值；第（4）步進行因開關動作而進行相應的網絡拓撲改變；第（4）步繼續尋找本步長內的其他開關動作點，并插值到新的開關動作時刻同時改變拓撲結構，直到找到最后一個開關動作點為止；第（5）步利用1%步長的向后歐拉法消除振蕩；第（6）、第（7）步利用外插回推過程得到開關動作后的變量值，即td（0+）時刻值；第（8）步利用td（0+）時刻的值向前積分，避免了td（0-）值的影響；第（9）、第（10）步利用內插技術插值到整步點上，完成下一步梯形法積分的準備。

（2）在上述算法第（4）步中，采用了試探檢測來尋找同一步長內的多次開關動作，且每次尋找到一個新動作點都進行內插計算和網絡拓撲改變，雖要多花去一部分計算時間，但通過該操作后能準確尋找到每一次開關動作時刻，從而為下一步積分提供正確的初始狀態，可明顯提高仿真精度。

（3）因電磁暫態仿真精度取決于數值積分方法和步長，梯形積分法的精度為二階，向后歐拉法精度為一階，因此在第（5）步中采用1%步長向后歐拉法，相比文獻[12，13]中采用的半步長向后歐拉法策略既能消除振蕩且精度更高。同時，本方法回到了整步點輸出，能與導納矩陣和暫態模型相匹配。由于本方法將計算得到的td（0+）時刻值用于下一步積分，相比文獻[14]獲得了更加精確的初始值，故總體仿真結果也具有更高精度。

（4）本文方法相對于非插值算法，相同步長下計算時間較長。但相對于插值類算法，雖計算量略有增加，但增加的不大。本文方法的突出優點是具有較高的仿真精度和穩定性，雖有一定的計算量增加，但相對于現在計算機的計算能力，這樣的增加幾乎可以忽略不計。

3 算例

3.1 算例1

如圖1所示，在開關斷開后，二極管導通為電感續流，當電感電流為零時二極管關斷。分別利用文獻[13，14]和本文插值法對該過程進行電磁暫態仿真，給出二極管關斷時電壓V的波形，并與10 μs步長仿真時的電壓波形對比，如圖4和圖5所示，其中圖4采用50 μs步長，圖5采用100 μs步長。

圖4 50 μs步長時電壓波形比較Fig.4 Comparison of the voltage waveforms at 50 μs time-step

通過圖4（a）和圖5（a）比較可知，50μs和100μs步長時，不利用插值算法的仿真結果誤差很大，而利用插值算法后的仿真波形接近10 μs仿真結果，說明插值技術對程序的精度有很大提高。

圖5 100 μs步長時電壓波形比較Fig.5 Comparison of the voltage waveforms at 100 μs time-step

由圖4（b）和圖5（b）可看出，本文提出的新插值方法相比文獻[13]，[14]的方法更接近于10 μs的仿真結果。從而也證明本文方法相較其他方法具有更高的仿真精度和更好的數值穩定性。

3.2 算例2

算例2為一完整整流逆變電路，如圖6所示。有4個換流器，24個電力電子開關元件，并包含分布式參數線路，濾波器等元件。在這個電路中，開關動作不僅包括自然換相也包括強制換相，而且用大步長仿真時會出現多重開關問題。本文采用頻率50 Hz，仿真步長分別為50 μs和100 μs，整流側觸發角α=15°，逆變側β=45°，仿真時間為1 s，測量長線路中靠近整流側端的電壓值。

圖7和圖8分別給出了50 μs和100 μs步長時本文插值方法、文獻[13]，[14]中的插值方法以及沒有插值時的電壓曲線。可看出利用插值方法后，直流電壓的振蕩較少、波形平穩，說明3種插值方法對于開關動作引起的振蕩有很好消除。然而，本文方法相對于文獻[13]，[14]方法得到的波形更加平穩，說明在包含多個開關動作的系統仿真中，本文方法更好地消除了數值振蕩和非特征諧波，具有更好的仿真精度和數值穩定性。同時從圖7和圖8的比較可知，步長從50 μs增加到100 μs時，對于本文插值方法得到的波形平穩度基本不受影響，說明本插值方法在大步長時仍適用，即對于大步長的電磁暫態仿真依然保證了一個較高精度，從而適用于大系統快速電磁暫態仿真。

圖6 整流逆變系統圖Fig.6 Rectifier inverter system

圖7 50 μs步長時直流電壓比較Fig.7 Comparison of DC voltage waveforms at 50 μs time-step

圖8 100 μs步長時直流電壓比較Fig.8 Comparison of DC voltage waveforms at 100 μs time-step

4 結語

針對快速仿真中多重開關的問題，本文提出一種適用于快速電磁暫態仿真的插值新算法。該方法的特點是通過內插技術尋找正確的開關動作點，采用試探法尋找一個步長內的多次開關動作，利用1%步長的向后歐拉法有效消除了數值振蕩；并利用外插回推技術得到正確的網絡拓撲結構和積分初始值。算例表明，本文算法相較于其他方法在處理多重開關時，既能更好地避免數值震蕩，又有更高的仿真精度。算例也表明，本文算法在較大積分步長時，仍具有較高的數值穩定性和精度，從而擁有實際工程應用的前景。

[1]Dommel H W．Digital computer solution of electromagnetic transient in single and multiphase networks[J]．IEEE Trans on Power Apparatus and Systems，1969，88（4）：388－399．

[2]De Kelper B，Dessaint L A，Al-Haddad K，et al．A comprehensive approach to fixed-step simulation of switched circuits[J]．IEEE Trans on Power Electronics，2002，17 （2）：216－224．

[3]Watson N，Arrillaga J．Power Systems Electromagnetic Transients Simulation [M]．London：The Institution of Electrical Engineers，2003．

[4]王成山，高毅，王丹，等（Wang Chengshan，Gao Yi，Wang Dan，et al）．考慮直流系統開關特性控制的變步長仿真算法 （Variable-step simulation method considering HVDC controls with switching characteristics）[J]．中國電機工程學報（Proceedings of the CSEE），2009，29（34）：16-21．

[5]劉文焯，湯涌，侯俊賢，等（Liu Wenzhuo，Tang Yong，Hou Junxian，et al）．考慮任意重事件發生的多步變步長電磁暫態仿真算法（Simulation algorithm for multi variable-step electromagnetic transient considering multiple events）[J]．中國電機工程學報（Proceedings of the CSEE），2009，29（34）：9-15．

[6]張益，周群（Zhang Yi，Zhou Qun）．電力系統數字仿真中的數值振蕩及對策（Numerical oscillation and its counter measures in power system simulation）[J]．上海交通大學學報（Journal of Shanghai Jiaotong University），1999，33 （12）：1545-1549．

[7]應迪生，張明（Ying Disheng，Zhang Ming）．動態相量法仿真中的數值振蕩及消除（Elimination of numerical oscillation in power system simulation using dynamic phasor）[J]．電力系統及其自動化學報（Proceedings of the CSU-EPSA），2007，19（4）：41-45，128.

[8]劉益青，陳超英（Liu Yiqing，Chen Chaoying）．用以消除數值振蕩的阻尼梯形法誤差分析與修正（Error analysis and correction of damping trapezoidal integration for eliminating numerical oscillation）[J]．中國電機工程學報（Proceedings of the CSEE），2003，23（7）：57-61．

[9]Lin Jiming，Marti J R．Implementation of the CDA procedure in the EMTP[J]．IEEE Trans on Power Systems，1990，5（2）：394-402．

[10]Marti J R，Lin J．Suppression of numerical oscillations in the EMTP power systems[J]．IEEE Trans on Power Systems，1989，4（2）：739-747．

[11]Kuffel P，Kent K，Irwin G．The implementation and effectiveness of linear interpolation within digital simulation[J]．International Journal of Electrical Power and Energy Systems，1997，19（4）：221-227．

[12]Zou Ming，Mahseredjian Jean，Joos Geza，et al．Interpolation and reinitialization in time-domain simulation of power electronics circuits[J]．Electric Power Systems Research，2006，76（8）：688-694．

[13]王成山，李鵬，黃碧斌，等（Wang Chengshan，Li Peng，Huang Bibin，et al）．一種計及多重開關的電力電子時域仿真插值算法（An interpolation algorithm for time-domain simulation of power electronics circuit considering multiple switching events）[J]．電工技術學報（Transactions of China Electrotechnical Society），2010，25（6）：83-88．

[14]Hong S N，Liu C R，Bo Z Q，et al．Elimination of numerical oscillation of dynamic phasor in HVDC system simulation [C]//IEEE Power&Energy Society General Meeting.Calgary，Canada：2009．

[15]王勇，趙錦成，解璞（Wang Yong，Zhao Jincheng，Xie Pu）．基于PSCAD/EMTDC的開關元件的暫態仿真算法研究（Study on transient simulation algorithm of switching element based on PSCAD/EMTDC）[J]．電氣開關（Electric Switchgear），2008，5（5）：17-18．

[16]Faruque M O，Dinavahi Venkata，Xu Wilsun．Algorithms for the accounting of multiple switching events in digital simulation of power-electronic systems[J]．IEEE Trans on Power Delivery，2005，20（21）：1157-1167．

Interpolation Algorithm for Electromagnetic Transient Simulation Considering Multiple Switching Events

LIU Tao1，YAN Yi-peng2，JIN Na3，WANG Xu-dong1，LIN Ji-keng4
（1.Electrical Power Research Institute of Electrical Company of Tianjin，Tianjin 300384，China；2．Qingpu Power Supply Branch of Shanghai Municipal Electric Power Company，Shanghai 201700，China；3.Jinshan Power Supply Branch of Shanghai Municipal Electric Power Company，Shanghai 201500，China；4.Department of Electrical Engineering，College of Electronic and Information Engineering，Tongji University，Shanghai 201804，China）

Due to the iteration oscillation and low precision frequently resulted from multiple switching events in electromagnetic transient program，a new interpolation algorithm to solve these problems is presented in the paper.Interpolation technology is applied to find the accurate switching time.The trial backward Euler integration method，considering multiple switching action during the simulation，is then utilized to integrate interval followed the latest action from switching action in current step.Extrapolation and back-ward technology are combined to obtain the variable value after switching points.Thus a correct topology and initial values for the next integration are achieved.Finally，a step of trapezoidal integration combined with interpolating is further proceeded to be back to the integer step point.The results of the samples demonstrate the validity of the proposed algorithm.

interpolation technology；electromagnetic transient simulation；numerical oscillation；multiple switching

TM711

A

1003-8930（2013）06-0143-05

劉 濤（1980—），碩士，工程師，從事電力系統仿真及繼電保護試驗工作。Email：net_liutao@126.com

2013-04-22；

2013-05-18

閻貽鵬（1987—），碩士，助理工程師，從事電力系統仿真分析工作。Email：yyp1125018@126.com

金 娜（1986—），碩士，助理工程師，從事電力系統仿真分析工作。Email：Jinna198777@163.com