基于深度強化學習算法的顆粒材料應力-應變關系數據驅動模擬研究1)

狄少丞 馮云田 瞿同明 于海龍

* (哈爾濱工程大學船舶工程學院,哈爾濱 150001)

? (斯旺西大學辛克維奇工程計算中心,英國斯旺西,SA1 8EP)

引言

在研究泥沙、巖土、粉體等顆粒材料宏觀力學行為時,通常在連續介質力學框架下建立其唯象本構模型,通過引入內變量來描述顆粒材料的路徑相關、歷史相關等特性[1-5].雖然基于連續假設的唯象本構模型在描述顆粒材料宏觀力學特性時具備很好的工程適用能力,但是顆粒材料由于其自身的離散特性,具有摩擦性,剪脹性,壓硬性,各向異性,荷載路徑和歷史相關性等典型特征,并且在外載作用下還會呈現出應變局部化等復雜的演化機制,連續模型要體現出顆粒材料精細的宏觀特性就不得不引入繁雜的唯象假設和眾多待標定的自由參數.

離散元方法[6](DEM)可跟蹤并計算每個散體單元的運動信息及單元間的接觸作用力,在顆粒材料力學行為研究中得到了廣泛的應用.隨著計算機硬件計算能力的提升,以及顆粒幾何特征的精準刻畫[7-9]和高精度接觸算法的提出[10-12],使得離散元方法的應用領域和計算精度得到大幅提升.然而離散元方法用于直接模擬較大空間尺度和時間尺度的工程問題仍需要巨大的計算成本.為此,人們提出了多尺度計算方法以提升顆粒材料問題的計算效率[13-14].在此框架中,通常基于均勻化思想采用低尺度的離散元計算在宏觀模擬中替代傳統的唯象本構關系.均勻化即從一定體積的顆粒聚集體內部顆粒間作用的微觀物理量獲得在該體積內連續性力學描述的過程[15].而均勻化的首要任務是創建代表性體積單元(RVE)[16],RVE 是包含若干顆粒的顆粒聚集體,是連接顆粒介質微結構尺度和宏觀工程尺度的橋梁.顆粒介質多尺度方法常采用DEM 與FEM 相結合的方式來實現,即在FEM 計算中提供施加給RVE 的變形條件,RVE 顆粒試樣在變形條件下產生與應變對應的應力值,隨后將應力-應變關系或者切向剛度矩陣參數映射到FEM 網格中高斯積分點上進行下一時步的求解[17-20].采用多尺度方法雖然避免了對顆粒介質全尺度的DEM 求解,但需要在FEM 網格每個高斯積分點上構建單獨的RVE[21-22],并且需要在計算中維護與FEM 網絡中高斯點數量相當的RVE 加載狀態信息,這又增加了計算機的存儲與計算開銷[19],對于跨多個尺度的問題計算量會更大.

為此,人們提出了“離線訓練(offline training)、在線模擬(online prediction)”的計算策略[23-24]來提升計算效率.相比于直接采用RVE 模擬結果的“在線”計算來說,“離線”策略下的RVE 模擬只產生用來訓練“代理模型”的數據集[25-26],即采用機器學習算法訓練并獲得RVE 的應力-應變關系,從而宏觀尺度FEM 的計算不再需要與細觀尺度RVE 的求解交替進行,在FEM 計算階段直接讀取RVE 的應力-應變關系,這就大大提高了計算效率[27].

可采用不同的神經網絡模型對顆粒介質的本構關系進行訓練與計算,常用的基礎模型有經典BP 神經網絡模型[28-29]和循環神經網絡模型[30],二者均可將應力、應變值分別作為輸入、輸出量進行訓練.循環神經網絡中的長短期記憶網絡(LSTM)和門控循環單元網絡(GRU)具有記憶功能,對于與加載路徑和加載歷史相關的顆粒材料而言,在應力-應變關系建立方面表現出了較強的適用性和較高的擬合精度[26].

在常規循環神經網絡模型基礎之上,也有一些學者結合其他理論模型進一步研究了顆粒材料力學行為.Wang 等[31-33]結合博弈論和強化學習模型研究了顆粒材料的界面位移與張力之間的力學行為以及本構關系.瞿同明等[34]結合細觀力學理論和循環神經網絡模型訓練了顆粒材料的應力-應變關系.Karapiperis 等[35]在有限元中采用了基于離散元數據訓練得到的循環神經網絡本構模型.Qu 等[36]基于坐標不變性原理討論了從主應力-應變空間到一般應力-應變空間的數據驅動本構模型訓練策略.

綜上所述,建立精度高、計算快的RVE 應力-應變關系對于顆粒介質的多尺度模擬是至關重要的.本文借助有向圖方法和AlphaGo Zero 深度強化學習算法,在顆粒物質力學的研究基礎之上,試圖建立一種自動尋找最優內變量間映射關系、完全依靠“數據驅動”的RVE 力學行為建模方法.

1 顆粒材料的本構關系

1.1 與顆粒材料應力-應變關系相關的重要內變量

離散元方法可有效地在細觀尺度下模擬顆粒材料的宏觀力學行為,本文將基于球形離散單元法線性接觸模型模擬顆粒材料的應力-應變關系.從細觀力學的角度,顆粒材料的宏觀應力-應變增量關系可表示為

其中,彈性剛度矩陣Cijmn在應變均質化的假設下可以表示為[37]

其中,kn和ks為顆粒之間的法向和切向接觸剛度;V為顆粒試樣的體積;δin為克羅內克函數;Lk為每個接觸對內兩個顆粒的中心距離;為接觸k的方向向量在i方向上的分量;NC為接觸對的數量,令

則式(2)可表示為

彈性剛度矩陣Cijmn可由顆粒間的法向接觸剛度kn和切向剛度ks,以及單位體積內的細觀結構張量表示.當顆粒材料承受外載荷作用時,顆粒間發生相對滑移,導致組構張量在外載作用下會不斷發生變化.

式(1)中,當切應力和切應變為0 時,采用增量形式表示的主應力和主應變關系為

其中,C1122=C2211,C1133=C3311,C2233=C3322,彈性剛度矩陣中共有6 個獨立分量表示主應力和主應變之間的關系.

除了上述推導所涉及的彈性剛度矩陣外,顆粒物質力學研究表明顆粒材料的孔隙率φ(對于單相顆粒材料,孔隙率定義為一定體積內非顆粒組分所占據的體積分數)和組構參數對于其宏觀力學行為具有重要的影響.其中,巖土材料作為一種典型顆粒材料,現代土力學理論的一個突破性進步就在于將孔隙率等狀態參數作為硬化法則的一部分,將其考慮到巖土材料的本構關系研究中[38-39].近年來巖土本構研究領域另一個重要進步在于認識到顆粒材料的組構對于其本構行為的重要影響[40],考慮組構參數的演化過程是近年來巖土本構理論的重要發展趨勢[41-43].通常顆粒材料的細觀組構有多種張量描述方式,本文采用式(3)和式(4)中的4 階組構張量來共同描述顆粒材料的組構演化過程.

1.2 基于深度強化學習的本構研究思路

顆粒材料的本構研究從不同的視角發展了不同的理論模型和方法,這些理論方法對問題本身的研究都取得了新的進展.如何將這些來自不同視角的重要發現整合到一個統一的框架內是一個思路合理但是實踐困難的挑戰.本文借助深度強化學習的思路,將從顆粒物質力學研究中得到的彈性剛度矩陣(式(2),用Cf表示),孔隙率參數φ和組構張量(式(3),用Af表示)結合在一起,試圖在這些重要變量的數據樣本基礎上,借助人工智能算法整合已有的知識資源,發掘新的認識,自動發現能描述顆粒材料應力應變行為的本構關系.

2 本構關系的有向圖表示

有向圖來源于數學上的圖論,主要指采用節點來表示物理量,采用邊線來描述物理量之間相互關系的一種方法.包含內變量的本構關系將通過有向圖來表達[31],有向圖為信息從本構關系輸入值開始,流經內變量,最終流向輸出值的信息流動網絡.

2.1 有向圖網絡

圖1 為以內變量和輸入輸出值構建的變量間信息流動的有向圖,圖中的箭頭代表了變量間的“源-目標”數據流動方向.有向圖中的變量稱為節點,變量間的鏈接稱為邊,對于每個鏈接而言,箭頭指出的節點稱為源節點,箭頭指向的節點稱為目標節點,整個有向圖的輸入稱為根節點,輸出稱為葉節點.

圖1 有向圖中的信息流動Fig.1 Information flow in a directed graph

一個有向圖代表一種本構關系,為保證本構關系中變量間信息流動的合理性,有向圖的構建需遵循以下規則:

(1)有向圖中僅存在唯一的葉節點(該節點不作為其他任何節點的源節點,僅作為目標節點);

(2)有向圖中僅存在唯一的根節點(該節點不作為其他任何節點的目標節點,僅作為源節點);

(3)有向圖中可存在孤立節點,該節點不參與到有向圖的建立中,即本構關系中不考慮該節點對應內變量對本構關系建立的貢獻;

(4)有向圖中不存在閉環信息流(即信息流從某一結點開始,最后又流回該節點);

(5)如果某一節點鏈接到了根節點或葉節點,那么此節點需包含在從根節點到葉節點的完整路徑中.

2.2 信息流動路徑

對于簡單的有向圖鏈接,如圖2(a)所示,信息流動方向可直接表示為ε→φ→Af→Cf→σ,而對于復雜的有向圖鏈接,如圖2(b)所示,需要建立合理的信息流動路徑.圖2(b)中的信息流動路徑需要從有向圖輸出變量σ開始以遞歸的方式向前尋找信息“源”,每一節點只尋找直接指向該節點的“源”節點,形成“一對一”或“多對一”的信息映射關系.圖2(b)中信息流動路徑對可表示為:{ε→Cf},{ε,Cf→Af},{ε,Af→φ}及{ε,φ,Cf→σ},即信息從ε流入有向圖后計算出Cf,隨后ε和Cf共同計算出Af,接著由ε和Af共同計算出φ,最后由ε,φ,Cf共同計算出σ,實現了ε到σ的信息流動.將采用循環神經網絡擬合出每對路徑中變量間的映射關系,用于后續“數據”本構模型的建立.

圖2 有向圖表示的兩種本構關系Fig.2 Two constitutive laws represented by two directed graphs

2.3 本構模型打分

對有向圖建立本構關系的預測精度進行打分,以評估通過數據建立的本構關系對于訓練數據的泛化能力.通過比較信息流路徑計算出的應力值與訓練數據中應力值的均方誤差(mean squared error,MSE),評價強化學習算法中不同有向圖代表的本構模型的優劣.

對于具有Ndata個數據樣本的樣本集,每個樣本的均方誤差可表示為

對集合{MSEi}內各均方誤差值進行升序排序,并取集合中第P%的MSE值(記為 εP%)計算本構模型打分值S

式中 εcrit?1,是需提前設定的MSE 誤差標準,以評估生成本構模型的精確度,其中P%取值90%,即集合中第90%的MSE值,以及 εcrit=1×10-6.

3 AlphaGo Zero 深度強化學習算法

本節將采用AlphaGo Zero 中的強化學習算法建立基于有向圖的本構關系.AlphaGo 是Google DeepMind 團隊開發的一個基于深度神經網絡的圍棋人工智能程序.2015 年10 月AlphaGo 擊敗了歐洲冠軍Fan Hui,2016 年3 月擊敗韓國棋手Lee Sedol.2017 年10 月,DeepMind 公布了最新的Alpha-Go Zero,在與AlphaGo 對戰中取得了100∶0 的戰績.與AlphaGo 不同的是,AlphaGo Zero 未使用任何圍棋領域的專業知識和人為干預,訓練數據全部來自于自我對弈(self-play),算法方面,AlphaGo Zero只采用了一個神經網絡,將棋盤狀態作為輸入特征,來指導蒙特卡洛樹搜索(Monte Carlo tree search,MCTS)探索最優的落子動作.

本構模型建立過程中,選定輸入變量、輸出變量和內變量后,可以通過“窮舉法”列出所有滿足信息流動規則的有向圖,然后對每一種有向圖打分,得分最高者即為最優本構關系.對于內變量較少的本構關系建立過程,可考慮采用“窮舉法”尋找出最優的本構關系,但如果引入的內變量的數量較多,采用直接尋找的方法將極大地增加計算成本,因此建立一種能夠自動尋找最優本構關系的算法是必要的.本構關系建立過程與AlphaGo Zero 實現過程類似,同樣可采用馬爾可夫決策過程(Markov decision process,MDP)對有向圖中變量鏈接關系進行選擇,以建立得分最高的有向圖鏈接關系.下面將給出采用AlphaGo Zero 算法建立本構關系有向圖的流程.

3.1 本構關系建立過程

本構關系建立過程可表述為馬爾可夫決策過程,本構關系的初始狀態如圖3(a)所示,本構關系建立的過程就是智能體在初始狀態下對有向圖中的邊進行激活,并獲得最大獎勵的過程.所有可能被激活的邊稱為“動作空間”,并遵循有向圖網絡建立準則,所有動作如圖3(b)所示,可描述為{ε→Af,ε→Cf,ε→φ,ε→σ,Af→Cf,Af→φ,Af→σ,Cf→Af,Cf→φ,Cf→σ,φ→Af,φ→Cf,ε→σ}.馬爾可夫決策過程中的“狀態s”為動作空間中有向圖邊的激活狀態;“動作a”為在當前狀態下準備激活的有向圖邊;“獎勵r”為根據本構關系打分情況確定的獎勵值,當打分高于給定的分值時取1,否則取-1.

圖3 MDP 中的初始狀態和可執行動作Fig.3 Initial state and all possible actions in MDP

在有向圖當前狀態st下,智能體執行動作at(激活某一條有向邊),當前狀態會轉移到t+1 時間步的新狀態st+1,其中動作at選自于st狀態下每一個有效動作概率的集合π(st),同時智能體獲得一個獎勵rt+1,該獎勵值的獲得是由于智能體在狀態st下執行了動作at,但是在實施過程中具體的獎勵值rt只有在本構關系完全建立后才能獲得.當本構關系得分高于給定分值時,本構建立過程中所有動作的獎勵rt≤T=1,反之,rt≤T=-1.圖4 為由馬爾可夫決策過程表示的一次完整的本構關系建立過程.

圖4 一次完整的有向圖本構關系建立過程Fig.4 A complete modeling process of generating a constitutive relationship

3.2 自主學習[44]

在AlphaGo Zero 算法中,主要通過訓練一個超參數為θ的策略/價值網絡fθ來指導本構關系有向圖建立過程中動作的選擇.策略/價值網絡將有向圖的鏈接狀態作為輸入,輸出一個向量p和標量v,即(p,v)=fθ(s),p表示在當前狀態下鏈接每一條有向邊的概率,v表示在當前狀態下本構關系獲得最高分的概率.

策略/價值網絡fθ通過自主學習來提升策略選擇的能力,使得有向圖鏈接得分最高.策略選擇能力的提升通過MCTS 算法來實現,而在每一狀態下蒙特卡洛樹擴展過程中的動作概率則是通過策略/價值網絡fθ進行計算,AlphaGo Zero 實現的自主學習過程如圖5 所示.自主學習過程中需運行多次MCTS算法使策略/價值網絡fθ變得更強,而MCTS 算法的實施均采用更新后的策略/價值網絡.

圖5 AlphaGo Zero 實現的自主學習過程Fig.5 Self-play reinforcement learning in AlphaGo Zero

蒙特卡洛搜索樹中的每個節點表示有向圖的鏈接狀態s,搜索樹中的每條邊(s,a)表示在狀態s下允許的動作a.搜索樹中每條邊維持一個統計量集合{N(s,a),W(s,a),Q(s,a),P(s,a)},其中N(s,a)為在樹搜索過程中節點被訪問的次數,W(s,a)為總動作價值,Q(s,a)為平均動作價值,P(s,a)為選擇邊的先驗概率.

(1)選擇(select):自主學習過程中從搜索樹的根節點s0開始,直至在時間步L遇到之前未被遍歷的葉節點sL,在t＜L時間步內的動作選擇則根據上置信邊界at=max[Q(st,a)+U(st,a)] 確定

其中,cpuct為MCTS 算法中控制探索程度的常數,在初始階段傾向于選擇具有較高先驗概率的動作,后續會逐步傾向于具有較高動作價值的動作;b的取值范圍為狀態s下所有可執行的動作.

(2)擴展和評價(expand and evaluate):對葉節點擴展時,葉節點狀態下的每條邊(sL,a)維持的統計量需初始化為{N(sL,a)=0,W(sL,a)=0,Q(sL,a)=0,P(sL,a)=pa},同時動作概率p與狀態價值v通過策略/價值網絡fθ預測得到.

(3)回溯(backup):對遍歷過的邊的統計量進行更新,其中訪問次數更新為N(st,at)=N(st,at)+1,動作價值更新為W(st,at)=W(st,at)+v,及Q(st,at)=

(4)執行動作(play):從根節點s0選擇動作的概率與節點被訪問次數相關,可表示為π(a|s0)=,其中τ為控制探索程度的系數.執行動作后該動作對應的子節點將更替為根節點,并且保留該子節點下的所有節點信息和相應的統計值,其它樹結構將被丟棄.

每一個時間步內,在每個有向圖當前狀態s下,每一次激活有向圖邊需進行若干次MCTS 模擬,最終MCTS 給出最優的動作策略π.MCTS 中策略π由 πt=αθi-1(st) 計算得到,采用了前一次訓練后的策略/價值網絡fθi-1.當建立起完整的本構有向圖后,可以獲得最終的獎勵值z∈{-1,1} .每一時間步內將狀態、動作策略以及獎勵存儲于數據集(st,πt,zt)中,作為策略/價值網絡的訓練數據,神經網絡的超參數θ通過梯度下降法使得損失函數的取值最小,其中損失函數定義為

損失函數綜合考慮了獎勵值的均方誤差和交叉熵損失.

表1 為采用AlphaGo Zero 強化學習算法建立應力-應變有向圖的計算流程,主要分為兩個階段:策略/價值網絡訓練階段和最優有向圖生成階段.訓練階段共執行Niter次策略/價值網絡的訓練,每次訓練執行Ncollect次完整的有向圖建立流程以收集所需的訓練數據,在鏈接每一條有向邊時執行NMCTS次蒙特卡洛樹搜索計算.有向圖生成階段,只需執行一次完整的有向圖建立流程,就可選定最終的應力-應變關系.

表1 利用強化學習算法建立最優應力-應變關系有向圖流程Table 1 Reinforcement learning of directed graph of the stressstrain relationship

4 建立“數據驅動”本構關系

4.1 數據制備

首先構建三軸數值試驗來生成顆粒材料的應力-應變關系數據,采用離散元方法建模的三軸試樣如圖6 所示,試樣包含4037 個球形顆粒,顆粒半徑在2 mm～ 4 mm 范圍內服從均勻分布.顆粒之間采用線彈性接觸模型,法向和切向接觸剛度分別為1 × 105N/m 和5 × 104N/m,顆粒密度為2600 kg/m3,局部阻尼系數為0.5.雖然顆粒之間的接觸作用是線彈性的,但試樣整體會由于細觀組構的不斷演化而發生復雜的塑性變形行為.

圖6 三軸數值試樣Fig.6 Numerical sample of triaxial test

試樣制備時,首先在一個較大空間內生成具有較小摩擦系數(0.05)的顆粒堆積體,采用伺服方法給墻體邊界施加速度直至試樣穩定地達到目標壓力值,然后再將制樣階段的接觸摩擦系數改成實際樣本摩擦系數0.5.這樣有利于減小試樣制備過程中的初始各向異性,生成較為密實的均勻試樣.

構造3 種類型的3 軸加載工況來生成數據,分別為:

(1)包含加卸載循環的常規三軸試驗;

(2)平均有效應力p不變的真三軸試驗(等p加載: σ11+σ22+σ33=p,σ11=σ22);

(3)中主應力b不變的真三軸試驗(等b加載:,σ11=常數),所有加載工況最大軸向應變為0.12.每組加載工況中,通過設定不同的卸載和再加載階段以生成不同的應力-應變關系數據,3 種加載方式共生成440 組應力-應變關系數據.

4.2 訓練變量間的映射關系

采用GRU (gated recurrent unit)循環神經網絡對有向圖內變量間的映射關系進行訓練,訓練后的映射關系用于后續有向圖本構關系的建立中.有向圖中共有5 個變量,包含1 個輸入變量ε,一個輸出變量σ,和3 個中間變量Af,Cf,φ.5 個變量間可形成36 組基礎映射關系,映射關系遵循2.2 節中的信息流動規則,這36 組基礎映射關系最終可搭建出種類繁多的有向圖鏈接關系.以映射關系{ε,φ→Cf}為例,應變值ε和孔隙率φ將共同作為GRU 網絡的輸入值,彈性剛度矩陣Cf作為輸出值,訓練時從440 組應力-應變關系數據中隨機選取340 組作為訓練集,剩余100 組作為測試集,訓練集中的34 組作為驗證集.

采用谷歌公司發布的基于TensorFlow 深度學習庫的Keras 搭建了GRU 循環神經網絡模型.本文選用的網絡層結構和超參數為:2 個神經元數量為50 的隱藏層,輸出層為采用linear 激活函數的全連接層.用于訓練的輸入和輸出數據采用最小最大歸一化處理,并且輸入變量包含當前和之前19 個歷史變量值,訓練中采用Adam 優化算法,batch size 取為256,模型在訓練400 次(epochs)后收斂.針對映射關系{ε,φ→Cf},訓練過程中GRU 模型在訓練集和驗證集上的性能表現如圖7 所示.

圖7 訓練集和驗證集上GRU 模型的訓練性能Fig.7 Training performance of GRU architecture on training and validation data

同樣地,可對其余35 組變量間映射關系訓練出GRU 神經網絡模型.值得注意的是,并不是每一組映射關系均能獲得圖7 所示的訓練性能,這主要取決于變量間是否存在內在的物理映射關系,直接影響到不同有向圖表征的本構模型得分的差異.

4.3 深度強化學習自動生成最佳“數據”本構關系

共執行Niter=10 次策略/價值網絡的訓練,每次訓練執行Ncollect=20 次完整的有向圖建立過程以收集所需的訓練數據,在激活每一條有向邊時執行NMCTS=30 次蒙特卡洛樹搜索.訓練過程中探索系數τ=1,以探索未知的使本構得分較高的動作;在利用訓練好的策略/價值網絡選擇最終本構鏈接時τ=0.01,選擇使本構得分較高的動作建立起最優的本構關系.采用卷積神經網絡訓練策略/價值網絡,卷積神經網絡的輸入值為本構有向圖的狀態空間.訓練時隨機改變策略/價值網絡的初始超參數值以研究訓練結果的收斂性,計算結果表明最終建立的有向圖鏈接是一致的.生成的最優本構關系鏈接如圖8所示,該鏈接關系在所有測試數據上打分為0.720,本構關系的信息流動路徑為:{ε→φ;ε,φ→Cf;ε,Cf→Af;ε,Af,Cf,φ→σ}.

圖8 采用強化學習算法建立的最優有向圖Fig.8 Optimal directed graph of stress-strain laws learned by deep reinforcement learning

為驗證“數據”本構模型的預測能力,在測試集中任選一組ε-σ數據,ε作為輸入值,由映射關系ε→φ計算出φ,再由ε,φ→Cf計算出Cf,ε,Cf→Af計算出Af,最后由ε,Af,Cf,φ→σ計算出最終的應力值σ.針對等p加載、等b加載以及常規加載3 種加載方式,幾組代表性工況下,由“數據”本構模型預測的應力-應變關系與離散元計算結果的比較如圖9～ 圖11 所示.可以看出,由強化學習算法生成的“數據”本構模型可很好地對測試集中的應力-應變關系及卸載、再加載過程進行預測.

圖9 等p 加載下“數據”本構模型預測比較Fig.9 Comparison between predictions and DEM simulation results for constant-p compression

需要指出的是,圖9(a)～ 圖9(c)分別是從3 個主應力和主應變的角度展示訓練模型對同一工況的預測表現.顆粒材料的主應力分量與主應變分量間并不存在唯一的對應關系.由于在等b加載工況下的最小主應力和常規三軸工況下的中主應力及最小主應力在加載過程中均保持初始圍壓不變,因此圖10 和圖11 中僅分別展示了部分主應力與主應變的對比結果.

圖10 等b 加載下“數據”本構模型預測比較Fig.10 Comparison between predictions and DEM simulation results for constant-b compression

圖11 常規三軸加載下“數據”本構模型預測比較Fig.11 Comparison between predictions and DEM simulation results for conventional triaxial compression

在強化學習訓練階段發現ε,Cf,σ3 個變量參與的本構關系可獲得最高0.702 的打分,對應的有向邊鏈接關系為{ε→Cf;ε,Cf→σ},有向圖如圖12(a)所示.該結果表明僅使用彈性剛度矩陣分量,可在數據驅動模型中最為準確的描述顆粒材料的應力-應變關系.

圖12 ε,Cf,σ 變量建立的有向圖Fig.12 Generated directed graph based on ε,Cf and σ

4.4 信息流動路徑的影響

基于有向圖的深度強化學習算法不僅能從數據中自動尋找與顆粒材料本構行為聯系最相關的物理量,還能定量地描述相關物理量之間的信息流動路徑,并且這種定量描述方式對于數據驅動模型的預測能力也是極為關鍵的.比如,與4.3 節中獲得的最佳本構訓練方式一樣,采用ε,Cf,σ3 個變量構建的鏈接關系{ε→Cf;Cf→σ},如圖12(b)所示,打分則為0.475.

本構關系{ε→Cf;Cf→σ}的預測能力如圖13 所示.從圖中可以看出,兩個模型在加載段均表現出了較高的預測能力,而在卸載、再加載階段的預測能力表現不佳.計算結果表明,本構關系的建立雖然采用了相同的內變量,但內變量間信息流動的路徑不同會生成預測精度不同的本構模型.將φ與Af作為內變量時同樣可以建立起不同打分的本構關系,作為內變量時的打分分別為0.664 (信息流動路徑:{ε→φ;ε,φ→σ})和0.453 (信息流動路徑:{ε→φ;φ→σ});Af作為內變量時的打分分別為0.657 (信息流動路徑:{ε→Af;ε,Af→σ})和0.522(信息流動路徑:{ε→Af;Af→σ}).以上結果表明,雖采用了相同的內變量來構造本構關系,但變量間不同的信息流動路徑會形成精度不同的本構關系.

圖13 有向圖鏈接{ε → Cf;Cf → σ}的預測精度Fig.13 Prediction performance of directed graph {ε → Cf;Cf → σ}

當直接將應變值ε和應力值σ鏈接成本構關系時,本構模型打分為0.708,得分高于圖12(a)中的有向圖.這表明內變量的參與會降低部分本構模型的預測精度,這是因為建模過程中激活的不同有向邊使得數據流動路徑發生了改變,而改變后的數據流動路徑需要用不同的變量間映射關系來表達.由于變量間內在的物理關系,并不能建立由GRU 神經網絡擬合的理想的映射關系,這就導致了某些映射關系并不能很好地用來建立本構關系,最終導致了本構模型預測精度降低.因此,為獲得更高的打分需要多激活精度較高的變量間映射關系,例如圖8 所示的信息流動路徑.

5 結論

本文以顆粒材料本構研究中所識別的重要相關內變量為基礎,采用有向圖和深度強化學習算法來對代表性體積單元的力學行為進行建模.一方面,采用有向圖來表征含有內變量的本構關系,采用變量間的鏈接關系來模擬數據在內變量間的流動路徑;另一方面,通過深度強化學習算法發掘數據之間的內在聯系,自主選擇鏈接有向邊,能夠使得有向圖表征的本構關系具有最優的預測精度.就計算結果討論如下:

(1)本構關系建模過程完全由數據驅動,即建模過程中不考慮內變量間是否存在理論層面的推演關系,只按照離散元模擬的應力-應變關系數據建立本構模型.

(2)建模過程采用了由葉節點向前遞歸尋找源節點,以及“一對一”、“多對一”的變量間信息流動規則,對于不同的建模需求可遵循不同的信息流動約束條件.此外,對于本構模型優劣也可采用不同的打分標準,例如,采用同時兼顧得分最高和有向邊數目最少的打分標準.

(3) 本構關系的建立是基于常規、等p、等b3 種圍壓,以及卸載后再加載的方式“激發”出的顆粒材料宏觀力學行為,除此之外,還可增加新的加載方式,產生新的訓練數據進一步提升本構模型的預測精度.

(4)有向圖與深度強化學習結合的方法,將有可能整合來自多重視角的理論模型,發展一個統一的、高質量的數據驅動本構模擬框架.

(5)借助深度強化學習算法在組合優化問題中的智能決策能力,可以在海量數據基礎上,尋找與顆粒材料本構行為最為相關的重要物理量和信息流動方式,這些發現將反過來提供新的信息,從數據驅動角度整合和促進理論本構模型的發展.