吳捷，郭峰，張素琴

（海軍航空大學(xué)青島校區(qū)，青島 266041）

引言

在現(xiàn)代高科技戰(zhàn)爭中，航空兵是決定戰(zhàn)爭勝負(fù)的重要因素之一。作戰(zhàn)飛機(jī)對航材的依賴性極大，因此做好航材需求預(yù)測工作具有重要的軍事和經(jīng)濟(jì)意義。現(xiàn)代戰(zhàn)爭具有很強(qiáng)的變異性、復(fù)雜性以及多樣性，加之作戰(zhàn)飛機(jī)的頻繁使用及戰(zhàn)場環(huán)境惡劣，飛機(jī)的戰(zhàn)斗損傷不可避免，戰(zhàn)時航材消耗較平時會呈現(xiàn)許多不同的特點(diǎn)和規(guī)律。但由于戰(zhàn)時樣本量極少，很難用一般的預(yù)測方法進(jìn)行需求預(yù)測，并且預(yù)測精度很低，這就給戰(zhàn)時航材保障人員在進(jìn)行航材保障時帶來很大困難。戰(zhàn)時航材需求預(yù)測是一個十分復(fù)雜的過程，構(gòu)建一個精確的模型對航材需求數(shù)量進(jìn)行準(zhǔn)確預(yù)測極其困難，戰(zhàn)爭不會重復(fù)發(fā)生，也不可能等到下一場戰(zhàn)爭爆發(fā)再用其來檢驗預(yù)測的準(zhǔn)確性。鑒于此，筆者以MATLAB為仿真平臺，結(jié)合蒙特卡洛方法以及馬爾科夫方法對戰(zhàn)時航材需求進(jìn)行計算機(jī)仿真模擬，以期預(yù)測戰(zhàn)時作戰(zhàn)飛機(jī)航材需求變化趨勢，為決策部門制定航材供應(yīng)數(shù)量提供一定的參考依據(jù)。

1 馬爾可夫模型

1.1 馬爾可夫預(yù)測法原理

馬爾可夫預(yù)測法是對事件下一時刻各種狀態(tài)發(fā)生的概率進(jìn)行預(yù)測的方法，通過對當(dāng)前狀態(tài)的分析來預(yù)測其未來各個階段狀態(tài)變化。假設(shè)一個事件所有可能的狀態(tài)分別記為E1，E2，…，En且該系統(tǒng)每次只能處于其中的一種狀態(tài)中，此事件下一步的狀態(tài)只有在時間流轉(zhuǎn)到固定且順序排列的時刻t上才會發(fā)生變化?？梢酝ㄟ^一系列的隨機(jī)狀態(tài)參數(shù)Xn（n=1，2，…）來描述，其中當(dāng)Xn=k，表示在t=n時刻，事件所處的狀態(tài)為

在事件各部分相互作用動態(tài)變化的過程中，從t時刻的狀態(tài)啟動，經(jīng)過k個時刻變化到下一狀態(tài)（可能為相同狀態(tài)也可能不同）的概率，稱作k步狀態(tài)轉(zhuǎn)移概率。經(jīng)過Ei變?yōu)镋j的狀態(tài)轉(zhuǎn)移概率簡記為：

假定事件所有可能存在的狀態(tài)共有n個，為E1，E2，…，En。把Pij作為矩陣的第i行第j列，則n個狀態(tài)（j=1，2，3，…，n）共有n行，則由狀態(tài)Ei轉(zhuǎn)變?yōu)闋顟B(tài)Ei的狀態(tài)轉(zhuǎn)移概率矩陣為：

1.2 戰(zhàn)時航材需求的馬爾科夫性分析

在戰(zhàn)爭爆發(fā)的時候，存在相當(dāng)數(shù)量的航材，其每個月的消耗數(shù)量是一個很小的正整數(shù)（其數(shù)量一般不超過5個），此種作戰(zhàn)所需航材的需求量是一個隨機(jī)的、不確定的量，并在戰(zhàn)爭持續(xù)時間內(nèi)表現(xiàn)出隨機(jī)波動的特點(diǎn)。時間序列預(yù)測法及灰色系統(tǒng)預(yù)測法均不能有效地預(yù)測其需求量。過去的戰(zhàn)爭數(shù)據(jù)和平時的訓(xùn)練數(shù)據(jù)只能作為一個參考，戰(zhàn)爭的演變速度很快，對戰(zhàn)時物資的需求日新月異，過去的數(shù)據(jù)信息對現(xiàn)在價值不大，因此我們可以假設(shè)此類航材未來的需求量的情況與“過去”的情況是無關(guān)的，只與當(dāng)前狀態(tài)有關(guān)，即該類航材需求量具有馬爾可夫性，馬爾科夫性也稱為無后效性[3]。事物第N次出現(xiàn)的狀態(tài)只與其第N-1次的狀態(tài)有關(guān)，它與再往前的歷史狀態(tài)無關(guān)[4]。

1.3 馬爾科夫模型分布研究

由該類航材馬爾可夫性可知，過程在時刻s處于狀態(tài)i條件下，在隨后的t個單位時間中過程仍不離開狀態(tài)i的概率，是它處于i至少t個單位的無條件概率[5-7]。若記hi為過程在轉(zhuǎn)移到另一個狀態(tài)之前停留在狀態(tài)i的時間，則對一切s，t≥0有：

由此可見，隨機(jī)變量hi具有無記憶性，對于一條馬爾可夫鏈，由于戰(zhàn)爭不會無限制地持續(xù)下去，我們僅截取馬爾可夫鏈中符合戰(zhàn)爭持續(xù)時間的一部分，每當(dāng)其演變到狀態(tài)i，在事件變化到下一個狀態(tài)的時間內(nèi)，馬爾科夫模型服從具有特定參數(shù)的指數(shù)分布[8，9]。

若已經(jīng)知道當(dāng)前月份的需求量，便可以利用前述矩陣P來計算未來有限月份內(nèi)航材需求量的各種狀態(tài)的發(fā)生概率。對未來3個月（即假定戰(zhàn)爭持續(xù)時間為三個月）的航材需求量各種狀態(tài)的概率累加就能夠得出戰(zhàn)爭持續(xù)時間內(nèi)航材需求量的概率分布情況。

1.4 馬爾科夫模型預(yù)測流程

第一步以時間t劃分階段，航材在某一時間段內(nèi)的消耗情況作為狀態(tài)量，該種低需求航材的月需求量為X，則該種低需求航材在狀態(tài)t1，t2，…，tn的取值分別為x1，x2，…，xn，該種航材需求數(shù)據(jù)中最大值記為Xmax，最小值記為Xmin，航材需求量為整數(shù)則，該種航材共有Xmax-Xmin+1種狀態(tài)。

第二步，根據(jù)航材樣本實際情況，用頻率近似替代概率，設(shè)Yi為處于狀態(tài)Ei的樣本數(shù)；Yij（k）為狀態(tài)Ei，經(jīng)過k步轉(zhuǎn)移到狀態(tài)Ej的數(shù)據(jù)個數(shù)，用Pij（k）表示k步轉(zhuǎn)移概率，如下所示：

則相應(yīng)的k步狀態(tài)轉(zhuǎn)移概率矩陣P（k）就可以表示為：

第三步，求解任意階段狀態(tài)概率向量，當(dāng)前時間航材需求個數(shù)即為當(dāng)前狀態(tài)，并通過之前的戰(zhàn)備訓(xùn)練記錄數(shù)據(jù)計算出轉(zhuǎn)移概率矩陣，就能夠預(yù)測出當(dāng)前時間往后k步航材數(shù)量所處于的各種狀態(tài)，用C（t）表示戰(zhàn)爭爆發(fā)后第t個月的航材各需求量出現(xiàn)的可能性矩陣，預(yù)測公式如下所示：

2 蒙特卡羅法

2.1 蒙特卡羅算法原理

蒙特卡羅算法是一種以大數(shù)定律為基礎(chǔ)，通過大量試驗來模擬隨機(jī)現(xiàn)象的方法，當(dāng)模擬次數(shù)達(dá)到一定數(shù)量后，模擬數(shù)據(jù)的特征就越接近現(xiàn)實情況。蒙特卡羅算法需要根據(jù)數(shù)據(jù)規(guī)律，把確定性問題與某個概率模型相聯(lián)系，用模擬試驗所得統(tǒng)計值代替無法精確計算的理論值，從而解決復(fù)雜問題[10-12]。

2.2 蒙特卡羅仿真基本框架

第一步，建立數(shù)據(jù)與目標(biāo)值的關(guān)系，構(gòu)造概率模型，對于本來不是隨機(jī)問題的確定性問題，需要人為地構(gòu)造一個概率分布模型。

第二步，實現(xiàn)從構(gòu)造概率分布的抽樣。為實現(xiàn)仿真，必須進(jìn)行隨機(jī)變量的抽樣，進(jìn)而得到符合構(gòu)造的概率分布的隨機(jī)序列，即樣本值。

第三步，建立各種統(tǒng)計量的估計。作為問題的解，可能是概率或期望。對于前者則用頻率代替，對于后者用樣本算術(shù)平均值代替，從而近似求解出預(yù)測值。

2.3 基于蒙特卡羅的航材狀態(tài)概率隨機(jī)模擬

首先分析戰(zhàn)時低需求航材的馬爾科夫性，通過MATLAB內(nèi)置的unifrnd (a，b，m，n)函數(shù)，產(chǎn)生m×n個在[a，b]上服從均勻分布的偽隨機(jī)數(shù)，將產(chǎn)生的結(jié)果與馬爾科夫模型產(chǎn)生的航材各需求量狀態(tài)概率進(jìn)行擬合。為了使蒙特卡羅模擬效果更具可靠性，更加符合實際，使用MATLAB內(nèi)置exprnd（MU，m，n）函數(shù)，先根據(jù)戰(zhàn)前積累的航材歷史數(shù)據(jù)，計算其統(tǒng)計平均值MU，接著計算機(jī)生成m×n個均值為MU的服從指數(shù)分布的隨機(jī)數(shù)序列。在剔除波動較大的數(shù)據(jù)序列后，將同一時間點(diǎn)的隨機(jī)數(shù)與真實數(shù)據(jù)相減，將其中方差最小的一組隨機(jī)模擬數(shù)據(jù)作為計算機(jī)模擬的戰(zhàn)時航材需求序列，與原始數(shù)據(jù)可視化對比，進(jìn)行誤差分析。具體步驟如下：

1）首先在MATLAB中以矩陣的形式輸入航材需求歷史數(shù)據(jù)，編寫程序找出其中最大值和最小值，計算出需求量可能存在的狀態(tài)。令t=0，記錄所需要預(yù)測的未來C個月的個數(shù)C。

2）統(tǒng)計每種狀態(tài)數(shù)量，以及由該種狀態(tài)下一步轉(zhuǎn)移狀態(tài)數(shù)量，計算出第1步狀態(tài)轉(zhuǎn)移概率矩，P（1），并以此類推計算出第2、3、…、C步狀態(tài)轉(zhuǎn)移矩陣。

3）根據(jù)最后一個狀態(tài)的樣本值，計算出未來C個月各個狀態(tài)對應(yīng)需求量的概率預(yù)測值。

4）令t=t+1，使用unifrnd函數(shù)生成C個在[0，1]上服從均勻分布的隨機(jī)數(shù)。

5）將概率擬合轉(zhuǎn)換為具體的航材需求數(shù)量，若生成的隨機(jī)數(shù)大于0且小于等于變?yōu)榈谝环N狀態(tài)的概率，這時該航材的需求量為0，若生成的隨機(jī)數(shù)大于變?yōu)榈谝环N狀態(tài)的概率且小于等于變?yōu)榍皟煞N狀態(tài)的概率之和，則該航材的需求量為1，以此類推求得第2、3、…、C個月的該航材需求量。

6）將這C個月的航材需求量求和，分別統(tǒng)計其中0、1、…、C×Xmax出現(xiàn)次數(shù)。

7）當(dāng)t小于規(guī)定的模擬次數(shù)M時，流程回到步驟2），當(dāng)t大于等于規(guī)定的模擬次數(shù)M時仿真停止，將步驟6）中各需求量出現(xiàn)的個數(shù)除以總仿真次數(shù)M，輸出結(jié)果并繪制圖形。

8）使用exprnd函數(shù)產(chǎn)生隨機(jī)航材需求數(shù)量序列，用該數(shù)據(jù)組與通過MATLAB產(chǎn)生的服從指數(shù)分布的航材需求數(shù)據(jù)組進(jìn)行逐項對比，檢驗?zāi)Ｐ皖A(yù)測效果。

馬爾科夫模型蒙特卡羅方法仿真流程如圖1所示。

圖1 馬爾科夫模型蒙特卡羅方法仿真流程圖

3 實例分析

已知某型飛機(jī)裝備實力與飛行任務(wù)量的變動不大，飛機(jī)上裝載的某型蓄壓器需求數(shù)據(jù)樣本統(tǒng)計時限為2019年4月~2020年7月，如表1所示。

由表1可知數(shù)據(jù)中最大值為2，最小值為0，該型蓄壓器需求數(shù)量共有三種狀態(tài)，對時間序列進(jìn)行統(tǒng)計，用頻率近似地表示一步轉(zhuǎn)移概率，由以上數(shù)據(jù)知，需求數(shù)為0的月份有2、6、9、10、12、14、15共7個，由需求數(shù)為0經(jīng)過一步轉(zhuǎn)移到需求數(shù)0的數(shù)量為2，則

表1 需求數(shù)據(jù)樣本

同理可以計算出一步概率轉(zhuǎn)移矩陣P（1）：

以此類推，利用計算MATLAB計算出P（2）、P（3）如下：

根據(jù)需求時間序列可知當(dāng)前此航材該月需求量為1，根據(jù)概率轉(zhuǎn)移矩陣P（1）、P（2）、P（3）可以計算出未來第1、2、3個月的航材需求狀態(tài)概率分布，利用MATLAB模擬1×104次，統(tǒng)計戰(zhàn)爭持續(xù)時間內(nèi)航材需求量為0的狀態(tài)出現(xiàn)次數(shù)為326，航材需求量為1的狀態(tài)出現(xiàn)次數(shù)為591，以此類推，除以仿真次數(shù)，得到航材總需求的累加概率，結(jié)果如表2所示。

表2 狀態(tài)概率對應(yīng)表

根據(jù)航材需求量和累計概率的對應(yīng)關(guān)系，可以繪制出在戰(zhàn)斗持續(xù)的三個月時間內(nèi)該航材的數(shù)量所對應(yīng)的戰(zhàn)時航材需求的滿足率，如圖2所示，由圖可知，該航材在戰(zhàn)時三個月內(nèi)需求量為7個的時候，對應(yīng)的需求滿足率達(dá)到了96 %以上，我們不必追求百分之百的滿足率而在戰(zhàn)前儲備10個該型蓄壓器，為了很小的滿足率的提高，需要投入巨大的人力成本和經(jīng)濟(jì)成本，雖然相較于經(jīng)濟(jì)成本來說，軍事效益是擺在第一位的，但是我們可以在戰(zhàn)爭開始的時候先儲備7個蓄壓器，再依據(jù)戰(zhàn)場態(tài)勢的變化靈活進(jìn)行補(bǔ)充。

圖2 航材個數(shù)與航材滿足率的關(guān)系

4 誤差分析

基于馬爾科夫模型和蒙特卡洛法的戰(zhàn)時航材需求預(yù)測，由于缺少真實的戰(zhàn)場需求數(shù)據(jù)作為誤差分析的依據(jù)，給誤差分析增添了許多的困難。筆者為了檢驗?zāi)M預(yù)測的可靠程度，根據(jù)現(xiàn)有蓄壓器需求數(shù)據(jù)，計算出其均值為0.75，其方差為0.60。使用MATLAB內(nèi)置exprnd函數(shù)，生成16×10 000個均值為0.75的服從指數(shù)分布的隨機(jī)數(shù)。將真實數(shù)據(jù)與仿真產(chǎn)生數(shù)據(jù)方差最小的序列作為隨機(jī)模擬蓄壓器需求量結(jié)果，繪制圖形如圖3所示。

圖3 原始數(shù)據(jù)與模擬數(shù)據(jù)對比圖

計算機(jī)模擬出的序列整體數(shù)據(jù)波動變小，大部分最值點(diǎn)在時間趨勢變化上與原需求序列擬合較好。依據(jù)指數(shù)分布所模擬出來的均值、方差等指標(biāo)與原始數(shù)據(jù)相差不是很大，說明該種蓄壓器需求量能較好的服從指數(shù)分布。根據(jù)新的蓄壓器需求序列模擬出未來第1、2、3個月的需求量分別為1.691 0、2.258 7、0.934 7，未來三個月的蓄壓器總需求量為4.884 4，可以得知在戰(zhàn)爭開始前儲備好的7個蓄壓器較為合理，能夠大概率的保證戰(zhàn)時該型蓄壓器的供應(yīng)。如果戰(zhàn)爭的持續(xù)時間超過三個月，可以繼續(xù)使用轉(zhuǎn)移概率矩陣進(jìn)行需求預(yù)測，但是會隨著預(yù)測時間的延長，導(dǎo)致轉(zhuǎn)移概率矩陣變化增大，從而導(dǎo)致誤差增大，后續(xù)還需要進(jìn)行誤差的修正。馬爾科夫模型和蒙特卡洛結(jié)合的預(yù)測方法整體模擬效果較好，雖然存在一些誤差，但是大體上的預(yù)測結(jié)果還是較為可靠的。

5 結(jié)語

傳統(tǒng)的戰(zhàn)時航材需求預(yù)測模式為了保證戰(zhàn)時航材供應(yīng)充足，往往會將所需航材大量儲備起來，從而造成航材積壓，消耗大量的資源。本文通過對馬爾科夫模型和蒙特卡羅仿真研究，表明了其可以為戰(zhàn)時航材需求量的決策提供一定參考，并以一個簡單的實例進(jìn)一步說明了其在戰(zhàn)時航材需求預(yù)測方面的應(yīng)用，最后通過誤差分析，對有限的歷史數(shù)據(jù)進(jìn)行仿真處理，以檢驗高預(yù)測的可靠性。同時，因為該預(yù)測方法簡單實用，只需輸入有關(guān)數(shù)據(jù)，就可以通過計算機(jī)軟件批量處理，且具有一定的可信度，所以還可以推廣到其它相似領(lǐng)域的需求預(yù)測。