加權隨機模擬的時變可靠性分析方法*

袁修開，鄭振軒，羅冬秀

(廈門大學 航空航天學院， 福建 廈門 361005)

在工程系統的可靠性分析中，時變不確定性廣泛存在，如材料強度會隨著工作時間發生變化。因此，系統性能的極限狀態函數是一個復雜的隨機過程。關于時變可靠性分析問題的研究早在20世紀40年代就開始進行了，目前的時變可靠性研究方法主要可歸納為三大類，分別為跨越率法[1-7]、極值方法[8-16]以及復合極限狀態函數求解法[17-21]等。

跨越率法是在20世紀40年代Rice提出了首次超越概率計算公式[1]。之后，Coleman給出了關于首次超越概率的可靠性分析方法[2]。Crandall又將數值模擬的方法引入跨越率法計算之中[3]。在近幾年，Hu和Du提出了一種聯合跨越率的方法求解時變可靠性問題[4]。Andrieu-Renaud等基于跨越率法提出一種計算方法，即PHI2方法[5]。Sudret在PHI2方法的基礎上，給出了PHI2方法的近似表達式來求解跨越率[6]。由于在實際工程中一般只能得知不確定性參數的區間范圍，故張德權等提出了關于區間的PHI2方法[7]。但跨越率法依賴于穿越過程假設，這會帶來精度上的誤差。

極值方法主要分為兩種，一種是將隨機過程近似為極值分布[8-10]，即對隨機過程在一段時間內的極大值進行統計，獲取極大值的統計參數，之后再運用靜態可靠性分析方法對其進行求解。如貢金鑫等[8]提出考慮抗力變化的結構可靠度分析方法，其將隨機過程近似為極值Ⅰ型分布，推導出關于時變可靠性的近似表達式。王新剛等[9]同樣也是采用將隨機過程近似為極值Ⅰ型分布進行求解。黃新萍等[10]基于貢金鑫的方法對非線性結構進行了時變可靠性分析。Li等[11-12]同樣也是通過將隨機過程轉化為極值問題來對時變可靠性進行求解。Hu和Du[13]提出一種通過統計隨機過程的極值分布規律，得出極值的分布參數，從而將時變可靠性問題轉化為靜態可靠性問題進行求解。近年來，有另外一種關于極值的方法是建立代理模型[14-16]。通過建立輸入隨機變量(隨機過程)與輸出響應極大值之間的代理模型來對時變可靠性進行求解。如Zhang等[16]通過將隨機過程離散成隨機變量，計算出結構響應并且取其極大值，建立起輸入參數和響應極大值之間的代理模型，最后通過一次二階矩進行計算。

復合極限狀態函數求解法[17-22]是首先將隨機過程離散化[20]，之后再運用時變可靠性中的瞬時失效概率和串聯系統可靠性將時變可靠性問題轉為靜態可靠性問題求解。Zissimos等[17-18]通過運用全概率定理和復合極限狀態的概念，針對輸入隨機變量和輸入隨機過程的時變可靠性分析，提出了一種新的可靠性分析方法。Li等[22]將單響應的時變可靠性通過復合極限狀態函數法轉化為多響應的靜態可靠性分析，再運用GSS(generalized subset simulation)方法計算得到累積失效概率，最后再通過插值的方法得到累積失效概率曲線。

本文在極值方法的基礎上，結合加權思想[23-25]，提出基于加權思想的時變可靠性分析方法，分別為加權蒙特卡洛(Weighted Monte Carlo Simulation， WMCS)算法和加權重要抽樣(Weighted Importance Sampling， WIS)算法，所提方法通過一次可靠性分析即可將時變失效概率表達成樣本的函數形式，從而大大提高時變可靠性分析效率。通過算例分析，將其結果與直接蒙特卡洛仿真的真實值對比來驗證基于加權思想的可靠性分析方法的準確性和高效性。

1 時變可靠性分析模型

在機械結構的可靠性分析中，依據應力-強度理論，時變機械結構響應Z(t)的功能函數(亦稱為極限狀態函數)g(·),形式可定義為：

Z(t)=g(R,Q(t),G)

(1)

其中：G是固定載荷或隨機載荷；Q(t)是可變載荷的隨機過程，為時間t的函數；R是機械結構的抗力值，本文中假設其為隨機變量。在設計基準期T內的累積失效概率可以表示為：

PF(T)={g(R,Q(t),G)≤0,?t∈[0,T]}

(2)

其中，?表示“至少存在一個”。

時變可靠性的極值方法是通過時間離散后，采用極值分布來將隨機過程等效為隨機變量，從而將時變問題轉化成靜態可靠性問題。以Q(t)轉化為例，首先將設計基準期T分為m個時段，每個時段為δ=T/m，則將Q(t)離散為m個參數為αδ和μδ的極值Ⅰ型分布隨機變量Qi(i=1,…,m)，然后求解Qi的極值，QT=max(Qi)，其對應的分布參數為：

(3)

則式(2)的時變累積失效概率可轉化為：

PF(T)={g(R,QT,G)≤0}

(4)

更進一步，任意時間段[0,t]下的時變累積失效概率函數為：

PF(t)={g(R,Qt,G)≤0}

(5)

其中，Qt為[0,t]上Qi的極值。可看出式(4)中的PF(T)為t=T時的累積失效概率。若要求解不同時間段內的失效概率，即失效概率與時間的函數關系，則一般需重復多次進行可靠性分析。

本文提出基于加權思想的時變可靠性分析方法。該方法通過一次常規可靠性分析模擬，將時變失效概率函數表達成樣本的加權函數形式，即可高效得到不同時段內的時變可靠性結果。

2 基于加權思想的可靠性分析方法

本文將采用加權思想，結合蒙特卡洛法和重要抽樣法提出加權蒙特卡洛法和加權重要抽樣法高效求解時變可靠性問題，所提方法的優勢在于僅需一次可靠性分析。

2.1 加權蒙特卡洛法

加權思想期望用一次可靠性分析來推斷不同時間段的失效概率，從而避免重復進行可靠性分析。

首先，構造一個基本隨機變量x=[R,Qt,G]的“輔助概率密度函數”f(x|t0)，其對應某一時刻t0變量的分布。式(5)的時變累積失效概率函數可轉化為：

=Et0[IF(x)rMCS(x,t,t0)]

(6)

其中：IF(·)為系統失效的指數函數；f(x|t)為對應任意t∈[0,T]時刻的變量分布；rMCS為“權重因子”，為概率密度函數的比值，即

(7)

通過式(6)可以看出t時刻的失效概率可以表示為該時刻權重因子在概率密度函數f(x|t0)的期望形式，時變失效概率的無偏估計可以表示為：

(8)

其中，x(j)為采用f(x|t0)抽取的樣本。進一步地，可得到失效概率函數估計的方差：

(9)

由于式(8)中對于不同的時間t用的是同一組隨機樣本x(j)(j=1,2,…,N)(來自f(x|t0))和相同的指示函數IF(x(j))，因此無須重新調用功能函數。在估計過程中，需要重復計算的僅是各時刻下的權重系數，而權重系數只是概率密度函數的比值，并不需要結構系統的分析，無須太大的計算代價，因此所提方法可以大大減小時變可靠性分析的工作量，提高計算效率。

前面提到f(x|t0)是輔助概率密度函數，它所產生的樣本將用于不同時刻失效概率的估計，因此需合理選擇f(x|t0)，盡可能地考慮所有時刻情況。在構建輔助概率密度函數時，首先把與時間無關的隨機變量歸為一組，記為xr=[x1,x2,…,xnr]，其余為xs=[x1,x2,…,xns]，則輔助概率密度函數可以表示為：

f(x|t0)=f(xr|θt0)f(xs)

(10)

其中，θt0表示與時間相關的隨機變量的分布參數在某一時刻t0的取值，如一般可以取t0=T/2。

2.2 加權重要抽樣法

加權蒙特卡洛法仍然受限于蒙特卡洛法的效率，為了提高計算效率，本文進一步提出加權重要抽樣法。采用加權的思想，可以將重要抽樣用于時變可靠性分析。

構造一個“全局”重要抽樣密度函數H(x)，式(5)的失效概率函數可以表示為：

(11)

在抽取H(x)樣本x(j)(j=1,2,…,N)后，失效概率函數的無偏估計可表示為：

(12)

類似地，可得到失效概率函數估計的方差：

(13)

本文在構建H(x)時借鑒基于設計點的構造方法。首先計算得到設計點x*(如可采用改進一次二階矩法)，然后選取正態分布密度作為抽樣密度，以設計點x*作為抽樣密度中心，則H(x)可以表示為：

H(x)=φ(x|x*,σH)

(14)

其中，φ(·)表示正態概率密度函數，σH表示隨機變量的標準差。

3 數值算例分析

3.1 管狀懸臂梁

算例來自參考文獻[10]，圖1為一管狀懸臂梁，結構基本參數與所受外力情況見表1，預期工作時間為10 a。其中承受外載P、F及U均為隨機變量，Q(t)為隨機過程，這里給出其在10 a內的極值QT=max(Qi)的分布。其他參數為確定值，即L1=120 mm，L2=60 mm,θ1=5°,θ2=10°,功能函數定義為強度R與固定端下表面最大應力值之差，即

G(x,t)=R-σmax(x,t)

(15)

圖1 管狀懸臂梁示意圖Fig.1 Tubular cantilever beam

表1 管狀懸臂梁結構隨機變量分布Tab.1 Random variable distribution of tubular cantilever beam structure

其中σmax(x)可通過以下計算獲得：

(16)

(17)

(18)

計算該懸臂梁在工作期內的失效概率隨時間變化的函數。

采用本文所提方法來求解結構時變可靠性，加權蒙特卡洛法的“輔助概率密度函數”f(x|t0)在t0為第5 a時構造，對應的變量均值為μx=[440,6.5,22,90,42,5,2.302]。加權重要抽樣法中“全局”重要抽樣密度函數H(x)為在上述均值點下計算出設計點x*=[245,6.948,2.235,90.04,4.19,5.0,2.435]后進行構建。

加權蒙特卡洛抽樣105次，加權重要抽樣20+103次(其中采用改進一次二階矩方法求解設計點調用極限狀態方程20次)。直接蒙特卡洛法對10個時間點進行計算，每個時間點抽取106個樣本點，共調用功能函數10×106次，作為精確結果進行對比。各方法結果如圖2所示(其中N為樣本數)。由圖2可以看出：該懸臂梁失效概率隨時間單調遞增，所提方法結果都接近于真實值。表2給出了所提兩種方法對應的估計值的誤差及變異系數。由此可看出，加權方法只需一次分析，具有較高的效率。進一步地，加權重要抽樣法比加權蒙特卡洛法的效率更高。

圖2 管狀懸臂梁1～10 a失效概率圖(取t0=5)Fig.2 Failure probability diagram of tubular cantilever beam in 1～10 years(t0=5)

表2 管狀懸臂梁1～10 a的失效概率(t0=5)Tab.2 Failure probability of tubular cantilever beams from 1 to 10 years (t0=5)

其次，為了進一步考察方法的穩健性，選取不同的t0，分別再用加權蒙特卡洛法和加權重要抽樣法進行分析，結果如圖3所示。由圖3可以看出，兩種加權方法都獲得了較為滿意的結果。

圖3 管狀懸臂梁1～10 a失效概率圖(取t0=1)Fig.3 Failure probability diagram of tubular cantilever beam in 1～10 years(t0=1)

3.2 十桿桁架

圖4所示為十桿桁架結構，左端是固定端，各桿的截面積分別為A1～A10，其中1、2、3、4、5和6號桿長度為L；彈性模量為E。該桁架承受3個外載荷，其中F1、F3為隨機變量載荷，F2(t)為隨機過程載荷，其在使用期限T為10 a的極值F2T的分布列于表3中。要求節點2處的豎直方向位移在使用期限T為10 a內且不得大于允許值d，計算該桁架在使用期限內的失效概率隨時間變化的函數。

圖4 十桿桁架結構Fig.4 Ten-bar truss structure

表3 十桿桁架隨機變量分布Tab 3 Distribution information of ten-bar truss random variable

由題意可知，十桿桁架的功能函數為：

G(t)=d-dy(t)

(19)

其中，dy(t)表示t時刻下2點位置的實際豎直方向位移，可以通過下列公式求得。

(20)

(21)

(22)

(24)

(25)

(26)

(27)

(28)

(29)

(30)

(31)

(32)

(33)

(34)

采用本文加權方法來進行分析，其中取t0為10 a。加權蒙特卡洛法的“輔助概率密度函數”f(x|t0)的均值為μx=[0.004,…,0.004, 9.1,70, 12, 20, 8, 0.02]。構造的“全局”重要抽樣密度函數H(x)為在設計點x*=[0.004,…，0.004,9.405，67.49，12.08，20.67，7.97，0.02]下構造。

采用所提方法分析得到十桿桁架的失效概率結果如圖5和表4所示。直接MCS方法抽取10×106個樣本點計算獲得的結果作為精確結果。WMCS抽取105和WIS抽取103個樣本點。從圖5可看出，失效概率隨時間增長逐漸增大，均得到較為精確的結果，證實了所提加權方法的正確性。結合表4可以看出，加權方法具有較高的效率和精度。

圖5 十桿桁架失效概率Fig.5 Failure probability diagram of ten-bar truss

表4 十桿桁架1～10 a的失效概率(t0=5)Tab.4 Failure probability of ten-bar truss from 1 to 10 years(t0=5)

4 結論

本文針對隨機模擬法在時變可靠性分析中計算效率低的問題，在極值方法的基礎上結合加權思想，提出了加權蒙特卡洛方法和加權重要抽樣方法。所提方法僅需一次常規可靠性分析模擬，即可得到時變失效概率函數結果，大大提升了分析的效率。對文中兩個算例進行了驗證。算例分析結果表明，加權方法能夠在一次可靠性分析下，得到較為滿意的可靠性結果，且加權重要抽樣方法比加權蒙特卡洛方法效率更高。