南京航空航天大学学报  2016, Vol. 48 Issue (5): 705-713   PDF    
可靠性及可靠性灵敏度分析的改进点估计方法
张永利, 吕震宙, 李生兰     
西北工业大学航空学院,西安,710072
摘要: 针对低维高度非线性问题,提出了一种可靠性及可靠性灵敏度分析的改进的点估计方法。基本思想是首先由空间分割来降低局部子空间中功能函数的非线性程度,然后用低精度的稀疏网格积分探索子空间中功能函数的概率响应特性,最后组合子空间中的信息来得到所需的可靠性及其灵敏度分析结果。方法的优点是适用于低维高度非线性功能函数的失效概率和可靠性灵敏度分析,由于无需求解功能函数的梯度函数,因此适用于复杂的隐式功能函数。另外,由于方法利用少量均匀抽样来估计对失效概率贡献最大的点,并依据其进行子空间的划分,从而使得子空间的划分更有利于提高可靠性及可靠性灵敏度分析的效率。用算例对所提方法进行了验证,结果表明:在低维高度非线性条件下,所提算法的精度和效率比类似的三点估计、直接稀疏网格积分方法有明显优势。
关键词: 可靠性     可靠性灵敏度     非线性     空间分割     稀疏网格     点估计    
Improved Point Estimation Method for Analyzing Reliability and Reliability Sensitivity
Zhang Yongli, Lü Zhenzhou, Li Shenglan     
School of Aeronautics, Northwestern Polytechnical University, Xi′an, 710072, China
Abstract: For low-dimensional and high-nonlinear performance function, an improved point estimation method is proposed for analyzing reliability and reliability sensitivity. In the proposed method, input space is firstly separated into some subspaces to reduce the nonlinearity of performance function in these subspaces. Secondly, low-level sparse grid integration method is used to estimate probabilistic response character in the subspaces. Finally, the probabilistic response characters in every subspace are combined to obtain the reliability and reliability sensitivity. The obvious advantage of the proposed method is its applicability for the reliability and reliability sensitivity of low-dimensional and high-nonlinear model, and it is also adaptive to complex implicit function because it is gradient free. Moreover, input space is partitioned according to the most probable point estimated by uniformly sampling few samples, which helps to improve the estimation efficiency of the reliability and reliability sensitivity. Several test examples demonstrate that for low-dimensional and high-nonlinear model the proposed method is more precise and efficient than existing point estimation methods, i.e., the three-point estimation and the direct sparse grid integral method.
Key words: reliability     reliability sensitivity     nonlinear     space partition     sparse grid     point estimation    

在可靠性及可靠性灵敏度分析过程中,点估计方法由于不需要求功能函数的导函数和寻找设计点而得到了广泛的应用。该类方法直接利用功能函数在一些特征点处的函数值及权函数来近似计算功能函数的低阶矩(主要是一到四阶矩),然后由功能函数的各阶矩来近似失效概率[1]。基于函数矩的可靠性灵敏度分析则依据灵敏度的定义,将失效概率对基本变量分布参数的偏导数转化成功能函数的各阶矩对基本变量分布参数的偏导数的形式,再将此偏导数转化成特征函数矩的形式。常用的点估计方法有两点估计、三点估计及稀疏网格积分等.Rosenblueth首先提出了两点估计法[2-3],用以计算单变量或多变量函数三阶(包括三阶)以下的概率矩;随后Gorman推导出了考虑输入随机变量前四阶矩的三点估计法(Three-point estimate,TPE)[4];Seo和Kwak利用试验设计技术也得到了相同的计算公式[5];然而Gorman和Zhao在实际使用中发现了TPE的一些不足,例如当功能函数的非线性程度较高时,TPE的精度会很低[4, 6],因此该方法不适合高度非线性的情况;文献[7]通过用Nataf变换替代传统方法中的Rosenblatt变换,将上述点估计方法推广到输入变量相关情况下的前四阶矩求解当中。近年来,以Smolyak准则为基础的稀疏网格方法(Sparse grid,SG)被广泛地应用到数值积分[8, 9]、插值[10-11]、微分方程的求解[12]、区域重要性测度估计[13]以及随机不确定性的传递中,并且已经被证明是一种特别适用于高维情况的有效离散化方法,但对于高度非线性问题,由于功能函数的复杂性,往往需要采用高精度的稀疏网格,这将导致计算量,即模型调用次数大幅度增加。由此可见对于高度非线性问题,SG方法的效率仍然有待提高。

如果能够降低功能函数的复杂性和非线性程度,则可以利用SG方法来进行有效的可靠性及相应的灵敏度分析,而空间分割则可以有效降低子空间中功能函数的复杂性和非线性程度,因此采用空间分割与低精度稀疏网格相结合的方法来进行相应的可靠性及可靠性灵敏度分析。空间分割的关键是分割策略,合理的分割策略才可以提高可靠性及其灵敏度分析的效率。由于实际工程问题中的失效概率一般非常小,使用输入变量的联合概率密度函数进行等概率空间分割时,会导致失效域中的子空间数量远远小于在安全域中的子空间数量,这将使得求解可靠性及其灵敏度的效率非常低下。本文将利用文献[14]中的方法来近似求得用于空间划分的重要抽样密度函数,其思想是在输入变量可能的取值区域进行均匀抽样,将落入失效域中的联合概率密度函数值最大的点作为重要抽样密度函数的密度中心,从而构造一个近似优化的重要抽样密度函数,并且在功能函数具有复杂性和高度非线性时,能够有效避免“多值”问题导致的寻找优化重要抽样密度函数的困难。本文按照此优化后的重要抽样密度函数对输入空间进行分割,使得对失效概率贡献较大的区域的子空间划分得较为密集,对失效概率影响小的区域的子空间划分得比较稀疏,从而达到了节约计算成本的目的。

在实际的工程应用中所面对的往往是高维高度非线性的复杂问题。目前已经存在多种降维技术可以有效地解决高维问题,比如渐进空间法[15]等可以将任意的高维模型降至低维,然而这种降维的代价是非线性程度的进一步增大。因此解决低纬高度非线性问题就显得非常必要,这类问题的解决再结合有效的降维技术就有望解决工程中的高维高度非线性难题,这也是本文研究的背景和所提方法的价值所在。

1 可靠性及可靠性灵敏度分析的矩方法

对于功能函数g(x),x=(x1, x2, …, xn)为不确定性输入变量,其失效概率可以表示成积分的形式

$ {P_f}=\int_F {{f_X}\left(\mathit{\boldsymbol{x}} \right)} {\rm{d}}\mathit{\boldsymbol{x=}}\int_{{R^n}} {{I_F}{f_X}\left(\mathit{\boldsymbol{x}} \right){\rm{d}}\mathit{\boldsymbol{x}}} $ (1)

式中:fX(x)表示X的联合概率密度函数;F={x:g(x)≤0}为g(x)定义的失效域;${I_F}=\left\{ \begin{array}{l} 1, \mathit{\boldsymbol{x}} \in F\\ 0, \mathit{\boldsymbol{x}} \notin F \end{array} \right.$为失效域F的指示函数;Rn表示整体n维输入变量空间。在数学上可靠性灵敏度是由失效概率Pf对基本变量xi的分布参数θxi的偏导数予以表达的

$ \frac{{\partial {P_f}}}{{\partial {\theta _{{x_i}}}}}=\int_{{{\bf{R}}^n}} {{I_F}\left(x \right)\frac{{\partial {f_X}\left(\mathit{\boldsymbol{x}} \right)}}{{\partial {\theta _{{x_i}}}}}{\rm{d}}\mathit{\boldsymbol{x}}} $ (2)
1.1 功能函数矩的计算

功能函数g(x)的前四阶概率矩αmg(m=1, 2, 3, 4)为

$ \begin{array}{l} {\alpha _{1g}}=\int {g\left(\mathit{\boldsymbol{x}} \right){f_X}\left(\mathit{\boldsymbol{x}} \right)} {\rm{d}}\mathit{\boldsymbol{x}}\\ {\alpha _{2g}}={\left[{\int {{{\left({g\left(\mathit{\boldsymbol{x}} \right)-{\alpha _{1g}}} \right)}^2}{f_X}\left(\mathit{\boldsymbol{x}} \right){\rm{d}}\mathit{\boldsymbol{x}}} } \right]^{\frac{1}{2}}}\\ {\alpha _{mg}}\alpha _{2g}^m=\int {{{\left({g\left(\mathit{\boldsymbol{x}} \right)- {\alpha _{1g}}} \right)}^m}{f_X}\left(\mathit{\boldsymbol{x}} \right){\rm{d}}\mathit{\boldsymbol{x}}} \;\;m=3, 4 \end{array} $ (3)

TPE和稀疏网格积分法是常用的计算前四阶矩的点估计方法,TPE的具体内容见文献[4, 5],SG方法由于其高效性而得到了广泛的应用[8-12, 16],它的思想及实现过程可归纳如下:

S1ijw1ij表示第j个变量一维空间中的第ij个积分点和权重,则n维空间中k精度下所有稀疏网格积分点的集合Sn(k)由以下Smolyak准则选取

$ \mathit{\boldsymbol{S}}_n^{\left(k \right)}=\mathop \cup \limits_{k+1 \le \left| i \right| \le q} \mathit{\boldsymbol{S}}_{{1^1}}^i \otimes \mathit{\boldsymbol{S}}_{{1^2}}^i \otimes \cdots \otimes \mathit{\boldsymbol{S}}_{{1^n}}^i $ (4)

式中:ⓧ表示张量积计算;q=k+n;|i|=i1+i2+…+in为多维指标之和。依据Smolyak准则,集合Snk中相应于第l个积分点ξl=[ξi1ji1, …, ξinjin]∈Sn(k)的权重w(l)

$ {w^{\left(l \right)}}={\left({ - 1} \right)^{q - \left| i \right|}}\left(\begin{array}{l} n - 1\\ q - \left| i \right| \end{array} \right)\left({w_{{j_{{i_1}}}}^{{i_1}} \cdots w_{{j_{{i_n}}}}^{{i_n}}} \right) $ (5)

则对非线性函数g(x)的积分可以由式(6)求得,并且能够达到2k+1阶多项式精度[8]

$ \begin{array}{l} \int {g\left(\mathit{\boldsymbol{x}} \right){f_X}\left(\mathit{\boldsymbol{x}} \right)} {\rm{d}}\mathit{\boldsymbol{x}} \approx \sum\limits_{l=1}^{M_{^n}^{\left(k \right)}} {{w^{\left(l \right)}}} g\left({{T^{ - 1}}\left({{{\bf{ \pmb{\mathsf{ ξ}} }}^{\left(l \right)}}} \right)} \right)=\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\sum\limits_{l=1}^{M_{^n}^{\left(k \right)}} {{w^{\left(l \right)}}} g\left({{s^{\left(l \right)}}} \right) \end{array} $ (6)

式中:Mn(k)表示k阶精度下相应的积分点个数;T-1(·)为任意分布的变量x向积分点空间变换函数的反函数,且在积分点ξl处的函数值为sl,则SG方法求解功能函数g(x)的前四阶矩的积分[16]可以表示为

$ \begin{array}{l} {\alpha _{1g}}=\int {g\left(\mathit{\boldsymbol{x}} \right){f_X}\left(\mathit{\boldsymbol{x}} \right)} {\rm{d}}\mathit{\boldsymbol{x}} \approx \sum\limits_{l=1}^{M_{^n}^{\left(k \right)}} {{w^{\left(l \right)}}} g\left({{s^{\left(l \right)}}} \right)\\ \;\;\;\;\;{\alpha _{2g}}={\left[{\int {{{\left({g\left(\mathit{\boldsymbol{x}} \right)-{\alpha _{1g}}} \right)}^2}{f_X}\left(\mathit{\boldsymbol{x}} \right){\rm{d}}\mathit{\boldsymbol{x}}} } \right]^{\frac{1}{2}}} \approx \\ \;\;\;\;\;\;\;\;\;\;\;\;\;{\left[{\sum\limits_{l=1}^{M_{^n}^{\left(k \right)}} {{w^{\left(l \right)}}} {{\left({g\left({{s^{\left(l \right)}}} \right)-{\alpha _{1g}}} \right)}^2}} \right]^{\frac{1}{2}}}\\ {\alpha _{mg}}\alpha _{2g}^m=\int {{{\left({g\left(\mathit{\boldsymbol{x}} \right)- {\alpha _{1g}}} \right)}^m}{f_X}\left(\mathit{\boldsymbol{x}} \right){\rm{d}}\mathit{\boldsymbol{x}}} \approx \\ \sum\limits_{l=1}^{M_{^n}^{\left(k \right)}} {{w^{\left(l \right)}}} {\left({g\left({{s^{\left(l \right)}}} \right)- {\alpha _{1g}}} \right)^k}\;\;\;m=3, 4 \end{array} $ (7)

基于Smolyak准则的稀疏网格能够在一定程度上克服传统数值积分计算量随变量维数呈指数级增长的维数灾难问题,对高维积分问题有很好的适用性,通过调整精度k可以方便地提高积分精度。但是对于非线性程度很高的模型,即使提高精度k也未必可以达到理想的精度。

1.2 基于前四阶矩的失效概率近似方法

基于功能函数概率矩近似失效概率的四阶矩方法包含高阶矩标准化技术法(High-order moment standardization technique,HO M ST)[17]、Edgeworth级数法[17]和频率曲线法[17-21]等。本文采用HOMST方法,可靠度指标β4M和失效概率Pf分别为

$ {\beta _{4M}}=\frac{{3\left({{\alpha _{4g}} - 1} \right){\beta _{2M}}+{\alpha _{3g}}\left({\beta _2^2M - 1} \right)}}{{\sqrt {\left({5\alpha _{3g}^2 - 9{\alpha _{4g}}+9} \right)\left({1 - {\alpha _{4g}}} \right)} }} $ (8)
$ {P_f}=\Phi \left({ - {\beta _{4M}}} \right) $ (9)

式中:${\beta _{2M}}=\frac{{{\alpha _{1g}}}}{{{\alpha _{2g}}}}$,且当α3g=0时,β4M=β2M

1.3 可靠性灵敏度分析的点估计方法

根据可靠性灵敏度的定义,由失效概率Pf与可靠度指标β4M的关系,以及β 4Mαmgm=1, …, 4的关系,可以采用函数求导法则推出考虑前四阶矩时Pf对基本变量Xi的第r个分布参数θ(r)xi的灵敏度为

$ \begin{array}{*{20}{c}} {\frac{{\partial {P_f}}}{{\partial \theta _{{x_i}}^{\left(r \right)}}}=\frac{{\partial {P_f}}}{{\partial {\beta _{4M}}}}\frac{{\partial {\beta _{4M}}}}{{\partial \theta _{{x_i}}^{\left(r \right)}}}=}\\ {\frac{{\partial {P_f}}}{{\partial {\beta _{4M}}}}\left({\frac{{\partial {\beta _{4M}}}}{{\partial {\beta _{2M}}}}\left({\frac{{\partial {\beta _{2M}}}}{{\partial {\alpha _{1g}}}}\frac{{\partial {\alpha _{1g}}}}{{\partial \theta _{{x_i}}^{\left(r \right)}}}+\frac{{\partial {\beta _{2M}}}}{{\partial {\alpha _{2g}}}}\frac{{\partial {\alpha _{2g}}}}{{\partial \theta _{{x_i}}^{\left(r \right)}}}} \right)} \right.}\\ {\left.{\frac{{\partial {\beta _{4M}}}}{{\partial {\alpha _{3g}}}}\frac{{\partial {\alpha _{3g}}}}{{\partial \theta _{{x_i}}^{\left(r \right)}}}+\frac{{\partial {\beta _{4M}}}}{{\partial {\alpha _{4g}}}}\frac{{\partial {\alpha _{4g}}}}{{\partial \theta _{{x_i}}^{\left(r \right)}}}} \right)} \end{array} $ (10)

式中:${\frac{{\partial {P_f}}}{{\partial {\beta _{4M}}}}}$${\frac{{\partial {\beta _{4M}}}}{{\partial {\beta _{2M}}}}}$${\frac{{\partial {\beta _{4M}}}}{{\partial {\alpha _{3g}}}}}$${\frac{{\partial {\beta _{4M}}}}{{\partial {\alpha _{4g}}}}}$以及αkgθkxi的偏导数${\frac{{\partial {\alpha _{mg}}}}{{\partial \theta _{{x_i}}^{\left(r \right)}}}}$的具体表达形式在文献[1]中有详细说明,此处不再赘述。

2 基于空间分割的改进点估计法 2.1 基本思想

基于空间分割的点估计方法通过空间分割后的子空间来降低功能函数的复杂性和非线性程度,从而达到提高计算效率的目的,其基本思想是:首先利用文献[14]中所介绍的方法使用少量样本近似最可能失效点(Most probable point,MPP),并以该点为密度中心构造近似优化的重要抽样密度函数hX(x);然后按照hX(x)对输入空间进行划分(Space parti tion,SP),进而在每个子空间中使用低精度的SG积分法探索功能函数的概率响应特性;最后组合所有子空间中的信息得到所需的可靠性及可靠性灵敏度结果。本文以下部分用SG-SP表示本文方法。

2.2 空间划分策略及划分后的失效率概率和可靠性灵敏度

空间划分策略总体上可分为两步:

步骤1    借鉴文献[14]中提出的MPP近似方法,构造重要抽样密度函数hX(x)。

(1)对系统的可靠度指标作出保守的假设,然后根据MPP估计中生成样本的最小区间(表 1)使用输入变量的累积分布函数的反函数确定“合适区间”的边界;

表 1 MPP估计中生成样本的最小区间 Table 1 The least required interval for sample generation in MPP estimation

(2)在这些“合适区间”内均匀抽取NMPP个样本xjj=1, …, NMPP);

(3)对于n维独立随机变量X,根据Wj=$\prod\limits_{i=1}^n {{f_{{X_i}}}\left({{x_{i \cdot j}}} \right)} $分别计算这些样本的权重系数W=(W1, …,WNMPP),其中fxi表示第i维变量的边缘概率密度函数;

(4)失效域中权重系数最大的样本点即为MPP的近似值;

(5)将输入变量联合密度函数fX(x)的中心平移至MPP处,形成重要抽样密度函数hX(x)。

步骤2    根据hX(x)等概率划分输入变量空间,生成子空间。

这种划分策略的结果是:对失效概率贡献较大的地方空间分割比较密集,在影响小的地方则较为稀疏,提高了后续失效概率计算的效率。

根据式(1, 2),失效概率和可靠性灵敏度在空间分割后可以表示为

$ \begin{array}{*{20}{c}} {{P_f} = \int {{I_F}\left( \mathit{\boldsymbol{x}} \right){f_X}\left( \mathit{\boldsymbol{x}} \right)} {\rm{d}}\mathit{\boldsymbol{x = }}}\\ {\;\int {{I_F}\left( \mathit{\boldsymbol{x}} \right)\frac{{{f_X}\left( \mathit{\boldsymbol{x}} \right)}}{{{h_X}\left( \mathit{\boldsymbol{x}} \right)}}} {h_X}\left( \mathit{\boldsymbol{x}} \right){\rm{d}}\mathit{\boldsymbol{x}} = }\\ {\sum\limits_{{i_1} = 1}^{{N_1}} {\sum\limits_{{i_2} = 1}^{{N_2}} { \cdots \sum\limits_{{i_n} = 1}^{{N_n}} {\int {_{{U_{{i_1}}}}\int {_{{U_{i2}}} \cdots \int {_{{U_{{i_n}}}}{I_F}\left( \mathit{\boldsymbol{x}} \right)\frac{{{f_X}\left( \mathit{\boldsymbol{x}} \right)}}{{{h_X}\left( \mathit{\boldsymbol{x}} \right)}}} } } } } } }\\ {\frac{{{h_X}\left( \mathit{\boldsymbol{x}} \right)}}{{\int {_{{U_{{i_1},}}{i_2}, \cdots ,{i_n}}} {h_X}\left( \mathit{\boldsymbol{t}} \right){\rm{d}}\mathit{\boldsymbol{t}}}}\left( {\int {_{{U_{{i_1},}}{i_2}, \cdots ,{i_n}}} {h_X}\left( \mathit{\boldsymbol{t}} \right){\rm{d}}\mathit{\boldsymbol{t}}} \right){\rm{d}}\mathit{\boldsymbol{x}} = }\\ {\sum\limits_{{i_1} = 1}^{{N_1}} {\sum\limits_{{i_2} = 1}^{{N_2}} { \cdots \sum\limits_{{i_n} = 1}^{{N_n}} {{P_{_{{i_1},}{i_2}, \cdots ,{i_n}}}{E_{_{{i_1},}{i_2}, \cdots ,{i_n}}}\left( {{I_F}\left( \mathit{\boldsymbol{x}} \right)\frac{{{f_X}\left( \mathit{\boldsymbol{x}} \right)}}{{{h_X}\left( \mathit{\boldsymbol{x}} \right)}}} \right)} } } } \end{array} $ (11)
$ \begin{array}{*{20}{c}} {\frac{{\partial {P_f}}}{{\partial \theta _{{x_i}}^{\left(r \right)}}}=\int {_{{R_n}}{I_F}\left(\mathit{\boldsymbol{x}} \right)\frac{{\partial {f_X}\left(\mathit{\boldsymbol{x}} \right)}}{{\partial \theta _{{x_i}}^{\left(r \right)}}}} {\rm{d}}\mathit{\boldsymbol{x=}}}\\ {\;\int {_{{R_n}}\frac{{{I_F}\left(\mathit{\boldsymbol{x}} \right)}}{{{f_X}\left(\mathit{\boldsymbol{x}} \right)}}\frac{{\partial {f_X}\left(\mathit{\boldsymbol{x}} \right)}}{{\partial \theta _{{x_i}}^{\left(r \right)}}}\frac{{{f_X}\left(\mathit{\boldsymbol{x}} \right)}}{{{h_X}\left(\mathit{\boldsymbol{x}} \right)}}} {h_X}\left(\mathit{\boldsymbol{x}} \right){\rm{d}}\mathit{\boldsymbol{x}}=}\\ {\sum\limits_{{i_1}=1}^{{N_1}} {\sum\limits_{{i_2}=1}^{{N_2}} { \cdots \sum\limits_{{i_n}=1}^{{N_n}} {\int {_{{U_{{i_1}}}}\int {_{{U_{i2}}} \cdots \int {_{{U_{{i_n}}}}\frac{{{I_F}\left(\mathit{\boldsymbol{x}} \right)}}{{{h_X}\left(\mathit{\boldsymbol{x}} \right)}}\frac{{\partial {f_X}\left(\mathit{\boldsymbol{x}} \right)}}{{\partial \theta _{{x_i}}^{\left(r \right)}}}} } } } } } }\\ {\frac{{{h_X}\left(\mathit{\boldsymbol{x}} \right)}}{{\int {_{{U_{{i_1}, }}{i_2}, \cdots, {i_n}}} {h_X}\left(\mathit{\boldsymbol{t}} \right){\rm{d}}\mathit{\boldsymbol{t}}}}\left({\int {_{{U_{{i_1}, }}{i_2}, \cdots, {i_n}}} {h_X}\left(\mathit{\boldsymbol{t}} \right){\rm{d}}\mathit{\boldsymbol{t}}} \right){\rm{d}}\mathit{\boldsymbol{x}}=}\\ {\sum\limits_{{i_1}=1}^{{N_1}} {\sum\limits_{{i_2}=1}^{{N_2}} { \cdots \sum\limits_{{i_n}=1}^{{N_n}} {{P_{_{{i_1}, }{i_2}, \cdots, {i_n}}}{E_{{h_{_{{i_1}, }{i_2}, \cdots, {i_n}}}}}\left({\frac{{{I_F}\left(\mathit{\boldsymbol{x}} \right)}}{{{h_X}\left(\mathit{\boldsymbol{x}} \right)}}\frac{{\partial {f_X}\left(\mathit{\boldsymbol{x}} \right)}}{{\partial \theta _{{x_i}}^{\left(r \right)}}}} \right)} } } } \end{array} $ (12)

式中:Uij(ij=1, 2, …, Nj, j=1, 2, …, n)表示第j个变量的第ij个小区间;Nj表示第j个变量被划分的总区间数;hi1, i2, …, in(x)=${\frac{{{h_X}\left(\mathit{\boldsymbol{x}} \right)}}{{\int {_{{U_{{i_1}, }}{i_2}, \cdots, {i_n}}} {h_X}\left(\mathit{\boldsymbol{t}} \right){\rm{d}}\mathit{\boldsymbol{t}}}}}$表示基本变量X1, X2, …,Xn的第i1, i2, …,in个小区间形成的输入变量子空间Ui1, i2, …, in=Ui1Ui2∩…∩Uin的联合概率密度函数;Ehi1, i2, …, in(·)表示以hi1, i2, …, in(x)为密度函数所求得的期望;Pi1, i2, …, in=∫U i1, i2,…,inhX(t) dt

2.3 SG-SP方法的实现步骤及计算量

结合前面介绍的SG方法,可以按照以下步骤估计失效概率和可靠性灵敏度:

步骤1    根据近似优化的重要抽样概率密度函数hX(x)将输入空间等概率划分成若干个子空间。

步骤2   在子空间Ui1, i2, …, in中,产生M(k)n个稀疏网格积分点s(l)i1, i2, …, in(l=1, 2, …, Mnk)和对应的权重系数w(l)i1, i2, …, in

步骤3   按照式(13, 14)计算Pf$\frac{{\partial {P_f}}}{{\partial \theta _{{x_i}}^{\left(k \right)}}}$

$ \begin{array}{*{20}{c}} {{P_f}=\sum\limits_{{i_1}=1}^{{N_1}} {\sum\limits_{{i_2}=1}^{{N_2}} { \cdots \sum\limits_{{i_n}=1}^{{N_n}} {{P_{_{{i_1}, }{i_2}, \cdots, {i_n}}} \cdot } } } }\\ {\left({\sum\limits_{l=1}^{M_{^n}^{\left(k \right)}} {w_{_{{i_1}, }{i_2}, \cdots, {i_n}}^{\left(l \right)}{I_F}\left({s_{_{{i_1}, }{i_2}, \cdots, {i_n}}^{\left(l \right)}} \right)\frac{{{f_X}\left({s_{_{{i_1}, }{i_2}, \cdots, {i_n}}^{\left(l \right)}} \right)}}{{{h_X}\left({s_{_{{i_1}, }{i_2}, \cdots, {i_n}}^{\left(l \right)}} \right)}}} } \right)} \end{array} $ (13)
$ \begin{array}{*{20}{c}} {\frac{{\partial {P_f}}}{{\partial \theta _{{x_i}}^{\left(r \right)}}}=\sum\limits_{{i_1}=1}^{{N_1}} {\sum\limits_{{i_2}=1}^{{N_2}} { \cdots \sum\limits_{{i_n}=1}^{{N_n}} {{P_{_{{i_1}, }{i_2}, \cdots, {i_n}}} \cdot } } } }\\ {\left({\sum\limits_{l=1}^{M_{^n}^{\left(k \right)}} {w_{_{{i_1}, }{i_2}, \cdots, {i_n}}^{\left(l \right)}\frac{{{I_F}\left({s_{_{{i_1}, }{i_2}, \cdots, {i_n}}^{\left(l \right)}} \right)}}{{_X\left({s_{_{{i_1}, }{i_2}, \cdots, {i_n}}^{\left(l \right)}} \right)}}\frac{{\partial {f_X}\left({s_{_{{i_1}, }{i_2}, \cdots, {i_n}}^{\left(l \right)}} \right)}}{{\partial \theta _{{x_i}}^{\left(r \right)}}}} } \right)} \end{array} $ (14)

从式(13, 14)可以看出,SG-SP方法的总计算量为$M_n^{\left(k \right)} \cdot \prod\limits_{i=1}^n {{N_i}+{N_{{\rm{MPP}}}}} $,式中NMPP表示用来估计MPP的样本量。尽管每个子空间中SG方法产生的积分点数M(n)k比较小,但是总的计算量却与维数成指数关系,所以该方法只适合低维的情况。

3 算例

在以下算例中,将蒙特卡洛(Monte Carlo,MC)数字模拟方法的结果作为参考解,分别使用TPE,SG积分方法和本文SG-SP方法来估计失效概率和可靠性灵敏度,从计算结果精度和总模型调用次数两方面来比较。

在使用SG-SP方法计算时,为了方便比较,对同一功能函数g(x),假定每一维输入变量Xj划分段数Nj相同,即N1=N2=…=Nn=N。同时为了充分证明本文所提方法的优势,在使用SG积分方法计算前四阶矩时,精度k1取结果稳定时的最小值。根据式(8~10),这些稳定的结果用来计算Pf${\frac{{\partial {P_f}}}{{\partial \theta _{{x_i}}^{\left(r \right)}}}}$,而${\frac{{\partial {P_f}}}{{\partial \theta _{{x_i}}^{\left(r \right)}}}}$包含的其他积分运算由精度为k2的稀疏网格估计, Pf不随k2变化。

3.1 算例1

假设功能函数为g(x1, x2)=a+0.004 63·(x1+x2-20)q-0.235 7(x1-x2),其中独立的输入随机变量x服从正态分布,即xi~N(10, 3)(i=1, 2)。功能函数中指数q,常数a和使用SG-SP方法估计MPP时均匀抽取的样本量NMPP表 2所示。

表 2 q,aNMPP的取值情况 Table 2 Values of q, a and NMPP

aq取不同值时前四阶中心矩和失效概率及可靠性灵敏度的计算结果见表 3~10,其中计算量表示总的模型函数调用次数。

表 3 q=1, N=17时输出量的前四阶中心矩 Table 3 First four central moments of output under q=1, N=17

表 4 q=1, N=17时可靠性及可靠性灵敏度分析结果 Table 4 Reliability and reliability sensitivity analysis results under q=1, N=17

表 5 q=2, N=13时输出量的前四阶中心矩 Table 5 First four central moments of output under q=2, N=1

表 6 q=2, N=13时可靠性及可靠性灵敏度分析结果 Table 6 Reliability and reliability sensitivity analysis results under q=2, N=13

表 7 q=3, N=13时输出量的前四阶中心矩 Table 7 First four central moments of output under q=3, N=13

表 8 q=3, N=13时可靠性及可靠性灵敏度分析结果 Table 8 Reliability and reliability sensitivity analysis results under q=3, N=13

表 9 q=4, N=12时输出量的前四阶中心矩 Table 9 First four central moments of output under q=4, N=12

表 10 q=4, N=12时可靠性及可靠性灵敏度分析结果 Table 10 Reliability and reliability sensitivity analysis results under q=4, N=12

结果分析:q=1和q=2时,功能函数的非线性程度比较低,SG方法在计算量小于SG-SP方法的最低计算量时已经得到了精度略高于后者的结果,在这种情况下SG方法无疑更具有优势。TPE方法在q=1时具有和SG方法相当的计算效率,但是在q=2时该方法的部分计算结果出现较大的偏差。

q=3和q=4时的计算结果来看,随着非线性程度的增加,在前四阶矩的计算中,TPE方法的计算结果已经部分错误,而SG方法的精度高于TPE方法。在表 810中,两种对照方法的结果几乎完全错误,即使提高SG方法的精度水平,也无法得到与MC解接近的结果。而SG-SP方法所得结果随着k的增加而逐渐收敛,且k=1时就已经达到了可用的精度。q=4时的情况同时也表明了所提方法在失效概率数量级较小时的良好适用性。需要说明的是,当q=4时,由于模型的可靠性及可靠性灵敏度过小,使用MC方法无法得到稳定的结果,并且前面3种情况已经证明了SG-SP方法的计算结果随着Nk的增加收敛于真实解,所以此情况下可以取N=40, k=8时的计算结果作为参考解。

q=4时的情况为例,简要说明SG-SP方法中Nk的选取方法。显然,为了获得要求的计算精度,可以通过增大Nk的取值来实现。N的取值越大,则在划分的子空间中功能函数的复杂性和非线性程度越低,这时就可以采用较小的k值来达到要求的精度。相反,N的取值越小,则要达到要求的精度就需要采用较大的k值来实现。一般来说k值较小时的结果较稳定。图 1给出了本算例在k取1,2,3和4时计算失效概率的均方误差随总计算次数变化的曲线,从图 1中可以看出:在相同计算量下,k=1时的均方误差最小,而相同均方误差情况下,k=1时的总计算量最小。因此一般情况下,在SG-SP方法中建议k取1或2,然后通过增加N来达到要求的计算精度。

图 1 均方误差与计算量的关系曲线 Figure 1 Mean square error of failure probability variying with computation

3.2 算例2

在汽车工程领域中,前轴承载着车辆前部的重力,所以必须保证其有足够的可靠性。工字梁由于具有良好的抗弯性能和质量较轻的特点而被广泛地应用于前轴的设计中。图 2为汽车前轴示意图,工字梁危险截面形状如图 3所示。已知危险截面的最大正应力为σ=M/Wxτ=T/Wρ,其中M=2.5×106 N·mm和T=2.1×106 N·mm分别为前轴所受的弯矩和扭矩,W xWρ分别为结构的截面系数和极截面系数,且有

$ {W_x}=\frac{{{{\left({h - 2t} \right)}^3}}}{{500h}}+\frac{b}{{6h}}\left[{{h^3}-{{\left({h-2t} \right)}^3}} \right] $ (15)
$ {W_\rho }=0.8b{t^2}+\frac{{100}}{3}{t^3}\left({h - 2t} \right) $ (16)
图 2 汽车前轴示意图 Figure 2 Diagram of front axle

图 3 工字梁截面 Figure 3 Cross section of front axle

根据前轴的材料特性,给定前轴的静强度极限为σs=600 MPa,故由静强度分析并构件模型的功能函数为g=${\sigma _s} - \sqrt {{\sigma ^2}+3{\tau ^2}} $,工字梁的几何参数b,t,h为相互独立的正态输入变量,其分布参数如表 11所示。采用TPE和SG方法作为对照,表 12为这两种方法对前四阶中心矩的估计结果。使用SG-SP方法对该结构进行可靠性及可靠性灵敏度分析时,将变量空间分割成83=512个子空间,均匀抽取400个样本估计MPP。失效概率及可靠性灵敏度的计算结果如表 13所示。

表 11 输入变量的分布参数 Table 11 Distribution parameters of input variables

表 12 输出量的前四阶中心矩 Table 12 First four central moments of output

表 13 可靠性及可靠性灵敏度分析结果 Table 13 Reliability and reliability sensitivity analysis results

结果分析:从以上计算结果可以发现两种对照方法SG和TPE估计的前四阶矩与参考值差别不大,但是基于前四阶矩的可靠性及可靠性灵敏度分析结果却几乎全部错误,而使用本文中所建立的SG-SP方法在k=2时已经具有良好精度,计算量为3 984。可靠性灵敏度反映了基本变量分布参数对失效概率的影响程度,无量纲正则化的可靠性灵敏度可以用于基本变量分布参数对可靠度的重要性排序。为了把失效概率控制在理想的范围,研究者可以依据这种重要性排序来决定哪些变量应该被重点考虑,这在工程上具有重要意义。

4 结束语

本文将空间分割思想,SG积分方法和最可能失效点估计方法结合在一起,提出了一种适合低维高度非线性模型可靠性及可靠性灵敏度分析的新方法。常用的传统点估计方法在非线性程度逐渐增加的情况下所计算的结果的误差也逐渐增大,甚至出现了严重的错误,而SG-SP方法却使用不大的计算量即可获得稳定的结果。需要再次强调的是在高维情况下,该方法的计算量会随着维数呈指数增加,在非线性程度低的情况下,与传统的点估计方法相比并不占有优势,但是在低维度、非线性程度较高的问题上本文方法体现出了明显的高效率和高精度的优势,与降维技术相结合将有望解决高维高度非线性问题。在SG-SP方法中,有关计算量随着维数呈指数增长问题将有望采用自适应空间剖分的方法加以解决,这将是作者后续的研究工作。

参考文献
[1] 吕震宙, 宋述芳, 李洪双, 等. 结构机构可靠性及可靠性灵敏度分析[M]. 北京: 科学出版社, 2009.
Lü Zhenzhou, Song Shufang, Li Hongshuang, et al. The reliability and reliability sensitivity analysis for structure and machine system[M]. Beijing: Science Press, 2009.
[2] Rosenblueth E. Point estimation for probability moments[J]. Proceedings of the National Academy of Science, 1975, 72(10): 3812–3814. DOI:10.1073/pnas.72.10.3812
[3] Rosenblueth E. Two-point estimates in probability[J]. Applied Mathematical Modeling, 1981, 5(5): 329–335. DOI:10.1016/S0307-904X(81)80054-6
[4] Gorman M. Reliability of structural systems [D]. Cleveland: Case Western Reserve University, 1980: 320-332.
[5] Seo H S, Kwak B M. Efficient statistical tolerance analysis for general distributions using three-point information[J]. International Journal of Production Research, 2002, 40(4): 931–944. DOI:10.1080/00207540110095709
[6] Zhao Yangang, Ono T. New point estimates for probability moments[J]. Journal of Engineering Mechanics, 2000, 126(4): 433–436. DOI:10.1061/(ASCE)0733-9399(2000)126:4(433)
[7] Li Hongshuang, Lü Zhenzhou, Yuan Xiukai. Nataf transformation based point estimate method[J]. Chinese Science Bulletin, 2008, 53(17): 2586–2592.
[8] Gerstner T, Griebel M. Numerical integration using sparse grids[J]. Number Algorithms, 1998, 18(01): 209–232.
[9] Novak E, Ritter K. High dimensional integration of smooth functions over cubes[J]. Numerische Mathematik, 1996, 75(1): 79–97. DOI:10.1007/s002110050231
[10] Klimke A. Sparse grid interpolation toolbox-User′s guide[R]. IANS Report 2006/001, 2006.
[11] Barthelmann V, Novak E, Ritter K. High dimensional polynomial interpolation on sparse grid[J]. Advances in Computational Mathematics, 2000, 12(4): 273–288. DOI:10.1023/A:1018977404843
[12] Xiu D B, Jan S H. High-order collocation methods for differential equations with random inputs[J]. SIAM Journal on Scientific Computing, 2005, 27(3): 1118–1139. DOI:10.1137/040615201
[13] 李璐祎, 吕震宙. 基本变量区域重要性测度及其稀疏网格解[J]. 力学学报, 2013, 45(4): 569–579.
Li Luyi, Lü Zhenzhou. Regional importance measure of the basic variable and its sparse grid solution[J]. Chinese Journal of Theoretical and Applied Mechanics, 2013, 45(4): 569–579.
[14] Mohsen R, Mahmoud M. A new efficient simulation method to approximate the probability of failure and most probable point[J]. Structural Safety, 2012, 39(4): 22–29.
[15] Luo Xiaopeng, Lü Zhenzhou, Xu Xin. Reproducing kernel technique for high dimensional model representations (HDMR)[J]. Computer Physics Communications, 2014, 185(12): 3099–3108. DOI:10.1016/j.cpc.2014.07.021
[16] Xiong Fenfen, Greene S, Wei Chen, et al. A new sparse grid based method for uncertainty propagation[J]. Structural and Multidisciplinary Optimization, 2010, 41(3): 335–349. DOI:10.1007/s00158-009-0441-x
[17] Zhao Yangang, Ono T. Moment methods for structural reliability[J]. Structural Safety, 2001, 23(1): 47–75. DOI:10.1016/S0167-4730(00)00027-8
[18] Zhao Yangang, Ono T. On the problems of the fourth moment method[J]. Structural Safety, 2004, 26(3): 343–347. DOI:10.1016/j.strusafe.2003.10.001
[19] Zhao Yangang, Lu Zhaohui. Applicable range of the fourth-moment method for structural reliability[J]. Journal of Asian Architecture and Building Engineering, 2007, 6(1): 151–158. DOI:10.3130/jaabe.6.151
[20] Zhao Yangang, Alfredo H S, Ang H M. System reliability assessment by method of moments[J]. Journal of Structural Engineering, 2003, 129(10): 1341–1349. DOI:10.1061/(ASCE)0733-9445(2003)129:10(1341)
[21] Hong H P. Point-estimate moment-based reliability analysis[J]. Civil Engineering System, 1996, 13(4): 281–294. DOI:10.1080/02630259608970204