可靠性靈敏度分析的一種偏倚蒙特卡羅方法

李靜輝 康 銳

(北京航空航天大學 可靠性與系統工程學院，北京 100191)

Ali Mosleh

(美國馬里蘭大學帕克分校 風險與可靠性中心，馬里蘭 20742)

在可靠性分析中，除獲得系統不可靠度等指標之外，通常還希望知道這些指標對相關參數(如元件的失效率和修復率等)的靈敏程度(相關導數).

可靠性靈敏度分析的一種傳統方法是進行解析地求導計算，但這種解析方法一般只對比較簡單的系統才可行，對于規模較大和較復雜的系統，通常需要借助于仿真.與系統不可靠度等非導數指標的估計相比，如何通過仿真獲得相關導數的估計目前研究相對較少.一種比較直接的方法是有限差分法［1］，但當感興趣的參數較多時，有限差分法需要運行很多輪仿真;此外，差分步長的選取涉及著名的“偏差－方差”難題.文獻［2］中探討了一種只需運行一輪仿真的方法，但仍存在差分步長的選取問題.而且對于高可靠度系統，原始蒙特卡羅仿真的效率一般會很低，需要采取降低方差技巧來加速仿真.文獻［3－5］等針對考慮修復的高可靠度馬爾可夫型系統的可靠性靈敏度分析進行了討論，并考慮了方差問題，但并沒有專門針對導數估計方差降低技巧的探討，而是直接利用降低非導數指標估計方差的技巧來同時減小導數估計的方差.

本文立足經典可靠性系統，探討應用仿真進行可靠性靈敏度分析的方法，并提出一種直接從降低導數估計方差出發的偏倚技巧來加速仿真.

1 導數估計方法選擇

本文主要考慮系統不可靠度(也稱為失效概率或累計失效概率)及其對相關參數導數的估計.一般地，系統不可靠度可表達成如下數學期望的形式:

式中，θ為參數向量;I{ω∈R}為指示函數;Ω為樣本空間;ω為樣本點;R為樣本空間中的系統失效域;f(θ，ω)為Ω上的概率密度函數.根據式(1)，系統不可靠度?(θ)可采用標準的蒙特卡羅仿真進行估計:

進一步，需要估計?(θ)對相關參數的導數:

θj為參數向量 θ =(θ1，θ2，…，θd)的第 j個分量;d為參數的個數.為簡化起見，下文有時將??(θ)/?θj簡記為.

本文關心的系統不可靠度對相關參數的導數估計具有如下幾個需求和特點:

1)感興趣的參數出現在概率分布中;

2)涉及的參數個數可能較多;

3)?(θ)為一指示函數的數學期望，而指示函數具有不光滑、難以變換處理等病態性質;

4)對于高可靠度系統，為獲得滿意精度的估計，可能需要運行大量次數的仿真.

目前文獻中已經提出一些可以用來估計如式(1)所示的數學期望對相關參數導數的仿真方法［6－8］，主要有:有限差分方法［1］、擾動分析方法［9－10］、似然比方法［11－12］和弱導數方法［13－15］.這些導數估計方法近年來在排隊系統、庫存模型和金融等領域中受到廣泛關注，但對于這些方法在可靠性分析中的應用，目前相關探討仍然較少.

在上述各種導數估計方法中，似然比方法較好地滿足本文關心的可靠性靈敏度分析的需求和特點.因此，本文選擇似然比方法作為基礎的導數估計方法.關于似然比方法的詳細介紹，可參考文獻［11－12］.

2 原始蒙特卡羅似然比導數估計

本節首先推導未采用重要抽樣的原始蒙特卡羅仿真環境下的似然比導數估計方法.在應用蒙特卡羅仿真進行可靠性分析時，可選擇兩種不同的基本形式:基于元件的直接蒙特卡羅方法和基于系統的間接蒙特卡羅方法［16］.本文主要討論基于元件的直接蒙特卡羅方法，并主要考慮如下類型的經典可靠性系統(模型):設系統由n個元件組成，各元件依次編號為1，2，…，n，初始時所有元件均處于正常工作狀態，假設各元件之間相互獨立，且壽命均服從指數分布，即 fj(t)=λje－λjt(λj為元件j的失效率)，任務時間為T，感興趣的量為系統不可靠度及其對各元件失效率的導數.

基于元件的直接蒙特卡羅方法從元件層次來模擬系統的行為，直接從相應的分布產生每個元件轉移時間的抽樣.記元件j抽樣到的轉移時間為tj(j=1，2，…，n)，由全部 tj可以確定一個具體的系統軌跡ω，對于每一個ω，可以判斷系統最終是成功還是失效，從而獲得I{ω∈R}的值(0或1).在整個仿真結束之后，計算式(2)可獲得系統不可靠度的估計.

為利用似然比方法獲得系統不可靠度對各元件失效率的導數，可取 ω =(t1，t2，…，tn)，ω 的概率密度為

將式(4)對λj求導，可得

從而

式(6)中Sλj(θ，ω)為統計學中著名的得分函數.根據似然比方法原理，可獲得?λj的無偏估計量:

3 一種利用結構函數的偏倚技巧

原始蒙特卡羅方法的一個主要缺陷是收斂速度慢，尤其對小概率事件，為獲得滿意精度的估計，需要運行大量次數的仿真.對比式(2)和式(7)可以看出，導數估計比系統不可靠度本身的估計更具有挑戰性.對于高可靠度系統，原始蒙特卡羅仿真采樣到的絕大多數系統軌跡ω都不會落入系統失效域R內，這些系統軌跡對各?λj貢獻的統計得分也都為零.而且，與估計不可靠度?的式(2)相比，估計?λj的式(7)還多了一個隨機項——得分函數，從而一般會帶來更大的方差.可見，導數估計比不可靠度估計本身一般更需要減小方差，以能在合理的時間內獲得滿意精度的估計.

為降低式(7)中系統不可靠度對各元件失效率導數估計的方差，本文在借鑒文獻［17－18］中方法的基礎上提出一種利用系統結構函數的偏倚技巧.

設Ij(tj)代表元件 j的狀態指示變量，Isys(tsys)代表系統的狀態指示變量，其中tj和tsys分別代表元件 j和系統的失效時間，Ij(tj)和Isys(tsys)均定義如下:

一般地，系統狀態指示變量可表達成各元件狀態指示變量的函數，該函數通常稱為系統結構函數.通過對結構函數進行邏輯化簡，可表達為

其中，m為結構函數化簡后的總項數;ci為第i項前面的系數;aji取0或1:如果第 i項中含有Ij(tj)，則 aji=1，否則 aji=0.

由于在基于元件的直接蒙特卡羅方法中，直接從相應的概率分布fj(t)抽樣每個元件的失效時間tj，應用重要抽樣的一個自然做法是改變tj的抽樣分布.設抽樣tj的新的重要抽樣分布為(t)，則相應的權重函數為

進一步，定義如下估計量:

值得一提的是，在上述證明過程中用到了以下兩個等式:

可以驗證，在 aji取0和1兩種情況下，均有式(16)和式(17)成立.

式(15)意味著可以用

作為?λj的近似估計值.

即Br代表式(12)中所有包含的項(亦即結構函數式(9)中所有包含Ir的項)，并可以將Br中的元素按從小到大的順序排列.進一步，對于r≠j，有

對于 r=j，有

對式(22)和式(23)作進一步整理，可得

由于式(24)和式(25)中第2項總大于0，故有

這意味著，式(19)中系統層次的優化問題可以分解為式(26)和式(27)中元件層次的優化問題，這是式(12)中所定義估計量的一個重要優勢，可以避免高維優化的難題.

對于壽命服從指數分布的情形，通過代入各相關變量(式(6)、式(10)、式(11))和一定的變形處理，式(26)和式(27)可分別轉化為如下2個非線性方程:

圖1 最優偏倚值隨自然值Λ的變化曲線

4 實例驗證

為驗證本文方法的有效性，本節考慮一個可以解析求解的簡單實例.系統的可靠性框圖如圖2所示.

圖2 一個簡單的串并聯混合系統

各元件的失效率為 λ1=1.0 ×10－6/h，λ2=1.5 × 10－6/h，λ3=2.0 × 10－6/h，λ4=8.0 ×10－4/h，λ5=9.0 ×10－4/h，任務時間為 T=1 h.對于該簡單系統，不難解析計算出系統不可靠度及其對各元件失效率的導數，相應的計算公式為

另外，采用MATLAB編程實現，應用第2節中描述的原始蒙特卡羅方法和第3節中描述的偏倚蒙特卡羅方法對該實例進行了求解.對于偏倚蒙特卡羅方法，除最優參數組外，還選取了其它幾組偏倚參數進行了試驗.解析計算和各組蒙特卡羅仿真試驗的結果見表1.

對于各組蒙特卡羅仿真試驗，表中第1行為估計值，第2行為估計方差和相對誤差(均為樣本值)，其中相對誤差根據下式計算:

5組偏倚蒙特卡羅仿真采用的偏倚參數分別為:

① (0.32，0.32，0.32，0.32，0.32)

② (1.60，1.60，1.60，1.60，1.60)

③ (8.00，8.00，8.00，8.00，8.00)

④ (0.32，0.80，1.60，3.20，8.00)

⑤ (0.01，0.015，0.02，8.00，9.00)

其中第②組為最優參數，第①，③，④組由在最優參數的基礎上乘上一定的因子(0.2，0.5，1，2，5)得到，第⑤組在各元件自然失效率的基礎上同時增大104倍.原始蒙特卡羅的仿真次數為108，5組偏倚蒙特卡羅的仿真次數均為106.

由表1可以看出，各組蒙特卡羅仿真試驗基本上都給出了較為合理的估計(各個量的估計值均比較接近其真值).此外，各組蒙特卡羅的估計性能也與預期相符:原始蒙特卡羅的效率較為低下，仿真108次獲得的結果精度仍不夠理想;5組偏倚蒙特卡羅的效率相對原始蒙特卡羅均有明顯的改善，其中對應最優參數組的偏倚蒙特卡羅符合預期地獲得了最小的方差和相對誤差，與原始蒙特卡羅相比，各個估計量都降低了至少6個數量級的方差.比較5組偏倚蒙特卡羅的估計方差和相對誤差可以看出偏倚參數的選取對于估計性能的影響.各組偏倚蒙特卡羅的估計結果基本上符合“偏倚參數離最優值越近，估計精度越高”的規律.特別地，如果一個偏倚參數組相對另一個偏倚參數組“絕對占優”(每個參數都離最優值更近)，仿真結果顯示其估計性能也“絕對占優”(對每個量的估計精度都更高).符合這個情況的有第①組相對第③組和第⑤組、第②組相對所有其余各組以及第④組相對第⑤組.為改善高可靠度系統仿真的效率，一般需要增大失效率，但第③組偏倚蒙特卡羅的試驗結果表明，選取太大的失效率(偏倚過度)獲得的估計性能可能并不理想(這一組中各個量的估計精度都比第①組和第④組要差至少一個數量級).另外，從第⑤組偏倚蒙特卡羅的仿真結果可以看到，將全部失效率統一擴大一定的倍數雖然是一種便捷的做法(而且在實際中常被仿真分析人員采用，尤其在缺少如何選取偏倚參數的有效指導的情況下)，但其估計性能可能并不是十分理想，特別是當各失效率的大小相差較大時.另一方面，雖然有些組中選取的偏倚參數離最優值不是很接近，但從整體上來看，各組偏倚蒙特卡羅仿真的估計性能都還較令人滿意(仿真106次獲得的相對誤差均小于5%).這主要源于所采用的利用系統結構函數的估計量，該估計量充分利用了系統結構信息，從而可以較為充分地利用仿真過程中產生的抽樣信息.另外，由表1可以看出的另一個值得一提的現象是，在同一次蒙特卡羅仿真試驗中，不同量的估計性能之間可能存在差異.對于這里考慮的實例，和的估計精度明顯比 ?，，和 ?λ3要低.

表1 解析計算和蒙特卡羅仿真估計結果

5 結論

1)似然比方法可以較好地滿足可靠性靈敏度分析(導數估計)的需求，能夠在一輪仿真中同時估計出系統不可靠度及其對相關參數的導數，并且僅需在估計系統不可靠度的基礎上增加較小的計算負擔;

2)對于高可靠度系統，似然比導數估計也面臨方差問題，需要采取降低方差技巧來加速仿真，事實上，導數量的估計通常比非導數量的估計更需要減小方差;

4)實例分析驗證了第3節中提出的偏倚蒙特卡羅方法的有效性，該方法對各個感興趣的量都給出了很好的估計，與原始蒙特卡羅方法相比，各個估計量的方差都降低了至少6個數量級.

References)

［1］Asmussen S，Glynn P W.Stochastic simulation:algorithms and analysis［M］.New York:Springer，2007

［2］Marseguerra M，Zio E.Monte Carlo estimation of the differential importance measure:application to the protection system of a nuclear reactor［J］.Reliability Engineering and System Safety，2004，86(1):11 －24

［3］Nakayama M K，Goyal A，Glynn P W.Likelihood ratio sensitivity analysis for Markovian models of highly dependable systems［J］.Operations Research，1994，42(1):137 －157

［4］Nakayama M K.Asymptotics of likelihood ratio derivative estimators in simulations of highly reliable Markovian systems［J］.Management Science，1995，41(3):524 －554

［5］Nakayama M K.On derivative estimation of the mean time to failure in simulations of highly reliable Markovian systems［J］.Operations Research，1998，46(2):285 －290

［6］L'Ecuyer P.An overview of derivative estimation［C］//Nelson B L，Kelton W D，Clark G M.Proceedings of the 1991 Winter Simulation Conference.Phoenix:IEEE，1991:207 －217

［7］Fu M C.Stochastic gradient estimation［C］//Henderson S G，Nelson B L.Handbooks in Operations Research and Management Science:Simulation.Boston:Elsevier，2006:575 －616

［8］Fu M C.What you should know about simulation and derivatives［J］.Naval Research Logistics，2008，55(8):723 －736

［9］Ho Y C.Performance evaluation and perturbation analysis of discrete event dynamic systems［J］.IEEE Transactions on Automatic Control，1987，32(7):563 －572

［10］Gong W B，Ho Y C.Smoothed(conditional)perturbation analysis of discrete event dynamical systems［J］.IEEE Transactions on Automatic Control，1987，32(10):858 －866

［11］Rubinstein R Y.Sensitivity analysis and performance extrapolation for computer simulation models［J］.Operations Research，1989，37(1):72 －81

［12］Rubinstein R Y，Kroese D.Simulation and the Monte Carlo method［M］.Hoboken:Wiley，2007

［13］Pflug G C.Sampling derivatives of probabilities［J］.Computing，1989，42(4):315 －328

［14］Heidergott B，Vázquez－Abad F J.Measure－valued differentiation for Markov chains［J］.Journal of Optimization Theory and Applications，2008，136(2):187 －209

［15］Heidergott B，Vázquez－Abad F J，Pflug G，et al.Gradient estimation for discrete－event systems by measure－valued differentiation［J］.ACM Transactions on Modeling and Computer Simulation，2010，20(1):1 －28

［16］Labeau P E，Zio E.Procedures of Monte Carlo transport simulation for applications in system engineering［J］.Reliability Engineering and System Safety，2002，77(3):217 －228

［17］Campioni L，Vestrucci P.Monte Carlo importance sampling optimization for system reliability applications［J］.Annals of Nuclear Energy，2004，31(9):1005 －1025

［18］Campioni L，Scardovelli R，Vestrucci P.Biased Monte Carlo optimization:the basic approach［J］.Reliability Engineering and System Safety，2005，87(3):387 －394