基于深度強化學習算法的顆粒材料應(yīng)力-應(yīng)變關(guān)系數(shù)據(jù)驅(qū)動模擬研究1)

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

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

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

引言

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

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

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

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

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

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

1 顆粒材料的本構(gòu)關(guān)系

1.1 與顆粒材料應(yīng)力-應(yīng)變關(guān)系相關(guān)的重要內(nèi)變量

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

其中,彈性剛度矩陣Cijmn在應(yīng)變均質(zhì)化的假設(shè)下可以表示為[37]

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

則式(2)可表示為

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

式(1)中,當切應(yīng)力和切應(yīng)變?yōu)? 時,采用增量形式表示的主應(yīng)力和主應(yīng)變關(guān)系為

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

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

1.2 基于深度強化學習的本構(gòu)研究思路

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

2 本構(gòu)關(guān)系的有向圖表示

有向圖來源于數(shù)學上的圖論,主要指采用節(jié)點來表示物理量,采用邊線來描述物理量之間相互關(guān)系的一種方法.包含內(nèi)變量的本構(gòu)關(guān)系將通過有向圖來表達[31],有向圖為信息從本構(gòu)關(guān)系輸入值開始,流經(jīng)內(nèi)變量,最終流向輸出值的信息流動網(wǎng)絡(luò).

2.1 有向圖網(wǎng)絡(luò)

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

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

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

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

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

(3)有向圖中可存在孤立節(jié)點,該節(jié)點不參與到有向圖的建立中,即本構(gòu)關(guān)系中不考慮該節(jié)點對應(yīng)內(nèi)變量對本構(gòu)關(guān)系建立的貢獻;

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

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

2.2 信息流動路徑

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

圖2 有向圖表示的兩種本構(gòu)關(guān)系Fig.2 Two constitutive laws represented by two directed graphs

2.3 本構(gòu)模型打分

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

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

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

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

3 AlphaGo Zero 深度強化學習算法

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

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

3.1 本構(gòu)關(guān)系建立過程

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

圖3 MDP 中的初始狀態(tài)和可執(zhí)行動作Fig.3 Initial state and all possible actions in MDP

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

圖4 一次完整的有向圖本構(gòu)關(guān)系建立過程Fig.4 A complete modeling process of generating a constitutive relationship

3.2 自主學習[44]

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

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

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

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

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

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

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

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

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

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

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

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

表1 利用強化學習算法建立最優(yōu)應(yīng)力-應(yīng)變關(guān)系有向圖流程Table 1 Reinforcement learning of directed graph of the stressstrain relationship

4 建立“數(shù)據(jù)驅(qū)動”本構(gòu)關(guān)系

4.1 數(shù)據(jù)制備

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

圖6 三軸數(shù)值試樣Fig.6 Numerical sample of triaxial test

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

構(gòu)造3 種類型的3 軸加載工況來生成數(shù)據(jù),分別為:

(1)包含加卸載循環(huán)的常規(guī)三軸試驗;

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

(3)中主應(yīng)力b不變的真三軸試驗(等b加載:,σ11=常數(shù)),所有加載工況最大軸向應(yīng)變?yōu)?.12.每組加載工況中,通過設(shè)定不同的卸載和再加載階段以生成不同的應(yīng)力-應(yīng)變關(guān)系數(shù)據(jù),3 種加載方式共生成440 組應(yīng)力-應(yīng)變關(guān)系數(shù)據(jù).

4.2 訓練變量間的映射關(guān)系

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

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

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

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

4.3 深度強化學習自動生成最佳“數(shù)據(jù)”本構(gòu)關(guān)系

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

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

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

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

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

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

圖11 常規(guī)三軸加載下“數(shù)據(jù)”本構(gòu)模型預測比較Fig.11 Comparison between predictions and DEM simulation results for conventional triaxial compression

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

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

4.4 信息流動路徑的影響

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

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

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

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

5 結(jié)論

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

(1)本構(gòu)關(guān)系建模過程完全由數(shù)據(jù)驅(qū)動,即建模過程中不考慮內(nèi)變量間是否存在理論層面的推演關(guān)系,只按照離散元模擬的應(yīng)力-應(yīng)變關(guān)系數(shù)據(jù)建立本構(gòu)模型.

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

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

(4)有向圖與深度強化學習結(jié)合的方法,將有可能整合來自多重視角的理論模型,發(fā)展一個統(tǒng)一的、高質(zhì)量的數(shù)據(jù)驅(qū)動本構(gòu)模擬框架.

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