馬永亮，曲先強(qiáng)，崔洪斌，石德新

(哈爾濱工程大學(xué) 多體船技術(shù)國(guó)防重點(diǎn)學(xué)科實(shí)驗(yàn)室，黑龍江 哈爾濱 150001)

結(jié)構(gòu)系統(tǒng)可靠性的最基本形式可以分為串聯(lián)系統(tǒng)和并聯(lián)系統(tǒng).根據(jù)一階可靠性方法(FORM)，Hohenbichler提出串聯(lián)系統(tǒng)和并聯(lián)系統(tǒng)可靠性都可以表示為高維正態(tài)積分的形式［1］.通過對(duì)高維正態(tài)積分的求解可以獲得結(jié)構(gòu)系統(tǒng)的失效概率.然而，由于高維正態(tài)概率函數(shù)的不可積性，該問題的求解必須借助于各種數(shù)值方法［2］.直接進(jìn)行數(shù)值積分是一種最基本的方法，Milton采用改進(jìn) Simpson方法［3］，Drezner采用一種特殊的 Gauss-Hermite 法進(jìn)行高維正態(tài)積分的求解［4］.Milton方法對(duì)于6維以上的積分，Drezner方法對(duì)于10維以上的積分很難實(shí)現(xiàn)［3，5］.因此，高維正態(tài)積分近似方法的研究愈加重要.Hohenbichler和Rackwitz提出了FOMN方法，F(xiàn)OMN方法受積分維數(shù)和相關(guān)系數(shù)的影響很大［1］;Tang和Melchers提出了I-FOMN方法和G-FOMN方法，對(duì)FOMN方法的計(jì)算精度進(jìn)行了改進(jìn)，其中GFOMN方法的計(jì)算精度較高［6］;Pandey提出了PCM方法［7］;Mori和Kato指出G-FOMN方法和PCM方法在進(jìn)行串聯(lián)系統(tǒng)可靠性計(jì)算時(shí)誤差可高達(dá)50%［8］;Yuan和Pandey后來提出了PCM方法的改進(jìn)方法IPCM方法［9］.這些近似方法雖然計(jì)算方便H，但計(jì)算誤差很大有時(shí)甚至是錯(cuò)誤的，所以仍然需要尋求計(jì)算效率高、計(jì)算精度好的方法.針對(duì)這一需要，本文采用基于 Halton 序列［10］和數(shù)論網(wǎng)格［11］的準(zhǔn) Monte Carlo方法求解高維正態(tài)積分，通過比較認(rèn)為基于數(shù)論網(wǎng)格的準(zhǔn)Monte Carlo方法計(jì)算精度高.從而提出一種串聯(lián)結(jié)構(gòu)系統(tǒng)可靠性計(jì)算方法，并將該方法和其他方法進(jìn)行了比較.

1 高維正態(tài)積分求解的準(zhǔn)Monte Carlo方法

1.1 高維正態(tài)積分表達(dá)式及積分變換

式中:

用于串聯(lián)結(jié)構(gòu)系統(tǒng)可靠性計(jì)算的高維正態(tài)積分表達(dá)式為

式中:m為失效模式的個(gè)數(shù)，也是積分的維數(shù);β=(β1，β2，…βm)為各個(gè)失效模式的可靠性指標(biāo);X=［x1x2… xm］T為積分變量，R(m×m)為失效模式之間的相關(guān)系數(shù)矩陣.

高維正態(tài)積分Φm(β，R)需要采用積分變換和數(shù)值方法進(jìn)行求解.由于Φm(β，R)的積分域?yàn)榘霟o限域，很難進(jìn)行積分求解，所以進(jìn)行合適的積分變換是必要的.通過以下3步變換可以將高維正態(tài)積分轉(zhuǎn)化為單位超立方域內(nèi)的高維積分［12］.

1)對(duì) R(m ×m)進(jìn)行 Cholesky分解［13］，得到下三角矩陣L，并令X=LY，式(1)可以表示為

式中:φ(y1)為標(biāo)準(zhǔn)正態(tài)分布概率密度函數(shù)，為積分上限 β＇i

2)使用變換，yi=Φ－1(Zi)，式(2)可以寫作:

3)使用變換Zi=wi·，式(3)可以寫為

1.2 高維積分的準(zhǔn)Monte Carlo求解方法

對(duì)于高維積分的求解，像Monte Carlo方法一樣，準(zhǔn)Monte Carlo方法的計(jì)算誤差與積分維數(shù)和積分區(qū)域無關(guān)，準(zhǔn)Monte Carlo方法可以得到真實(shí)誤差，且已經(jīng)證明這個(gè)誤差要比概率誤差?。?1］.所以選用準(zhǔn)Monte Carlo方法進(jìn)行高維正態(tài)積分的求解.式(4)可以寫為

采用準(zhǔn)Monte Carlo方法，式(5)的估計(jì)值為［14］

估計(jì)值的方差為［14］

式中:(xi1、xi2、…、xim)表示積分區(qū)域上的m維確定性點(diǎn)，采用哈頓序列和數(shù)論網(wǎng)格點(diǎn)來表示.

1.2.1 哈頓序列

哈頓(Halton)序列［10］于 1960年提出，最初用于一維均勻序列的產(chǎn)生，后來被擴(kuò)展到了高維情況.

一維Halton序列的產(chǎn)生可以分成2步:

1)將非負(fù)整數(shù)i轉(zhuǎn)化為Pk進(jìn)制的整數(shù)，表示為

2)將Pk進(jìn)制的整數(shù)反轉(zhuǎn)變?yōu)镻k進(jìn)制的小數(shù)，再轉(zhuǎn)化為十進(jìn)制的小數(shù)，表示為

多維Halton序列可以通過產(chǎn)生多個(gè)一維序列獲得，隨著積分維數(shù)的變化Pk需要取不同的素?cái)?shù).對(duì)于m個(gè)隨機(jī)變量的積分，取N個(gè)積分點(diǎn)的情況，首先選取前m個(gè)素?cái)?shù)，產(chǎn)生一維Halton序列，然后再形成m維Halton序列，表示為

式中:i=0，1，…N －1.

1.2.2 數(shù)論網(wǎng)格法(Lattice method)

數(shù)論網(wǎng)格法［11］最早由柯洛波夫(Korobov)于1959年提出.對(duì)于m個(gè)隨機(jī)變量的積分，取N個(gè)網(wǎng)格點(diǎn)的情況，最簡(jiǎn)單的單位超立方網(wǎng)格定義為

式中:{·}表示取小數(shù)部分，(a1，a2，…，am)稱為最優(yōu)系數(shù).

按照定義，網(wǎng)格點(diǎn)的產(chǎn)生主要與最優(yōu)系數(shù)(a1，a2，…，am)有關(guān)，最優(yōu)系數(shù)(a1，a2，…，am)的選取與積分計(jì)算結(jié)果的誤差有關(guān)，式(6)的數(shù)論網(wǎng)格法(Lattice method)積分可以表示為

式中:E為誤差項(xiàng).

數(shù)論網(wǎng)格法(Lattice method)特別適用于被積函數(shù)為周期函數(shù)的情形，對(duì)于被積函數(shù)不是周期函數(shù)的情況，需要通過Fourier展開表示為周期函數(shù).所以誤差項(xiàng)是由Fourier級(jí)數(shù)形式表示的，最優(yōu)系數(shù)(a1，a2，…，am)要滿足誤差項(xiàng) E 最小的原則，最優(yōu)系數(shù)可以參考文獻(xiàn)［15］來選取，本文按照J(rèn)oe提出的準(zhǔn)則進(jìn)行最優(yōu)系數(shù)的選擇［16］.

1.3 計(jì)算精度比較

1.3.1 Benchmark 解

在各個(gè)失效模式的可靠性指標(biāo)相等以及失效模式之間的相關(guān)系數(shù)相等的情況下，即βi=β，rij=r，高維正態(tài)積分函數(shù)Φm(β，R)可以被化簡(jiǎn)為一維積分，表達(dá)式為［17］

將這個(gè)特殊的精確解作為Benchmark解，用來評(píng)估本文所提出的方法的準(zhǔn)確性.

1.3.2 兩種方法的比較

根據(jù)準(zhǔn)Monte Carlo方法，高維正態(tài)積分的求解主要與3個(gè)因素有關(guān):樣本數(shù)量、可靠性指標(biāo)和相關(guān)系數(shù).為了研究本文提出方法的準(zhǔn)確性，采用2個(gè)指標(biāo)作為衡量標(biāo)準(zhǔn):1)結(jié)果的變異系數(shù)(coefficient of variation，Cov);2)結(jié)果的誤差.其中，變異系數(shù)表示結(jié)果的收斂快慢，變異系數(shù)越小表示結(jié)果的收斂程度越好;誤差表示結(jié)算結(jié)果的準(zhǔn)確性，誤差越小說明結(jié)果越準(zhǔn)確.

根據(jù)式(6)和式(7)，結(jié)果的變異系數(shù)定義為

根據(jù)式(6)和式(13)，結(jié)果的誤差定義為

采用1 000 000個(gè)樣本，按照式(14)、(15)2種標(biāo)準(zhǔn)進(jìn)行了計(jì)算，計(jì)算結(jié)果如圖1～2所示.

從結(jié)果的收斂速度來說，3種方法的收斂速度隨可靠性指標(biāo)的增加而變快，隨相關(guān)系數(shù)的增大而減慢，Halton序列方法對(duì)結(jié)果的收斂速度影響不大，數(shù)論網(wǎng)格方法(Lattice method)的收斂速度遠(yuǎn)大于隨機(jī)抽樣.

圖1 不同β和r時(shí)的Cov計(jì)算結(jié)果Fig.1 The result of Cov for differentβ and r

從結(jié)果的準(zhǔn)確性來看，3種方法的準(zhǔn)確性都隨可靠性指標(biāo)的增加和隨相關(guān)系數(shù)的增大而降低.本文提出的2種準(zhǔn)Monte Carlo方法的準(zhǔn)確性遠(yuǎn)高于隨機(jī)抽樣的一般Monte Carlo方法，其中數(shù)論網(wǎng)格方法(Lattice method)的精度可以控制在3%以內(nèi)，Halton序列的精度可以控制在10%以內(nèi).

圖2 不同β和r時(shí)的Error計(jì)算結(jié)果Fig.2 The result of Error for different β and r

2 串聯(lián)結(jié)構(gòu)系統(tǒng)的可靠性計(jì)算

在結(jié)構(gòu)系統(tǒng)中任何一個(gè)模式失效都有可能導(dǎo)致整個(gè)結(jié)構(gòu)系統(tǒng)失效，則稱這種體系為串聯(lián)系統(tǒng).一般地，由m個(gè)失效模式組成的串聯(lián)系統(tǒng)，通常用圖3所示的符號(hào)表示［17］.

圖3 串聯(lián)結(jié)構(gòu)系統(tǒng)示意圖Fig.3 Sketch of the series structural system

對(duì)于串聯(lián)體系，設(shè) F1，F(xiàn)2，…，F(xiàn)m，表示各個(gè)模式的失效事件，則該系統(tǒng)的失效概率為

本文提出的2種準(zhǔn)Monte Carlo方法中，數(shù)論網(wǎng)格法(Lattice method)無論在收斂速度還是在計(jì)算精度上都比Halton序列方法好，在以后的計(jì)算中都采用數(shù)論網(wǎng)格法.

根據(jù)式(13)和式(16)，串聯(lián)結(jié)構(gòu)系統(tǒng)的失效概率可以表示為

式(17)可以用來評(píng)估串聯(lián)結(jié)構(gòu)系統(tǒng)可靠性數(shù)值方法計(jì)算的準(zhǔn)確性.采用式(18)來計(jì)算近似方法的誤差.

采用以上評(píng)估方法，將本文提出數(shù)論網(wǎng)格方法(Lattice method)和計(jì)算串聯(lián)結(jié)構(gòu)系統(tǒng)可靠性方法中的 G-FOMN 法［6］、PCM 法［7］、IPCM 法［9］以 及Equivalent components法［18］進(jìn)行了比較.為了研究可靠性指標(biāo) β和相關(guān)系數(shù) r的影響，取 m=20，1 000 000個(gè)樣本進(jìn)行計(jì)算，計(jì)算結(jié)果和其他方法的比較如圖4所示.

圖4 不同β和r時(shí)的計(jì)算結(jié)果比較Fig.4 Comparison of the results for different β and r

為了研究積分維數(shù)m和相關(guān)系數(shù)的影響r，取β=4，1 000 000個(gè)樣本進(jìn)行計(jì)算，計(jì)算結(jié)果和其他方法的比較如圖5所示.為了研究可靠性指標(biāo)β的影響，取r=0.9，1 000 000個(gè)樣本進(jìn)行計(jì)算，計(jì)算結(jié)果和其他方法的比較如圖5(c)和圖6所示.

從圖4可以看出，在 r=0.1和 r=0.5時(shí) GFOMN法、PCM法、IPCM法的計(jì)算誤差隨可靠性指標(biāo)的增大經(jīng)歷了一個(gè)由小到大再變小的過程，在r=0.9時(shí)這3種方法的計(jì)算誤差隨可靠性指標(biāo)的增大而增大;在r=0.5時(shí)Equivalent components法的計(jì)算誤差隨可靠性指標(biāo)的增大經(jīng)歷了一個(gè)由小到大再變小的過程，在r=0.1和r=0.9時(shí)Equivalent components法的誤差隨可靠性指標(biāo)的增大而減小;本文提出方法的計(jì)算誤差隨可靠性指標(biāo)的增大而增大.

從圖5可以看出，在r=0.1時(shí)本文方法和其他4種方法的計(jì)算誤差都隨積分維數(shù)的增大而增大;在r=0.5和在r=0.9時(shí)本文提出的方法、G-FOMN法、IPCM法的計(jì)算誤差隨積分維數(shù)的增大而增大，而PCM法和Equivalent components法的計(jì)算誤差隨積分維數(shù)的增大而減小.

從圖5(c)和圖6可以看出，在r=0.9時(shí)，不論積分維數(shù)取值是多少，本文提出的方法和其他4種方法的計(jì)算誤差都隨可靠性指標(biāo)的增大而增大.

綜合所有的計(jì)算結(jié)果來看，本文提出的數(shù)論網(wǎng)格方法(Lattice method)的計(jì)算誤差較小，計(jì)算精度遠(yuǎn)高于其他4種方法.

圖5 不同m和r時(shí)的計(jì)算結(jié)果比較Fig.5 Comparison of the results for different m and r

圖6 β=4，r=0.9時(shí)的計(jì)算結(jié)果比較Fig.6 Comparison of the calculation results for β =4，r=0.9

3 計(jì)算實(shí)例

我國(guó)現(xiàn)行的《潛水系統(tǒng)和潛水器入級(jí)與建造規(guī)范》，給出了耐壓環(huán)肋圓柱殼結(jié)構(gòu)的5種校核公式［19］.大量的潛水器結(jié)構(gòu)可靠性研究文獻(xiàn)將這5種校核公式作為失效模式的極限狀態(tài)方程.

失效模式1 肋骨跨中殼板縱剖面內(nèi)中面周向應(yīng)力強(qiáng)度不足，表達(dá)式為:G1=σs－.

失效模式2 肋骨跨端板殼內(nèi)表面縱向應(yīng)力強(qiáng)度不足，表達(dá)式為:G1=σs－σ1.

失效模式3 肋骨應(yīng)力強(qiáng)度不足，表達(dá)式為:G3=σs－σf.

失效模式4 相鄰肋骨間殼板失穩(wěn).表達(dá)式為:G4=Pcr－P.

失效模式5 艙段整體失穩(wěn)，表達(dá)式為:G5=－P.其中，符號(hào)的具體涵義以及隨機(jī)變量的分布和特征參數(shù)見文獻(xiàn)［13］.首先通過一次二階矩方法，計(jì)算出各個(gè)失效模式的可靠性指標(biāo)，再計(jì)算出失效模式之間的相關(guān)系數(shù)，然后通過式(15)求解系統(tǒng)的失效概率.各個(gè)失效模式的可靠性指標(biāo)如表1所示，失效模式之間的相關(guān)系數(shù)如表2所示.

表1 各失效模式的可靠性指標(biāo)Table 1 Reliability index of every failure modes

表2 各個(gè)失效模式之間的相關(guān)系數(shù)Table 2 The correlation among failure modes

將數(shù)論網(wǎng)格方法(Lattice method)的計(jì)算結(jié)果和Drezner方法［4］進(jìn)行了比較，比較結(jié)果如表3所示.從表3可以看出本文方法和Drezner方法的計(jì)算結(jié)果很接近，說明該方法的計(jì)算精度較高.

表3 計(jì)算結(jié)果的比較Table 3 Comparison of the calculation results

4 結(jié)束語(yǔ)

本文提出一種求解高維正態(tài)積分的準(zhǔn)Monte Carlo方法，并將此方法用于串聯(lián)結(jié)構(gòu)系統(tǒng)可靠性的計(jì)算中.與其他近似方法相比較，本文提出的方法具有與積分維數(shù)無關(guān)、受可靠性指標(biāo)和相關(guān)系數(shù)影響不大的特點(diǎn).計(jì)算實(shí)例表明，本文提出的方法計(jì)算精度高，計(jì)算結(jié)果穩(wěn)定，適應(yīng)性廣.同時(shí)，本文提出的高維正態(tài)積分計(jì)算方法也適用于并聯(lián)系統(tǒng)可靠性的計(jì)算.

［1］HOHENBICHLER M，RACKWITZ M.First-order concepts in system reliability［J］.Structural Safety，1983，1:177-188.

［2］貢金鑫.一類多維正態(tài)積分的近似解法［J］.西安建筑科技大學(xué)學(xué)報(bào)，1996，28(1):52-56.GONG Jinxin.Approximate solution for a class of multi-normal integral［J］.Journal of Xi’an University of Architecture and Technology，1996，28(1):52-56.

［3］MILTON R C.Computer evaluation of the multivariate normal integral［J］.Technometrics，1972，14(4):881-889.

［4］DREZNER Z.Computation of the multivariate normal integral［J］.ACM Transactions on mathematical Software，1992，18(4):881-889.

［5］馬永亮，曲先強(qiáng)，崔洪斌，等.基于數(shù)值積分方法的潛水器耐壓圓柱殼結(jié)構(gòu)系統(tǒng)可靠性研究［C］//黑龍江省造船工程學(xué)會(huì)學(xué)術(shù)年會(huì)，哈爾濱，中國(guó)，2010:92-96.

［6］TANG L K，MELCHERSR E.Improved approximation for multi-normal integral［J］.Structural Safety，1987，4:81-93.

［7］PANDEY M D.An effective approximation to evaluate multinormal integrals［J］.Structural Safety，1998，20:51-67.

［8］YASUHIRO M，TERUYUKI K.Multi-normal integrals by importance sampling for series system reliability［J］.Structural Safety，2003，25:363-378.

［9］YUAN X X，PANDEY M D.Analysis of approximations for multi-normal integration in system reliability computation［J］.Structural Safety，2006，28:361-377.

［10］GLASSERMAN P.Monte Carlomethods in financial engineering［M］.Beijing:Higher Education Press，2008:293-297.

［11］宮野.計(jì)算物理［M］.大連:大連工學(xué)院出版社，1987:275-292.

［12］GENZ A.Numerical computation of multivariate normal probabilities［J］.Journal of Comput Graph Stat，1992，1:141-150.

［13］馬永亮.考慮腐蝕影響的潛艇結(jié)構(gòu)可靠性研究［D］.哈爾濱:哈爾濱工程大學(xué)，2009:14-18.MA Yongliang.Reliability assessment of submarine structure considering corrosion［D］.Harbin:Harbin Engineering University，2009:14-18.

［14］DAVISP，RABINOWITZP.數(shù)值積分法［M］.馮振興，伍富良，譯.北京:高等教育出版社，1986:235-238.

［15］馮康.數(shù)值計(jì)算方法［M］.北京:國(guó)防工業(yè)出版社，1978:82-84.

［16］JOE S，SLOAN IH.Implementation ofa lattice method for numerical multiple integration［J］.ACM Transactions on Mathematical Software，1993，19(4):523-545.

［17］何水清，王善.結(jié)構(gòu)可靠性分析與設(shè)計(jì)［M］.北京:國(guó)防工業(yè)出版社，1993:208-216.

［18］GOLLWITZER S，RACKWITZ R.Equivalent components in first-order system reliability［J］.Structural Safety，1983，5(2):359-366.

［19］中國(guó)船級(jí)社.潛水系統(tǒng)和潛水器入級(jí)與建造規(guī)范［S］.北京:人民交通出版社，1996:32-42.