可重構(gòu)的微分方程通用解算器研究和實(shí)現(xiàn)

張多利,魏 可,胡永陽(yáng),聶言碩,侯 寧,宋宇鯤

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

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

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

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

綜上所述,為了滿足各種應(yīng)用對(duì)ODEs求解的靈活性、高效性、實(shí)時(shí)性等多方面的需求,本文設(shè)計(jì)了一種高性能通用可重構(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計(jì)算流程,有效滿足不同求解收斂速度的要求。

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

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

1.1 RK4算法原理

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

(1)

(2)

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

(3)

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

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

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

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

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

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

1.2.1 RK模塊

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

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

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

1.2.2 RCA模塊

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

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

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

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

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

1.2.3 ILC模塊

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

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

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

2 應(yīng)用映射與實(shí)驗(yàn)分析

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

2.1 范德波(Van der Pol)方程

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

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

(4)

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

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

(5)

2.1.1 方程映射

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

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

2.1.2 實(shí)驗(yàn)與結(jié)果分析

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

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

圖5 幾組Van der Pol方程軟硬件對(duì)比圖像

表1 軟硬件計(jì)算Van der Pol方程時(shí)間比較

2.2 阿波羅(Apollo)衛(wèi)星運(yùn)動(dòng)軌跡方程

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

(6)

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

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

(7)

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

圖6 幾組Apollo運(yùn)動(dòng)軌跡軟硬件對(duì)比圖像

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

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

表2 軟硬件計(jì)算Apollo方程時(shí)間比較

2.3 炮彈彈道方程

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

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

圖7 彈道方程計(jì)算部分DFG和映射圖

2.3.2 實(shí)驗(yàn)與結(jié)果分析

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

圖8 幾組彈道軟硬件對(duì)比圖像

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

表3 軟硬件計(jì)算彈道方程時(shí)間比較

3 結(jié) 論

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