2. 北京航空航天大学 可靠性与环境工程技术重点实验室, 北京 100083
2. Science and Technology on Reliability and Environment Engineering Laboratory, Beihang University, Beijing 100083, China
当现有的可靠性分析方法应用于复杂的工程结构时，往往面临巨大的挑战[1]，功能函数通常是高度非线性甚至是隐式的，而且需要借助有限元分析(Finite Element Analysis，FEA)进行评估，计算量大，计算时间长[2-3]。一阶可靠度算法(First Order Reliability Methodology，FORM)、二阶可靠度算法(Second Order Reliability Methodology, SORM)已被广泛用于估算结构系统失效的概率[4]，但是只能适用于功能函数表达式已知的情况。蒙特卡罗模拟(Monte Carlo Simulation, MCS)方法是结构可靠性分析中最简单也是最广泛使用的方法，计算精度高，适用于结构可靠性的隐式问题。但由于工程结构的失效概率通常较小，MCS须生成大量样本以确保准确性，所带来的计算成本可能会变得非常高，从而极大地限制了MCS在工程中的应用。
为了降低可靠性分析的计算成本，代理模型越来越多地被用来实现结构可靠性。代理模型作为一个近似模型可以找到输入和输出之间的潜在关系。如果与功能函数足够接近，则代理模型可以完全表示功能函数，以准确评估性能函数的值。常用的方法有响应面方法(Response Surface Methodology，RSM)[5]、支持向量机[6]、神经网络[7]、Kriging模型[8]。其中，Kriging模型是一种方差估计最小的无偏估计的随机过程算法，由于其对非线性函数的良好近似能力和独特的误差估计功能[9]，越来越多地用于可靠度的求解过程。Zhang等[10]运用Kriging模型进行结构可靠度分析，Liu等[11]运用Kriging模型和MCS方法进行齿轮接触可靠性的分析。刘瞻[12]和陈志英[13]等将人工智能优化算法融入到Kriging模型的构造中，进一步提高了预测精度。为了降低复杂结构的可靠性分析的计算成本并实现快速可靠性分析，尽可能减少样本数量，Ma等[14]在采样过程中应用学习函数，从而构造Kriging模型，Li等[15]提出一种基于最可能失效点(Most Probable Point, MPP)的局部抽样方法，提高了Kriging的精度。Echard等[16]提出动态Kriging联合MCS的可靠度计算方法。Xia等[17]提出了一种基于自适应动态泰勒Kriging的高效可靠性优化算法。通过二元粒子群算法优化选择基函数，并确定具有期望拟合误差的最小数量的采样数据。
由于初始训练的高质量样本总是难以获得，本文在训练小样本的基础上，采用修改代理模型拟合误差的动态更新机制。针对复杂结构隐式功能函数的可靠性分析问题，本文将混合粒子群-模拟退火(Particle Swarm Optimization-Simulated Annealing, PSOSA)算法应用到动态Kriging模型中相关参数的寻优过程，并将改进的Kriging模型应用于可靠度的分析计算问题。该方法运用PSOSA算法高效求解最优相关参数，构造Kriging模型；同时基于动态更新机制，逐步加入最佳样本点，从而不断修正Kriging模型的精度，并结合一次二阶矩法进行可靠度计算。
1 动态Kriging模型 1.1 Kriging模型近年来, Kriging模型作为一种新型的响应面模型技术在航天等工程优化领域中得到了广泛应用。对于k个n维样本点及其对应的响应值Y=[G(x1), G(x2), …, G(xk)]T，Kriging模型采用高斯随机过程模型，其插值结果定义为已知样本函数响应值的线性加权[9]，即
|
(1) |
式中：x=[x1, x2, …, xk]T；f(x)=[f1(x), f2(x), …, fp(x)]T为回归多项式基函数向量，p为回归多项式的数量；β=[β1, β2, …, βp]T为多项式参数向量；z(x)为服从正态分布N(0, σ2)的随机过程，σ为标准差，协方差方程为
|
(2) |
式中：R(xi, xj)为样本点中任意2个输入向量xi和xj的空间相关方程，通常采用高斯相关方程，其形式为
|
(3) |
式中：xil为输入向量xi的第l个分量；θl为相关性参数，可通过最大似然估计得到
|
(4) |
式中：R为相关矩阵，其中元素Rij=R(xi, xj)(i, j=1, 2, …, k)，相关参数直接控制着Kriging模型的输入输出特性。
基于给定的样本点，多项式参数向量β与随机过程方差σ2的估计值计算式分别为
|
(5) |
|
(6) |
式中：F为由回归多项式函数值组成的回归矩阵[18]。
对于任意未知点x0，其函数预测值
|
(7) |
|
(8) |
|
(9) |
式中：r0=[R(x0, x1), R(x0, x2), …, R(x0, xk)]T。
1.2 动态更新机制为了提高可靠性分析的效率，动态代理模型应符合以下条件：①初始样本数量小；②用于小样本学习的替代模型拟合；③每迭代一次，对FEA进行一次调用，以修正代理模型。
本文选择动态Kriging模型来减少调用FEA的次数。代理模型的计算精度主要取决于对极限状态面的逼近程度。一般来说，使用一些预设的训练样本很难构建具有高拟合精度的响应曲面，而响应曲面的动态更新可用于校正拟合误差。动态Kriging可靠度具体算法流程如下：
步骤1 建立初始的样本集。通过拉丁超立方抽样生成初始样本点，并调用FEA计算其响应值。
步骤2 构造动态Kriging。首先，通过优化技术寻找最优相关参数；然后，设置最优基函数，并构造代理模型来近似目标函数。
步骤3 利用DACE工具箱可求出模型对各个变量的梯度，结合FORM法计算MPP[12]。判断MPP是否收敛。若是，求解可靠度等指标。若否，进入步骤4。
步骤4 上述MPP作为新的样本点加入初始样本集，继续迭代。
上述过程实现了可靠度计算过程中，Kriging模型在实际极限状态响应面附近拟合误差的自适应修正。减少了FEA的调用次数，这使计算成本大大降低，有利于减少模型对样本的依赖程度，加快了收敛速度。
2 改进的动态Kriging模型 2.1 优化算法 2.1.1 PSO算法假设搜索空间
|
(10) |
|
(11) |
式中：vk为第k个粒子的速度；R1′和R2′为0~1之间的随机数；非负常数c1和c2为学习因子，迭代终止条件一般选为最大迭代次数或粒子群迄今为止搜索到的最优位置满足预定最小适应阈值；系数w为权重因子，实际上控制了搜索空间的搜索，是跳出目前最小点的随机机制，其作用是允许粒子在整个搜索区域内移动，让搜索在局部区域内开展。
本文中系数w由式(12)确定：
|
(12) |
式中：C=c1+c2。
PSO的工作原理如下：首先，粒子群由搜索区域内随机分布的粒子组成，然后，PSO通过更新粒子的速度和位置来迭代寻找最优解。当满足收敛标准时，循环终止。但是，运用公式Pb和Pg来搜索最优解，收敛速度快。然而所有粒子有移向最优经验点的趋势，并在局部区域内搜索，所以可能会导致收敛过早的问题。因此，尽管PSO比其他方法收敛速度快，但结果可能精度不高。
2.1.2 SA算法模拟退火算法起源于固体退火原理，模拟了高温金属降温的热力学过程。在搜索过程中具有概率突跳的能力，能够有效地避免搜索过程中陷入局部极小解。试验点从xi到xi+1的接受概率为
|
(13) |
式中：y为目标函数，Δy=y(xi+1)－y(xi)和T为控制参数，可类比为温度，温度的降低实际上控制了接受概率。
设φj为每个变量接受的移动次数与试验次数的比值，则扰动可写为
|
(14) |
式中：R′(－1, 1)为-1~1之间的随机数；V为步长向量。Vbi是变量上下界之差的一半，可表示为
|
(15) |
式中：Vi为V中的元素。
理论上证明，当采用非常缓慢的温度下降率时，找到全局最小值的概率可以接近1。
最常见的温度下降规律为
|
(16) |
式中：T0为初始温度；RT为0.8~0.999 9之间的常数。
2.1.3 混合PSOSA算法在混合PSOSA算法中，搜索最优解时，速度矢量控制粒子在搜索空间内的移动方式。速度矢量有3部分组成：
1) 惯性项(wvki)，使得粒子保持原来的速度向量，以防粒子的速度急剧变化。
2) 认知行为(c1R1′(Pkb－Xki))，包含粒子经历过的最优位置。
3) 社会行为(c2R2′(Pg－Xki))，给予粒子朝群体最优位置运动的趋势。
本文PSOSA算法的主要思想是通过更新群体的群体最佳位置来改善群体的社会行为。当PSO循环中最优解没有改进时，旧的群体最佳位置将被替换为使用SA算法计算的新位置。新的最佳位置实际上向群体的社会领导者(即前一个全局最优解所属的粒子)发出信号，用于更新其方向。PSOSA算法充分结合了PSO算法全局搜索能力和SA算法的局部搜索能力，能达到很好的收敛结果。算法示意图如图 1所示, 具体算法如下：
|
| 图 1 PSOSA算法示意图 Fig. 1 Sketch map of PSOSA algorithm |
步骤1 初始化。m为粒子群粒子的数；imax为允许最大迭代次数；w为权重因子；c1和c2为学习因子；vmax为速度阈值，T0为初始温度；RT为温度降低参数。
步骤2 粒子的位置。可行域均匀分布m个粒子；确定每个粒子的适应度值。
步骤3 初始最佳位置的确定。确定每个粒子Pkb和群体的最佳位置Pg。
步骤4 更新每个粒子的速度和位置：
|
步骤5 更新每个粒子的Pkb和群体的Pg：若f(Xki+1)≤f(Pkb)，则Pkb=f(Xki+1)；若f(Xki+1)≤f(Pg)，则Pg=f(Xki+1)。
步骤6 判断最优解是否有所改进。如果第i+1次迭代的Pg没有优于第i次迭代的Pg，则进入步骤7；否则，更新温度Ti+1=TiRT，返回步骤4。
步骤7 SA算法。计算一个测试位置，Ptrialg=Pg+R′(－1, 1)V，R′取随机数。如果f(Ptrialg)≤f(Pg)，则Pg=Ptrialg，更新温度，Ti+1=TiRT并继续迭代；如果f(Ptrialg)>f(Pg)，检测Metropolis Monte Carlo接受概率，如果e－Δf/T≥R′(0, 1)，则Pg=Ptrialg并更新温度；如果e－Δf/T < R′(0, 1)，取目前的Pg为最优解。返回步骤4。
步骤8 重复整个过程直到收敛。
该算法主要有以下几个特点：
1) 通过改进速度矢量的组成促使群体向整体最优位置所在的搜索空间移动，这加快了收敛过程。
2) PSO算法和SA算法的结合提高了整个寻优机制的质量和稳定性。
3) 在SA迭代中PSO速度公式中的认知成分得到了保留。
2.2 Kriging建模流程目前Kriging建模过程普遍采用的模式搜索法对初始点十分敏感，容易造成求解θ*无法收敛，Kriging预测精度很低。主要原因是模式搜索所采用的单点序列搜索方式决定了它只能从一个指定的起点出发，并收敛于一个距起点较近的局部最优解。为了寻找最优θ*，只能将起点接近于θ*，然而这很难确定。不仅如此，模式搜索法搜索路径单一，不能根据似然函数的变化而变化，因此即使根据经验信息，使起点接近θ*，搜索路径也有可能偏离导致无法收敛[13]。
本文一方面运用PSOSA算法进行搜索，种群中各粒子在搜索最优解的过程中实现信息共享，每个粒子不断根据个体极值与群体极值的变化情况动态调整着自己的搜索方向，另外，当PSO循环没有改善，原来全局最优位置被通过SA计算得到的新位置代替，从而在缺乏先验信息的前提下，也能搜索到最优解θ*。另一方面, 通过自适应动态更新机制减少调用FEA的次数，从而提高计算效率。
其中，运用PSOSA算法搜索极大似然意义下的最优相关参数时：
1) 似然函数作为适应度函数。
2) 为转化为无约束优化，粒子采用对数形式。
综上所述，本文的总体计算流程如图 2所示。
|
| 图 2 改进的动态Kriging模型流程图 Fig. 2 Flowchart of improved dynamic Kriging model |
假设某结构的极限状态函数[12]为
|
(19) |
式中：x1和x2的分布类型为正态分布，x1~N(38, 3.8)，x2~N(54, 2.7)，并且x1和x2相互独立，具体分布参数见表 1。本案例中，c1=1.5，c2=1.5，RT=0.8，w=2。
| 方法 | 样本点 | βr | 失效概率/10-3 |
| MCS | 108 | 6.3 | |
| RSM | 65 | 2.392 7 | 8.3 |
| 经典Kriging | 40 | 2.475 1 | 6.7 |
| PSO-Kriging | 40 | 2.480 5 | 6.56 |
| 本文 | 40 | 2.489 3 | 6.4 |
表 1中，βr为可靠度系数，βr越大，可靠度越大。
对比MCS、RSM、经典Kriging和PSO-Kriging[13]方法，由表 1结果可知，MCS方法选取样本点108，计算失效概率结果为6.3×10－3；RSM方法处理非线性问题有一定的局限性，计算结果误差较大；经典Kriging方法、PSO-Kriging方法与本文算法选取同样数量的样本数量，失效概率与MCS方法接近，但是经典Kriging方法与MCS方法相比，误差为4.76%，PSO-Kriging方法误差为4.13%，本文算法误差为1.58%。由此看出本文算法精度相比于其他2种方法有所提高。
将式(19)转化为标准正态分布变量空间下的功能函数，即
|
(20) |
式中：u1和u2相互独立，均服从标准正态分布。转换后的极限状态函数G(u1, u2)=0为一曲线，G(u1, u2)>0一侧为安全域，G(u1, u2) < 0一侧为失效域，如图 3所示。根据本文所提算法，首先在给定随机变量空间内，生成初始设计点，建立Kriging模型，并通过后续的动态更新机制逐渐增加样本点。初始设计点和拟合过程的样本点如图 4所示。由图可知，通过更新机制可以有效选择功能函数附近的样本点，可以充分利用少量样本点的信息拟合待求功能函数。
|
| 图 3 极限状态函数 Fig. 3 Limit state function |
|
| 图 4 算法说明图 Fig. 4 Illustration of algorithm |
极限状态函数为G(x)=x13+x23-4，x1和x2均服从正态分布且相互独立，此极限状态函数高度非线性。x1~N(3, 1), x2~N(2.9, 1)。c1=1.2，c2=1.2，RT=0.8，w=2。
将本文计算结果与其他方法进行对比，结果如表 2所示。表 2说明RSM在处理高度非线性问题时无法准确收敛和获得正确结果。与其他方法相比，本文算法的结果与MCS结果最接近。其中MCS计算时间长达60.68 s，本文算法需时12.18 s。
| 方法 | 函数调用次数 | βr | 失效概率/10-3 |
| MCS | 106 | 4.16 | |
| RSM | 5 | 1.472 9 | 70.39 |
| 经典Kriging | 30 | 2.624 5 | 4.30 |
| PSO-Kriging | 30 | 2.630 7 | 4.26 |
| 本文 | 30 | 2.639 6 | 4.15 |
3.3 工程案例
涡轮盘是现代航空发动机的核心部件，它在发动机燃烧室内受到高温燃气的推动，将燃气的热能转化为机械能，驱动发动机的运转。其可靠性水平的高低直接影响发动机的技术水平。
本文选择某航空发动机涡轮盘为例，盘上销钉沿周向均匀分布。参数不确定性的主要来源主要包括载荷、材料的随机性。本文取转速均值为9 550 r/min，销轴的材料为3Cr13，轮盘的材料为TC11，叶片的载荷施加在销轴上，合力均值为24 925 N。具体参数信息见表 3。
| 参数 | 涡轮盘的转速ω | 弹性模量(轮盘)E1 | 泊松比(轮盘)ε1 | 密度(轮盘)ρ1 | 弹性模量(销轴)E2 | 泊松比(销轴)ε2 | 密度(销轴)ρ2 | 均布载荷P |
| 均值 | 9 550 r/min | 123 GPa | 0.33 | 4.48 g/cm3 | 219 GPa | 0.3 | 7.76 g/cm3 | 24 925 N |
| 变异系数 | 0.1 | 0.015 | 0.01 | 0.02 | 0.015 | 0.01 | 0.002 | 0.1 |
涡轮盘的可靠性水平的高低主要取决于结构危险点处的应力应变分布。由于结构形状、载荷为完全对称的，本文选取轮盘的1/37为研究模型。首先以各变量的均值为设置参数，借助ANSYS 18.1进行仿真计算，得到结果如图 5所示，轮盘与销轴交界处应力应变水平最高，是结构的危险点，其中最大应力为0.834 GPa，最大应变为0.010 17。
|
| 图 5 应力和应变计算结果 Fig. 5 Calculation results of stress and strain |
当危险点的最大应力超过允许值Rmax=0.97 GPa，则结构发生失效，相应的功能函数可表示为
|
(21) |
上述功能函数为隐式函数。
采用本文所提算法进行可靠度计算，首先建立Kriging模型，其中相关参数θ*的迭代过程见图 6。
|
| 图 6 θ*的迭代过程 Fig. 6 Iteration process of θ* |
由图 6可以看出，随着迭代次数的增加，本文所提算法在第15次迭代时基本寻得最优相关参数θ*。
利用本文算法与ANSYS 18.1、Isight进行联合仿真计算，最终求得可靠度。为了验证本文算法的有效性和工程实用性，将其与MCS方法进行对比，结果如表 4所示。
| 方法 | 函数调用次数 | 样本点 | 失效概率/10-3 |
| MCS方法 | 105 | 3.300 | |
| 本文算法 | 15 | 40 | 3.352 |
由表 4可知，本文算法的结果与MCS结果十分接近，误差仅为1.576%，验证了本文算法的有效性以及准确性，并适用于工程实践。
4 结论1) 将PSOSA算法引入Kriging的构造过程，PSO算法的全局搜索能力及SA算法跳出局部极小值的能力相结合，克服了经典Kriging过程的局限性和对初始样本点的依赖，保证了极大似然意义下的最优相关参数的求解精度，使寻优过程更加高效和精确，从而有效保证了Kriging模型的预测精度。
2) 本文同时通过动态更新机制，尽可能地减少样本点和调用有限元等数值计算的次数，能够较好地解决功能函数隐式和高度非线性的问题，通过数值案例分析，相比于其他可靠度算法，能够有效提高可靠度计算的效率和精度。通过工程案例分析，本文所提算法可应用于工程实际问题，尤其是功能函数隐式且复杂的问题。有一定的工程实用价值。
| [1] |
DITLEVSEN O, MADSEN H O. Structural reliability methods[M]. Chichester: John Wiley & Sons Ltd, 2007.
|
| [2] |
PELLISSETTI M F, SCHU LLER G I. On general purpose software in structural reliability-An overview[J]. Structural Safety, 2006, 28(1-2): 3-16. DOI:10.1016/j.strusafe.2005.03.004 |
| [3] |
REH S, BELEY J D, MUKHERJEE S, et al. Probabilistic finite element analysis using ANSYS[J]. Steel Construction, 2008, 28(1): 17-43. |
| [4] |
SCHU LLER G I, PRADLWARTER H J, KOUTSOURELAKIS P S. A critical appraisal of reliability estimation procedures for high dimensions[J]. Probabilistic Engineering Mechanics, 2004, 19(4): 463-474. DOI:10.1016/j.probengmech.2004.05.004 |
| [5] |
ZHAO W, QIU Z. An efficient response surface method and its application to structural reliability and reliability-based optimization[J]. Finite Elements in Analysis & Design, 2013, 67(5): 34-42. |
| [6] |
CHERKASSKY V, MA Y. Practical selection of SVM parameters and noise estimation for SVM regression[J]. Neural Networks the Official Journal of the International Neural Network Society, 2004, 17(1): 113-126. DOI:10.1016/S0893-6080(03)00169-2 |
| [7] |
GOMES H M, AWRUCH A M. Comparison of response surface and neural network with other methods for structural reliability analysis[J]. Structural Safety, 2004, 26(1): 49-67. DOI:10.1016/S0167-4730(03)00022-5 |
| [8] |
MATHERON G. Principles of geostatistics[J]. Economic Geology, 1963, 58(8): 1246-1266. DOI:10.2113/gsecongeo.58.8.1246 |
| [9] |
韩忠华. Kriging模型及代理优化算法研究进展[J]. 航空学报, 2016, 37(11): 3197-3225. HAN Z H. Kriging surrogate model and its application to design optimization:A review of recent progress[J]. Acta Aeronautica et Astronautica Sinica, 2016, 37(11): 3197-3225. (in Chinese) |
| [10] |
ZHANG J, ZHANG L M, TANG W H. Kriging numerical models for geotechnical reliability analysis[J]. Soils and Foundations, 2011, 51(6): 1169-1177. DOI:10.3208/sandf.51.1169 |
| [11] |
LIU H, CUI W M, HE B, et al.The contact reliability analysis of gear rack based on the Kriging and RSM method[C]//Proceedings of the 2015 International Conference on Applied Science and Engineering Innovation.Paris: Atlantis Press, 2015: 279-284.
|
| [12] |
刘瞻, 张建国, 王灿灿, 等. 基于优化Kriging模型和重要抽样法的结构可靠度混合算法[J]. 航空学报, 2013, 34(6): 1347-1355. LIU Z, ZHANG J G, WANG C C, et al. Hybrid structure reliability method combining optimized Kriging model and importance sampling[J]. Acta Aeronautica et Astronautica Sinica, 2013, 34(6): 1347-1355. (in Chinese) |
| [13] |
陈志英, 任远, 白广忱, 等. 粒子群优化的Kriging近似模型及其在可靠性分析中的应用[J]. 航空动力学报, 2011, 26(7): 1522-1530. CHEN Z Y, REN Y, BAI G C, et al. Particle swarm optimized Kriging approximate model and its application to reliability analysis[J]. Journal of Aerospace Power, 2011, 26(7): 1522-1530. (in Chinese) |
| [14] |
MA J G, REN Z Y, ZHAO G X, et al. A new reliability analysis method combining adaptive Kriging with weight index Monte Carlo simulation[J]. IEEE Transactions on Magnetics, 2018, 54(3): 7001904. |
| [15] |
LI X K, QIU H B, CHEN Z Z, et al. A local Kriging approximation method using MPP for reliability-based design optimization[J]. Computers & Structures, 2016, 162: 102-115. |
| [16] |
ECHARD B, GAYTON N, LEMAIRE M. A Kriging improvement of the Monte Carlo simulation:AK-MCS method[M]. Boca Raton: CRC Press-Taylor & Francis Group, 2011: 679-686.
|
| [17] |
XIA B, REN Z Y, KOH C S. A novel reliability-based optimal design of electromagnetic devices based on adaptive dynamic Taylor Kriging[J]. IEEE Transactions on Magnetics, 2017, 53(6): 7201504. |
| [18] |
佟操, 孙志礼, 杨丽, 等. 一种基于Kriging和Monte Carlo的主动学习可靠度算法[J]. 航空学报, 2015, 36(9): 2992-3001. TONG C, SUN Z L, YANG L, et al. An active learning reliability method based on Kriging and Monte Carlo[J]. Acta Aeronautica et Astronautica Sinica, 2015, 36(9): 2992-3001. (in Chinese) |