基于最大熵和仿真的導(dǎo)彈裝備備件需求預(yù)測

趙建忠，徐廷學(xué)，李海軍，徐衡博

(1．海軍航空工程學(xué)院 兵器科學(xué)與技術(shù)系，山東 煙臺，264001； 2．中國人民解放軍91423部隊 裝備部，山東 萊陽 265200)

0 引言

備件是導(dǎo)彈裝備使用和維修的重要物質(zhì)基礎(chǔ)，其保障水平對于保持導(dǎo)彈裝備的戰(zhàn)備完好性有著直接的影響。準(zhǔn)確地預(yù)測備件的消耗數(shù)目才能在有限的經(jīng)費條件下，最大限度地滿足導(dǎo)彈裝備保障的備件需求。目前，大部分的備件預(yù)測方法或模型將備件視為在冷儲備條件，即不考慮備件在儲備期間發(fā)生故障的情況，而事實上導(dǎo)彈裝備經(jīng)常會出現(xiàn)零部件或備件在儲備期間損壞而導(dǎo)致短缺的問題，所以它們的預(yù)測精度往往不高。因此，建立間斷工作導(dǎo)彈裝備備件需求預(yù)測模型，更準(zhǔn)確地確定備件需求量成為亟待解決的問題。

在建立模型之前，首先需要確定導(dǎo)彈裝備的壽命分布。導(dǎo)彈裝備故障機(jī)理十分復(fù)雜，而且其備件消耗大都呈現(xiàn)“短周期、小樣本”的特點，故采用通過大量數(shù)據(jù)樣本進(jìn)行建模的方法計算導(dǎo)彈備件的壽命分布參數(shù)，效果往往不太理想。實踐表明，對于常見統(tǒng)計分布函數(shù)，壽命分布假定有時在0.05的水平被拒絕[1]。為了解決這一問題，本文將導(dǎo)彈壽命分布形式看作未知，提出用最大熵法確定壽命分布，并利用遺傳算法求解系統(tǒng)建模過程中的非線性規(guī)劃問題。

文獻(xiàn)[2]利用最大熵法來確定統(tǒng)計分布函數(shù)參數(shù)。由于傳統(tǒng)最大熵模型中的約束條件是基于經(jīng)典矩構(gòu)建的，因此，基于經(jīng)典矩的最大熵法能夠有效擬合大樣本數(shù)據(jù)的分布函數(shù)，然而導(dǎo)彈裝備故障的數(shù)據(jù)往往不能滿足大樣本條件[3]。本文結(jié)合最大熵原理與概率加權(quán)矩，提出了基于最大熵原理和概率加權(quán)矩的導(dǎo)彈裝備壽命分布確定方法。它直接從樣本信息出發(fā)，不需要對待估隨機(jī)變量的統(tǒng)計分布類型作任何假定，從而為建立統(tǒng)計分布函數(shù)提供了便利。

1 基于最大熵原理與概率加權(quán)矩的部件導(dǎo)彈單元壽命分布確定

為了克服基于不同樣本量估計的模型對于系統(tǒng)統(tǒng)計特性的描述偏差相差極大的缺陷,本文引入概率加權(quán)矩(probability weighted moments, PWMs)方法。根據(jù)最大熵原理[4]，對于定義域R上的連續(xù)隨機(jī)變量x(x的概率密度函數(shù)為f(x))，假設(shè)其累積分布函數(shù)F(x)=P(X≤x)及其逆函數(shù)x=x(F)存在，則其概率加權(quán)矩可以表示為

(1)

式中：變量i,n,k為實數(shù)。

目前常用的概率加權(quán)矩有[5]

(2)

(3)

并設(shè)αk與βn對應(yīng)的無偏估計為ak和bn。

由式(1)可得

(4)

在上述分析的基礎(chǔ)上，可得最大熵逆累積分布函數(shù)模型

(5)

求解得

(6)

式中：N為樣本矩的最高階數(shù)；λi,i=0,1,…,N為與第i階矩約束相對應(yīng)的拉格朗日算子。

如果式(5)中約束等于各階經(jīng)典矩，式(6)概率密度函數(shù)的約束等于各階概率加權(quán)矩，則式(6)為逆累積分布函數(shù)。

在上述分析的基礎(chǔ)上，bn可表示為[7-8]

(7)

式(7)是關(guān)于λi(i=0,1,…,N)的N個方程。

為便于數(shù)值求解,將式(7)改寫為如下形式：

(8)

確定如下優(yōu)化問題:

(9)

通過遺傳算法求解式(9),可得λi,i=0,1,…,N的值。具體的算法步驟為：①對樣本數(shù)據(jù)進(jìn)行排序；②計算概率加權(quán)矩；③建立殘差優(yōu)化模型，確定優(yōu)化初始值；④應(yīng)用遺傳算法對模型進(jìn)行求解；⑤判斷是否收斂，若收斂則轉(zhuǎn)⑦，否則轉(zhuǎn)⑥；⑥重新設(shè)定初始值，轉(zhuǎn)④；⑦輸出計算結(jié)果；⑧計算結(jié)束。

2 基于保障度的間斷工作導(dǎo)彈裝備備件需求模型的建立

為建立導(dǎo)彈裝備維修備件需求量模型作如下假設(shè)：

(1) 某導(dǎo)彈備件供應(yīng)保障系統(tǒng)由S+1個單元組成。當(dāng)導(dǎo)彈裝備處于工作狀態(tài)時，認(rèn)為其中1個單元工作，其余單元作為溫儲備。當(dāng)工作單發(fā)生故障時，由其中一個沒有發(fā)生故障的溫儲備單元替換；如果導(dǎo)彈裝備處于停機(jī)狀態(tài)，所有未發(fā)生故障的單元均作為溫儲備。

(2) 初始狀態(tài)，所有零部件及其備件均為新品，系統(tǒng)中1個零部件在工作，S個備件處于溫儲備用于替換故障件。

(3) 在每個階段狀態(tài)中，系統(tǒng)部件和備件只有正常/故障2種狀態(tài)，且相同零部件或備件壽命分布相同，零部件在工作期間和備件在儲備期間壽命服從同一分布。

(4) 部件只有2種狀態(tài)：“工作”或“不工作”，即維修保障周期為1年(365天)，裝備每周工作tw(h)，其余時間ts=168-tw(h)停機(jī)。導(dǎo)彈裝備工作時，其零部件處于工作狀態(tài)，相應(yīng)的備件處于儲備狀態(tài)。導(dǎo)彈裝備停機(jī)時，零部件的故障率與備件的故障率相同。

(5) 零部件及其備件均不可修，導(dǎo)彈裝備發(fā)生故障時進(jìn)行換件維修，且修復(fù)后性能如初，換件時間忽略不計。

(6) 導(dǎo)彈裝備如果在工作期間發(fā)生故障立即實施維修，如果在停機(jī)期間發(fā)生故障等到再次開機(jī)時實施維修。

(7) 裝備系統(tǒng)各組成器件的失效相互獨立，其失效不會發(fā)生在同一時刻。

(8) 任一器件在發(fā)現(xiàn)其不能工作之前總是完好的，即不能工作時間從故障發(fā)現(xiàn)時開始。

(9) 裝備數(shù)量及其所屬的零部件數(shù)量一定。

可推導(dǎo)出

(10)

(11)

式中:t1為第1個零部件發(fā)生故障前最后一次開機(jī)保持良好狀態(tài)的時間;i表示天數(shù)。

則第1個零部件的壽命可表示為

T1=24n1+t1.

(12)

如果處于儲備狀態(tài)的備件沒有發(fā)生故障，處于工作狀態(tài)的第1個零部件故障后直接用其替換，則第2個備件的累計故障分布函數(shù)可表示為

(13)

壽命T2=24(n1+n2)+t2。以此類推，第N個備件的累計故障分布函數(shù)為

(14)

通過迭代計算，可以推導(dǎo)出單元的壽命為TN=TN-1-(tw-tN-1)+24n+tN。其中tN為該零部件最后一次開機(jī)后良好運行的時間。

這樣依次循環(huán)直到所以備件都用完，就認(rèn)為導(dǎo)彈裝備不能工作，則最后1個沒有在儲備期間發(fā)生故障的零部件壽命就是備件的可保障時間。

3 仿真思路規(guī)劃

根據(jù)MC直接抽樣法[11]，可得到每個單元的累計故障概率抽樣值ηj，即[0,1]區(qū)間均勻分布的隨機(jī)數(shù)(j=1,2,…,S+1)。由于處于溫儲備狀態(tài)的備件也可能發(fā)生故障，所以當(dāng)?shù)?個零部件在工作期間發(fā)生故障需要進(jìn)行換件維修時，首先要判斷第2個單元是否已經(jīng)發(fā)生故障。如果沒有故障，則立即進(jìn)行換件維修；否則按以上方法對第3個單元進(jìn)行判斷。以此類推，直到所以備件都用完，就認(rèn)為導(dǎo)彈裝備不能工作。當(dāng)最后一個零部件的壽命不小于備件供應(yīng)保障周期時，就認(rèn)為備件的數(shù)量滿足需要，否則就是不滿足需要。根據(jù)Monte-Carlo思想，將上述過程重復(fù)Nmn次，記備件滿足需要的次數(shù)為Nsuc，當(dāng)Nmn足夠大時，Nsuc/Nmn即為備件的保障概率。針對故障概率和平均壽命的計算，進(jìn)行1 000次仿真運行即可，本文取Nmn=1 000次。

3.1 實例分析

根據(jù)調(diào)研分析，某型導(dǎo)彈裝備及其零部件平時主要處于停機(jī)和使用2種狀態(tài)，即工作和儲備2種模式，裝備平均每周使用大約12 h，其余時間處于停機(jī)儲備狀態(tài)，且裝備和備件的儲備條件一樣。裝備發(fā)生故障后采用換件維修的方式。對于該型導(dǎo)彈裝備某電子模塊工作和儲備期間的故障情況，這里統(tǒng)計分析了有連續(xù)故障記錄的20個樣本，以該模塊的故障間隔時間表示，單位為h，具體如表1所示。要求備件保障度為0.9。

表1 故障統(tǒng)計表

由于該電子模塊故障數(shù)據(jù)較少，適合采用基于概率加權(quán)矩的最大熵方法確定壽命分布函數(shù)。其中，所用矩的階數(shù)m取為2，具體的m選取參見文獻(xiàn)[7-8]。在構(gòu)造遺傳算法時，由于目標(biāo)函數(shù)是求最小值，同時為增大適應(yīng)度函數(shù)間的差距，以提高優(yōu)秀個體被選中的概率，因此選擇目標(biāo)函數(shù)值倒數(shù)的1010倍作為遺傳算法的適應(yīng)度函數(shù)，求得其參數(shù)為：λ0=6.25×10-3，λ1=4.49×10-3，λ2=2.63×10-3，其算法尋優(yōu)過程如圖1所示。經(jīng)擬合優(yōu)度檢驗知，該電子模塊工作壽命服從參數(shù)為λ=6.25×10-3的指數(shù)分布。在此基礎(chǔ)上，運用隨機(jī)生成服從相同參數(shù)分布的50組數(shù)據(jù)樣本，分別基于概率密度函數(shù)和逆累積分布函數(shù)來計算1階矩，并用在定義域內(nèi)的積分對于理論總體均值的相對誤差來評判不同方法的優(yōu)劣，具體計算結(jié)果如圖2所示。圖2表明，當(dāng)樣本量少于10時，采用經(jīng)典矩方法相對誤差達(dá)到50.033 7%，而采用概率加權(quán)矩方法僅為3.969 4%。因此，當(dāng)樣本≤30時，基于經(jīng)典矩的最大熵方法不太適用；若樣本量>30，兩種方法均可行，但概率加權(quán)矩方法精確度相對高些。同理，可求得儲備壽命為服從μ=1.851 8×10-4的指數(shù)分布。

圖1 適應(yīng)度函數(shù)圖Fig.1 Fitness function figure

圖2 不同樣本下2種方法計算相對誤差Fig.2 Bias of two different methods obtained from various sample

3.2 需求仿真

利用Matlab編程進(jìn)行1 000次仿真，該模塊1年的備件需求量部分計算結(jié)果如表2所示。其中，該導(dǎo)彈裝備每周工作20 h時，1年所需備件需求量如圖3所示。

3.3 結(jié)果分析與評價

表2中，tw=0表示備件一直處于儲備狀態(tài)，可作為溫儲備看待，即所有單元均以μ=1.851 8×10-4h的故障率進(jìn)行儲備，可利用并聯(lián)系統(tǒng)公式[12]進(jìn)行計算，結(jié)果為S=1，P=0.928 6。tw=24可視為連續(xù)工作的熱儲備系統(tǒng)，利用熱儲備系統(tǒng)公式[12]可解得S=6，P=0.928 1。由此可見，2種方法的計算結(jié)果是一致的，這說明了本文仿真方法的正確性。

表2 備件需求量計算結(jié)果

圖3 備件需求量仿真圖Fig.3 Simulation figure of spare parts demand

4 結(jié)束語

針對導(dǎo)彈裝備故障數(shù)據(jù)樣本量較小的情況，提出了基于最大熵原理與概率加權(quán)矩的分析方法，能夠在確定概率分布的同時估計分布參數(shù)。該方法的優(yōu)勢在于適用范圍較廣，而且與基于經(jīng)典矩的最大熵方法相比，對于樣本量小于30的情形擬合效果更為理想。備件需求模型既分析了導(dǎo)彈裝備工作期間發(fā)生情況，也考慮了導(dǎo)彈裝備在停機(jī)期間及備件在儲備期間也會發(fā)生故障的情況，更加貼近實際。

參考文獻(xiàn)：

[1] PANDEY M D. Direct Estimation of Quantile Functions Using the Maximum Entropy Principle [J]. Structural Safety, 2000,22(1): 61-79.

[2] 俞禮軍,嚴(yán)海,嚴(yán)寶杰. 最大熵原理在交通流統(tǒng)計分布模型中的應(yīng)用[J].交通運輸工程學(xué)報, 2001, 1(3): 91-94.

YU Li-jun, YAN Hai, YAN Bao-jie. Maximum Entropy Principle and Its Application in Probability Density Function of Traffic Flow [J].Journal of Traffic and Transportation Engineering, 2001, 1 (3): 91-94.

[3] HOSKING R M, WALLIS J R. Regional Frequency Analysis[M].New York: Cambridge University Press, 1997.

[4] TAN L, TANIAR D. Adaptive Estimated Maximum Entropy Distribution Model[J].Information Sciences, 2007,15(177):3110-3128.

[5] 俞禮軍,徐建閩. 出行時間價值最大熵分布估計模型[J].交通運輸工程學(xué)報, 2008, 8(1):83-96.

YU Li-jun, XU Jian-min. Estimation Models of Distribution Functions for Travel Time Value with Maximum Entropy Principle[J]. Journal of Traffic and Transportation Engineering, 2008, 8(1):83-96.

[6] 俞禮軍, BenHeydecker. 基于概率加權(quán)矩與最大熵原理的交通流統(tǒng)計分布估計方法[J].公路交通科技, 2008, 25(2):113-133.

YU Li-jun, BenHeydecker. Probability Density method of Traffic Flow Based on Probability Weighted Moments and Maximum Entropy Principle[J]. Science and Technology of Road Transportation, 2008, 25(2):113-133.

[7] 朱堅民，郭冰菁，王中宇，等.基于最大熵方法的測量結(jié)果估計及測量不確定度評定[J].超重運輸機(jī)械，2005，42(8)：5-8.

ZHU Jian-min,GUO Bing-jing,WANG Zhong-yu. Study on Evaluation of Measurement Result and Uncertainty based on Maximum Entropy Method[J].Journal of Electronic Measurement and Instrument,2005,42(8):5-8.

[8] YU X M,TONG L. Evaluation of Uncertainty in Measurement Based on the Maximum Entropy Method[J].Journal of Electronic Measurement and Instrument,2004,17(6):98-103.

[9] 葛陽, 寧劍平. 間斷工作導(dǎo)彈裝備維修備件需求量仿真模型[J].起重運輸機(jī)械, 2010,22(12):42-48.

GE Yang, NING Jian-ping. Demand Simulation Model of the Missile Equipment Maintenance Parts in Discontinuous Work[J].Lift Transportation Mechanism, 2010,22(12):42-48.

[10] 金星, 文明, 李俊美. 壽命服從指數(shù)分布產(chǎn)品相關(guān)失效解析分析[J].裝備指揮技術(shù)學(xué)院學(xué)報, 2002, 13(4):37-39.

JIN Xing, WEN Ming, LI Jun-mei. Correlative Failure Analysis on Exponential Distribution Manufacture[J].Journal of Academy of Equipment, 2002, 13(4):37-39.

[11] 程文鑫, 陳立強(qiáng), 龔沈光. 基于蒙特卡洛法的艦船裝備戰(zhàn)備完好性仿真[J].兵工學(xué)報, 2006, 27(6):1090-1094.

CHENG Wen-xin, CHEN Li-qiang, GONG Shen-guang. Ship Equipment Operational Readiness Simulation Based on Monte Carlo Method[J].Acta Armamentarii, 2006, 27(6):1090-1094.

[12] 曹晉華, 程侃. 可靠性數(shù)學(xué)引論[M].北京:高等教育出版社, 2006.

CAO Jin-hua, CHENG Kan. Reliability Mathematics Generality Induction[M].Beijing: Higher Education Press, 2006.