施 文，遲 頌，郭 亮

（1.河北工業(yè)大學(xué)電氣工程學(xué)院省部共建電工裝備可靠性與智能化國家重點(diǎn)實(shí)驗(yàn)室，天津300132；2.河北工業(yè)大學(xué)電氣工程學(xué)院河北省電磁場與電器可靠性重點(diǎn)實(shí)驗(yàn)室，天津300132）

以實(shí)時(shí)仿真為核心的數(shù)字雙胞胎技術(shù)，能夠幫助企業(yè)在實(shí)際投入生產(chǎn)之前在虛擬環(huán)境中優(yōu)化、仿真和測試，實(shí)現(xiàn)高效的柔性生產(chǎn)，是企業(yè)邁入工業(yè)4.0 的解決方案之一， 然而由于寬禁帶器件的高頻特性使得電力電子裝置與數(shù)字雙胞胎結(jié)合的難度增加。 寬禁帶器件開關(guān)頻率可達(dá)百kHz，最小穩(wěn)態(tài)時(shí)間可低于1 μs，而常見的基于現(xiàn)場可編程門陣列FPGA（field programmable gate array）實(shí)時(shí)仿真平臺RT-LAB、PXI 平臺等仿真步長在250 ns 到1 μs 之間[1-2]。過高的開關(guān)頻率會使得仿真中出現(xiàn)開關(guān)時(shí)刻偏移等問題，一般采用變步長計(jì)算。

文獻(xiàn)[3-5]提出了實(shí)時(shí)仿真并行化的方法，證明了在并行化實(shí)時(shí)仿真中，可以異步仿真，無須考慮與通訊部分同步， 每部分仿真模型可以將仿真步長減小至算法極限。因此，減少仿真計(jì)算時(shí)間可以減小仿真步長，有助于在高頻下簡化開關(guān)動(dòng)作情況，提高仿真精度。算法方面，文獻(xiàn)[6-7]利用多步插值的變步長算法，解決了多重開關(guān)動(dòng)作的問題，但算法復(fù)雜難以實(shí)時(shí)計(jì)算；文獻(xiàn)[8]比較了常用離散迭代算法仿真精度，證明了梯形法綜合性能更優(yōu)；文獻(xiàn)[9]使用權(quán)重?cái)?shù)值積分的方法解決了變步長仿真中的數(shù)值發(fā)散問題，但都沒有涉及實(shí)時(shí)仿真；文獻(xiàn)[10]通過單步插值實(shí)現(xiàn)了實(shí)時(shí)仿真變步長算法；文獻(xiàn)[11]通過二階段積分和自校正單步插值的方法解決了變步長的相關(guān)問題，但都未考慮計(jì)算時(shí)間對仿真步長的影響。

針對現(xiàn)有實(shí)時(shí)仿真中使用的變步長算法計(jì)算速度和數(shù)值發(fā)散的問題，提出一種基于FPGA 程序結(jié)構(gòu)的實(shí)時(shí)仿真數(shù)值算法，通過改進(jìn)梯形迭代法以提高算法收斂程度， 在此基礎(chǔ)上改進(jìn)仿真算法結(jié)構(gòu)，使計(jì)算更加快速穩(wěn)定。以三相逆變電路為例，以離線仿真軟件精度為標(biāo)準(zhǔn)，對改進(jìn)算法進(jìn)行實(shí)時(shí)仿真，證明算法的有效性。

1 高頻電力電子實(shí)時(shí)仿真存在問題

隨著電力電子系統(tǒng)開關(guān)頻率的不斷提高，仿真中會出現(xiàn)開關(guān)動(dòng)作時(shí)刻與采樣點(diǎn)時(shí)刻不一致的情況， 這時(shí)開關(guān)動(dòng)作時(shí)刻會被歸算到采樣點(diǎn)時(shí)刻，出現(xiàn)開關(guān)動(dòng)作偏移的問題。一般采用變步長算法進(jìn)行處理，但隨之帶來了關(guān)于數(shù)值發(fā)散[11]與實(shí)時(shí)性的問題。 文獻(xiàn)[12]總結(jié)了變步長算法及適用情況，文獻(xiàn)[13]對開關(guān)動(dòng)作偏移問題進(jìn)行了詳細(xì)的描述。

（1）數(shù)值發(fā)散的問題。 由于步長過大或者算法精度不足導(dǎo)致仿真結(jié)果發(fā)散，通常可以使用高階迭代算法、多步迭代等方法解決，但在實(shí)時(shí)仿真中會增加計(jì)算時(shí)間，進(jìn)而增大步長，使問題更加復(fù)雜，甚至無法保證實(shí)時(shí)。

（2）實(shí)時(shí)性問題。 實(shí)時(shí)仿真最重要的是滿足實(shí)時(shí)性，即：①實(shí)時(shí)性，步長設(shè)定時(shí)間等于數(shù)據(jù)傳輸間隔時(shí)間；②穩(wěn)定性，最長仿真計(jì)算時(shí)間小于步長設(shè)定時(shí)間。

常用的高階迭代算法會增加總體仿真計(jì)算時(shí)間，變步長仿真中采用半步長歐拉法仿真增加最長仿真計(jì)算時(shí)間，從而變相增大了步長。

本文在實(shí)時(shí)仿真變步長算法的基礎(chǔ)上，針對問題（1）提出改進(jìn)的梯形法以提高收斂性，在此基礎(chǔ)上針對問題（2）改進(jìn)算法計(jì)算結(jié)構(gòu)，使計(jì)算更加穩(wěn)定快速，可使仿真步長進(jìn)一步減小。

2 基于系數(shù)控制的數(shù)值解算方法

2.1 改進(jìn)的梯形迭代算法

一般實(shí)時(shí)仿真主要使用歐拉法迭代配合狀態(tài)方程對離散電路模型進(jìn)行求解，公式表示為

其中：第1 式為狀態(tài)方程，第2 式為歐拉法公式。式中：h 為仿真步長；X 為狀態(tài)變量；為狀態(tài)變量的導(dǎo)數(shù)；A 和B 為狀態(tài)方程的數(shù)值矩陣； 下標(biāo)n 代表離散時(shí)間點(diǎn)。

歐拉法精度較低， 仿真結(jié)果易出現(xiàn)數(shù)值發(fā)散，可以使用梯形法來提高精度[8-9]。 同樣對于狀態(tài)變量的導(dǎo)數(shù)求解使用狀態(tài)變量法，梯形迭代法公式為

式中，X'為對下一時(shí)刻狀態(tài)變量的預(yù)測值。 在實(shí)時(shí)仿真中，預(yù)測計(jì)算會延長仿真計(jì)算時(shí)間，導(dǎo)致仿真步長增大。 通過圖1 進(jìn)行說明。

圖1 梯形迭代法計(jì)算流程Fig. 1 Calculation flow diagram of trapezoidal iterative algorithm

圖1 為梯形迭代法的簡化計(jì)算流程，其中細(xì)箭頭表示常數(shù)乘法，1/2 表示系數(shù)數(shù)值，未寫數(shù)值則系數(shù)為1； 而粗箭頭代表狀態(tài)方程AX+B 的計(jì)算，最為消耗計(jì)算時(shí)間，省略了步長h 乘法計(jì)算。 梯形迭代法會進(jìn)行2 次狀態(tài)方程計(jì)算， 存在時(shí)序問題，無法并行計(jì)算。相比于歐拉法，計(jì)算時(shí)間增長近1 倍。

為了減少計(jì)算時(shí)間，對梯形法進(jìn)行一定改進(jìn)。圖1 中，tn到tn+1中第一次狀態(tài)方程的計(jì)算結(jié)果，與tn-1到tn的基于歐拉法的預(yù)測計(jì)算結(jié)果在一定條件下近似等效。 將替換為，可以節(jié)省一次狀態(tài)方程計(jì)算，改進(jìn)的梯形迭代法計(jì)算流程如圖2 所示。

圖2 改進(jìn)的梯形迭代法計(jì)算流程Fig. 2 Calculation flow diagram of improved trapezoidal iterative algorithm

梯形迭代法從tn-1到tn的仿真計(jì)算中，存在使用歐拉法對tn的狀態(tài)變量導(dǎo)數(shù)進(jìn)行預(yù)測計(jì)算結(jié)果，將其運(yùn)用在tn對tn+1的仿真計(jì)算中，如圖2 虛線部分所示，整理后得到

相對于傳統(tǒng)梯形迭代法，改進(jìn)迭代法計(jì)算省去一次狀態(tài)方程的計(jì)算，速度更快。

精度方面， 以buck 電路為例對比算法誤差進(jìn)行說明。 電路拓?fù)淙鐖D3 所示。 其中，Es為電源電壓，取100 V；電感L 為0.1 mH；負(fù)載電阻Rload為2.5 Ω；占空比取0.5；開關(guān)頻率為10 kHz。對比歐拉法、 改進(jìn)梯形法與梯形法低步長的仿真結(jié)果，以PLECS 仿真結(jié)果為標(biāo)準(zhǔn)，計(jì)算誤差[14]公式為

式中：δerror為相對范數(shù)誤差；x 為標(biāo)準(zhǔn)結(jié)果；y 為對比數(shù)據(jù)。 誤差曲線如圖4 所示。

圖3 Buck 電路拓?fù)銯ig. 3 Topology of buck circuit

圖4 buck 電路輸出電壓誤差曲線Fig. 4 Error curves of buck circuit output voltage

表1 給出了本算例中不同算法預(yù)計(jì)計(jì)算時(shí)間，其中計(jì)算按1 個(gè)時(shí)鐘周期（clk）計(jì)時(shí)（使用定點(diǎn)數(shù)計(jì)算），開關(guān)判斷為1 clk，并行計(jì)算時(shí)間為最長路徑計(jì)算時(shí)間。

仿真精度方面， 改進(jìn)梯形法與梯形法相近，都遠(yuǎn)小于歐拉法；在計(jì)算時(shí)間方面，相較于梯形法，改進(jìn)梯形法時(shí)間約減少了45.5%。

表1 不同算法預(yù)計(jì)計(jì)算時(shí)間Tab. 1 Computation time predicted using different algorithms

2.2 變步長算法計(jì)算步驟

文獻(xiàn)[10]中的變步長算法，在完成插值后，將步長調(diào)整為原步長的一半， 使用歐拉法迭代計(jì)算2次，以此提高收斂性與計(jì)算精度，但在變步長時(shí)刻需要進(jìn)行2～3 次迭代計(jì)算，仿真時(shí)間較長。 為加快計(jì)算速度，改進(jìn)了變步長算法計(jì)算步驟[12]，計(jì)算流程如圖5 所示。

圖5 改進(jìn)的變步長算法計(jì)算流程Fig. 5 Calculation flow diagram of improved variablestep algorithm

步驟1從tn-1到tn時(shí)刻進(jìn)行計(jì)算， 并檢測到開關(guān)動(dòng)作；

步驟2對狀態(tài)變量插值，將狀態(tài)還原到開關(guān)時(shí)刻tsw；

步驟3步長變?yōu)閔+ΔT 進(jìn)行仿真計(jì)算， 仿真時(shí)間與計(jì)算時(shí)間相等，保證實(shí)時(shí)性。

步驟2 和步驟3 對應(yīng)的計(jì)算公式為

式中，ΔT 為圖5 中對應(yīng)的時(shí)間間隔。

在圖5 計(jì)算過程中， 經(jīng)過了2 次仿真計(jì)算，完成了2 個(gè)步長的仿真計(jì)算，可以保證每個(gè)步長只進(jìn)行一次迭代計(jì)算，仿真計(jì)算時(shí)間穩(wěn)定，單步長內(nèi)計(jì)算快速。

2.3 改進(jìn)的實(shí)時(shí)仿真數(shù)值計(jì)算方法

同時(shí)使用改進(jìn)的梯形法、歐拉法、改進(jìn)的變步長算法可以完成變步長計(jì)算，保證變步長計(jì)算中仿真結(jié)果的收斂性。 3 種方法都只進(jìn)行一次矩陣乘法，計(jì)算速度相近，計(jì)算快速。

但由于FPGA 屬于分布式處理器，每多一種情況都要預(yù)留相應(yīng)的計(jì)算資源，同時(shí)使用3 種算法使FPGA 程序?qū)崿F(xiàn)的難度提升， 因此需要對3 種算法（式（1）、式（3）、式（5））的計(jì)算結(jié)構(gòu)進(jìn)行改進(jìn)，減少不同的計(jì)算分支。參考權(quán)重?cái)?shù)值積分[9]的方法，將權(quán)重轉(zhuǎn)變?yōu)閹讉€(gè)確定的系數(shù)，通過控制系數(shù)進(jìn)行計(jì)算[10]，將算法的選擇轉(zhuǎn)變?yōu)橄禂?shù)的選擇，提升FPGA 程序效率。

以式（5）變步長算法的計(jì)算結(jié)構(gòu)為原型，進(jìn)行算法結(jié)構(gòu)的改進(jìn)；式（3）的改進(jìn)梯形法與變步長算法結(jié)構(gòu)相近，不做調(diào)整；式（1）歐拉迭代法向式（5）計(jì)算結(jié)構(gòu)進(jìn)行靠近，對式（1）進(jìn)行擴(kuò)充，增加2 個(gè)0系數(shù)的部分， 并將過程量等效為梯形法中的相等量，使計(jì)算結(jié)構(gòu)靠近式（5），最終得到

對式（6）、式（3）、式（5）相同部分進(jìn)行提取，不同部分使用系數(shù)a、b、c 代替得到

系數(shù)a、b、c 的不同取值及其對應(yīng)的算法如表2所示。 計(jì)算步驟如下。

步驟1通過開關(guān)動(dòng)作采樣程序， 得到上一個(gè)步長中的開關(guān)動(dòng)作情況，即是否開關(guān)動(dòng)作和動(dòng)作時(shí)間點(diǎn)ΔT。

步驟2根據(jù)步驟1 中開關(guān)動(dòng)作情況進(jìn)行系數(shù)的選擇，選擇原則如下：

（1）判斷開關(guān)是否動(dòng)作。若開關(guān)未動(dòng)作，使用改進(jìn)梯形法；若開關(guān)動(dòng)作，進(jìn)入下一步判斷。

（2）判斷開關(guān)動(dòng)作時(shí)間點(diǎn)ΔT。若ΔT>0，使用變步長算法；若ΔT=0，使用歐拉法。

步驟3根據(jù)選擇的系數(shù)進(jìn)行之后的計(jì)算。

除了只進(jìn)行一次矩陣乘法計(jì)算快速外，與切換算法相比， 切換系數(shù)對FPGA 的資源占用減小很多，便于程序編寫與仿真系統(tǒng)的建立；與使用多種算法相比，改進(jìn)的算法計(jì)算結(jié)構(gòu)在每個(gè)步長中的計(jì)算量一定，計(jì)算時(shí)間穩(wěn)定，易于與通信程序和采樣程序配合。

表2 不同算法對應(yīng)系數(shù)設(shè)置Tab. 2 Coefficient settings of different algorithms

3 仿真算法FPGA 程序結(jié)構(gòu)

3.1 實(shí)時(shí)仿真算法的FPGA 程序結(jié)構(gòu)

FPGA 程序硬件設(shè)計(jì)結(jié)構(gòu)如圖6 所示，圖中，實(shí)線部分為原有歐拉法定步長算法，虛線部分為本文算法增加部分，最下方為時(shí)間軸；圖6 為一個(gè)步長中涉及的計(jì)算，計(jì)算模塊是考慮并行運(yùn)算的基礎(chǔ)上根據(jù)計(jì)算的時(shí)間先后順序進(jìn)行排列，仿真計(jì)算時(shí)根據(jù)從左到右的順序進(jìn)行計(jì)算（對應(yīng)時(shí)間軸相同即為并行運(yùn)算）。

虛線部分改進(jìn)主要體現(xiàn)在以下3 個(gè)方面：

（2）與系數(shù)a 相關(guān)的算法切換與數(shù)值補(bǔ)償部分，可以與矩陣計(jì)算同時(shí)進(jìn)行，不增加計(jì)算時(shí)間；

（3）與系數(shù)b 相關(guān)的控制步長變換部分。

與傳統(tǒng)的定步長歐拉法計(jì)算對比，核心的矩陣計(jì)算時(shí)間不變，最長計(jì)算路徑增加2 個(gè)乘法與1 個(gè)加法。僅需提取開關(guān)信號并對狀態(tài)變量進(jìn)行處理即可，不需要對算法進(jìn)行根本的改動(dòng)，便于對已有程序的結(jié)構(gòu)進(jìn)行升級與優(yōu)化。

如圖6 所示，在資源使用方面，增加計(jì)算資源主要針對狀態(tài)變量的數(shù)量，正比于狀態(tài)方程矩陣階數(shù)；而整體計(jì)算資源近似正比于狀態(tài)方程矩陣階數(shù)二次方，即隨著矩陣階數(shù)增大，增加計(jì)算資源所占比例將逐漸減小。

圖6 FPGA 程序硬件結(jié)構(gòu)設(shè)計(jì)示意Fig. 6 Schematic of FPGA program hardware structure design

3.2 實(shí)時(shí)仿真算法的FPGA 程序結(jié)構(gòu)

由于仿真所選拓?fù)洳煌?、硬件不同，仿真時(shí)間不確定。 為保證對比簡潔，以第2.1 節(jié)中的算例為例進(jìn)行分析。 如圖5 所示，變步長的插值算法需要增加1 個(gè)乘法與1 個(gè)加法。以算例中的并行計(jì)算方式統(tǒng)計(jì)預(yù)計(jì)計(jì)算時(shí)間，如表3 所示。

在本算例中，計(jì)算時(shí)間可以節(jié)省72%，由于本文算法所占據(jù)的時(shí)間可并行，不隨矩陣階數(shù)增加而增加，節(jié)省時(shí)間比率相對穩(wěn)定。

表3 不同算法的計(jì)算時(shí)間Tab. 3 Computation time of different algorithms

4 實(shí)時(shí)仿真平臺實(shí)驗(yàn)驗(yàn)證

4.1 仿真程序?qū)嶒?yàn)平臺

本文使用的FPGA 開發(fā)板為Arria II GX 系列開發(fā)板，芯片型號為EP2AGX125EF35。 該FPGA 開發(fā)板提供124 100 個(gè)邏輯單元、576 個(gè)18×18 DSP乘法器，并配置頻率100 MHz 晶振器，可以提供足量的計(jì)算資源。

實(shí)驗(yàn)數(shù)據(jù)通過Altera 公司的開發(fā)套件Quartus Prime 中SignalTap II 程序， 設(shè)置探針程序?qū)PGA仿真結(jié)果進(jìn)行數(shù)據(jù)提取。

4.2 實(shí)驗(yàn)用仿真模型

仿真電路采用三相逆變器接LC 濾波器帶對稱阻感負(fù)載，電路拓?fù)淙鐖D7 所示，為方便程序計(jì)算規(guī)定點(diǎn)o 為參考節(jié)點(diǎn)。 電路模型使用狀態(tài)變量法。開關(guān)模型方面， 開通等效為電壓源串聯(lián)電阻模型（見圖7 中虛線），關(guān)斷等效為開路。 仿真電路參數(shù)如表4 所示。

圖7 三相逆變器電路拓?fù)銯ig. 7 Topology of three-phase inverter circuit

表4 仿真電路參數(shù)數(shù)據(jù)Tab. 4 Data of simulation circuit parameters

控制采用雙極性正弦脈寬調(diào)制SPWM（sinusoidal pulse width modulation）控制，正弦波通過直接數(shù)字頻率合成DDS（direct digital frequency synthesis）技術(shù)[15]將正弦波離散化后存入存儲器中，根據(jù)仿真步長調(diào)用數(shù)據(jù)。 受開發(fā)板存儲限制，為保證控制信號足夠精確，排除控制不同帶來的誤差，輸出設(shè)定為基頻250 Hz 的三相正弦波，LC 濾波器的截止頻率按照10 倍基頻整定。

算例主要使用定點(diǎn)數(shù)的數(shù)值進(jìn)行計(jì)算，本文算法計(jì)算使用時(shí)間為130 ns，步長設(shè)定在130 ns 以上即可保證實(shí)時(shí)性。 控制信號周期與步長同步，開關(guān)頻率取110 kHz。

4.3 仿真結(jié)果對比與分析

4.3.1 算法收斂性驗(yàn)證

為便于觀測發(fā)散現(xiàn)象，避免算法諧波誤差與數(shù)值發(fā)散的結(jié)果疊加，步長采用本文算法500 ns 與歐拉法500 ns 進(jìn)行仿真，仿真結(jié)果如圖8 所示。 本文在計(jì)算時(shí)間方面與歐拉法相當(dāng)，而在同樣步長500 ns 情況下， 本文算法的輸出波形是正常的正弦波形，而歐拉法出現(xiàn)了明顯的數(shù)值發(fā)散情況。

圖8 不同迭代算法的仿真結(jié)果比較Fig. 8 Comparison of simulation results between different iterative algorithms

4.3.2 算法仿真驗(yàn)證

由于并行仿真無需預(yù)留過多通信時(shí)間，仿真步長可以接近計(jì)算時(shí)間， 算例中設(shè)置步長為150 ns，大于計(jì)算時(shí)間130 ns，滿足實(shí)時(shí)性。 圖9 為本文算法步長150 ns 下，三相逆變器的三相輸出電壓（圖7 中為A 相uao、B 相ubo、C 相uco）與三相輸出電流（圖7 中為A 相iao、B 相ibo、C 相ico）。

由圖9 可見， 逆變器輸出電壓和輸出電流為平滑的三相正弦曲線，無諧波影響，本文算法的變步長部分計(jì)算正常，算法收斂性正常。仿真結(jié)果證明了算法的有效性。 將圖9 中的計(jì)算結(jié)果與PLECS 軟件變步長計(jì)算結(jié)果作對比，以A 相輸出電壓和輸出電流為例，結(jié)果如圖10 所示。使用式（4）對圖9 的三相輸出結(jié)果進(jìn)行誤差計(jì)算，計(jì)算結(jié)果如表5 所示。

圖9 三相逆變器輸出波形Fig. 9 Output waveforms of three-phase inverter

圖10 A 相輸出波形對比Fig. 10 Comparison of output waveform in phase A

表5 實(shí)時(shí)仿真結(jié)果誤差Tab. 5 Errors of real-time simulation results

結(jié)果表示誤差較小，仿真結(jié)果可用。 算法計(jì)算時(shí)間為130 ns，與一次歐拉法計(jì)算時(shí)間相近，而使用半步長歐拉法進(jìn)行變步長計(jì)算，需要約3 倍的歐拉法計(jì)算時(shí)間。與之相比，本文算法可以節(jié)省約2/3的計(jì)算時(shí)間，而由于計(jì)算快速，步長較小，誤差可以進(jìn)一步減小。

5 結(jié)論

針對變步長實(shí)時(shí)仿真給嵌入式實(shí)時(shí)仿真帶來的數(shù)值發(fā)散和計(jì)算時(shí)間不穩(wěn)定的問題，提出一種基于FPGA 程序結(jié)構(gòu)的實(shí)時(shí)仿真算法，優(yōu)勢為：①計(jì)算時(shí)間穩(wěn)定；②收斂快且計(jì)算快速，相較于梯形法，計(jì)算時(shí)間減少近一半；③仿真精度方面，與梯形法相近，但遠(yuǎn)小于歐拉法。

穩(wěn)定的計(jì)算時(shí)間有助于通信程序的設(shè)置與實(shí)時(shí)仿真模塊之間的數(shù)據(jù)交互，而快速計(jì)算可以減小步長，有助于降低多重開關(guān)事件的出現(xiàn)，降低仿真計(jì)算難度。

算例中使用100 MHz 晶振的FPGA 開發(fā)板，如使用RT-LAB 等實(shí)驗(yàn)級FPGA 開發(fā)板，仿真速度還可進(jìn)一步提升。后續(xù)將研究實(shí)驗(yàn)級FPGA 開發(fā)板對仿真速度的提升情況，進(jìn)一步探究步長減小與電力電子仿真開關(guān)頻率提升程度的具體關(guān)系。