在系统随机不确定性分析与设计研究领域, 为全面衡量输入变量对输出参数性能的影响, 众多学者进行了大量的学术与工程应用研究, 建立了以重要性测度(重要度)来衡量这种影响的理论体系^[1]。而重要度(灵敏度)分析又可以分为局部重要度分析和全局重要度分析。局部重要度定义为输入变量分布参数的变化引起输出性能统计特征变化的比率^[2]。这一比率可以用输出性能统计特征对分布参数的偏导数来描述, 对此, 已有很多国内外学者进行了细致的研究。但是局部重要度存在明显的局限性, 即只能反映输入变量在名义值处对输出性能统计特征的影响, 而不能反映输入变量的完整不确定性对输出性能统计特征的影响。而全局重要度分析可以从输入变量的整个分布范围来衡量其不确定性对我们所感兴趣的输出统计特征的影响^[1]。

在全局重要度研究领域, Iman和Hora^[3]于1990年提出了理想的重要性测度指标的计算方法, 但是该方法需要同时求解三种不同的重要度指标, 计算量大。Sobol^[4]于1993年建立了一系列基于方差的基本变量重要性测度指标。Saltelli等^[5]在Sobol重要性测度指标的基础上总结出基本变量重要性测度指标的要求, 即“全局性, 可量化性, 通用性”, 并提出了相应的指标。此外, Homma和Saltelli^[6]还引入了总效应指标表征变量自身及其与其他变量相互作用对输出性能指标方差的影响, 并给出有效的计算方法^[7]。

失效概率是不确定性系统的一种重要性能指标, 而基于失效概率的全局重要度研究还未全面展开, 研究成果较少。Li和Lu等^[8]提出了基于失效概率的全局重要度, 并推导了其与基于方差的重要度指标之间的关系, 为可靠性评估与全局重要性测度分析建立了新的桥梁。基本变量对失效概率的重要性测度的求解涉及到无条件失效概率和有条件失效概率的计算对于求解失效概率, 通常有近似解析法和数字模拟法^[9]。近似模拟法包括一次二阶矩法和改进的一次二阶矩法等, 而数字模拟法的经典代表就是直接蒙特卡罗法(Monte Carlo, MC)。由于数字模拟法在失效概率评估和重要度分析上具有较好的相容性, 因此在基于失效概率的全局重要度研究中, 数字模拟法更为通用。李贵杰等^[10]建立了鞍点线抽样方法。Wei和Lu等^[11]给出了3种计算基于失效概率的全局重要度的数字模拟法, 包括直接蒙特卡罗法、重要抽样法(importance sampling, IS)和截断重要抽样法(truncated importance sampling, TIS)。Cheng等人研究了拒绝抽样技术在基于失效概率的全局重要度分析中应用^[12]。对比而言, 重要抽样法和截断重要抽样法的计算效率较高, 但是存在设计点位置及个数难以确定的困难, 进一步引起重要抽样函数选择困难, 降低了其应用范围。

交叉熵方法是一种重要抽样法, 其本质是对重要抽样函数选择过程的改进, 可以部分避免确定设计点位及个数。因此本文将建立求解基本变量对失效概率影响的重要性测度评估的交叉熵方法, 并对所建立方法进行计算精度和效率方面的验证, 为失效概率重要性测度评估提供一种新的可选方法。

1 基于失效概率的全局重要性测度

在包含n维不确定性变量的可靠性模型Y=g(x)=g(x₁, x₂, …, x_n)中, x=[x₁, x₂, …, x_n]^T为基本输入随机变量, g(x)为功能函数。记无条件失效概率值为P_fY, 即P_fY=P{g(x)≤0}, 当基本输入变量x_i取其实现值x_i^*时, 记条件失效概率值为P_fY|xi。基本输入变量x_i不确定性的消失会对失效概率造成一定的影响, 这种影响可以用P_fY和P_fY|xi的差异来度量^[10]。当考虑x_i在其整个分布域内变化时对失效概率的影响, 可以建立一个反映输入变量x_i对失效概率影响程度的失效概率重要度指标。定义基本输入随机变量x_i基于失效概率的矩独立重要性测度为^[13]

(1)

式中，F表示由功能函数g(x)定义的失效域, 即F={X:g(x)≤0}, f_Xi(x_i)为输入变量x_i的概率密度函数。

(1) 式中的绝对值符号不利于运算, 所以将其等价转化为平方运算, 等价转换后的输入变量x_i对失效概率的重要性测度记为δ_i^P, 其表达式如下^[8]

(2)

而上式中的2个失效概率, 其精确表达式为输入变量的联合概率密度函数在失效域中的积分。所以无条件失效概率P_fY可以改写为下式所示的失效域指示函数I_F的数学期望形式

(3)

式中，Rⁿ为n维变量空间, I_F为指示函数, 其定义为

(4)

类似的, 条件失效概率P_fY|xi可以改写为

(5)

将(3) 式和(5) 式带入(2) 式中, 就可以得到失效概率矩独立重要性测度与相应的失效域指示函数基于方差的重要性测度之间的关系

(6)

将基于失效概率的重要性测度转化为指示函数的方差形式将会给数值模拟带来极大便利。在(6) 式的基础上, 定义一组响应量为I_F的基于方差的重要度指标, 即基于失效概率的矩独立重要度指标的方差形式^[11]

(7)

(8)

式中，x_~i为输入变量x中除去x_i的其他分量的集合。定义S_i为主效应, S_i^T为总效应。

由全方差公式和(3) 式得

(10)

2 基于交叉熵的全局重要性测度

在上一节中, 将(2) 式转化为方差形式, 进而可以通过一些高效的抽样方法来求解原本比较难以求解的积分问题。

传统的重要抽样法是在蒙特卡罗法的基础上引入重要抽样密度函数, 从而使得有较多的样本点落入失效域^[14]。

对于(3) 式的求解, 传统的重要抽样法的公式为

(11)

式中，h_X(x)为重要抽样密度函数。通常将抽样中心选择在设计点处, 以有效地增加样本点落入失效域的概率, 从而减少抽样次数, 提高计算效率。

重要抽样法虽然拥有计算效率高的优点, 但是在实际使用过程中也存在2个困难:① 随着维数的增加, 失效概率的估计值通常会与真实值相差较远, 通常是低估失效概率; ② 当设计点难以求得或者设计点数目较多时, 构造重要抽样密度函数将变得十分困难。而事实上重要抽样密度函数的好坏将直接影响计算结果精度和效率。

相比较传统的重要抽样法, 交叉熵方法(cross-entropy method)可以一定程度上避免设计点求解的困难。交叉熵的概念由Rubinstein^[15]于1997年最先提出用于稀疏事件的模拟, 随后Rubinstein又将其拓展为随机优化方法^[16]。

交叉熵方法实际上是一种自适应重要抽样法。与传统的重要抽样法相同, 交叉熵方法的思想也是将抽样中心移至对失效概率贡献较大的区域, 不同点在于:使用了序贯方法确定重要抽样函数, 避免了难以确定重要抽样函数的困难。

交叉熵定义为2个概率密度函数比值对数的期望

(12)

从(12) 式可以看出, 当f_X(x)=h_X(x)时, D(f, h)=0。由此可以给求解最优重要抽样密度函数提供一个思路:设重要抽样密度为h_X(x; v), 其中v为未知的分布参数, 通过选择恰当的参数v的值, 使得重要抽样密度h_X(x; v)与最优重要抽样密度h_X^*(x)之间的差异最小, 即转化为D(h^*, h)的最小化问题。

最优重要抽样密度函数的理论值为

(13)

式中，θ为分布参数。将(13) 式和h_X(x; v)带入(11) 式, 即可将最小化问题转化为最大化问题

(14)

对于(14) 式, 可以首先采用Monte Carlo法, 从原分布f_X(x; θ)中抽取N个随机样本点{x₁, x₂, …, x_N}则(14) 式可以转化为

(15)

对于稀疏事件, 仍然存在极少数点落入失效域的问题, 因此可以再一次进行重要抽样, 其重要抽样密度为h_X(x; w), 则(14) 式可以等价为

(16)

定义一个权函数

(17)

从重要密度函数中抽取N个随机样本点{x₁, x₂, …, x_N}, 可以得到

(18)

可以证明(18) 式的优化问题关于v是凸的且可微, 并且对于常用的正态分布、指数分布和威布尔分布等自然指数分布族, (18) 式存在解析解, 其最优问题与下式的最小化问题等价

(19)

可以得出最小化问题的最优解

(20)

(21)

当失效概率非常小时, 直接应用此结论依然十分困难, 因为大部分的随机样本还是落在失效域之外。因此Rubinstein等人提出了多层交叉熵方法(multi-level cross entropy method)。这种方法是通过增加一个控制参数ρ——中间失效概率, 产生中间失效事件, 并使得中间失效事件逐渐向失效域靠近。由此构建交叉熵算法:

1) 初始化待定分布参数v₀=θ, 迭代计数k=0。

2) 从h_X(x; v_k)产生N个随机样本{x₁, x₂, …, x_N}, 计算功能函数值{g(x_i)}, 并且对{g(x_i)}从小到大排序, 得到{g_i}, 求得功能函数值次序统计量的ρ分位数b_k=g_[ρN]。

3) 使用相同的N个随机样本{x₁, x₂, …, x_N}, 解(19) 式，得到新的参数v_k+1的估计值。

4) 若b_k>0, 设k=k+1, 转步骤2)。

5) 基于最终的重要抽样密度, 生成所需的随机样本。

通过引入交叉熵的概念, 改进传统的重要抽样法, 可以得到求解重要性测度的高效抽样求解法。交叉熵方法求解主效应S_i和总效应S_i^T的步骤为^[11]:

(1) 由交叉熵法求得重要抽样密度函数h_X(x), 由h_X(x)生成2N个样本{x₁, x₂, …, x_2N}, 将这2N个样本分配给2组N×n的样本矩阵A和B(n为输入变量的维度), 由此得到的矩阵A和B是独立同分布的。这组样本可以通过简单随机抽样、拉丁方抽样、Sobol序列等抽样技术得到, 一般来说, Sobol序列产生的样本收敛速度优于其他方法。得到的A和B矩阵如下所示

(22)

(23)

(2) 求得另一组N×n的样本矩阵Cⁱ(i=1, 2, …, n), 其中Cⁱ的第i列来自于矩阵A, 其余的n-1列来自于矩阵B。

(24)

(3) 对样本矩阵A、B和Cⁱ, 计算(25) 式中的4个N维向量。

(25)

式中，ω(A)和ω(B)中的元素ω_A(x_i)(i=1, 2, …, N)和ω_B(x_i)(i=N+1, …, 2N)为指示函数与权函数的乘积。

(26)

函数F_A是一个2n-1维的函数

(27)

函数F_B是一个k+1维的函数

(28)

由此可得4个N维向量, 即公式(25)。

(4) 主效应S_i和总效应S_i^T的值由(29) 式和(30) 式进行计算。

(29)

(30)

式中

(31)

(32)

(33)

(34)

(35)

式中上标j表示向量中第j个元素。

(5) Homma和Saltelli^[6]还给出了计算全局重要度的误差估计方法, 结合Wei等^[11]的研究, 由(36) 式和(37) 式计算估计值主效应S_i和总效应S_i^T的概率误差。

(36)

(37)

式中

(38)

(39)

(40)

3 算例

为了验证本文所提方法的性能, 本节对2个数值算例和1个工程算例进行分析。分析时将蒙特卡罗法的结果作为精确解。通过对比失效概率和重要性测度的计算结果, 以及所使用的样本量来衡量方法的计算精度和效率。

算例1  考虑1个包含2个设计点的问题, 即

(41)

式中，x₁和x₂服从标准正态分布。该问题存在2个设计点x₁^*={2.55, 0.50}^T和x₂^*={-2.55, 0.50}^T。

首先采用直接蒙特卡罗法在标准正态空间内求解失效概率和重要度。取样本量N=10⁶, 可得失效概率估计值为。

在使用重要抽样法时, 重要抽样密度需要很好地覆盖所有对失效概率贡献大的区域, 为了模拟工程中不知道设计点的情况, 假设不知道设计点2的存在。在仅知道设计点1的情况下, 构建如下的重要抽样密度

(42)

样本量为N=3 000, 失效概率估计值为=4.660×10^-3, 相比较真实值大约低估了一半。由此可见设计点的信息对重要抽样法至关重要, 因此重新构造一个混合高斯分布作为重要抽样密度

(43)

在重新构建重要抽要密度之后, 失效概率估计值为=9.077×10^-3。

使用多层交叉熵方法求解此问题, 设定初始参数ρ=0.1, 单层样本量N₀=1 000, 重要抽样样本量N=3 000。在1 000个样本中, 有498个样本点落入失效域, 迭代次数2次。这说明该方法的抽样效率非常高。得到重要抽样密度的均值和方差分别为={0.199, 0.499}^T和={3.799, 0.922}^T。最终失效概率估计值为=9.129×10^-3。

在验证了传统的重要抽样法和多层交叉熵方法计算失效概率的精度后, 采用第二节提到的方法, 使用交叉熵方法计算输入参数的重要性测度, 并与直接蒙特卡罗法和重要抽样法的结果作比较, 其结果如表 1所示。表中括号内数值表示概率误差。

从表中数据可以看到一个有意思的现象, 就本算例而言, 缺失一个设计点信息导致传统重要抽样法错误地估计了失效概率, 而重要度的计算结果却是正确的。这是由于重要性测度指标中的归一化处理和极限状态函数轴对称形态引起的, 所以分布在其中一个设计点附近的样本点对其重要度的贡献与所有样本点对总体的贡献是一样的。但是这仅仅是针对这个算例的特殊情况, 并不能说明传统重要抽样法在错误计算失效概率的情况下能够得到正确的重要度。相反, 说明在全局重要性测度评估中存在陷阱。

如果是设计点已知的情况, 采用传统重要抽样法和交叉熵方法都可以得出较为精确的结果, 既有输入变量x₁的重要度要大于x₂。而且这2种方法相比较直接蒙特卡罗法, 所需样本量明显少, 有效地提高了计算效率。对比后2组数据, 发现交叉熵方法在确定重要抽样密度的参数时, 迭代次数为2次, 单层样本量1 000, 加上最后使用重要抽样所用的3 000个样本, 一共使用了5 000个样本。因此可以看出, 相比于已知设计点的传统重要抽样法, 交叉熵方法的计算效率略有降低。这是由于在设计点已知时, 重要抽样法的样本点可以很好地覆盖设计点, 而交叉熵方法则是通过扩大方差的方法增大样本点的覆盖范围, 以较好地覆盖所有设计点。图 1分别是设计点未知的重要抽样法、设计点已知的重要抽样法和交叉熵方法所使用的随机样本在平面上分布情况。

从图中可以看出, 在已知设计点的前提下, 所构造的重要抽样密度的样本很集中地分布在2个设计点的附近, 因此抽样效率极高; 而多层交叉熵方法的抽样中心位于2个设计点的中间, 由于x₁的方差比较大, 所以样本点也可以比较好地覆盖2个设计点, 而相应的缺点就是计算效率和精度有所下降。不过相比较直接蒙特卡罗法, 交叉熵方法依然是大大提高了计算效率。

算例2  考虑如下的数值算例

(44)

其中3个输入变量x_i(i=1, 2, 3) 均服从标准正态分布。使用直接蒙特卡罗法, 重要度排序为:x₃>x₂>x₁。将此结果作为标准值。

使用传统重要抽样法, 通过Matlab中的fminsearch函数求设计点。设计点为极限状态方程上到原点距离最小的点, 即转化为(45) 式的最小化问题

(45)

式中

(46)

即可转化为关于x₂和x₃的最优化问题。当x₂和x₃取不同的初值时, 优化结果也不尽相同。表 2给出了若干不同初值下求得的设计点以及相应的失效概率和重要度计算结果。

失效概率和重要性测度的计算结果和直接蒙特卡罗法差异较大。实际上, 在增加样本量降低估计值变异系数后, 传统重要抽样法的计算结果并没有显著变好。这说明计算结果的不准确并不是样本量不足造成的, 而是因为所求得的设计点不甚理想。事实上像这种带有三角函数的极限状态函数, 由于三角函数的周期性, 设计点的计算会比较困难, 而设计点的计算误差将很大程度上影响失效概率和重要性的评估。

采用交叉熵方法, 设定初始参数ρ=0.07, 单层样本量N₀=1 000, 重要抽样样本量N=5 000。其中控制参数是通过参数寻优求得。计算过程中, 多层交叉熵的迭代次数仅为1次, 总计使用6 000个样本点。抽样中心的计算结果为:均值={0.447, 0.309, 1.486}^T, 方差={1.040, 1.749, 1.066}^T。使用交叉熵方法之后, 不论是失效概率还是重要度, 结果都与标准值十分吻合, 并且相比蒙特卡罗法效率大大提高。

由MC法、IS法和CE法求得的3个输入变量的主效应和总效应以及相应的失效概率误差见表 3。

算例3  本算例来自文献[17], 考虑如图 2所示的工字梁问题。该问题包含8个输入变量, x=[P, L, a, S, d, b_f, t_w, t_f]^T, 各变量的分布信息如表 4所示。

该问题的极限状态函数为

(47)

式中，S为强度, σ_max为最大应力, 最大应力的计算表达式为

(48)

式中

(49)

使用直接蒙特卡罗法和交叉熵方法的计算结果如表 5所示。

根据交叉熵方法的计算结果, 输入变量重要性排序为: t_f, L, a, P, S, d, t_w, b_f (由大到小)。这与使用蒙特卡罗法得到的结果相吻合, 也与文献中使用鞍点线抽样方法得到的结果完全一致, 证明了引入交叉熵方法计算重要性测度在工程中的可行性。值得注意的是:本算例中, 交叉熵方法相比蒙特卡罗法, 计算效率的提高并没有算例1和算例2中明显。其原因在于本算例中失效概率较大, 即使是使用直接蒙特卡罗法, 也有大约20%的样本点落入失效域。也就是说, 通过转移抽样中心来提高抽样效率的方法, 在小失效概率事件中会有更为显著的优势。

4 结论

本文针对失效概率的重要性测度分析, 提出了一种基于交叉熵方法的高效评估方法。与传统重要抽样方法相比, 避免了设计点位置和个数的求解。所提方法通过自适应移动抽样中心的方法, 逐步确定一个合理的重要抽样函数, 回避多设计点问题, 使得更多的样本点落入失效域, 以提高在失效概率和重要性测度的计算效率。相比于传统的重要抽样法, 交叉熵方法可能略微降低了计算效率, 但是免去了设计点求解的困难, 适用于非线性极限状态函数和多设计点的复杂情况, 为重要性测度研究提供了一个新的思路和方向。文中算例不但验证了该方法的优越性, 同时也提供了传统重要抽样法明显不足的典型例证。