一種改進的基于深度神經(jīng)網(wǎng)絡的偏微分方程求解方法*

陳新海,劉 杰,萬 仟,龔春葉

(1.國防科技大學并行與分布處理國家重點實驗室，湖南 長沙 410073；2.復雜系統(tǒng)軟件工程湖南省重點實驗室，湖南 長沙 410073)

1 引言

隨著計算能力的提升，數(shù)值分析在諸多科學與工程領域發(fā)揮著越來越重要的作用[1]。大部分數(shù)值分析問題的底層數(shù)學模型都可以歸結(jié)為對偏微分方程的數(shù)值求解。傳統(tǒng)的數(shù)值求解方法主要包括有限差分法、有限元法和有限體積法[2]。這些方法首先將計算域離散為獨立的網(wǎng)格單元，然后在單元子域上迭代求解偏微分方程系統(tǒng)，以此獲得未知函數(shù)的數(shù)值近似解。然而，為了保證數(shù)值分析的準確性和有效性，傳統(tǒng)方法通常十分耗時且強依賴人工經(jīng)驗：一方面，迭代求解大型線性/非線性方程組需要昂貴的計算開銷；另一方面，為了避免計算失效，通常在網(wǎng)格剖分階段還需要頻繁的人機交互來判別和優(yōu)化計算網(wǎng)格的質(zhì)量，以滿足求解器和模擬精度的要求[3,4]。隨著所分析的問題越來越復雜，繁重的交互和開銷限制了傳統(tǒng)數(shù)值求解方法在實時模擬、氣動參數(shù)設計和優(yōu)化空間探索等問題上的應用效率。

近年來，以深度神經(jīng)網(wǎng)絡為代表的人工智能技術取得了令人矚目的突破和成果。以往的研究已證明，深度神經(jīng)網(wǎng)絡是復雜知識系統(tǒng)和函數(shù)擬合的有效替代[5,6]。通過訓練由神經(jīng)元組成的深度神經(jīng)網(wǎng)絡，其模型能夠從高維的非線性空間中學習到潛在的規(guī)則和特征，最終擬合得到輸入與輸出之間的函數(shù)映射關系。深度神經(jīng)網(wǎng)絡在各個領域的應用被認為是“解放勞動力”的過程[2]。受啟發(fā)于深度神經(jīng)網(wǎng)絡在圖形圖像、自然語言處理等領域的廣泛應用，深度學習方法越來越多地應用到受限于高額計算開銷和領域人工經(jīng)驗的數(shù)值分析中。

Chen等[7,8]證明了神經(jīng)網(wǎng)絡的萬能逼近定理，即由包含非線性激活函數(shù)(如Tanh、ReLU、Sigmoid激活函數(shù))的隱藏層構(gòu)成的前饋神經(jīng)網(wǎng)絡能夠近似任何從一個有限維空間到另一個有限維空間的Borel可測函數(shù)。該定理隨后被進一步推廣到了微分算子的近似，為使用神經(jīng)網(wǎng)絡求解偏微分方程奠定了理論基礎?；谠摾碚摚琈alek等[9]提出了一種求解高階微分方程的混合神經(jīng)網(wǎng)絡模型，該模型使用了Nelder-Mead方法來構(gòu)造常微分方程的封閉解，并在微分方程的初始/邊界值問題上取得了很好的效果。Mall等[10]通過構(gòu)造滿足邊界條件的損失項和神經(jīng)網(wǎng)絡模型，提出了一套求解常微分方程的回歸算法，同時還研究了神經(jīng)元數(shù)量和初值對預測解精度的影響。但是，這些方法都只適用于常微分方程求解。

Lagaris等[11]提出了適用于不規(guī)則邊界的基于神經(jīng)網(wǎng)絡的偏微分方程求解模型，該模型使用多層感知機和徑向基神經(jīng)網(wǎng)絡來擬合復雜邊界，并在一些較為簡單的二維、三維偏微分方程上取得了較高的預測準確率。Sirignano等[12]提出了基于長短時記憶LSTM(Long Short-Term Memory)網(wǎng)絡的偏微分方程求解方法，用于求解高維的偏微分方程。但是，以上方法都只適用于特定類型的偏微分方程，缺乏一定的通用性。2019年，Raissi等[13]提出了物理信息神經(jīng)網(wǎng)絡PINN(Physics-Informed Neural Networks),用于近似偏微分方程中的待求解函數(shù)。該方法利用深度神經(jīng)網(wǎng)絡的函數(shù)擬合能力，用控制方程組和邊界函數(shù)來約束網(wǎng)絡中可調(diào)參數(shù)(權(quán)重和偏置)的迭代更新過程，最終得到輸入的時空點坐標到待求函數(shù)的端到端映射模型。整個求解過程不需要進行網(wǎng)格離散以及先驗的數(shù)值求解，且訓練后的模型能夠在極小計算開銷下自動預測計算域內(nèi)任意點的觀測值，這大大提高了數(shù)值分析的效率。但是，傳統(tǒng)的物理信息神經(jīng)網(wǎng)絡在求解許多基于偏微分方程的物理問題時存在預測精度低的情況。

為了進一步提高物理信息神經(jīng)網(wǎng)絡方法在求解實際問題時的準確性，本文提出了一種改進的基于深度神經(jīng)網(wǎng)絡的偏微分方程求解方法TaylorPINN。與傳統(tǒng)PINN中直接使用回歸模型輸出離散的預測值不同，TaylorPINN利用深度神經(jīng)網(wǎng)絡的函數(shù)逼近能力，并結(jié)合泰勒公式，通過構(gòu)造解的形式將神經(jīng)網(wǎng)絡輸出與待求函數(shù)聯(lián)系起來，實現(xiàn)了無網(wǎng)格的數(shù)值求解過程。在Klein-Gordon方程[14]、 Helmholtz方程[15]和Navier-Stokes方程[16,17]上的數(shù)值實驗結(jié)果表明，TaylorPINN能夠準確擬合計算域內(nèi)時空點坐標與待求函數(shù)之間的映射關系，并提供了準確的數(shù)值預測結(jié)果。與常用的基于物理信息神經(jīng)網(wǎng)絡方法相比，對于不同數(shù)值問題,TaylorPINN將預測精度提升3～20倍。

2 理論背景介紹

2.1 偏微分方程

為方便討論，設偏微分方程的一般形式如式(1)所示：

(1)

其中,D表示微分算子，u是待求解的函數(shù)，計算域Ω?Rd。方程(1)的邊界條件定義如式(2)所示：

B[u(x,t)]=h(x,t),x∈?Ω

(2)

初始條件如式(3)所示：

(3)

2.2 物理信息神經(jīng)網(wǎng)絡PINN

2019年，Raissi等[13]提出了物理信息神經(jīng)網(wǎng)絡PINN的概念，即用深度神經(jīng)網(wǎng)絡模型來求解偏微分方程控制的物理問題。在PINN中，偏微分方程的數(shù)值解問題被轉(zhuǎn)化為基于神經(jīng)網(wǎng)絡的最優(yōu)化問題，并通過搭建深度神經(jīng)網(wǎng)絡u(θ)作為黑盒回歸模型來擬合待求的函數(shù)u。在網(wǎng)絡訓練優(yōu)化過程中，PINN只使用方程組本身的信息來約束神經(jīng)網(wǎng)絡參數(shù)的優(yōu)化。具體地，用于指導優(yōu)化的損失函數(shù)Loss由3部分組成，分別是方程殘差項MSEu以及條件項MSEb和MSEI：

Loss=MSEu+MSEb+MSEI

(4)

(5)

(6)

(7)

通過最小化損失函數(shù)Loss，PINN能利用神經(jīng)網(wǎng)絡的函數(shù)擬合能力得到同時滿足方程項和條件項的待求解函數(shù)u的近似解u(θ)，即：

u≈u(θ)

(8)

該近似解由一個前饋神經(jīng)網(wǎng)絡表示，對于計算域中的任意點坐標，訓練后的網(wǎng)絡模型都能夠通過幾次矩陣乘高效地預測出對應的函數(shù)值，從而實現(xiàn)點坐標到物理場的端到端映射?；赑INN的訓練和預測過程消除了對計算網(wǎng)格的依賴，大幅提高了數(shù)值分析的效率。但是，現(xiàn)有的物理信息神經(jīng)網(wǎng)絡在求解許多基于偏微分方程的物理問題時存在預測精度低的情況，從而限制了該方法的實用性和有效性。

3 基于深度神經(jīng)網(wǎng)絡的偏微分方程求解方法TaylorPINN

為了進一步提升物理信息神經(jīng)網(wǎng)絡的預測精度，本文提出了一種改進的基于深度神經(jīng)網(wǎng)絡的偏微分方程求解方法TaylorPINN。該方法同樣將微分方程的求解問題轉(zhuǎn)化為函數(shù)的最優(yōu)化問題。但是，與傳統(tǒng)的物理信息神經(jīng)網(wǎng)絡不同的是，TaylorPINN在預測過程中加入了對待求解函數(shù)的形式構(gòu)造過程，即不再直接使用回歸模型輸出離散的預測值，而是通過構(gòu)造解的方式將神經(jīng)網(wǎng)絡輸出與待求解函數(shù)聯(lián)系起來。在這個過程中，如何構(gòu)造解的形式成為了問題轉(zhuǎn)化后的關鍵。

泰勒公式是一種廣泛用于函數(shù)逼近的方法，定理1是一元函數(shù)的n階泰勒展開定理：

定理1設函數(shù)f(x)在x=x0的鄰域內(nèi)n+1階可導，則對于位于此鄰域內(nèi)的任一點x有：

f(x)=f(x0)+f′(x0)(x-x0)+…+

(9)

根據(jù)上述定理，本文采用基于泰勒公式的多項式函數(shù)空間來構(gòu)造解的形式，即令偏微分方程中的每個待求解函數(shù)如式(10)所示：

(10)

相比于其他函數(shù)空間(如傅里葉級數(shù)或各類徑向基函數(shù))，由時空坐標(x,t)的線性加權(quán)及其乘積構(gòu)成的多項式空間擁有靈活的擴展性，在擁有強的函數(shù)表達能力的同時也易于實現(xiàn)，且計算開銷相對較小。

Figure 1 Architecture of TaylorPINN圖1 TaylorPINN的整體結(jié)構(gòu)

在模型架構(gòu)方面，TaylorPINN除了使用全連接層外，還引入了一個輸入增強層作為一種增強手段來擴大輸入點坐標的維度。該增強層實現(xiàn)了一個基于時空點坐標(x,t)的離散線性映射φ，φ如式(11)所示：

(11)

此外，在損失函數(shù)構(gòu)造方面，TaylorPINN引入了一個靜態(tài)權(quán)重參數(shù)α1來平衡不同損失項對迭代收斂的貢獻，使神經(jīng)網(wǎng)絡模型在訓練優(yōu)化階段就能夠更好地滿足各損失項，從而更快地收斂到最優(yōu)值。最終，TaylorPINN的損失函數(shù)如式(12)所示：

Loss=MSEu+α1(MSEb+MSEI)

(12)

圖1顯示了TaylorPINN的整體結(jié)構(gòu)，其結(jié)構(gòu)包含輸入層、輸入增強層、進行高維函數(shù)空間特征學習的隱藏層和最后的輸出層。其中，TaylorPINN的輸入層以計算域的時空點坐標為輸入，在擴維后傳播到由全連接層組成的隱藏層。前向傳播得到輸出值λ，作為基于泰勒公式的解形式中的系數(shù)來計算待求解函數(shù)的值。在訓練階段，得到的函數(shù)值用于計算損失函數(shù)的值及反向傳播的梯度。在預測階段，該輸出值將直接用于計算待求解函數(shù)的預測解。TaylorPINN的訓練過程如算法1所示。

算法1TaylorPINN訓練過程

輸入:計算域內(nèi)的時空點坐標。

輸出:訓練后的TaylorPINN模型。

Step1分別從計算域內(nèi)和邊界隨機選取時空點坐標，組成訓練樣本集。

Step2搭建TaylorPINN神經(jīng)網(wǎng)絡并根據(jù)待求解的偏微分方程構(gòu)建損失函數(shù)Loss。

Step3初始化網(wǎng)絡參數(shù)。

Step4隨機從樣本集中選取時空點坐標作為訓練樣本輸入TaylorPINN。

Step5使用優(yōu)化算法最小化損失函數(shù)Loss，并更新神經(jīng)網(wǎng)絡中的優(yōu)化參數(shù)。

Step6重復步驟4和步驟5，直到完成預定的迭代輪數(shù)。

4 數(shù)值實驗與分析

本節(jié)使用Klein-Gordon方程[14]、Helmholtz方程[15]和Navier-Stokes方程[16,17]3個數(shù)值案例來評估TaylorPINN的偏微分方程求解能力，并在相同網(wǎng)絡結(jié)構(gòu)、超參數(shù)及訓練迭代輪數(shù)的條件下與傳統(tǒng)PINN進行預測性能對比。實驗中使用的初始學習率均為0.001，衰減率為每1 000輪衰減0.9，優(yōu)化器選用Adam[18]，batchsize為128。實驗的測試平臺是Intel Core i7-6700和NVIDIA Tesla P100，實驗環(huán)境是TensorFlow1.15.0[19]。偏微分方程中的微分項統(tǒng)一用TensorFlow中的自動微分函數(shù)“tf.gradient()”[20]。神經(jīng)網(wǎng)絡模型的預測精度使用L2誤差來衡量，其計算方式如式(13)所示：

(13)

其中,uref是待求方程的參考解，通常由解析解或高精度數(shù)值解表示;upred是神經(jīng)網(wǎng)絡提供的預測解。為了保證實驗結(jié)果的可復現(xiàn)性，本文將訓練時的隨機種子固定為666。

4.1 Klein-Gordon方程

在第1個測試案例中，本文使用Klein-Gordon方程來評估所提出方法的待求函數(shù)擬合能力。Klein-Gordon方程是量子物理、固體物理等領域的基礎偏微分方程之一[14]。一維Klein-Gordon方程的形式如式(14)所示：

utt+αuxx+βuγ=q(x,t),

(x,t)∈Ω×[0,T],Ω=[0,1]

(14)

邊界條件如式(15)所示：

u(x,t)=h(x,t),(x,t)∈?Ω×[0,T]

(15)

初始條件如式(16)和式(17)所示：

u(x,0)=g1(x),x∈Ω

(16)

ut(x,0)=g2(x),x∈Ω

(17)

當設g1(x)=g2(x)=0,α=β=1,γ=3時，本文可以構(gòu)造出方程的解析解作為神經(jīng)網(wǎng)絡預測的參考解，如式(18)所示：

uref(x,t)=xcos(5πt)+(xt)3

(18)

本文使用該解析解反推得到h(x,t)和q(x,t),用于后續(xù)計算損失函數(shù)中的殘差項和條件項。對于偏微分方程(式(14)～式(17))，本文構(gòu)建如式(19)所示的損失函數(shù)用于后續(xù)的神經(jīng)網(wǎng)絡訓練：

(19)

Figure 3 Prediction results of TaylorPINN and PINN on Klein-Gordon equation圖3 TaylorPINN和PINN在Klein-Gordon方程上的預測結(jié)果

在本案例中，本文搭建了一個包含3個隱藏層的TaylorPINN網(wǎng)絡模型，每層神經(jīng)元個數(shù)為100。在訓練階段，本文通過最小化損失函數(shù)來尋找最優(yōu)的網(wǎng)絡參數(shù)，以擬合Klein-Gordon方程的封閉預測解。本文令訓練的總輪數(shù)為40 000輪，靜態(tài)權(quán)重α1為10。測試集為計算域內(nèi)均勻分布的10 000個時空點。不同泰勒公式階數(shù)下，TaylorPINN的預測精度(L2誤差)如圖2所示。

Figure 2 Prediction performance of TaylorPINN on Klein-Gordon equation with different orders圖2 不同階數(shù)下TaylorPINN在Klein-Gordon方程上的預測性能

從圖2中可以看出，隨著階數(shù)的增加，預測精度也逐漸增加，當階數(shù)k=5時，模型實現(xiàn)了最佳的預測性能，L2誤差為6.63E-03;但是，當繼續(xù)增加泰勒公式階數(shù)時，精度開始下降，其主要原因是由于階數(shù)越大，輸出層的神經(jīng)元個數(shù)也就越多，對應模型收斂所需要的訓練輪數(shù)、隱藏層寬度和深度也隨之增大。因此，當固定模型結(jié)構(gòu)和訓練輪數(shù)后，TaylorPINN會在階數(shù)較大時出現(xiàn)欠擬合的情況，從而影響模型收斂，導致預測精度下降。

圖3和表1對比了PINN和TaylorPINN(k=5)在Klein-Gordon方程上的預測結(jié)果。其中圖3a是由式(18)給出的參考解的可視化結(jié)果，圖3b和圖3d是由2個模型給出的預測解的可視化結(jié)果，圖3c和圖3e是參考解與預測解之間的絕對誤差。從圖3可以看出，TaylorPINN預測解與參考解基本吻合，且預測精度明顯高于PINN的。從表1可以看到，在L2誤差方面，雖然TaylorPINN在構(gòu)建解的形式時會帶來額外的計算開銷，但相比于PINN,TaylorPINN將預測精度提升大約20倍。

Table 1 Comparison of the performance of TaylorPINN with PINN on Klein-Gordon equation

4.2 Helmholtz方程

在第2個測試案例中，本文使用TaylorPINN來求解二維Helmholtz方程。該方程在聲學、電磁學和彈性力學等領域有著廣泛的應用[15]。二維Helmholtz方程如式(20)所示：

Δu(x,y)+k2u(x,y)=q(x,y),(x,y)∈Ω

(20)

方程的邊界條件如式(21)所示：

u(x,y)=h(x,y),(x,y)∈?Ω

(21)

其中，Δ是拉普拉斯算子，計算域Ω為[0,1]×[0,1]。當k=1，源項q(x,y)如式(22)所示時，

q(x,y)=-π2sin(πx)sin(4πy)-

k2sin(πx)sin(4πy)

(22)

可以構(gòu)建出方程的解析解,如式(23)所示:

uref=sin(πx)sin(4πy)

(23)

在TaylorPINN中，Helmholtz方程對應的損失函數(shù)如式(24)所示：

(24)

在本案例中，本文通過改變網(wǎng)絡層數(shù)和單層神經(jīng)元個數(shù)，測試了不同網(wǎng)絡規(guī)模對TaylorPINN預測性能的影響。為了更好地分析不同規(guī)模下的性能結(jié)果，本文將TaylorPINN與相同規(guī)模下的PINN進行了比較。圖4顯示了當階數(shù)為3，固定單層神經(jīng)元數(shù)為50時，TaylorPINN在Helmholtz方程案例上的預測結(jié)果。

Figure 4 Prediction performance of TaylorPINN on Helmholtz equation with different numbers of network layers圖4 不同網(wǎng)絡層數(shù)的TaylorPINN在Helmholtz方程上的預測性能

從圖4可以看出，當網(wǎng)絡層數(shù)為3時，TaylorPINN實現(xiàn)了最佳的預測性能，其L2誤差為1.86E-02。當繼續(xù)增加網(wǎng)絡層數(shù)時，預測精度反而會有所退化。但在不同網(wǎng)絡層數(shù)下，TaylorPINN都獲得了比PINN更高的預測性能，且預測精度最大可提升12倍。

表2顯示了當固定網(wǎng)絡層數(shù)為3，變化單層神經(jīng)元數(shù)時PINN和TaylorPINN的L2誤差。從表2可以看出，隨著神經(jīng)元個數(shù)的增加，2個網(wǎng)絡模型都能得到預測精度提升。當單層神經(jīng)元個數(shù)為50時，相比于PINN,TaylorPINN的預測精度提升了大約8倍，其可視化結(jié)果如圖5所示。

Table 2 Influence of different numbers of neurons on the prediction performance of TaylorPINN

4.3 Navier-Stokes方程

在最后一個案例中，本文使用一個經(jīng)典的流體力學案例——二維頂蓋驅(qū)動流(如圖6所示)[16,17]，來分析TaylorPINN的偏微分方程求解能力。

Figure 5 Prediction results of TaylorPINN and PINN on Helmholtz equation圖5 TaylorPINN和PINN在Helmholtz方程上的預測結(jié)果

Figure 6 Lid-driven cavity flow case圖6 頂蓋驅(qū)動流案例

該案例的控制方程為不可壓Navier-Stokes方程，如式(25)所示：

(25)

(26)

u(x,y)=(1,0),(x,y)∈Γ1

(27)

u(x,y)=(0,0),(x,y)∈Γ0

(28)

其中，u表示待求解的速度矢量場;p是標量壓力場;Re是流體的雷諾數(shù);計算域Ω為[0,0.1]×[0,0.1];Γ1是方腔的頂蓋邊界，其x方向速度分量為1 m/s;Γ0是其他3個邊界，速度分量均為0。本文使用OpenFOAM[21]計算的該案例的數(shù)值解作為參考解。

在TaylorPINN的訓練過程中，本文使用Navier-Stokes方程的渦量-速度形式VV(Vorticity-Velocity)來構(gòu)建損失函數(shù)，在該形式下，u的速度分量(u,v)可以使用函數(shù)ψ來計算,如式(29)所示：

(29)

根據(jù)式(29)，方程自然滿足不可壓流體的連續(xù)性。此外，本文還去掉了時間項來進一步簡化模型。最終構(gòu)造的用于訓練的Navier-Stokes方程損失函數(shù)如式(30)所示：

Loss=

(30)

根據(jù)邊界條件的定義，對于Γ1上的采樣點，h1(x,y)=(1,0),對于Γ0上的采樣點，h0(x,y)=(0,0)。值得注意的是，該案例給定的邊界條件在方腔的2個邊界交界處(見圖6中的黑色標記點)形成了不連續(xù)的奇點，即左右邊界條件給定的交界點處的速度值為0，而上頂蓋將交界點處的速度賦值為1。這種不連續(xù)性會影響神經(jīng)網(wǎng)絡的訓練優(yōu)化。因此，在本案例的樣本點隨機選取過程中，本文避免使用這2個間斷點用于網(wǎng)絡訓練。

最終，本文搭建了一個包含3個隱藏層(單層神經(jīng)元個數(shù)為100)的TaylorPINN。圖7展示了當Re為100，泰勒展開式階數(shù)為5時，TaylorPINN和PINN在頂蓋驅(qū)動流案例上的預測結(jié)果。從圖7可以看出，相比于PINN，TaylorPINN能夠更準確地捕捉到流體流動的狀態(tài)。表3中的結(jié)果顯示，TaylorPINN在該案例上實現(xiàn)了7.76E-02的預測誤差，而相同網(wǎng)絡規(guī)模下PINN的L2誤差為4.17E-01。

Figure 7 Prediction results of TaylorPINN and PINN on Navier-Stokes equation圖7 TaylorPINN和PINN在Navier-Stokes方程上的預測結(jié)果

5 結(jié)束語

隨著人工智能理論和技術的發(fā)展，使用深度神經(jīng)網(wǎng)絡技術提升傳統(tǒng)數(shù)值分析的效率成為了當前學術界的研究熱點。偏微分方程求解是數(shù)值分析的計算核心，本文提出了一種基于深度神經(jīng)網(wǎng)絡的偏微分方程求解方法TaylorPINN。該方法結(jié)合了神經(jīng)網(wǎng)絡和泰勒公式的函數(shù)擬合能力，實現(xiàn)了無網(wǎng)格的高效數(shù)值求解過程。在不同偏微分方程案例上的數(shù)值實驗結(jié)果表明，TaylorPINN能夠很好地擬合計算域內(nèi)時空點坐標與待求函數(shù)值之間的映射關系，并提供準確的數(shù)值預測結(jié)果。在未來的工作中，考慮在網(wǎng)絡模型上做進一步研究，例如將注意力機制引入網(wǎng)絡模型，從而提高模型的準確率；也會改進提出的方法并將其應用到更復雜的物理問題場景中，以進一步增強該類方法的實用性。

Table 3 Comparison of the performance of TaylorPINN with PINN on Navier-Stokes equation