結(jié)構(gòu)可靠度求解的自適應(yīng)細(xì)分-重要抽樣法

王新愿 周金宇 謝里陽(yáng) 程錦翔

1.江蘇理工學(xué)院機(jī)械工程學(xué)院，常州，213001 2.金陵科技學(xué)院機(jī)電工程學(xué)院，南京，211169 3.東北大學(xué)機(jī)械工程與自動(dòng)化學(xué)院，沈陽(yáng)，110819

0 引言

迄今為止，研究者們已提出很多可靠度求解方法。其中，基于設(shè)計(jì)點(diǎn)的方法是最為常見(jiàn)的方法[1-2]，該類(lèi)方法將功能函數(shù)在設(shè)計(jì)點(diǎn)處進(jìn)行泰勒展開(kāi)，保留展開(kāi)式一階項(xiàng)或高階項(xiàng)，進(jìn)而計(jì)算結(jié)構(gòu)可靠度，當(dāng)隨機(jī)變量呈明顯非正態(tài)或功能函數(shù)高度非線性時(shí)，計(jì)算誤差較大。數(shù)值模擬方法中最基本的是蒙特卡羅模擬(Monte Carlo simulation，MCS)方法，該方法具有較高的精度和魯棒性，主要用于驗(yàn)證其他方法的有效性。這種方法在小失效概率的情況下，需要大量樣本保證可靠度求解的精度，計(jì)算效率較低。

針對(duì)涉及非正態(tài)隨機(jī)變量的可靠性問(wèn)題，通常對(duì)隨機(jī)變量當(dāng)量正態(tài)化，這種轉(zhuǎn)換會(huì)增加功能函數(shù)的非線性程度，降低可靠度求解的精度。鞍點(diǎn)近似(saddlepoint approximation，SA)方法可以直接評(píng)估具有非正態(tài)隨機(jī)變量的功能函數(shù)的概率分布[3-4]，得到了廣泛應(yīng)用和發(fā)展。ZHANG等[5]將SA法與線抽樣(line sampling，LS)法相結(jié)合，提出鞍點(diǎn)近似線采樣(SA-LS)法，具有準(zhǔn)確性高、效率高的優(yōu)點(diǎn)。LU等[6]提出一種改進(jìn)的基于高階矩的SA法，提高了可靠度求解的精度。目前改進(jìn)方法僅對(duì)單峰低偏度的類(lèi)正態(tài)分布有效，而對(duì)高偏度、多峰分布或均勻分布存在較大誤差或?qū)е率?。針?duì)涉及小失效概率的可靠性問(wèn)題，國(guó)內(nèi)外學(xué)者引入重要抽樣(importance sampling，IS)和子集模擬(subset simulation，SS)方法[7-8]，取得很多研究成果。KURTZ等[9]提出交叉熵重要抽樣(cross entro-py importance sampling，CE-IS)法，該方法與IS法相比，具有更高的抽樣效率。史朝印等[10]結(jié)合自適應(yīng)Kriging(adaptive Kriging，AK)代理模型，提出改進(jìn)的CE-IS法，顯著減小CE-IS法的計(jì)算量，提高了可靠度求解的效率。BOURINET等[11]結(jié)合支持向量機(jī)(support vector machine，SVM)和SS法，提高了可靠度求解的效率和準(zhǔn)確性。HUANG等[12]將SS法與AK代理模型方法相結(jié)合，提出AK-SS方法，可以有效處理涉及小失效概率和高維功能函數(shù)的可靠性問(wèn)題。針對(duì)涉及高度非線性功能函數(shù)的可靠性問(wèn)題，傳統(tǒng)的基于梯度的方法可能無(wú)法得到收斂解。ELEGBEDE[13]應(yīng)用粒子群算法求解非線性結(jié)構(gòu)失效概率，可得到較精確的結(jié)果。ZOU等[14]受粒子群算法的啟發(fā)，提出一種全局協(xié)調(diào)搜索算法(NGHS)，具有更高的效率。YANG[15]引入混沌控制方法，提高了非線性結(jié)構(gòu)可靠度求解的精度，但效率大幅降低。KESHTEGAR[16]引入共軛混沌控制方法，與混沌控制方法相比，提高了可靠度求解的效率。然而，對(duì)于同時(shí)涉及非正態(tài)隨機(jī)變量、小失效概率以及高度非線性功能函數(shù)的可靠性問(wèn)題，現(xiàn)有方法尚不能提供有效的解決方案。

通用生成函數(shù)(universal generating function，UGF)是建立在廣泛使用概率論原理的生成函數(shù)基礎(chǔ)上，用于枚舉分析系統(tǒng)狀態(tài)的一種方法。20世紀(jì)80年代，Ushakov首次提出UGF概念，隨后LISNIANSKI[17]、LEVITIN[18]將UGF應(yīng)用于多狀態(tài)系統(tǒng)可靠性分析領(lǐng)域，UGF逐漸成為多狀態(tài)離散系統(tǒng)可靠性分析的重要工具[19-21]。自LISNIANSKI[17]提出利用UGF計(jì)算連續(xù)狀態(tài)系統(tǒng)可靠度指標(biāo)的界以來(lái)，該方法逐漸被擴(kuò)展到具有連續(xù)型隨機(jī)變量的結(jié)構(gòu)失效概率分析中[22-25]。連續(xù)變量結(jié)構(gòu)系統(tǒng)經(jīng)變量離散化后，狀態(tài)空間通常較為龐大，對(duì)這類(lèi)系統(tǒng)進(jìn)行可靠性分析時(shí)會(huì)產(chǎn)生海量的狀態(tài)組合，因此該方法未能廣泛應(yīng)用于結(jié)構(gòu)可靠性分析中。

由于現(xiàn)有的可靠度求解方法無(wú)法同時(shí)滿(mǎn)足復(fù)雜結(jié)構(gòu)可靠度求解的精度和效率要求，故本文將UGF工具、自適應(yīng)細(xì)分原理以及重要抽樣技術(shù)有機(jī)結(jié)合，提出結(jié)構(gòu)可靠度求解的自適應(yīng)細(xì)分-重要抽樣法。該方法針對(duì)涉及非正態(tài)隨機(jī)變量、小失效概率以及高度非線性功能函數(shù)的可靠性問(wèn)題，在可控的計(jì)算成本內(nèi)保證了可靠度求解的精度。

1 結(jié)構(gòu)可靠性分析的UGF法

UGF的描述對(duì)象為離散型隨機(jī)變量，因此在建立連續(xù)型隨機(jī)變量UGF時(shí)，需將隨機(jī)變量離散化。設(shè)連續(xù)型隨機(jī)變量X的累積分布函數(shù)及概率密度函數(shù)分別為FX(x)和fX(x)，將其在定義域(xmin，xmax)內(nèi)均勻離散成m個(gè)點(diǎn)，由下式計(jì)算各離散點(diǎn)xi對(duì)應(yīng)的概率值：

(1)

式中，離散步長(zhǎng)δ=(xmax-xmin)/m。

因此，根據(jù)離散數(shù)據(jù)集{(xi，pi)|i=1，2，…，m}定義連續(xù)型隨機(jī)變量X的UGF為

(2)

設(shè)n維連續(xù)型隨機(jī)向量X=(X1,X2,…,Xn)，根據(jù)式(1)、式(2)獲得X各分量的UGF，記為

(3)

式中，xiji、piji分別為隨機(jī)變量Xi的第ji個(gè)狀態(tài)值和對(duì)應(yīng)的概率值；mi為Xi的離散狀態(tài)數(shù)。

設(shè)結(jié)構(gòu)功能函數(shù)為G(X)，當(dāng)G(X)≤0時(shí)，結(jié)構(gòu)失效；當(dāng)G(X)>0時(shí)，結(jié)構(gòu)可靠。對(duì)各分量的UGF進(jìn)行復(fù)合運(yùn)算，獲得結(jié)構(gòu)系統(tǒng)的UGF，即

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

(4)

式中，UG(z)為功能函數(shù)G(X)的結(jié)構(gòu)UGF；ΦG為復(fù)合算子，用于描述各變量UGF間的運(yùn)算規(guī)則。

結(jié)構(gòu)UGF由式(4)可進(jìn)一步簡(jiǎn)化為

(5)

依據(jù)該UGF的概率項(xiàng)進(jìn)行條件求和，得到結(jié)構(gòu)的可靠度為

(6)

式中，λ(·)為條件求和算子；I(·)為示性函數(shù)，當(dāng)Gk<0時(shí)I(Gk)取0，否則取1。

顯然，UGF法通過(guò)枚舉結(jié)構(gòu)系統(tǒng)各隨機(jī)變量離散狀態(tài)，并通過(guò)遞推運(yùn)算實(shí)現(xiàn)結(jié)構(gòu)的可靠性分析。該方法適用于涉及非正態(tài)隨機(jī)變量以及高度非線性功能函數(shù)的可靠性問(wèn)題。

2 臨界域的自適應(yīng)細(xì)分原理

根據(jù)式(4)定義隨機(jī)變量狀態(tài)組(x1j1,x2j2,…,xnjn)對(duì)應(yīng)的n維空間為焦元Q(j1,j2,…,jn)。通過(guò)各變量的復(fù)合運(yùn)算可得到若干個(gè)n維空間焦元Qk及對(duì)應(yīng)的概率值pk和功能函數(shù)值Gk。設(shè)各焦元功能函數(shù)極小值和極大值分別為GkL、GkU。如果GkU<0，焦元Qk在失效域；如果GkL>0，焦元Qk在可靠域；如果GkL≤0≤GkU，焦元Qk在臨界域。

計(jì)算焦元功能函數(shù)極值，可采取兩種簡(jiǎn)化方法。一種是分別計(jì)算焦元Qk超立方體各頂點(diǎn)的功能函數(shù)值，取最小值為GkL，最大值為GkU。該方法需要調(diào)用2n次功能函數(shù)，當(dāng)隨機(jī)變量個(gè)數(shù)較多時(shí)，計(jì)算成本偏高。另一種方法是引入?yún)^(qū)間分析技術(shù)高效求解焦元功能函數(shù)極值[26]，基本原理如下。

將功能函數(shù)G(X)在焦元Qk的中心點(diǎn)Xk處進(jìn)行一階泰勒展開(kāi)：

(7)

式中，Xi為隨機(jī)向量X的分量；xki為焦元中心點(diǎn)Xk的分量，i=1,2,…,n。

式(7)中的偏導(dǎo)數(shù)可用前向差分法求解，即

(8)

式中，Γi為n維向量，該向量的第i個(gè)分量為1，其余分量均為0；δi為隨機(jī)變量xi的離散步長(zhǎng)。

將式(8)代入式(7)，則焦元子空間Qk中功能函數(shù)的極小值GkL和極大值GkU可通過(guò)下式求解：

(9)

(10)

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

利用式(9)、式(10)求解焦元功能函數(shù)極值時(shí)，由于焦元中心點(diǎn)Xk的功能函數(shù)值G(Xk)等于Gk，直接由結(jié)構(gòu)UGF相關(guān)信息獲知，所以計(jì)算每個(gè)焦元的功能函數(shù)極值時(shí)僅需再調(diào)用n次功能函數(shù)。

基于UGF的可靠性分析方法通常使用遞推操作實(shí)現(xiàn)隨機(jī)變量的復(fù)合運(yùn)算，并借助同類(lèi)項(xiàng)合并技術(shù)，以獲得結(jié)構(gòu)性能的概率分布。然而，離散步長(zhǎng)較大時(shí)，會(huì)導(dǎo)致分析精度不足；離散步長(zhǎng)較小時(shí)，會(huì)增加復(fù)合運(yùn)算的成本。因此，基于自適應(yīng)細(xì)分原理[27]，對(duì)臨界域進(jìn)行細(xì)分，減少離散區(qū)間長(zhǎng)度，并通過(guò)遞歸操作對(duì)多變量進(jìn)行非均勻自適應(yīng)細(xì)分(圖1)。

圖1 臨界域的自適應(yīng)細(xì)分

根據(jù)隨機(jī)變量的定義域確定結(jié)構(gòu)初始狀態(tài)的UGF，即

(11)

式中，UAB(z)為超立方體A≤X≤B的結(jié)構(gòu)UGF；A和B分別為由X的定義域確定的超立方體的下界和上界；K為超立方體A≤X≤B的焦元總數(shù)；pk為超立方體A≤X≤B中第k個(gè)焦元對(duì)應(yīng)的概率值。

當(dāng)焦元功能函數(shù)極值滿(mǎn)足GkL≤0≤GkU時(shí)，焦元Qk處于臨界域，相應(yīng)的變量空間會(huì)被細(xì)分。通過(guò)二分法細(xì)分超立方體A≤X≤B，得到2n個(gè)超立方體，超立方體A′≤X≤B′的結(jié)構(gòu)UGF可表示為

(12)

(13)

(14)

通過(guò)遞歸運(yùn)算重復(fù)細(xì)分后可得細(xì)分空間對(duì)應(yīng)的結(jié)構(gòu)UGF：

(15)

式中，“?”為遞歸運(yùn)算符。

根據(jù)細(xì)分空間對(duì)應(yīng)的結(jié)構(gòu)UGF，可得結(jié)構(gòu)可靠度

(16)

式中，M為細(xì)分空間的焦元總數(shù)；φ(·)為示性函數(shù)，當(dāng)GkL>0時(shí)取1，當(dāng)GkU<0時(shí)取0。

根據(jù)式(15)對(duì)變量空間進(jìn)行重復(fù)細(xì)分后，包含極限狀態(tài)的超曲面的區(qū)域變小，結(jié)構(gòu)狀態(tài)UGF的項(xiàng)數(shù)增加，概率分析的精度提高。但是隨著空間維數(shù)的增加，計(jì)算量呈指數(shù)增長(zhǎng)，即使遞歸操作只針對(duì)不斷縮小的臨界域，也會(huì)產(chǎn)生較大的計(jì)算成本。因此，在上述方法的基礎(chǔ)上，考慮將變量空間細(xì)分到一定程度，得到失效域概率以及細(xì)分后的臨界域，再對(duì)臨界域焦元進(jìn)行重要抽樣，最終得總體失效概率。遞歸終止準(zhǔn)則可定義為變量子空間大小的體積比小于某一閾值，即

(17)

式中，ε為給定的閾值。

3 面向臨界域焦元的重要抽樣技術(shù)

重要抽樣法是一種應(yīng)用廣泛的方差縮減方法，它通過(guò)構(gòu)建重要抽樣密度函數(shù)替代原來(lái)的概率密度函數(shù)，使樣本盡可能多地落入失效域，從而提高計(jì)算效率。

設(shè)隨機(jī)向量X的概率密度函數(shù)為fX(X)，結(jié)構(gòu)的失效域?yàn)镕={X|G(X)≤0}，結(jié)構(gòu)失效概率可由以下積分求得：

(18)

設(shè)細(xì)分后的臨界域?yàn)棣浮鋍，構(gòu)建重要抽樣密度函數(shù)hX(X)為臨界域Ω′c上的均勻分布，即

(19)

式中，X∈∪(Q′c,1,Q′c,2,…,Q′c,r)；Q′c,k為第k個(gè)焦元，k=1,2,…,r；Vc為r個(gè)臨界域焦元體積之和。

(20)

式中，Eh(·)為按hX(X)抽樣的樣本數(shù)學(xué)期望；Xi(i=1,2,…,N)為按hX(X)抽取的N個(gè)樣本；IF(·)為樣本失效示性函數(shù)，當(dāng)樣本點(diǎn)處的功能函數(shù)值小于0時(shí)取1，否則取0。

設(shè)臨界域Ω′c占概率總量為pc(由臨界域焦元概率求和得出)，可表示為

(21)

將式(19)、式(21)代入式(20)，進(jìn)一步整理可得

(22)

4 自適應(yīng)細(xì)分-重要抽樣法的基本步驟

首先對(duì)臨界域進(jìn)行自適應(yīng)細(xì)分，失效域焦元概率直接累加求取失效域概率pf1，再對(duì)細(xì)分后的臨界域焦元采用重要抽樣技術(shù)得到臨界域失效概率pf2，兩者相加得結(jié)構(gòu)失效概率pf。

綜上可得，自適應(yīng)細(xì)分-重要抽樣法的基本步驟如下：

(1)由隨機(jī)向量X的定義域確定超立方體的上下界B和A。根據(jù)結(jié)構(gòu)功能函數(shù)G(X)及焦元的概率值pk，由式(11)得UAB(z)。

(2)由式(9)、式(10)計(jì)算焦元Qk的功能函數(shù)極小值GkL和極大值GkU。當(dāng)GkL≤0≤GkU時(shí)，由式(12)～式(14)細(xì)分生成子空間，對(duì)每個(gè)子空間執(zhí)行由式(17)定義的遞歸終止準(zhǔn)則，如果滿(mǎn)足條件，則停止循環(huán)操作。否則，當(dāng)前子空間被更新為初始空間，再次執(zhí)行以上遞歸操作。

(3)遞歸運(yùn)算完成后，由式(16)定義的加權(quán)求和算子λ獲得失效域概率pf1。

(4)運(yùn)用式(4)、式(5)對(duì)臨界域內(nèi)各變量UGF進(jìn)行復(fù)合運(yùn)算，得

(5)將臨界域焦元中概率大于足夠小閾值e(約低于失效概率預(yù)估值兩個(gè)數(shù)量等級(jí))的焦元?jiǎng)潪闊狳c(diǎn)焦元，其余為邊際焦元。對(duì)邊際焦元做近似處理，將其在UG(z)中的對(duì)應(yīng)項(xiàng)概率系數(shù)求和后取半得邊際焦元失效概率。

(7)根據(jù)步驟(3)、(5)、(6)的計(jì)算結(jié)果，求和得總體結(jié)構(gòu)失效概率pf。

5 算例分析及討論

5.1 數(shù)學(xué)算例

考慮如下3個(gè)功能函數(shù)：

(23)

(24)

(25)

其中，隨機(jī)向量X=(X1，X2，X3)。對(duì)于功能函數(shù)式(23)，X1服從均值為3.8、標(biāo)準(zhǔn)差為0.6的正態(tài)分布，X2服從均值為2.7、標(biāo)準(zhǔn)差為0.6的正態(tài)分布，X3服從均值為6.4、標(biāo)準(zhǔn)差為0.6的正態(tài)分布。對(duì)于功能函數(shù)式(24)和功能函數(shù)式(25)，X1服從區(qū)間[2，5.6]上的均勻分布，X2服從均值為2.7、標(biāo)準(zhǔn)差為0.6的正態(tài)分布，X3服從均值為6.4、標(biāo)準(zhǔn)差為0.6的對(duì)數(shù)正態(tài)分布。

首先針對(duì)功能函數(shù)式(23)計(jì)算結(jié)構(gòu)失效概率。根據(jù)以上提出的自適應(yīng)細(xì)分-重要抽樣法，確定初始超立方體上界B=[6.2 5.1 8.8]和下界A=[1.4 0.3 4]。

運(yùn)用式(11)得到UAB(z)，進(jìn)一步由式(12)～式(14)細(xì)分生成子空間，遞歸運(yùn)算直至滿(mǎn)足遞歸終止準(zhǔn)則，得到失效域概率pf1=0.0116。運(yùn)用式(4)、式(5)對(duì)臨界域內(nèi)各變量UGF進(jìn)行復(fù)合運(yùn)算，得

表1 不同算法計(jì)算結(jié)果對(duì)比(針對(duì)功能函數(shù)式(23))

針對(duì)式(24)、式(25)給出的功能函數(shù)，運(yùn)用FORM法、JC法、AS-IS法和MCS法計(jì)算結(jié)構(gòu)的失效概率，相應(yīng)的計(jì)算結(jié)果如表2、表3所示。

表2 不同算法計(jì)算結(jié)果對(duì)比(針對(duì)功能函數(shù)式(24))

表3 不同算法計(jì)算結(jié)果對(duì)比(針對(duì)功能函數(shù)式(25))

由表1～表3可知，對(duì)于涉及小失效概率和計(jì)算成本高的功能函數(shù)的可靠性問(wèn)題，MCS方法需要大量的樣本才能獲得精確解，計(jì)算效率較低[10]。例如，針對(duì)功能函數(shù)式(25)，MCS方法計(jì)算結(jié)構(gòu)失效概率估計(jì)值的變異系數(shù)達(dá)到0.1904，需要顯著提高抽樣次數(shù)才能縮小估計(jì)值的分散性，提高可靠度求解精度。FORM法將設(shè)計(jì)點(diǎn)附近的結(jié)構(gòu)功能函數(shù)近似為線性函數(shù)，對(duì)于線性或弱非線性結(jié)構(gòu)可靠度求解而言，計(jì)算精度尚可，但是，對(duì)于算例所示的高度非線性結(jié)構(gòu)，計(jì)算精度不符合工程分析精度要求[28]。JC法于涉及單峰低偏度的類(lèi)正態(tài)分布變量以及非線性程度不高的結(jié)構(gòu)可靠性問(wèn)題有效，但是，對(duì)于涉及均勻分布變量以及高度非線性功能函數(shù)的結(jié)構(gòu)可靠性問(wèn)題，存在較大誤差或失效[29]。本文提出的結(jié)構(gòu)可靠度求解的自適應(yīng)細(xì)分-重要抽樣法，能在計(jì)算成本可控的條件下獲得滿(mǎn)意的計(jì)算精度，方法可行。

5.2 工程算例

圖2所示為各向同性材料制作的懸臂梁[30]，橫向分布載荷q服從均值為10 N/mm、標(biāo)準(zhǔn)差為2 N/mm的對(duì)數(shù)正態(tài)分布，懸臂梁長(zhǎng)度l服從區(qū)間[3000，3800] mm上的均勻分布，彈性模量E服從均值為73000 Pa、標(biāo)準(zhǔn)差為1000 N/mm2的正態(tài)分布，慣性矩I服從均值為1.067×109mm4、標(biāo)準(zhǔn)差為100 000 mm4的正態(tài)分布，懸臂梁自由端最大位移Δmax=4 mm。

圖2 懸臂梁結(jié)構(gòu)

根據(jù)懸臂梁自由端的位移公式，定義功能函數(shù)為

(26)

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

針對(duì)式(26)給出的功能函數(shù)，運(yùn)用FORM法、JC法、AS-IS法以及MCS法計(jì)算結(jié)構(gòu)失效概率，相應(yīng)的計(jì)算結(jié)果如表4所示。

表4 不同算法計(jì)算結(jié)果對(duì)比(針對(duì)功能函數(shù)式(26))

6 結(jié)論

(1)針對(duì)涉及非正態(tài)隨機(jī)變量、小失效概率以及高度非線性功能函數(shù)的可靠性問(wèn)題，傳統(tǒng)方法存在求解精度低或無(wú)法求解的問(wèn)題，而本文方法可以獲得較高的求解精度，且計(jì)算效率遠(yuǎn)高于MCS方法。

(2)實(shí)際工程問(wèn)題中的功能函數(shù)通常為隱式函數(shù)，每次調(diào)用都帶來(lái)極大的計(jì)算量。后續(xù)工作將在本文方法的基礎(chǔ)上，進(jìn)一步研究涉及隱式功能函數(shù)的可靠性問(wèn)題的求解方法。