基于代理模型的工程结构可靠性分析

易 平， 杨 潍 宁， 谢 东 赤

( 大连理工大学 建设工程学部, 辽宁 大连 116024 )

0 引 言

工程实践中往往无法得到结构功能与结构随机变量的显式表达式，因而常用的可靠性分析方法不能得到有效运用．基于工程结构有限元软件和计算机随机数生成，计算工程结构可靠度的蒙特卡罗法得到重要运用．但由于随机模拟次数多，计算成本大，蒙特卡罗法计算结构可靠度也有很大的局限性．鉴于这种情况，张哲等引入代理模型方法计算结构可靠度[1]．

代理模型是指通过数学手段构造计算量小，但计算结果与数值分析或物理实验结果相近的近似数学模型，以代替原分析模型．传统二次多项式响应面法(response surface method，RSM)和Kriging方法是目前常用的两种代理模型．Kriging 方法最早于1951年由南非地质学家Krige提出；1989年Sacks等[2]建立了用于计算机试验设计(design and analysis of computer experiments，DACE)的Kriging方法；Lophaven等[3]采用Matlab编程实现了Kriging方法的DACE工具箱．Kriging代理模型在结构领域中多用于结构优化设计[4-5]．王红等[6]基于Kriging代理模型提出了混合优化算法，成功提高了高维参数优化问题的计算效率．

Kaymaz[7]首先将Kriging方法应用于结构可靠性分析，并与传统响应面法做了比较．2011年，陈志英等[8]采用粒子群优化(particle swarm optimization，PSO)算法搜索Kriging近似模型参数的最优值，分析了涡轮盘低循环疲劳的可靠度．2013年，冯欢欢等[9]应用样本重复使用的策略，充分利用迭代过程中的样本点避免了资源信息的浪费，有效地提高了可靠度的计算效率和精度．2015年，Yi等[10]提出了样本选择累积策略即优先选择更接近极限状态曲面的点，并利用该策略基于传统响应面模型和改进Kriging模型进行了边坡稳定可靠度研究．

目前代理模型在结构优化中得到了广泛应用，但其在结构可靠性分析中仍然没有普及．代理模型方法作为一种计算精度高、近似误差小的可靠性分析方法拥有广阔的前景．本文采用基于粒子群优化的Kriging方法(PSO-Kriging)计算功能函数的可靠指标并与传统响应面法(RSM)做对比，同时研究多次拟合近似代理模型过程中的样本全部累积和选择累积策略，最后将基于粒子群优化的Kriging方法运用于实际结构的可靠性分析．

1 结构可靠性理论

结构在规定时间内和规定条件下，完成预定功能的概率称为结构可靠度，该预定功能可以通过如下功能函数表达：

Y=g(x1，x2，…，xn)

(1)

其中x1，x2，…，xn为n个随机变量；Y>0表示结构可靠，Y<0表示结构失效．

工程中一般多用失效概率Pf反映结构可靠度，Pf为失效域内联合概率密度函数的多维积分．当随机变量增多时，通过数值积分求解失效概率非常困难甚至不可行，因此引入可靠指标β的概念．可靠指标β和失效概率Pf通过下式进行转换：

β=-Φ-1(Pf)

(2)

由于结构的复杂性，一般很难写出其显式功能函数，代理模型方法是求解此类问题的一种有效方法．代理模型方法利用有限次数的结构分析结果拟合一个显式表达式近似代替未知的功能函数进行可靠指标计算．

2 代理模型方法

基于代理模型方法进行可靠性分析，首先要通过试验设计方法选取样本点．拉丁超立方抽样(Latin hypercube sampling，LHS)是一种多维分层抽样方法，其基本思想是在试验空间内等概率抽样，将每个变量的抽样空间划分为N个(样本点数)等概率的子空间，从而整个试验空间被划分为Nn个等概率的子空间，在这些子空间中进行随机变量的取值配对，最终得到N个样本点[11]．前期的一些研究工作表明拉丁超立方抽样试验均衡性较好，实现简单，所以本文采用LHS方法选取样本点．利用LHS选取样本点后，分别基于传统响应面模型和Kriging模型进行可靠性分析．

2.1 传统响应面法(RSM)

传统响应面法又称多项式响应面法，是一种应用广泛的代理模型方法．近似模型通常采用不含交叉项的二次多项式形式：

g-(x)=a0+∑ni=1aixi+∑ni=1aiix2i

(3)

a=(a0aiaii)T，为2n+1个待定系数．进行m(m≥2n+1)次独立试验，得到其对应的响应值，组成用于确定系数a的方程组，然后利用最小二乘QR分解算法(least squares QR，LSQR)[12]确定待定系数得到近似函数．

2.2 Kriging方法

Kriging模型由回归项和随机误差项两部分组成：

g-(x)=∑pj=1βjfj(x)+z(x)

(4)

式中：βj为回归系数，fj(x)表示基函数，z(x)为随机误差．随机误差z(x)是Kriging方法与RSM最主要的不同点，以下是随机误差z(x)的主要性质：

E[z(x)]=0； var[z(x)]=σ2cov [z(xi)，z(xj)]=σ2R(xi，xj)

(5)

R为含有参数θ的关于样本点xi和xj的相关函数，通常采用如式(6)所示的高斯函数[13]：

(6)

其中θk为参数向量θ的第k个元素．参数θ的取值影响着Kriging模型的精度，需要获取极大似然意义下的最优相关参数θ，从而确保Kriging预测结果的最优无偏性．

Lophaven等编写了Kriging方法的Matlab程序——DACE工具箱[3]．在DACE工具箱中，利用输入的样本点和响应值得到相应的参数值，进而根据式(4)得到功能函数的显式表达．DACE工具箱采用模式搜索法寻找最优参数θ，这种方法需要给定初始点θ0，θ0对Kriging模型的精度有很大影响，θ0取值不当会使得参数θ陷入局部最优，从而降低Kriging模型的精度．为了避免初值选择对Kriging模型的影响，本文采用粒子群优化(PSO)算法[14]搜索最优参数θ．有别于模式搜索法，PSO算法不需要给定初始点θ0，采用多点并行搜索，种群中各粒子在搜索过程中实现信息共享，每个粒子追随自身最优位置和种群最优位置，通过不断调整自身的坐标和速度在解空间中进行搜索，最终趋近最优解．

为了对比Kriging模型和传统响应面模型的差异，选取二维非线性测试函数，函数表达式如下：

2sin 2x1sin 2x2

变量x1和x2的取值范围为[0.5，3.5]，采用LHS方法选取55个样本点(图1(b)、(c)所示的黑色散点)．分别使用PSO-Kriging和RSM对函数拟合，并画出近似函数图像，如图1所示．

(a) 原函数图像

(b) PSO-Kriging

(c) RSM

图1 原函数、Kriging模型和RSM模型近似曲面对比

Fig.1 Comparison of original function image with the approximate surface of Kriging and RSM models

从图1可以看出，PSO-Kriging较好地拟合了原函数，样本点都落在近似曲面上．而基于RSM的近似函数图像和原函数图像出入很大，样本点也多不在近似曲面上．这表明当拟合非线性程度较高的函数时，Kriging模型的拟合精度要高于传统响应面模型．这也很容易理解，由于RSM中近似功能函数的表达形式完全预设，当真实功能函数的非线性程度增加时，RSM很难构造精度很高的近似函数．

3 样本累积

使用代理模型方法求解结构可靠度时，由于其近似性需多次拟合功能函数，多次求得近似可靠指标直至收敛．首次拟合功能函数时通常选择随机变量的均值点作为抽样中心，利用LHS方法产生2n+1个样本点；后续迭代中则以前一次求得的验算点作为抽样中心选取新的样本点．传统方法是每次迭代生成2n+1个样本点，只利用新生成的样本点构造功能函数．下章算例表明这种摈弃旧样本点信息的方法在求解可靠指标的精度和效率上都不高．

样本累积即样本重复利用，将历史样本点和新的样本点组合在一起共同构造代理模型，计算可靠指标，不断循环直至收敛．其方法流程如图2所示．

图2 基于样本累积的代理模型方法流程图Fig.2 Flow chart of surrogate model based on sample accumulation

样本累积可分为样本全部累积和样本选择累积．样本全部累积是指在第k(k>1)次迭代中，将之前k-1次迭代用到的所有样本点累积到本次迭代产生的2n+1个样本点中，利用k×(2n+1)个样本点构造代理模型计算可靠指标．而样本选择累积是指从上次迭代使用的样本点中选择距离极限状态面最近的l个样本点，即使得响应Y绝对值最小的l个样本点，累积到本次迭代产生的2n+1个样本点中，利用2n+1+l个样本点构造代理模型计算可靠指标．在接下来的算例中分别使用样本不累积、样本全部累积以及样本选择累积进行可靠性分析，并进行比较．

4 算例分析

下面的两个显式算例均使用LHS方法每次迭代选取2n+1个样本点，其中选择样本累积时l=2n+1．近似功能函数则分别采用基于粒子群优化的Kriging方法(PSO-Kriging)和传统响应面法(RSM)构造，最后根据可靠指标β的几何意义[15]利用Matlab优化工具箱求解可靠指标β．同时使用一次二阶矩计算可靠指标，对比计算效率和计算精度．收敛标准均设为βi-βi-1≤0.001．使用蒙特卡罗法计算可靠指标作为参考解，其中蒙特卡罗法所必需的样本点数M通过下式进行估算：

(7)

基于PSO-Kriging和RSM，在迭代过程中采用3种不同的样本处理方法，计算功能函数g(x)的可靠指标，并与蒙特卡罗法计算的参考解相对照，计算结果如表1所示．

表1 算例1中不同方法的计算结果

由表1可以看出，综合考虑3种样本使用情况，在保证计算精度前提下，PSO-Kriging的计算效率要高于RSM．不论是PSO-Kriging还是RSM，迭代过程采用样本累积都提高了计算效率，其中样本选择累积在计算效率上要优于样本全部累积．针对该显式算例，传统的一次二阶矩具有较高的计算效率和计算精度．

表2 算例2中不同方法的计算结果

由表2可以看出，无论是使用PSO-Kriging还是RSM，样本不累积时都不能得到最终的收敛解．若迭代过程中样本全部累积，PSO-Kriging经过5次迭代25次函数计算得到最终的可靠指标值；RSM的计算效率略低，经过11次迭代55次函数计算得到收敛解．对于RSM，样本选择累积在计算效率和计算精度上均要高于样本全部累积，但对于PSO-Kriging，样本选择累积较样本全部累积只稍提高了计算精度．此外样本累积的PSO-Kriging的计算效率和计算精度要高于传统的一次二阶矩．

由以上两个数值算例可知，迭代过程中累积之前的样本点避免了资源的浪费，使迭代更加高效稳定．对于RSM，样本选择累积较样本全部累积提高了计算精度和效率；而对于PSO-Kriging，样本选择累积较样本全部累积的优势不明显．但整体来看，PSO-Kriging明显比RSM效率高，所以接下来的3个结构可靠性分析隐式算例采用PSO-Kriging．

图3 多层框架结构Fig.3 Multi-story frame structure

表3 算例3中随机变量信息

表4 算例3中不同方法的计算结果

算例4一板柱结构如图4所示[18]，其宽度为6 m，柱间距为5 m，层高为3 m．板的厚度为0.2 m，方柱的名义边长为0.5 m，材料弹性模量为200 GPa，泊松比为0.3，密度为7 800 kg/m3．取结构每层方柱实际边长与名义边长的比值ki为随机变量，各个随机变量相互独立，其统计特征值列于表5．进行结构模态分析，当结构基频超过4.107 2 Hz时认为结构失效，求其可靠指标．计算结果如表6所示．

图4 板柱结构Fig.4 The slab-column structure

表5 算例4中随机变量信息

表6 算例4中不同方法的计算结果

算例5图5所示矩形薄板[19]的尺寸为0.254 0 m×0.304 8 m，在位于角点到中点距离的2/7处设置4个支撑．板的边缘和对称线上分别施加Q1=26.27 kN/m和Q2=35.02 kN/m的分布荷载．根据对称性，选择板的1/4建立如图6所示的有限元模型，共49个四节点板单元，所有单元采用密度为2 134 kg/m3的同一种材料．板厚采用变量联接技术(如图6所示，颜色相同的单元的板厚用同一个变量表示)．7个板厚和材料弹性模量是相互独立的随机变量，其统计信息列于表7．功能函数取为板中点的垂直位移要小于0.152 4 mm．计算结果如表8所示．

图5 矩形薄板Fig.5 A rectangular plate

由以上3个工程算例可知，基于PSO-Kriging 代理模型进行结构可靠性分析，若未进行样本累积，有可能迭代过程振荡不能收敛到最终的可靠指标值；采用样本累积能显著改善收敛性能，样本选择累积在计算效率和计算精度上都要优于样本全部累积．在分析结构算例时，常用的一次二阶矩的计算效率较低，且误差较大．

图6 矩形薄板的有限元模型Fig.6 Finite element model for the rectangular plate

表7 算例5中随机变量信息

表8 算例5中不同方法的计算结果

5 结论与展望

(1)RSM难以对非线性程度高的极限状态曲面作出较好的拟合，PSO-Kriging有良好的预测能力，无论是在计算效率还是精度上都要优于RSM．

(2)迭代过程中摒弃之前迭代样本点信息的做法造成了资源的浪费，增加了迭代次数和计算误差，甚至造成迭代过程振荡不收敛．

(3)迭代过程中采用样本累积能显著改善收敛性能，提高收敛速度，得到较为精确的可靠指标值；样本选择累积在计算效率和计算精度上通常都优于样本全部累积，但有时对效率的改善不明显，后续工作应继续研究样本选择方法，尝试利用Kriging方法不仅能提供最优无偏预测，而且能对预测结果误差作定量估计的特性，构造学习函数进行样本选择．

[1] 张 哲,李生勇,滕启杰. 一种改进的结构可靠度分析中响应面法[J]. 大连理工大学学报, 2007,47(1):57-60.

ZHANG Zhe, LI Shengyong, TENG Qijie. An improved response surface method for structural reliability analysis [J].JournalofDalianUniversityofTechnology, 2007,47(1):57-60. (in Chinese)

[2] SACKS J, WELCH W J, MITCHELL T J,etal. Design and analysis of computer experiments [J].StatisticalScience, 1989,4(4):409-423.

[3] LOPHAVEN S N, NIELSEN H B, S∅NDERGAARD J. DACE - a Matlab kriging toolbox, version 2.0: IMM-TR-2002-12 [R]. Copenhagen:Technical University of Denmark, 2002.

[4] SIMPSON T W, MAUERY T M, KORTE J J,etal. Kriging models for global approximation in simulation-based multidisciplinary design optimization [J].AIAAJournal, 2001,39(12):2233-2241.

[5] LAURENCEAU J, SAGAUT P. Building efficient response surfaces of aerodynamic functions with kriging and cokriging [J].AIAAJournal, 2008,46(2):498-507.

[6] 王 红,王希诚,李克秋. 混合Kriging代理模型的高维参数估计优化算法[J]. 大连理工大学学报, 2015,55(2):215-222.

WANG Hong, WANG Xicheng, LI Keqiu. Hybrid Kriging surrogate model optimization algorithm for high-dimension parameter estimation [J].JournalofDalianUniversityofTechnology, 2015,55(2):215-222. (in Chinese)

[7] KAYMAZ I. Application of kriging method to structural reliability problems [J].StructuralSafety, 2005,27(2):133-151.

[8] 陈志英,任 远,白广忱,等. 粒子群优化的Kriging近似模型及其在可靠性分析中的应用[J]. 航空动力学报, 2011,26(7):1522-1531.

CHEN Zhiying, REN Yuan, BAI Guangchen,etal. Particle swarm optimized Kriging approximate model and its application to reliability analysis [J].JournalofAerospacePower, 2011,26(7):1522-1531. (in Chinese)

[9] 冯欢欢,蒋向华. 样本重复使用失效响应曲线分析结构可靠度方法[J]. 航空动力学报, 2013,28(10):2228-2234.

FENG Huanhuan, JIANG Xianghua. Sample points reused failure response curve method for structural reliability analysis [J].JournalofAerospacePower, 2013,28(10):2228-2234. (in Chinese)

[10] YI Ping, WEI Kaitao, KONG Xianjing,etal. Cumulative PSO-Kriging model for slope reliability analysis [J].ProbabilisticEngineeringMechanics, 2015,39(2015):39-45.

[11] MCKAY M D, BECKMAN R J, CONOVER W J. Comparison of three methods for selecting values of input variables in the analysis of output from a computer code [J].Technometrics, 1979,21(2):239-245.

[12] PAIGE C C, SAUNDERS M A. LSQR: An algorithm for sparse linear equations and sparse least squares [J].ACMTransactionsonMathematicalSoftware, 1982,8(1):43-71.

[13] MATHERON G. Principles of geostatistics [J].EconomicGeology, 1963,58(8):1246-1266.

[14] PARSOPOULOS K E, VRAHATIS M N. Recent approaches to global optimization problems through Particle Swarm Optimization [J].NaturalComputing, 2002,1(2/3):235-306.

[15] HASOFER A M, LIND N C. Exact and invariant second-moment code format [J].JournaloftheEngineeringMechanicsDivision, 1974,100(1):111-121.

[16] GONG Jinxin, YI Ping, ZHAO Na. Non-gradient - based algorithm for structural reliability analysis [J].JournalofEngineeringMechanics, 2014,140(6):04014029.

[17] 张 明. 结构可靠度分析——方法与程序[M]. 北京:科学出版社, 2009.

ZHANG Ming.StructuralReliabilityAnalysis-MethodsandProcedures[ M]. Beijing: Science Press, 2009. (in Chinese)

[18] YI Ping, ZHU Zuo. Step length adjustment iterative algorithm for inverse reliability analysis [J].StructuralandMultidisciplinaryOptimization, 2016,54(4):999-1009.

[19] YI Ping, CHENG Gengdong, JIANG Lei. A sequential approximate programming strategy for performance-measure-based probabilistic structural design optimization [J].StructuralSafety, 2008,30(2):91-109.