结构可靠度求解的自适应细分-重要抽样法

时间：2024-07-28

王新愿周金宇谢里阳程锦翔

1.江苏理工学院机械工程学院，常州，213001 2.金陵科技学院机电工程学院，南京，211169 3.东北大学机械工程与自动化学院，沈阳，110819

0 引言

迄今为止，研究者们已提出很多可靠度求解方法。其中，基于设计点的方法是最为常见的方法[1-2]，该类方法将功能函数在设计点处进行泰勒展开，保留展开式一阶项或高阶项，进而计算结构可靠度，当随机变量呈明显非正态或功能函数高度非线性时，计算误差较大。数值模拟方法中最基本的是蒙特卡罗模拟(Monte Carlo simulation，MCS)方法，该方法具有较高的精度和鲁棒性，主要用于验证其他方法的有效性。这种方法在小失效概率的情况下，需要大量样本保证可靠度求解的精度，计算效率较低。

针对涉及非正态随机变量的可靠性问题，通常对随机变量当量正态化，这种转换会增加功能函数的非线性程度，降低可靠度求解的精度。鞍点近似(saddlepoint approximation，SA)方法可以直接评估具有非正态随机变量的功能函数的概率分布[3-4]，得到了广泛应用和发展。ZHANG等[5]将SA法与线抽样(line sampling，LS)法相结合，提出鞍点近似线采样(SA-LS)法，具有准确性高、效率高的优点。LU等[6]提出一种改进的基于高阶矩的SA法，提高了可靠度求解的精度。目前改进方法仅对单峰低偏度的类正态分布有效，而对高偏度、多峰分布或均匀分布存在较大误差或导致失效。针对涉及小失效概率的可靠性问题，国内外学者引入重要抽样(importance sampling，IS)和子集模拟(subset simulation，SS)方法[7-8]，取得很多研究成果。KURTZ等[9]提出交叉熵重要抽样(cross entro-py importance sampling，CE-IS)法，该方法与IS法相比，具有更高的抽样效率。史朝印等[10]结合自适应Kriging(adaptive Kriging，AK)代理模型，提出改进的CE-IS法，显著减小CE-IS法的计算量，提高了可靠度求解的效率。BOURINET等[11]结合支持向量机(support vector machine，SVM)和SS法，提高了可靠度求解的效率和准确性。HUANG等[12]将SS法与AK代理模型方法相结合，提出AK-SS方法，可以有效处理涉及小失效概率和高维功能函数的可靠性问题。针对涉及高度非线性功能函数的可靠性问题，传统的基于梯度的方法可能无法得到收敛解。ELEGBEDE[13]应用粒子群算法求解非线性结构失效概率，可得到较精确的结果。ZOU等[14]受粒子群算法的启发，提出一种全局协调搜索算法(NGHS)，具有更高的效率。YANG[15]引入混沌控制方法，提高了非线性结构可靠度求解的精度，但效率大幅降低。KESHTEGAR[16]引入共轭混沌控制方法，与混沌控制方法相比，提高了可靠度求解的效率。然而，对于同时涉及非正态随机变量、小失效概率以及高度非线性功能函数的可靠性问题，现有方法尚不能提供有效的解决方案。

通用生成函数(universal generating function，UGF)是建立在广泛使用概率论原理的生成函数基础上，用于枚举分析系统状态的一种方法。20世纪80年代，Ushakov首次提出UGF概念，随后LISNIANSKI[17]、LEVITIN[18]将UGF应用于多状态系统可靠性分析领域，UGF逐渐成为多状态离散系统可靠性分析的重要工具[19-21]。自LISNIANSKI[17]提出利用UGF计算连续状态系统可靠度指标的界以来，该方法逐渐被扩展到具有连续型随机变量的结构失效概率分析中[22-25]。连续变量结构系统经变量离散化后，状态空间通常较为庞大，对这类系统进行可靠性分析时会产生海量的状态组合，因此该方法未能广泛应用于结构可靠性分析中。

由于现有的可靠度求解方法无法同时满足复杂结构可靠度求解的精度和效率要求，故本文将UGF工具、自适应细分原理以及重要抽样技术有机结合，提出结构可靠度求解的自适应细分-重要抽样法。该方法针对涉及非正态随机变量、小失效概率以及高度非线性功能函数的可靠性问题，在可控的计算成本内保证了可靠度求解的精度。

1 结构可靠性分析的UGF法

UGF的描述对象为离散型随机变量，因此在建立连续型随机变量UGF时，需将随机变量离散化。设连续型随机变量X的累积分布函数及概率密度函数分别为FX(x)和fX(x)，将其在定义域(xmin，xmax)内均匀离散成m个点，由下式计算各离散点xi对应的概率值：

(1)

式中，离散步长δ=(xmax-xmin)/m。

因此，根据离散数据集{(xi，pi)|i=1，2，…，m}定义连续型随机变量X的UGF为

(2)

设n维连续型随机向量X=(X1,X2,…,Xn)，根据式(1)、式(2)获得X各分量的UGF，记为

(3)

式中，xiji、piji分别为随机变量Xi的第ji个状态值和对应的概率值；mi为Xi的离散状态数。

设结构功能函数为G(X)，当G(X)≤0时，结构失效；当G(X)>0时，结构可靠。对各分量的UGF进行复合运算，获得结构系统的UGF，即

UG(z)=ΦG(UX1(z),UX2(z),…,UXn(z))=

(4)

式中，UG(z)为功能函数G(X)的结构UGF；ΦG为复合算子，用于描述各变量UGF间的运算规则。

结构UGF由式(4)可进一步简化为

(5)

依据该UGF的概率项进行条件求和，得到结构的可靠度为

(6)

式中，λ(·)为条件求和算子；I(·)为示性函数，当Gk<0时I(Gk)取0，否则取1。

显然，UGF法通过枚举结构系统各随机变量离散状态，并通过递推运算实现结构的可靠性分析。该方法适用于涉及非正态随机变量以及高度非线性功能函数的可靠性问题。

2 临界域的自适应细分原理

根据式(4)定义随机变量状态组(x1j1,x2j2,…,xnjn)对应的n维空间为焦元Q(j1,j2,…,jn)。通过各变量的复合运算可得到若干个n维空间焦元Qk及对应的概率值pk和功能函数值Gk。设各焦元功能函数极小值和极大值分别为GkL、GkU。如果GkU<0，焦元Qk在失效域；如果GkL>0，焦元Qk在可靠域；如果GkL≤0≤GkU，焦元Qk在临界域。

计算焦元功能函数极值，可采取两种简化方法。一种是分别计算焦元Qk超立方体各顶点的功能函数值，取最小值为GkL，最大值为GkU。该方法需要调用2n次功能函数，当随机变量个数较多时，计算成本偏高。另一种方法是引入区间分析技术高效求解焦元功能函数极值[26]，基本原理如下。

将功能函数G(X)在焦元Qk的中心点Xk处进行一阶泰勒展开：

(7)

式中，Xi为随机向量X的分量；xki为焦元中心点Xk的分量，i=1,2,…,n。

式(7)中的偏导数可用前向差分法求解，即

(8)

式中，Γi为n维向量，该向量的第i个分量为1，其余分量均为0；δi为随机变量xi的离散步长。

将式(8)代入式(7)，则焦元子空间Qk中功能函数的极小值GkL和极大值GkU可通过下式求解：

(9)

(10)

GkL=min(G1，G2)GkU=max(G1，G2)

利用式(9)、式(10)求解焦元功能函数极值时，由于焦元中心点Xk的功能函数值G(Xk)等于Gk，直接由结构UGF相关信息获知，所以计算每个焦元的功能函数极值时仅需再调用n次功能函数。

基于UGF的可靠性分析方法通常使用递推操作实现随机变量的复合运算，并借助同类项合并技术，以获得结构性能的概率分布。然而，离散步长较大时，会导致分析精度不足；离散步长较小时，会增加复合运算的成本。因此，基于自适应细分原理[27]，对临界域进行细分，减少离散区间长度，并通过递归操作对多变量进行非均匀自适应细分(图1)。

图1 临界域的自适应细分

根据随机变量的定义域确定结构初始状态的UGF，即

(11)

式中，UAB(z)为超立方体A≤X≤B的结构UGF；A和B分别为由X的定义域确定的超立方体的下界和上界；K为超立方体A≤X≤B的焦元总数；pk为超立方体A≤X≤B中第k个焦元对应的概率值。

当焦元功能函数极值满足GkL≤0≤GkU时，焦元Qk处于临界域，相应的变量空间会被细分。通过二分法细分超立方体A≤X≤B，得到2n个超立方体，超立方体A′≤X≤B′的结构UGF可表示为

(12)

(13)

(14)

通过递归运算重复细分后可得细分空间对应的结构UGF：

(15)

式中，“⨁”为递归运算符。

根据细分空间对应的结构UGF，可得结构可靠度

(16)

式中，M为细分空间的焦元总数；φ(·)为示性函数，当GkL>0时取1，当GkU<0时取0。

根据式(15)对变量空间进行重复细分后，包含极限状态的超曲面的区域变小，结构状态UGF的项数增加，概率分析的精度提高。但是随着空间维数的增加，计算量呈指数增长，即使递归操作只针对不断缩小的临界域，也会产生较大的计算成本。因此，在上述方法的基础上，考虑将变量空间细分到一定程度，得到失效域概率以及细分后的临界域，再对临界域焦元进行重要抽样，最终得总体失效概率。递归终止准则可定义为变量子空间大小的体积比小于某一阈值，即

(17)

式中，ε为给定的阈值。

3 面向临界域焦元的重要抽样技术

重要抽样法是一种应用广泛的方差缩减方法，它通过构建重要抽样密度函数替代原来的概率密度函数，使样本尽可能多地落入失效域，从而提高计算效率。

设随机向量X的概率密度函数为fX(X)，结构的失效域为F={X|G(X)≤0}，结构失效概率可由以下积分求得：

(18)

设细分后的临界域为Ω′c，构建重要抽样密度函数hX(X)为临界域Ω′c上的均匀分布，即

(19)

式中，X∈∪(Q′c,1,Q′c,2,…,Q′c,r)；Q′c,k为第k个焦元，k=1,2,…,r；Vc为r个临界域焦元体积之和。

(20)

式中，Eh(·)为按hX(X)抽样的样本数学期望；Xi(i=1,2,…,N)为按hX(X)抽取的N个样本；IF(·)为样本失效示性函数，当样本点处的功能函数值小于0时取1，否则取0。

设临界域Ω′c占概率总量为pc(由临界域焦元概率求和得出)，可表示为

(21)

将式(19)、式(21)代入式(20)，进一步整理可得

(22)

4 自适应细分-重要抽样法的基本步骤

首先对临界域进行自适应细分，失效域焦元概率直接累加求取失效域概率pf1，再对细分后的临界域焦元采用重要抽样技术得到临界域失效概率pf2，两者相加得结构失效概率pf。

综上可得，自适应细分-重要抽样法的基本步骤如下：

(1)由随机向量X的定义域确定超立方体的上下界B和A。根据结构功能函数G(X)及焦元的概率值pk，由式(11)得UAB(z)。

(2)由式(9)、式(10)计算焦元Qk的功能函数极小值GkL和极大值GkU。当GkL≤0≤GkU时，由式(12)～式(14)细分生成子空间，对每个子空间执行由式(17)定义的递归终止准则，如果满足条件，则停止循环操作。否则，当前子空间被更新为初始空间，再次执行以上递归操作。

(3)递归运算完成后，由式(16)定义的加权求和算子λ获得失效域概率pf1。

(4)运用式(4)、式(5)对临界域内各变量UGF进行复合运算，得

(5)将临界域焦元中概率大于足够小阈值e(约低于失效概率预估值两个数量等级)的焦元划为热点焦元，其余为边际焦元。对边际焦元做近似处理，将其在UG(z)中的对应项概率系数求和后取半得边际焦元失效概率。

(7)根据步骤(3)、(5)、(6)的计算结果，求和得总体结构失效概率pf。

5 算例分析及讨论

5.1 数学算例

考虑如下3个功能函数：

(23)

(24)

(25)

其中，随机向量X=(X1，X2，X3)。对于功能函数式(23)，X1服从均值为3.8、标准差为0.6的正态分布，X2服从均值为2.7、标准差为0.6的正态分布，X3服从均值为6.4、标准差为0.6的正态分布。对于功能函数式(24)和功能函数式(25)，X1服从区间[2，5.6]上的均匀分布，X2服从均值为2.7、标准差为0.6的正态分布，X3服从均值为6.4、标准差为0.6的对数正态分布。

首先针对功能函数式(23)计算结构失效概率。根据以上提出的自适应细分-重要抽样法，确定初始超立方体上界B=[6.2 5.1 8.8]和下界A=[1.4 0.3 4]。

运用式(11)得到UAB(z)，进一步由式(12)～式(14)细分生成子空间，递归运算直至满足递归终止准则，得到失效域概率pf1=0.0116。运用式(4)、式(5)对临界域内各变量UGF进行复合运算，得

表1 不同算法计算结果对比(针对功能函数式(23))

针对式(24)、式(25)给出的功能函数，运用FORM法、JC法、AS-IS法和MCS法计算结构的失效概率，相应的计算结果如表2、表3所示。

表2 不同算法计算结果对比(针对功能函数式(24))

表3 不同算法计算结果对比(针对功能函数式(25))

由表1～表3可知，对于涉及小失效概率和计算成本高的功能函数的可靠性问题，MCS方法需要大量的样本才能获得精确解，计算效率较低[10]。例如，针对功能函数式(25)，MCS方法计算结构失效概率估计值的变异系数达到0.1904，需要显著提高抽样次数才能缩小估计值的分散性，提高可靠度求解精度。FORM法将设计点附近的结构功能函数近似为线性函数，对于线性或弱非线性结构可靠度求解而言，计算精度尚可，但是，对于算例所示的高度非线性结构，计算精度不符合工程分析精度要求[28]。JC法于涉及单峰低偏度的类正态分布变量以及非线性程度不高的结构可靠性问题有效，但是，对于涉及均匀分布变量以及高度非线性功能函数的结构可靠性问题，存在较大误差或失效[29]。本文提出的结构可靠度求解的自适应细分-重要抽样法，能在计算成本可控的条件下获得满意的计算精度，方法可行。

5.2 工程算例

图2所示为各向同性材料制作的悬臂梁[30]，横向分布载荷q服从均值为10 N/mm、标准差为2 N/mm的对数正态分布，悬臂梁长度l服从区间[3000，3800] mm上的均匀分布，弹性模量E服从均值为73000 Pa、标准差为1000 N/mm2的正态分布，惯性矩I服从均值为1.067×109mm4、标准差为100 000 mm4的正态分布，悬臂梁自由端最大位移Δmax=4 mm。

图2 悬臂梁结构

根据悬臂梁自由端的位移公式，定义功能函数为

(26)

式中，X=(q，l，E，I)。

针对式(26)给出的功能函数，运用FORM法、JC法、AS-IS法以及MCS法计算结构失效概率，相应的计算结果如表4所示。

表4 不同算法计算结果对比(针对功能函数式(26))