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

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

0 引言

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

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

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

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

為了克服基于不同樣本量估計(jì)的模型對(duì)于系統(tǒng)統(tǒng)計(jì)特性的描述偏差相差極大的缺陷,本文引入概率加權(quán)矩(probability weighted moments, PWMs)方法。根據(jù)最大熵原理[4]，對(duì)于定義域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í)數(shù)。

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

(2)

(3)

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

由式(1)可得

(4)

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

(5)

求解得

(6)

式中：N為樣本矩的最高階數(shù)；λi,i=0,1,…,N為與第i階矩約束相對(duì)應(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個(gè)方程。

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

(8)

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

(9)

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

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

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

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

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

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

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

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

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

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

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

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

可推導(dǎo)出

(10)

(11)

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

則第1個(gè)零部件的壽命可表示為

T1=24n1+t1.

(12)

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

(13)

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

(14)

通過(guò)迭代計(jì)算，可以推導(dǎo)出單元的壽命為T(mén)N=TN-1-(tw-tN-1)+24n+tN。其中tN為該零部件最后一次開(kāi)機(jī)后良好運(yùn)行的時(shí)間。

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

3 仿真思路規(guī)劃

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

3.1 實(shí)例分析

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

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

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

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

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

3.2 需求仿真

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

3.3 結(jié)果分析與評(píng)價(jià)

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

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

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

4 結(jié)束語(yǔ)

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

參考文獻(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)計(jì)分布模型中的應(yīng)用[J].交通運(yùn)輸工程學(xué)報(bào), 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] 俞禮軍,徐建閩. 出行時(shí)間價(jià)值最大熵分布估計(jì)模型[J].交通運(yùn)輸工程學(xué)報(bào), 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ì)分布估計(jì)方法[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ān)民，郭冰菁，王中宇，等.基于最大熵方法的測(cè)量結(jié)果估計(jì)及測(cè)量不確定度評(píng)定[J].超重運(yùn)輸機(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] 葛陽(yáng), 寧劍平. 間斷工作導(dǎo)彈裝備維修備件需求量仿真模型[J].起重運(yùn)輸機(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é)報(bào), 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é)報(bào), 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.