高 军,黄再兴

(1.上海民航职业技术学院，上海 200232；2.南京航空航天大学机械结构力学及控制国家重点实验室，江苏 南京210016)

多种群遗传算法在PBX本构模型参数识别中的应用*

高 军1,黄再兴2

(1.上海民航职业技术学院，上海 200232；2.南京航空航天大学机械结构力学及控制国家重点实验室，江苏 南京210016)

利用多种群并行结构对标准遗传算法SGA进行并行化处理，引入移民算子和精华种群形成多种群遗传算法MPGA，并设计了自适应交叉和变异概率对算法的收敛速度进行改进。结合ABAQUS软件和改进的多种群遗传算法，建立了材料本构模型参数识别方法。采用该方法对PBX炸药黏弹性损伤本构模型参数进行了模拟识别，并同基于标准遗传算法的参数识别方法进行了比较。结果证明，基于改进多种群遗传算法IMPGA的方法对克服算法未成熟收敛有显著的效果，识别结果更稳定。同时该方法的收敛速度更快，寻优能力更强，适合复杂非线性问题的优化，此方法可以被应用到其他材料本构模型的参数识别中。

固体力学；参数识别；多种群遗传算法；PBX炸药；ABAQUS；本构模型

高聚物黏结炸药(polymer bonded explosive, PBX)是以高能单质炸药为主体，加入黏结剂等物质的固体炸药，在常规武器战斗部中具有广泛应用[1]。PBX在制造、运输及发射等过程中经历的复杂应力过程，影响着武器弹药的安全和使用。因此，对PBX的力学行为进行有限元模拟研究具有重要意义[2]。PBX力学行为的有限元模拟依赖于本构模型的正确性，一个合适本构模型中参数的准确性对于PBX的有限元模拟至关重要。材料本构的参数识别属于带有约束条件寻找最优解的优化问题，针对该问题，目前已有一系列的研究报道[3-4]。从研究现状来看，PBX炸药本构模型具有强非线性特点，其参数无法通过对标准拉伸等实验数据直接处理获得。另外，参数空间具有非凸性，使参数搜索可能陷入局部最优值。所以，传统优化方法不能保证获得待识别参数的全局最优解。

随着智能分析方法的发展，遗传算法、蚁群算法等被不断引入参数识别中。其中，遗传算法在处理复杂非线性问题具有独特的优势，得到了广泛的应用。B.M.Chaparro等[5]、G.H.Majzoobi等[6]通过遗传算法，对多种本构模型参数进行了识别；P.A.Muoz-Rojas等[7]结合自编有限元程序和遗传算法，进行参数识别；陈炳瑞等[8]将遗传算法和神经网络结合，进行岩体流变参数的识别；高军等[9]通过遗传算法，对PBX本构模型参数进行了识别。然而遗传算法虽然简单、通用，但算法本身受遗传操作参数、群体规模等条件影响，客观存在着早熟收敛等问题，最终可能无法获得最优解。解决这些问题的一个有效途径是调整遗传算法的搜索方式。随着大规模并行机的快速发展，结合遗传算法固有的并行性，遗传算法的多种群并行化成为解决这些问题的有效方法[10]。

本文中，采用一种改进的多种群并行遗传算法来取代传统的遗传算法，并结合有限元软件ABAQUS，建立材料本构模型参数识别方法。并利用该方法，对PBX炸药黏弹性损伤本构模型参数进行模拟识别，与标准遗传算法识别结果进行比较。

1 数学模型

本构模型的参数识别是，在已知材料本构模型的基础上，结合同种材料试件的相关实验数据，合理地选择材料本构模型参数，使材料模型可以真实地模拟材料的物理特性。本构模型中的待识别参数，可以表示为[11]：

X=(x1,x2,…,xm)T

(1)

式中：m是待识别参数的个数。

通过实验测量可以获得选定测点处的实验值，如应变、位移等。测点处实验值可表示为：

(2)

测点处值是材料参数变量X=(x1,x2,…,xm)T的函数，则有：

(3)

测点处的模拟计算值，可以在给定待识别参数X=(x1,x2,…,xm)T的初始值时通过数值计算获得，表示为：

(4)

式中：u1、u2、…、un为测点处模拟值。

参数识别就是寻找一组材料参数，使模拟计算结果和测量结果之间的差值最小。假设计算两者之间差值的目标函数为：

(5)

则参数识别的过程就是寻找一组材料参数，使下式成立：

(6)

这组材料参数X*即为本构模型参数的最优解。通常，材料参数在一定的允许范围内约束条件为：xi,min≤xi≤xi,max(i=1,2,…,m)，xi,min和xi,max为材料参数的允许范围。

2 多种群遗传算法

遗传算法(standard genetic algorithm, SGA)是根据生物界自然选择和自然遗传机制发展而来的随机搜索算法。它基于群体运算的方法，对每个个体进行的各种运算都具有一定的相互独立性，所以具有一种隐含的并行性[12]。根据这种并行性，遗传算法的并行结构分为单种群主从并行结构、细粒度并行结构、多种群并行结构和分级混合并行结构。其中多种群并行结构计算效率高，结构实际操作简单，因此该结构被广泛应用[13-14]。通过多种群并行结构对遗传算法进行并行化处理，形成多种群遗传算法(multiple-population genetic algorithm, MPGA)。

多种群遗传算法的算法结构如图1所示。图中种群1～N均是基本的遗传算法，N个种群均运行基本的遗传算法，同时进行优化搜索。进化过程中，通过移民操作进行各种群之间的信息交换，将源种群中的最优个体利用移民算子定期地迁移到目标种群，并同目标种群的最劣个体比较后淘汰较差个体。移民操作将各种群相互联系，否则MPGA相当于独立进行多次SGA计算。

为了更好地模拟生物进化过程，并行进化时对各种群选择不同的遗传操作参数，实现不同的搜索目的。在SGA中，交叉和变异概率的恰当选择决定着全局和局部搜索能力的均衡。一般建议取较大的交叉概率(0.7～0.9)和较小的变异概率(0.001～0.1)。对于不同的取值组合，优化结果也具有较大的差异。MPGA通过多个赋以不同遗传操作参数的种群协同进化,在取值范围内随机产生多个遗传参数组合，兼顾算法的全局搜索和局部搜索。

进化过程中，利用人工选择算子选出各种群的最优个体放入精华种群进行保存，并从中选取全局最优个体作为算法的最优解。精华种群同时作为判断算法终止的条件，采用最优个体最少保持代数作为终止判据[15]。

图1 MPGA算法结构图Fig.1 Structure of MPGA

3 算法收敛速度改进策略

基本遗传算法中，交叉和变异概率均为一个固定常量，不同适应度的个体被选择进行交叉和变异操作的概率相同，使子代个体的产生随机性较强，父代的优良基因易被破坏，对算法的收敛和计算速度产生不利的影响。本文中，设计了一种与父代种群整体适应度相关联的自适应交叉和变异概率：

(7)

(8)

式中：Pc0和Pm0为各种群中预先给定的交叉和变异概率；Pc,min和Pm,min为最小交叉和变异概率，这里分别取0.7和0.001；Favg为父代种群中个体的平均适应度，Fmax为父代种群中个体的最大适应度。Favg和Fmax计算公式如下：

(9)

(10)

这种自适应变化概率能够根据遗传进化过程中父代种群的整体适应度不断地自动调整交叉和变异概率，可以增大遗传算法获得全局最优解的可能性，同时可以提高遗传算法的收敛速度。将自适应交叉和变异概率引入多种群遗传算法中，形成改进的多种群遗传算法(improved multiple-population genetic algorithm, IMPGA)。

4 参数识别的实现

基于多种群遗传算法的本构模型参数识别方法由正演和反演两个部分组成。正演部分是在参数识别过程中，计算一组参数所对应的模拟值。本文中，通过有限元软件ABAQUS的二次开发接口UMAT，将本构模型嵌入到ABAQUS进行计算。反演部分由多种群遗传算法组成，本文中采用MATLAB语言进行计算程序的编写。具体实现步骤如下。

(1) 初始化。种群中个体采用实数编码，在各待识别参数范围内，随机产生种群数N组初始种群。种群数和种群中群体规模根据实际情况选择。每个染色体个体(x1,x2,…,xm)为m维向量，其中xi(i=1,2,…,m)为待识别参数。随机产生N组交叉和变异概率，分别赋予各组种群。

(2) 正演计算。调用ABAQUS命令程序流，利用真实模拟试验件几何尺寸和实验条件的有限元模型，分别计算N组种群中每个参数个体对应的模拟计算值。

(3) 个体适应度计算。利用测点处实验值和模拟计算值，通过式(5)计算每个个体的目标函数值。本文中参数识别为求最小值问题，因此适应度亦可取式(5)，则适应度为：

(11)

图2 参数识别的计算流程Fig.2 Calculation process of parameter identification

(4) 移民操作。根据适应度大小判断个体优劣，将源种群中的最优个体迁移到目标种群，并同目标种群的最劣个体比较，再淘汰较差个体。同时根据适应度，通过人工选择算子选出各种群的最优个体放入精华种群MaxChrom中。

(5) 终止判断。选出精华种群MaxChrom中的最优个体，同上一代中选出的最优个体进行比较。若两者相同，最优个体保持代数gen=gen+1，否则令gen=0，同时将本次选出的最优个体代替全局最优个体进行保存。最优个体保持代数gen同设定的最少保持代数MINGEN比较，若gen>MINGEN则识别过程终止，此时的全局最优个体为参数识别的最优解；否则，继续下一步。

(6) 对各种群中的个体分别执行选择、交叉、变异操作，生成下一代种群，再返回到步骤(2)进行新一轮的优化。遗传操作中采用轮盘赌选择、实数交叉和随机变异方法。

具体计算流程如图2所示。

5 参数识别方法的比较

通过基于多种群遗传算法MPGA和IMPGA的参数识别方法，对PBX黏弹性损伤本构模型参数进行模拟识别，验证该方法的正确性和可靠性，并同基于基本遗传算法SGA的参数识别方法的结果进行比较。

5.1 PBX黏弹性损伤本构模型

Visco-scram方程被广泛应用于PBX炸药的力学行为描述。在该方程中，引入表示内部细观缺陷的损伤变量c，表示微裂纹的特征尺寸(平均尺寸)，则PBX炸药黏弹性损伤本构模型[9]为：

(12)

式中：σm为体积应力，S为偏应力，η为黏性系数，G为剪切模量，K为体积模量，βm与βs是和微裂纹扩展相关的两个常数。基于文献[16]，损伤变量c的演化方程具有幂律形式：

(13)

式中：χ与n是两个材料常数。本构模型中，共涉及7个材料参数，分别为杨氏模量E、泊松比ν、黏性系数η及4个常数βm、βs、χ与n。通过ABAQUS二次开发用户材料子程序UMAT，将PBX黏弹性损伤本构模型嵌入ABAQUS中，应用于PBX材料。

参数识别正演部分分析模型采用PBX炸药圆柱体试件模拟压缩，试件长度6 mm，直径10 mm，单元类型为C3D10M四面体实体单元。一端固定，另一端通过位移加载形式进行压缩。压缩加载最大位移0.2 mm，加载时间1 s。取待识别参数的一组合理值为真实值，如表1。利用这组参数进行模拟压缩，得到选定测点的数据，将此数据作为实验数据进行本构模型的参数识别。通过比较参数识别结果和参数真实值的差距，分析识别方法的优劣。

表1 参数识别结果Table 1 Results of parameter identification

图3 测点处的载荷-位移曲线Fig.3 Load-displacement curve of measurement point

参数识别中，避免多解的一般方法是，预先确定待识别参数的取值范围，将参数取值约束在变化范围内。但是，在实际操作中取值范围难于直接确定，只能根据理论分析或工程经验来预估。通常可以在符合理论前提下暂定较大的取值范围，通过参数识别确定一组参数，再进行其他工况下的验证。PBX本构模型参数取值范围，见表1。

参数识别中，实验数据的类型和质量直接影响参数识别问题的适定性，即解的存在性、稳定性和唯一性问题。因此，实验数据需要尽量选择对本构参数敏感度较大的，另外测点位置也要根据实际问题进行选择。本文中，选取圆柱体加载端表面一点为测点，取模拟载荷值作为实验数据用于参数识别，该点处载荷-位移曲线如图3所示。

5.2 算例结果比较

分别采用基于SGA的参数识别方法、本文中建立的基于MPGA和IMPGA的参数识别方法，对PBX黏弹性损伤本构模型进行参数识别。

基于SGA参数识别方法中的操作参数取值分别为：种群规模20，迭代次数20次，交叉概率0.8，变异概率0.1。图4为SGA方法运行5次得到的种群最优解适应度随进化代数的变化曲线。从图4可以看出，5次得到的适应度均不相同，相应的识别得到的参数也各不相同。其中一次识别过程中适应度较大时已停止进化，说明算法陷入某个局部收敛区。更改遗传操作交叉和变异概率后进行多次识别计算，最终适应度均不能收敛到一个稳定的范围。所以，对于待识别参数较多且具有强非线性的本构模型，通过标准遗传算法进行参数识别时结果并不稳定，需要计算多次取相对最优解。

采用MPGA和IMPGA方法进行参数识别时,取种群数N为10, 每个种群的规模均为10。各种群的初始交叉和变异概率分别在区间(0.7～0.9)和(0.001～0.1)中随机选取，IMPGA方法中的交叉和变异概率在进化过程中进行自适应调整，算法终止条件是最优个体最少保持代数为10。采用两种方法计算运行5次得到最优解的适应度变化曲线，如图5～6所示。从图5～6可以看出，两种方法运行5次的适应度结果都能趋近于稳定值,说明稳定性很好。和SGA方法相比，多种群遗传算法在更少的进化次数下收敛到一个稳定值，收敛速度快而且寻优能力更强。IMPGA方法和MPGA方法相比，适应度从进化初期快速下降，进化代数也更少，基本在10代内已趋于稳定。IMPGA方法收敛速度更快，寻优能力更强，计算效率也更高。

图4 SGA方法适应度Fig.4 Fitness value with SGA

图5 MPGA方法适应度Fig.5 Fitness value with MPGA

图6 IMPGA方法适应度Fig.6 Fitness value with IMPGA

图7 载荷-位移曲线对比Fig.7 Contrast of load-displacement curves

对5次采用SGA、MPGA和IMPGA方法的参数识别结果中最优解进行对比，见表1。从表中可以看出，采用3种方法的结果均满足工程计算基本要求，但是采用多种群遗传方法的结果明显优于采用SGA方法的，MPGA和IMPGA方法的结果较接近。将通过3种方法获得的参数代入正演分析模型中，计算参数对应的测点载荷和位移，图7为计算结果和实验数据的载荷-位移曲线对比。从图7可以看出，3种方法识别得到的参数模拟曲线和实验曲线整体上都比较吻合, 但MPGA和IMPGA方法的结果优于SGA方法的。采用SGA的参数识别方法在优化结果上不稳定，但是在进行多次识别计算后，也可以获得较优解。而采用MPGA和IMPGA的参数识别方法，运行5次的结果最优值近乎相等，算法稳定性更好，识别效率更高。

可见，多种群遗传算法MPGA中，通过多个种群同时对参数空间进行搜索，兼顾了全局和局部搜索能力，降低了算法对遗传操作参数的敏感性，对克服未成熟收敛有显著的效果，同时算法寻优结果比标准遗传算法SGA的稳定性更好，计算效率更高。多种群遗传算法中，IMPGA比MPGA算法收敛速度更快、寻优能力更强。

6 结 论

通过多种群并行结构对标准遗传算法进行并行化处理，并引入移民算子和精华种群形成多种群遗传算法。利用移民算子实现各种群之间的信息交换，采用精华种群作为最优个体储存空间和算法终止判据。针对算法收敛速度的改进，设计了与父代种群整体适应度相关联的自适应交叉和变异概率。自适应交叉和变异率随着遗传进化过程不断的自动调整，可以增大遗传算法获得全局最优解的可能性。

基于改进的多种群遗传算法和ABAQUS软件，建立了材料本构模型参数识别方法。通过对PBX黏弹性损伤本构模型参数进行模拟识别，验证了该方法的正确性和可靠性，并同基于标准遗传算法的参数识别方法的结果进行比较。算例证明，基于改进多种群遗传算法的参数识别方法识别结果的稳定性更好，对克服未成熟收敛有显著的效果。同时，该方法收敛速度更快，寻优能力更强，计算效率更高，适用于复杂非线性材料参数的识别。

[1] 梁增友．炸药冲击损伤与起爆特性[M]．北京：电子工业出版社，2009.

[2] 郭虎，罗景润．循环载荷下 PBX力学行为研究[J]．爆炸与冲击，2013，33(增刊1)：105-110. Guo Hu, Luo Jingrun. Mechanical behavior of PBX under cyclic loadings[J]. Explosion and Shock Waves, 2013，33(Suppl 1):105-110.

[3] Rauchs G, Bardon J. Identification of elasto-viscoplastic material parameters by indentation testing and combined finite element modelling and numerical optimization[J]. Finite Elements in Analysis and Design, 2011,47(7):653-667.

[4] Springmann M, Kuna M. Identification of material parameters of the Gurson-Tvergaard-Needleman model by combined experimental and numerical techniques[J]. Computational Materials Science, 2005,32(3/4):544-552.

[5] Chaparro B M, Thuillier S, Menezes L F, et al. Material parameters identification: Gradient-based, genetic and hybrid optimization algorithms[J]. Computational Materials Science, 2008,44(2):339-346.

[6] Majzoobi G H, Dehgolan F R. Determination of the constants of damage models[J]. Procedia Engineering, 2011,10:764-773.

[8] 陈炳瑞，冯夏庭，丁秀丽，等．基于模式-遗传-神经网络的流变参数反演[J]．岩石力学与工程学报，2005，24(4)：553-558． Chen Bingrui, Feng Xiating, Ding Xiuli, et al. Back analysis on rheological parameters based on pattern-genetic-neural network[J]. Chinese Journal of Rock Mechanics and Engineering, 2005,24(4):553-558.

[9] 高军，黄再兴．PBX 炸药粘弹性损伤本构模型的参数识别[J]．工程力学，2013，30(7)：299-304． Gao Jun, Huang Zaixing. Parameter identification for viscoelastic damage constitutive model of PBX[J]. Engineering Mechanics, 2013，30(7)：299-304.

[10] Solano G J, Rodriguez V K, Garcia N D. Model-based spectral estimation of Doppler signals using parallel genetic algorithms[J]. Artificial Intelligence in Medicine, 2000,19(1):75-89.

[11] 李守巨，刘迎曦，孙伟．智能计算与参数反演[M]．北京：科学出版社，2008.

[12] Cantu-Paz E. Designing efficient and accurate parallel genetic algorithms (parallel algorithms)[D]. University of Illinois at Urbana-Champaign, 1999.

[13] Potts J C, Giddens T D, Yadav S B. The development and evaluation of an improved genetic algorithm based on migration and artificial selection[J]. IEEE Transactions on Systems, Man, and Cybernetics, 1994,24(1):73-86.

[14] 刘桂萍．基于微型遗传算法的多目标优化方法及应用研究[D]．长沙：湖南大学，2008.

[15] Sheblé G B, Brittig K. Refined genetic algorithm-economic dispatch example[J]. IEEE Transactions on Power Systems, 1995,10(1):117-124.

[16] Paris P C, Erdogan F. A critical analysis of crack propagation laws[J]. Journal of Fluids Engineering, 1963,85(4):528-533.

(责任编辑 丁 峰)

Application of multiple-population genetic algorithm in parameter identification for PBX constitutive model

Gao Jun1, Huang Zaixing2

(1.ShanghaiCivilAviationCollege,Shanghai200232,China;2.StateKeyLaboratoryofMechanicsandControlofMechanicalStructures,NanjingUniversityofAeronautics&Astronautics,Nanjing210016,Jiangsu,China)

In this work, the standard genetic algorithm (SGA) was parallel processed using multiple parallel structures. Based on the structures, a multiple-population genetic algorithm (MPGA) was established by introducing the immigration operator and quintessence population. Self-adaptive operators of crossover probability and mutation probability were designed to improve the convergence speed of the MPGA. Combining ABAQUS with the improved multiple-population genetic optimized algorithm, a parameter identification method of constitutive model was built. Using the proposed method, a simulation example of parameter identification for PBX viscoelastic damage constitutive model was carried out. Comparison was made between methods based on SGA and MPGA. The results show that the MPGA method can effectively overcome the difficulty of the premature convergence and the identification result is more robust. The method is suitable for the optimization of complex nonlinear systems due to its superiority in the convergence speed and searching ability, and it can be applied to the parameter identification of other models.

solid mechanics; parameter identification; multiple-population genetic algorithm; PBX explosive; ABAQUS; constitutive model

10.11883/1001-1455(2016)06-0861-08

2015-01-27； < class="emphasis_bold">修回日期：2015-05-28

2015-05-28

高 军(1984— )，男，博士研究生,gaojun.nuaa@foxmail.com。

O343;TJ55 <国标学科代码：1301520 class="emphasis_bold"> 国标学科代码：1301520 文献标志码：A国标学科代码：1301520

A