可靠性全局灵敏度求解的Meta-IS-AK算法

周苏婷， 吕震宙， 凌春燕，王燕萍

西北工业大学 航空学院，西安 710072

全局灵敏度分析可以从输入变量的整个分布范围来衡量输入变量的不确定性对工程设计中所感兴趣的输出性能统计特征的贡献程度[1-2]，通常在可靠性分析中，更关注结构是否失效以及输入变量对于结构失效概率的影响程度。为了衡量输入变量对于结构失效概率的影响并对输入变量的重要性进行排序，Cui等[3]提出了可靠性全局灵敏度指标。可靠性全局灵敏度指标定义为某一输入变量固定时，结构的无条件失效概率与条件失效概率的差异在输入变量的整个分布范围内的平均值，该指标反映了输入变量固定时对失效概率的平均影响。与可靠性局部灵敏度指标[4-5]不同，可靠性全局灵敏度指标可以从输入变量的整个不确定性范围来衡量其对于失效概率的平均影响。

本文提出了利用贝叶斯公式转换以及元重要抽样算法[6-7]结合自适应Kriging代理的可靠性全局灵敏度，并组织了求解可靠性全局灵敏度指标的高效算法。基于贝叶斯公式可以将可靠性全局灵敏度指标的原始定义等价表示为输入变量的无条件概率密度函数(Probability Density Fanction, PDF)与条件概率密度函数之间的差异与失效概率的乘积形式，从而避免条件失效概率的求解[8-9]。基于此等价形式，只需要一组样本，便可以得到各输入变量的可靠性全局灵敏度指标。利用贝叶斯算法进行可靠性全局灵敏度分析的研究已有很多，Wang等[8]利用贝叶斯公式，将可靠性全局灵敏度指标的定义式进行了等价转换，并采用蒙特卡罗法来计算可靠性全局灵敏度指标。该方法适用范围虽广，但在计算工程实际中的小失效概率问题时计算量过于庞大。Wang等[9]则结合重要抽样方法来计算可靠性全局灵敏度指标的贝叶斯等价形式。该方法只需一组样本便可以计算得到所有输入变量的可靠性全局灵敏度指标，且计算量独立于输入变量的维数，明显地提高了计算的效率。但是重要抽样方法[10-11]需要构造重要抽样密度函数，一般来说需要将重要抽样密度函数的中心放在设计点处，这就使得该方法依赖于其他方法来寻找设计点，并且重要抽样方法对于多设计点的情况和多失效域的情况并不适用。Yun等[12]基于可靠性全局灵敏度指标的贝叶斯公式转换后的等价形式，利用子集模拟方法[13-16]结合重要抽样的方法进行可靠性全局灵敏度的计算。该方法在利用子集模拟结合重要抽样的方法计算得到无条件失效概率的同时，通过重复利用子集模拟重要抽样方法抽取的样本，并利用Metropolis-Hastings(M-H)[17]准则来实现失效样本的转换，进而估计各输入变量的失效条件概率密度函数。该方法同样只需要一组样本便可计算得到可靠性全局灵敏度指标，且对于小失效概率问题是有效可行的。但该方法抽样过程中所结合的重要抽样方法仍然依赖于设计点的选取，因而该方法依旧很难适用于多设计点和多失效域的情况。

对于实际工程中普遍存在大量复杂的功能函数，自适应Kriging(Adaptive Kriging, AK)代理模型法是处理这类问题的高效算法[18]。可靠性分析中涉及到的自适应Kriging算法有自适应Kriging代理模型结合Monte Carlo模拟(Monte Carlo Simulation, MCS)法(AK-MCS)[19]、自适应Kriging代理模型结合重要抽样(Important Sampling, IS)法(AK-IS)[20]以及元重要抽样算法(Meta-IS)[6-7]等。AK-MCS算法将自适应Kriging代理模型与Monte Carlo模拟法相结合，其首先由输入变量的联合概率密度函数产生MCS备选样本池，然后利用学习函数在备选样本池中逐步挑选对失效面拟合贡献较大的点来更新Kriging模型，最终能够确保Kriging模型能在一定的置信水平下识别样本池内样本的功能函数值的正负号，从而较准确地求得失效概率。但对于工程上常见的小失效概率问题，AK-MCS算法的样本池容量庞大，进而使得代理过程十分耗时。将重要抽样法与自适应Kriging过程相结合形成的AK-IS算法可以大大降低备选样本池的规模，从而提高自适应学习的效率，减少自适应学习过程所消耗的时间。但是对于多设计点及多失效域问题，基于一次二阶矩法构造重要抽样密度函数的AK-IS算法将失效，因而对于多设计点及多失效域的问题需要另寻解决途径。元模型重要抽样算法(Meta-IS)通过元模型构造重要抽样密度函数，进而来抽取重要抽样样本点。该方法在降低AK-MCS样本池规模的同时避免了设计点的求解，适用于多设计点和多失效域问题的可靠性分析。然而，Meta-IS算法在抽取重要抽样样本后仍需要求解重要抽样样本点处真实的功能函数值来估计失效概率，所以该算法仍有较大的计算量。如果在元重要抽样的基础上嵌入自适应Kriging模型，形成Meta-IS-AK算法，则有可能极大程度地提高失效概率求解的效率，本文将采用这种思路求解失效概率，并在此基础上组织可靠性全局灵敏度的算法。

利用可靠性全局灵敏指标的贝叶斯转换形式，只要能求得无条件失效概率和失效域条件下输入变量的概率密度函数即可直接求得可靠性全局灵敏度，而求解无条件失效概率的数字模拟法可以同时完成失效概率的计算以及失效域条件下输入变量概率密度函数的求解，为此本文将Meta-IS与自适应Kriging模型结合起来，以便完成可靠性全局灵敏度的高效求解。本文所组织的Meta-IS-AK算法主要分3步来执行。第1步是由元重要抽样的迭代策略得到逐步逼近最优重要抽样函数的样本点，在迭代收敛后便可以得到第1步的Kriging模型以及相应的重要抽样样本点。第2步是基于重要抽样样本构建自适应Kriging模型，该步的Kriging模型是在第1步得到的Kriging模型的基础上逐步更新得到的，其目的是为了构建能够对第1步中重要抽样样本点失效与否做出准确预测的Kriging模型，从而高效地求解出无条件失效概率。然而Meta-IS-AK算法抽取的失效域内的样本点的密度函数为最优重要抽样密度函数，因而无法直接用于估计输入变量失效域内的条件概率密度函数，为此第3步则对失效域内服从于重要抽样密度函数的样本点进行转换，该转换将利用Metropolis-Hastings准则由重要抽样密度函数的失效样本点来得到原概率密度函数在失效域内的样本点，进而求得各个输入变量在失效域中的条件概率密度函数。

本文主要由以下几部分组成。第1节简要地回顾了已有的利用贝叶斯公式转换的可靠性全局灵敏度的形式，第2节详细介绍了本文所提的计算可靠性全局灵敏度的高效算法。第3节分别给出了数值与工程算例，验证了所提算法的高效性，第4节则得出结论。

1 基于贝叶斯公式的可靠性全局灵敏度指标的转换形式

(1)

利用贝叶斯公式，可将条件失效概率Pf|Xi进行如下转换[8-9,12,23-24]：

(2)

将式(2)代入到式(1)中，即可得到基于贝叶斯公式的可靠性全局灵敏度指标表达式：

(3)

2 求解可靠性全局灵敏度指标的高效算法

2.1 失效概率Pf的估计

(4)

步骤1以gK1(x)作为gK2(x)的初始模型，即令gK2(x)=gK1(x)。

步骤3计算g(xu)，将(xu,g(xu))添加到训练集T中，形成更新的训练集T={T∪(xu,g(xu))}。

步骤4由T更新gK2(x)，并返回步骤2。

(5)

(6)

(7)

2.2 失效域条件下输入变量的概率密度函数估计

步骤2计算下列比值r,并由r确定下一个马尔可夫链样本

i=1,2,…,NF-1

(8)

i=1,2,…,NF-1

(9)

式中：u为[0,1]区间上服从均匀分布的随机数。

2.3 Meta-IS-AK算法求解可靠性全局灵敏度指标流程图

由2.1和2.2节的内容可给出Meta-IS-AK算法求解可靠性全局灵敏度的流程如图1所示。

图1 Meta-IS-AK算法求解的流程图Fig.1 Flowchart of estimating Meta-IS-AK algorithm

3 算例分析

为了验证所提算法在求解可靠性全局灵敏度的效率和精度，本节给出了3个算例，分别采用准Monte Carlo算法(QMC)、AK-MCS算法、Meta-IS算法及Meta-IS-AK算法求解可靠性全局灵敏度指标。其中QMC算法结果将作为参照解。

3.1 不连续多失效域算例

考虑一个多失效域的算例，其功能函数为

(10)

(11)

图2 算例3.1失效边界Fig.2 Failure boundary of Example 3.1

表1 算例3.1可靠性全局灵敏度指标的计算结果

Table 1Results of reliability global sensitivity index for Example 3.1

算法X1/10-3X2/10-3NcallQMC14.5[0.30]14.5[0.46]1×105AK-MCS14.2[0.23]14.1[0.35]440Meta-IS14.2[0.41]14.0[0.65]42+5 000Meta-IS-AK14.2[0.43]13.9[0.65]42+254

注：数据右上标为可靠性全局灵敏度指标估计值的标准差,Ncall为功能函数的调用次数。

3.2 钢架结构算例

对于一钢架结构，有4个失效模式：

(12)

因为是串联系统，所以系统的功能函数g为g=min{g1,g2,g3,g4}。在各个失效模式的功能函数中，Mi(i=1,2,3)与S是独立的，均值和标准差分别为:μMi=2(i=1,2,3),μS=1,σMi=2(i=1,2,3),σS=0.25。各输入变量所对应的可靠性全局灵敏度指标的计算结果如表2所示。

表2 算例3.2可靠性全局灵敏度指标的计算结果

由表2中结果可以看出，Meta-IS-AK算法得到的计算结果与QMC算法计算所得结果基本一致，且模型的调用次数低于QMC算法与AK-MCS算法以及Meta-IS算法。这说明Meta-IS-AK算法能够很好地适用于多个失效模式的情况。对于此算例，输入变量的重要性排序为S>M3>M1>M2。此外，输入变量S对失效概率的影响较其他输入变量的影响大很多。因此在进行可靠性设计时通过调整输入变量S的不确定性可以更有效地满足可靠性要求。此外，M2对失效概率影响的灵敏度指标较其他的输入变量小很多，因此可以忽略M2的不确定性对于失效概率的影响，可以通过将M2固定在均值处来简化失效概率的求解模型。

3.3 简化的翼盒结构算例(多设计点算例)

图3(a)是简化翼盒模型的示意图。这种翼盒结构由64个杆和42个板组成。64个杆根据它们的方向分为3组，x、y和z方向上杆的长度分别是2L、L和 3L。所有杆的截面积均为A，所有板的厚度均为TH，E和P分别是所有板和杆的弹性模量和外部载荷。泊松比为0.3。假设输入变量是独立的正态变量，其分布参数如表3所示。

图3 翼盒结构Fig.3 Wing box structure

表3 算例3.3输入变量的分布参数

Table 3Distribution parameters of input variables of Example 3.3

变量分布均值变异系数A/m2正态1×10-40.1L/m正态0.20.1E/Pa正态7.1×10100.12P/N正态1 5000.1TH/m正态2.5×10-30.15

由于翼盒结构承担了来自机翼的大部分垂直方向的载荷，其自由端的位移可能相对较大。大的位移可能会引起杆的变形并导致杆或板的破坏，这在飞行过程中是不允许的。因此计算翼盒结构的位移是很有必要的，且位移的计算可以通过有限元分析完成。图3(b)显示了在ANSYS 14.0中构建的翼盒结构的有限元模型，其中输入变量固定在它们的均值处，此外，翼盒的变形如图3(c)所示。

本文中，该翼盒的隐式功能函数是根据翼盒的最大位移Dmax不超过临界值构造的(这里临界值等于0.01)，其表达式为

g(X)=0.01-|Dmax|

(13)

式中：Dmax=D(L,A,E,P,TH)，是关于输入变量的隐式函数，由有限元分析确定。算例3.3的可靠性全局灵敏度计算结果如图4所示，计算结果的标准差如表4所示。

图4 算例3.3可靠性全局灵敏度指标的计算结果Fig.4 Results of reliability global sensitivity index for Example 3.3

表4 算例3.3可靠性全局灵敏度指标结果

Table 4Reliability global sensitivity index forExample 3.3

变量SDi/10-5QMCAK-MCSMeta-ISMeta-IS-AKA/m20.270.200.170.18L/m0.260.110.140.11E/Pa5.363.875.045.05P/N0.350.420.330.32TH/m8.178.856.376.30Ncall5×10512033+5 00033+59

由表4中的结果可以看出，Meta-IS-AK算法得到的计算结果精度较高且模型的调用次数低于QMC算法与AK-MCS算法以及Meta-IS算法，这说明了Meta-IS-AK算法能够很好地适用于隐式功能函数的情况。对于算例3.3，输入变量的的不确定性对于失效概率影响的重要性排序为TH>E>P>L>A。从图4的结果可以看出，输入变量TH对失效概率的影响较其他输入变量的影响大很多，因此在进行可靠性设计时调整输入变量TH的不确定性可以更有效地满足可靠性要求。从图4的结果还可以看出，L和A的不确定性对于失效概率的影响较其他的输入变量小很多，因此忽略L和A的不确定性对于失效概率的影响可以简化失效概率的求解模型。

4 结 论

1) 可靠性全局灵敏度指标可以有效地衡量输入变量对于结构失效概率的影响。为了高效地计算可靠性全局灵敏度指标，本文提出了可靠性全局灵敏度指标计算的Meta-IS-AK算法。

2) Meta-IS-AK算法在计算可靠性全局灵敏度指标时充分利用了可靠性全局灵敏度指标贝叶斯算法的维度独立性，也使得该指标的求解转化为输出样本的分类问题，这就使可靠性全局灵敏度指标可以通过嵌入式代理模型结合数字模拟法来求解。相比于传统的重要抽样函数，元重要抽样法构造的准最优抽样密度函数能够很好地适用于多设计点和多失效域的情况。在此基础上构造重构的功能函数的代理模型能够准确地对准重要抽样密度函数抽取的样本进行分类，进一步地提高了失效概率的计算效率。M-H准则无需额外调用功能函数便可以将准重要抽样密度函数的失效样本转化成原概率密度函数的失效样本，大大提高了求解输入变量的失效条件概率密度函数的效率。

3) 所提方法能够很好地适用于多设计点和多失效域的情况，也适用于隐函数问题的可靠性灵敏度分析。所提方法的计算效率远远高于QMC算法，同时较AK-MCS算法以及Meta-IS算法也有提高。所给算例充分验证了以上结论。