随着现代机械技术的发展, 其结构复杂度也随之增加, 导致产品极限状态函数(limited state function, LSF)呈现高度非线性且为隐函数。传统的可靠性解析手段如一阶可靠性方法(first order reliability method, FORM)、二阶可靠性方法(second order reliability method, SORM)等, 是在最可能失效点处对功能函数进行一次或二次函数展开。当极限状态函数为隐函数时很难得到每一次迭代计算过程中的梯度向量。蒙特卡罗抽样(monte carlo sampling, MCS)有效解决极限状态函数为隐函的问题, 且计算精度高, 被广泛应用于结构可靠性计算。然而, MCS方法会造成巨量计算成本, 所以极大限制了其在实际工程应用中的推广^[1-2]。

代理模型可靠性分析手段从20世纪开始发展到现在, 出现了各种各样的模型结构, 主流的代理模型方法有多项式响应面、支持向量机、Kriging模型和神经网络，其中Kriging模型是一种包含多项式和随机过程函数的模型结构。国内外为了提高Kriging在可靠性指标求解过程中的计算精度, 提出了各种不同的Kriging模型构建方法, 总体可分为3类:①通过优化Kriging模型中相关函数的参数θ来提升模型精度。魏娟等利用混合粒子群-模拟退火算法进行参数求解, 提高模型预测精度; 刘瞻等利用人工蜂群算法对Kriging模型参数进行优化, 再用优化后的模型模拟隐式极限状态函数^[3-4]。②通过优化求解最佳样本点, 即通过主动学习方法序列选择最佳样本点, 然后重新拟合Kriging模型, 逐渐提高Kriging模型精度。从该方法发展到现在出现了各式各样的主动学习选点策略。Jones等提出一种EGO(efficient global optimization)方法, 其中采用EI(expected improvement)函数来选择优化过程中的最佳样本点; Bichon在EGO方法的基础上, 提出了一种EFF(expected feasibility function)函数选择最佳样本点, 着重考虑在极限状态函数附近且方差较大的样本点; Echard提出一种U主动学习方法, 考虑样本点符号被错误估计概率最大的点; LYU等提出H主动学习函数, 结合线抽样方法进行可靠性求解; Sun等提出一种最小改进函数(LIF), 同时考虑Kriging方差和联合概率密度函数值对失效概率估计精度的影响, 避免对概率密度函数值过小的区域抽样造成样本浪费^[5-9]。③选择不同的数字仿真方法求解结构失效概率, 在候选点集的构造和失效概率计算的过程中都采用了数字仿真方法产生样本点。在一般情况下, 使用拉丁超立方采样(latin hypercube sampling, LHS)和MCS方法来构造初始点集, 成为初始设计点, 并结合MCS、IS、子集模拟法(subset simulation, SS)、线抽样(line sampling, LS)等不同抽样方法和Kriging模型结合计算结构可靠性指标^{[8, 10-13]}。

本文基于Kriging模型可靠性计算框架, 考虑Kriging模型主动学习过程中对全局和局部搜索能力的综合要求, 针对传统EGO主动学习方法选点位置不对等问题, 提出一种IEGO(improved efficient global optimization)主动学习方法, 该方法既能在极限状态面附近局部搜索最佳样本点, 又能根据模型预测方差的大小全局寻找最佳样本点, 进而不断修正Kriging模型的精度, 最后结合蒙特卡罗方法完成结构可靠性指标计算。

1 Kriging模型

Kriging模型包含2个部分:多项式回归和随机部分^[14]。随机部分通常有高斯过程、指数过程等, 在本文里面统一使用高斯过程进行计算。Kriging模型需k个初始实验点X=[x₁, x₂, …, x_k]^T, 其中第i个实验点x_i∈Rⁿ(i=1, 2, …, k), 并根据相关分析手段得到这k个实验点的真实响应值G=[G(x₁), G(x₂), …, G(x_k)]^T, 其中第i个实验点对应的响应值为G(x_i)∈R(i=1, 2, …, k)。则通过Kriging代理模型的结构形式如下

(1)

式中, x为任意输入变量，为对应预测响应值，β^T=[β₁, β₂, …, β_p]是一个回归系数向量, f^T(x)=[f₁(x), f₂(x), …, f_p(x)]是多项式的基函数向量。研究表明, 第一部分F(β, x)的具体形式对整个Kriging模型的拟合精度影响不大, 因此将多项式f^T(x)取为固定常向量I, 此时F(β, x)=β。z(x)是服从N(0, σ²)的高斯随机过程, 样本点中任意2个样本点x_i和x_j的协方差方程为

(2)

式中：R(θ, x_i, x_j)是表征样本点x_i和x_j空间相关性的相关函数；θ=[θ₁, θ₂, …, θ_n]^T为相关性参数。本文相关性函数采用高斯相关函数, 其具体形式为

(3)

式中, x_i^e和x_j^e分别表示样本点x_i和x_j的第e个方向分量; θ_e是对应着第e方向分量的相关性参数, 表示设计样本点之间的距离和相关性程度的大小关系, 其大小直接控制着Kriging模型的预测结果。

基于给定的k个初始样本点X=[x₁, x₂, …, x_k]和与之相对应的k个真实极限状态函数响应值G=[G(x₁), G(x₂), …, G(x_k)], 对于给定的相关参数θ=[θ₁, θ₂, …, θ_n]^T, 回归系数β和高斯随机过程方差σ²的估计值分别为

(4)

(5)

式中: 为估计得到的多项式参数向量; 为方差σ²的估计值; R为相关系数矩阵, 且R=[R_ij]_k×k, R_ij=R(θ, x_i, x_j)。

为了保证在某一预测点x处的预测值无偏和预测均方误差最小, 在该点处的Kriging预测均值和方差分别为

(6)

(7)

式中：r(x)=[R(θ, x_i, x), R(θ, x_i, x), …, R(θ, x_i, x)]; u(x)=F^TR^-1r-f; Kriging方差表征着模型在预测点x处的预测误差, 是主动学习方法的重要依据。

2 基于主动学习的可靠度求解算法

主动学习方法的基本原理就是通过学习函数, 在大量Monte Carlo模拟样本点中选择最佳样本点, 并将该样本点加入初始样本实验点, 对Kriging模型进行更新, 如此循环迭代直至达到相关停止要求。通过主动学习方法, 仅需要少量的初始样本点和主动学习筛选策略就可以构建高精度的代理模型, 较传统方法较好提升了计算效率。

2.1 EGO主动学习方法

EGO主动学习方法是由Jones等初次提出的用于复杂函数的代理模型优化, 基于Kriging代理模型, 采用EI函数指标计算样本点处对模型的期望改善程度, 最大EI函数值对应的样本点用于模型的再次拟合进行模型本身的序列更新, 从而逐步提高模型精度, 得到函数极值的优化过程^[15]。

对于空间任何一点x, 其预测响应值服从高斯分布, 即

(8)

采用Monte Carlo抽样, 选择m个抽样样本点, 对应的预测响应值为G₁, G₂, …, G_m, 选择当前响应值中的最小值G_min=min(G₁, G₂, …, G_m), 定义在样本点x处对代理模型的改进程度为

(9)

式中, 是一个不确定量, 代表着模型在点x处响应值的不确定性, 所以I(x)也是一个随机变量。则在样本点x处对模型的改进程度函数的期望值为

(10)

式中:φ(·)和Φ(·)分别是标准正态分布的概率密度函数和累积分布函数。定义期望值E(I(x))最大所对应的样本点为当前Kriging模型和此次抽样样本点集下的最佳样本点, 即

(11)

计算最佳样本点x_best处的真实响应值G(x_best), 并加入到最初的k个样本集, 更新Kriging模型。然后, 再次在当前模型下进行大样本抽样, 计算样本点的E(I(x)), 直到满足迭代停止条件:max(E(I(x))) < 0.001。

2.2 可靠性计算中EGO方法的不适用性和IEGO主动学习方法

由2.1节可知EGO算法的选点策略是基于表征模型期望改进程度的函数, 它是从代理模型优化问题引申出来的一套主动选点策略。首先, 它选取当前Kriging模型的最小值作为最佳函数值G_min, 当抽样得到点x的时候, 并不知道其真实响应值G(x)大小, 因此我们假设G服从正态分布, 其分布均值和标准差由Kriging模型预测给出。未知点处真实响应值G < G_min的概率为

(12)

随机变量G处于直线G(x)=G_min下面概率刚好处于正态分布的尾端, 称直线G(x)=G_min下的线段长度为在点x处的改善程度。为了衡量该改善值的期望值, 根据在该点处真实响应值的概率密度函数求取改善值的期望函数值, 即

(13)

由(13)式可以看出, 在某一点x处影响期望改善值E(I(x))大小的主要因素有2个:①当前响应值大小, 即在点x处正态分布的均值大小。②当前响应值标准差大小, 即直接影响点x处离散程度的高低。

经典EGO方法的选点策略不符合可靠性计算求解理念, 在可靠性求解过程中, 主要着重于极限状态面处的拟合精度, 而不是在响应函数最小值处的拟合情况。因此, 基于经典EGO方法提出一种改进的EGO方法(IEGO)。通过改变期望改善值E(I(x)的定义方式, 使其最优点选取策略能优先考虑极限状态面附近的点和方差较大的点。

首先, 在原有Kriging模型预测值的基础上, 将部分按照直线进行翻转, 即在期望改善值的计算过程中, 点x处随机变量的分布均值变为

(14)

此时, 假设点x处的分布函数标准差保持不变, 仍为。可以将修正后的Kriging模型在点x处的预测值当做是服从均值为、标准差为的正态随机变量。而且, 当前响应值的最小值重新定义为

(15)

此时对函数的改善程度定义为

(16)

则期望改善值变为

(17)

式中：等式右边第一项是点x处改善程度和小于当前最小目标函数值概率的乘积, 搜索的重点放在 (即极限状态面)附近, 强化了极限状态面周围的局部搜索能力; 等式右边第二项是点x处预测标准差和正态密度函数的乘积, 强化了IEGO方法的全局搜索能力。所以, 采用改进的IEGO算法, 具备良好的平衡极限状态面区域局部和全局的搜索能力。

2.3 基于IEGO算法的可靠性计算流程

基于Kriging代理模型, 结合2.2节所提出的IEGO主动学习方法, 给出结构可靠性算法计算步骤如下:

1) 生成候选样本点。基于Monte Carlo抽样原理进行仿真, 根据基本随机变量的联合概率密度函数生成N_MC个随机样本S_MC=[x₁, x₂, …, x_NMC]^T。这些样本点作为后续主动学习过程中的学习样本, 称为候选样本点。

2) 生成初始设计样本点。采用拉丁超立方采样, 在基本随机变量空间(μ_i-3σ_i, μ_i+3σ_i) i=1, 2, …, n中生成初步拟合Kriging模型所需的试验设计样本点S_initial=[x₁, x₂, …, x_N]^T, 初始样本点的数量选择要适中, 并在后续的模型学习过程中不断去增加初始实验设计样本S_initial。

3) 初步拟合Kriging模型。计算初始实验设计样本点S_initial的真实响应函数值Y_initial, 由S_initial和Y_initial建立Kriging预测模型, 并计算候选样本点集中所有点x_i∈S_MC的预测值和标准差。

4) 根据(14)和(15)式计算样本点集内所有点新的预测值和当前最小值G_min^*, 并基于(17)式得到改进的期望改善值E(I^*(x_i))。

5) 对所有候选样本点期望改善值E^*(I(x_i))进行排序, 寻找最大值max(E^*(I(x_i)))

6) 若max(E(I^*(x_i))) < 0.001, 学习结束, 到第7)步进行相关可靠性计算。否则, 选择最大值max(E(I^*(x_i)))函数所对应的样本点为最佳样本点x^*, 计算其极限状态函数响应值G(x^*), 并将新的数据(x^*, G(x^*))加入到初始设计样本集, 回到步骤3);

7) 根据Monte Carlo抽样模拟, 采用当前的Kriging模型来计算失效概率和变异系数

(18)

式中, N_MC为Monte Carlo抽样数量。当≤0.03时, 则可靠度计算收敛。

3 实际算例验证

为了验证本文算法的精确性和有效性, 提出3个实际算例验证, 其中包含可靠性方面的经典数值算例和一个工程算例。在每一个算例中, 基于Kriging代理模型, 将本文算法(IEGO)分别与EGO算法、EFF算法和U算法进行比较, 同时还和经典二次响应面作比较, 验证Kriging模型和IEGO算法组合在可靠性领域的优异性。以Monte Carlo模拟抽样计算结果P_f^MC为标准解, 分别从以下2个方面作比较:①调用极限状态函数次数N_call来评价方法的效率; ②计算各可靠性计算方法得到的失效率与P_f^MC的相对误差

(19)

来直观评价方法的计算精确性。

算例1   选择一个四分支串联系统, 其极限状态函数形式如(20)式所示^[16], 其中基本变量x₁和x₂为独立同分布标准正态变量。

(20)

首先在各基本变量均值附近±3σ范围采用拉丁超立方采样获得初始设计样本, 求取真实响应值用来构建初始Kriging模型, 并逐次加入最佳样本点进行模型更新。针对不同的主动学习函数, 最终得到其最佳Kriging模型和最佳样本点的分布如图 1所示。并添加MC计算方法和经典二项式响应面方法与Kriging代理模型进行对比, 最终的结果如表 1所示。

如图 1所示, 可知：①EGO方法本身就是适用于最小值优化过程中主动选点策略, 所以该主动学习过程, 在极限状态面附近的拟和精度较差; ②本文所提出的IEGO主动学习方法, 在极限状态面附近的拟合精度大大提升, 且每次模型更新的最佳样本点不过分集中, 基本均匀处在极限状态面上; ③U和EFF主动学习方法在极限状态面附近的拟合效果较好, 然而, U方法学习点相对比较靠拢, 可能造成区域的过度选点。相对于其他主动学习方法, EFF方法选点数量相对最少。

以MC方法的抽样模拟计算结果作为精确解, 由表 1可知, 针对这种串联型的可靠性计算, 模型拟合精度较差。其中EGO方法和EFF方法的调用极限状态函数次数最少, 然而由于EGO的选点策略不适用于可靠性计算等原因, 导致其失效概率计算结果误差很大; EFF方法由于选点过少, 模型在极限状态面附近拟合精度不高, 失效概率计算结果误差较大。IEGO和U方法, 在失效概率计算精度上, IEGO占有微弱优势, 但是由于U方法在极限状态面上选点过于集中, 导致调用极限状态函数次数较多, 大大增加计算成本。

算例2  选择Rastrigin非线性函数, 其函数具体形式如(21)式所示^[20]。

(21)

式中基本变量x₁和x₂为独立同分布的标准正态变量。针对不同学习方法的最终Kriging模型和更新过程中的最佳样本点分布如图 2所示。

由图 2可知:①EGO方法在极限状态面附近的拟合情况依然十分不理想, 选点主要在极限状态函数最小值附近, 导致最佳样本点集中在图 2a)的外围; ②本文所提出的IEGO主动学习方法, 对于这种复杂的极限状态面情况, 拟合效果较令人满意; ③U方法整体拟合效果不错, 然而可以看到在图 2c)外围的极限状态面拟合效果不佳, 甚至有些区域选点策略并没有考虑到; ④EFF方法在整个极限状态面的拟合效果也较好, 然而该方法选点相对较均匀, 比较图 2b)和图 2d)2个图中心部分, 可以看到EFF相对于本文提出的IEGO方法出现空白部分选点情况较多。

在算例2中, 对不同的模型可靠性计算结果进行对比如表 2所示, 可以看到经典二次响应面针对这种复杂极限状态函数的模型拟合精度不高, 导致计算结果很差。传统EGO方法, 虽然调用极限状态函数次数少, 然而由于选点位置问题, 导致相对误差较大; 本文所提出的IEGO主动学习方法, 和U方法及EFF方法相比, 调用极限状态函数次数差不多, 然而本文针对失效概率的计算误差相对较小。

算例3   图 3所示为十杆桁架结构^[17], 该结构在优化设计和可靠性计算领域得到广泛应用。由于桁架结构的水平杆、垂直杆以及倾斜杆采用3种不同的铝合金杆件, 分布将3种杆件结构的横截面积A₁, A₂和A₃考虑为基本随机变量, 且基本变量A₁~N(13, 1.3), A₂~N(2, 0.2), A₁~N(9, 0.9)(in²)。

假设结构其他参数均为确定性变量, 如弹性模量E=10⁷ psi(1 psi=6.895 kPa)、每根杆件长度L=360 in(1 in=2.54 cm), 顶点(2)处的真实载荷P=100 000 lb(1 lb=0.453 6 kg), 顶点(2)的许用位移值为d_allow=4 in。

所以结构的极限状态函数为

(22)

式中, A=(A₁, A₂, A₃)为不同杆的横截面积; d_real为顶点(2)的实际位移值。

针对算例3, 不同方法计算得到的失效概率和调用极限状态函数次数如表 3所示。其中二项式响应面方法调用极限状态函数次数较多, 计算失效概率相对误差也较大。IEGO方法相对于传统EGO方法, 不仅在计算效率(N_call少了5次)上有所提高, 同时计算精度上也提升了一个量级; U方法计算精度也较好, 然而其调用极限状态函数次数较多, 因此计算效率较差; EFF方法在几种方法中计算效率最高, 然而其计算精度相对较低。所以, 权衡计算精度(ε的大小)和计算效率(N_call)两方面特性, 本文所提出的IEGO方法具有较好的优势。

4 结论

本文基于模型结构的Kriging代理可靠度求解体系, 提出一种改进的主动学习函数IEGO, 使得Kriging预测模型能在有限的样本点信息下, 最大限度提高拟合精度, 保证计算得到的可靠性指标的正确性。

1) 所提出的IEGO方法能够使主动学习函数在极限状态面附近进行局部选点, 保证Kriging模型在极限状态面附近的拟合精度。同时, 还能兼顾全局方差过大地方的选点, 使代理模型不会有局部过拟合状况发生。

2) 2个经典数值算例和1个工程实例计算结果表明, 本文所提出的的IEGO主动学习方法, 在可靠性工程应用中, 能同时兼顾计算效率和计算精度2个方面问题, 具有较高的工程实用价值。