李 珂，顧 偉，柳 偉，曹 陽，樓冠男，鄒德虎

（1. 東南大學電氣工程學院，江蘇省 南京市 210096；2. 南京理工大學自動化學院，江蘇省 南京市 210094）

0 引言

隨著分布式發(fā)電和交直流混聯(lián)等技術的發(fā)展，電網(wǎng)的電力電子化程度不斷提高，網(wǎng)絡中功率流動和暫態(tài)過程的分析更加復雜［1］，需要模型更完善、計算精度更高的仿真系統(tǒng)對網(wǎng)絡特性進行研究。以變流器為代表的電力電子設備越來越多地融入電力系統(tǒng)中，其中部分開關器件的工作頻率可高達104Hz以上，需要微秒級的仿真步長反映其暫態(tài)過程，大大增加了仿真系統(tǒng)的計算負擔。

在仿真計算效率和計算機計算能力有限的前提下，實時化電磁暫態(tài)仿真的仿真規(guī)模、仿真精度和計算速度難以得到兼顧。為了滿足變流器模型中快動態(tài)的開關元件對微秒級仿真步長的需求，RT-LAB等商業(yè)仿真平臺從底層硬件的角度對仿真系統(tǒng)進行了 改 進［2-3］，通 過 引 入 現(xiàn) 場 可 編 程 門 陣 列（field programmable gate array，F(xiàn)PGA）進一步優(yōu)化系統(tǒng)的計算性能，以實現(xiàn)電磁暫態(tài)仿真的實時化。文獻［4-5］采用一種基于伴隨離散電路（associated discrete circuit，ADC）開關模型的電磁暫態(tài)實時仿真計算方法，利用該模型下導納矩陣不隨開關狀態(tài)變化的特性在FPGA 中進行實時仿真計算。文獻［6］基于ADC 開關模型實現(xiàn)了光伏發(fā)電系統(tǒng)的暫態(tài)仿真，并提出了控制器、光伏電池等元件在FPGA 中的建模方法，對仿真步長與結果誤差的關系進行了研究。文獻［7-8］通過計算所有開關可能的狀態(tài)組合并預存每種狀態(tài)下的系統(tǒng)導納矩陣，反映了開關狀態(tài)變化對系統(tǒng)的影響。以上基于FPGA 獨立實現(xiàn)的電磁暫態(tài)仿真存在以下問題：

1）采用ADC 開關模型帶來了虛擬功率損耗問題［9］，通過預存導納矩陣反映開關狀態(tài)時程序占用的存儲空間隨開關元件的數(shù)量迅速上升［10］；

2）FPGA 成本較高且硬件資源有限，一定程度上限制了仿真網(wǎng)絡的規(guī)模；

3）網(wǎng)絡控制策略和系統(tǒng)參數(shù)的調整過程較為復雜，對實時仿真的靈活性造成了限制。

文獻［11］提出一種基于CPU-FPGA 異構平臺的逆變器實時仿真方案。該方案結合不同計算單元的優(yōu)勢，對變流器模型進行解耦后將控制系統(tǒng)的計算由基于CPU 的服務器實現(xiàn)，其余計算由FPGA 實現(xiàn)，在保證實時性的前提下降低了FPGA 的硬件資源消耗，但在不同步長時序配合的過程中，直接以對側最后一個仿真步長的結果作為輸入，沒有進一步考慮多速率解耦并行仿真過程中的時序誤差。本文提出一種變流器電磁暫態(tài)硬件加速實時仿真方法，通過對并行仿真時序下控制系統(tǒng)解耦的多速率仿真誤差進行分析，采用結合拉格朗日插值的改進線性插值算法得到更精確的計算結果，并基于FPGA 實現(xiàn)了電磁暫態(tài)實時化硬件加速仿真計算，建立了計算能力更強的實時仿真系統(tǒng)。最后，本文搭建了三相橋式逆變器模型算例對所提出仿真方法的性能進行分析，并將仿真結果與PSCAD 離線仿真結果進行對比，驗證了實時仿真結果的準確性。

1 控制系統(tǒng)解耦的多速率并行計算方法

與工作狀態(tài)高頻率變化的開關元件相比，變流器控制系統(tǒng)的動態(tài)特性更接近于變流器以外的慢動態(tài)網(wǎng)絡，實現(xiàn)控制系統(tǒng)的多速率解耦仿真可以在反映高頻暫態(tài)過程的同時有效提高仿真計算效率。本章重點介紹變流器控制系統(tǒng)與電氣狀態(tài)解耦求解的多速率并行計算方法。

1.1 多速率并行計算時序分析

實現(xiàn)變流器控制系統(tǒng)與電氣狀態(tài)解耦計算時，需要將離線仿真中的控制系統(tǒng)與電氣狀態(tài)從以單一步長串行求解的仿真時序轉變?yōu)橐圆煌介L并行求解的仿真時序。圖1 反映了不同電磁暫態(tài)計算時序的對比，圖中的箭頭表示依次發(fā)生的數(shù)據(jù)計算或傳輸過程。

圖1 電磁暫態(tài)計算時序示意圖Fig.1 Schematic diagram of time sequence of electromagnetic transient calculation

設變流器控制系統(tǒng)采用較大仿真步長ΔT進行仿真，而電氣狀態(tài)采用小步長Δt求解，其中ΔT為Δt的整數(shù)倍，定義其倍數(shù)δ=ΔT/Δt為步長倍率（以δ=4 為例），將大小步長之間的數(shù)據(jù)交互設置在每個大步長的起始時刻進行?？梢钥闯?，仿真時序的轉變將在兩方面產(chǎn)生計算誤差：

1）控制系統(tǒng)由串行計算時序轉變?yōu)椴⑿杏嬎銜r序增加的延遲誤差［12］。對比圖1（a）和（b）中t+Δt到t+2Δt時刻的電氣狀態(tài)求解過程可以發(fā)現(xiàn)，單一速率串行時序下，該時間段求解時使用的控制參數(shù)以t+Δt時刻的電氣狀態(tài)為基礎在過程5 中計算得到，經(jīng)過程6 傳遞到電氣系統(tǒng)側；而在單一速率并行時序下使用的控制參數(shù)根據(jù)t時刻的電氣狀態(tài)在過程1 中得到，經(jīng)過程2 傳遞到電氣系統(tǒng)側，由此控制系統(tǒng)的計算產(chǎn)生一個步長的額外延遲，對結果的準確性造成影響。

2）單一速率仿真轉變?yōu)槎嗨俾史抡娈a(chǎn)生的延遲誤差。對比圖1（a）和（c）可知，采用多速率串行時序替代單一速率串行時序時，大步長仿真計算的頻次低于小步長仿真，因此，并非每個小步長仿真計算采用的大步長側輸入量都是對應時刻的實際值。即多速率仿真過程中，變流器控制系統(tǒng)的變化對于電氣系統(tǒng)計算的影響將無法得到完全反映，由此產(chǎn)生了相應的計算誤差。當不同時間常數(shù)的系統(tǒng)之間動態(tài)特性差異更小時，或者大小步長之間的倍數(shù)更大時，該問題產(chǎn)生的誤差將更加顯著。

1.2 結合拉格朗日插值的改進線性插值算法

產(chǎn)生以上誤差的根源在于，并行時序下多速率仿真計算中對側的實時數(shù)據(jù)難以得到準確表示，需要結合插值算法對對側實時數(shù)據(jù)進行有效預測，補償控制系統(tǒng)在一個大步長對應的多個小步長下的變化，減小多速率并行計算過程中的誤差。文獻［13］采用線性內插計算一個大步長所對應的多個小步長時刻慢動態(tài)網(wǎng)絡的歷史狀態(tài)，解決了多速率網(wǎng)絡解耦計算時，小步長側歷史電流源計算所需的前一時刻慢動態(tài)網(wǎng)絡狀態(tài)未知的問題，但內插法只能用于已知區(qū)間內，無法將插值范圍擴大到未計算的時刻。外插法是一種根據(jù)已知的離散數(shù)據(jù)集合對范圍外的數(shù)據(jù)進行預測的插值方法，利用先前時刻的已知結果通過插值對未來時刻的結果進行估計可達到減小延遲誤差的效果。文獻［14］采用單次線性外插預測電氣狀態(tài)，降低了計算時序由串行轉變?yōu)椴⑿兴黾拥难舆t誤差。這種方法具有形式簡潔、計算量小的優(yōu)勢，在步長極小的情況下具有較好的數(shù)值穩(wěn)定性［15］。但文獻［14］沒有對多速率仿真場景下的誤差補償做進一步分析，在多速率仿真步長倍率較大時準確性有所下降。

本文提出一種結合拉格朗日插值的改進線性插值算法。首先，采用兩次迭代外插擴大已知區(qū)間，再使用拉格朗日內插法得到下一大步長所對應的各小步長時刻下控制系統(tǒng)的狀態(tài)，并以此為依據(jù)計算開關狀態(tài)，在下一大步長開始時將結果傳遞至電氣狀態(tài)求解側，對控制系統(tǒng)解耦后的多速率仿真及并行仿真計算過程中的兩方面誤差做進一步補償，最后得到更準確的計算結果。

圖2 反映了本文所提出改進線性插值算法的計算時序。圖中，黃色部分和紅色部分分別代表t―ΔT到t時刻，以及t到t+ΔT時刻變流器控制系統(tǒng)的計算過程，藍色部分代表t到t+ΔT時刻電氣狀態(tài)的求解過程。

圖2 改進線性插值算法時序示意圖Fig.2 Schematic diagram of time sequence of improved linear extrapolation algorithm

從t到t+ΔT時刻，控制系統(tǒng)計算過程具體如下：

步驟1：接收t-ΔT到t時刻電氣系統(tǒng)的計算結果。

步驟2：根據(jù)t時刻的電氣狀態(tài)計算該時刻下控制系統(tǒng)的輸出vmod（t）。

步驟3：基于vmod（t）以及先前t-ΔT和t-2ΔT時 刻 的 控 制 系 統(tǒng) 輸 出 量vmod(t-ΔT)、vmod(t-2ΔT)，采 用 兩 次 迭 代 插 值 對t+ΔT和t+2ΔT時刻的控制系統(tǒng)輸出量vmod（t+ΔT）、vmod（t+2ΔT）進行預測，如式（1）所示。

步驟4：結合預測結果vmod（t+ΔT）和vmod（t+2ΔT），利用拉格朗日插值法對t+ΔT到t+2ΔT時刻之間每個小步長時刻下控制系統(tǒng)的輸出量進行預測，即

式中：ti為t+ΔT到t+2ΔT時刻之間的小步長時刻，ti=t+ΔT+nΔt，n=0，1，…，δ-1；lk（ti）為ti時刻插值節(jié)點k對應插值基函數(shù)的值；tλ為插值區(qū)間的終點t+2ΔT；K為插值節(jié)點的個數(shù)。

步驟5：根據(jù)步驟4 的結果計算t+ΔT到t+2ΔT時刻的開關狀態(tài)，并于t+ΔT時刻向小步長側傳遞。

與此同時，小步長側在t時刻接收大步長側于t-ΔT到t時刻計算出的t到t+ΔT時刻的開關狀態(tài)，并以此為依據(jù)對t到t+ΔT時刻的電氣狀態(tài)進行求解，實現(xiàn)與大步長側控制系統(tǒng)的并行計算。

將以上方法用于變流器控制系統(tǒng)與電氣系統(tǒng)進行多速率解耦并行計算的實時仿真場景?？紤]到控制系統(tǒng)的計算步長遠大于電氣狀態(tài)求解，因此，改進插值算法增加的計算量不會對仿真系統(tǒng)的實時性產(chǎn)生明顯負擔。

2 基于FPGA 的硬件加速實時仿真實現(xiàn)

在大規(guī)模系統(tǒng)中，全部采用微秒級小步長進行電磁暫態(tài)仿真的效率是難以接受的［16］。因此，基于前文所介紹的改進插值算法，本文將電磁暫態(tài)仿真模型劃分到不同的計算單元中，利用FPGA 并行性更強的優(yōu)勢處理小步長電氣狀態(tài)求解，使用基于CPU 的上位機處理規(guī)模更大、邏輯更復雜的大步長仿真計算，通過PCIe 總線完成不同步長仿真之間的數(shù)據(jù)交互，實現(xiàn)基于FPGA 的硬件加速實時仿真。具體結構見附錄A 圖A1。

仿真系統(tǒng)的計算時序及FPGA 部分的主要計算過程見附錄A 圖A2。本文采用文獻［17］中的算法實現(xiàn)電磁暫態(tài)仿真?？紤]到開關不同運行狀態(tài)下導納矩陣的變化將產(chǎn)生巨大的運算負擔［18-19］，對于狀態(tài)高頻率變化的開關元件采用基于交叉初始化的變流器參數(shù)化恒導納模型［20］進行建模求解。該模型下系統(tǒng)導納矩陣恒定不變，可以將變流器導納矩陣及相應逆矩陣在初始化過程中提前計算并存儲。大步長開始時刻，上位機首先通過PCIe 總線發(fā)送上一步長下計算得到的接口狀態(tài)、控制參數(shù)等信息到FPGA 端，然后，啟動本時刻的變流器控制系統(tǒng)計算。與此同時，F(xiàn)PGA 端的電氣狀態(tài)求解過程主要包括以下步驟：

步驟1：在初始化階段，利用導納矩陣不隨開關狀態(tài)改變的特性分塊預存變流器模型導納矩陣及相應的逆矩陣。

步驟2：讀取控制系統(tǒng)前一步長計算得到的各開關狀態(tài)及接口狀態(tài)。

步驟3：結合開關狀態(tài)對各支路歷史電流源在t時刻的電流Ihis（t）進行計算。

式中：Geq為支路等效電導；ubr(t-Δt)和ibr(t-Δt)分別為t-Δt時刻支路兩端的電壓和流過支路的電流；對于電阻、電感、電容等狀態(tài)固定的支路，α、β均為定值［21］，對于開關支路，α、β的值取決于本時刻的開關狀態(tài)［20］。

步驟4：根據(jù)系統(tǒng)電源以及各歷史電流源的值計算t時刻各節(jié)點的注入電流，例如對于節(jié)點i及其相鄰節(jié)點的集合w（i），有

步驟5：計算并更新各節(jié)點電壓、各支路電流等電氣狀態(tài)量。

式中：U（t）和Iinj（t）分別為t時刻各節(jié)點電壓和注入電流組成的向量；Y為系統(tǒng)導納矩陣。

步驟6：更新仿真時間，若未達到結束時刻則重復步驟2—步驟6；若達到則將計算結果存儲到指定地址等待上位機讀取。

在與以上過程并行的大步長結束前，上位機將讀取變流器電氣狀態(tài)的計算結果，完成基于FPGA的硬件加速實時仿真計算。

3 仿真分析

3.1 仿真精度分析

基于以上仿真方法，本文搭建了基于FPGA 的變流器硬件加速實時仿真系統(tǒng)。其中，CPU 端處理器型號為Intel Core i7-4770，核心數(shù)為8，主頻為3.40 GHz；FPGA 芯 片 型 號 為Kintex Ultrascale XCKU060，采用PCIe3.0 接口與CPU 端上位機進行連接。同時，建立圖3 所示的變流器模型作為算例，設CPU 端控制系統(tǒng)求解步長為100 μs，F(xiàn)PGA 端電氣狀態(tài)求解步長為1 μs，控制系統(tǒng)采用恒功率控制，電氣及控制參數(shù)如表1 所示。圖中：S1—S6為變流器開關元件；uabc、udq和iabc、idq分別為dq變換前、后的變流器交流側電壓和電流；Udc為變流器直流側電壓；f和θ分別為交流側電壓頻率和相位；Pgrid和Qgrid分別為變流器輸出的有功和無功功率；Pref和Qref分別為恒功率控制的有功和無功功率參考值。

表1 變流器算例參數(shù)設置Table 1 Parameter setting of converter example

圖3 三相橋式變流器算例系統(tǒng)示意圖Fig.3 Schematic diagram of example system of three-phase bridge-type converter

3.1.1 仿真結果

為了驗證本文所提出仿真方法的有效性，分別采用基于改進插值法的多速率并行時序（EP-ON）、無插值預測的多速率并行時序（EP-OFF）以及單一速率并行時序（SP），在以下2 種場景下對仿真平臺進行測試，并在采用單一速率串行時序計算的離線仿真平臺PSCAD 中搭建等效模型進行步長為1 μs的仿真作為對照。

1）場 景A：Pref和Qref分 別 為50 kW 和0 kvar，5.5 s 時交流側BC 兩相發(fā)生接地短路故障，0.1 s 后故障消除，得到仿真結果如圖4 所示。

2）場景B：控制系統(tǒng)Pref初始設定為50 kW，Qref初始設定為0 kvar。5 s 時Pref突變?yōu)?0 kW，Qref突變?yōu)?0 kvar，得到仿真結果如圖5 所示。

圖4、圖5 分別反映了在發(fā)生短路故障與控制參數(shù)變化的情況下，不同仿真計算方法所得到系統(tǒng)暫態(tài)過程的對比情況。圖中，P、Q和Ia和Udc分別代表變流器交流側輸出的有功功率、無功功率A 相交流電流和直流側電壓；EP-ON、EP-OFF 和SP 分別代表改進插值方法和無插值方法下的多速率并行仿真和單一速率并行仿真結果。表2 給出了幾種算法相對于離線仿真結果中變流器輸出有功功率Pout和A相交流電流IacA的均方根誤差，其中分別以發(fā)生故障或參數(shù)調整時刻前后0.2 s 內的仿真結果為樣本，反映穩(wěn)態(tài)及暫態(tài)情況下幾種算法的誤差大小?？梢钥闯觯捎脝我凰俾什⑿袝r序得到的仿真結果與PSCAD 得到的結果極為接近，證明了本文所采用仿真模型的準確性。在穩(wěn)態(tài)情況下，幾種算法得到的仿真結果基本一致，但系統(tǒng)參數(shù)調整或發(fā)生故障后，無插值預測的多速率并行仿真結果與PSCAD離線仿真結果相比出現(xiàn)了一定偏差，而采用改進插值法的多速率并行仿真得到的結果與離線仿真結果基本一致，保持了較高的仿真精度，驗證了本文所提出算法對于降低計算時序由單一速率串行轉變?yōu)槎嗨俾什⑿袝r所產(chǎn)生時序誤差的有效性。

表2 不同算法的均方根誤差對比Table 2 Comparison of root mean square errors of different methods

圖4 場景A 下的仿真波形對比Fig.4 Comparison of simulation waveforms in scenario A

圖5 場景B 下的仿真波形對比Fig.5 Comparison of simulation waveforms in scenario B

3.1.2 步長變化影響分析

本文所提出的改進插值算法用于降低多速率實時仿真場景下計算時序轉變產(chǎn)生的時序誤差，仿真步長的變化將對算法效果產(chǎn)生一定的影響。當仿真步長擴大時，場景A 下變流器交流側輸出電流的有效值相對于PSCAD 離線仿真結果的相對誤差如附錄B 圖B1 所示。其中，SE 代表文獻［14］所提出外插方法下的計算結果；設步長倍率δ固定為50，電氣狀態(tài)求解步長分別為1、2、3 μs。可以看出，在穩(wěn)態(tài)條件下采用插值算法與無插值預測得到的仿真結果較為接近，但故障發(fā)生后的暫態(tài)過程中，無插值算法的多速率并行仿真時序下計算誤差隨步長的增加快速上升，而采用插值算法后誤差的增長得到明顯抑制，且本文所提出的改進插值算法對于計算精度的改善效果更好，驗證了本文所提出仿真方法的有效性。

當電氣狀態(tài)求解步長固定時，時序變化引起的誤差隨步長倍率δ的增加而上升，而步長倍率取決于仿真場景的實際需求。在電力電子化程度較高的仿真場景下，可以采用微秒甚至納秒級的仿真步長反映開關的高頻特性，而配電網(wǎng)采用微秒或毫秒級的步長模擬［22］，需要根據(jù)網(wǎng)絡模型的動態(tài)特性合理選擇步長倍率。以場景B 為例，討論δ變化對于本文所提出仿真方法的影響。設電氣狀態(tài)求解步長為3 μs，步長倍率δ分別為20、50 和100 時變流器交流側B 相輸出電流Ib與PSCAD 離線仿真結果的對比見附錄B 圖B2。可以看出，控制參數(shù)Pref與Qref調整前的穩(wěn)態(tài)情況下，無論是否采用插值算法得到的結果都與離線仿真接近，但參數(shù)調整后的暫態(tài)過程中，無插值算法的多速率并行計算產(chǎn)生明顯誤差，誤差收斂的速度與δ負相關。與文獻［14］中的方法相比，改進插值算法對該誤差表現(xiàn)出更好的抑制效果，使計算精度得到了顯著改善。

3.2 計算性能與硬件資源消耗分析

實時仿真要求在每個步長下，仿真計算消耗的總時間低于仿真步長。在本文所建立的仿真平臺上，上位機向FPGA 寫入控制信息的平均耗時為10.391 μs，從FPGA 讀 取 電 氣 狀 態(tài) 的 平 均 耗 時 為10.446 μs，數(shù)據(jù)傳輸?shù)目倳r間為20.837 μs。由于控制系統(tǒng)與電氣狀態(tài)并行求解且控制系統(tǒng)的求解次數(shù)遠小于電氣狀態(tài)求解，制約仿真平臺實時性能的主要因素在于FPGA 側的求解速度。以上算例中，F(xiàn)PGA 計算模塊時鐘頻率為142.8 MHz，變流器電氣狀態(tài)求解過程中每個環(huán)節(jié)的計算耗時如表3 所示。其中，100 次小步長計算的總耗時47.747 μs 與數(shù)據(jù)傳輸?shù)目倳r間20.837 μs 之和為68.584 μs，小于控制系統(tǒng)計算步長100 μs 且仍有充足裕度，有實現(xiàn)更大規(guī)模實時仿真計算的條件。在本文提出的硬件加速仿真方法中，CPU 端無須承擔變流器電氣狀態(tài)求解任務，且解耦后的控制系統(tǒng)的求解次數(shù)相比串行時序顯著減少，大大降低了CPU 端的仿真計算負擔，提高了系統(tǒng)的實時仿真能力。

表3 仿真實時性能分析Table 3 Analysis on real-time simulation performance

在計算不變的情況下，高并行度計算通常需要占用更多的硬件資源。表4 所示為本文提出的硬件加速仿真結構與FPGA 獨立實現(xiàn)仿真時，核心計算模塊的主要硬件資源消耗情況對比。

表4 硬件資源消耗分析Table 4 Analysis on hardware resource consumption

表4 中，F(xiàn)PGA-S 表示FPGA 獨立以單一速率并行時序實現(xiàn)仿真時的資源消耗；COSIM 表示在本文建立的硬件加速仿真平臺中FPGA 的資源消耗情況。與本文所述硬件加速結構相比，F(xiàn)PGA 獨立實現(xiàn)仿真時需計算控制系統(tǒng)輸出，增加了并行程度較低、乘加運算較多的計算流程，故DSP48E 的消耗上升顯著。與FPGA 獨立實現(xiàn)的電磁暫態(tài)仿真相比，本文提出的硬件加速計算結構在保持高計算精度的情況下，有效降低了FPGA 的硬件資源消耗。

4 結語

本文提出了一種通過控制系統(tǒng)解耦實現(xiàn)多速率仿真的變流器電磁暫態(tài)實時仿真方法，分析了控制系統(tǒng)解耦計算過程中的時序和誤差問題，采用一種結合拉格朗日插值的改進線性插值算法進行誤差補償，有效提高了控制器多速率解耦仿真的計算精度，并以該算法為依據(jù)搭建了利用FPGA 實現(xiàn)硬件加速計算的實時化仿真平臺。將CPU 端與FPGA 端仿真步長分別設為100 μs 與1 μs，搭建采用功率控制的逆變器模型在系統(tǒng)故障和控制參數(shù)調整的場景下，分別與無插值算法的多速率仿真和PSCAD 離線仿真的結果進行對比。

仿真結果表明，本文提出的改進線性插值算法在精度上優(yōu)于無插值算法的控制器多速率解耦仿真，更接近離線仿真得到的計算結果。此外，通過仿真計算耗時和硬件資源消耗情況分析，驗證了本文所提出的變流器電磁暫態(tài)硬件加速實時仿真方法在保持高計算精度和實時性能的前提下，實現(xiàn)了更少的硬件資源消耗，具有良好的可行性和有效性。后續(xù)研究考慮建立種類更豐富的電磁暫態(tài)模型，實現(xiàn)多類型器件的大規(guī)模仿真。

附錄見本刊網(wǎng)絡版（http：//www.aeps-info.com/aeps/ch/index.aspx），掃英文摘要后二維碼可以閱讀網(wǎng)絡全文。