王新愿 周金宇 謝里陽 程錦翔

1.江蘇理工學院機械工程學院，常州，213001 2.金陵科技學院機電工程學院，南京，211169 3.東北大學機械工程與自動化學院，沈陽，110819

0 引言

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

針對涉及非正態(tài)隨機變量的可靠性問題，通常對隨機變量當量正態(tài)化，這種轉換會增加功能函數(shù)的非線性程度，降低可靠度求解的精度。鞍點近似(saddlepoint approximation，SA)方法可以直接評估具有非正態(tài)隨機變量的功能函數(shù)的概率分布[3-4]，得到了廣泛應用和發(fā)展。ZHANG等[5]將SA法與線抽樣(line sampling，LS)法相結合，提出鞍點近似線采樣(SA-LS)法，具有準確性高、效率高的優(yōu)點。LU等[6]提出一種改進的基于高階矩的SA法，提高了可靠度求解的精度。目前改進方法僅對單峰低偏度的類正態(tài)分布有效，而對高偏度、多峰分布或均勻分布存在較大誤差或導致失效。針對涉及小失效概率的可靠性問題，國內外學者引入重要抽樣(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方法，可以有效處理涉及小失效概率和高維功能函數(shù)的可靠性問題。針對涉及高度非線性功能函數(shù)的可靠性問題，傳統(tǒng)的基于梯度的方法可能無法得到收斂解。ELEGBEDE[13]應用粒子群算法求解非線性結構失效概率，可得到較精確的結果。ZOU等[14]受粒子群算法的啟發(fā)，提出一種全局協(xié)調搜索算法(NGHS)，具有更高的效率。YANG[15]引入混沌控制方法，提高了非線性結構可靠度求解的精度，但效率大幅降低。KESHTEGAR[16]引入共軛混沌控制方法，與混沌控制方法相比，提高了可靠度求解的效率。然而，對于同時涉及非正態(tài)隨機變量、小失效概率以及高度非線性功能函數(shù)的可靠性問題，現(xiàn)有方法尚不能提供有效的解決方案。

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

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

1 結構可靠性分析的UGF法

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

(1)

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

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

(2)

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

(3)

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

設結構功能函數(shù)為G(X)，當G(X)≤0時，結構失效；當G(X)>0時，結構可靠。對各分量的UGF進行復合運算，獲得結構系統(tǒng)的UGF，即

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

(4)

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

結構UGF由式(4)可進一步簡化為

(5)

依據(jù)該UGF的概率項進行條件求和，得到結構的可靠度為

(6)

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

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

2 臨界域的自適應細分原理

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

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

將功能函數(shù)G(X)在焦元Qk的中心點Xk處進行一階泰勒展開：

(7)

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

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

(8)

式中，Γi為n維向量，該向量的第i個分量為1，其余分量均為0；δi為隨機變量xi的離散步長。

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

(9)

(10)

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

利用式(9)、式(10)求解焦元功能函數(shù)極值時，由于焦元中心點Xk的功能函數(shù)值G(Xk)等于Gk，直接由結構UGF相關信息獲知，所以計算每個焦元的功能函數(shù)極值時僅需再調用n次功能函數(shù)。

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

圖1 臨界域的自適應細分

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

(11)

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

當焦元功能函數(shù)極值滿足GkL≤0≤GkU時，焦元Qk處于臨界域，相應的變量空間會被細分。通過二分法細分超立方體A≤X≤B，得到2n個超立方體，超立方體A′≤X≤B′的結構UGF可表示為

(12)

(13)

(14)

通過遞歸運算重復細分后可得細分空間對應的結構UGF：

(15)

式中，“?”為遞歸運算符。

根據(jù)細分空間對應的結構UGF，可得結構可靠度

(16)

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

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

(17)

式中，ε為給定的閾值。

3 面向臨界域焦元的重要抽樣技術

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

設隨機向量X的概率密度函數(shù)為fX(X)，結構的失效域為F={X|G(X)≤0}，結構失效概率可由以下積分求得：

(18)

設細分后的臨界域為Ω′c，構建重要抽樣密度函數(shù)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)抽樣的樣本數(shù)學期望；Xi(i=1,2,…,N)為按hX(X)抽取的N個樣本；IF(·)為樣本失效示性函數(shù)，當樣本點處的功能函數(shù)值小于0時取1，否則取0。

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

(21)

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

(22)

4 自適應細分-重要抽樣法的基本步驟

首先對臨界域進行自適應細分，失效域焦元概率直接累加求取失效域概率pf1，再對細分后的臨界域焦元采用重要抽樣技術得到臨界域失效概率pf2，兩者相加得結構失效概率pf。

綜上可得，自適應細分-重要抽樣法的基本步驟如下：

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

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

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

(4)運用式(4)、式(5)對臨界域內各變量UGF進行復合運算，得

(5)將臨界域焦元中概率大于足夠小閾值e(約低于失效概率預估值兩個數(shù)量等級)的焦元劃為熱點焦元，其余為邊際焦元。對邊際焦元做近似處理，將其在UG(z)中的對應項概率系數(shù)求和后取半得邊際焦元失效概率。

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

5 算例分析及討論

5.1 數(shù)學算例

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

(23)

(24)

(25)

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

首先針對功能函數(shù)式(23)計算結構失效概率。根據(jù)以上提出的自適應細分-重要抽樣法，確定初始超立方體上界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 不同算法計算結果對比(針對功能函數(shù)式(23))

針對式(24)、式(25)給出的功能函數(shù)，運用FORM法、JC法、AS-IS法和MCS法計算結構的失效概率，相應的計算結果如表2、表3所示。

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

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

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

5.2 工程算例

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

圖2 懸臂梁結構

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

(26)

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

針對式(26)給出的功能函數(shù)，運用FORM法、JC法、AS-IS法以及MCS法計算結構失效概率，相應的計算結果如表4所示。

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

6 結論

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

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