曹汝男, 孙志礼, 郭凡逸, 王健
东北大学 机械工程与自动化学院, 辽宁 沈阳 110819
收稿日期：2020-07-06
基金项目：国家自然科学基金资助项目(51775097); 中央高校基本科研业务费专项资金资助项目(N180303031)。
作者简介：曹汝男(1992-)，男，辽宁锦州人, 东北大学博士研究生;
孙志礼(1957-)，男，山东巨野人，东北大学教授，博士生导师。

摘要：首次穿越法可以得到失效概率随时间的变化, 但是计算复杂提高计算成本.拟静态法计算量较小但是只能得到给定时间段对应的失效概率值, 不能得到其变化趋势.针对以上问题, 本文结合Kriging模型和Monte Carlo方法提出一种动态可靠性分析方法AK-MCS-T.AK-MCS-T即拥有拟静态法计算量小的优点又可以给出失效概率随时间变化情况.为提高Kriging模型精度, AK-MCS-T提出一种新的选点准则.新的选点准则不仅考虑到样本点处基于Kriging模型的响应估值符号判断错误的概率, 还考虑了该样本点对应的概率密度.同时, 为保证方法计算结果的准确性, AK-MCS-T提出了新的学习停止条件.最后通过腐蚀梁结构和悬臂管结构证明了新方法的高效性和准确性.
关键词：Monte Carlo方法Kriging模型选点准则学习停止条件动态可靠性分析
Time-Dependent Reliability Algorithm Based on Kriging and Monte Carlo
CAO Ru-nan, SUN Zhi-li, GUO Fan-yi, WANG Jian
School of Mechanical Engineering & Automation, Northeastern University, Shenyang 110819, China
Corresponding author: CAO Ru-nan, E-mail: Caorunan123@163.com.

Abstract: The first outcrossing method can get the change of failure probability with time, but the calculation is complex and may call for more cost. The quasi-static method requires a smaller amount of calculation, but it can only get the failure probability value corresponding to a given period of time, and cannot obtain its change trend.In view of the above problems, AK-MCS-T, a time-dependent reliability analysis method is proposed by combining the Kriging model and the Monte Carlo method. It not only has the advantage of small calculation amount of the quasi-static method, but also can give the change of failure probability with time. AK-MCS-T proposes a new selection criterion for the Kriging model which not only considers the probability of wrong judgment of response evaluation value based on the Kriging model at a sample point, but also considers the probability density corresponding to this sample point. At the same time, a new learning stopping condition is proposed to ensure the accuracy of AK-MCS-T. Finally, the effectiveness and accuracy of the new method are proved by two examples: acorroded beam structure and a cantilever tube structure.
Key words: Monte Carlo methodKriging modelselection criterionlearning stopping conditiontime-dependent reliability analysis
近年来, 机械结构的动态可靠性问题受到广泛关注.顾名思义, 动态可靠性指的是机械结构在给定时间区间[t₀, t_e]内能够完成规定功能而不失效的概率.相较于传统的静态可靠性分析, 动态可靠性的影响因素包含两个部分: ①确定性随机变量，如零部件的初始尺寸、安装误差等不受时间影响的因素；②随机过程变量, 如受到腐蚀、磨损导致的尺寸误差, 循环应力作用下的零件剩余疲劳强度等.动态可靠性分析方法主要包括: 动态Monte Carlo法、首次穿越法和拟静态法.
动态Monte Carlo法^[1]的原理跟经典的Monte Carlo法相同, 即基于大数定律, 当抽取样本数足够大时失效发生的频率近似代替失效概率.但由于Monte Carlo法需要对大量样本点的响应进行计算, 导致计算成本高, 因此通常作为新方法结果的检验标准, 或者配合其他代理模型法使用, 如响应面^[2], Kriging等^[3-4].
相较于动态Monte Carlo法, 基于Rice^[5]公式的首次穿越法在动态可靠性分析中应用广泛.此类方法重点研究功能函数首次超过给定安全阈值的上限或下限(穿越率)的情况, 并基于穿越率进行失效概率的计算.Andrieu等^[6]提出的PHI2可靠性分析方法通过采用一次二阶矩推导各时间点对应的穿越率; Mourefzadeh等^[7]利用并行系统可靠性公式的基本思想来计算过程连续情况下的穿越率; 王正等^[8]推导出在一定假设条件下承受随机载荷作用的机械零件动态可靠性模型.首次穿越法的优势在于可以近似得到机械结构寿命的概率密度函数, 但是计算量大, 效率较低, 这也是首次穿越法急需解决的问题.
拟静态法又称为极值法, 通过计算功能函数随时间变化的极值将动态可靠性问题转换成静态问题.Wang等^[9]使用全局优化EGO算法计算给定随机变量的响应值随时间变化的极小值, 再建立功能函数极小值关于随机变量的函数, 从而计算失效概率; Hu等^[10]对文献[9]进行了改进, 提出了基于Mix-EGO的分析方法, 提高了计算效率.以上两种算法均属于双环结构, 即内环计算响应值的极值, 外环构建关于极值响应模型.为了简化结构, Hu等^[11]提出单环结构(SILK方法), 摒弃了寻找极值点的过程, 大幅度提高计算效率.然而SILK等方法只能计算时间区域[t₀, t_e]对应的失效概率, 不能得到其他子区间的失效概率情况.因此本文提出一种基于Monte Carlo和Kriging的单环动态可靠性分析方法, 该方法不仅可以计算给定时间区域的失效概率, 还可以给出失效概率随时间的变化情况.
1 理论基础1.1 动态Monte Carlo法已知m维随机变量x=[x₁, x₂, …, x_m]^T和n维随机过程v(t)=[v₁(t), v₂(t), …, v_n(t)]^T, 功能函数G(x, v(t), t)在时间区域[t₀, t_e]的失效概率可以表示为

(1)

动态Monte Carlo方法的基本思路是根据随机抽样的样本推导出母体的某些统计概率特性.首先将时间区间进行离散化处理得到K个时间段(对应K+1个时刻), 根据变量的分布信息抽取N_MC个随机样本点构成如下样本集:

(2)

各点对应的功能函数值为
对应区间[t₀, t_e]内的失效概率为

(4)

1.2 AK-MCS法Echard等^[12]结合Kriging模型和Monte Carlo法提出了一种主动学习方法(AK-MCS法), 并证明了该方法的高效性和准确性.在AK-MCS法中, 为提高效率, 使用Kriging模型替代了真实功能函数.
Kriging模型是一种半参数化的差值方法, 通过部分已知的样本点[x₁, x₂, …, x_N]及其对应的真实响应[y₁, y₂, …, y_N]对真实功能函数进行拟合, 并估算出其他样本点x对应的响应均值和标准差.其真实响应值G(x)满足正态分布:

(5)

AK-MCS法的核心思想是通过不断学习, 使用有限个样本点建立精确的Kriging模型代替真实功能函数, 并应用MCS法对Kriging模型进行可靠性分析.针对如何高效选择构成Kriging模型的最佳样本, AK-MCS法中定义了U函数.

(6)

AK-MCS法属于静态可靠性分析方法, 基于Kriging的失效概率表示为

(7)

其中, 是指示函数.当时, 否则.
根据式(7)可知, 失效概率的取值仅与的符号正负有关而与的大小无关, U(x)值越小, 对应的x符号判断错误的概率也就越大.当U(x)>2时, 对应的符号判断正确的概率在0.977以上, 认为是可信的.所以AK-MCS中选取U_min(x)对应的点, 用于构造Kriging模型.当U_min(x)>2, 模型精度满足要求, 学习停止.
2 动态可靠性分析方法2.1 基本原理本文基于AK-MCS法的思想提出一种动态可靠性分析方法, 简称为AK-MCS-T.该方法的主要思路如下:
为阐述方便以功能函数y=G(x, t)为例, 首先使用MCS法生成N_MC个随机样本点[x₁, x₂, …, x_NMC]^T, 同时将时间区间[t₀, t_e]等分为K个时间段, 最终得到N_MC×(K+1)个样本点, 记为样本集S.(注: 如果功能函数包含随机过程可参考文献[11]进行离散化处理.)

(8)

这些样本点用于失效概率的计算以及候补样本用于提高Kriging模型精度.基于Kriging模型计算样本集S中样本点对应的响应值和标准差:

(9)

(10)

根据式(3)、式(4), 可以得到基于Kriging模型在子区间[t₀, t_k]内失效概率估计值：

(11)

(12)

其中k≤K.
2.2 最佳样本选点准则当Kriging模型精度没有满足要求时, 会导致失效概率的误差, 本文方法结合U函数和样本点的分布信息提出一种新的动态可靠性选点准则, 根据式(9)~式(10)样本点x_i在t时刻对应的U函数为, 而该点的联合概率密度函数为f(x_i).Φ(－U(x_i, t))越大说明对应的符号判断错误的概率越大, f(x_i)越大说明该样本点出现的频率越高.因此新的选点准则表示为

(13)

(14)

2.3 学习停止条件根据式(11)~式(12), 的取值直接影响失效概率的大小.分两种情况进行讨论.
1) 当时(即样本失效), 如果在区间[t₀, t_k], 存在某一点(x_i, t_j)满足≤0, 且对应的U(x_i, t_j)>2可以认为=1是可信的.
2) 如果, 需要证明在区间[t₀, t_k]内所有的样本点(x_i, t_j), j=1, …, k，对应的且U(x_i, t_j)>2才能认为=0是可信的.
定义某区间[t₀, t_k]中不准确的数量与(失效)的数量之比为误差系数, 用δ(t₀, t_k)表示.为使失效概率误差在0.05以内, 则δ(t₀, t_k) < 0.05, ?k=0, …, K.因此, AK-MCS-T方法的停止条件为

(15)

2.4 AK-MCS-T方法计算流程步骤1 ?将区间[t₀, t_e]离散化为K+1个时刻, 使用Monte Carlo方法建立公式(2)所示的样本集合S;
步骤2 ?应用拉丁超立方抽样选取N₀个初始样本点x_i, v_i(t_i), t_i, i=1, …, N₀;
步骤3 ?计算真实响应y=[y₁, y₂, …, y_N0], 并建立Kriging模型;
步骤4 ?使用当前Kriging模型计算集合S中样本对应的和U(x, v(t), t);
步骤5 ?根据式(11)计算各时间区间[t₀, t_k]的失效概率估计值和δ(t₀, t_k), k=0, …, K;
步骤6 ?根据学习停止条件判断Kriging模型精度是否满足要求, 如果满足要求执行步骤9, 不满足则执行步骤7;
步骤7 ?根据公式(13)选取最佳样本点;
步骤8? 将最佳样本点加入到DoE中, 计算其对应的真实响应, 并更新Kriging模型, 返回步骤3;
步骤9 ?输出结果.
3 动态可靠性算例3.1 算例1:腐蚀梁结构动态可靠性分析图 1为腐蚀梁结构^[10], 该梁长度为L=5 m, 受到作用点为中点处大小随时间发生变化的随机过程载荷F(t)作用, 及分布均匀的重力作用(ρ_st=78.5 kN/m³), 根据文献[11, 13], F(t)可表示为

(16)

图 1(Fig. 1)

图 1 腐蚀梁结构Fig.1 The corroded beam

式中：ξ_i (i=1, 2, …, 7)服从正态分布, 分布参数如表 1所示；a_ij, b_ij, c_ij数值见矩阵a, b, c.
表 1(Table 1)

表 1 算例1变量参数Table 1 Variable parameters of example 1 变量 分布类型 均值 标准差 σu/Pa 正态 2.4×108 2.4×107 a0/m 正态 0.2 0.01 b0/m 正态 0.04 4×10-3 ξ1 正态 0 100 ξ2 正态 0 50 ξ3 正态 0 98 ξ4 正态 0 121 ξ5 正态 0 227 ξ6 正态 0 98 ξ7 正态 0 121

表 1 算例1变量参数 Table 1 Variable parameters of example 1

梁受到腐蚀作用的影响, 其横截面积随时间变化满足:

(17)

式中: a₀, b₀是梁的初始截面尺寸; k是腐蚀速度, k=5×10^-5m/年.
t时刻, 梁中点处对应的弯矩以及发生屈服时对应的弯矩分别表示为

(18)

(19)

式中, σ_u是材料的屈服极限.
腐蚀梁的功能函数为

(20)

式中，X=[σ_u, a₀, b₀]^T，F(t), t是输入变量.
使用AK-MCS-T方法对腐蚀梁结构进行动态可靠性分析, 将时间区间[0, 35]等分成100个时间段, 即K=100.设定拉丁超立方抽取样本数为10.按照2.4节所示流程进行分析.
图 2中显示了各子区间对应的失效概率估计值和误差系数δ(0, t_k), k=0, …, K, 随样本点数变化情况, 其中样本点数用N_call表示, 失效概率标准值通过MCS得到, 其中PHI2、独立EGO、混合EGO方法的结果数据来自文献[10].通过图 2和图 3可以发现, 使用AK-MCS-T方法可以使失效概率估值很快地收敛到标准值.
图 2(Fig. 2)

图 2 算例1的和δ(0, t)曲线Fig.2 The curves of and δ(0, t) for example 1 (a)—Ncall=15；(b)—Ncall=30；(c)—Ncall=45.

图 3(Fig. 3)

图 3 算例1的δmax曲线Fig.3 The curves of δmax for example 1

表 2基于区间[0, 35]对应的失效概率对不同方法进行比较.通过样本点数和失效概率的对比可以发现, AK-MCS-T方法在保证精度的前提下, 调用样本点数远小于其他算法, 具备更高的计算效率.
表 2(Table 2)

表 2 算例1不同可靠性方法比较Table 2 Comparison of different methods for example 1 方法 样本点数 相对误差/% MCS 108 3.02×10-2 — PHI2 6501 2.85×10-2 5.63 独立EGO 496 3.27×10-2 8.28 混合EGO 283 2.99×10-2 0.99 AK-MCS-T 41 3.01×10-2 0.33

表 2 算例1不同可靠性方法比较 Table 2 Comparison of different methods for example 1

3.2 算例2:悬臂管结构动态可靠性分析悬臂管结构^[14]如图 4所示, 悬臂管结构承受2个定载荷(F₂和P)以及2个时变载荷(F₁(t)和T(t)), 材料的屈服强度随时间减小.

(21)

图 4(Fig. 4)

图 4 悬臂管结构Fig.4 Structure of cantilever tube

式中, σ₀是初始屈服强度.
悬臂管功能函数表示为

(22)

(23)

(24)

(25)

(26)

(27)

其中: L₁=60 mm; L₂=120 mm; θ₁=10°; θ₂=5°; 输入变量X=[d, h, σ0, F₂, P]^T; Z(t)=[F₁(t), T(t)]^T.参数分布信息如表 3所示, 时间区间为[0, 5]，算例2的变化曲线见图 5.
表 3(Table 3)

表 3 算例2变量参数Table 3 Variable parameters of example 2 变量 分布 均值 标准差 自相关系数 d/mm 正态 42 0.499 8 — h/mm 正态 5 0.1 — σ0/MPa 正态 560 56 — F1(t)/N 高斯过程 1 800 180 exp(-|Δt|/4) F2/N 正态 1 800 180 — P/N Gumbel 1 000 100 — T(t)/(N·mm) 高斯 1.9×106 1.9×105 exp(-4|Δt|2)

表 3 算例2变量参数 Table 3 Variable parameters of example 2

图 5(Fig. 5)

图 5 算例2的和δ(0, t)变化曲线Fig.5 The curves of and δ(0, t) for example 2 (a)—Ncall=15；(b)—Ncall=30；(c)—Ncall=45.

使用AK-MCS-T方法进行可靠性分析, 将区间[0, 5]分成80个时间段, K=80.已知功能函数的输入变量为7, 设定拉丁超立方抽取样本数为15.悬臂管结构失效概率估计值和误差系数随样本点数变化趋势如图 6所示.不同方法的结果对比如表 4所示.其中PHI2和t-PCE方法的数据结果来自文献[14].
图 6(Fig. 6)

图 6 算例2的δmax曲线Fig.6 The curve of δmax for example 2

表 4(Table 4)

表 4 算例2不同可靠性方法比较Table 4 Comparison of different methods for example 2 方法 样本点数 相对误差/% MCS 108 1.01×10-2 — PHI2 22 468 1.39×10-2 37.62 t-PCE 11 070 9.39×10-3 7.03 AK-MCS-T 37 0.99×10-2 1.98

表 4 算例2不同可靠性方法比较 Table 4 Comparison of different methods for example 2

通过图 5，图 6可以发现, 通过AK-MCS-T方法得到的各区间失效概率可以很快地收敛到标准值.在表 4中, 将多种动态可靠性方法进行比较, 可以发现AK-MCS-T方法明显优于其他几种方法, 且具有更高的计算效率.
4 结论1) 为提高计算效率, AK-MCS-T方法提出一种新的选点准则, 同时考虑到样本点对应的U函数和概率密度, 选取对失效概率精度影响最大的点更新Kriging模型.
2) AK-MCS-T方法可以计算各时间区间内的失效概率估计值, 得到失效概率随时间的变化趋势.
3) 为保证计算结果的精度, 提出了新的学习停止条件.
4) 通过2个算例, 与其他方法对比, 可以发现AK-MCS-T方法具有更高的计算效率, 且适用于高维非线性问题, 具有很高的工程价值.
参考文献

[1]	Au S K, Beck J L. A new adaptive importance sampling scheme for reliability calculations[J]. Structural Safety, 1999, 21(2): 135-158. DOI:10.1016/S0167-4730(99)00014-4
[2]	Zhang D, Han X, Jiang C, et al. Time-dependent reliability analysis through response surface method[J]. Journal of Mechanical Design, 2017, 139(4): 1-12.
[3]	Wei P, Wang Y, Tang C. Time-variant global reliability sensitivity analysis of structures with both input random variables and stochastic processes[J]. Structural and Multidisciplinary Optimization, 2017, 55(5): 1883-98. DOI:10.1007/s00158-016-1598-8
[4]	曹汝男, 孙志礼, 张毅博, 等. 基于主动学习的复杂机械结构的可靠性分析[J]. 东北大学学报(自然科学版), 2020, 41(2): 223-228. (Cao Ru-nan, Sun Zhi-li, Zhang Yi-bo, et al. Reliability analysis of complex mechanical structures based on active learning[J]. Journal of Northeastern University (Natural Science), 2020, 41(2): 223-228.)
[5]	Rice S O. Mathematical analysis of random noise[J]. Bell Labs Technical Journal, 1944, 23(3): 282-332. DOI:10.1002/j.1538-7305.1944.tb00874.x
[6]	Andrieu-Renaud C, Sudret B, Lemaire M. The PHI2 method: a way to compute time-variant reliability[J]. Reliability Engineering & System Safety, 2004, 84(1): 75-86.
[7]	Moarefzadeh M R, Sudret B. Implementation of directional simulation to estimate outcrossing rates in time-variant reliability analysis of structures[J]. Quality and Reliability Engineering, 2018, 34(8): 1818-1827. DOI:10.1002/qre.2374
[8]	王正, 谢里阳, 李兵. 随机载荷作用下的零件动态可靠性模型[J]. 机械工程学报, 2007, 43(12): 20-25. (Wang Zheng, Xie Li-yang, Li Bing. Part dynamic reliability model under random load[J]. Journal of Mechanical Engineering, 2007, 43(12): 20-25. DOI:10.3321/j.issn:0577-6686.2007.12.004)
[9]	Wang Z, Wang P. A nested extreme response surface approach for time-dependent reliability-based design optimization[J]. Journal of Mechanical Design, 2012, 134(12): 121007. DOI:10.1115/1.4007931
[10]	Hu Z, Du X P. Mixed efficient global optimization for time-dependent reliability analysis[J]. Journal of Mechanical Design, 2015, 137(5): 051401. DOI:10.1115/1.4029520
[11]	Hu Z, Mahadevan S. A single-loop Kriging surrogate modeling for time-dependent reliability analysis[J]. Journal of Mechanical Design, 2016, 138(6): 061406. DOI:10.1115/1.4033428
[12]	Echard B, Gayton N, Lemaire M. AK-MCS: an active learning reliability method combining Kriging and Monte Carlo simulation[J]. Structural Safety, 2011, 33(2): 145-154. DOI:10.1016/j.strusafe.2011.01.002
[13]	Hu Z, Du X. A sampling approach to extreme value distribution for time-dependent reliability analysis[J]. Journal of Mechanical Design, 2013, 135(7): 071003. DOI:10.1115/1.4023925
[14]	Hawchar L, El Soueidy C P, Schoefs F. Principal component analysis and polynomial chaos expansion for time-variant reliability problems[J]. Reliability Engineering & System Safety, 2017, 167: 406-416.