基于深度學習長時預測的非線性系統建模方法研究

李聰，時宏偉

（四川大學計算機學院，成都610065）

0 引言

非線性動力系統辨識在現代控制領域內一直是致力于解決的難題之一。它基于這樣的假設，系統能夠由一組確定的常微分方程所確定，然而由于實際的建模過程中，這樣的方程往往難以建?；蛘呓⒌哪Ｐ头浅碗s，于是就想通過機器學習的方法來解決該問題。文獻[1]通過從數據中建模出動態化線性模型，對比自適應控制模型，該模型性能優異。文獻[2]利用電壓數據預測電路故障。文獻[3]使用高斯過程回歸建立預測模型，文獻[4]從數據當中建立導彈退化故障模型，從而預測導彈的故障。這些模型都是從歷史的數據中建模出系統動力模型，然后用該模型預測下一時刻的系統變化。即從系統采集到的數據發掘系統的動態變化規律。對于一般的動力系統，我們往往能獲取系統中各個時刻的狀態參數，想通過這些數據來揭示出系統模型。現在常用的方法有基于動態模式分解（Dynamic mode Decomposition，DMD）[4]，試圖在局部時間關系建立線性關系，通過奇異值分解的方法求解偽逆的方法來求解動力系統中的參數。也有通過加入先驗知識通過搜索的方法揭示出動力學模型，事先在一個備用的高次多項式的字典中搜索[5]，直到找到滿足條件的動力學系統，但是文中所示實驗僅僅以Lorenz 三維系統，方法新奇但不適合一般的較復雜的動力系統。也有基于神經網絡結合Jacobian 正則的方法來建模動力系統[7]，能夠很好的長時建模和預測動力系統，但是模型本身沒有建模長時動力系統，隨著迭代增長動力系統誤差必然增大。文獻[8]建模了多個時間步動力系統，但是該方法在Lorenz 系統上表現不佳，隨著迭代時間步的增長，一點點誤差導致系統狀態參數誤差增長不可控制。基于深度Koopman 算子建模非線性動力系統，然而系統僅建模了單步動力系統[11]，所以本文考慮在單步模型上考慮到了減少模型誤差，所以引入了總均方誤差（Total Least Square，TLS）[12]的方法求解單步動力模型。文獻[15]首次使用TLS 在非線性動力系統辨識中，文中模型首次將輸入數據和輸出數據都考慮為含噪數據，且噪聲性質相同，即輸入數據和輸出數據是對稱的，此外在模型上擴展了均方誤差，文中使用廣義線性模型來建模數據。由于擴展了模型，使得廣義線性模型具有更一般的表達能力。

非線性動力系統的數據也可以看作是時間序列數據，時間序列數據往往具有周期性，但往往不具有標準的周期，在一定的時間內具有某種非規范的周期性。為解決這個問題，文獻[13]首次提出將時間序列數據分解為多個內蘊函數（Intrinsic Mode Function，IMF），非線性動力系統所表示的數據也可以表征為某種形式的內蘊函數的加權疊加，也有一些在這方面的嘗試[20]。還有一些模型考慮從過去的多個時間步中提取信息，建模非線性動力學，借助長短時記憶網絡（Long Short Memory Network，LSTM）[23]的記憶能力來表達系統的非線性特征。例如線性循環自編碼網絡[24]考慮了前面多個時刻的狀態參數，將多個時刻的數據作為系統輸入，再利用自編碼機將數據映射到某一線性空間，相鄰時間步的系統狀態數據也建模為線性關系，最后在幾個數據集上驗證，獲得不錯的結果。

與之類似的方法在文獻[25]提出，文章試圖用自編碼機從數據中學習出Koopman 不變子空間，文中方法在編碼后的數據中建模時序依賴關系，這仍能建模出數據的局部線性關系?；谏疃葘W習和Koopman 算子的結合的方法還有很多[26-27]，大致和文獻[14]類似，都是試圖學習出不變子空間，在線性空間上做狀態參數演化，這種方法均是僅考慮狀態參數在局部空間上的演化，沒有考慮模型隨時間演化誤差較大的問題。

為了解決上述問題，本文在深度動態模式分解的算法框架之下，引入了Jacobian 正則化方法增加神經網絡的魯棒性；其次在線性模型的部分廣義化了原模型，引入了Total Least Square 作為求解廣義化模型的方法；最后在幾個常用的測試數據集上進行實驗驗證，得到非常不錯的效果。算法模型

本文使用模型包含三部分。第一，本文中的基礎模型采用深度動態模式分解[14]，該模型在單步預測中取得較好效果，但在長期（時間步在1000～4000）預測上表現不佳，因此需要引入額外的約束。第二，在基礎模型之上，通過在損失函數中加入Jacobian 正則，使得模型能夠在保持單步預測不失真的條件下，又能讓模型在循環預測中抑制誤差的傳播，相比較于基礎模型，加入正則后的模型能夠在理論上保證誤差減小，得在后續的預測中更加準確。第三，總體最小二乘法（Total Least Square，TLS）方法減小基礎模型中線性模型的誤差，TLS 也是減小誤差的方法，只不過在輸入上也考慮到誤差，且TLS 是廣義的最小二乘法，如果輸入數據中沒有誤差，TLS 會退化到普通最小二乘法（Ordinary Least Square，OLS）。

1 建模方法

1.1 非線性動力模型

本文中考慮的非線性動力系統的形式為：

其中系統狀態向量xt∈Rn為系統t 時刻系統狀態參數，xt+1∈Rn為下一時刻的系統狀態參數，函數f ∈Rn表示當前時刻系統參數到下一時刻系統狀態參數的映射，P 為其參數。在模型訓練完成之后，可以使用單步預測。對于長時預測，使用模型需要循環預測，初始化系統參數x0已知，公式化表達如下：

1.2 深度動態模式分解

文獻[14]基于擴展動態系統（Extend Dynamic Mode Decomposition）提出深度動態模式分解（Deep Dynamic Mode Decomposition，D-DMD），首次使用Autoencoder將系統狀態向量映射到線性空間，模型假設經過編碼后的數據，在時刻t 與t+1 時刻的zt，zt+1之間是線性關系。系統示意圖如圖1 所示。

圖1 深度動態模式分解模型框架

需要注意的是在模型僅考慮單個時刻間的動力學性質時。在實際使用模型時，t 時刻的系統參數參數xt總是被假設為已知，且被用來預測下一時刻的系統參數xt+1，這樣模型本身只具備有單步預測的能力，然而在實際應用中往往需要系統具備更長時間的預測能力。在文獻[14]中，上述可以被公式化描述為如下：

其 中 xt,xt+1∈Rn，zt,zt+1∈Rl，其 中 系 統 參 數P={θenc,K,θdec} ，n 是系統原始參數，l 是隱空間參數的維度。 fenc，fdec自編碼機的編解碼器，參數為Θ={θenc,θdec}，可以是全連接神經網絡，對圖像數據也可以是卷積神經網絡。記優化目標函數為J( K,θ )，上述算法的優化目標如下：

其中超參數λ1，λ2通常選取較小的值，例如0.01，0.001，||*||2表示l2范數，其中神經網絡參數Θ 通過梯度下降法訓練可以得到。矩陣參數K 通過求解偽逆得到[15]。為了敘述的連貫性，下面給出求解過程。對于經過神經網絡編碼的數據：Z=[z0,z1,z2,…,zT]，矩陣Zl×T表示隱空間中的數據集合。為了滿足線性關系，不妨記: X=[z0,z1,z2,…,zT-1] ，矩 陣 形 式 如（5）所 示。Y=[ z1,z1,z2,…,zT]，如（6）所示，矩陣X，Y ∈Rl×T，參數矩陣K ∈Rl×l。

線性表達式：

則：

其中X?表示矩陣X 的Moore-Penrose 偽逆運算[12，15]。根據（7）式要求解矩陣K，則需要求解矩陣X的偽逆，一般使用奇異值分解[5，15，28，29]的方法來求解矩陣K，首先對矩陣Y 進行奇異值分解，得到X=UΣV*，其中U 為l×l 正交陣，V 為T×T 的正交陣，Σ 為l×T 的奇異矩陣。則有：

通常而言在DMD 算法中需要多次左乘矩陣K，所以文獻[5，15，28，29]使特征值特征向量的方法計算。求解線性矩陣的過程中隱含了求解多變量最小均方誤差模型[5，15，28，29]：

即優化目標：

整個實際優化過程中，神經網絡參數的反向傳播和矩陣求解（11）式交替進行。

雖然模型（3）能夠很好地表達單步動力系統的關系，但是由于深度學習往往只能找到局部最優解，模型總存在一個較小的誤差，導致模型在實際使用時每一步的預測都存在著較小的誤差。假設上一步t=i 時刻預測的值為，t 時刻真實值xi，則t 時刻誤差為。如果采用長時預測，則下一步t+1 時刻的模型計算輸出為：

由于輸出數據中誤差εi的存在，導致模型輸出的不可控制，隨著迭代次數的增加，后續模型預測的值與真實值誤差越來越大。是要想表達長期的關系，需要在迭代過程中控制單步的誤差εi，以減小誤差對后續迭代的影響，讓單步誤差在傳遞過程中被減小，因此引入Jacobian 正則化[7]的方法來減小在單步迭代中的誤差。

1.3 Jacobian正則化方法

正則化方法在多個領域中有引入，在文獻[16]中，作者首次使用Jacobian 正則化方法提高卷積神經網絡識別物體的魯棒性，使得模型在對抗攻擊上表現良好，模型具備有抗噪的能力。在文獻[7]中作者首次使用Jacobian 正則化在長時動力系統預測中，文中給出了誤差在下一次預測中如何得到抑制，并給出了完整的分析過程。本文簡述了Jacobian 正則如何在非線性動力預測系統中抑制噪聲。

由（12）式可知，我們要想對含有噪聲的輸入進行預測進而獲得輸出：

于是對f( xi+ εi)在xi處進行泰勒展開，得到：

要想只保留f( xi)項，則需要在迭代的過程中使得除f( xi)以外的項盡可能的小，考慮到εi較小，高階項忽略，所以有：

也即是Jacobian 正則項。通過（15）式可知，當Jacobian 矩陣的值為0 矩陣時，系動力系統中的誤差則不被傳遞到下一次預測當中，所以在原模型的基礎上加上Jacobian 正則?？紤]到（16）是矩陣，通常使用Frobenius 范式，則優化目標函數J( K,θ )為：

1.4 Total Least-Square線性模型

為更進一步減小模型誤差，需要對線性模型（10）進行擴展，文獻[15]首次將線性模型進行廣義化，文章指出公式（10）中的模型是不對稱的，由于zt，zt+1,t=0,1,2…,T 在時間上具有連貫性，所以zt，zt+1都有相同性質的誤差，且指出（10）式是（18）式的一種特殊情況，所以本文將原框架中線性模型（10）進行廣義化：

注意到（18）式在輸入數據X 上引入了誤差?X，這是因為在多步預測時，當前時刻的預測值（輸出值）又作為下一時刻的輸入，而當前時刻的數據又包含有誤差，所以對于多步預測數據時，模型必須考慮到輸入數據中包含誤差。其次在單步模型中，輸入也是包含有誤差的，因為本文算法架構中，線性模型是對經過encoder 編碼之后的數據來建模的，所以單步模型也能適用。對于求解模型（18）中的方程即求解總體均方誤差（Total Least Square，TLS）的方法在文獻[15，17]中給出詳細介紹，利用奇異值分解的方法求解該問題。TLS模型作為最小均方誤差的擴展模型在文獻[9，10，11]中有廣泛應用。在文獻[15]中首次應用到動態模式求解算法中。對于求解（18）式，一般而言是先是對（18）式中右邊展開移項后得到：

對增廣矩陣Z ∈R2l×T奇異值分解，得到分解矩陣U2l×2l,Σ2l×T,VT×T，只保留矩陣V 的前l 個列向量Vl，計算得到投影后的時序數據：

再將（21）式的按照（8）式求解，轉化成求解偽逆問題即可[15]。

1.5 算法框架

本文算法在實際計算步驟分為兩部分，一是Autoencoder 的將數據的編解碼，二是使用TLS 求壓縮后求解壓縮后數據的線性表達式。所以本文在文獻[7]的基礎上加上Jacobian 正則和TLS，提出精確深度模式分解算法（Exact Dynamic Mode Decomposition，ED-DMD）來建模非線性動力系統。

算法1：精確深度模式分解算法

輸入：時序數據X=[x0,x1,x2,x3,…,xT]

輸出：模型參數P={ }

θenc,K,θdec

初始化：模型參數θ={ }θenc,θdec，迭代次數Iter，學習率η，因子λ，編碼維度l

重復直到迭代次數為Iter:

階段1：固定參數θ，更新矩陣K

1） 使用參數θ，使用當前encoder 編碼時序數據數據，Z=fenc( X )；則Z=[ z0,z1,z2,z3,…,zT]

2） 組織數據：

3） 對Zt奇異值分解，得到U,Σ,V

4） 更新矩陣K:K=Zt+1VΣ-1U*

階段2：固定參數K，更新網絡參數θ

1） 計算損失函數fl(x,θ)

2） 計算Jacobian 矩陣的F 范式JF(x,θ)

3） θ=θ-η?θ(fl+λJF)

算法1 結束。

由整個算法的過程可知，更新神經網絡參數θ 和K 是交替迭代的，更新其中一項參數時，另一項參數保持不變。參數θ 最先是給定初始值，算法交替更新參數時首先更新的是參數K，然后更新神經網絡參數θ，交替迭代的過程也可以由圖2 看出。為對比較時間復雜，假設輸數據維度n，數據量T，m 層encoder 神經網絡，每層神元數依次l1,l2…lm，且li?T(i=1,2…m)，對應的decoder 的m 層網絡神經元數為lm…l2,l1。在訓練數據相同的情況下，對于原算法[14]提出模型，第i(i=1,2…m)層神經元計算復雜度為O(T*li-1*li)，對于求解線性部分是求解lm個最小二乘回歸，其時間復雜度O(lm*lm2*T)，所以總的時間復雜度O(lm3*T)，其中lm個最小人可以并行求解。本文中模型計算量大的部分是求解總體最小二乘線性模型和計算Jacobian 矩陣。求解TLS 主要是對矩陣Z2lm×T求解SVD，選取l(l ≤lm)個特征向量，時間復雜度為O(T2*2l)，計算Jacobian 矩陣時間復雜度O(T*lm2)，總的時間復雜度為O（T2*2l）。

圖2 算法流程圖

2 實驗結果

本文中采兩個常用的常微分方程組，分別從這方程組中采取一定長度的軌跡數據，然后用數據來訓練模型，直至模型收斂完成訓練；在使用模型時只給出初始值x0，作為訓練模型的輸入數據，后續迭代按照（2）式所示進行。用表示模型循環預測的結果，xi表示采集第i 步的真實數據。最后對比訓練好的模型預測數據與真實數據。特別地，在使用模型時，單個數據的誤差為：

所有數據的總體誤差為：

系統維度n，數據長度為T。在本文模型中，模型參數還需要滿足總體誤差Esum盡可能的小。由于僅從減小Jacobian 誤差和總體均方誤差最小化的角度來建模，以期望達到長時間精準建模動力系統。所以并沒有明確的將長期誤差考慮到優化目標函數（17）之內，且（23）式難以優化。于是采取在迭代過程中取出滿足（23）式的模型參數作為最終模型，以達到最好的實驗效果。

2.1 范德波（Van der Pol）方程

第一個實驗采用的一個電路系統中常用的方程，可以使用如下方程組描述：

方程（24）描述的即是一個二維微分動力系統。參數u 取值u=4.0，在采集方程（24）的數據時，使用系統初始參數為[x0,y0]T=[-1,-4]T作為系統的初始狀態。時間間隔Δt=0.01，從t=0 時刻開始采集數據，時間步為3000 步。神經網絡采用自編碼機，編碼部分的神經網絡采用三層神經網絡，分別使用256，64，32 個神經元。解碼器部分使用與之對應的神經元數目。實驗結果如圖3 所示。

圖3

左側系統狀態參數隨時間變化，右側為系統的相圖?？梢钥吹侥Ｐ皖A測得到的數據在整體上系統還是大致接近真實系統數據，作為第一個闡述性的例子，既描述出了系統參數圖，系統參數隨時間變化的值（真實值和預測值）；又描述出了系統相圖。模型在整體上預測還是較好，但是在局部還是有誤差，如圖3。

2.2 Glycolysis動力系統

Glycolysis 動力系統方程組是在生物化學領域用于描述糖酵過程，該系統描述了7 個系統狀態參數的演化過程。這一組方程常常被用來做基準測試動力模型算法。在文獻[14，21，22]都用方程（25）來測試模型，所以：

本文也用它來測試模型的表達能力。數據采集時間間隔0.01，從t=0 時刻開始采集數據，采集4000步。神經網絡編碼器網絡層數為：7，256，128，64，32，32，32。解碼器有對應的網絡層數。從實驗結果可以看出，長期預測結果得到不錯性能；整體上，本文模型更接近真實數據，由圖4 的誤差對比可以看出，本文模型誤差在整體上還是比較小，模型性能有所提升。需要注意的是在算法在某些局部空間誤差較大，如圖4。

圖4 模型預測結果誤差對比

但是模型在整體上比原模型更接近真實數據如圖5。以上兩個實驗結果初步證明的本文提出的改進的算法在兩個微分方程組所確定的數據集上取得不錯效果。尤其是對Glycolysis 動力系統的改進效果更加明顯。原文算法[14]在僅僅能預測100～500 步，本文算法減小了前一步誤差對后續的影響，所以如實驗中的數據所示本文提出的算法達到了3000～4000 步的精確預測。

圖5 兩模型系統參數圖和真實數據對比

在模型的訓練過程中，模型的損失函數隨著時間變化如圖6 所示。然而本文算法對三維Lorenz 系統卻表現不佳，本文模型能預測500 步左右的狀態參數，有很多文獻使用Lorenz 作為驗證算法模型[6，8]，均不能取得較好的結果，這是由于該系統的Lyapunov 系數為正，微小的系統參數的誤差將在后續模型預測中導致巨大誤差，誤差呈現指數級的增長。相比較于D-DMD算法，雖然本文提出的ED-DMD 算法收斂速度較慢，這是由于新增了Jacobian 正則和總體最小二乘求解運算，增加了計算時間，導致收斂速度較慢；但是隨著迭代次數的增加，兩個算法都趨于收斂。本文實驗在NVIDIA GTX 1060 6GB 顯卡和TensorFlow 框架完成。

圖6 模型loss值隨時間（單位：秒）的變化

本文算法在兩個模型中取得較好的結果，誤差曲線總體來說也比較低，分析了算法的時間復雜度和模型loss 值隨時間的變化，由圖6 可以看出本文算法時間消耗沒有太大差距。實驗結果也證明了算法的優越性。

3 結語

本文提出的算法旨在揭示一般系統的演化規律，即非線性動力學。這些一般的演化規律能夠由一組確定的微分方程所表示，系統中的數據是由確定的但未知的規律所確定的。本文從系統數據入手，采用機器學習方法來揭示系統內在的規律，提出的方法能夠在一定的時間步中學習到數據的內在演化規律，并在長時預測商保持較低的模型誤差，實驗結果也證明這點。然而模型還有不足的地方需要改進，例如模型沒有明確的建模長期動力系統的關系[8]，模型沒有學習到數據長時的內蘊模式[13]，模型沒有使用循環神經網絡來維持數據之間的時序關系[17-18]等；這些問題希望能在后續的工作中得到解決。