基于蒙特卡羅模擬的航空發(fā)動機故障風險預測

趙洪利，劉宇文

(1.中國民航大學 航空工程學院，天津300300;2.中國民航大學 中歐航空工程師學院，天津300300)

航空發(fā)動機的風險水平指在規(guī)定的時間內(nèi)、正常飛行的條件下，航空發(fā)動機出現(xiàn)故障的概率或次數(shù)［1］.故障風險預測方法能夠給出航空發(fā)動機在未來一段時間周期內(nèi)可能發(fā)生故障的次數(shù).航空運營企業(yè)為了保證飛機運行的可靠性，需要計劃好發(fā)動機維修，并儲備相應數(shù)量的需備件，這就需要對發(fā)動機的故障風險水平做出預測，估計發(fā)動機在未來一段時間內(nèi)各部件及發(fā)動機系統(tǒng)出現(xiàn)故障的次數(shù)，從而為生產(chǎn)或運營計劃做出安排.

航空發(fā)動機部件多，結(jié)構(gòu)復雜，不同部件都有各自的故障概率分布.經(jīng)過多年發(fā)展，航空發(fā)動機部件故障概率的分析方法有很多［2-3］.其中，經(jīng)典的解析方法雖然能夠給出單個部件統(tǒng)計特性的精確表達式，但是用于確定多部件組成的復雜系統(tǒng)的故障風險等可靠性參數(shù)是十分困難的.而其他系統(tǒng)分析方法，如故障樹分析法、馬爾可夫狀態(tài)轉(zhuǎn)移法等，其問題求解的規(guī)模往往會隨著部件數(shù)目的增多而呈指數(shù)量級增大［4-5］，且對所求解的系統(tǒng)往往有一定限制(如部件故障概率分布形式等).蒙特卡羅方法作為一種隨機模擬方法，通過對模型或系統(tǒng)的觀察或抽樣試驗來計算所求參數(shù)的統(tǒng)計特征，能夠很好地避免以上問題.

1 蒙特卡羅方法

蒙特卡羅方法的主要步驟［6］包括:針對具體問題構(gòu)造相應的概率過程;實現(xiàn)已知概率分布的隨機抽樣;利用隨機抽樣的結(jié)果計算具體的估計量.

1.1 構(gòu)造概率過程

構(gòu)造概率過程即找到具體問題所服從的概率分布，如二項分布、指數(shù)分布、正態(tài)分布等.準確描述問題的概率過程是蒙特卡羅求解方法的關(guān)鍵，整個隨機模擬過程都是建立在已知的概率分布之上.如何確定問題所服從的概率分布，這需要根據(jù)已有的數(shù)據(jù)并結(jié)合相應的數(shù)據(jù)處理方法來確定.

1.2 隨機抽樣與統(tǒng)計

構(gòu)造完概率過程后，就知道了系統(tǒng)各個部分所服從的概率分布，如一輛汽車中儀表盤等電子元器件壽命可能服從指數(shù)分布，發(fā)動機故障時間可能服從威布爾分布等.在已知概率分布的基礎(chǔ)上，產(chǎn)生滿足相應分布的隨機變量即為隨機抽樣過程.隨機抽樣的核心在于隨機數(shù)的生成，常用的方法有工程上的物理方法和數(shù)學公式遞推法.工程上的物理方法無法重復便捷地生成隨機數(shù)，且價格昂貴.數(shù)學遞推公式法利用計算機就可重復產(chǎn)生大量的隨機數(shù)，雖然無法做到真正意義上的隨機性，但只要滿足相應的隨機數(shù)檢驗，即可滿足大部分要求［7］.

蒙特卡羅模擬方法能否得到好的結(jié)果關(guān)鍵在于整個隨機過程的構(gòu)造和計算過程中隨機數(shù)的性質(zhì).在隨機抽樣的基礎(chǔ)上，根據(jù)所求解的問題，對隨機抽樣過程中產(chǎn)生的數(shù)據(jù)進行記錄、統(tǒng)計，進而確定問題的估計量，得到問題的解.

2 故障風險預測算法

2.1 故障概率模型

在航空發(fā)動機故障數(shù)據(jù)處理中，威布爾分布是最常見的，且與數(shù)據(jù)符合程度相對較好的一種分布［8］.航空發(fā)動機是一種高可靠性的產(chǎn)品，通常只有小樣本的故障數(shù)據(jù).威布爾分析在處理小樣本數(shù)據(jù)時，與其他分析方法相比通常都具有較好的效果［9-10］.因此，在航空發(fā)動機故障風險預測中，本文采用了威布爾分布模型.

威布爾累積概率分布函數(shù)為

式中，F(xiàn)(t;β，η)為累積概率，t為故障時間，β 為形狀參數(shù)，η為尺度參數(shù);t0為起始點.

對于小樣本航空發(fā)動機故障數(shù)據(jù)，為了得到合適的尺度參數(shù)和形狀參數(shù)，采用中位秩回歸參數(shù)估計法［11］.

1)將故障時間數(shù)據(jù)T=［T1T2… Tn］T從小到大排列，并利用式(2)計算各個數(shù)據(jù)的中位秩，當故障時間數(shù)據(jù)中存在刪失數(shù)據(jù)時，需利用式(3)調(diào)整數(shù)據(jù)排序值.

式中，i為數(shù)據(jù)序號;Ri為第i個數(shù)據(jù)的中位秩;N為數(shù)據(jù)總數(shù).

式中，i′為調(diào)整后的排序值;j為前一個數(shù)據(jù)調(diào)整后的排序值;p為當前刪失數(shù)據(jù)之后的數(shù)據(jù)個數(shù).

2)計算故障時間的自然對數(shù)，記為Y=ln T，并令

式中，I為單位矩陣;RY為Y的中位秩.

式中，xi，yi為向量 X，Y 的分量.

2.2 隨機變量的抽樣方法

滿足威布爾分布的隨機數(shù)生成算法有很多，反變換法是其中一種便捷、有效的方法［12-13］.設需產(chǎn)生分布函數(shù)F(x)的連續(xù)隨機數(shù)x，若已有［0，1］區(qū)間均勻隨機數(shù) r，則產(chǎn)生x的反變換公式［14］為

則威布爾的反變換公式為

式中，T為隨機故障時間.

反變換法的關(guān)鍵在于［0，1］區(qū)間上高品質(zhì)的均勻隨機數(shù).采用目前廣泛使用的乘同余組合發(fā)生器來產(chǎn)生［0，1］區(qū)間上的均勻隨機數(shù).遞推公式［15］為

式中，i=0，1，…;j=1，2，…，m;ri為［0，1］區(qū)間上隨機數(shù);mod為求余運算符;m為組合數(shù).

選取組合數(shù)m=2，且令

2.3 蒙特卡羅模擬流程

設發(fā)動機有n種故障模式，記為F1，F(xiàn)2，…，F(xiàn)n.這n種故障模式相互獨立，均服從威布爾分布，參數(shù)記為(η1，β1)，(η2，β2)，…，(ηn，βn)，故障時間抽樣為

式中，rk∈R，R為乘同余組合發(fā)生器生成的隨機數(shù)序列.

假定發(fā)動機定期檢查時間為T，故障修復后，發(fā)動機使用時間重新計，累積到時間T時，再做下次定檢，則蒙特卡羅模擬運行流程如圖1所示.

圖1 蒙特卡羅模擬流程Fig.1 Monte Carlo simulation process

由圖1可知，模擬的基本步驟如下:

1)輸入機隊的原始數(shù)據(jù)，包括機隊總數(shù)，使用率，機隊年齡分布;

2)利用乘同余組合發(fā)生器構(gòu)建［0，1］區(qū)間均勻隨機數(shù)表;

3)針對每臺發(fā)動機，從隨機數(shù)表中順序選取隨機數(shù)，利用式(13)計算各個故障模式的隨機故障時間 F1，F(xiàn)2，…，F(xiàn)n，首次故障時間需大于發(fā)動機的已安全運行時間，否則需重新產(chǎn)生一組隨機故障時間，直到均大于發(fā)動機已安全運行時間;

4)判斷 min{F1，F(xiàn)2，…，F(xiàn)n}是否小于定期檢查時間T，若是，則記錄該故障Fk，該故障模式發(fā)生次數(shù)累加1;

5)針對k故障模式，重新抽樣故障時間Fk;

6)重復步驟4)和步驟5)，直到min{F1，F(xiàn)2，…，F(xiàn)n}超過發(fā)動機定期檢查時間T，則一次模擬結(jié)束;

7)遍歷機隊內(nèi)所有發(fā)動機，每臺發(fā)動機均進行N次模擬，并記錄結(jié)果;

8)計算發(fā)動機各個故障模式發(fā)生的概率;

9)結(jié)合機隊原始數(shù)據(jù)，計算機隊未來一段時間內(nèi)發(fā)動機各個故障模式可能發(fā)生的次數(shù).

3 算例及結(jié)果分析

3.1 算例描述

為了驗證該模型的準確性與適用性，本文采用了某發(fā)動機公司手冊中的數(shù)據(jù)作為輸入條件，利用上述故障風險預測算法進行模擬仿真，并將結(jié)果與發(fā)動機公司給出的軟件仿真結(jié)果進行分析對比.手冊中提供了某噴氣式發(fā)動機4種相互獨立的故障模式:F1表示過熱，F(xiàn)2表示葉片裂紋，F(xiàn)3表示油管裂紋，F(xiàn)4表示燃燒室裂紋.4種故障模式均服從不同參數(shù)的威布爾分布，具體參數(shù)見表1.同時，發(fā)動機的定期檢查時間為1 000 h，整個機隊發(fā)動機運行時間分布如圖2所示.假定發(fā)動機使用率為25 h/月.

表1 4種故障模式參數(shù)列表Table1 Parameters list of 4 failure modes

3.2 模擬流程

針對初始時間為0的發(fā)動機，詳細闡述蒙特卡羅風險預測方法，模擬流程見圖3.

1)為4種模式生成隨機故障時間.利用隨機數(shù)表，順序選取前4個［0，1］區(qū)間上的隨機數(shù):0.007，0.028，0.517，0.603.根據(jù)式(13)，則有:

圖2 機隊發(fā)動機運行時間分布Fig.2 Operation time distribution of engine fleet

圖3 模擬流程Fig.3 Simulation process

2)4種故障模式中，最小的故障時間為951 h，并未達到定期檢查時間1000 h.

3)最先發(fā)生的為F1故障，記錄該故障發(fā)生時間為951 h，相當于在未來第38個月發(fā)生故障.從隨機數(shù)表中選取下一個隨機數(shù)0.442，為F1故障重新生成隨機故障時間，F(xiàn)1=8827 h.此時，min{F1，F(xiàn)2，F(xiàn)3，F(xiàn)4}已超過定期檢查時間 1 000 h.至此，一次模擬結(jié)束.

針對每臺發(fā)動機重復模擬N次，并遍歷所有發(fā)動機.根據(jù)大數(shù)定理，則發(fā)動機故障發(fā)生概率近似為

式中，i=1，2，3，4 為故障模式序號;j=1，2，…，12為發(fā)生故障的月份;Ni，j為第i種故障在j月份發(fā)生的次數(shù);Pi，j為第i種故障模式在j月份發(fā)生的概率;Na為額外抽樣次數(shù);Ne為發(fā)動機數(shù)量.

則發(fā)動機故障風險預測結(jié)果為

式中，F(xiàn)i，j為第 i種故障在 j月份發(fā)生的預測次數(shù).

發(fā)動機一年內(nèi)的故障風險預測結(jié)果見表2，某發(fā)動機公司提供的內(nèi)部軟件預測結(jié)果見表3，圖4給出了12月份發(fā)動機故障次數(shù)的收斂結(jié)果.

表2 發(fā)動機累積故障次數(shù)預測結(jié)果Table2 Forecasting results of cumulative failure number of engine

表3 某發(fā)動機公司預測結(jié)果Table3 Forecasting results provided by an engine company

圖4 12月份發(fā)動機故障次數(shù)收斂結(jié)果Fig.4 Convergence results of engine failure number on December

3.3 結(jié)果分析

在模擬過程中，為了得到穩(wěn)定的結(jié)果，需要實時觀察模擬輸出的數(shù)值，以便確定輸出結(jié)果是否穩(wěn)定.以12月份發(fā)動機的故障次數(shù)預測結(jié)果為例，每10000次輸出一次模擬數(shù)值，并做100次抽樣，則總模擬次數(shù)為106.由圖4可知，最終結(jié)果已收斂.一般來說，模擬次數(shù)越大，預測結(jié)果越逼近真實情況.在試驗過程中，通過對比總模擬次數(shù)為106，107和108時的收斂圖形，發(fā)現(xiàn)結(jié)果差別不大.因此，總模擬次數(shù)設定為106已滿足收斂性要求.同時，由表2和表3可以看出，機隊內(nèi)發(fā)動機一年內(nèi)各部件累積故障次數(shù)的預測結(jié)果與發(fā)動機公司提供的數(shù)據(jù)相差不大，故障總次數(shù)僅相差0.97.預測結(jié)果準確，驗證了該算法的有效性與適用性.

4 結(jié)論

1)在確定了失效分布規(guī)律后，利用蒙特卡羅模擬算法進行仿真，能夠比較準確地估算出整個機隊發(fā)動機在未來不同時間段內(nèi)的故障風險水平，從而為發(fā)動機維修管理提供可靠性指標參考;

2)蒙特卡羅模擬是多故障模式下風險預測的一個有效方法，并且在多故障模式下，該算法不但可以預測機隊整體風險水平，而且還可確定何種故障模式最易發(fā)生，從而可以有針對性地對該故障制定相應的改進措施，降低機隊的故障風險;

3)當某個易發(fā)故障得到改進后，可重新進行仿真模擬來找到下一個易發(fā)故障，不斷迭代，從而實現(xiàn)對機隊內(nèi)發(fā)動機的動態(tài)管理.

致謝 中國民航大學郭慶副教授提供了某發(fā)動機廠商的技術(shù)資料，并為文章寫作提出了建設性意見，在此表示感謝!

References)

［1］ 趙宇.可靠性數(shù)據(jù)分析［M］.北京:國防工業(yè)出版社，2011:14-20.Zhao Y.Data analysis of reliability［M］.Beijing:National Defense Industry Press，2011:14-20(in Chinese).

［2］ 楊宇航，馮允成.基于仿真的復雜系統(tǒng)可靠性、可用性和MTBF評估文獻綜述［J］.系統(tǒng)工程理論與實踐，2003，42(2):80-85.Yang Y H，F(xiàn)eng Y C.Survey of reliability and availability evaluation of complex system using Monte Carlo techniques［J］.Systems Engineering-Theory ＆ Practice，2003，42(2):80-85(in Chinese).

［3］ 劉曉平，唐益明，鄭利平.復雜系統(tǒng)與復雜系統(tǒng)仿真研究綜述［J］.系統(tǒng)仿真學報，2008，20(23):6303-6315.Liu X P，Tang Y M，Zheng L P.Survey of complex system and complex system simulation［J］.Journal of System Simulation，2008，20(23):6303-6315(in Chinese).

［4］ 邵偉.蒙特卡洛方法及在一些統(tǒng)計模型中的應用［D］.濟南:山東大學，2012.Shao W.Monte Carlo methods theory and their applications in some statistical model［D］.Jinan:Shandong University，2012(in Chinese).

［5］ 袁明偉.復雜系統(tǒng)可靠性分析［D］.馬鞍山:安徽工業(yè)大學，2013.Yuan M W.Reliability analysis of complex system［D］.Maanshan:Anhui University of Technology，2013(in Chinese).

［6］ Dickman B H，Gilman M J.Technical note Monte Carlo optimization［J］.Journal of Optimization Theory and Applications，1989，60(1):149-157.

［7］ 周文彬.組合式偽隨機數(shù)發(fā)生器的研究與設計［D］.哈爾濱:哈爾濱工程大學，2013.Zhou W B.Design and research of the combined pseudo-random number generator［D］.Harbin:Harbin Engineering University，2013(in Chinese).

［8］ Baumshtein M V，Prokopenko A V，Ezhov V N.Probabilistic prediction of the fatigue life of gas-turbine engine compressor blades under two-level programmed loading［J］.Strength of Ma-terials，1985，17(5):587-592.

［9］ Nee A Y C，Song B，Ong S K.Re-engineering manufacturing for sustainability［M］.Singapore:Springer Singapore，2013:699-703.

［10］ Saralees N，Samuel K.Strength modeling using Weibull distributions［J］.Journal of Mechanical Science and Technology，2008，22(7):1247-1254.

［11］ 麻曉敏，張士杰，胡麗琴，等.可靠性數(shù)據(jù)威布爾分析中秩評定算法改進研究［J］.核科學與工程，2007，27(2):152-155.Ma X M，Zhang S J，Hu L Q，et al.An improved rank assessment method for weibull analysis of reliability data［J］.Nuclear Science and Engineering，2007，27(2):152-155(in Chinese).

［12］ Yadolah D.The concise encyclopedia of statistics［M］.New York:Springer New York，2008:446-447.

［13］ Johnson G E.Constructions of particular random processes［J］.Proceedings of the IEEE，1994，82(2):270-285.

［14］ 金暢.蒙特卡羅方法中隨機數(shù)發(fā)生器和隨機抽樣方法的研究［D］.大連:大連理工大學，2005.Jin C.Study on random number generator and random sampling in Monte Carlo method［D］.Dalian:Dalian University of Technology，2005(in Chinese).

［15］ 楊自強，魏公毅.常見隨機數(shù)發(fā)生器的缺陷及組合隨機數(shù)發(fā)生器的理論與實踐［J］.數(shù)理統(tǒng)計與管理，2001，20(1):45-51.Yang Z Q，Wei G Y.Drawbacks in classical random number generators-theory and practice of combined generator［J］.Journal of Applied Statistics and Management，2001，20(1):45-51(in Chinese).