基于內嵌物理機理神經網絡的熱傳導方程的正問題及逆問題求解

趙 暾，周 宇，程艷青，錢煒祺，*

(1. 空氣動力學國家重點實驗室，綿陽 621000；2. 中國空氣動力研究與發展中心，綿陽 621000)

0 引言

自從Hinton[1]在2006年提出深度學習的概念之后，最近幾年深度學習在科學領域和商業領域已經取得了巨大的成功。一般來說，深度學習需要大量的訓練樣本，通常應用在大數據領域，而且由于過度依賴訓練數據，導致泛化能力差，同時還具有以下缺點：缺乏推理和因果關系的表達能力而導致表征能力弱；無法解釋形成決策的推理過程；很難在計算過程中導入選擇性注意力機制。

在實際工程領域，受成本限制，所獲取到的樣本數據量是有限的。因此，很難將深度學習技術直接應用于小樣本問題領域。但是這些數據通常都包含由物理定律、經驗規則或者領域知識構成的先驗信息。如何將深度學習方法有效運用到包含先驗信息的小樣本數據領域，相關學者做了大量研究。其中，美國布朗大學Karniadakis教授團隊做出了較為出色的工作。他們通過將蘊藏在小樣本數據中的先驗信息編碼到深度網絡中，提出了一種新的深度學習算法—內嵌物理機理神經網絡[2]（Physics-informed neural networks，PINNs），該網絡可以由神經網絡對輸入及模型參數進行自動微分所獲取。PINNs通常受到數據中包含的物理定律的約束，而這些物理定律一般可以建模成時間相關的非線性偏微分方程。這樣就可以利用PINNs發展出一類新的偏微分方程數值求解方法以及一類新的系統辨識和模型反演的數據驅動方法。PINNs方法的核心思想是將上述原問題（即求解偏微分方程的正問題或者逆問題）轉換為優化問題求解。與傳統的偏微分方程數值求解方法以及系統辨識方法相比，PINNs具有以下幾個特點：

1）PINNs可以作為訓練辨識一體化的方法，即在求解偏微分方程的同時，根據附加信息直接辨識出方程中的未知參數。

2）與傳統的基于網格的偏微分方程數值求解方法相比，PINNs是無網格的，這樣可以節省大量的網格繪制與生成時間[3]，而且還可以避免求解高維PDEs時會遇到的維度災難問題。

3）與傳統的基于變分原理求解偏微分方程方法相比，PINNs可以求解不滿足變分原理的偏微分方程[4]。

4）由于PINNs采用了自動微分，因此可以避免傳統數值方法中由差分格式帶來的誤差影響。而且自動微分是直接對神經網絡的輸出而不是對數據進行微分，因此，PINNs對包含噪聲的輸入數據具有較強的魯棒性。

5）PINNs可以有效地解決不適定的偏微分方程逆問題。

6）由于PINNs中的損失函數是非凸的，會存在多個局部極小值，因此基于梯度下降的優化算法幾乎一定會陷入某個極小值，從而導致PINNs的求解結果不唯一。

7）PINNs對于每個新的實例（例如，方程參數變化）就需要重新訓練網絡。后續，Lu等[5]提出了DeepONet去解決這個問題。該網絡可以直接對數據中的非線性算子進行近似，例如，如果數據是從偏微分方程獲得，則可以直接學習出積分算子，從而近似出偏微分方程的解。

Jagtap等[6]將PINNs中的激活函數修改為自適應激活函數，大大加速了訓練收斂速度。具體做法是在激活函數中引入一個超參數，并將該參數作為神經網絡參數進行優化。在PINNs的理論基礎層面，Shin等[4]針對橢圓型和拋物線型偏微分方程，證明了PINNs訓練得到的網絡可以一致收斂到偏微分方程的解。

經過最近兩年的發展，PINNs發展出了一大類家族—cPINNs[7]、fPINNs[8]、nPINNs[9]、sPINNs[10-11]、vPINNs[12]以及XPINNs[13]等。前面四種可以分別用來求解離散域非線性守恒律方程、分數階偏微分方程、非局部橢圓型偏微分方程、隨機偏微分方程，vPINNs則是以Petrov-Galerkin方法為基礎求解偏微分方程，而XPINNs可以求解任意復雜幾何域上的非線性偏微分方程。

Lu等[14]基于Tensorflow開發出了一個Python工具包—DeepXDE，該工具包對上述算法進行了整合，并且還可以用來進行多保真度數據建模，以及從數據中直接學習非線性算子和函數關系。更進一步地，PINNs由于其原理簡單，靈活度高，效果良好，已被英偉達、ANSYS、西門子等公司采用。英偉達已經開發出了SimNet[15]工具包，用來求解偏微分方程的正問題、逆問題以及數據融合問題。

雖然PINNs方法目前發展迅速，但是對求解正問題時的精度以及求解逆問題時的魯棒性還需要進一步分析。本文則針對一維熱傳導方程對上述問題進行了研究。

1 PINNs的基本結構

下面以求解非線性偏微分方程的正問題和逆問題為例，闡釋PINNs的基本網絡結構及其訓練方式。

非線性偏微分方程的一般形式可以表示為:

其中，u是上述偏微分方程的解， N[u,λ]為帶有參數λ的非線性算子， Ω 是RD的 子集，tT為終端時刻，h(x)為方程的初始條件，g(x,t)為方程的邊界條件。

首先，建立一個深度網絡n(x,t;θ)去近似u，并定義偏微分方程的殘差為PINN的輸出:

上式中的偏導數都可以通過自動微分得到。相應地，損失函數可以給定為:

損失函數的訓練集分為三部分，第一部分為初始條件和邊界條件的采樣點集，第二部分為對偏微分方程殘差進行采樣的配置點集。配置點集一般有兩種選取方式，第一種為類網格點選取方式，第二種為偽隨機離散取點方式，如圖1所示。當針對高維非常規區域上偏微分求解問題，偽隨機離散取點方式具有更強的靈活性。第三部分為已知的數據樣本集，該部分標簽數據用來對方程中的參數進行辨識。

圖1 配置點離散方式Fig. 1 The distribution of collocation points

從損失函數式（3）的組成可以看出，偏微分殘差項起到了正則化作用，這樣可以增強網絡n(x,t;θ)的泛化能力。而從另一個角度看，通過自動微分構建偏微分方程殘差并將其引入到損失函數的過程，可以看作是將該偏微分方程表示的物理規律嵌入到深度網絡的過程。

最后，基于該損失函數，通過梯度下降方法對深度網絡參數θ和方程參數λ進行訓練。

最終，通過訓練就可以求得使損失函數（3）最小的 參 數 θ?和 λ?，而 方 程（1）的解可以用神經網絡=n(x,t;θ?)來表示。PINNs的網絡結構及訓練過程可以用圖2表示。

圖2 PINNs基本原理圖Fig. 2 A sketch of PINNs

2 PINNs方法的應用

2.1 一維熱傳導方程求解及其精度分析

在物體邊界給定為第二類邊界條件后，可得到如下的一維熱傳導控制方程[16]：

邊界條件：

初始條件：

本節的目的是在已知邊界條件和初始條件的情況下，利用PINN求解導熱系數k= 1時方程（5）的解。

首先設計一個深度網絡n(x,t;θ)去近似方程（5）的解，即：

該網絡包含2個輸入，1個輸出，3個隱層且每個隱層包含40個節點，隱層節點的激活函數設置為sigmoid。隨后對n(x,t;θ)的輸出進行自動微分，定義PINN為：

接著，選擇如下均方誤差的和作為損失函數：

其中，Nb和N0分別為邊界點和初始點訓練集大小，Nf為方程配置點訓練集大小。在本文中，首先給定Nb和N0分 別為80和100，以及在時間空間域 Ω×R1內采用類網格點方式選擇Nf= 8000個配置點。接著，采用Adam優化器并設置學習率為0.001對所建的網絡n(x,t;θ)進 行訓練。即在訓練過程中，通過給定邊界條件數據集、初始條件數據集{ti,xi;Ti}以及方程殘差點集 {ti,xi;fi}，通過優化尋找最接近觀測數據集以及偏微分方程殘差的參數 θ?所確定的深度網絡。最后，利用訓練完的網絡預測區域(x,t)∈[0,1]×(0,1]內的溫度值。

還需要特別注意的是，由于PINN會從不同的網絡初值收斂到不同的解，因此，在本文中，隨機初始化20次訓練PINN，損失函數隨迭代步數變化情況如圖3所示。

圖3 20次不同初始權重網絡的損失函數隨迭代步數變化Fig. 3 Loss functions with 20 different initial weights versus the iteration step

從圖3可以看出，從不同的初始權重出發，網絡經過大約10000次左右的迭代訓練后，損失函數都可以收斂。

然后利用訓練完的網絡分別預測溫度場，最后將平均值作為最終結果，并將其與溫度場的解析解進行對比，如圖4所示。與解析解相比，利用PINN預測整體溫度場的均方誤差約為2.2865× 10?7。

圖4 基于PINN的一維熱傳導方程求解結果與解析解對比Fig. 4 Comparison of the solution obtained by PINN and the analytical solution of the one-dimensional heat conduction equation

通過對PINN方法分析可知，其求解誤差主要包括有限訓練點集導致的采樣誤差、神經網絡的逼近誤差以及陷入局部最優解導致的優化誤差等三部分組成。下面在神經網絡結構固定的情況下，即逼近誤差一定，分別從空間和時間兩個尺度分析采樣誤差以及優化誤差對精度的影響，并且與控制體積方法的求解結果進行對比分析。為分析空間尺度影響，采用以下配置點選取方式：固定100個時間點，分別選取20、40、80、160、320、640和1280個空間點；而為分析時間尺度影響，則采用以下配置點選取方式：固定80個空間點，分別選取25、50、100、200、400、800和1600個時間點，具體如表1和表2所示。

表1 分析空間尺度離散影響的采樣點選取方式Table 1 The selection method of sampling points to analyze the influence of spatial discretization

表2 分析時間尺度離散影響的采樣點選取方式Table 2 The selection method of sampling points to analyze the influence of temporal discretization

本文訓練環境為英偉達RTX2060顯卡，顯存6 GB，運行頻率1800 MHz。在該環境中隨機初始化20次訓練網絡，訓練時間隨采樣點數的變化如圖5所示?？梢钥闯觯S著采樣點數的增加，訓練時間呈線性增長的趨勢。

圖5 訓練時間隨采樣點數的變化Fig. 5 The training time as a function of the sampling number

利 用 訓 練 完 的 網 絡 計 算 其 在x= 0 m、0.1 m、0.25 m、0.6 m及1 m處的平均預測值及其對應的均方誤差，并與控制體積的求解結果相對比，如圖6所示，圖中的黑色虛線表示的是O(N?2)收斂速率，其中N為采樣點數。

圖6 PINN方法與控制體積方法求解精度比較Fig. 6 CA comparison of the solution accuracy of between the PINN and the control volume method

通過對圖6的分析，可以得出以下幾個結論：

1）當采樣點數較少時，時間和空間尺度采樣點數對PINN方法的求解精度影響與控制體積方法類似，即在空間采樣點固定情況下，求解誤差按照O(N?2)速率收斂，而且時間尺度的采樣誤差相較空間尺度占主要地位。

這是由于在損失函數（10）中，時間尺度采樣點對應的有左邊界殘差和右邊界殘差兩項，而空間尺度采樣點只對應一項初始條件殘差。

2）PINN方法在采樣點數較少時，采樣誤差為主要誤差來源，隨著采樣點數增加，其求解精度隨之提高，但是當總采樣點數約大于8000時，兩種配置點采樣方式的求解精度幾乎收斂到同一水平。

這是由于過多的采樣點數會導致過于復雜的損失函數，從而使優化誤差增大，這讓優化誤差成為主要誤差來源。

下面，固定8000個配置點，分別分析初始和邊界條件的采樣點對求解精度的影響。左右邊界分別選取100個采樣點，初始條件分別選取20、40、80、160、320、640和1280個采樣點；初始條件選取80采樣點，左右邊界分別選取25、50、100、200、400、800、1600個采樣點時，方程（5）求解精度如圖7所示。

從圖7可以得出與圖6一致的結論，即采樣點數少時，采樣誤差主導了收斂誤差，當采樣點數增加到一定程度后，優化誤差成為了收斂誤差的主要組成部分。

圖7 配置點固定時PINN方法求解精度分析Fig. 7 The solution accuracy of the PINN with fixed collocation points

2.2 一維熱傳導方程參數辨識及其魯棒性分析

本節針對熱傳導方程（5），在已知初始條件、混合邊界條件：

以及Nt個測點溫度的情況下，對方程中的導熱系數k進行辨識，其中，T0為基于2.1節方程的解析解。通過對上節方法稍加改動，就可以對參數k進行辨識。具體實現步驟如下：

首先，在構建PINN時，將參數k作為待優化變量，即PINN的輸出為：

接著，將Nt個測點溫度作為標簽數據，定義如下損失函數：

最后，通過Adam算法對所構建的網絡進行優化，得到使損失函數最小的k?。

本文選取網絡結構如下：2個輸入及1個輸出，6個隱層且每個隱層包含20個節點，隱層節點的激活函數設置為tanh。在固定網絡初始參數的情況下，同時給定無噪聲、10%噪聲以及50%噪聲測點數據，經過10次蒙特卡洛仿真，PINN訓練過程中的損失函數變化如圖8所示。可以看出，隨著數據中噪聲的增大，損失函數的收斂值也變大。

圖8 給定不同噪聲數據時損失函數隨迭代步數變化Fig. 8 Loss functions versus the iteration step for training data with different noisy levels

隨后，將PINN方法與共軛梯度法（CGM）的辨識結果進行對比，如表3所示。從表3可以看出，基于PINN的辨識結果與CGM結果精度相當。對于50%的噪聲數據，基于PINN的辨識結果有1.253%的相對誤差，這說明利用PINN的參數辨識方法對噪聲數據具有較強的魯棒性。這是因為，第一，損失函數中包含初始、邊界條件以及偏微分方程殘差構成的正則化項，可以防止網絡對噪聲過擬合；第二，由于采用了自動微分，一方面在計算溫度的偏導數時，不會引入標簽溫度的噪聲，另一方面，也可以避免傳統數值方法中低精度差分格式對辨識結果的影響[17]。在計算過程中注意到，網絡的初始權重對結果影響較大，當選取不恰當時，會導致訓練發散。

表3 導熱系數k辨識結果Table 3 Identification results of the thermal conductivity k

下面研究小導熱系數以及邊界條件對基于PINN的辨識方法的影響。

已知導熱系數k= 0.01時x= 0.1 m處的溫度場測點數據，以及邊界條件（11），利用PINN方法的導熱系數辨識結果如表4所示。從表4可以看出，在給定較小導熱系數的測點數據時，基于PINN方法可以對參數k進行高精度的辨識。

表4 導熱系數k辨識結果Table 4 Identification results of the thermal conductivity k

當給定第二類邊界條件（6）時，將導熱系數k= 1時x= 1 m處的Nt個溫度場測點數據作為樣本值，定義如下損失函數：

基于PINN方法的導熱系數辨識結果如表5所示。從表5可以看出，當熱傳導方程采用第二類邊界條件時，基于PINN方法的導熱系數的辨識精度與采用混合邊界條件熱傳導方程的辨識精度相當。

表5 導熱系數k辨識結果Table 5 Identification results of the thermal conductivity k

3 結論

針對一維熱傳導方程的正問題和逆問題，本文介紹了一種內嵌物理機理神經網絡的求解方法。該方法將神經網絡作為解的近似，并在給定數據和偏微分方程的約束下，對網絡進行訓練。通過將偏微分方程作為正則化部分加入到損失函數中，使得上述深度網絡具有較強的泛化性能以及一定程度上的可解釋性。

通過分析熱傳導方程求解的收斂誤差得到，采樣點較少時，采樣誤差是收斂誤差的主要組成部分，而當采樣點較多時，主導誤差變為優化誤差。再通過對方程中的未知參數辨識的魯棒性進行分析得知，由于損失函數中包含了初始條件、邊界條件以及方程殘差構成的正則項，因此，可以防止網絡對噪聲出現過擬合，從而增強辨識方法的魯棒性。但是，基于PINN的導熱系數辨識方法在給定不恰當的超參數時，可能會導致訓練發散。