李靜輝 康 銳

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

Ali Mosleh

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

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

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

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

1 導數(shù)估計方法選擇

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

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

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

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

本文關(guān)心的系統(tǒng)不可靠度對相關(guān)參數(shù)的導數(shù)估計具有如下幾個需求和特點:

1)感興趣的參數(shù)出現(xiàn)在概率分布中;

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

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

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

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

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

2 原始蒙特卡羅似然比導數(shù)估計

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

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

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

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

從而

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

3 一種利用結(jié)構(gòu)函數(shù)的偏倚技巧

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

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

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

一般地，系統(tǒng)狀態(tài)指示變量可表達成各元件狀態(tài)指示變量的函數(shù)，該函數(shù)通常稱為系統(tǒng)結(jié)構(gòu)函數(shù).通過對結(jié)構(gòu)函數(shù)進行邏輯化簡，可表達為

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

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

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

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

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

式(15)意味著可以用

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

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

對于 r=j，有

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

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

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

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

圖1 最優(yōu)偏倚值隨自然值Λ的變化曲線

4 實例驗證

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

圖2 一個簡單的串并聯(lián)混合系統(tǒng)

各元件的失效率為 λ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.對于該簡單系統(tǒng)，不難解析計算出系統(tǒng)不可靠度及其對各元件失效率的導數(shù)，相應的計算公式為

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

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

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

① (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)

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

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

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

5 結(jié)論

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

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

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

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