基于深度强化学习算法的颗粒材料应力-应变关系数据驱动模拟研究

1、基于深度强化学习算法的颗粒材料应力-应变关系数据驱动模拟研究摘要 颗粒材料的宏观力学行为受颗粒组分等材料参数，孔隙率、配位数等状态参数的影响，同时又具备复杂 的加载路径和加载历史相关性,建立包含多个内变量以及各变量间相互关联的颗粒材料本构模型是一个重要 的科学难题.不同于传统的基于屈服面、流动法则和硬化函数框架下的唯象本构模型，本文基于颗粒物质力学 的研究基础，以颗粒材料平均孔隙率、细观组构参数和弹性刚度参数作为内变量，结合深度学习方法建立以有 向图表征的数据本构模型.有向图中以不同的链接网络表示不同的内变量信息流动方向，各个内变量间的映射 关系采用循环神经网络来建立，将各个神经网络相互组合，

2、形成包含不同内变量且具有不同预测能力的本构模 型.该本构模型的建立过程等价于在众多可能的内变量链接关系空间中寻找最能描述实际材料宏观应力应变 行为的优化问题.因此，可将有向图本构模型的建立过程看作“马尔可夫决策过程”，采用深度强化学习算法构 建有向图的内变量链接组合优化过程，具体采用AlphaGo Zero算法自动寻找最优的颗粒材料数据驱动本构模 型建模路径.研究结果表明，采用有向图和深度强化学习算法可建立起完全依靠“数据驱动”的颗粒材料应力 -应变关系.此外，本方法提供了一种将不同理论模型从数据角度统一起来,且基于人工智能算法发展更优模型 的研究思路，可为相似问题的研究提供借鉴.关键词 颗粒

3、材料，应力-应变关系，离散元法,深度强化学习,有向图DATA-DRIVEN STRESS-STRAIN MODELING FOR GRANULAR MATERIALSTHROUGH DEEP REINFORCEMENT LEARNING1)Abstract The macroscopic mechanical behaviours of granular materials are affected by not just material properties such as particle composition but also state parameters of granular

4、 assembly, like porosity and coordination number, etc. Meanwhile, granular materials are of complex loading path- and loading history-dependent features. Establishing a constitutive model incorporating multiple internal variables and their inherent relations for granular materials is an important sc

5、ientific challenge. Different from the traditional phenomenological constitutive model based on the framework of yield surface, flow rule and hardening function, this study establishes a directed graph-based data-driven constitutive model with the average porosity, fabric tensor and equivalent elast

6、ic stiffness tensor being considered as internal variables, which are critical to the constitutive behaviour of granular materials from the perspectives of the particulate mechanics. The constitutive models containing different internal variables and having different predictive capabilities are repr

7、esented by different directed graphs with various internal variables linking networks. The recurrent neural network is trained to represent the source-target mapping relationships of the information flows between internal variables. The process of establishing constitutive models is simplified as a

8、sequence of forming graph edges with the goal of finding the optimal combination of internal variables. Therefore, the modelling of the constitutive model can be formulated as a Markov decision process and implemented by the deep reinforcement learning algorithm. Specifically, the well-known AlphaGo

9、 Zero algorithm is used to automatically discover the optimal data-driven constitutive modelling path for granular materials. Our numerical examples show that this modeling framework can produce constitutive models with higher prediction accuracy. Furthermore, this study provides a new research para

10、digm by integrating different theoretical models from the point of data and leveraging the algorithms in artificial intelligence to develop a superior model. The same idea can be extended to seek new insights for similar scientific problems.Key words granular materials, strain-stress relationship, d

11、iscrete element method, deep reinforcement learning, directed graph引言在研究泥沙、岩土、粉体等颗粒材料宏观力学 行为时，通常在连续介质力学框架下建立其唯象本 构模型，通过引入内变量来描述颗粒材料的路径相 关、历史相关等特性虽然基于连续假设的唯象 本构模型在描述颗粒材料宏观力学特性时具备很好 的工程适用能力，但是颗粒材料由于其自身的离散 特性，具有摩擦性，剪胀性,压硬性，各向异性，荷载 路径和历史相关性等典型特征，并且在外载作用下 还会呈现出应变局部化等复杂的演化机制，连续模 型要体现出颗粒材料精细的宏观特性就不得不引入 繁杂的

12、唯象假设和众多待标定的自由参数.离散元方法6(DEM)可跟踪并计算每个散体单 元的运动信息及单元间的接触作用力，在颗粒材料 力学行为研究中得到了广泛的应用.随着计算机硬 件计算能力的提升，以及颗粒几何特征的精准刻 画7-9和高精度接触算法的提出1-12,使得离散元方 法的应用领域和计算精度得到大幅提升.然而离散 元方法用于直接模拟较大空间尺度和时间尺度的工 程问题仍需要巨大的计算成本.为此，人们提出了多 尺度计算方法以提升颗粒材料问题的计算效率13-14. 在此框架中，通常基于均匀化思想采用低尺度的离 散元计算在宏观模拟中替代传统的唯象本构关系. 均匀化即从一定体积的颗粒聚集体内部颗粒间作用

13、的微观物理量获得在该体积内连续性力学描述的过 程15.而均匀化的首要任务是创建代表性体积单元 (RVE)16, RVE是包含若干颗粒的颗粒聚集体，是连 接颗粒介质微结构尺度和宏观工程尺度的桥梁.颗 粒介质多尺度方法常采用DEM与FEM相结合的 方式来实现，即在FEM计算中提供施加给RVE的 变形条件,RVE颗粒试样在变形条件下产生与应变 对应的应力值，随后将应力-应变关系或者切向刚度 矩阵参数映射到FEM网格中高斯积分点上进行下 一时步的求解17-20.采用多尺度方法虽然避免了对 颗粒介质全尺度的DEM求解，但需要在FEM网格 每个高斯积分点上构建单独的RVE21-22,并且需要 在计算中维护

14、与FEM网络中高斯点数量相当的 RVE加载状态信息，这又增加了计算机的存储与计 算开销19,对于跨多个尺度的问题计算量会更大.为此，人们提出了“离线训练(offline training) 在线模拟(online prediction)”的计算策略23-24来提 升计算效率.相比于直接采用RVE模拟结果的“在 线计算来说，“离线策略下的RVE模拟只产生用 来训练“代理模型”的数据集25-26,即采用机器学习 算法训练并获得RVE的应力-应变关系，从而宏观 尺度FEM的计算不再需要与细观尺度RVE的求解 交替进行，在FEM计算阶段直接读取RVE的应力- 应变关系，这就大大提高了计算效率27.可米

15、用不同的神经网络模型对颗粒介质的本构 关系进行训练与计算，常用的基础模型有经典BP神 经网络模型29和循环神经网络模型3,二者均可 将应力、应变值分别作为输入、输出量进行训练. 循环神经网络中的长短期记忆网络(LSTM)和门控 循环单元网络(GRU)具有记忆功能，对于与加载路 径和加载历史相关的颗粒材料而言，在应力-应变关 系建立方面表现出了较强的适用性和较高的拟合 精度26.在常规循环神经网络模型基础之上，也有一些 学者结合其他理论模型进一步研究了颗粒材料力学 行为.Wang等容结合博弈论和强化学习模型研 究了颗粒材料的界面位移与张力之间的力学行为以 及本构关系.瞿同明等34结合细观力学理论

16、和循环 神经网络模型训练了颗粒材料的应力-应变关系. Karapiperis等35在有限元中采用了基于离散元数据 训练得到的循环神经网络本构模型.QU等36基于 坐标不变性原理讨论了从主应力-应变空间到一般 应力-应变空间的数据驱动本构模型训练策略.综上所述，建立精度高、计算快的RVE应力- 应变关系对于颗粒介质的多尺度模拟是至关重要 的.本文借助有向图方法和AlphaGo Zero深度强化 学习算法，在颗粒物质力学的研究基础之上，试图建 立一种自动寻找最优内变量间映射关系、完全依靠 “数据驱动的RVE力学行为建模方法.1颗粒材料的本构关系1.1与颗粒材料应力-应变关系相关的重要内变量离散元方

17、法可有效地在细观尺度下模拟颗粒材 料的宏观力学行为，本文将基于球形离散单元法线 性接触模型模拟颗粒材料的应力-应变关系.从细观 力学的角度，颗粒材料的宏观应力-应变增量关系可 表示为Ab =j,n 1; 2; 3(1)其中，弹性刚度矩阵CIjmn在应变均质化的假设下可 以表示为37NcCijmn =(Lk 洒球 AVk=1其中，kn和ks为颗粒之间的法向和切向接触刚度; V为颗粒试样的体积;&为克罗内克函数;时为每个 接触对内两个颗粒的中心距离；奸为接触k的方向 向量在I方向上的分量;N为接触对的数量，令1 Ncaijmn V(Lk契鸟殊由V k=1则式(2)可表示为Cijmn (kn ks)

18、aij mn+ ksbijmn(4)(5)弹性刚度矩阵Cjmn可由颗粒间的法向接触刚 度kn和切向刚度ks,以及单位体积内的细观结构张 量表示.当颗粒材料承受外载荷作用时，颗粒间发生 相对滑移，导致组构张量在外载作用下会不断发生 变化.式(1)中，当切应力和切应变为0时，采用增量 形式表示的主应力和主应变关系为AbiiAb22Ab332112CCCC1 C2 C222212 312CCC3333121理路AAA其中，C1122 = C2211, C1133 = C3311, C2233 = C3322, 弹性 刚度矩阵中共有6个独立分量表示主应力和主应变 之间的关系.除了上述推导所涉及的弹性刚

19、度矩阵外，颗粒 物质力学研究表明颗粒材料的孔隙率0 (对于单相 颗粒材料，孔隙率定义为一定体积内非颗粒组分所 占据的体积分数)和组构参数对于其宏观力学行为 具有重要的影响.其中，岩土材料作为一种典型颗粒 材料，现代土力学理论的一个突破性进步就在于将 孔隙率等状态参数作为硬化法则的一部分，将其考 虑到岩土材料的本构关系研究中取-39.近年来岩土 本构研究领域另一个重要进步在于认识到颗粒材料 的组构对于其本构行为的重要影响40,考虑组构参 数的演化过程是近年来岩土本构理论的重要发展趋 势41-43.通常颗粒材料的细观组构有多种张量描述 方式，本文采用式(3)和式(4)中的4阶组构张量来 共同描述颗

20、粒材料的组构演化过程.1.2基于深度强化学习的本构研究思路颗粒材料的本构研究从不同的视角发展了不同的理论模型和方法，这些理论方法对问题本身的研 究都取得了新的进展.如何将这些来自不同视角的 重要发现整合到一个统一的框架内是一个思路合理 但是实践困难的挑战.本文借助深度强化学习的思 路，将从颗粒物质力学研究中得到的弹性刚度矩阵 (式(2),用Cf表示)，孔隙率参数0和组构张量(式,用4f表示)结合在一起，试图在这些重要变量的 数据样本基础上，借助人工智能算法整合已有的知 识资源，发掘新的认识，自动发现能描述颗粒材料应 力应变行为的本构关系.2本构关系的有向图表示有向图来源于数学上的图论，主要指采

21、用节点 来表示物理量，采用边线来描述物理量之间相互关 系的一种方法.包含内变量的本构关系将通过有向 图来表达El,有向图为信息从本构关系输入值开始, 流经内变量,最终流向输出值的信息流动网络.2.1有向图网络图1为以内变量和输入输出值构建的变量间信 息流动的有向图，图中的箭头代表了变量间的“源- 目标数据流动方向.有向图中的变量称为节点，变 量间的链接称为边，对于每个链接而言，箭头指出的 节点称为源节点，箭头指向的节点称为目标节点，整 个有向图的输入称为根节点，输出称为叶节点.一个有向图代表一种本构关系，为保证本构关 系中变量间信息流动的合理性，有向图的构建需遵 循以下规则：有向图中仅存在唯一

22、的叶节点(该节点不作 为其他任何节点的源节点,仅作为目标节点)；有向图中仅存在唯一的根节点(该节点不作 为其他任何节点的目标节点，仅作为源节点)；有向图中可存在孤立节点，该节点不参与到 有向图的建立中，即本构关系中不考虑该节点对应有向图中不存在闭环信息流(即信息流从某 一结点开始，最后又流回该节点)；如果某一节点链接到了根节点或叶节点，那 么此节点需包含在从根节点到叶节点的完整路径中. 2.2信息流动路径对于简单的有向图链接，如图2(a)所示,信息流 动方向可直接表示为 - 0 -Af - Cf C,而对于 复杂的有向图链接，如图2(b)所示，需要建立合理的 信息流动路径.图2(b)中的信息流

23、动路径需要从有 向图输出变量a开始以递归的方式向前寻找信息 “源每一节点只寻找直接指向该节点的“源节点, 形成“一对一或“多对一的信息映射关系.图2(b) 中信息流动路径对可表示为：( - Cf, (, Cf Af, (, Af 一 0及, 0, Cf 一 a,即信息从流入有向图 后计算出Cf,随后和Cf共同计算出Af,接着由 和Af共同计算出0,最后由, 0, Cf共同计算出a, 实现了 到a的信息流动.将采用循环神经网络拟 合出每对路径中变量间的映射关系，用于后续“数 据本构模型的建立.(a)(b)图2有向图表示的两种本构关系Fig. 2 Two constitutive laws rep

24、resented by two directed graphs2.3本构模型打分对有向图建立本构关系的预测精度进行打分, 以评估通过数据建立的本构关系对于训练数据的泛 化能力.通过比较信息流路径计算出的应力值与训 练数据中应力值的均方误差(mean squared error, MSE),评价强化学习算法中不同有向图代表的本构 模型的优劣.对于具有Mata个数据样本的样本集，每个样本 的均方误差可表示为i Nf2MSEi =万 Sj(Y#ta)_Vf 7=1其中，Yjta和匕? 1是第j个数据样本(应力-应变曲 线)的第j个特征值，样本分别来自训练数据值和本 构模型预测值;Nf是输出特征的数量

25、(对应于本构 模型应力、应变的3个分量值),Sj是数据标准化中 的缩放因子(采用最小最大原则对数据进行归一化 处理).对集合MSEj内各均方误差值进行升序排序, 并取集合中第P%的MSE值(记为P%)计算本构 模型打分值SS = max lgm/p%E(8)lg crit式中crit 1 ,是需提前设定的MSE误差标准，以评 估生成本构模型的精确度，其中P%取值90%,即集 合中第90%的MSE值，以及crit = 1X 10-6 .3 AlphaGo Zero深度强化学习算法本节将采用AlphaGo Zero中的强化学习算法 建立基于有向图的本构关系.AlphaGo是Google DeepM

26、ind团队开发的一个基于深度神经网络的围 棋人工智能程序.2015年10月AlphaGo击败了欧 洲冠军Fan Hui, 2016年3月击败韩国棋手Lee Sedol. 2017 年 10 月，DeepMind 公布了最新的 AlphaGo Zero,在与AlphaGo对战中取得了 100 : 0的战 绩.与AlphaGo不同的是,AlphaGo Zero未使用任何 围棋领域的专业知识和人为干预，训练数据全部来 自于自我对弈(self-play),算法方面，AlphaGo Zero 只采用了一个神经网络,将棋盘状态作为输入特征, 来指导蒙特卡洛树搜索(Monte Carlo tree sear

27、ch, MCTS)探索最优的落子动作.本构模型建立过程中, 选定输入变量、输出变 量和内变量后，可以通过“穷举法列出所有满足信 息流动规则的有向图，然后对每一种有向图打分，得 分最高者即为最优本构关系.对于内变量较少的本 构关系建立过程，可考虑采用“穷举法寻找出最优 的本构关系，但如果引入的内变量的数量较多，采用 直接寻找的方法将极大地增加计算成本，因此建立 一种能够自动寻找最优本构关系的算法是必要的. 本构关系建立过程与AlphaGo Zero实现过程类似, 同样可采用马尔可夫决策过程(Markov decision process, MDP)对有向图中变量链接关系进行选择, 以建立得分最高

28、的有向图链接关系.下面将给出采 用AlphaGo Zero算法建立本构关系有向图的流程.3.1本构关系建立过程本构关系建立过程可表述为马尔可夫决策过 程，本构关系的初始状态如图3(a)所示，本构关系建 立的过程就是智能体在初始状态下对有向图中的边 进行激活，并获得最大奖励的过程.所有可能被激活 的边称为“动作空间，并遵循有向图网络建立准则, 所有动作如图3(b)所示,可描述为 一 Af, 一 Cf, 8 如 c, Af Cf, Af 如 Af c, Cf Af, Cf 0, Cf c,放 Af,放 Cf, 8 c.马尔可夫决策过 程中的状态S”为动作空间中有向图边的激活状态; “动作/为在当前

29、状态下准备激活的有向图边;“奖 励,为根据本构关系打分情况确定的奖励值，当打 分高于给定的分值时取1,否则取-1.(a)初始状态(a) Initial state(b)动作空间(b) All actions space图3 MDP中的初始状态和可执行动作Fig. 3 Initial state and all possible actions in MDP在有向图当前状态St下,智能体执行动作at (激 活某一条有向边)，当前状态会转移到t + 1时间步的 新状态&+1,其中动作at选自于&状态下每一个有 效动作概率的集合n(st),同时智能体获得一个奖励 S1,该奖励值的获得是由于智能体在状

30、态st下执行 了动作牛但是在实施过程中具体的奖励值rt只有在本构关系完全建立后才能获得.当本构关系得分ry T = 1,反之,W T = I-图4为由马尔可夫决策过高于给定分值时，本构建立过程中所有动作的奖励程表示的一次完整的本构关系建立过程.initial state s0state S1state %0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 00, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 00, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0saction % =action % =e-今J final stat

31、e s5state &state &0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 1, 1, 00, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 1, 00, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0state SiC.score/rewardpredict policy/valuetrain policy/value network f (s)actionstate SiC.score/rewardpredict policy/valuetrain policy/value network f (s)actionn -initial

32、state soaction%nfinal state sT action action图5 AlphaGo Zero实现的自主学习过程Fig. 5 Self-play reinforcement learning in AlphaGo Zeromonte carlo treesearches% Vo3.2自主学习3在AlphaGo Zero算法中，主要通过训练一个超 参数为9的策略/价值网络fg来指导本构关系有向 图建立过程中动作的选择.策略/价值网络将有向图 的链接状态作为输入，输出一个向量p和标量v,即 (p, v) = fg (s), p表示在当前状态下链接每一条有向 边的概率,v表示

33、在当前状态下本构关系获得最高分 的概率.策略/价值网络通过自主学习来提升策略选 择的能力，使得有向图链接得分最高.策略选择能力 的提升通过MCTS算法来实现，而在每一状态下蒙 特卡洛树扩展过程中的动作概率则是通过策略/价 值网络fg进行计算,AlphaGo Zero实现的自主学习 过程如图5所示.自主学习过程中需运行多次MCTS 算法使策略/价值网络fg变得更强，而MCTS算法的 实施均采用更新后的策略/价值网络.蒙特卡洛搜索树中的每个节点表示有向图的链 接状态s,搜索树中的每条边(s, a)表示在状态s下 允许的动作a.搜索树中每条边维持一个统计量集合 N(s, a), W(s, a), Q

34、(s, a), P(s, a),其中 N(s, a)为在树 搜索过程中节点被访问的次数,W(s, a)为总动作价值, Q(s, a)为平均动作价值,P(s, a)为选择边的先验概率.选择(select):自主学习过程中从搜索树的根 节点s开始,直至在时间步L遇到之前未被遍历的 叶节点辎在t 敢；,敢- Cf； , Cf 一4f； , 4f, C& 敢 一 e.为验证“数据”本构模型的预测能力，在测试集 中任选一组2数据,作为输入值，由映射关系 敢计算出敢,再由,敢 Cf计算出Cf, , Cf Af计算出Af,最后由s, Af, Cf,放一 g-计算出最终的 应力值g.针对等p加载、等b加载以及

35、常规加载 3种加载方式,几组代表性工况下，由“数据”本构模 型预测的应力-应变关系与离散元计算结果的比较 如图9图11所示.可以看出，由强化学习算法生图8采用强化学习算法建立的最优有向图Fig. 8 Optimal directed graph of stress-strain laws learned by deepreinforcement learning图9等p加载下“数据”本构模型预测比较Fig. 9 Comparison between predictions and DEM simulation results forconstant-p compression(b)易3 - 3

36、3图10等b加载下“数据”本构模型预测比较Fig. 10 Comparison between predictions and DEM simulation resultsfor constant-b compression图11常规三轴加载下“数据”本构模型预测比较Fig. 11 Comparison between predictions and DEM simulation resultsfor conventional triaxial compression成的“数据”本构模型可很好地对测试集中的应力- 应变关系及卸载、再加载过程进行预测.需要指出的是，图9(a)图9(c)分别是从3

37、个 主应力和主应变的角度展示训练模型对同一工况的 预测表现.颗粒材料的主应力分量与主应变分量间 并不存在唯一的对应关系.由于在等b加载工况下 的最小主应力和常规三轴工况下的中主应力及最小 主应力在加载过程中均保持初始围压不变，因此 图10和图11中仅分别展示了部分主应力与主应变250250的对比结果.0.020.040.060.080.020.040.060.08在强化学习训练阶段发现, Cf, a 3个变量参与 的本构关系可获得最高0.702的打分，对应的有向边 链接关系为一 Cf; , Cf 一町,有向图如图12(a)所 示.该结果表明仅使用弹性刚度矩阵分量，可在数据 驱动模型中最为准确的

38、描述颗粒材料的应力-应变 关系.(b)打分 0.475(b) Score = 0.475图12 e, Cf, a变量建立的有向图Fig. 12 Generated directed graph based on e, Cf and a4.4信息流动路径的影响基于有向图的深度强化学习算法不仅能从数据 中自动寻找与颗粒材料本构行为联系最相关的物理 量，还能定量地描述相关物理量之间的信息流动路 径，并且这种定量描述方式对于数据驱动模型的预 测能力也是极为关键的.比如，与4.3节中获得的最 佳本构训练方式一样，采用e, Cf, a 3个变量构建的 链接关系e 一 Cf； Cf 一 a,如图12(b)所示,打分则 为 0.475.本构关系e 一 Cf； Cf 一 a的预测能力如 图13所示.从图中可以看出，两个模型在加载段均 表现出了较高的预测能力，而在卸载、再加载阶段 的预测能力表现不佳.计算结果表明,本构关系的建 立虽然采用了相同的内变