張多利,魏 可,胡永陽,聶言碩,侯 寧,宋宇鯤

(1.合肥工業(yè)大學 電子科學與應(yīng)用物理學院,安徽 合肥 230601; 2.教育部IC設(shè)計網(wǎng)上合作研究中心,安徽 合肥 230601; 3.河南城建學院 電氣與控制工程學院,河南 平頂山 467041)

常微分方程(ordinary differential equation,ODE)是對世界認知的一種基礎(chǔ)數(shù)學方法。常微分方程的初值問題(initial-value problem)在滿足李普希茨(Lipschitz)條件時存在唯一的連續(xù)可微解[1],因為其精確解析解難以獲得,所以通常使用數(shù)值方法計算求解。在實時檢測宏觀和微觀事物、自然現(xiàn)象、運動規(guī)律和連續(xù)變化過程的工程應(yīng)用中,為高效求解需求量大、運算量大、精度要求高的ODEs,往往需要采用硬件代替通用處理器實現(xiàn)。

四階龍格-庫塔法(fourth-order Runge-Kutta,RK4)是一類高精度迭代法,常用于非剛性O(shè)DE數(shù)值解。目前,基于現(xiàn)場可編程邏輯門陣列(field programmable gatearray, FPGA)的RK4微分方程數(shù)值解算實現(xiàn)已得到廣泛應(yīng)用[2-3]。文獻[4-5]介紹了實時求解洛倫茲系統(tǒng)的FPGA實現(xiàn),分別應(yīng)用于數(shù)據(jù)流加密和混沌動力系統(tǒng)計算中;文獻[6-7]采用RK4結(jié)構(gòu)在FPGA上分別實現(xiàn)了風力渦輪機和感應(yīng)電機的動態(tài)仿真;文獻[8]介紹了系數(shù)可配置的衛(wèi)星飛輪RK4動力模型,并在FPGA完成了硬件實現(xiàn);文獻[9-13]在不同的微弱信號檢測應(yīng)用中實現(xiàn)了Duffing微分方程RK4的FPGA實現(xiàn);文獻[14]完成了使用RK4求解彈道方程的FPGA實現(xiàn);文獻[15]介紹了在FPGA中采用RK4求解醫(yī)學領(lǐng)域PK/PD模型。

上述研究均針對確定ODE展開,其求解電路結(jié)構(gòu)的通用性考慮不足。其中,對于RK4計算結(jié)構(gòu),文獻[4]采用計算全展開結(jié)構(gòu)以獲得高性能;文獻[5-7]采用迭代復(fù)用結(jié)構(gòu),減少了所需電路面積。對于ODE計算結(jié)構(gòu),文獻[9,12-13]采用非流水線結(jié)構(gòu),在同一時間內(nèi)只能完成單個ODE求解;文獻[10]采用流水線結(jié)構(gòu),大批量數(shù)據(jù)處理速度優(yōu)勢。對于ODE求解靈活性方面,文獻[12-13]通過多次求解得到不同系數(shù)條件下的求解結(jié)果;文獻[8]給出的求解方程可改變相關(guān)系數(shù)數(shù)值,能夠完成不同系數(shù)條件下的方程求解。

綜上所述,為了滿足各種應(yīng)用對ODEs求解的靈活性、高效性、實時性等多方面的需求,本文設(shè)計了一種高性能通用可重構(gòu)微分方程解算器(reconfigurable differential equation solver,RDES),RDES的主要特征如下:

(1) 采用可重構(gòu)的原理,支持在線重構(gòu)完成不同結(jié)構(gòu)、階數(shù)、變量等條件的ODEs求解任務(wù),增加系統(tǒng)通用性。

(2) 可根據(jù)需要重構(gòu)為流水線或循環(huán)迭代結(jié)構(gòu)完成RK4計算流程,有效滿足不同求解收斂速度的要求。

(3) 確定ODEs下,支持不同初始條件下的批量求解,顯著提高初始值、系數(shù)等因素的收斂速度。

1 系統(tǒng)架構(gòu)設(shè)計

1.1 RK4算法原理

RK4是由歐拉方法改進得到的一種在工程上應(yīng)用廣泛的高精度ODE數(shù)值解法,可求解一階微分方程或微分方程組。對于包含多個變量的ODE,對應(yīng)的RK4求解公式為:

(1)

(2)

其中,h為步長。對于變量x,4組Kx值對應(yīng)的表達式為:

(3)

y、z等變量表達式與(3)式相似。

1.2 系統(tǒng)架構(gòu)

RDES架構(gòu)和系統(tǒng)工作流程如圖1所示。RDES基于可重構(gòu)方法實現(xiàn),整體架構(gòu)包括RK模塊、可重構(gòu)計算陣列(reconfigurable computing array,RCA)模塊、迭代循環(huán)控制(iterative loop controller,ILC)模塊、配置信息和變量數(shù)據(jù)存儲器(memory)及與外界交互信息的接口(interface)。

由圖1b可知,通過編譯工具整理待求解ODE的結(jié)構(gòu),合并分支,確定其包含的基礎(chǔ)算子。將整體計算結(jié)構(gòu)拆解為由基本運算單元組成的單步或連續(xù)計算?；具\算單元為RCA包含的加、減、比較和選擇、乘、除、開方、三角函數(shù)計算等。通過映射工具將拆解出的各組運算映射到RCA中,并生成相應(yīng)的配置信息。完成后RDES就能開始工作。

圖1 RDES架構(gòu)和系統(tǒng)工作流程

在RDES的工作流程中,通過interface寫入配置信息和變量數(shù)據(jù)信息,存儲到memory中,變量數(shù)據(jù)可以是大量不同初值條件ODEs的變量信息。首先取出配置信息對RCA進行配置,配置完成后開始ODE求解計算。對應(yīng)(2)式計算流程,從memory中逐個取出各ODE的變量數(shù)據(jù)x0,y0,…,z0,輸入RK模塊,經(jīng)RCA計算得到新一組變量x1,y1,…,z1。不滿足停止計算條件的變量數(shù)據(jù)寫回memory,作為下一組變量數(shù)據(jù)輸入RK模塊計算。由此迭代循環(huán)直到輸入系統(tǒng)的所有ODEs都滿足條件,批量求解完成。系統(tǒng)的數(shù)據(jù)交互、迭代循環(huán)計算流程都由ILC控制實現(xiàn)。

1.2.1 RK模塊

RK模塊結(jié)構(gòu)如圖2所示,該模塊是按Runge-Kutta計算流程設(shè)計的系統(tǒng)核心計算模塊,主要由RCA模塊、f(Ki)運算模塊、斜率生成(slope generator)模塊組成。RCA模塊根據(jù)輸入的變量值計算各個變量的Ki值;f(Ki)運算模塊根據(jù)需要對已獲得的K值進行處理;斜率生成根據(jù)所有RCA輸出的K值計算變化斜率,進而得到輸出結(jié)果。

圖2 RK模塊結(jié)構(gòu)

RK模塊可實現(xiàn)任意階數(shù)的Runge-Kutta計算。本文已完成的RK模塊取i值為4,適用于四階以內(nèi)的Runge-Kutta計算。以RK4為例,計算步驟由(2)式、(3)式得到,即輸入變量初值x0,y0,…,z0逐次經(jīng)過4個RCA模塊得到K1～K4的值,3個f(Ki)運算模塊分別計算K1h/2、K2h/2、K3h的值,斜率生成計算(K1+2K2+2K3+K4)h/6,最后與初值求和,得到輸出新值x1,y1,…,z1。

1.2.2 RCA模塊

RCA是RDES的基本計算模塊,是由每行4個處理單元(processor elements,PE)形成的固定互連結(jié)構(gòu)與行間數(shù)據(jù)傳輸網(wǎng)(data transmission network,DTN)交替排布組成的陣列,模塊結(jié)構(gòu)如圖3所示。通過重構(gòu)實現(xiàn)不同ODEs的計算結(jié)構(gòu),即(1)式所有輸入變量到輸出變量一階導(dǎo)數(shù)值的數(shù)據(jù)路徑。

圖3 RCA模塊結(jié)構(gòu)

本文已完成的RCA結(jié)構(gòu)行數(shù)為5行,適用于包含4個以內(nèi)變量數(shù)的ODEs。能夠?qū)崿F(xiàn)同行相鄰PE的數(shù)據(jù)交互,行間數(shù)據(jù)經(jīng)DTN的單向傳遞,能夠處理ODE中的計算分支、變化系數(shù)、變量復(fù)制等情況。對于超過4變量或超過單個RCA結(jié)構(gòu)允許計算量的ODEs,可通過多次重構(gòu)計算。

RCA模塊中3種PE為固定互連關(guān)系的不同運算實現(xiàn):① 比較選擇計算單元,如圖3中PE1,包含浮點比較器和選擇器,可以處理ODE計算中的分支情況;② 乘加計算單元,如圖3中PE2,包含浮點乘法器和加減法器,實現(xiàn)最常用的乘法和加減法計算;③ 復(fù)雜計算單元,如圖3中PE3包含浮點除法器、開方器、CORDIC計算單元,主要實現(xiàn)三角函數(shù)和多種超越函數(shù)的計算。

RCA模塊主要由控制路徑和數(shù)據(jù)路徑2個部分組成。外層RK模塊通過控制路徑寫入配置信息,廣播到所有PE和DTN進行配置。配置完成后,RCA模塊形成從變量輸入到K值輸出的完整數(shù)據(jù)路徑,即可根據(jù)輸入變量值計算得到對應(yīng)K值,輸出至RK模塊。

1.2.3 ILC模塊

ILC模塊負責RDES系統(tǒng)從信息配置、數(shù)據(jù)輸入到迭代循環(huán)停止的完整流程的控制。ILC模塊控制配置信息和變量數(shù)據(jù)的接收。檢測到配置完成信號后,控制取出初始變量數(shù)據(jù)x0,y0,…,z0,寫入RK模塊計算。RK模塊計算輸出x1,y1,…,z1。ILC模塊在數(shù)據(jù)計算輸出后,根據(jù)停止迭代信號iteration-finish的值判斷下一步操作。若信號值為0,則保留該組數(shù)據(jù),作為新的輸入變量數(shù)據(jù)作下一輪迭代計算;若信號值為1,則ODE結(jié)束計算,該組數(shù)據(jù)不再存儲。ILC模塊對結(jié)束計算的ODE計數(shù),finish-num計數(shù)到輸入的N個方程都結(jié)束則完成計算任務(wù)。ILC模塊控制所有變量遵循同樣的迭代循環(huán)流程,其中對變量x的迭代循環(huán)控制過程的偽代碼如下:

while finish-num

if configuration done then

if iteration-time==0 then

x←x0;

data valid sign is decided by input data

valid sign.

else

x←x1;

if iteration is not done then

data valid sign is decided by iteration data valid sign.

else

data valid sign is set to 0;

finish-num←finish-num+1;

end while

通過配置信息可以得到停止迭代條件,設(shè)計支持包括迭代次數(shù)限制、停止時間限制、變量數(shù)值達標等多種條件,滿足實際應(yīng)用中不同ODEs的需求。ILC模塊在每輪迭代后作相應(yīng)的條件判斷,結(jié)果寫入iteration-finish信號,以完成對變量數(shù)據(jù)的控制。

2 應(yīng)用映射與實驗分析

軟件實現(xiàn)在Intel(R) Core(TM) i7-4790 CPU@ 3.60 GHz,32 GiB內(nèi)存、64位操作系統(tǒng)的PC機上,使用Microsoft Visual Studio 2017軟件的x86編譯環(huán)境,用C++語言作為編程工具,設(shè)定數(shù)據(jù)格式為雙精度浮點數(shù);硬件實現(xiàn)采用Xilinx Vivado 2018.3作為編譯工具,并在Virtex-7 XC7V2000TFLG上實現(xiàn)整個設(shè)計,設(shè)定數(shù)據(jù)格式為單精度浮點數(shù)。使用Matlab軟件完成計算數(shù)據(jù)圖像繪制,并對比軟硬件執(zhí)行時間。

2.1 范德波(Van der Pol)方程

Van der Pol方程是描述振蕩的一種基礎(chǔ)模型。其往復(fù)的周期運動易受擾動影響,且對初值變化敏感。因此實際研究時往往需要實時求解。方程常用形式為:

x″+μ(x2-1)x′+x=0

(4)

其中,μ為可變的阻尼系數(shù)。

將方程轉(zhuǎn)化為一階微分方程組。設(shè)y=x′,則有:

(5)

2.1.1 方程映射

該方程計算數(shù)據(jù)流圖(data flow graph,DFG)及其在圖3 RCA中的映射結(jié)構(gòu)如圖4所示。

DFG通過拆解分組(set),映射到RCA的PE中,形成圖4的計算路徑,并作為配置信息存儲下來。這里1個Van der Pol方程使用了2個PE2完成映射。

2.1.2 實驗與結(jié)果分析

實驗100組不同初值或系數(shù)條件的Van der Pol方程求解。分別在PC和FPGA板上完成計算。

硬件實現(xiàn)結(jié)果能夠滿足Van der Pol方程實時求解的需求,其中使用Matlab軟件作出的隨機選取的幾組FPGA實測結(jié)果和軟件結(jié)果的對比圖像如圖5所示。軟硬件計算結(jié)果最大絕對誤差為2.043×10-4,因此軟硬件計算結(jié)果可視為一致,不同初始條件下的方程可能求解出不同的振蕩關(guān)系,據(jù)此可使初始值、系數(shù)條件的選擇范圍快速收斂。Van der Pol方程實驗軟硬件實現(xiàn)時間對比見表1所列。由表1可知,硬件實現(xiàn)相對于軟件實現(xiàn)約有16.93倍的加速比。

圖5 幾組Van der Pol方程軟硬件對比圖像

表1 軟硬件計算Van der Pol方程時間比較

2.2 阿波羅(Apollo)衛(wèi)星運動軌跡方程

人造衛(wèi)星的運動軌跡與三體引力、太陽光壓、當前位置和速度、調(diào)整軌道的推力等多方面因素有關(guān)。因其復(fù)雜多變的干擾因素,對批量求解衛(wèi)星在不同條件下的運動軌跡就有了需求。采用簡化的Apollo衛(wèi)星運動軌跡方程作為測試例,求解不同初值條件下Apollo衛(wèi)星運動軌跡和狀態(tài)。衛(wèi)星運動軌跡(x,y)滿足:

(6)

其中:μ為變化系數(shù);λ=1-μ。

將方程轉(zhuǎn)化為一階微分方程組。設(shè)x1=x′,y1=y′,則有:

(7)

實驗150組不同初值或系數(shù)條件的Apollo運動方程求解。硬件實現(xiàn)結(jié)果能夠滿足Apollo方程批量求解的需求,其中隨機選取的幾組FPGA實測結(jié)果和軟件結(jié)果的對比軌跡圖像如圖6所示。

圖6 幾組Apollo運動軌跡軟硬件對比圖像

軟硬件計算結(jié)構(gòu)最大絕對誤差為6.524×10-4,因此軟硬件計算結(jié)果可視為一致,不同初始條件下的方程求解可能產(chǎn)生不同的運動軌跡,據(jù)此可對滿足需求的初始條件快速篩選,從而實現(xiàn)快速收斂。

Apollo運動方程實驗軟硬件實現(xiàn)時間對比見表2所列,由表2可知，硬件實現(xiàn)相對于軟件實現(xiàn)約有112.64倍的加速比。

表2 軟硬件計算Apollo方程時間比較

2.3 炮彈彈道方程

炮彈的飛行軌跡受大氣狀況、初始條件、炮彈自身性質(zhì)等多方面因素影響。為在當前大氣條件下精準擊中目標點,需要大批量求解炮彈發(fā)射初始條件、變化系數(shù)條件有細微調(diào)整的彈道方程，因此批量快速求解彈道方程就有重要的意義。簡化彈道方程模型為包含速度v、飛行方向與水平面的夾角θ、當前位置距離發(fā)射點的垂直距離y、當前位置距離發(fā)射點的水平距離x4個變量,微分方程組表達式[16-17]為:

v′=-cH(y)F(Vτ)-gsinθ,

y′=vsinθ,

x′=vcosθ

(8)

其中:c為彈形系數(shù);g為重力加速度,是變化系數(shù);H(y)為空氣密度函數(shù);F(Vτ)為阻力系數(shù)，是中間變量[18]。

2.3.1 方程映射

在彈道方程計算流程中,變量v的K值求解流程較復(fù)雜,需要求解出6個中間變量,包括虛擬溫度τ、空氣密度函數(shù)H(y)、虛擬速度vτ、阻力系數(shù)F(Vτ)。將方程中求解τ和F(Vτ)的條件分支結(jié)構(gòu)通過比較和選擇操作合并計算后,得到彈道方程的部分DFG及對應(yīng)映射圖像如圖7所示。

圖7 彈道方程計算部分DFG和映射圖

2.3.2 實驗與結(jié)果分析

實驗200組不同初值或系數(shù)條件的炮彈彈道方程求解。硬件實現(xiàn)結(jié)果能夠滿足彈道方程批量求解的需求,其中隨機選取的幾組FPGA實測結(jié)果和軟件結(jié)果的對比彈道圖像如圖8所示。

圖8 幾組彈道軟硬件對比圖像

軟硬件計算結(jié)構(gòu)最大絕對誤差為1.482×10-6,因此軟硬件計算結(jié)果可視為一致,不同初始條件下的炮彈能夠產(chǎn)生不同的飛行軌跡,產(chǎn)生不同的落點,據(jù)此可對靠近目標落點的初始條件快速篩選,從而實現(xiàn)快速收斂。炮彈彈道方程實驗軟硬件實現(xiàn)時間對比見表3所列,硬件實現(xiàn)相對于軟件實現(xiàn)約有45.23倍的加速比。

表3 軟硬件計算彈道方程時間比較

3 結(jié) 論

本文提出了一種可重構(gòu)的微分方程通用解算器RDES,能夠重構(gòu)完成不同結(jié)構(gòu)、階數(shù)、變量等條件的ODEs求解任務(wù),拓展了系統(tǒng)通用性。幾組實際應(yīng)用的微分方程驗證實驗結(jié)果表明,RDES在完成不同初始條件的ODEs批量求解計算時較通用處理器具有很大的速度優(yōu)勢。