需求隨機依賴庫存環境下的訂貨仿真優化模型

歐 劍，閔 杰

OU Jian,MIN Jie

安徽建筑大學 數理系，合肥 230601

Department of Mathematics&Physics,Anhui Jianzhu University,Hefei 230601,China

1 引言

庫存控制是供應鏈管理的一個關鍵因素，庫存控制策略的優化可以極大地優化整個供應鏈的成本和利潤水平。現實環境中需求通常不是確定和靜態的，而是具有隨機性和動態性。目前關于隨機需求庫存模型的研究方法按照優化目標形式和求解方式總體上分為三類。第一類是在需求服從某種隨機分布的假設下，推導出贏利或成本的數學期望的函數表達式，然后利用解析的方式求出極值點，亦即最優訂貨策略。這類方法在對訂貨周期內受隨機需求影響的庫存成本和缺貨成本建模計算時通常采用兩種方法：（1）用某一時間點（期初或期末）的庫存量或者某幾個時間點庫存量的函數來近似整個周期的平均庫存量[1-2]，這類方法簡單易行但誤差較大。（2）應用微元法精確計算實時的庫存和缺貨成本[3]，這類方法結果精確但優化目標函數形式復雜、求解困難。第二類研究方法針對解析推導的目標函數難以求解的問題，應用進化計算、模擬退火等智能算法對目標函數進行尋優[3]。這兩類研究方法都是用數學方法推導出目標函數的解析表達式。但是，一旦隨機需求量服從的分布稍加改變，原有的解析模型就不再適用，而且，考慮更多復雜因素的訂貨模型往往難以推導出顯式的目標函數。因此，前兩類解析方法在求解效果和效率上往往不是十分令人滿意，在實際的供應鏈決策中應用起來具有較大的缺陷。

為了解決上述問題，很多研究者提出了通過系統仿真評價訂貨模型、采用智能算法優化訂貨策略的第三類方法，如胡幼華、潘蔭榮等[4]建立了隨機性的批量制造庫存系統的仿真模型，Burcu B.Keskin和Sharif H.Melouk[5]設計了隨機庫存系統的基于分散搜索的啟發式優化方法，陰艷超、郭成[6]提出了一種啟發式動態響應算法實現對隨機需求提前期環境下的庫存系統的優化，Ming-Feng Yang，Yi Lin[7]提出用非線性粒子群優化算法求解多級庫存系統的最優庫存策略。

但是，在目前關于隨機庫存的研究文獻中，均沒有考慮需求隨庫存水平動態變化這一客觀現象。自從20世紀70年代以來，許多供應鏈管理學家都發現，在超市之類的零售系統中庫存商品展示得越多，顧客的購買欲望就越強，需求量就越大。這種商品庫存促進銷售的模式稱之為“需求依賴于庫存”。關于需求受庫存水平影響的形式主要有兩大類。第一類是將需求率看成初始庫存水平的函數，Gupta和Vrat[8]首次在庫存模型中考慮并用需求函數反映了這一規律。第二類是將需求看成是當前庫存水平的函數。Mandal[9]首次將需求率假設為庫存水平的線性函數，而Baker和Urban[10]提出了需求率是庫存水平的多項式函數，該模型能夠較準確地表現需求與庫存的關系，因此在其后被其他模型廣泛采用。Karabi等人[11]又假設了需求率是分階段的，開始時的需求率是當時庫存水平的多項式函數，當庫存水平下降到一定值時需求率就變成常數。閔杰等人[12]研究了需求依賴庫存和通貨膨脹環境下考慮不同缺貨策略的庫存模型，但是以上這些模型都沒有考慮到隨機需求的情形。

本文建立了需求隨機依賴庫存的訂貨模型并通過基于仿真的進化算法進行求解。首先，假設需求是當前庫存水平的單調非減函數，即庫存水平越高，需求率相應越高。然后，用計算機仿真的方法觀察庫存模型的平均贏利水平，評價訂貨策略。最后，由于仿真方法不存在顯式的目標函數，因此將其與進化計算相結合，對訂貨策略進行優化。

2 需求隨機依賴庫存的（Q，T）型訂貨模型

2.1 模型假設與符號

（1）補貨是瞬時完成的；

（2）最大庫存水平為Q，訂貨周期為T；

（3）t時刻的庫存水平記為I（t），需求率記為D（t）；

（4）單位購買成本為c，銷售價格為p，單位時間（一天）單位物品的存儲成本為h，一次訂貨固定成本為A；

（5）允許缺貨，單位時間（一天）單位物品的缺貨成本為d；

（6）第t0天庫存降為零。

2.2 允許短缺的隨機（Q，T）訂貨策略

（Q，T）庫存控制策略是一種定期訂貨策略，即每隔T時刻訂貨一次，將庫存瞬時補充到最大庫存水平Q。如果存在缺貨銷售，則還需將累積的短缺需求補足。如果需求量是隨機的，則庫存水平的下降速度是動態隨機變化的。一個訂貨周期結束時的庫存可能存在三種情況：（1）庫存為零，即庫存Q在一個周期內恰好銷售完；（2）庫存大于零，即存在剩余庫存；（3）庫存為負，即存在由缺貨銷售累積的短缺需求。整個系統的庫存水平變化過程如圖1所示。

圖1 瞬時補貨的庫存水平隨機變化示意圖

2.3 需求是庫存函數的訂貨模型及改進

傳統的隨機庫存模型假設需求是隨機變量，而其所服從的分布在整個訂貨周期是固定不變的。但在實際環境中，某一時刻的需求率通常是受到當時庫存水平的影響，庫存越多，需求越大。1988年，Baker和Urban[10]在庫存模型中重新考慮了庫存水平對需求的影響方式，修正了之前需求是庫存的線性函數的簡單假設，而是假定需求率是瞬時庫存水平的多項式函數。這類多項式需求函數的形式為：

其中，α＞0，稱為規模參數；0＜β＜1，稱為形狀參數。式（1）表明，隨著庫存水平的降低，需求率加速減小；當庫存水平降至零庫存，需求率趨于零，即沒有需求。

但是，由于允許短缺，即在缺貨的情況下仍然可以銷售，因此可以假設零庫存時的需求率為D0＞0。顯然，庫存水平大于0時的需求率應大于D0，因此，可將需求函數修改為如下形式：

2.4 需求依賴庫存的隨機（Q，T）模型

假設需求是隨機的，且依賴庫存水平，即庫存水平越高，需求率總體上就越大，則式（2）修改為 D(t)=αI(t)β+D0+ε，其中ε為隨機波動。

因此，隨機需求的數學期望值隨著庫存水平的增高而變大。不妨假設單位時間（一天）的隨機需求量服從參數為λ的泊松分布，因為泊松分布的期望即為參數λ，所以有如下規律：

其中，λ(t)表示第t天的需求量D(t)所服從的Poisson分布的參數，λ0表示零庫存時一天的隨機需求量D0所服從的Poisson分布的參數。式（3）表明隨著銷售的進行，庫存I不斷降低，隨機需求量的均值λ隨著庫存I的降低而減小，直至庫存降為零后，在缺貨期間隨機需求的均值恒定為λ0。不妨假設第t0天庫存降為零。

在上述假設下，（Q，T）型訂貨模型的單位時間（一天）內平均利潤的數學期望ATP具有如下計算形式：

本文的目標是確定最優的訂貨周期T和最大庫存水平Q，使該系統單位時間（一天）內的平均利潤的數學期望達到最大。

3 離散事件系統仿真

3.1 計算機仿真

近年來，人們所研究的系統越來越復雜，系統構成因素復雜且相互關聯，因此很難對系統建立解析模型，即使建立了數學模型也難以求解。而計算機仿真的方法由于不需要進行數學解析，可以考慮各種復雜因素，隨意地調整系統結構和參數的隨機分布，從而能更貼切地、動態地反映出系統運行的實際特征。因此，通過計算機仿真結合統計的方法可以對系統性能進行有效的評價。

在本文研究的訂貨模型中，由于隨機需求D(t)依賴庫存水平 I(t)，而I(t)=I(t-1)-D(t-1)，即t時刻之前的需求影響到t時刻的庫存和需求，因此隨機變量D(t)和I(t)(t=1，2，…，T)相互不獨立，所以訂貨模型的平均利潤的期望式（4）很難以解析形式表示出來。因此，用計算機仿真的方法對（Q，T）型訂貨策略下的平均利潤水平進行評價是一種合適的手段。

3.2 訂貨模型的離散事件系統仿真

很多實際的庫存系統中，需求量是離散的隨機變量，庫存水平不是連續變化的，而是在離散時間點上改變的。因此可用離散事件系統仿真的方法模擬隨機庫存系統的運行。離散事件系統仿真有事件步長法和時間步長法兩種仿真策略[13]。時間步長法就是按照時間流逝的順序，一步一步地對系統的活動進行仿真，在整個仿真過程中時間步長保持不變。仿真時鐘每步進一次，就對系統中所有的實體、屬性和活動進行一次全面的掃描，按照預定的目標記錄、計算系統狀態的變化，直至仿真時鐘結束為止。

本文研究的（Q，T）訂貨模型假設需求是離散型隨機變量，且訂貨周期是固定的，因此可以建立時間步長法的離散系統仿真模型，以一天為時間步長進行仿真。首先檢查當天是否為訂貨日期，如果是則將庫存量刷新為最大庫存Q，累計上一次訂貨費；如果不是，則庫存量不變。接著按照式（3）確定隨機需求量服從的隨機分布的參數，產生當天的需求值，計算當天的收益。若庫存量大于需求量，則將庫存量減去需求量，計算當天的存貯費用；反之，則累加上缺貨量，計算當天的缺貨費用。當仿真完成N（充分大）個周期后，用累計的收益減去累計的訂貨費、存貯費、缺貨費，得到凈利潤，然后除以仿真天數得到平均利潤的仿真值。這一仿真過程的流程如圖2所示。

圖2 （Q，T）訂貨模型的仿真流程框圖

4 基于仿真的進化計算

4.1 進化計算

進化計算[14]（Evolutionary Algorithm，EA）是一種模仿生物進化機制的智能優化算法，它模擬生物群體的自然選擇、雜交、突變等進化形式，采用選擇（select）、交叉（crossover）、變異（mutation）三種遺傳算子，對編碼的解群體進行進化操作，指導種群向最優解方向搜索。進化算法具有并行計算、自組織、自適應、自學習等良好的特性，更重要的是，進化計算不會受到搜索空間限制性條件（如連續、可導、單峰）的約束，不需要其他輔助信息（如梯度、海塞矩陣）。因此，進化計算不僅尋優效率較高，而且相較傳統優化算法更加簡單、易行、通用。

本文所研究的需求依賴庫存的（Q，T）訂貨模型，其目標函數是平均利潤的數學期望，由于期望的解析式難以顯式表達，因此采用時間步長法的離散事件系統仿真得到在某一具體的（Q，T）參數下平均利潤期望的近似值。所以，既然沒有顯式的目標函數表達式，那么，傳統的依賴梯度的優化方法就不再適用，而基于仿真的進化計算[15]則不會受到目標函數的局限，可以很好地解決該模型的優化問題。

4.2 編碼和適應度函數

需求隨機依賴庫存的（Q，T）訂貨模型的決策變量是最大庫存Q和訂貨周期T。在進化計算中采用實數編碼方案，取控制變量（Q，T）作為進化個體的編碼。

適用度函數表明解個體的優劣性，一般是由目標函數轉化而成。因為（Q，T）訂貨模型優化的目標是平均利潤的期望最大化，所以，將個體（Q，T）作為系統仿真控制參數，運行離散事件系統仿真模型得到平均利潤期望的仿真值，以該仿真值直接作為個體（Q，T）的適用度。

4.3 選擇算子

選擇操作的目的是從當前群體中以較大概率選出適用度較高的優良個體，進行交叉、變異產生新個體，個體的適用度越高，選擇概率就應越大。輪盤賭選擇法是一種常見的選擇算法，但是當種群中個體適應度（仿真結果）差不多時，不同個體的選擇概率相差不大，此時輪盤賭選擇就近似于隨機選擇，起不到優勝劣汰的效果。因此，本文采用Michalewicz的線性排序的選擇概率計算公式：

其中i表示個體按適應度從大到小排序的序號。這種選擇概率計算公式使得個體之間無論適應度多么接近，都會按優劣順序拉開被選中概率的差距。q越大，不同個體的選擇概率差距就越大，排名靠前的個體被選中的概率就越大，因此，q不應過大，以免較優的少數個體充斥整個種群，造成“早熟”。

4.4 交叉算子與變異算子

交叉操作采用中間重組法，子個體按照如下公式產生：

當s∈(0，1)時，子個體是一對父體的凸組合。如果父代種群是凸集，那么通過交叉產生的子代個體都是這個凸集的內點。

變異操作采用高斯變異算子，變異公式為：

σ2為高斯分布方差，σ2越大，則變異幅度越大。當進化出現“早熟”，種群陷入局部最優解，可通過加大變異幅度使個體跳出局部最優解。

4.5 基于仿真的進化計算結構

進化計算的種群可分為重疊和非重疊兩種組織形式，遺傳算子可分為重疊和非重疊兩種執行方式。本文針對（Q，T）訂貨模型的特性，設計種群重疊、遺傳操作非重疊的進化計算結構如下：

5 基于仿真的模擬退火與粒子群優化

與基于仿真的進化算法類似，基于仿真的模擬退火（Simulated Annealing，SA）和粒子群優化（Particle Swarm Optimization，PSO）等其他基于仿真的優化算法（Simulation Based Optimization，SBO）也可用于隨機（Q，T）訂貨模型的優化。

基于仿真的SA算法如下：

步驟1設定一個初始狀態 X=(Q，T)，較高的起始溫度t。

步驟2以（Q，T）作為控制參數運行離散事件系統仿真模型，得到平均利潤期望的仿真值，取其負值作為X的能量E。

步驟3根據式（7）得到新狀態Xnew，并計算ΔE=E(Xnew)-E(X)。

步驟4如果ΔE＜0，則置 X=Xnew；否則，以概率exp(-ΔE/t)置 X=Xnew 。

步驟5 降低溫度，t=η×t，0＜η＜1。

步驟6若達到最大迭代次數，則停止，輸出X；否則轉步驟3。

基于仿真的PSO算法如下：

步驟1初始化N個粒子，每個粒子位置即為一個庫存控制變量 X=(Q，T)，在可行域內隨機設定初始值和初始速度，k=1，2，…，N 。

步驟2將種群中每個粒子的位置（Q，T）輸入離散事件系統仿真模型，得到平均利潤期望的仿真值作為粒子的當前評價值。

步驟3更新粒子的當前最好位置Pk和全局最好位置Gbest。

步驟4根據式（8）和式（9）改變粒子的速度和位置。

步驟5如果達到最大迭代次數，則停止，輸出Gbest；否則轉到步驟2。

6 仿真實例及其分析

6.1 進化計算的實例分析

假設單位貨物售價 p=10，購買成本c=5，單位貨物單位時間存貯成本h=0.6，缺貨成本d=0.7，一次訂貨成本 A=80，規模參數 α=1.5，形狀參數 β=0.4，單位時間需求量 D～Poisson(λ)，缺貨時的需求率λ0=20。

設置仿真長度=100個周期，選擇概率系數q=0.1，交叉概率=0.5，變異概率=0.1，交叉點s=0.618，高斯變異標準差σ=1，種群規模N=20，進化代數=20，得到庫存控制策略的進化曲線如圖3所示。由圖3易見，從第1代到第3代種群的整體適應度迅速提升，說明進化算法在進化初期的收斂速度很快，種群迅速向最優個體集中。當進化到第14代，最優適應度不再提高，種群整體適應度則始終緩慢提升。

圖3 進化過程中適應度曲線圖

進化完成后得到平均利潤最優的5個訂貨策略如表1所示。從表1中可以看到，如果要提高服務水平，降低缺貨率，則必須提高最大庫存Q值，或縮短訂貨周期T，從而導致平均利潤降低。因此，在實際的需求隨機依賴庫存的（Q，T）訂貨模型中要綜合考慮盈利水平和服務水平，作出最佳的訂貨策略。

為考察模型參數對最優解的影響，通過分別改變參數λ0、α、β的數值來進行模型靈敏度分析，結果列于表2。

表1 仿真結果（每天平均利潤）

表2 模型的靈敏度分析

由表2易得如下結論：（1）缺貨期平均需求率λ0、規模參數α和形狀參數β越大，則最優的平均利潤水平越高，最優的最大庫存水平Q越大，最優的訂貨周期T越小，缺貨率越低，這說明需求率越高或需求依賴庫存效應越明顯時，利潤水平越高，最優訂貨策略越傾向于高庫存、短周期、現貨銷售。（2）最優訂貨周期T對α和 β的改變不是非常敏感，對λ0較敏感。（3）最大庫存水平Q對 β的改變較λ0和α更敏感。

6.2 不同SBO算法對比分析

用模擬退火算法優化該實例，取初始退火溫度為1000，冷卻系數η=0.95，得到能量下降曲線如圖4所示。退火算法在迭代到495次時收斂到最小值－52.8592。

用粒子群算法優化該實例，設置粒子群規模N=20，最大迭代次數=20，慣性權重w=0.3，學習因子c1=2，c2=1.5。當迭代到11次時平均利潤收斂到50.95。

圖4 模擬退火的能量下降曲線

由于SBO算法的計算量主要是仿真的計算量，所以可用仿真的次數來度量算法的計算量。因為EA、PSO和SA算法在每次產生新解時進行一次仿真評價，所以算法的計算量為種群規模和迭代次數的乘積。三種算法的性能對比如表3所示。由表3可見，EA和PSO算法的收斂速度較快，而SA算法為了避免陷入局部最優則需迭代較多次數。因此，針對本文的隨機庫存優化模型，基于仿真的EA和PSO算法均為效率較高的優化方法。

表3 三種算法的計算量對比表

7 結束語

本文建立的需求隨機依賴庫存的訂貨模型反映了實際庫存系統中需求的隨機性和動態性。由于該模型中平均利潤的數學期望難以解析表示，所以無法通過傳統的優化方法求得最優訂貨策略，而基于離散事件系統仿真的進化算法則可以較好地解決這一優化問題，該算法不僅具有較高的運算效率和靈活性、通用性，而且可以獲得多個候選策略以供決策者綜合權衡選擇。但本文的進化計算僅將單位時間的平均利潤最大化作為進化目標，即只把平均利潤作為適應度函數，而實際決策時不僅要考慮利潤最大化，還要考慮服務質量，即缺貨率盡可能低，因此如何設計高利潤和低缺貨率的多目標優化算法是下一步的研究目標。

[1]Feng Cheng，Sethi S P.A periodic review inventory model with demand influenced by promotion decisions[J].Management Science，1999，45（11）：1510-1523.

[2]Bookbinder J H，Cakanyildirim M.Random lead times and expedited orders in（Q，r）inventory systems[J].European Journal of Operational Research，1999，115（2）：300-313.

[3]紀鵬程，宋士吉，吳澄.隨機需求和提前期環境下的精確庫存成本建模[J].計算機集成制造系統，2008，14（11）：2129-2133.

[4]胡幼華，潘蔭榮，房愛蓮.隨機性的批量制造庫存系統的仿真模型[J].計算機工程，2003，29（8）：70-71.

[5]Keskin B B，Melouk S H，Meye I L.A simulation-optimization approach for integrated sourcing and inventory decisions[J].Computers&Operations Research，2010，37（9）：1648-1661.

[6]陰艷超，郭成.不確定環境下的一種動態響應庫存優化方法[J].計算機工程與應用，2012，48（17）：238-242.

[7]Yang Mingfeng，Lin Yi.Applying the linear particle swarm optimization to a serial multi-echelon inventory model[J].Expert Systems with Applications，2010，37（3）：2599-2608.

[8]Gupta R，Vrat P.Inventory model for stock-dependent consumption rate[J].Opsearch，1986，23（1）：19-24.

[9]Mandal B N，Phaujdar S.An inventory model for deteriorating item and stock-dependent consumption rate[J].Journal of the Operational Research Society，1989，40（5）：483-488.

[10]Baker R C，Urban T L.A deterministic inventory system with an inventory-level-dependent demand rate[J].Journal of the Operational Research Society，1988，39（9）：823-831.

[11]Karabi P，Datta T K，Chaudhuri K S，et al.An inventory model with two-component demand rate and shortages[J].Journal of the Operational Research Society，1996，47：1020-1036.

[12]閔杰，歐劍，周永務.通貨膨脹環境下需求依賴庫存的經濟訂貨批量模型[J].計算機集成制造系統，2010，16（2）：299-309.

[13]肖田元，張燕云，陳加棟.系統仿真導論[M].北京：清華大學出版社，2000.

[14]Jones M T.Artificial intelligence[M].北京：電子工業出版社，2009：116-117.

[15]李東，汪定偉.基于仿真的優化方法綜述[J].控制工程，2008，15（6）：674-677.