基于FPGA 的高頻電力電子裝置快速實時仿真數(shù)值算法

施 文，遲 頌，郭 亮

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

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

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

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

1 高頻電力電子實時仿真存在問題

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

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

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

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

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

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

2.1 改進的梯形迭代算法

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

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

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

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

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

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

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

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

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

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

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

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

圖3 Buck 電路拓撲Fig. 3 Topology of buck circuit

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

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

仿真精度方面， 改進梯形法與梯形法相近，都遠小于歐拉法；在計算時間方面，相較于梯形法，改進梯形法時間約減少了45.5%。

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

2.2 變步長算法計算步驟

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

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

步驟1從tn-1到tn時刻進行計算， 并檢測到開關(guān)動作；

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

表3 不同算法的計算時間Tab. 3 Computation time of different algorithms

4 實時仿真平臺實驗驗證

4.1 仿真程序?qū)嶒炂脚_

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

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

4.2 實驗用仿真模型

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

圖7 三相逆變器電路拓撲Fig. 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 倍基頻整定。

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

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

4.3.1 算法收斂性驗證

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

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

4.3.2 算法仿真驗證

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

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

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

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

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

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

5 結(jié)論

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

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

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