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

李靜輝 康 銳

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

Ali Mosleh

(美國(guó)馬里蘭大學(xué)帕克分校 風(fēng)險(xiǎn)與可靠性中心，馬里蘭 20742)

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

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

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

1 導(dǎo)數(shù)估計(jì)方法選擇

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

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

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

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

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

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

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

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

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

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

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

2 原始蒙特卡羅似然比導(dǎo)數(shù)估計(jì)

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

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

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

將式(4)對(duì)λj求導(dǎo)，可得

從而

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

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

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

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

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

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

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

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

進(jìn)一步，定義如下估計(jì)量:

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

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

式(15)意味著可以用

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

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

對(duì)于 r=j，有

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

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

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

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

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

4 實(shí)例驗(yàn)證

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

圖2 一個(gè)簡(jiǎn)單的串并聯(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，任務(wù)時(shí)間為 T=1 h.對(duì)于該簡(jiǎn)單系統(tǒng)，不難解析計(jì)算出系統(tǒng)不可靠度及其對(duì)各元件失效率的導(dǎo)數(shù)，相應(yīng)的計(jì)算公式為

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

對(duì)于各組蒙特卡羅仿真試驗(yàn)，表中第1行為估計(jì)值，第2行為估計(jì)方差和相對(duì)誤差(均為樣本值)，其中相對(duì)誤差根據(jù)下式計(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ù)的基礎(chǔ)上乘上一定的因子(0.2，0.5，1，2，5)得到，第⑤組在各元件自然失效率的基礎(chǔ)上同時(shí)增大104倍.原始蒙特卡羅的仿真次數(shù)為108，5組偏倚蒙特卡羅的仿真次數(shù)均為106.

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

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

5 結(jié)論

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

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

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

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