基于馬爾可夫到達過程的兩級可修備件（S-1，S）庫存優(yōu)化模型

陳　童，黎　放，狄　鵬

（海軍工程大學管理工程系，武漢430033）

基于馬爾可夫到達過程的兩級可修備件（S-1，S）庫存優(yōu)化模型

陳童，黎放，狄鵬

（海軍工程大學管理工程系，武漢430033）

本文以兩級可修備件庫存系統(tǒng)為研究對象，采用馬爾可夫到達過程（MAP）描述備件需求規(guī)律，考慮有限維修設施的情況，假設故障件維修時間、備件運輸時間以及采購時間均服從phase-type（PH）分布，建立了一種描述能力更強、解析計算性更好的（S-1，S）庫存優(yōu)化模型，并推導出系統(tǒng)缺貨量分布函數(shù)；然后通過算例演示了模型的優(yōu)化效果，驗證了模型的正確性和適用性。

（S-1，S）庫存策略；兩級庫存；可修備件；馬爾可夫到達過程

1　前言

在軍事裝備、民用高技術(shù)產(chǎn)品等領(lǐng)域，對于具有較高可靠性，并且流通量較低的昂貴可修備件，管理者通常選擇多級（S-1，S）庫存策略。Sherbrooke[1]最早根據(jù)這種庫存策略開發(fā)了METRIC模型，他選擇費用約束下的期望缺貨量最小為優(yōu)化目標，假設各級維修能力無限。Grave[2]針對Sherbrooke研究中設計了一種便于理解的近似計算方法。文獻[3，4]則針對METRIC模型無限修理能力的假設，研究了維修臺有限情況下的多級可修備件庫存模型。付興方等[5]研究了兩級庫存系統(tǒng)備件分配優(yōu)化方法。Lau和Song[6]采用非平穩(wěn)Poisson過程研究了維修能力有限的兩級可修庫存系統(tǒng)。這些研究均采用Poisson過程描述部件的失效過程，假設維修時間等服從指數(shù)分布，這樣的假設具有一定的不合理性[7]。Kim等[8]為了放寬假設條件，進一步假設維修時間服從一般分布，維修臺有限；但該研究仍是建立在需求到達是Poisson流的基礎上，因此假設條件仍略顯嚴格。

本文針對已有研究假設條件過于嚴格的情況，提出采用馬爾可夫到達過程（MAP）描述備件的需求情況，考慮有限維修設施的情況，假設故障件維修時間、備件運輸時間以及采購時間均服從一般分布，采用phase-type（PH）分布進行表示；建立了一個描述能力更強的兩級可修備件（S-1，S）庫存模型。采用MAP描述需求到達過程時，（S-1，S）策略通常并不是最優(yōu)策略，這是因為當需求過程比較平穩(wěn)時，適當增加訂貨延遲時間是有助于降低成本的。但在軍事和工業(yè)領(lǐng)域，備件缺貨除了對費用的影響，往往還會導致任務失敗或增加設備的安全風險，因此采用（S-1，S）策略更加穩(wěn)妥，也更符合實際情況。

利用MAP和PH分布開展研究具有如下好處。

1）將備件需求過程表示為MAP形式便可以利用“MAP類在[0，+∞）的概率空間上，對平穩(wěn)點過程類似稠密[9]”這一性質(zhì)。該性質(zhì)使得MAP類涵蓋了許多常用的點過程[10]，如泊松過程、PH更新過程等；并且實際的備件需求往往具有突發(fā)性，而MAP本身結(jié)構(gòu)的特點[11]使它能夠很好的描述和分析突發(fā)現(xiàn)象。因此，MAP非常適合描述備件的需求發(fā)生規(guī)律。

2）采用PH分布表示一般分布是因為：任何分布總可以選擇一個適當?shù)腜H分布把它擬合到任意精確的程度[12]，因此PH分布可以作為許多假設分布的一般表述，采用PH分布描述維修時間等能有效避免指數(shù)分布等經(jīng)典分布的局限性。

2　預備知識

首先介紹PH分布和MAP的有關(guān)概念。

定義1[7][0,+∞)上的概率分布函數(shù)F(?)稱為PH分布，當且僅當它是一個有限狀態(tài)馬爾科夫過程的吸收時間分布，有分布函數(shù)

式（1）中，T為m階方陣；α=(α1，α2，...，αm)為其瞬態(tài)的初始概率向量；e為元素均為1的m階列向量；(α，T)稱為該PH分布的m階表示。PH分布中的每一個瞬態(tài)稱為相位。

定義2[11]一個有限不可約馬爾可夫鏈具有狀態(tài)空間S={1，2，...，m}，設D為該馬爾可夫鏈的無窮小生成元。狀態(tài)i的逗留時間服從參數(shù)為λi的指數(shù)分布，在該狀態(tài)即將結(jié)束時，有下列兩個事件之一發(fā)生：

1）存在一個馬爾可夫鏈，使得從狀態(tài)i到 j（1≤i，j≤m）的轉(zhuǎn)移中，有一個事件到達的概率為pij(1)；

2）存在一個馬爾可夫鏈，使得從狀態(tài)i到 j（1≤i，j≤m；i≠j）的轉(zhuǎn)移中，沒有事件到達的概率為pij(0)。

可以看到只有通過一次事件到達，才能使得狀態(tài)i返回到自身。并且有：

D0和D1分別是控制沒有事件到達和有單個事件到達的狀態(tài)轉(zhuǎn)換矩陣，它們均為m階子隨機矩陣，可以給出無窮小生成元D=D0+D1。

然后，定義A(t)為在(0，t]內(nèi)事件到達的次數(shù)；J(t)為在時刻t嵌入馬爾可夫鏈所處的狀態(tài)，其狀態(tài)空間為{i:1≤i≤m}。則{A(t)，J(t)}為一個兩維馬爾可夫過程，其狀態(tài)空間為{(n，i):n≥0，1≤i≤m}。而{A(t)，J(t)}被稱為馬爾可夫到達過程。

MAP是半馬氏過程的特例，其狀態(tài)轉(zhuǎn)移概率矩陣為：

設π是無窮小生成元為D的穩(wěn)態(tài)概率向量，則有：πD=0,?πe=1。并且πD1e表示在MAP進入穩(wěn)態(tài)后，單位時間內(nèi)事件到達的期望個數(shù)。

下面給出Kronecker乘積的定義，它在后面的模型中被大量使用。

定義3[13]如果 A和B分別為m1×m2和n1×n2的矩陣，它們的 Kronecker乘積 A?B則為m1n1×m2n2的矩陣，且：

3　問題描述

本文假設如下。

1）一個中心倉庫通常負責N個基地的備件供應和維修，在中心倉庫和基地分別存有S0和Si（i=1，...，N）個某型備件，并且在中心倉庫和每個基地均有一個維修站。

2）基地i（i=1，...，N）的故障到達是一個MAP，具有m階表示（E0(i)，E1(i)）。

3）當基地i中出現(xiàn)一個備件需求時，若該基地中仍有庫存，則立刻滿足該需求?；豬備件需求處理過程示意圖見圖1。該更換下來的故障件有r的概率可在基地進行維修，否則需交給中心倉庫處理。

4）故障件交還中心倉庫得同時，基地向中心倉庫申領(lǐng)一個備件，若中心倉庫有庫存則立即下發(fā)，否則需要等待維修完成或訂貨到達；

5）在中心倉庫，故障件有p的概率可修，1-p的概率報廢；當發(fā)生報廢時，中心倉庫需向生產(chǎn)廠家訂購；

6）基地i（i=1，...，N）的維修站維修時間服從PH分布，具有ki階表示(αi，Ti)；中心倉庫維修站維修時間服從PH分布，具有k0階表示(α0，T0)；所有維修站均采用先到先服務原則；

7）采購時間服從PH分布，具有 kC階表示(αC，TC)；

圖1　基地i備件需求處理過程示意圖Fig.1　The disposal process for spare parts demand of base i

8）從中心倉庫向各基地的運輸時間服從PH分布，具有kiY階表示(αiY，TiY)；而從基地向中心倉庫運送故障件的時間則可以計入中心倉庫的維修時間或采購時間中，不單獨列出。

9）所有節(jié)點均選擇（S-1，S）補充策略。

4　系統(tǒng)特性

下面首先對整個系統(tǒng)的備件供應流程進行分析。

4.1備件供應流程

根據(jù)假設條件，可以將整個備件供應流程劃分為4個部分，如圖2所示。

圖2　備件供應流程劃分示意圖Fig.2　The partition of spare parts supply flow

令OPi、DIi和BOi分別表示基地i的工作部件數(shù)、維修隊長和缺貨量；和表示在中心倉庫的等待維修隊長和采購隊列隊長，則表示中心倉庫不可用件數(shù)量；而中心倉庫的缺貨量為BO0，其中基地i在中心倉庫未滿足訂單的比例為 fi；OH0和OHi分別表示中心倉庫和各基地的可用件數(shù)量；在運往基地i途中的備件數(shù)為TRi。下面研究這些參數(shù)之間的關(guān)系：

由上述關(guān)系式可知，BO0和BOi可以表示為多個隨機變量的卷積形式。而又由假設條件可知和TRi均是MAP/PH/1排隊系統(tǒng)的隊長，則通過研究MAP/PH/1排隊系統(tǒng)可知它們所服

從的分布。

4.2MAP/PH/1排隊系統(tǒng)

令ψ(t)=(S(t)，I(t)，J(t))，其中S(t)、I(t)和J(t)分別為在時刻t的系統(tǒng)顧客數(shù)，顧客到達過程相位和服務臺工作相位；另假設服務時間服從k階PH (a，T)分布。則為連續(xù)時間馬爾可夫鏈，其狀態(tài)空間 Ω可以表示為：Ω=h1?h2。其中，h1={(0，i)：1≤i≤m}表示系統(tǒng)內(nèi)沒有顧客；h2={(s，i，j)：1≤s，1≤i≤m，1≤j≤k}表示系統(tǒng)內(nèi)有s個顧客，其中一個顧客正在接受服務。

通常稱S(t)為狀態(tài)空間的宏狀態(tài)。對宏狀態(tài)按字典序排列，可以給出該馬爾可夫鏈的無窮小生成元矩陣：

式（9）中，B00=D0；B01=D1?a；B10=D0?T0；A0=D1?T；A1=D0?T+D1?T0a；A2=D0?T0a。

對維修臺工作情況進行分析，發(fā)現(xiàn)維修臺的狀態(tài)轉(zhuǎn)移過程亦為一個MAP，其中無窮小生成元矩陣為C=T+T0α；因此，存在一組穩(wěn)態(tài)概率向量πs使得：πsC=0，πse=1。

定理1矩陣A和C均是不可約隨機陣，并且矩陣A的穩(wěn)態(tài)概率向量πA=π?πs。

證明由A和C的定義已知它們均為不可約隨機陣，而：

因此矩陣A的穩(wěn)態(tài)概率向量πA=π?πs。定理得證。

由MAP和PH分布性質(zhì)知，顧客到達率和服務率分別為：λ=πD1e，u=-aT-1e。

Neuts[14，15]給出了求矩陣Q穩(wěn)態(tài)概率向量的一個重要矩陣——率陣R，R是矩陣方程R2A2+RA1+A0=0的最小非負解。

證明根據(jù)文獻[7]定理2.4知若該馬氏鏈正常返，則必須{} ψ(t)的 R的譜半徑 sp(R)＜1，即πAA2＞πAA0，該式正是經(jīng)典排隊系統(tǒng) μ＞λ條件向矩陣空間的推廣。定理得證。

在定理1、定理2的基礎上，運用矩陣幾何解理論，不難得出：

推論1若 ρ＜1，則矩陣Q的穩(wěn)態(tài)概率向量θ=(θ0，θ1，θ2，...)存在，且θk=θ1Rk-1，k＞1；而 θ0和θ1滿足下面的矩陣方程：

并且由下式歸一化后可得θ：

θ就是系統(tǒng)在任意時刻的隊長分布函數(shù)，因此容易給出任意時刻隊長的各階矩參數(shù)。根據(jù)推論1可以獲得各維修隊列、采購隊列的平均隊長及方差等參數(shù)。

此外，備件從中心倉庫向基地i運送的過程也可以看作是一個MAP/PH/1排隊系統(tǒng)，這是因為備件進入運輸渠道是一個MAP，而運輸時間可以視為服務臺服務時間，因此根據(jù)推論1容易確定運輸中的備件個數(shù)分布。

4.3缺貨量分布

記x0為中心倉庫不可用備件數(shù)，則中心倉庫缺貨量的期望值可以表示為：

在到達過程為MAP的情況下，由MAP性質(zhì)[15]可知：

令xi為基地i的不可用備件數(shù)，則基地i的期望缺貨量則為：

而BOi的兩階矩為：

在獲得各基地缺貨量的期望和方差后，Grave[2]通過大量數(shù)據(jù)分析認為缺貨量分布可以近似表示為參數(shù)a和b的負二項分布，因此本文亦直接采用該結(jié)論，由式（20）、式（21）可確定出a、b的值。

5　庫存優(yōu)化模型

定義PBO(S0，Si)為在給定S0和Si值后，基地i發(fā)生缺貨的概率，由4.3節(jié)獲得的缺貨量分布容易得到：

因為各基地是一個并聯(lián)系統(tǒng)，因此該多級庫存系統(tǒng)的備件滿足率可以表示為：

以總費用ctotal最小為目標函數(shù)，確定合適的S*=(S0，S1，...，SN)使得整個系統(tǒng)滿足備件滿足率約束。令c1、c2、c3、c4和c5分別為備件單價，單位備件庫存費用，單位備件基地維修費用，單位備件中心倉庫維修費用和單位缺貨費用，因此該優(yōu)化模型可以表示為如下形式：

式（23）中，Li和Wi分別表示各庫存點庫存量和維修隊列長。

顯然，備件滿足率與單個庫存點備件數(shù)之間存在單調(diào)遞增關(guān)系；但A(S0，S1，...，SN)與{(S0，S1，...，SN)}的關(guān)系由于涉及多個相互獨立的變量，就非常復雜了。對這類非線性整數(shù)規(guī)劃模型而言，許多成熟的智能算法均能有效獲得最優(yōu)解。本文對智能算法不做深入研究，而直接采用遺傳算法進行求解。

下面通過示例演示MAP和PH分布能夠有效描述兩級可修備件的（S-1，S）庫存模型。

6　算例

某中心倉庫負責4個基地的備件供應和維修保障任務?；豬的某型部件備件需求到達過程可以表示為MAP（D0(i)，D1(i)）；而各基地的維修時間分布、向中心倉庫運輸故障部件的時間分布可以分別表示為PH（αi，Ti）和PH（αiY，TiY），表1、表2是各基地的需求、維修及運輸數(shù)據(jù)。

中心倉庫的維修時間分布為PH(α0，T0)，采購備件的時間分布為PH(αC，TC)，其中：

而在中心倉庫可修的概率則為p=0.79。c1、c2、c3、c4和c5分別為備件單價，單位備件年庫存費用，單位備件基地維修費用，單位備件中心倉庫維修費用和單位缺貨費用，c1=26 800元，c2=90元，c3=440元，c4=1 900元，c5=7 000元。

表1　各基地備件需求到達過程Table 1　The spare parts demand process of each base

表2　各基地備件維修時間分布及運輸時間分布Table 2　The distributions of repair time and transit time for each base

下面利用4.3節(jié)獲得的缺貨量分布，研究中心倉庫和各基地備件數(shù)對整個系統(tǒng)的備件滿足率的影響。由于本優(yōu)化模型包含5個自變量，因此為了便于更直觀地了解備件數(shù)與備件滿足率之間的關(guān)系，以中心倉庫備件數(shù)S0、基地1備件數(shù)S1同時對系統(tǒng)備件滿足率的影響為例進行說明。

如圖3所示，令其他基地備件數(shù)均為0時，可以看到整個系統(tǒng)的備件滿足率早期隨著S0、S1增加而增加幅度較大；但當S0≥8并且S1≥2時，備件滿足率就保持在0.711 2附近，說明要使?jié)M足率繼續(xù)增加，必須調(diào)整其他基地的備件量。當其他基地的備件數(shù)均設為1時，系統(tǒng)備件滿足率有大幅提升，迅速達到了0.937 4附近。而當其他基地的備件數(shù)均為5時，滿足率上升則沒有早期那么明顯了，穩(wěn)定在0.999 9附近。從圖3可以看到，當給定其他基地的備件數(shù)時，S0、S1對滿足率的影響十分相似。同時也表明無論中心倉庫還是基地的備件數(shù)增加，均會使得備件滿足率增加，只是增加的幅度與系統(tǒng)內(nèi)其他倉庫的備件數(shù)量密切相關(guān)；如果只增加S0,…,S4中任意一個，整個系統(tǒng)的備件滿足率增長則十分有限。而當S0,...,S4均較大時，系統(tǒng)的備件滿足率非常接近1，但此時僅購置費就會大幅上升。

圖3　中心倉庫和基地1備件數(shù)對備件滿足率的影響Fig.3　The influence of the spare amount of central depot and base 1 on the system spare fill rate

在對備件滿足率的分析過程中，人工對各庫存點的庫存量進行調(diào)整，發(fā)現(xiàn)當(S0，...，S4)=(6，0，2，2，1)時，滿足率為0.9092；而(S0，...，S4)=(5，1，1，1，1)時，滿足率達到了0.921 5。這說明為了提高系統(tǒng)的備件滿足率，除了增加總的備件量外，還必須合理分配各單位的備件數(shù)。

7　結(jié)語

本文分析了以往多級可修備件庫存模型中存在的問題，用MAP描述各基地的故障到達流，將維修與采購時間、運輸時間均假設服從一般分布，建立了一個適應性更好的多級可修備件（S-1，S）模型，并給出了該系統(tǒng)缺貨量分布。該類模型涉及到的庫存優(yōu)化算法是庫存模型研究的一個重要方向，結(jié)合該類庫存問題的實際背景和特點，開發(fā)更合理高效的啟發(fā)式算法對本文模型在實際工作中獲得有效應用有著重要的意義。本文建立的模型在計算時涉及到矩陣的運算，而目前高性能計算機和矩陣解析方法的應用能對大型矩陣的運算提供良好的支持，所以本文建立的模型具有很好的實用價值。

[1]Sherbrooke C C.METRIC：a multi-echelon technique for recoverable item control[J].Operations Research，1968，16（1）：122-141.

[2]Graves S C.A multi-echelon inventory model for a repairable item with one-for-one replenishment[J].Management Science，1985，31：1247-1256.

[3]Díaz A，F(xiàn)u M C.Models for multi-echelon repairable item inventory systems with limited repair capacity[J].European Journal of Operational Research，1997，97：480-492.

[4]Sleptchenko A，Heijden M C，Harten A.Effects of finite repair capacity in multi-echelon，multi-indenture service part supply systems[J].International Journal of Production Economics，2002，79（3）：209-230.

[5]付興方，李繼軍，李宗植.基于兩級供應關(guān)系的可修復航材存儲策略模型研究[J].系統(tǒng)工程理論與實踐，2004，24（1）：111-115.

[6]Lau H C，Song H.Multi-echelon repairable item inventory system with limited repair capacity under nonstationary demands[J]. International Journal of Inventory Research，2008，1（1）：67-92.

[7]田乃碩，岳德權(quán).擬生滅過程與矩陣幾何解[M].北京：科學出版社，2002.

[8]Kim J S，Hur S，Kim T Y.A multiechelon repairable item inventory system with lateral transshipment and a general repair time distribution[J].Informatica，2006，17（3）：381-392.

[9]Asmussen S，Koole G.Marked point processes as limits of Markovian arrival streams[J].Journal Application of Probability，1993，30（2）：365-372.

[10]Asmussen S，Bladt M.Point processes with finite-dimensional conditional probabilities[J].Stochastic Processes and Their Applications，1999，82（1）：127-142.

[11]Chakravarthy S R.The batch Markovian arrival process：A review and future work，in advances in probability theory and stochastic processes[Z].A.Krishnamoorthy.New Jersey，2001，21-49.

[12]Huang Z，Liang L，Guo B.A general repairable spare part demand model based on quasi birth and death process[J].Acta Automacita Sinica，2006，32（2）：200-206.

[13]田乃碩.休假隨機服務系統(tǒng)[M].北京：北京大學出版社，2001.

[14]Neuts M F.Matrix-Geometric Solutions in Stochastic Models-An Algorithmic Approach[M].Baltimore：Johns Hopkins University Press，1981.

[15]Breuer L，Baum D.An Introduction to Queueing Theory and Matrix-Analytic Methods[M].Dordrecht：Springer，2005.

A two-echelon（S-1，S）inventory model for repairable items based on markovian arrival process

Chen Tong，Li Fang，Di Peng
（Department of Management Engineering，Naval University of Engineering，Wuhan 430033，China）

This paper investigates a two-echolon inventory system with（S-1，S）policy that consists of several same repairable items and single repair facility，and assumes that the item demand occur according to a markovian arrival process（MAP），the repair time，ship time and procurement time follow the general distribution which is represented by phase-type（PH）distribution.Then a inventory optimization model with better description ability and analytical performance is given，and the probability distribution of backorder is obtained.Finally，a numerical example was given to illustrate the effectiveness of the model.

（S-1，S）policy；two-echelon inventory；repairable item；markovian arrival process

F253.4

A

1009-1742（2015）05-0113-07

2015-03-15

陳童，1980年出生，男，河南駐馬店市人，博士，研究方向為裝備綜合保障，系統(tǒng)可靠性；E-mail：chentong@nudt.edu.cn