基于Kriging模型和模拟退火粒子群算法的结构有限元模型修正

康俊涛 柯志涵 胡 佳

(武汉理工大学土木工程与建筑学院 武汉 430070)

0 引 言

工程中，实际结构现场试验响应值与建立的结构有限元分析模型计算值往往存在误差，需要对建立的初始有限元模型进行修正.模型修正在结构的性能评价[1]、健康监测[2]、损伤识别[3],以及验证结构设计等领域中广泛应用.有限元模型修正，是建立的初始有限元分析模型利用结构现场试验响应值来修正模型计算值，使得修正后模型计算值与结构试验值误差最小.直接调用有限元软件进行模型修正存在计算效率低，不易于二次开发等缺点[4-5].采用代理模型结合智能算法进行模型修正，克服了修正过程中大量调用有限元模型计算的不足.代理模型是利用显式函数式拟合构造出结构有限元模型设计参数与其响应之间复杂的隐式关系的一种数学模型，其避免了复杂的有限元计算.任伟新等[6]在模型修正中采用代理模型，提高了修正的效率.

常用的代理模型有响应面方法[7](response surface method,RSM)、径向基函数[8](radical basis function,RBF)、Kriging模型[9]等.其中Kriging模型由于存在随机过程，更适合于复杂非线性的结构输入、输出模型，具有非常高的精度.

模型修正问题归根结底可以看成优化问题，常结合优化算法求解.基本粒子群算法(particle swarm optimization，PSO)的思想源于模拟鸟群的觅食行为，是由Kennedy等[10]提出的一种群体协作随机搜索算法.其实现容易、概念简单、鲁棒性强、收敛速度快，在结构模型修正[11]、优化设计[12]等领域得到广泛应用.但PSO在搜索过程中容易陷入局部极值解，存在早熟等缺点[13].而模拟退火粒子群算法[14](simulated annealing particle swarm optimization，SimuAPSO)在搜索过程中受温度参数的控制而避免陷入局部极值解，能够较快的收敛到全局最优解,提高计算效率.

1 Kriging函数模型及试验设计

1.1 Kriging函数模型

Kriging函数模型是由回归模型与随机过程相加而成.其形式为

y(X)=fT(X)β+z(X)

(1)

式中：X为维度为d的设计参数[x1,x2,…,xd]；y(X)为模型响应值；fT(X)β为回归模型，fT(X)=[f1(X),f2(X),…,fm(X)]，β=[β1,β2,…,βm]T为回归矩阵，常用二次无交叉项多项式模型为

(2)

式中：z(X)为均值为0协方差为σ2的随机过程，提供模拟的局部近似.z(X)中任意两个参数的协方差表示为

Cov(z(X)z(ω))=σ2R(X,ω)

(3)

式中：R(X,ω)为任意两个样本点X和ω之间的变异函数,其一般形式为

(4)

式中：Rj(dj)为变异函数的核函数，通常采用高斯函数,则变异的函数形式为

(5)

式中：θj为计算空间相关函数中大于0的系数；n为变异的维数.假设有N个样本点，响应阵和多项式阵分别为Y=[y(X1)…y(XN)]，F=[fT(X1)…fT(XN)]T，得到β和σ2值为

(6)

(7)

式中：空间相关矩阵为

(8)

(9)

式中：|R|为矩阵R的行列式值；θ值可以通过遗传算法、模式搜索、粒子群算法等优化方法得到.确定了θ1,θ2,…,θd的值，最后利用已得到的模型系数和变异函数参数，得到最好线性无偏估计值为

(10)

1.2 试验设计

试验设计是构建Kriging函数的基础，对修正参数进行试验设计抽样得到样本点作为输入集，并将样本点代入有限元计算得到结构响应作为输出集，根据输入和输出集拟合得到Kriging函数.因此试验设计方法直接影响构造Kriging模型的精度，同时也是决定结构有限元计算量的主要因素.

本文样本抽样方法选用的是由Timothy 等[15]提出的一种基于分层抽样思想的拉丁超立方设计(LHD).LHD将变量空间分成n个子区间，每个子区间概率相同且互不重叠,然后分别独立随机地在每个子区间内采样一次，以保证样本点的均匀性.拉丁超立方抽样的样本为

(11)

式中：n为样本点数；i为水平数；j为维数；dij为1到n的独立随机排列数；eij为[0,1]之间的随机数.

2 模拟退火粒子群算法

模拟退火粒子群算法是一种源于固体退火原理的启发式搜索算法，其步骤如下.

步骤1随机初始化各粒子的位置与速度.

步骤2评价每个微粒的适应度，将第i个微粒所经历的历史最好点储存在pi中，将所有pi中的最优值储存在pg中.

(12)

(13)

vi(t+1)=φ{vi(t)+c1r1[pi-xi(t)]+

c2r2[pg-xi(t)]}

(14)

xi(t+1)=xi(t)+vi(t+1)

(15)

(16)

式中：vi(t)为粒子温度t时的速度；xi(t)为粒子温度t时的位置；ω为惯性权重；pi为每个粒子经过的历史最优点；pg为所有粒子经过的全局最优点；c1、c2为学习因子；r1,r2为扰动因子，r1,r2∈[0,1].

步骤6计算各微粒新的目标值,更新各微粒的pi值及群体的pg值.

步骤7进行退温操作，tk+1=λtk，其中0<λ<1.

步骤8若满足运算精度或迭代次数等停止条件时，输出当前最优值，搜索结束，否则返回步骤4继续搜索.

3 桁架数值模拟算例

通过一空间桁架进行数值模拟，对设计参数人为给定摄动量，构造Kriging函数模型代替空间桁架Midas结构，结合SimuAPSO进行寻优，完成对桁架结构的修正.

具体修正流程图见图1.

图1 模型修正流程图

3.1 桁架模型

利用一个空间桁架结构作为算例验证构建Kriging法用于模型修正的可行性，见图2.桁架竖直和水平方向每根杆单元长度均为0.5 m，总长度为2 m、宽度为0.5 m、高度0.5 m，弹性模量为2.06×105MPa，采用管形截面,外径为16 mm，管壁厚2.5 mm，截面面积为7.85×10-4m2.桁架材料的初始密度(用ρ表示)为7 852 kg/m3,用Midas建立空间桁架模型见图3.

图2 空间桁架结构

图3 Midas建立的空间桁架模型

桁架由34个杆件单元和16个节点组成，节点1约束住三个平动自由度，节点5，6约束住z方向平动自由度，节点10约束住y，z方向平动自由度，且该四个节点均约束住绕x，z轴方向的转动自由度.对34根杆件的密度参数进行扰动，代入有限元模型计算其响应，经过敏感性分析，选取ρ1，ρ5，ρ9，ρ20四个敏感参数，人为地扰动这四个参数数值，ρ1增加20%至9 420 kg/m3，ρ5增加15%至9 028 kg/m3，ρ9增加30%至10 205 kg/m3，ρ20增加10%至8 635 kg/m3.未损伤模型代表初始模型，扰动后的损伤模型代表实际模型，经过Midas空间桁架模型计算得到实际模型与初始模型的前6阶频率见表1.

表1 初始模型和真实模型的前6阶频率

3.2 试验设计

采用LHD抽样方法，对选取的ρ1，ρ5，ρ9，ρ20这四个参数进行抽样40次得到试验点并代入有限元模型中计算得到对应的频率数据.设计参数的变化范围为±50%，即抽样点密度下界为3 925 kg/m3，上界为11 775 kg/m3.将ρ1，ρ9抽样的样本散点图绘于图4，试验点均匀地散布于输入参数空间中.

图4 利用LHD得到的ρ1，ρ9样本散点图

3.3 构建Kriging模型及精度检验

LHD抽样得到样本点作输入集，利用Midas计算得到的前六阶频率作输出集，构建Kriging模型.选取杆件单元1和9的密度参数作为自变量，杆件单元5和20的参数值为初始值不变，做Kriging函数拟合图，见图5.由图5可知,构造的Kriging函数响应面平滑，图6是构建Kriging模型的均方误差图，其MSE值数量级在10-4，表明拟合构建的Kriging函数设计参数估计值与真值误差较小，精确度较高，能代替有限元模型.

图5 频率f1和密度ρ1，ρ9之间的函数拟合图

图6 Kriging函数的均方误差图

3.4 构建目标函数

以构建的Kriging代理模型得到的前六阶频率与实际结构的前六阶频率相对误差的平方和作为目标函数.

(17)

式中：fkrigi为构建的Kriging代理模型得到的第i阶频率响应；fri为实际损伤结构的第i阶频率响应.

3.5 SimuAPSO优化

SimuAPSO算法参数设置如下：初始粒子总数为100，迭代次数200次，粒子维数四维，个体认知部分权重与社会认知部分权重均取2.05.优化结果及误差见图7.由图7可知，四个参数的初始值为7 850 kg/m3,经模拟退火粒子群优化得到的修正值分别为9 807，8 521,10 147,8 901 kg/m3,损伤结构的真实值分别为9 420,9 028，10 205,8 635 kg/m3.图7中条形图上数字表示修正值与真实值的相对误差，可见最大相对误差为5.61%，最小相对误差为0.57%，满足误差要求.

图7 模拟退火粒子群算法参数修正图

图 8为模拟退火粒子群算法和经典粒子群算法对目标函数的最优解进化的曲线图.

图8 两种算法对目标函数的最优解进化曲线图

由图8可知，两种算法修正后的目标函数值基本一致.但SimuAPSO优化效率较高，可以避免陷入局部极值解，更快的收敛到全局最优解.标准粒子群算法得到最优解的迭代次数需要20次以上，而本文采取的模拟退火粒子群将迭代次数缩小到10次左右，提高了优化效率.

4 结 论

1) 基于Kriging代理模型的结构有限元模型修正避免了重复调用有限元进行计算，缩短了运算时间.

2) 基于MATLAB编程的Kriging模型修正过程易于与优化算法相结合，采用模拟退火粒子群算法，较快收敛到全局最优解，显著提高优化效率.

3) 通过空间桁架算例表明，基于Kriging代理模型，结合模拟退火粒子群算法实现有限元模型修正是一个行之有效的方法.在结构优化设计、模型修正与确认、损伤识别等工程领域有很好的应用价值.