基于子集模擬法的橋式起重機(jī)主梁失效概率分析

李耀宗，范小寧

（太原科技大學(xué) 機(jī)械工程學(xué)院，太原 030024）

0 引言

橋式起重機(jī)，在起重運(yùn)輸行業(yè)中占有重要的地位，其主要由橋架、小車和電氣部分三大部分構(gòu)成，其中作為橋架主要組成部分的主梁，是整機(jī)的主要承載構(gòu)件，一旦發(fā)生結(jié)構(gòu)破壞，輕則造成重大財(cái)產(chǎn)損失，重則引起人員傷亡，因此必須保證其具有非常高的可靠性，即很小的失效概率。目前起重機(jī)金屬結(jié)構(gòu)設(shè)計(jì)計(jì)算中均采用許用應(yīng)力法[1]，該方法通過設(shè)置某一固定的安全系數(shù),在結(jié)構(gòu)構(gòu)件服役期間，截面上任何一點(diǎn)的應(yīng)力值都不應(yīng)大于給定的許用應(yīng)力值，雖然許用應(yīng)力法使用簡單，但其無法充分考慮結(jié)構(gòu)的載荷，材料性質(zhì)，結(jié)構(gòu)尺寸等的不確定性對其可靠性的影響，不能定量的描述可靠性的大小。

針對金屬結(jié)構(gòu)的概率可靠性計(jì)算方法，如直接蒙特卡洛模擬法[2,3]，盡管具有適應(yīng)性強(qiáng)，對變量維數(shù)不敏感，且不需要事先知道失效域的任何信息等優(yōu)點(diǎn)，但由于抽樣數(shù)目過于龐大，導(dǎo)致其在模擬小失效概率時(shí)效率非常低下。其他縮減方差的數(shù)字模擬法，如重要抽樣法[4,5]，分層抽樣法等，盡管比直接蒙特卡洛抽樣法在抽樣次數(shù)上有所降低，但在計(jì)算大型復(fù)雜結(jié)構(gòu)的小失效概率時(shí)其計(jì)算效率仍難以滿足實(shí)際要求。

針對以上問題，本文介紹了一種能極大提高多維小失效概率問題計(jì)算效率的方法-子集模擬方法[6]，它的主要思想是通過選擇合適的中間失效事件，進(jìn)而可以將概率空間劃分成一系列的子集，小的失效概率則采用一系列較大的中間失效事件條件概率的乘積來表示，該方法最早由Au和Beck于2001年提出，并在多個(gè)領(lǐng)域取得廣泛而成功的應(yīng)用。文獻(xiàn)[7,8]將子集模擬法用于邊坡可靠性分析中，且文獻(xiàn)[8]解決了含相關(guān)非正態(tài)變量的邊坡可靠性問題。文獻(xiàn)[9,10]將子集模擬法用于結(jié)構(gòu)的動力可靠度計(jì)算，并取得了可喜的研究成果, 可見其在求解小概率失效問題時(shí)的有效性和可行性, 因而成為當(dāng)前進(jìn)行結(jié)構(gòu)可靠度研究的熱點(diǎn)問題。鑒于此, 本文將子集模擬法引入到起重機(jī)金屬結(jié)構(gòu)的可靠度計(jì)算中，通過對簡單數(shù)值算例和工程算例可靠度進(jìn)行子集模擬，并與直接蒙特卡洛法計(jì)算結(jié)果進(jìn)行對比，表明所提方法不受極限狀態(tài)方程非線性程度的影響，可有效減少計(jì)算成本，具有較高的計(jì)算效率和精度。

1 子集模擬法

1.1 子集模擬法的基本原理

子集模擬法[6,11,12]的主要思想是通過選擇合適的中間失效事件，將小概率問題轉(zhuǎn)化為一系列較大條件概率的乘積。若結(jié)構(gòu)的功能函數(shù)為Z(x)，失效域定義為F={x|Z(x)＞lim}，通過引入一系列的臨界值t1＜t2＜...＜tm=lim，可以形成一系列具有包含關(guān)系的中間失效事件Fk，F(xiàn)k={x|Z(x)＞tk},(k=1,2,...,m)，即F1? F2?...?Fm=F ，且有。根據(jù)條件概率中乘法定理的定義，失效概率可以進(jìn)一步表述為：

令P1=P(F1),Pj=P(Fj|Fj?1)(,j=2,3,...,m)

若m=4，Pj為0.1量級，則Pf可以達(dá)到10-4量級，在其他情況都相同的情況下，不管是從效率還是精度方面來說，數(shù)字模擬法計(jì)算0.1量級的失效概率，都比計(jì)算10-4量級更有優(yōu)勢。對于10-4量級的失效概率，若采用直接蒙特卡羅法，則需要100/Pf[13,14]即1000000次抽樣才能達(dá)到需要的精度，若采用子集模擬法，則只需大概4×100/Pj次即4000次抽樣就可以達(dá)到需要的精度，效率提高了99.6%。通過以上分析可知，子集模擬法可以把小的失效概率表示成一連串容易求解的相對較大的失效概率的乘積，采用數(shù)字模擬法則可以高效的求解相對較大的失效概率，從而顯著提高小失效概率的計(jì)算效率。

1.2 條件樣本點(diǎn)的 MCMC 模擬

式(2)中P1的估計(jì)值 P?1可以通過直接蒙特卡羅方法獲得，而 Pj, (j = 2,3,...,m)的估計(jì)值P?j則可以通過抽取條件樣本點(diǎn)來進(jìn)行計(jì)算。

子集模擬法的關(guān)鍵在于構(gòu)造合適的馬爾科夫鏈[6,15]，并通過馬爾科夫鏈獲得一系列的樣本，隨著馬爾科夫鏈鏈長的不斷擴(kuò)展，鏈上樣本相對應(yīng)的工程系統(tǒng)響應(yīng)值的概率分布會越來越接近響應(yīng)的真實(shí)分布。馬爾可夫鏈蒙特卡羅法中出現(xiàn)較早且目前應(yīng)用較為廣泛的方法是M-H算法[11,15]，該算法要求隨機(jī)向量的各個(gè)分量相互獨(dú)立[16]，此法通過在當(dāng)前樣本點(diǎn)附近的隨機(jī)游動來得到下一點(diǎn)，進(jìn)而可以構(gòu)建出一條馬爾科夫鏈，簡稱馬氏鏈。設(shè)馬氏鏈的當(dāng)前狀態(tài)點(diǎn)為x，x∈F，F(xiàn)為失效域，M-H算法的目的是得到下一個(gè)狀態(tài)點(diǎn)y∈F，其具體操作過程如下：

1）選取具有對稱性的建議分布函數(shù)p(y|x)[6,11]，因其對稱性，有 p (y|x)= p(x|y)。從理論上講, 建議分布函數(shù)p(y|x)可以隨意選取，但實(shí)際應(yīng)用中，算法效率會受到建議分布函數(shù)的影響。目前普遍認(rèn)為建議分布函數(shù)的形式越接近目標(biāo)函數(shù)分布[6,15]，則模擬的效果越好。實(shí)際應(yīng)用中多取一維均勻分布或正態(tài)分布作為“建議分布函數(shù)”，本文選取一維均勻分布為建議分布函數(shù)，建議分布函數(shù)的標(biāo)準(zhǔn)差取為對應(yīng)隨機(jī)變量的標(biāo)準(zhǔn)差。

2）以當(dāng)前樣本點(diǎn)作為建議分布函數(shù)的中心，通過建議分布函數(shù)隨機(jī)產(chǎn)生一點(diǎn)ζ，計(jì)算接受概率r=min(1 ,f(ζ)/f(x))，其中f為功能函數(shù)。以概率r接受ζ，即=ζ，以概率1-r拒絕ζ，即=x。具體可以此操作：從均勻分布[0,1]隨機(jī)產(chǎn)生一個(gè)數(shù)字u，若u＜r，則=ζ，若u＞r，則=x。

第二步：接受或拒絕階段

第三步：產(chǎn)生新的狀態(tài)點(diǎn)

繼續(xù)以新的y作為新的當(dāng)前狀態(tài)點(diǎn)，重復(fù)第一步和第二步，即可產(chǎn)生下一個(gè)新的狀態(tài)點(diǎn)。重復(fù)一、二和三步，即可生成一條完整的馬氏鏈。

以上過程僅為一維狀態(tài)下產(chǎn)生馬氏鏈的過程，在本文的研究背景下，各個(gè)設(shè)計(jì)變量是獨(dú)立同分布的，所以在具有多個(gè)設(shè)計(jì)變量的情況下，只需將當(dāng)前樣本點(diǎn)的各個(gè)分量作為馬氏鏈的起始狀態(tài)點(diǎn)，分別對每個(gè)分量重復(fù)步驟一、二和三，即可產(chǎn)生一條多維狀態(tài)下的馬氏鏈。

1.3 中間失效事件的選擇

中間失效事件[6,11]的選擇在子集模擬過程中占據(jù)重要地位。一種情況下，如果中間失效事件非常多(m值很大)，即 tj, (j = 1,2,...,m)減小變慢，則相應(yīng)的條件失效概率可以很大，用較少的條件樣本點(diǎn)即可以模擬對應(yīng)的條件失效概率，但總的樣本點(diǎn)數(shù)將會增多；在另一種情況下，中間事件比較少，則相應(yīng)的條件失效概率會很小，此時(shí)需要較多的條件樣本點(diǎn)才能精確的模擬很小的條件失效概率，同樣會導(dǎo)致總的抽樣點(diǎn)數(shù)的增多。故目前中間失效事件的選擇方法為：在模擬條件失效概率的樣本點(diǎn)數(shù)Nj和中間失效事件的數(shù)量m上折衷[6]，即通過提前選定條件概率p0，根據(jù)實(shí)際需要自適應(yīng)地進(jìn)行分層，本文亦采用該法進(jìn)行中間失效事件的選擇。

1.4 子集模擬法的具體步驟

1）用直接蒙特卡羅模擬法產(chǎn)生N（本文中N=1000）個(gè)服從聯(lián)合概率密度函數(shù)為fX(x)的相互獨(dú)立的樣本{:k= 1,2,...,N}。

3）從落在 Fj-1(j =2,3,...,m)域內(nèi)的p0N個(gè)樣本點(diǎn)出發(fā)，以每個(gè)樣本點(diǎn)作為馬氏鏈的種子，共產(chǎn)生p0N條馬氏鏈，每條馬氏鏈可產(chǎn)生(1-p0)/p0個(gè)新的樣本點(diǎn)，p0N條馬爾科夫鏈共產(chǎn)生(1-p0)/N個(gè)條件樣本點(diǎn)，這(1-p0)/N個(gè)新產(chǎn)生的樣本點(diǎn)與作為種子的p0N個(gè)樣本點(diǎn)總數(shù)共為N，且它們均為服從密度函數(shù) q (x| Fj-1)的條件樣本點(diǎn)，即{,k = 1,2,...,N}。

4）與步驟（2）相似，用功能函數(shù)Z(x)計(jì)算這N個(gè)條件樣本相應(yīng)的系統(tǒng)響應(yīng)值，{Z(),k=1,2,...,N;j=2,3,...,m}，同樣對響應(yīng)值按照由小到大的順序排序，以第(1-p0)/N個(gè)值作為中間失效事 件 Fj={x|Z(X )＞tj}的 臨 界 值tj， 同 時(shí) 得 到P（Fj|Fj-1)=p0。

5）重復(fù)步驟(3)和(4)，直到第(1-p0)/N個(gè)響應(yīng)值的值大于lim，令tm=1im，自動分層至此結(jié)束。計(jì)算落入失效域F中的服從密度函數(shù) q (x| Fm-1)的條件樣本點(diǎn)的個(gè)數(shù)Nf，可得 P(Fm|Fm-1)=Nf/N。

2 算例分析

為證明所提方法的可行性和有效性，以下給出數(shù)值算例和工程算例，并將子集模擬法的計(jì)算結(jié)果與直接蒙特卡羅法進(jìn)行對比。為研究方便假定算例中各隨機(jī)變量之間相互獨(dú)立，且均為正態(tài)分布（非正態(tài)分布的隨機(jī)變量可通過Rosenblatt變換轉(zhuǎn)換成正態(tài)分布[17]）。

2.1 數(shù)值算例

數(shù)值算例1

g(x1,x2,x3)=x31+2x2-5x3，界限值為5，其中x1～N(5,0.1)，x2～N(36,0.2)，x3～N(39,0.1)。

數(shù)值算例2

g(x1,x2,x3,x4)=x31+2.x22-2.x3.X4，界限值為0，其中x1～N(2,0,01)，x2～N(5,0.02)，x3～N(5,0.01)，x4～N(6,0.03)。

數(shù)值算例3

g(x1,x2,x3,x4,x5)=x51-x22/(x3.x4)-2.x25，界限值為-2，其中x1～N(3,0.003)，x2～N(15,0.02)，x3～N(2,0.01)，x4～N(10,0.5)，x5～N(11,0.01)。

三個(gè)數(shù)值算例的直接蒙特卡羅法與子集模擬法的計(jì)算結(jié)果如表1所示。

表1 數(shù)值算例計(jì)算結(jié)果

三個(gè)數(shù)值算例分析結(jié)果顯示, 盡管直接蒙特卡羅法具有抽樣量大，計(jì)算效率低下的缺點(diǎn)，但是由于其計(jì)算精度高，適應(yīng)性強(qiáng)，通常將其計(jì)算結(jié)果作為精確值與其他計(jì)算方法進(jìn)行對比。從表1的計(jì)算結(jié)果可見, 三個(gè)數(shù)值案例的子集模擬法與直接蒙特卡羅法的計(jì)算誤差分別為8.3%，3.7%，4.3%，而子集模擬法的抽樣次數(shù)分別僅為直接蒙特卡洛法的3%，0.05%，0.000009%，可以看出，在計(jì)算精度近似相同的情況下，與直接蒙特卡洛法相比，子集模擬法所需的抽樣次數(shù)顯著降低，從而提高了計(jì)算效率。

2.2 工程算例

已知某偏軌雙梁橋式起重機(jī)，其性能參數(shù)為：額定起重量mQ=32t，小車質(zhì)量mx=1t，吊具質(zhì)量m0=1t，跨度25.5m，工作級別為A6。雙梁橋式起重機(jī)的兩根主梁對稱布置，載荷通過小車車輪作用于主梁，兩根主梁具有相同的的受力特征，以一根主梁進(jìn)行分析即可。圖1是主梁受力圖，圖2為主梁的截面圖。滿載小車位于主梁跨中是主梁最危險(xiǎn)的工況，因此以滿載小車位于主梁跨中時(shí)跨中的靜剛度和跨中截面的第3個(gè)危險(xiǎn)點(diǎn)的應(yīng)力[18]建立功能函數(shù)如下：

靜剛度功能函數(shù)：

許用撓度為L/800。

3點(diǎn)應(yīng)力功能函數(shù)：

許用應(yīng)力為175MPa。

式中：

以上各式中：x2,y2：主梁形心坐標(biāo)；Ix,Iy：主梁慣性矩；I2：端梁慣性矩；L：主梁跨度；K：小車軌距；B0：大車軸距；B：端梁全長；b：小車輪距；d1：大車運(yùn)行機(jī)構(gòu)與大車軌道距離；d2：司機(jī)室與大車軌道距離；mg：小車軌道質(zhì)量；ml：欄桿質(zhì)量；mj：大車運(yùn)行機(jī)構(gòu)質(zhì)量；ms：司機(jī)室質(zhì)量。

圖1 主梁受力圖

圖2 主梁截面圖

表2 基本隨機(jī)變量分布參數(shù)表

分別用直接蒙特卡羅法和子集模擬法計(jì)算本例的失效概率，結(jié)果如表3所示。

表3 工程算例計(jì)算結(jié)果

該工程算例結(jié)果分析顯示，在2.2的工程算例中,以靜剛度的失效概率為例，當(dāng)每層的抽樣次數(shù)N=1000時(shí)，子集模擬法的抽樣次數(shù)為直接蒙特卡羅法的0.05%，而計(jì)算誤差為22.2%，此時(shí)的計(jì)算誤差還是比較大的；當(dāng)N=2000，N=3000時(shí)，子集模擬法的抽樣次數(shù)分別為直接蒙特卡羅法的0.1%和0.15%，而計(jì)算誤差卻逐步降低為8.9%和3.7%；由此可見，隨著子集模擬法抽樣次數(shù)的增加，子集模擬法的計(jì)算結(jié)果越來越接近直接蒙特卡羅法的計(jì)算結(jié)果，然而子集模擬法的抽樣次數(shù)仍遠(yuǎn)遠(yuǎn)低于直接蒙特卡羅法的抽樣次數(shù)，由此可見，子集模擬法在計(jì)算小失效概率時(shí)的計(jì)算效率遠(yuǎn)遠(yuǎn)高于直接蒙特卡羅法；點(diǎn)3的計(jì)算結(jié)果具有類似的規(guī)律。數(shù)值算例和工程算例的結(jié)果表明，與直接蒙特卡羅法相比，子集模擬法幾乎具有同樣的計(jì)算精度，且抽樣次數(shù)大為減少，從而提高了計(jì)算效率。

3 結(jié)論

1）本文主要介紹了基于M-H抽樣方法的子集模擬法，并將子集模擬法應(yīng)用于橋式起重機(jī)主梁的失效概率計(jì)算，通過與直接蒙特卡羅法進(jìn)行對比，驗(yàn)證了所提方法在計(jì)算高維小失效概率時(shí)的可行性與計(jì)算精確度。

2）通過工程算例對比可知，當(dāng)每層的抽樣次數(shù)N=1000時(shí)，在失效概率降低1個(gè)數(shù)量級的情況下，蒙特卡羅法的抽樣次數(shù)擴(kuò)大為原來的10倍，而子集模擬法僅需增加1個(gè)分層的樣本數(shù)即N=1000即可，此時(shí)抽樣次數(shù)僅為原來的1.2倍，抽樣次數(shù)的減少是顯而易見的。在失效概率為10-6量級，每層的抽樣次數(shù)N=3000時(shí)，子集模擬法的抽樣次數(shù)僅為蒙特卡羅法的0.018%，效率提高了99.98%。根據(jù)大數(shù)定律，直接蒙特卡羅法和子集模擬法的計(jì)算精度都會隨著抽樣點(diǎn)數(shù)的增加而提高，兩種方法的區(qū)別是后者與前者相比具有更快的收斂速度，換言之就是：在精度近似一樣的情況下，子集模擬法所需抽樣點(diǎn)數(shù)遠(yuǎn)小于直接蒙特卡羅法。

在本文的研究背景下，隨機(jī)變量均為正態(tài)分布，若隨機(jī)變量為非正態(tài)分布，可以采用等概率變換原則將非正態(tài)變量等效為正態(tài)變量即可求解。

[1] 國家標(biāo)準(zhǔn) GB/T3811-2008.起重機(jī)設(shè)計(jì)規(guī)范[M].北京:中國標(biāo)準(zhǔn)出版社,2008.

[2] 鄔錢涌,程文明,崔寶祥.岸邊集裝箱起重機(jī)結(jié)構(gòu)輕量化設(shè)計(jì)及可靠性分析[J].機(jī)械設(shè)計(jì),2012,29(6):40-43,47.

[3] 葉黔元,陳錦輝,朱勝波.堆垛起重機(jī)鋼結(jié)構(gòu)體系可靠性分析的蒙特卡羅模擬[J].上海理工大學(xué)學(xué)報(bào),2006,28(1):83-86,90.

[4] Au S K, Beck J L.A New Adaptive Important Sampling Scheme for Reliability Calculations[J].Structural Safety,1999,21(2):139-163.

[5] Au S K. Probabilistic Failure Analysis by Importance Sampling Markov Chain Simulation[J].Journal of Engineering Mechanics, 2004,130(3):303-311.

[6] Au S K,Beck J L.Estimation of small failure probabilities in high dimensions by simulation[J].Probabilistic Engineering Mechanics, 2001,16(4):263-277.

[7] 曹子君.子集模擬在邊坡可靠性分析中的應(yīng)用[D].成都:西南交通大學(xué),2009.

[8] 張曼,唐小松,李典慶.含相關(guān)非正態(tài)變量邊坡可靠度分析的子集模擬方法[J].武漢:武漢大學(xué)學(xué)報(bào)（工學(xué)版）,2012,45(1):41-45,53.

[9] 徐瑞,張家興,蘇成.非平穩(wěn)隨機(jī)激勵下結(jié)構(gòu)動力可靠度時(shí)域顯式子集模擬法[J].工程力學(xué),2013,30(7):28-33,39.

[10] 劉佩,姚謙峰.基于子集模擬法的非線性結(jié)構(gòu)動力可靠度計(jì)算[J].工程力學(xué),2010,27(12):72-76.

[11] Au S K, Beck J L.Subset simulation and its application to seismic risk based on dynamic analysis. Journal of engineering mechanics. 2003,129(8):901-917.

[12] Au S K,Ching J,Beck J L.Application of subset simulation methods to reliability benchmark problems[J].Structural Safety, 2007,29(3):183-193.

[13] Rubinstein RY.Simulation and the Monte-Carlo method.New York: Wiley,1981.

[14] Fishman GS.Monte Carlo:concepts, algorithm,and applications. New York:Springer,1996.

[15] 茆詩松,王靜龍,濮曉龍.高等數(shù)理統(tǒng)計(jì)[M].北京:高等教育出版社,1998:450-454.

[16] Haario H, Saksman E, Tamminen J.Componentwise adaptation for high dimensional MCMC. Computational Statistics,2005,20(2):265-273.

[17] Rosenblatt M.Remarks on a multivariate transformation. Annals of Mathematical Statistics 1952,23:470-472.

[18] 徐格寧.金屬裝備結(jié)構(gòu)設(shè)計(jì)（第二版）[M].機(jī)械工業(yè)出版社，2009.09.