廣義乘子法求解構造變分問題的神經網絡方法

歐陽曄，江 巍，吳 怡，馮 強，鄭 宏,3

(1.三峽庫區地質災害教育部重點實驗室，湖北，宜昌 443002；2.三峽大學土木與建筑學院，湖北，宜昌 443002；3.北京工業大學建筑工程學院，北京 100124)

在工程領域，力學問題的數學模型往往都是以微分方程的形式提出的，微分方程與初邊值條件一起構成了微分方程的定解問題[1]。求解偏微分方程的途徑可分為解析法與數值解法兩大類。解析法適用于形狀規則的簡單問題，以逆解法、半逆解法為代表。實際問題的邊界條件往往比較復雜，解析法求解困難，于是發展出了以有限差分法、有限單元法等為代表的數值解法，這些算法的求解效果嚴重依賴事先剖分的網格。此外，有限差分法在不同的方程上需使用不同的差分格式；有限單元法針對不同的問題需選取不同的單元類型，這些局限性大大影響了其應用范圍[2]。

神經網絡算法的飛速發展[3]和Tensorflow[4]、Pytorch[4-5]等深度神經網絡框架的發布，使得近年來神經網絡的應用研究呈現井噴式發展，如葉繼紅和楊振宇[6]采用生成式對抗神經網絡生成的風場數據與目標數據之間具有相同的特性；許澤坤和陳雋[7]利用神經網絡準確預測了地震作用下框架結構的非線性結構響應；程詩焱等[8]基于BP 神經網絡提出的地震易損性曲面分析方法，利用該方法可以得到可信的損傷概率結果；趙林鑫等[9]將比例邊界有限元法與神經網絡相結合，建立適用于薄板結構裂紋狀缺陷識別的反演模型；鄭秋怡等[10]建立了基于長短時記憶神經網絡的大跨拱橋溫度-位移相關模型，可準確描述多元溫度與位移的非線性映射關系。為克服傳統數值解法求解偏微分方程時的局限性，國內外學者嘗試將神經網絡引入偏微分方程的求解，深度Ritz 法[11]、深度伽遼金法[12]、深度能量法[13]等相繼被提出，并取得了初步成功應用。RAISSI 等[14-15]提出的物理信息神經網絡(Physics Informed Neural Networks,PINN)，將原問題的物理信息嵌入到神經網絡中，通過訓練神經網絡獲得近似解；LU 等[16]提出了一種改善PINN訓練效率的自適應細化算法，并開發了DeepXDE 求解器；郭宏偉和莊曉瑩[17]提出了一種使用Adam 與LBFGS 兩種優化器訓練神經網絡的策略，克服了傳統神經網絡訓練的結果易陷入局部最優解的問題；黃鐘民等[18]發展了一種耦合神經網絡方法用于求解面內變剛度薄板的位移與撓度；唐明健和唐和生[19]使用神經網絡求解了矩形薄板的力學正反問題；HE 等[20]從逼近論的角度證明了使用ReLU 激活函數的神經網絡與有限元是等價的。

上述研究中，邊界條件的施加與傳統有限元略有不同。現有研究中神經網絡的損失函數是一個包含控制方程與邊界條件的泛函，通常使用經典罰函數法進行構造，將原約束優化問題轉化為無約束優化問題，如郭宏偉等[17]、黃鐘民等[18,21]。理論上罰函數法中罰因子的數值應為一個無窮大的數[22]，在實際計算中，罰因子取得過大會造成方程的病態導致無法求解，取得過小則無法起到懲罰的作用，其值不易把控。以L1精確罰函數法[23]為代表的改進方法，可有效克服經典罰函數法中罰因子需趨近于無窮的局限，并在非線性規劃問題的神經網絡求解中取得成功應用[24-25]，罰因子的具體取值仍依賴于人為確定。

Lagrange 乘子法[26]也可將約束優化問題轉化為無約束優化問題，但采用Lagrange 乘子法構建的神經網絡同樣存在一定不足。當神經網絡損失函數對應的Lagrange 函數在平衡點處其Hessian 矩陣為非正定矩陣時，該平衡點將偏離原問題的最優解[24]。理論上，廣義乘子法可克服Lagrange 乘子法的這一局限[22]。鑒于此，針對邊界條件復雜的偏微分方程組，本文提出一種采用廣義乘子法施加邊界條件的神經網絡方法。該方法首先通過神經網絡獲得預測解，再使用廣義乘子法構建損失函數并計算損失值，最后利用梯度下降法進行參數尋優，判斷損失值是否滿足要求；若不滿足，則更新罰因子與乘子后再進行求解，直至損失滿足要求。為驗證新方法的求解精度和計算效率，本文求解了Kirchhoff 薄板和簡支梁問題，并與使用經典罰函數法、L1精確罰函數法、Lagrange乘子法施加邊界條件的神經網絡求解進行比較。此外，采用不同拓撲結構的神經網絡和域內配點求解Kirchhoff 薄板問題，探討了神經網絡拓撲結構與域內配點數對新方法求解精度的影響。

1 微分方程及其泛函

工程技術領域，許多自然規律是以偏微分方程(組)的形式提出的，常見的偏微分方程分為三大類：橢圓型方程、拋物型方程、雙曲型方程，上述方程的共同點在于方程內部包含待求解的未知場量及其導數，即：

式中：A為偏微分方程組；u為物理量；Ω 為問題求解區域。靜態場問題邊值條件通常分為三類：第一類條件為指定未知場量在邊界上的值；第二類條件為指定未知場量在邊界外法線方向上的值；第三類條件為上述兩類邊界條件的線性組合，即：

式中：B為邊界條件組；Γ 為求解區域的邊界。

式(1)所示的控制方程與式(2)所示的邊值條件一起構成了微分方程的定解問題。通常來說，微分方程對解的要求較高，若直接對微分方程進行求解較為困難，因此常使用加權余量法獲得微分方程的等效積分弱形式，弱形式會降低對解的要求。對滿足線性、自伴隨的微分方程而言，采用伽遼金加權余量法所獲得的等效積分弱形式與使用自然變分原理獲得的泛函是相同的，弱形式將具備明確的物理意義，例如，彈性力學中平衡微分方程的等效積分弱形式就是虛功方程。對不滿足線性、自伴隨的微分方程而言，需使用構造變分原理來獲得泛函，使用罰函數法構造的泛函表達式如下：

式中：v為一組與A(u)方程個數相同的任意函數；M為罰因子，是一個標量。

2 神經網絡方法

2.1 神經網絡基本原理

神經網絡是一個由權重w、偏置b、激活函數σ 等定義的非線性函數，神經網絡的輸入x，輸出為u(x)，其前向傳播可表達為：

式中，下標l表示神經網絡的層數，Fl(l=1, 2, 3, … ,l)為第l層的映射規則，定義為：

初始狀態下，神經網絡中參數θ={w;b}是隨機生成的，不能滿足實際要求，需進行參數更新，現階段常使用BP 算法進行參數更新。

萬能近似定理[20]指出，神經網絡具有逼近任意復雜非線性函數的能力。神經網絡具有強大非線性函數擬合能力的關鍵在于其激活函數是非線性的。現有研究中常使用ReLU 函數和Tanh 函數兩種激活函數。ReLU 函數[27]的表達式為：

該函數及其導數圖形如圖1(a)和圖1(c)所示。ReLU 函數形式簡單、計算成本低，可避免出現梯度消失現象。Tanh 函數具有C∞性質，其表達式為：

圖1 兩種常用的神經網絡激活函數及其導數Fig.1 Two widely used activation function in neural network and their derivatives

當式(3)涉及未知場量的高階導數時，Tanh函數更為適用，因此在目前神經網絡求解力學問題中廣泛使用[17-18]，該函數及其導數圖形如圖1(b)、圖1(d)和圖1(e)所示，本文采用Tanh 函數為激活函數。

2.2 深度配點法

深度配點法是由RAISSI 等[14-15]提出的一種由物理信息驅動的神經網絡模型。根據萬能近似定理，偏微分方程組的解u可使用神經網絡來逼近。對于深度配點法，需首先在求解域內布置一系列的散點xΩ，另外邊界上也需布置一系列的散點xΓ。圖2 給出了一矩形區域的配點示意，點★代表xΩ，點●代表xΓ，散點x={xΩ;xΓ}構成了神經網絡的訓練數據。

圖2 配點示意圖Fig.2 Schematic diagram of collocation points

從偏微分方程出發構建如圖3 所示的物理信息網絡。左邊的神經網絡作為偏微分方程的近似解，輸入為散點坐標，輸出 ^u(x;θ)由式(4)計算。右邊的物理信息則是將近似解^u帶入控制方程與邊界條件中所產生的余量。

圖3 物理信息網絡示意圖Fig.3 Diagram of Physics Informed Neutral Network

深度配點法一般采用經典罰函數法施加邊界條件，建立的神經網絡損失函數表達式為：

式中：NA、NB分別為微分方程與邊界條件方程的個數；PΩ、PΓi分別為域內、相應邊界上的配點數；M為罰因子；其中A(x;θ)與B(x;θ)共享同一套超參數θ，并且都是坐標的函數。當loss 取得最小值時，^u(x;θ)就是偏微分方程的近似解。

將經典罰函數法替換為L1精確罰函數法，則神經網絡損失函數的表達式變化為：

式中，|·|為絕對值。

如果使用Lagrange 乘子法施加邊界條件，建立的神經網絡損失函數為：

式中，ki為Lagrange 乘子。

2.3 采用廣義乘子法施加邊界條件的深度配點法

如引言中所述，使用經典罰函數法、L1精確罰函數法、Lagrange 乘子法施加邊界條件的神經網絡在求解時均存在一定的局限性。廣義乘子法[20]結合了罰函數與Lagrange 乘子法的優點，使得罰因子在適當大的情況下也能求得原約束問題的解。因此，本文將廣義乘子法引入神經網絡之中，提出一種基于廣義乘子法施加邊界條件的神經網絡方法。該方法左邊的神經網絡仍由式(4)計算，右邊的物理信息則基于廣義乘子法構建，最終建立的神經網絡損失函數表達式為：

式中：M為罰因子；ki為相應邊界方程對應的乘子，當loss 取得最小值時，^u(x;θ)就是偏微分方程的近似解。使用梯度下降法訓練參數θ，從而使得loss 取得最小值，具體做法如下：

S1 在求解域內與邊界上布置一系列的配點；

S2 利用廣義Lagrange 乘子法將約束問題轉化無約束的待優化函數，得到式(11)；

S3 構建前饋神經網絡，并將其作為試函數；

S4 給定初始值θ0，初始乘子初始罰因子M0，放大系數α>1，誤差閾值ε，參數γ∈(0,1)。這里，各量的上標代表迭代次數；

S5 問題求解：以θn-1為起始點，利用梯度下降法求解當前乘子、罰因子下的最優解θn，具體做法為：

① 計算損失函數loss 的值，判斷損失函數的值是否達到精度要求；若達到要求，則進入S6，若未達到要求，則進入②。

② 將損失標量loss 對神經網絡參數θ 求偏導，并更新神經網絡參數，其更新公式為：

式中：Δθn-1為θn-1的更新量；η 為神經網絡學習率。參數更新完畢后進入①。

S6 判斷迭代終止條件：根據下式計算邊界上的誤差h(θn)，若h(θn)<ε 成立，則停止迭代并輸出θn；否則，進入步驟S7；

否則：

S8 更新乘子項：使用如下公式更新乘子項及迭代次數，并返回步驟S5：

3 算例分析

本文求解程序均采用Python 的Tensorflow2.6版本進行編寫，運行平臺為Windows10，硬件配置為CPU：Intel Xeon E-2224 @ 3.40 GHz，GPU：NIVIDIA Quadro RTX 4000。

3.1 周邊簡支矩形Kirchhoff 薄板

如圖4 所示的周邊簡支矩形Kirchhoff 薄板，邊長a=b=1 m，薄板厚度h=0.01 m，彈性模量E=10 MPa，泊松比ν=0.2，板中面的均布荷載為q=10 N/m2，該問題的重三角級數解[28]為：

圖4 Kirchhoff 矩形薄板示意圖Fig.4 Diagram of a Kirchhoff thin rectangular plate

該問題的微分控制方程為：

其邊界條件為：

式中：Γ={Γ1，Γ2，Γ3，Γ4}為邊界，分別代表{x=0，x=a，y=0，y=a}四條邊；Mn為彎矩。

根據廣義乘子法構建的損失函數為：

式中：PΩ、PΓi分別為域內、相應邊界上的配點數；ki為相應邊界方程對應的乘子。

初始設定罰因子M=1，乘子ki=1(i=1, 2, …, 8)，放大系數α=2，誤差閾值ε=1×10-3，參數γ=0.01，采用本文提出方法進行求解。按照圖2 的方式在求解域與邊界上進行配點，其中域內配點20×20 個，邊界上每邊配點40 個。神經網絡的拓撲結構設置為2-1*30-1，表示輸入層神經元數量為2，分別為配點的橫縱坐標；隱藏層層數為1，隱藏層神經元數量設置為30 個；輸出層神經元數量為1，為配點的撓度預測解，優化器選取Adam，學習率設置為0.01。在神經網絡訓練過程中，若損失值連續500 輪沒下降則自動更新罰因子、乘子后再進行計算，且在每輪罰因子、乘子下的訓練不超過10000 輪。

經典罰函數法、L1精確罰函數法、Lagrange乘子法施加邊界條件建立的神經網絡求解該算例時，配點方式、神經網絡拓撲結構、激活函數、優化器和學習率等設定均與本文算法求解時的設定保持一致，訓練完畢后取損失值最小的一輪作為輸出結果。為探討罰因子取值對結果的影響，經典罰函數法和L1精確罰函數法的罰因子分別取M=1、10、50、100、500 和1000。Lagrange 乘子法的乘子初始值設定與本文算法求解時的設定相同。

以配點處的神經網絡預測解與重三角級數解[28]的絕對誤差平均值作為參考，對上述神經網絡求解結果的精度進行比較。如圖5 所示，使用經典罰函數法時，計算誤差隨著罰因子的增加呈先減小后增加的趨勢，當M=50 時其數值精度最高，絕對誤差平均值為5×10-4m。這表明在使用基于經典罰函數法構建的神經網絡時，盲目的增大罰因子其效果可能適得其反。當罰因子M=1 時，L1精確罰函數法具有最好的求解精度，此時絕對誤差平均值為5.2×10-4m。隨著罰因子增大，其誤差反而快速增加至與Lagrange 乘子法同數量級。對于本算例而言，Lagrange 乘子法的求解精度最低，絕對誤差平均值達到1.4×10-2m，其原因正如引言中所提及的，Lagrange 乘子法建立的神經網絡損失函數對應的Lagrange 函數在平衡點處的Hessian矩陣可能為非正定矩陣。本文建立的廣義乘子法經過5 次參數更新之后，罰因子M=32 時絕對誤差平均值減小至4.6×10-4m，5 次參數更新的詳細情況參見表1。

表1 Kirchhoff 薄板算例罰因子、乘子更新值Table 1 The updating of penalty factor and multipliers in Kirchhoff thin plate example

圖5 四種神經網絡模型撓度結果的絕對誤差平均值Fig.5 Mean value of the absolute deflection errors resulted by four neural network models

以矩形簡支薄板邊界處的撓度計算結果為參考，進一步對經典罰函數法、L1精確罰函數法和本文算法進行比較。邊界處撓度理論上為零，經典罰函數法罰因子M=50 時，撓度絕對值最小為6.1×10-4m；L1精確罰函數法罰因子M=1 時，撓度絕對值最小為4.2×10-4m；本文算法經過5 次參數更新后，撓度絕對值最小達到3.2×10-4m。本文算法求解出的薄板撓度云圖如圖6 所示，經典罰函數法、L1精確罰函數法和Lagrange 乘子法均可以獲得與圖6 接近的撓度分布，結果的主要區別在于計算精度。

圖6 本文算法求解的薄板撓度云圖Fig.6 Deflect distribution of thin plate resulted by the proposed algorithm

基于上述比較分析，可以發現本文方法的求解精度整體上要優于其他方法。經典罰函數法和L1精確罰函數法，當罰因子取值合適時，求解精度與本文方法具有一定的可比性；若罰因子取值不合理，則計算結果誤差較大。由于避開了罰因子的取值問題，因此本文方法的實用性強于罰函數法。

3.2 均布荷載作用下的簡支梁

設有一受均布荷載作用的簡支梁，如圖7 所示，梁長l=1 m、梁寬b=0.01 m，梁高h=0.01 m，均布荷載q=10 N/m2，彈性模量E=1.0 GPa，求梁的撓度方程w(x)。

圖7 平面簡支梁示意圖Fig.7 Diagram of a simple supported beam

該問題的控制方程為：

邊界條件為：

使用廣義乘子法構造的損失函數為：

使用本文涉及的四種不同神經網絡進行求解時，神經網絡拓撲結構均設定為1-1*10-1，優化器為Adam，學習率設定為0.01，域內配點數為10，邊界配點數為2。基于3.1 節中算例求解的經驗，使用罰函數時應避免罰因子取值過大，本算例執行時，經典罰函數法罰因子采用50 和100 兩種取值，L1精確罰函數法的罰因子采用1 和50 兩種取值。廣義乘子法求解時初始設定為：罰因子M=5，乘子k1=k2=k3=k4=1，放大系數α=5，誤差閾值ε=1×10-5，參數γ=0.25。Lagrange 乘子法的乘子初始值設定與廣義乘子法保持一致。

四種神經網絡求解得到的撓度曲線如圖8 所示。經典罰函數法、Lagrange 乘子法和廣義乘子法計算得到的撓度結果與解析解近乎完全重合，其絕對誤差量級為10-4m。L1精確罰函數法計算所得的撓度結果誤差則相對較大，其絕對誤差量級為10-2m，這進一步證明了罰因子取值對基于罰函數法構建的神經網絡求解結果的影響。

圖8 采用不同神經網絡計算的平面彎曲梁撓度曲線Fig.8 Deflection curve of plane bending beam by different neural network models

圖9 給出了不同神經網絡求解時的誤差-耗時曲線。可以發現，由于數學理論的不同，四種神經網絡求解過程呈現顯著的差異。罰函數法的單次訓練過程耗時較短，在向最終解逼近的過程中呈現反復的來回震蕩現象，L1精確罰函數法的震蕩幅度略小于經典罰函數法，但在本算例中由于罰因子取值不太合理所以誤差無法繼續降低。Lagrange 乘子法初期誤差下降較快且無明顯震蕩，當誤差下降至10-2m 量級后，其向最終解逼近時呈現逐步向下的來回震蕩。本文采用的廣義乘子法，經過4 次迭代后誤差達到收斂標準。第一次迭代時耗時較長，原因在于初始值距離最終解較遠，計算參數可能不合理。隨迭代輪數增加，其計算參數將根據最終解的預測位置實時調整，后續迭代明顯耗時減少。在向最終解逼近過程中，誤差雖然也存在震蕩現象，但每次震蕩后誤差均會大幅降低。比較四種神經網絡求解該問題的耗時可發現，本文提出方法總耗時最少，與Lagrange 乘子法相比甚至呈現數量級的區別，因此在計算效率方面本文提出的方法明顯優于其它的神經網絡方法。

圖9 采用不同神經網絡計算的誤差-耗時曲線Fig.9 Error and time-consuming curves by different neural network models

此算例求解時，廣義乘子法神經網絡求解4 次迭代過程中的參數更新情況如表2 所示。

表2 簡支梁算例罰因子、乘子更新值Table 2 The updating of penalty factor and multipliers in simple supported beam example

4 廣義乘子法神經網絡求解的影響因素分析

4.1 神經網絡拓撲結構對求解過程的影響

針對3.1 節中的Kirchhoff 薄板問題，構造隱藏層層數為1 層、2 層和3 層，每層神經元數量分別為10 個、20 個和30 個的神經網絡，對其進行求解，觀察神經網絡拓撲結構對求解過程的影響。為簡單起見將拓撲結構表示為2-L*N-1，其中，L為隱藏層層數，N為每層隱藏層的神經元數量，其他參數設置與3.1 節中保存一致。使用不同拓撲結構的神經網絡進行求解，最終得到的誤差耗時曲線如圖10 所示。

圖10 不同結構神經網絡求解的誤差-耗時曲線Fig.10 Error and time-consuming curves by different neural network structures

可以發現，隱藏層層數為1 層時，神經網絡結果的絕對誤差平均值均可收斂至10-3m 左右，增加隱藏層神經元數量可顯著降低神經網絡結果的誤差。隱藏層層數為2 層時，隱藏層神經元數量增加對于降低神經網絡結果的誤差有一定作用，但是誤差曲線在收斂過程中可能會產生較大振蕩。隱藏層層數為3 層時，神經網絡結果始終存在一定誤差，而且該誤差不會隨著神經元數量的增加而減小。增加隱藏層層數與增加神經元數量都會導致求解耗時增加，但相對來說，隱藏層層數對計算效率影響更大。基于上述比較，可認為本文建立的廣義乘子法神經網絡求解偏微分方程時，過多的隱藏層層數反而會導致計算效率和求解精度降低。

4.2 配點數與隱藏層神經元數量

理論上，神經網絡的計算精度與訓練集的大小和隱藏層神經元數量密切相關。基于前面的分析，使用2-1*N-1 的神經網絡拓撲結構，其中N=5、10、15、20、30、50，其他參數設置與3.1 節中保持一致，選取5×5、10×10、15×15、20×20、30×30 和50×50 等六種不同的配點方式，探討不同配點數量在不同隱藏層神經元數量下的收斂情況，并將結果繪制于圖11。

圖11 不同配點數與隱藏層神經元數量組合求解的誤差Fig.11 Error resulted by different combinations of collocation points and neurons in the hidden layer

可以發現，配點數為5×5 時，不同隱藏層神經元數量的神經網絡求解結果的絕對誤差平均值均較大；使用其他配點方式，當隱藏層神經元數量超過20 個時，神經網絡求解結果的誤差均能降低至10-3m 以下。因此，使用本文建立的廣義乘子法神經網絡求解偏微分方程時，配點數量不宜過少，且隱藏層神經元數量應達到一定標準。

5 結論

針對以往物理信息神經網絡采用經典罰函數法進行求解時罰因子取值難以確定的問題，本文結合廣義乘子法提出一種改進方法，經算例驗證，可得出如下結論：

(1)與使用經典罰函數法、L1精確罰函數法和Lagrange 乘子法施加邊界條件構造的神經網絡相比，采用廣義乘子法施加邊界條件構造的神經網絡具有更好的數值精度和更高的求解效率，而且求解過程更加穩定。

(2)使用本文提出的方法求解偏微分方程組時，神經網絡隱藏層層數不宜過多。隱藏層層數設置合理的前提下，隱藏層神經元數量和配點數量均應達到一定標準。