蒙特卡羅方法在維修過程仿真中的應(yīng)用研究

趙思橋，劉井泉

（清華大學(xué)工程物理系，北京100084）

鑒于核電站對安全的特殊要求，對核電站設(shè)備進(jìn)行預(yù)防性維修已經(jīng)成為必不可少的一個環(huán)節(jié)。目前在有關(guān)核電站維修策略制定過程中，主要應(yīng)用解析計算的方法，然而鑒于電站設(shè)備壽命過程的隨機性，解析方法已經(jīng)難以滿足電站的經(jīng)濟性、安全性對維修計劃的要求。因此采用仿真方法求解維修策略已經(jīng)逐漸成為研究的熱點。目前并沒有針對于存在維修條件下設(shè)備壽命過程的仿真方法相關(guān)的研究。本文對各類維修條件下設(shè)備的壽命過程模型進(jìn)行了分析，對完全維修與最小維修的仿真方法進(jìn)行了總結(jié)并闡述了仿真方法的特點，本文中同樣對不完全維修條件下設(shè)備的壽命過程進(jìn)行了建模，并對其仿真方法進(jìn)行了構(gòu)造。

根據(jù)傳統(tǒng)分類可將維修過程分為三種類型：完全維修、最小維修和不完全維修。完全維修是指設(shè)備維修后恢復(fù)最初時刻的狀態(tài)，完全維修條件下設(shè)備逐次的故障間隔時間分布相互獨立同分布；最小維修是指設(shè)備在維修后僅恢復(fù)到維修前一時刻的狀態(tài)，設(shè)備的壽命過程服從非齊次泊松過程；不完全維修是介于完全維修與最小維修之間的一種維修，也就是設(shè)備維修后恢復(fù)到“修復(fù)如新”與“修復(fù)如舊”之間的某一個狀態(tài)。

對于維修過程的仿真，其主要任務(wù)在于對各種維修條件下設(shè)備的逐次失效間隔時間進(jìn)行抽樣。對于完全維修，鑒于每次的故障間隔時間分布相互獨立同分布，只需重復(fù)對相同的壽命分布進(jìn)行抽樣即可。對于最小維修，其抽樣就是對確定強度函數(shù)的非齊次泊松過程進(jìn)行抽樣，目前業(yè)界對非齊次泊松過程的抽樣方法有一定的研究，其主要思路是根據(jù)逐次失效間隔時間的聯(lián)合密度分布函數(shù)或非齊次泊松過程的特性構(gòu)造相應(yīng)的抽樣方法，關(guān)于最小維修的仿真在第2節(jié)會有詳細(xì)的說明。不完全維修不屬于泊松過程，目前有多種描述不完全維修的模型，在核電行業(yè)中應(yīng)用最多的為等效役齡回退因子模型所描述的不完全維修，本文中所提出的有關(guān)不完全維修仿真方法是在最小維修仿真方法的基礎(chǔ)上進(jìn)行改進(jìn)而構(gòu)造的，有關(guān)這方面的內(nèi)容在第3節(jié)會有詳細(xì)的闡述。

1 最小維修過程仿真

1．1 最小維修

最小維修的定義為：維修使得系統(tǒng)的故障率恢復(fù)到失效前的狀態(tài)。在維修決策優(yōu)化中，為簡化模型通?？紤]監(jiān)測維修為最小維修。在核電站的維修中最小維修的維修時間相對一般較短，在不考慮維修時間的情況下，最小維修條件下的壽命過程是一種非齊次泊松過程，其存在每次故障間隔時間相互不獨立的特點。

1．2 非齊次泊松過程

非齊次泊松過程是用來描述可修復(fù)系統(tǒng)的相鄰故障間隔呈某種趨勢的壽命過程，是齊次泊松過程的一種推廣。

非齊次泊松過程的數(shù)學(xué)的定義為：

｛N（t）；t≥0｝若滿足以下條件稱作是一個非齊次泊松過程。

（1）N（0）＝0；

（2）N（t）具有獨立增量；

（3）對任意的0≤t＜s，在（t，s］中的故障數(shù)N（s）－N（t）有參數(shù)的泊松分布：

非齊次泊松過程與齊次泊松過程的區(qū)別在于其相鄰故障間隔時間既不獨立亦不同分布，在某一區(qū)間上其故障參數(shù)還依賴于區(qū)間的起點。

泊松過程的強度函數(shù)決定了泊松過程，電廠的設(shè)備壽命分布多為威布爾分布，因此當(dāng)累計強度函數(shù)呈威布爾分布形式時將作為主要的研究對象，這類非齊次泊松過程稱為威布爾過程。

1．3 非齊次泊松過程仿真方法

對于最小維修過程的仿真主要即為對服從某一非齊次泊松過程的設(shè)備失效間隔時間進(jìn)行抽樣。目前業(yè)界存在幾種關(guān)于非齊次泊松過程的仿真方法，但沒有有關(guān)各種方法的適用范圍和使用特點的說明，事實上其中有很多種是難以在實際應(yīng)用中得以實施的。根據(jù)對目前主要應(yīng)用的仿真方法進(jìn)行比較，本人篩選出適合于核電設(shè)備的兩種：間隔時間產(chǎn)生方法、尺度變換方法。

方法一：間隔時間產(chǎn)生方法

間隔時間產(chǎn)生方法是根據(jù)設(shè)備的失效密度函數(shù)直接抽樣從而得到抽樣樣本，根據(jù)非齊次泊松過程的特點，設(shè)備第n次失效間隔時間的分布函數(shù)是與第n－1次失效的發(fā)生時間Sn－1相關(guān)，其關(guān)系式如定理1所述。

定理1［1］：設(shè)S0＝s0＝0，在Sn－1＝sn－1的條件下，tn＝Sn－Sn－1的條件密度函數(shù)為：

式中：Sn代表發(fā)生第n次失效的時間點，tn代表第n次失效發(fā)生的失效間隔時間，λ為非齊次泊松過程的強度函數(shù)。

根據(jù)定理1的敘述，在Sn－1＝sn－1條件下第n次的失效間隔時間的條件分布函數(shù)為：

采用間隔時時間產(chǎn)生方法其最大的限制在于需要對條件分布函數(shù)進(jìn)行抽樣，無論采用何種蒙卡的抽樣方法都需要分布形式不能過于復(fù)雜，因此間隔時間產(chǎn)生方法更多適用于失效率形式相對簡單的情況，例如威布爾分布（Λ（t）＝（t／α）β）等。

在設(shè)備壽命服從威布爾過程的條件下，其條件壽命分布函數(shù)的形式為：

式中：α是威布爾分布中的尺度參數(shù)，β為形狀參數(shù)。

該方法是通過依次產(chǎn)生設(shè)備失效間隔時間來對設(shè)備的壽命過程進(jìn)行仿真，在仿真過程中需要給出一個截止條件，可以選擇為觀測總時間。

方法二：尺度變換方法

尺度變換方法主要是根據(jù)非齊次泊松過程與齊次泊松過程之間的關(guān)系而構(gòu)造出的一種抽樣方法，其主要依據(jù)如定理2所述：

定理2［1］：設(shè)T1，…，Tn為強度函數(shù)為λ（t）的非齊次泊松過程頭n個事件發(fā)生時刻，記，則Z1，…，Zn獨立同分布，服從參數(shù)為1的指數(shù)分布。

根據(jù)定理2所構(gòu)造的尺度變換方法的仿真步驟為：

1．產(chǎn)生λ＝1的齊次泊松過程前n個事件的發(fā)生時刻Z1，…，Zn；

2．計算Si＝Λ－1（Zi），則S1，…，Sn滿足強度為λ（t）的非齊次泊松過程。

根據(jù)上面的步驟可以看出尺度變換方法是利用已經(jīng)產(chǎn)生的齊次泊松過程的失效時間通過一定的時間尺度變換映射，得出非齊次泊松過程的失效時間。在這一過程之中最主要的計算步驟在于計算累計強度函數(shù)的逆函數(shù)，這一點就可以降低該方法相對于時間間隔產(chǎn)生方法對失效率形式復(fù)雜度的要求。對于形式較為復(fù)雜的情況，例如常用于描述電器類設(shè)備失效的對數(shù)正態(tài)分布，很難構(gòu)造出適用于它的時間間隔產(chǎn)生方法的抽樣方法，而采用尺度變換方法，可以通過數(shù)值積分的方法計算出變量t與累計強度函數(shù)Λ（t）之間的映射關(guān)系，在計算Λ－1的時候，只需從之前得到的映射中尋找對應(yīng)值即可。

綜上，尺度變換方法主要適用于失效率分布形勢復(fù)雜的情況。在仿真的過程中需要給出一定的截止條件，由于在做尺度變換的前提是已經(jīng)仿真得出所有的Z1，…，Zn，因此截止條件可設(shè)置為總的失效次數(shù)n。

1．4 實例計算

為了驗證上面敘述的兩種仿真方法是否正確，對一假設(shè)參數(shù)的算例進(jìn)行抽樣，并對結(jié)果進(jìn)行擬合優(yōu)度檢驗。為方便起見選擇失效率為威布爾分布：λ（t）＝2×10－8t。適用于威布爾分布的檢驗方法可以采用Cramer Von－Mises檢驗。

1．4．1 Cramer Von－Mises擬合優(yōu)度檢驗

Cramer Von－Mises檢驗方法［2］為：記n次故障的故障時間分別為T1，T2，…，Tn；

對定數(shù)截尾情況，統(tǒng)計量C2M計算表達(dá)式：

若統(tǒng)計量計算值大于某一置信度下的統(tǒng)計量臨界值，則說明在該置信度下接受樣本服從該過程，反之則舍棄該樣本。

1．4．2 實例計算結(jié)果

分別采用時間間隔產(chǎn)生方法以及尺度變換方法對設(shè)備最小維修條件下的5次壽命過程進(jìn)行了仿真計算，得到每一次壽命過程中的設(shè)備逐次失效時間點Ti的結(jié)果如表1所示。

表1 最小維修抽樣結(jié)果Table 1 Sampled data of least maintenance

壽命過程1壽命過程2壽命過程3壽命過程4壽命過程5 9 947．92 7 671．81 16 130．96 12 770．2 9 597．52 T218 193．34 9 815．34 19 353．89 13 047．81 9 714．52 T322 337．03 13 172．31 19 371．01 14 229．92 19 895．92 T418 722．35 21 422．28 14 990．78 20 079．97 T525 777．95 23 853．56 19 294．04 20 080．48 T628 050．23 29 158．68 20 633．39 T721 010．94 T821 070．96 T924 942．03 T1025 388．5 T1 T11 25 409．91

采用Cramer Von－Mises檢驗：

通過擬合優(yōu)度檢驗可以在較大的置信度下接受抽樣樣本服從威布爾分布，因此兩種方法得到的結(jié)果是可信的，由此證明了兩種抽樣方法的可行性。

2 不完全維修過程仿真

不完全維修是最重要也是最難研究的一類維修過程，在電站的管理中，考慮設(shè)備的預(yù)防性維修為不完全維修。不完全維修的難點在于設(shè)備在該條件下維修后狀態(tài)的無法確定，設(shè)備的壽命過程無法繼續(xù)采用非齊次泊松過程進(jìn)行分析。因此在解決不完全維修過程仿真方法的過程中，建立設(shè)備壽命過程模型是其中一個關(guān)鍵的問題。

在目前研究中，Pham和Wang將描述不完全維修的各種方法歸納為七種類型［3］：（p，q）法則、（p（t），q（t））法則、復(fù)合（p，q）法則、（α，β）法則、沖擊模型、等效役齡、回退因子。在諸多描述方法中等效役齡與回退因子具有較強的靈活性和實用性，在核電站的數(shù)據(jù)分析中最多考慮。本章中將主要根據(jù)等效役齡與回退因子建立壽命過程模型并提出仿真方法。

2．1 等效役齡與回退因子

在Chan和Shaw［4］的研究中認(rèn)為預(yù)防性維修過程改善了設(shè)備的可靠性，設(shè)備可靠性的改善體現(xiàn)在失效率與可靠度的減少，其減少的程度用回退因子η來描述，而回退因子通常是由設(shè)備的役齡與維修的次數(shù)等參數(shù)決定的。根據(jù)之前的研究結(jié)果，η可以寫成預(yù)防性維修成本與維修次數(shù)的函數(shù)：

式中：a，b是回退因子的計算參數(shù)；Cpmi是第i次預(yù)防性維修的總成本，它的具體數(shù)值是與預(yù)防性維修固定成本Cpmf、預(yù)防性維修變動成本Cpmv以及設(shè)備役齡Ti均存在關(guān)系，計算關(guān)系式可以寫為：

因此回退因子的表達(dá)式可以寫為：

2．2 不完全維修過程仿真方法

目前還沒有任何有關(guān)不完全維修仿真方法的研究，通過分析引入回退因子后設(shè)備壽命過程與非齊次泊松過程之間的關(guān)系，在以非齊次泊松過程的仿真方法為基礎(chǔ)之上提出了不完全維修過程的仿真方法，填補了這方面的空白，為應(yīng)用仿真方法求解維修策略提供了支持。

本文中提出的不完全維修過程仿真方法是在最小維修的兩種仿真方法：間隔時間產(chǎn)生方法與尺度變換方法基礎(chǔ)之上進(jìn)行針對于回退因子的修改構(gòu)造得出的。

方法一：間隔時間產(chǎn)生方法

若設(shè)備第n－1次維修后認(rèn)為達(dá)到等效役齡Yn－1，則設(shè)備的第n次失效間隔時間分布函數(shù)可以寫為：

因此在考慮設(shè)備的等效役齡后，其失效間隔時間分布密度函數(shù)只需要將非齊次泊松過程的壽命分布函數(shù)中的第n－1次失效發(fā)生時間sn－1替換為第n次失效發(fā)生時設(shè)備的等效役齡Yn－1。

同樣，這種方法只適用于失效率分布形式相對簡單的情況，例如對于失效率服從威布爾分布情況，條件密度函數(shù)服從：

方法二：時間尺度變換方法

根據(jù)定理3的描述，對于服從非齊次泊松過程的前n個失效時刻S1，…，Sn，對其進(jìn)行尺度變換后的Z1，…，Zn服從強度為1的齊次泊松過程。在存在回退因子的情況下，顯然不能簡單地將Zn進(jìn)行尺度變換的逆變換，這是由于第n－1次維修后設(shè)備的等效役齡不同于維修前設(shè)備的等效役齡。因此需要對n－1次維修前的役齡進(jìn)行一定的變換使其與維修后的役齡相協(xié)調(diào)才可以對Z1，…，Zn進(jìn)行時間尺度的逆變換。

這種方法的本質(zhì)在于構(gòu)造出符合于存在回退因子情況下經(jīng)過時間尺度變換后服從于齊次泊松過程的逐次失效發(fā)生時間點Z＊1，…，Z＊n，再對其進(jìn)行尺度逆變換從而得到需要的結(jié)果。

其具體的步驟為：

1）對Z1進(jìn)行尺度變換產(chǎn)生第一次維修前的設(shè)備壽命；

2）設(shè)備壽命與回退因子η相乘產(chǎn)生第一次維修后的等效役齡；

5）以S2作為步驟二的輸入，重復(fù)上述步驟，直到產(chǎn)生所有的樣本S1，…，Sn。

其流程如圖1所示。

圖1 尺度變換方法流程圖Fig．1 Flowchart of time scale conversion method

實例計算：

為了驗證本文提出的兩種不完全維修仿真方法是否正確，需要對實例進(jìn)行計算，并對結(jié)果進(jìn)行驗證。需要說明的是目前并沒有針對役齡回退因子模型適合的檢驗方法，為此在實例計算中采用與最小維修一節(jié)相同的算例，并對兩者之間的結(jié)果進(jìn)行比較，從而給出分析與結(jié)論。采用的抽樣分布為威布爾分布：λ（t）＝2×10－8t；役齡回退因子與設(shè)備役齡之間的關(guān)系設(shè)定為η＝（0．1＋1×10－5T）2，對該條件下設(shè)備壽命過程進(jìn)行仿真，逐次的失效時點抽樣樣本如表2所示。

表2 不完全維修抽樣結(jié)果Table 2 Sampled data of incomplete maintenance

將上表中的結(jié)果與表1中的結(jié)果進(jìn)行比較，由于采用的隨機數(shù)完全一致，因此兩種維修模型設(shè)備出現(xiàn)故障的時間應(yīng)處在相同的分布密度函數(shù)區(qū)間內(nèi)。通過比較可以看出對于每一次壽命過程除首次故障時間外設(shè)備的壽命均有所增加，出現(xiàn)這種情況恰是由于引入了役齡回退因子，改善了維修后的設(shè)備可靠性，增加了設(shè)備的維修間隔時間。

3 結(jié)論

在核電站的維修策略優(yōu)化過程中，預(yù)防性維修是計算的主要對象。由于設(shè)備的失效存在概率性的特點，因此選用蒙特卡羅方法對存在維修的設(shè)備壽命過程進(jìn)行模擬所得到的結(jié)果相對數(shù)值方法而言存在應(yīng)用簡便、易于擴展的特點。

在本文中分別針對最小維修與不完全維修提出了設(shè)備壽命過程的抽樣方法。最小維修模型即非齊次泊松過程，在該模型下本文總結(jié)了兩種設(shè)備壽命過程的抽樣方法：間隔時間產(chǎn)生法、尺度變換方法。通過對兩種方法的使用特點進(jìn)行分析提出間隔時間產(chǎn)生方法局限于失效率形式相對簡單的情況，尺度變換方法步驟相對復(fù)雜但可以處理失效率復(fù)雜的情況。通過實例計算，對設(shè)備的壽命過程進(jìn)行仿真并對設(shè)備逐次的失效時間點進(jìn)行擬合優(yōu)度檢驗，可以證明兩種抽樣方法的正確性。不完全維修是一種復(fù)雜的過程，它不屬于非齊次泊松過程但與其存在一定的聯(lián)系。在通過對最小維修模型的兩種抽樣方法進(jìn)行修改后本文構(gòu)造出適用于等效役齡回退因子描述的不完全維修模型下設(shè)備壽命過程的仿真方法。通過分析實例計算得到的數(shù)據(jù)可以證明兩種方法的正確性不完全維修過程仿真方法的提出完善了核電站設(shè)備壽命過程的仿真方法，為應(yīng)用蒙卡仿真方法計算維修策略提供了支持。

［1］ 曹晉華，程侃．可靠性數(shù)學(xué)引論［M］．北京：高等教育出版社，2006．

［2］ Darling D A．The Kolmogorov－Smirnov，Cramer－von Mises Tests，Ann［J］．Math．Statist，1957，284，28（4）：823．

［3］ Chan J K，Shaw L．Modeling repairable systems with failure rates that depend on age and maintenance［J］．IEEE Trans．on Reliability，1993，42：566－570．

［4］ H Wang，H Pham．Some maintenance models and availability with imperfect maintenance in production systems［J］．Annals of Operations Research，1999，91，（0）．

［5］ 徐文靜．不完全維修條件下的可用度與維修策略分析［D］．湖南：國防科技大學(xué)，2008．

［6］ 程志君，高大化．不完全維修條件下的視情維修優(yōu)化模型［A］．系統(tǒng)工程與電子技術(shù)，2006，28（7）：1106－1108．

［7］ 徐文靜，郭波．基于回復(fù)因子的不完全維修模型分析［A］．科技創(chuàng)新導(dǎo)報，2008（24）：74－75．

［8］ 張寧平，陳躍良．基于可靠度經(jīng)濟性約束的不完全維修周期研究［A］．新技術(shù)新工藝，2010，3：17－19．

［9］ Ebastian Martorell，Ana Sanchez，Vicente Serradell．Age－dependent reliability model considering effects of maintenance and working conditions［J］．Reliability Engineering and System Safety，1999，64，64（1）：19－31．

［10］ Pham H，Wang H．Optimal maintenance policies for several imperfect maintenance models［J］．International Journal of Systems Science，1996，12，12（2）：87－101．