基于高維正態積分的串聯結構系統可靠性研究

馬永亮，曲先強，崔洪斌，石德新

(哈爾濱工程大學 多體船技術國防重點學科實驗室，黑龍江 哈爾濱 150001)

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

1 高維正態積分求解的準Monte Carlo方法

1.1 高維正態積分表達式及積分變換

式中:

用于串聯結構系統可靠性計算的高維正態積分表達式為

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

高維正態積分Φm(β，R)需要采用積分變換和數值方法進行求解.由于Φm(β，R)的積分域為半無限域，很難進行積分求解，所以進行合適的積分變換是必要的.通過以下3步變換可以將高維正態積分轉化為單位超立方域內的高維積分［12］.

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

式中:φ(y1)為標準正態分布概率密度函數，為積分上限 β＇i

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

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

1.2 高維積分的準Monte Carlo求解方法

對于高維積分的求解，像Monte Carlo方法一樣，準Monte Carlo方法的計算誤差與積分維數和積分區域無關，準Monte Carlo方法可以得到真實誤差，且已經證明這個誤差要比概率誤差?。?1］.所以選用準Monte Carlo方法進行高維正態積分的求解.式(4)可以寫為

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

估計值的方差為［14］

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

1.2.1 哈頓序列

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

一維Halton序列的產生可以分成2步:

1)將非負整數i轉化為Pk進制的整數，表示為

2)將Pk進制的整數反轉變為Pk進制的小數，再轉化為十進制的小數，表示為

多維Halton序列可以通過產生多個一維序列獲得，隨著積分維數的變化Pk需要取不同的素數.對于m個隨機變量的積分，取N個積分點的情況，首先選取前m個素數，產生一維Halton序列，然后再形成m維Halton序列，表示為

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

1.2.2 數論網格法(Lattice method)

數論網格法［11］最早由柯洛波夫(Korobov)于1959年提出.對于m個隨機變量的積分，取N個網格點的情況，最簡單的單位超立方網格定義為

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

按照定義，網格點的產生主要與最優系數(a1，a2，…，am)有關，最優系數(a1，a2，…，am)的選取與積分計算結果的誤差有關，式(6)的數論網格法(Lattice method)積分可以表示為

式中:E為誤差項.

數論網格法(Lattice method)特別適用于被積函數為周期函數的情形，對于被積函數不是周期函數的情況，需要通過Fourier展開表示為周期函數.所以誤差項是由Fourier級數形式表示的，最優系數(a1，a2，…，am)要滿足誤差項 E 最小的原則，最優系數可以參考文獻［15］來選取，本文按照Joe提出的準則進行最優系數的選擇［16］.

1.3 計算精度比較

1.3.1 Benchmark 解

在各個失效模式的可靠性指標相等以及失效模式之間的相關系數相等的情況下，即βi=β，rij=r，高維正態積分函數Φm(β，R)可以被化簡為一維積分，表達式為［17］

將這個特殊的精確解作為Benchmark解，用來評估本文所提出的方法的準確性.

1.3.2 兩種方法的比較

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

根據式(6)和式(7)，結果的變異系數定義為

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

采用1 000 000個樣本，按照式(14)、(15)2種標準進行了計算，計算結果如圖1～2所示.

從結果的收斂速度來說，3種方法的收斂速度隨可靠性指標的增加而變快，隨相關系數的增大而減慢，Halton序列方法對結果的收斂速度影響不大，數論網格方法(Lattice method)的收斂速度遠大于隨機抽樣.

圖1 不同β和r時的Cov計算結果Fig.1 The result of Cov for differentβ and r

從結果的準確性來看，3種方法的準確性都隨可靠性指標的增加和隨相關系數的增大而降低.本文提出的2種準Monte Carlo方法的準確性遠高于隨機抽樣的一般Monte Carlo方法，其中數論網格方法(Lattice method)的精度可以控制在3%以內，Halton序列的精度可以控制在10%以內.

圖2 不同β和r時的Error計算結果Fig.2 The result of Error for different β and r

2 串聯結構系統的可靠性計算

在結構系統中任何一個模式失效都有可能導致整個結構系統失效，則稱這種體系為串聯系統.一般地，由m個失效模式組成的串聯系統，通常用圖3所示的符號表示［17］.

圖3 串聯結構系統示意圖Fig.3 Sketch of the series structural system

對于串聯體系，設 F1，F2，…，Fm，表示各個模式的失效事件，則該系統的失效概率為

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

根據式(13)和式(16)，串聯結構系統的失效概率可以表示為

式(17)可以用來評估串聯結構系統可靠性數值方法計算的準確性.采用式(18)來計算近似方法的誤差.

采用以上評估方法，將本文提出數論網格方法(Lattice method)和計算串聯結構系統可靠性方法中的 G-FOMN 法［6］、PCM 法［7］、IPCM 法［9］以 及Equivalent components法［18］進行了比較.為了研究可靠性指標 β和相關系數 r的影響，取 m=20，1 000 000個樣本進行計算，計算結果和其他方法的比較如圖4所示.

圖4 不同β和r時的計算結果比較Fig.4 Comparison of the results for different β and r

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

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

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

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

綜合所有的計算結果來看，本文提出的數論網格方法(Lattice method)的計算誤差較小，計算精度遠高于其他4種方法.

圖5 不同m和r時的計算結果比較Fig.5 Comparison of the results for different m and r

圖6 β=4，r=0.9時的計算結果比較Fig.6 Comparison of the calculation results for β =4，r=0.9

3 計算實例

我國現行的《潛水系統和潛水器入級與建造規范》，給出了耐壓環肋圓柱殼結構的5種校核公式［19］.大量的潛水器結構可靠性研究文獻將這5種校核公式作為失效模式的極限狀態方程.

失效模式1 肋骨跨中殼板縱剖面內中面周向應力強度不足，表達式為:G1=σs－.

失效模式2 肋骨跨端板殼內表面縱向應力強度不足，表達式為:G1=σs－σ1.

失效模式3 肋骨應力強度不足，表達式為:G3=σs－σf.

失效模式4 相鄰肋骨間殼板失穩.表達式為:G4=Pcr－P.

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

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

表2 各個失效模式之間的相關系數Table 2 The correlation among failure modes

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

表3 計算結果的比較Table 3 Comparison of the calculation results

4 結束語

本文提出一種求解高維正態積分的準Monte Carlo方法，并將此方法用于串聯結構系統可靠性的計算中.與其他近似方法相比較，本文提出的方法具有與積分維數無關、受可靠性指標和相關系數影響不大的特點.計算實例表明，本文提出的方法計算精度高，計算結果穩定，適應性廣.同時，本文提出的高維正態積分計算方法也適用于并聯系統可靠性的計算.

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

［2］貢金鑫.一類多維正態積分的近似解法［J］.西安建筑科技大學學報，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］馬永亮，曲先強，崔洪斌，等.基于數值積分方法的潛水器耐壓圓柱殼結構系統可靠性研究［C］//黑龍江省造船工程學會學術年會，哈爾濱，中國，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］宮野.計算物理［M］.大連:大連工學院出版社，1987:275-292.

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

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

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

［15］馮康.數值計算方法［M］.北京:國防工業出版社，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］何水清，王善.結構可靠性分析與設計［M］.北京:國防工業出版社，1993:208-216.

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

［19］中國船級社.潛水系統和潛水器入級與建造規范［S］.北京:人民交通出版社，1996:32-42.