基于GEP算法的高階常微分方程預測模型

崔 未，王衛華，黃樟燦，談 慶

武漢理工大學 理學院，武漢 430070

1 引言

隨時間變化的動態系統在日常生活中廣泛存在，如溫度變化、降水量變化、金融數據變化等。如何對隨時間變化的動態系統建模完成預測，一直是研究的熱點。合適的時間序列模型對投資風險把控、投入產出評估等方面具有重要意義。

時間序列預測是一種根據現有數據的規律信息構建模型，并將模型外推來完成預測的方法。預測效果主要受模型選擇影響[1]，這是由于時間序列數據是一個無顯著規律的動態系統，數據隨時間波動變化規律復雜，不同模型對數據的處理、構建聯系和規律發現的方法有很大區別，使得模型對歷史數據的描述有不同程度的偏差，進而對預測結果有直接影響。針對隨時間變化的動態系統采樣數據，通常采用ARIMA模型與人工神經網絡（ANN）模型進行預測。ARIMA模型是一類線性模型，對非線性問題處理能力較差。為提高它的非線性能力，文獻[2-5]將ARIMA模型與深度信念網絡、支持向量機和GARCH等模型結合，在紅潮預測、鈾礦價格預測、網絡流量預測和地鐵乘客量短期預測等方面取得了一定的效果。雖然上述方法使得ARIMA模型在許多非線性系統預測問題中得到了運用，但是ARIMA模型對數據平穩性要求高，對非線性數據規律的捕捉能力較弱的問題依舊存在。與ARIMA模型及其改進相比，ANN有更強的非線性能力和自適應能力，對復雜非線性數據的規律發現能力更強，但它存在對數據數量要求高、參數設定復雜和模型解釋性差的缺陷[6]。對此文獻[7-11]采用PSO、GA等優化算法對BP神經網絡和RBF神經網絡進行參數優化；文獻[12]將前饋神經網絡（FNN）嵌入到EMD預測框架中，利用EMD框架的加權重組策略對FNN進行賦權，實現了FNN參數的自適應；文獻[13]中將改進小世界網絡應用到Leaky ESN中，在實現稀疏連接的同時，改善了神經節點的連接方式，提升了Leaky ESN的適應性，取得了更高的預測精度。雖然這些混合算法對參數設定復雜、數據量要求高的問題進行了解決，在一類時間序列預測問題上取得了良好效果，但是它們本質上是一個無顯式模型的“黑盒”，對模型輸出缺少合理解釋，對模型的進一步分析研究存在阻礙。

為了彌補這些模型具有“黑盒”性質的缺陷，研究人員提出通過符號回歸，以代數形式構建微分方程模型。基于演化算法的常微分方程建模方法早在2000年就被提出[14]，該方法在建模過程中分為兩步，首先使用遺傳規劃（GP）構建模型，接下來通過遺傳算法（GA）對模型進行參數優化。由于GP的遺傳操作是基于樹結構進行的，因此在建模過程中，存在容易陷入局部最優的問題。除此之外，受制于兩步建模以及當時的計算條件，進行一次成功建模需要50 min，因此該方法雖然在預測精度上表現較好，但是未得到廣泛應用和推廣。

基因表達式編程（GEP）是一種融合了GA與GP優點的新型演化算法，尋優效率比二者高出100～60 000倍[15]。這是因為GEP對表達式樹結構的更改通過對線性編碼的修改完成，相比GP對樹結構直接進行修改有更高的靈活性，種群多樣性更豐富，能夠有效地在保證精度的前提下提升建模效率。

綜上，為了提升演化算法在常微分方程建模中的效率，改善人工神經網絡等非線性方法在時間序列預測建模中無法獲得顯式模型和精度較低的缺陷，本文提出基于GEP的高階常微分方程建模方法，并對太陽黑子年平均數進行建模預測。

2 基因表達式編程（GEP）

葡萄牙科學家C.Ferreira結合了遺傳算法（GA）和遺傳程序設計（GP）的優點，于2001年提出了基因表達式編程（GEP）[16]。GEP是一種擁有“基因型-表現型”特性的遺傳算法，即采用固定長度的線性編碼，對線性編碼使用遺傳算子進行操作，通過“Karva”語言將線性編碼轉化為表達式樹（ET），樹結構最終轉化為表達式。下面對GEP編碼及操作進行詳細介紹。

2.1 GEP編碼說明

GEP對染色體通過遺傳操作完成進化。染色體由基因通過連接函數構成。染色體是由函數集F和終止集T構成。函數集包含基本初等函數、邏輯運算符和變量；終止集T中有且只有變量及參數。基因分為頭部和尾部兩部分，頭部由函數集F與終止集T構成，尾部由終止集T構成。尾長度由頭長度決定，關系式為：

其中t為尾長度，n為運算符目數，h為頭長度。下面以一個例子對GEP的編碼進行說明。

形如“Q+*/aaaba//+//sin+abababa”的線性編碼對應的表達式為：

式（2）對應的表達式樹和染色體編碼如圖1所示。

圖1 式（2）對應的染色體編碼與表達式樹

從圖1可以看出，該表達式通過兩個等長基因構成的染色體表示，連接函數“+”將兩個基因連接。基因的頭結點長度為4，尾結點長度通過式（7）計算可知為5，總長度為9。基因1將尾部進行了部分表達，基因2對尾部沒有進行表達。這種編碼方式根據基因中符號集元素的分布情況，靈活地選擇是否對基因全部進行表達，與GP直接對樹結構操作相比更為靈活，與GA的二進制編碼相比有更豐富的表現形式。

2.2 GEP基本流程

GEP在種群在進化過程通過對染色體進行遺傳操作，使種群向符合適應度函數的方向進化[17]。GEP的操作流程如圖2所示。

對圖2進行總結，得到GEP的操作流程如下：

步驟1根據設定參數進行種群初始化。

步驟2計算種群適應度，判斷是否達到結束條件，是則跳出循環，否則進入步驟3。

步驟3進行遺傳操作。

步驟4得到一組新解，返回步驟2。

圖2 GEP操作流程圖

通過上述流程完成模型建立。下面對GEP的遺傳算子進行簡單介紹，遺傳算子包括選擇復制算子、轉座算子、重組算子和變異算子，四種算子及操作方法如下：

（1）選擇復制算子：GEP中采用基于精英策略的輪盤賭選擇，即對適應度最高的個體進行保留并復制，其余個體的選擇通過輪盤賭進行。

（2）轉座算子：轉座是指染色體中的基因片段能夠轉移到該染色體或其他染色體的各個位置。轉座根據轉移目標的不同，分為IS轉座、RIS轉座和基因轉座。

①IS轉座：選擇帶有運算符的基因片段，將該片段插入到基因的頭部中（除了起始位置以外的任意節點）。

②RIS轉座：從頭部開始選擇一個基因片段，從片段中選取第一個可以表示完整函數的部分作為RIS片段，將該片段插入到基因的尾部中。

③基因轉座：將整個基因插入到染色體的起始位置，染色體長度保持一致。

（3）重組算子：重組根據基因片段截取方式不同分為單點重組，兩點重組和基因重組。

①單點重組：父代染色體配對，確定截斷點，將截斷點后的部分進行交換。

②兩點重組：父代染色體配對，通過兩個階段點確定截取片段，對截取片段進行交換。

③基因重組：父代染色體配對，對選中基因進行交換。

（4）變異算子：變異算子是最高效的遺傳算子，變異可以發生在染色體的任何位置。在頭部中的運算符和終止點均可以通過變異相互轉換。尾部的終止點只能進行終止點變異。

染色體經過遺傳操作，不斷迭代，最終得到適應度最高的個體，即為最優解。下面對基于GEP的高階常微分方程建模進行詳細介紹。

3 基于GEP的高階常微分方程預測模型

3.1 數學方法說明

針對一組動態系統觀測樣本：

對這組數據進行預測的方法是通過不同方法對這組數據進行規律挖掘，建立模型 f如下：

通過 f進行n步外推，對后n個時刻的數據進行預測。常見的動態系統大多為非線性系統，如股票價格、溫度變化和太陽黑子變化等。它們都是非線性、非平穩、復雜度高的時間序列數據，對數據規律的發現和總結存在較大難度，因此以人工神經網絡為代表的非線性建模方法在這類數據的預測上得到了廣泛應用。但是這類方法的問題在于它們所建的模型 f是一個隱式模型，只能通過對參數的修改完成對模型的改進，無法對模型進行更為深入的研究。為了改善這個缺陷，本文利用常微分方程對動態系統準確描述的特性，構建高階常微分方程模型 f并完成預測。對X構建高階常微分方程的主要步驟如下所述：

對于規模為n的動態系統X進行m階差分近似，得到矩陣 Mn×m：

通過矩陣M中的元素，對X進行高階常微分方程建模得到：

使得：

成立。其中左邊部分具體為：

其中，n代表數據規模。對式（6）進行數值求解，得到接下來τ步的預測值：

式（6）中，f是由符號回歸而成的。基于GEP的高階常微分方程演化建模，即通過對規模為n的連續觀察數據進行差分處理，利用GEP的自組織、自學習的特點對已處理數據進行建模，得到形如式（6）的表達式。下面對上述步驟分別進行詳細說明。

3.2 有限差分

為了得到數據波動規律，首先需要對數據進行有限差分處理，計算導數近似值[18]，構建矩陣M。以一階導數為例，通過下面的數值差分公式得到導數近似值。

對t0時刻的數據進行如下處理：

對t1～tn-3時刻的數據進行如下處理：

對tn-2～tn-1時刻的數據進行如下處理：

其中，h代表單位時間間隔，xi代表i時刻的數據。更高階的計算以此類推，保證不同時刻的數據誤差始終為O(h2)。

3.3 GEP-HODE算法

通過有限差分得到矩陣式（5），使用GEP對M 構建如式（6）所示函數關系 f，整體流程如圖3所示。

圖3 GEP-HODE算法流程圖

根據圖3，對算法步驟進行歸納，具體流程如下：

輸入：時序數據X，時間t。

輸出：高階常微分方程模型 f。

步驟1輸入HODE模型階數m。

步驟2計算m階導數近似值，得到矩陣M。

步驟3確定迭代次數并對GEP進行參數設定。

步驟4進行種群初始化，并計算適應度。

步驟5對種群進行遺傳操作，形成新種群。

步驟6根據迭代次數判斷是否結束運算，是則結束，否則返回步驟2。

步驟7根據GEP建模得到模型 f。

步驟8利用數值方法對 f進行求解，得到最終預測值。

3.4 適應度函數選擇和模型求解

3.4.1 適應度函數選擇

模型的建立是通過GEP以適應度函數為目標演化進行的，因此適應度函數的選擇對模型的構建有直接影響。GEP-HODE模型是以式（7）為評判標準的，即向模型擬合值與實際值誤差最小的方向進化。標準絕對誤差函數的作用與式（7）相同，故本文選擇該函數作為適應度函數，如式（13）所示：

其中K為選擇范圍，一般情況下選擇1 000。CY為樣本數，Cj為模型所得，Tj為目標值。

3.4.2 模型求解

由GEP-HODE算法得到的高階常微分方程模型是一類非線性模型，因此常采用數值解法進行求解。這種方法的思路是將高階方程轉換為多個一階方程構成的方程組，利用數值方法對一階方程逐個求解，利用所得結果進行遞推，進而得到高階方程的數值解。針對式（6）所示模型，首先進行如下變換：

根據變量替換，將式（6）所示的高階常微分方程轉化為式（15）所示的常微分方程組：

其中t表示t0～tn的所有時間節點。常微分方程組的初值為：

通過數值方法，對式（15）所示的方程組從下至上遞歸求解，最終得到預測值y0( )t+τ，其中τ為預測步長。

4 計算機建模與仿真

4.1 數據選擇與參數設定

本文使用太陽黑子年平均數據[19]，設計實驗進行對比分析。實驗使用1749—1919年太陽黑子年平均數據，與文獻[14]中所用數據一致。其中訓練集為1749—1863年的數據，樣本總數為115，測試集為1864—1919年的數據，樣本總數為56。對1920—1924五年的太陽黑子年平均數進行預測，將GEP-HODE模型的預測結果與文獻[14]中方法（下面稱為GPGA方法）、SVM-ARIMA模型[20]、PSO-BP神經網絡模型[21]和單因素GEP預測模型進行結果對比分析。數值解法采取步長為0.1的四階龍格庫塔算法，初值為最后一個數據點對應的值。參數設定表如表1所示。

4.2 誤差評判標準

本文選取的誤差評判標準為平均絕對誤差MAE、平均絕對誤差百分比MAPE和均方誤差RMSE。平均絕對誤差用來衡量預測值與實際值的直接差距；平均絕對誤差百分比用來衡量預測值相對實際值的變化情況；均方誤差用于判斷預測值中是否有偏離實際值較大的點。它們的計算公式分別為：

表1 參數設定表

4.3 實驗結果與對比分析

1749—1919年的太陽黑子年平均數量樣本共171個，樣本按年分布如圖4所示。

圖4 1749—1919年太陽黑子年平均數目分布圖

從圖4中可以看出太陽黑子年平均數據是一組非線性、非平穩的時間序列數據，沒有明顯的變化規律。采取GEP-HODE方法進行常微分方程建模，得到1～4階模型。每階模型進行五次建模，取預測誤差平均值進行分析，各階模型預測誤差對比如表2所示。

表2 各階模型預測誤差對比表

從表2可以看出在三種誤差評判標準下，三階常微分方程模型（ODE（3））在對太陽黑子年平均數據建模預測效果最好。效果最好的ODE（3）模型如下所示：

模型參數如下：

為了體現GEP在建模中的高效性和模型的精確性，下面將式（20）所示模型（下面稱為GEP-HODE模型）與GPGA等四種方法進行對比。預測精度結果如表3所示。

表3 預測結果誤差對比表

從表3中能夠看出GEP-HODE模型的預測誤差在三種評判標準上均小于后三種預測模型，在MAE與RMSE上低于GPGA方法，在MAPE上高于GPGA，總體預測結果較為接近，下面對造成這種結果的原因進行分析。太陽黑子年平均數是一個隨時間劇烈變化的非線性動態系統，ARIMA本質上是一個線性模型，即使通過SVM將數據映射到高維空間試圖削弱非線性造成的影響，也無法獲得更為準確的預測結果，同時SVM的參數設定以及徑向基函數的選擇也較為主觀，因此SVMARIMA模型預測結果在五種方法中最差；單因素GEP預測模型的問題在于其本質上是一高階復合函數擬合模型，由于數據維度為1，因此該模型對數據的規律總結能力與微分方程模型相比較弱，但與ARIMA相比，它的非線性能力更強，因此該模型的預測結果相較SVMARIMA模型更好；BP神經網絡盡管通過PSO算法對各連接層的權值進行了優化，但數據量較少、隱藏層數人為設定等問題使得模型訓練效果較差，導致預測精度較低，但是該模型出色的非線性能力使得它在常微分方程模型之外的三種對比方法中精度最高；GEP-HODE與GAGP預測誤差相近是因為二者同屬于對常微分方程的演化建模，因此它們本質上是一致的。由于使用的演化算法不同，導致了它們的結果雖然相近，但是仍然存在著一定的差異。因此需要對造成這種差異的原因進行進一步分析。

首先對兩種方法的預測結果進行對比，如圖5所示。

圖5 兩種方法預測結果對比圖

從圖5可以看出GPGA方法在整體數據點的數值預測上精度要略高于GEP-HODE方法，從表3中則可以發現二者在預測誤差上相差不大。作為一種預測手段，對數據趨勢的預測也屬于模型預測能力的一部分，從圖5中可以看出GPGA方法在1923—1924年這一段的預測與實際趨勢完全相反，而GEP-HODE方法則準確地反映了這段時間的實際情況。因此，雖然二者在精度上是比較接近的，但從預測圖像上來看，GEP-HODE方法相較GPGA方法表現更好。從理論上看，二者均屬于遺傳算法，在模型表達方面和遺傳操作方面十分類似，因此得出的數值解精度也十分相近。導致預測趨勢產生偏差的原因在于GPGA方法通過GP進行模型建立，利用GA進行參數優化，這種方法受制于GP的編碼結構，導致模型結構存在一定的缺陷。詳細的說，由于GP在進化過程中通過對樹結構進行變換和修飾，實現遺傳操作獲得新解，這種操作方式會導致解的多樣性不足，因此GP容易陷入局部最優。GEP算法改進了編碼方式，將遺傳操作應用在線性編碼中，通過表達式樹對遺傳操作進行體現。這種“基因型-表現型”的操作方式豐富了解的多樣性，相較GPGA方法跳出局部最優的能力更強，因此在模型的建立方面相較GPGA存在一定優勢，使得預測趨勢更為精確。

接下來對兩種方法的建模效率進行對比。兩種方法在函數集相同、連接函數相同、迭代次數相同以及同一環境下（intel core i5 4200h）重復進行五次，將兩者用時進行對比。對比結果如表4所示。

表4 建模用時對比表

從表4中可以看出五次實驗GEP算法的效率遠高于GPGA方法。一方面因為GEP通過改進編碼方式，從本質上比GP的進化效率高；另一方面由于GPGA屬于混合算法，在建模過程中首先需要通過GP進行模型建立，隨后需要通過GA進行參數優化，“兩步走”的建模過程導致了運算時間的增加。與之相比，GEP一方面將模型建立與參數優化同時進行，另一方面采取了更為高效的編碼方式，因此用時更短，效率更高。

綜上，GEP-HODE方法在具有代表性的太陽黑子年平均數據預測上取得了一定的成果，與改進人工神經網絡模型和混合ARIMA模型相比具有更高的精度，能夠在保證預測精度的同時能夠獲得顯式模型的優勢；與GPGA混合建模方法相比擁有更高的效率，具備更優秀的趨勢預測能力。如何通過減小導數近似帶來的誤差以及對GEP算法進行改進，提高模型對歷史數據的描述能力，從而進一步提升模型精度是接下來要研究的問題。

5 結語

本文針對動態系統預測建模問題中存在的建模效率低和無顯式模型的缺陷，提出了一種基于GEP的高階常微分方程建模方法。該方法兼容了高階微分能夠對短期波動數據變化具有較強描述能力，能夠深入挖掘數據本征的特點和GEP算法的自組織、自適應特性。具備高效率、高精度和擁有顯式模型的優點。以不同時間段的太陽黑子年平均數為例構建模型，通過與GAGP、SVM-ARIMA、PSO-BP等方法進行對比，分析說明GEPHODE方法的優勢以及造成其他方法出現誤差的原因，最終結果表明GEP-HODE方法在精度和效率上分別具有一定的優勢，可以對動態時間系統進行建模并完成預測。