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

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

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

引言

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

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

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

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

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

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

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

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

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

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

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

則式(2)可表示為

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

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

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

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

1.2 基于深度強(qiáng)化學(xué)習(xí)的本構(gòu)研究思路

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

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

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

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

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

圖1 有向圖中的信息流動(dòng)Fig.1 Information flow in a directed graph

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

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

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

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

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

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

2.2 信息流動(dòng)路徑

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

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

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

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

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

3 AlphaGo Zero 深度強(qiáng)化學(xué)習(xí)算法

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

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

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

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

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

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

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

3.2 自主學(xué)習(xí)[44]

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

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

圖5 AlphaGo Zero 實(shí)現(xiàn)的自主學(xué)習(xí)過(guò)程Fig.5 Self-play reinforcement learning in AlphaGo Zero

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

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

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

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

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

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

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

損失函數(shù)綜合考慮了獎(jiǎng)勵(lì)值的均方誤差和交叉熵?fù)p失.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖7 訓(xùn)練集和驗(yàn)證集上GRU 模型的訓(xùn)練性能Fig.7 Training performance of GRU architecture on training and validation data

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

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

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

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

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

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

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

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

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

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

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

4.4 信息流動(dòng)路徑的影響

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

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

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

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

5 結(jié)論

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

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

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

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

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

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