孔子慶， 劉白楊， 劉 濟(jì)

（ 1. 華東理工大學(xué)信息科學(xué)與工程學(xué)院，上海 200237；2. 國家管網(wǎng)集團(tuán)北方管道有限責(zé)任公司壓縮機(jī)組維檢修中心，河北廊坊 065000）

由于需求數(shù)據(jù)的隨機(jī)性和非平穩(wěn)性，不常用備件的需求預(yù)測是庫存控制研究領(lǐng)域的一大難題。對于不常用備件的間斷型需求數(shù)據(jù)，基于時(shí)間序列分析的指數(shù)平滑法[1]、Croston法[2]、神經(jīng)網(wǎng)絡(luò)[3]等點(diǎn)預(yù)測方法精度較低，因此，基于概率統(tǒng)計(jì)理論的預(yù)測方法受到研究者的青睞。傳統(tǒng)的統(tǒng)計(jì)預(yù)測方法選擇常見的概率分布函數(shù)來描述備件的需求規(guī)律，例如Negative exponential分布、Normal分布及Poisson分布等[4]。但是，需求分布函數(shù)的選擇缺乏有效的理論依據(jù)，而且需求數(shù)據(jù)常常并不滿足這些常見的概率分布。對此，許多研究者選擇通過一種增廣樣本統(tǒng)計(jì)方法——Bootstrap算法來估計(jì)需求量的概率分布情況。如Willemain等[5]將Markov模型與Bootstrap法相結(jié)合，得到提前期內(nèi)需求量的概率分布。Aris等[6]對文獻(xiàn)[5]方法與指數(shù)平滑法進(jìn)行了比較，肯定了文獻(xiàn)[5]方法的優(yōu)越性。但文獻(xiàn)[5]僅考慮需求發(fā)生與不發(fā)生兩狀態(tài)之間的轉(zhuǎn)移概率，基于此概率預(yù)測提前期需求發(fā)生次數(shù)，預(yù)測精度還有提升空間。

另一方面，一些學(xué)者將研究重點(diǎn)聚焦于庫存優(yōu)化上，通過設(shè)置訂購點(diǎn)、安全庫存等方法來提高需求預(yù)測的容錯性和庫存控制的魯棒性。如Yonit 等[7]通過連續(xù)檢查庫存(Q，R) 模型，最大程度地降低庫存總成本。Wang等[8]則通過構(gòu)建(T，S)模型來應(yīng)對備件的供應(yīng)風(fēng)險(xiǎn)。這些庫存優(yōu)化模型都是建立在需求服從某一已知概率分布的假設(shè)上，缺乏合理性。針對間斷型需求預(yù)測的難點(diǎn)，另有一些學(xué)者將需求預(yù)測與庫存優(yōu)化結(jié)合起來。如Ehrenthal等[9]通過研究需求數(shù)據(jù)中的季節(jié)性因素，證明了在庫存管理中考慮非平穩(wěn)需求可有效降低庫存成本和缺貨成本。劉琴等[10]將需求分布預(yù)測與(s，S)模型相結(jié)合，構(gòu)建備件庫存優(yōu)化模型，但該模型中的需求分布是基于文獻(xiàn)[5]方法得到的，模型的準(zhǔn)確性有待提高。

本文提出了一種將需求預(yù)測與庫存優(yōu)化相結(jié)合的新方法，采用高斯過程回歸 (Gaussian Process Regression，GPR) 方法對需求發(fā)生間隔進(jìn)行預(yù)測，結(jié)合Bootstrap算法對需求量序列進(jìn)行重采樣，預(yù)測提前期內(nèi)需求量的概率分布，基于此預(yù)測結(jié)果建立庫存總成本最小的優(yōu)化模型，使用粒子群算法尋優(yōu)得到最優(yōu)庫存決策結(jié)果。

1 高斯過程回歸方法

不常用備件的庫存決策流程如圖1所示。其中，高斯過程回歸方法用于需求間隔預(yù)測。

圖1 不常用備件庫存決策流程Fig. 1 Rarely used spare parts inventory decision process

高斯過程（Gaussian Process，GP）是指一個隨機(jī)過程，其中任意隨機(jī)變量的有限子集均服從聯(lián)合高斯分布[11]。

對于一個給定樣本數(shù)據(jù)集D=(xi，yi)，i=1，2，…，n，其中輸入數(shù)據(jù)xi∈Rd，為d維輸入量，輸出數(shù)據(jù)yi∈R。在所給數(shù)據(jù)集D的有限集合中，f(x1)，…，f(xn)構(gòu)成一個具有聯(lián)合高斯分布的隨機(jī)變量集，此時(shí)高斯過程的所有統(tǒng)計(jì)特性由均值函數(shù)m(x)與協(xié)方差函數(shù)k(xi,xj)構(gòu)成，即

其中：xi，xj∈Rd；f(x)表示均值函數(shù)為m(x)、協(xié)方差函數(shù)為k(xi，xj)的高斯過程函數(shù)。

建立一個包含噪聲的高斯過程模型

當(dāng)m(x)=0時(shí)，有限觀測值y的先驗(yàn)分布為

其中：δij是Kronecker delta函數(shù)。

將觀測值y的協(xié)方差以矩陣的形式表示：

其中：X=(x1,…,xn)T；C(X，X)表示n×n的協(xié)方差矩陣；In表示n×n的單位矩陣；K(X，X)為n×n的對稱正定協(xié)方差矩陣，其中的矩陣元素Kij=k(xi，xj)。

給定測試樣本x?,則觀測值y與其預(yù)測值y*的聯(lián)合先驗(yàn)分布為

由此可以計(jì)算出y*的條件概率分布為

其中

GPR的預(yù)測結(jié)果與協(xié)方差函數(shù)(高斯核函數(shù)）的選取密切相關(guān)。一般選取核函數(shù)為平方指數(shù)協(xié)方差：

其中：λ為距離尺度；為協(xié)方差函數(shù)的信號方差。通過極大似然法求取超參數(shù)集合 θ={λ,σf,σn} 的最優(yōu)值。使用共軛梯度法、牛頓法等優(yōu)化方法求得令似然函數(shù)L(θ)極小的超參數(shù)最優(yōu)值。L(θ)的形式如下：

2 基于GPR的需求間隔預(yù)測

需求間隔是指相鄰兩次需求發(fā)生的時(shí)間間隔長度?；跁r(shí)間序列的需求預(yù)測算法大都將需求量和需求間隔分離開各自進(jìn)行預(yù)測[2]，這樣，充滿了許多零值的間斷型原始數(shù)據(jù)便剝離為兩個連續(xù)時(shí)間序列?需求量序列和需求間隔序列，可以用各種常規(guī)的時(shí)間序列回歸和預(yù)測方法進(jìn)行分析。

根據(jù)Kolmogorov定理可知，任何一個時(shí)間序列都可以看成是由一個非線性機(jī)制確定的輸入輸出系統(tǒng)，因此，可以將任何時(shí)間序列預(yù)測問題當(dāng)作一個非線性系統(tǒng)回歸問題來解決。需求間隔序列為連續(xù)時(shí)間序列，具有一定統(tǒng)計(jì)規(guī)律，可以近似為常見的高斯過程。因此，本文通過GPR進(jìn)行需求間隔序列的預(yù)測。

針對間斷型需求時(shí)間序列Z=(z1，z2，···，zn)，可直接抽取得到需求間隔序列G=(g1，g2，···，gt)。其中zi表示第i期的需求量，i=1，2，···，n；gj表示相鄰兩期非零需求的時(shí)間間隔，j=1，2，···，t。

首先選取需求間隔序列的預(yù)測滑動窗數(shù)據(jù)gt?β,gt?(β?1)，···，gt，建立下一個需求間隔值gt+1與該輸入序列之間的高斯回歸模型

其中：β為參與預(yù)測的數(shù)據(jù)滑動窗參數(shù)，則β+1為滑動窗大小，即預(yù)測維數(shù)。恰當(dāng)?shù)剡x取預(yù)測維數(shù)是數(shù)據(jù)預(yù)測方法的關(guān)鍵。

然后，以需求間隔樣本序列G=(g1，g2，···，gt)作為學(xué)習(xí)樣本，訓(xùn)練高斯過程回歸模型，獲得模型超參數(shù)最優(yōu)解，便可以獲得需求間隔預(yù)測值t+1及其方差。

3 基于GPR和Bootstrap的需求概率分布預(yù)測

Bootstrap算法是一種增廣樣本統(tǒng)計(jì)推斷技術(shù)，其基本思想是在原始數(shù)據(jù)集中進(jìn)行多次有放回的隨機(jī)抽樣，每次抽取N個數(shù)據(jù)，構(gòu)成觀測樣本，然后根據(jù)觀測樣本生成反映原始數(shù)據(jù)的抽樣分布。Williamain等[5]根據(jù)需求發(fā)生與不發(fā)生的歷史轉(zhuǎn)移概率，對整個需求數(shù)據(jù)進(jìn)行抽樣并放回，相較于前人的工作有較大的提升，但是對歷史數(shù)據(jù)的挖掘并不充分，其后相應(yīng)的改進(jìn)研究非常稀少[10]。

本文提出了一種新的GPR與Bootstrap結(jié)合的需求概率預(yù)測方法，其基本思想是：利用GPR循環(huán)預(yù)測需求發(fā)生間隔均值及其方差，由此進(jìn)行隨機(jī)抽樣，推斷提前期內(nèi)需求發(fā)生的次數(shù)。基于提前期需求次數(shù)，采用Bootstrap算法對歷史需求量序列進(jìn)行重采樣，獲得多組提前期需求量預(yù)測序列，根據(jù)這些預(yù)測序列便可估計(jì)提前期內(nèi)需求量的概率密度和累積分布。具體步驟如下：

(1)假定提前期為L，初始化需求間隔預(yù)測序列M=[]，令i，j=0。

(2)生成預(yù)測滑動窗(gt?(β?i)，gt?(β?1?i)，···，gt+i)，通過GPR方法預(yù)測需求間隔分布N(t+i+1,t+i+12)。

(3)對需求間隔分布進(jìn)行隨機(jī)抽樣，獲得下一時(shí)刻需求間隔預(yù)測值rj，將rj存入需求間隔序列M中，此時(shí)M=(r0，r1，…，rj)。判斷預(yù)測需求間隔之和是否大于提前期L，若r0+r1+…+rj>L，則將rj從序列M中剔除，并跳轉(zhuǎn)到步驟(4)。若r0+r1+…+rj＜L，則令i=i+1，j=j+1，gt+i=rj，跳轉(zhuǎn)到步驟（2）。

(5)重復(fù)步驟（1）～（4）共m次，得到m組DL，由其擬合得到DL的概率密度函數(shù)f(DL)和累積概率分布函數(shù)F(DL)。

具體算法流程如圖2所示。

圖2 需求量概率分布預(yù)測流程圖Fig. 2 Flow chart of probabilistic distribution forecast

4 基于需求概率分布的隨機(jī)庫存優(yōu)化建模

由于間斷型需求數(shù)據(jù)的隨機(jī)性，任何一種基于時(shí)間點(diǎn)預(yù)測的方法都難以保證其準(zhǔn)確性，因此本文避開需求量時(shí)間序列預(yù)測的難題，從庫存控制決策的角度出發(fā)，尋求提前期內(nèi)需求量的統(tǒng)計(jì)概率[12]，并以此為條件，從統(tǒng)計(jì)學(xué)的一般意義上來建立庫存成本優(yōu)化模型，讓提前期需求量的累積概率作為庫存成本控制的一個約束條件，以此獲得有效的庫存決策。

以連續(xù)檢查庫存的定點(diǎn)訂購策略(Q，R)為例[7]，定量化隨機(jī)性需求所帶來的缺貨損失，建立缺貨不補(bǔ)情況下的庫存總成本優(yōu)化模型。假定提前期L不變、備件到貨入庫時(shí)間很短、備件的單次訂購費(fèi)A、儲存費(fèi)率h和缺貨費(fèi)率s基本不變。基于上述條件，定義DL為提前期內(nèi)需求量(隨機(jī)變量）；f(DL)為提前期需求量的概率密度函數(shù)；F(DL)為提前期需求量的累積概率分布函數(shù)；R為訂購點(diǎn)；Q為單次訂購量；p為備件的單價(jià)；h為單位時(shí)間單位價(jià)值的儲存費(fèi)率（%）；s為單位價(jià)值備件缺貨一次的缺貨費(fèi)率(%)。

假定在一個訂購周期ti內(nèi)，庫存的總成本由儲存成本、訂購成本、缺貨成本構(gòu)成。

(1)訂購成本CO

一般考慮一次訂購所產(chǎn)生的物流、運(yùn)輸、人員、稅務(wù)等費(fèi)用基本固定，即不考慮訂購費(fèi)因訂購量而打折促銷的因素。假設(shè)單次訂購發(fā)生的訂購成本為常數(shù)。

(2)缺貨成本CL

若單位價(jià)值的備件一次缺貨造成的損失成本為ps，則缺貨總成本為

(3)儲存成本CS

CS為備件儲存在倉庫中所產(chǎn)生的管理費(fèi)、場地費(fèi)等，與庫存量及其儲存時(shí)間有關(guān)。若發(fā)生缺貨，則平均庫存量的計(jì)算應(yīng)排除缺貨量，庫存時(shí)間也發(fā)生相應(yīng)變化[13]。訂貨周期ti內(nèi)的平均儲存成本為

這樣，一個訂購周期ti內(nèi)總成本為Ctotal(ti)為

將式(15)等式兩邊同時(shí)除以ti，可以得到單位時(shí)間內(nèi)總庫存成本為C，因此，以單位時(shí)間總庫存成本C最小為控制目標(biāo)的優(yōu)化問題描述為

其中：CSL為客戶要求的備件服務(wù)水平，代表備件庫存滿足需求的最低概率[14]；F?1表示F(DL)的反函數(shù)；Qmax、Rmax分別表示訂購量和訂購點(diǎn)的優(yōu)化邊界值。

庫存優(yōu)化模型的主要思想是：在決策變量Q和R的可行解空間(約束條件s.t.給出)，尋找令由訂購成本CO、儲存成本CS、缺貨成本CL之和構(gòu)成的庫存總成本最小的最優(yōu)解 (Q*,R*)，該最優(yōu)解代表了最優(yōu)訂購點(diǎn)和最優(yōu)訂購批量。

5 實(shí)例分析應(yīng)用

針對不常用備件的間斷型需求數(shù)據(jù)，本文提出了一種基于GPR和Bootstrap的需求概率預(yù)測方法，并以此為條件來建立隨機(jī)庫存優(yōu)化模型，以提高庫存控制決策的準(zhǔn)確性和魯棒性。借鑒文獻(xiàn)[15]的核電設(shè)備備件以及文獻(xiàn)[16]的動車組齒輪箱歷史需求數(shù)據(jù)，對本文算法進(jìn)行比較驗(yàn)證。以文獻(xiàn)[15]中核電設(shè)備的備件需求數(shù)據(jù)為例，初始化提前期L為5個月，其需求數(shù)據(jù)以月為單位進(jìn)行統(tǒng)計(jì)，共有132條數(shù)據(jù)，其中需求量為0的數(shù)據(jù)共91條，以第一次發(fā)生需求為起點(diǎn)進(jìn)行統(tǒng)計(jì)，平均需求間隔長度為3.075個月，需求量均值為0.47件，標(biāo)準(zhǔn)差為0.604件。

5.1 需求間隔預(yù)測實(shí)例驗(yàn)證

從第一次發(fā)生需求開始計(jì)算需求間隔，將需求序列Z轉(zhuǎn)換為需求間隔序列G=(1, 6, 5, 3, 1, 2, 1, 1, 1,18, 3, 5, 2, 2, 2, 1, 4, 2, 3, 2, 3, 1, 1, 1, 12, 3, 3, 2, 4, 2, 4,2, 2, 2, 3, 3, 1, 2, 4, 3)。將前37條數(shù)據(jù)作為訓(xùn)練數(shù)據(jù)，第38、39、40條數(shù)據(jù)作為測試數(shù)據(jù)，采用本文提出的需求間隔預(yù)測算法，得到需求間隔的預(yù)測結(jié)果如圖3所示。其中超參數(shù)集合 θ 的優(yōu)化值為λ=0.112,σf=1.156，σn=1.162。

圖3 基于GPR的備件需求間隔預(yù)測結(jié)果Fig. 3 Spare parts demand interval prediction results based on GPR

表1 本文方法與其他算法的比較Table 1 Comparison of the proposed method and other algorithms

從表1可以看出，本文方法的平均絕對誤差及均方根誤差均在0.5以下，預(yù)測誤差小于ES與SVM，預(yù)測效果較好。高斯過程回歸在進(jìn)行需求間隔預(yù)測的同時(shí)，還給出了預(yù)測值的置信區(qū)間（方差），可提高后續(xù)Bootstrap抽樣的可靠性，降低預(yù)測誤差帶來的影響。

5.2 需求量概率分布預(yù)測實(shí)例驗(yàn)證

將GPR與Bootstrap相結(jié)合，預(yù)測備件在提前期內(nèi)的需求量概率分布。仍然使用文獻(xiàn)[15]的備件需求數(shù)據(jù)，根據(jù)其需求序列Z的前123個月數(shù)據(jù)，可得到未來提前期內(nèi)需求量的概率分布及累積分布。將本文方法與文獻(xiàn)[10]中基于Markov與Bootstrap的方法進(jìn)行對比，需求量概率分布的對比如圖4所示，累積分布的對比如圖5所示。假定備件服務(wù)水平要求為90%以上，則提前期內(nèi)需求量累積概率達(dá)到90%的點(diǎn)設(shè)置為安全庫存量。分別設(shè)置提前期L=5、8、11，對本文方法與文獻(xiàn)[10]方法的預(yù)測結(jié)果進(jìn)行對比，結(jié)果如表2所示，其中，置信區(qū)間表示概率為95%以內(nèi)的需求量范圍。

圖4 提前期內(nèi)需求量概率直方圖Fig. 4 Probability histogram of quantity demanded in lead period

圖5 提前期內(nèi)需求量累積分布直方圖Fig. 5 Histogram of cumulative distribution of demand in lead period

表2 兩種方法的評估結(jié)果Table 2 Evaluation results of two methods

從實(shí)驗(yàn)結(jié)果可以看到，不同提前期設(shè)置下，需求量真實(shí)值均落在兩種方法的95%置信區(qū)間內(nèi)，但本文方法的置信區(qū)間比文獻(xiàn)[10]的置信區(qū)間更窄，精度更高。并且，安全庫存量比文獻(xiàn)[10]方法得到的更小，可顯著降低備件的庫存總費(fèi)用。

5.3 隨機(jī)庫存決策實(shí)例驗(yàn)證

采用GPR與Bootstrap方法獲得提前期內(nèi)需求量概率函數(shù)f(DL)和F(DL)，由此建立隨機(jī)庫存最優(yōu)化模型式(16)，求解得到最優(yōu)決策變量作為庫存決策依據(jù)。依舊采用文獻(xiàn)[15]中核電設(shè)備備件的需求數(shù)據(jù)，備件單價(jià)為23 881元，儲存費(fèi)率設(shè)置為1.5%/a。單次訂購的固定費(fèi)用設(shè)置A=94元，缺貨損失費(fèi)率s=120%/次。服務(wù)水平CSL=90%，即須保證提前期內(nèi)備件需求得到滿足的概率至少為90%。邊界值Rmax與Qmax均設(shè)置為50件。

粒子群算法作為一種啟發(fā)式算法，具有操作簡單、尋優(yōu)速度快的特點(diǎn)[18]，因此本文采用粒子群算法求解優(yōu)化模型式(16)。在種群大小為100，最大迭代次數(shù)為200，慣性權(quán)重為0.6，自我學(xué)習(xí)因子為0.5，群體學(xué)習(xí)因子為0.5條件下，得到最優(yōu)尋優(yōu)結(jié)果為Q*=5，R*=4。庫存總成本的迭代收斂過程如圖6所示，可以看出，按照所求得的決策量進(jìn)行庫存控制，可使庫存總成本最小，有效地控制了庫存總費(fèi)用。

圖6 月最小總成本優(yōu)化過程圖Fig. 6 Optimization process diagram of minimum total monthly cost

為了檢驗(yàn)?zāi)Ｐ颓蠼獾膸齑鏇Q策數(shù)據(jù)能否滿足實(shí)際需求，對歷史需求序列Z采用滑動窗的形式統(tǒng)計(jì)該備件提前期內(nèi)需求量的大小以及發(fā)生頻次，窗口大小根據(jù)提前期設(shè)定為5，每次滑動步數(shù)為1，統(tǒng)計(jì)信息如表3所示。

從表3可以得到，提前期內(nèi)需求量最大值為9，最小值為0，其中需求量小于或等于最優(yōu)訂購點(diǎn)R*=4的累積概率為91.41%，服務(wù)水平達(dá)到90%以上，滿足模型預(yù)設(shè)的服務(wù)水平要求。

表3 歷史需求量統(tǒng)計(jì)信息Table 3 Historical demand statistics

6 結(jié)束語

針對不常用備件需求不穩(wěn)定、波動性強(qiáng)的特點(diǎn)，本文提出了一種需求概率預(yù)測和庫存優(yōu)化結(jié)合的新方法，以提高庫存決策的準(zhǔn)確性和魯棒性。首先采用GPR預(yù)測需求發(fā)生間隔，接著結(jié)合Bootstrap進(jìn)行需求概率分布預(yù)測，最后以此為條件建立隨機(jī)庫存優(yōu)化模型并求解。實(shí)際備件的需求數(shù)據(jù)數(shù)值分析結(jié)果表明：基于GPR的需求間隔預(yù)測精度相比于指數(shù)平滑法、SVM算法都有一定的提高，同時(shí)還獲得了預(yù)測值的置信區(qū)間信息；GPR與Bootstrap結(jié)合的需求概率分布預(yù)測結(jié)果比文獻(xiàn)[10]的置信區(qū)間更窄，精度更高，預(yù)測結(jié)果下的庫存總費(fèi)用相比顯著減??；基于上述預(yù)測結(jié)果所建立的庫存優(yōu)化模型，在保證服務(wù)水平的前提下，庫存總成本收斂至最小值，且最優(yōu)決策變量數(shù)值完全符合歷史需求統(tǒng)計(jì)結(jié)果。