動態故障樹頂事件發生概率的蒙特卡洛近似算法＊

楊恒占，張 可，錢富才

（1．西安工業大學 電子信息工程學院，西安710021；2．西安理工大學 自動化與信息工程學院，西安710048）

動態故障樹（Dynamic Fault Tree，DFT）是采用動態邏輯門進行實際物理系統故障診斷的一種方法．自1992年動態故障樹模型首次提出起［1］，已在不同領域得到深入的研究．由于動態故障樹建模方便、高效，其實際應用也越來越廣泛．例如容錯計算機、航空電子設備系統、FTPP（容錯并行處理器）配置、航天系統等系統的可靠性分析中，大量應用動態故障樹計算故障事件發生概率［2-3，10-11］．目前，在動態故障樹算法理論研究方面，主要形成了三種算法思路［4］：基于Markov狀態轉移過程的分析算法、基于貝葉斯網絡（Bayesian Network，BN）的分析算法和基于梯形公式的近似分析算法［5-7］．這些方法在不同領域都得到了一定的應用，但普遍存在計算復雜、計算效率低問題，阻礙了在線分析應用，尤其在故障樹的規模較大時更為明顯．

在動態故障樹中，故障事件位于故障樹頂部，引起它發生的隨機因素很多，而其發生概率就是眾多隨機變量的概率密度函數在特定區域的積分．這是一個多重積分的計算問題，除非該積分能轉化為多個單重積分，否則這類問題的傳統算法都會遇到計算爆炸問題，而且隨積分維數的增加計算量呈指數增長．這也是傳統算法效率低下的根本原因．因此，迫切需要尋找一種簡單、高效的求解動態故障樹頂事件故障概率的算法．

本文在動態故障樹分析過程中應用Monte Carlo方法，將多重積分等價為具有特定分布隨機函數的期望值，設計隨機實驗以獲得隨機被積函數的樣本．

1 動態邏輯模塊頂事件發生概率的計算公式

動態故障樹由多種動態邏輯門根據實際物理關系組合而成．動態邏輯門包括順序相關門、優先與門、冷備件門、熱備件門和功能相關門．已知一系列隨機變量的分布就可以求得這些隨機變量函數的分布，因此通過分析動態故障樹子樹的失效分布就能夠得到動態故障樹頂事件的失效分布．各種動態邏輯模塊失效分布的運算機理為

1．1 順序相關門

順序相關門的邏輯功能是當底事件以從左到右的順序依次發生時頂事件才發生．假設底事件xi（i＝1，2，…，m-1，m）的發生時間Ti的概率分布函數依次為Gi（xi），那么頂事件的發生時間T的概率分布函數GS（t）為

1．2 優先與門

優先與門是順序相關門輸入為兩個時的特例，在邏輯上相當于“與門”，但是要求輸入事件必須以指定的順序發生．假設底事件x1和x2的發生時間T1和T2的概率分布函數分別為G1（x1）和G2（x2），則頂事件的發生時間 的概率分布函數GY（t）為

1．3 冷備件門

冷備件門有一個主件和一個備件，且二者都是基本事件．主件一開始就進入工作狀態；備件則是指一開始不工作只是作為主件的替代備件，其失效率為零．假設主件x1和備件x2的發生時間T1和T2的分布函數分別為G1（x1）和G2（x2），則頂事件的發生時間的分布函數GL（t）為

1．4 熱備件門

熱備件門有一個主件和一個備件，且二者都是基本事件．主件與備件同時處于工作狀態，負荷分擔，主件與備件的失效分布函數基本相同；當主件和備件均失效時，頂事件發生．一般應用在可靠性要求較高系統中．假設主件x1和備件x2的發生時間T1和T2的分布函數分別為G1（x1）和G2（x2），則頂事件的發生時間T的分布函數GR（t）為

1．5 功能相關門

功能相關門的邏輯功能為當觸發事件發生時就會引起基本事件發生，則頂事件也會發生；或者基本事件單獨發生時，頂事件也會發生．假設觸發事件x1和底事件x2的發生時間T1和T2的分布函數分別為G1（X1）和G2（X2），則頂事件的發生時間T的分布GG（t）為

綜上分析，多重積分的計算成為動態故障樹分析方法的核心問題，而這些積分的計算量過大問題嚴重制約了動態故障樹方法的推廣應用．

2 Monte Carlo方法

Monte Carlo方法是指任何利用隨機數列來做隨機模擬的數值方法，該方法一般通過構造符合一定規則的隨機數來解決數學上的各種問題．對于計算過于復雜而難以得到解析解或者根本沒有解析解的問題，Monte Carlo方法是一種有效的求出數值解的方法，Monte Carlo方法在數學中最常見的應用就是Monte Carlo積分．許多物理問題都會涉及到多重積分的計算問題，對于這類問題，計算量隨維數的增加呈指數增長，傳統的數值算法難以求解，Monte Carlo積分方法能很好的求解多重積分問題，實現多重積分的快速計算，因為該方法的計算不依賴于維數．

Monte Carlo方法的解題過程可以歸結為三個主要步驟：構造或描述概率過程；實現從已知概率分布抽樣；建立各種估計量．Monte Carlo求解多重積分發展較成熟的方法主要有兩種［8-9］：隨機投點法和均值估計法．因為隨機投點法的誤差較大且做法繁瑣，通常不被采用．均值估計法比隨機投點法更好操作，精度更高（若選擇的密度函數合適）．文中采用均值估計法求解動態故障樹中的多重積分，最終得到頂事件發生概率．均值估計法原理．

考慮k重積分

其中Ω為k維積分區域．

選取Ω上的一個概率密度函數g（x1，x2，…，xk-1，xk），使它滿足當（x1，x2，…，xk-1，xk）∈Ω時g（x1，x2，…，xk）≠0．

其中m是隨機點落在Ω區域中的投點數．

Monte Carlo積分本質在于將積分平均轉化為統計平均，因此大大簡化了積分計算過程．本文結合某衛星太陽翼驅動機構的步進電機實例，將此方法應用于實際系統中．

3 衛星太陽翼驅動機構步進電機實例分析

該步進電機由定子、轉子、軸承、定子繞組等組成，而其易損零部件主要為軸承與定子繞組．為了提高系統的可靠性，航天器機構中的步進電機通常采用雙繞組作為備份手段．其中軸承模塊包括軸承1和軸承2兩個，它們與定子繞組之間為“或門”關系，即其中一個故障，系統失效．定子繞組模塊由定子繞組主件和定子繞組備件組成，它們之間是冷備件關系，即主件失效，則啟用備件，當備件再失效，則宣布系統失效．

由以上分析可得“步進電機故障”的動態故障樹模型如圖1所示．

圖1 “步進電機故障”DFT模型Fig．1 DFT model for“Stepping motor fault”

圖1中A，C，X1和X2為模型的底事件，分別代表軸承1，軸承2，定子繞組主件和定子繞組備件；B為模型的中間事件，代表定子繞組；T為模型的輸出事件，代表步進電機失效．

由以上DFT模型可知，頂事件T發生概率為

假設軸承1和軸承2的壽命均服從指數分布，其失效率為1．0×10-6h-1，其失效率函數分別為

假設定子繞組的主件與備件的壽命服從w（100 000，3）分布，其失效概率函數分別為

根據式（3）可得定子繞組的失效概率函數為

選取g（t1，t2）為積分區域Ω上的均勻分布，設Ω的體積為Vs，則

由式（6）和式（10）可轉換為

將定子繞組的失效概率函數以及軸承1和軸承2作為頂事件的輸入，即可求得“步進電機故障”動態故障樹模型頂事件在t時刻的發生概率．

采用Monte Carlo算法在不同時刻編程求解式（11），結合式（7）～（8）求得頂事件T的發生概率．

為了說明本文方法的優越性，針對“步進電機故障”動態故障樹模型與梯形公式法進行了仿真結果對比，得出頂事件不可靠度和計算效率的比較結果見表1．

可見，Monte Carlor方法可以應用于實際航天系統中重要的數值積分計算，實現航天系統的故障診斷．應用Monte Carlo方法計算動態故障樹頂事件故障概率，簡化了計算過程，顯著提高了計算效率，并且隨著積分維數的增加，Monte Carlo方法性能優勢更加明顯．

表1 頂事件不可靠度和仿真效率結果比較Tab．1 Comparison between the top event reliability and the computational efficiency

4 結 論

1）本文將Monte Carlo方法應用于動態故障樹頂事件發生概率的求解中，將頂事件概率轉換為不同的多重積分，把多重積分等價為具有特定分布隨機函數的期望值，該樣本均值就是多重積分的近似值，設計隨機實驗以獲得隨機被積函數的樣本，用Monte Carlo方法求出樣本均值，實現了動態故障樹的定量分析．該樣本均值就是多重積分的近似值．

2）通“衛星太陽翼驅動機構步進電機”實例分析與梯形公式法進行比較，本文方法與梯形公式法相比較在計算結果上基本一致，仿真效率比梯形公式法提高了50%，滿足動態故障樹頂事件發生概率的計算要求．

3）采用Monte Carlo方法進行動態故障樹的定量分析，不依賴于積分維數，尤其適用于復雜多維積分問題的求解．能夠有效解決動態故障樹方法固有的計算困難問題，有效避免了傳統故障樹分析中“組合爆炸”問題的發生，滿足航天系統在線分析的要求．

［1］ DUGAN J B，BAVUSO S J，Boyd M A．Dynamic Fault-Tree Models for Fault Tolerant［J］．IEEE Transactions on Reliability，1992，41（3）：363．

［2］ 范長征，蔣凡，曾凡平．綜合故障樹分析方法在容錯計算機系統中的應用［J］．計算機仿真，2006（4）：63．FAN Chang-zheng，JIANG Fan，ZENG Fan-ping．Aplication of Integrated Fault Tree Analysis to Fault Tree Analysis to Fault-tolerant Computer Systems［J］．The Computer Simulation，2006（4）：63．（in Chinese）

［3］ 邢曉辰，蔡遠文，程龍，等．基于故障樹的目標航天器在軌故障定位［J］．兵工自動化，2013（10）：71．XING Xiao-chen，CAI Yuan-wen，CHENG Long，et al．On-Orbit Fault Positioning of Tanget Space Craft Based on Fault Free［J］．These Automation，2013（10）：71．（in Chinese）

［4］ 段凌昊，郭愛民，潘勇．動態故障樹分析算法研究綜述［J］．電子產品可靠性與環境試驗，2013（4）：59．DUAN Ling-hao，GUO Ai-min，PAN Yong．Review on Dynamic Fault Tree Analysis［J］．Electronic Product Reliability and Environmental Testing，2013（4）：59．（in Chinese）

［5］ BOUDALI H，DUGAN J B．A Discrete-Time Bayesian Network Reliability Modeling and Analysis Framework［J］．Reli-ability Engineering and System Safety，2005，873：37．

［6］ AMARI S．A New Approach to Solve Dynamic Fault Tree［C］／／Proceeding Annual Reliability and Maintainability Symposium．Alexander：RAMS，2005：374．

［7］ GULATI R，DUGAN J B．A Modular Approach to Analyzing Static and Dynamic Fault Tree［C］／／Proceedings Annual Reliability and Maintainability Symposium．Fort Worth：RAMS，2009：57．

［8］ 馬海峰，劉宇熹．基于 Monte-Carlo方法計算定積分的算法研究［J］．計算機與現代化，2011（11）：15．MA Hai-feng，LIU Yu-xi．Research on Calculating Definite Integrals Method Based on Monte-Carlo［J］．Computer and Modern，2011 （11）：15．（in Chinese）

［9］ 方再根．計算機模擬和蒙特卡洛方法［M］．北京：北京工業學院出版社，1988．FANG Zai-gen．Computer Simulation and Monte Carlo Method［M］．Beijing：Beijing Industrial Institute Publishing House，1988．（in Chinese）

［10］ 梅啟智．系統可靠性工程基礎［M］．北京：科學出版社，1992．MEI Qi-zhi．System Reliability Engineering Foundation［M］．Beijing：Science Press，1992．（in Chinese）

［11］ 吳嘉寧，閻紹澤，謝里陽．基于RBD和FTA的航天器太陽翼可靠性分析［J］．清華大學學報：自然科學版，2012（1）：15．WU Jia-ning，YAN Shao-ze，XIE Li-yang．Reliability Analysis of Spaceraft Solar Arrays Based on RBD and FTA［J］．Journal of Tsinghua University：Science and Technology，2012（1）：15．（in Chinese）