基于馬爾科夫鏈-多目標(biāo)模型的應(yīng)急供應(yīng)鏈決策優(yōu)化研究*

陳偉炯，董雯玉，李咪靜，張善杰，李曉戀,，康與濤

(1.上海海事大學(xué) 物流科學(xué)與工程研究院,上海 201306；2.上海海事大學(xué) 物流供應(yīng)鏈風(fēng)險(xiǎn)控制研究中心，上海 201306；3.上海海事大學(xué) 海洋科學(xué)與工程學(xué)院，上海 201306)

0 引言

突發(fā)災(zāi)難事件頻發(fā)，嚴(yán)重影響人類生活、社會(huì)穩(wěn)定和經(jīng)濟(jì)發(fā)展。國(guó)家應(yīng)急管理部數(shù)據(jù)顯示，2020年我國(guó)自然災(zāi)害以地質(zhì)災(zāi)害、洪澇為主，全年受災(zāi)人群達(dá)1.38億人，死亡及失蹤人數(shù)共591人,造成直接經(jīng)濟(jì)損失3 701.5億元[1]。2020年新冠疫情曾造成口罩、防護(hù)服等相關(guān)應(yīng)急醫(yī)療物資短缺，嚴(yán)重威脅公眾生命安全。應(yīng)急救援亟待低風(fēng)險(xiǎn)、高效率的應(yīng)急供應(yīng)鏈。

目前，應(yīng)急供應(yīng)鏈風(fēng)險(xiǎn)分析與應(yīng)急物資優(yōu)化調(diào)度成為國(guó)內(nèi)外研究熱點(diǎn)。在應(yīng)急供應(yīng)鏈風(fēng)險(xiǎn)分析方面，研究趨勢(shì)逐漸由靜態(tài)向動(dòng)態(tài)轉(zhuǎn)變，動(dòng)態(tài)研究能夠更加真實(shí)、客觀地反應(yīng)供應(yīng)鏈風(fēng)險(xiǎn)及其變化。一方面,學(xué)者基于馬爾可夫過程、混合數(shù)學(xué)規(guī)劃、動(dòng)態(tài)貝葉斯網(wǎng)絡(luò)、魯棒控制優(yōu)化等[2-4]方法構(gòu)建動(dòng)態(tài)應(yīng)急供應(yīng)鏈模型，在災(zāi)情更新情況下分析供應(yīng)鏈中的潛在風(fēng)險(xiǎn)，為應(yīng)急管理提供實(shí)時(shí)有效的政策調(diào)整建議。另一方面，單一目標(biāo)的應(yīng)急物資調(diào)度規(guī)劃模型已無法滿足應(yīng)急管理實(shí)際需求，因此,部分學(xué)者聚焦于多目標(biāo)規(guī)劃模型的研究，包括以最小化延遲時(shí)間、最小化運(yùn)輸成本或系統(tǒng)損失的雙目標(biāo)規(guī)劃模型[5-6]，在此基礎(chǔ)上考慮最大化線路可靠性[7]、最大化供應(yīng)鏈可視性[8]、最小化物資短缺[9]等因素的3目標(biāo)優(yōu)化模型。為有效提高應(yīng)急供應(yīng)鏈中斷響應(yīng)效率，學(xué)者展開應(yīng)急供應(yīng)鏈決策優(yōu)化研究，基于可持續(xù)性、魯棒性、韌性等，設(shè)計(jì)優(yōu)化應(yīng)急組織配置策略[10]、協(xié)同優(yōu)化方法[11]、韌性優(yōu)化策略[12]，增強(qiáng)應(yīng)急供應(yīng)鏈恢復(fù)力，高效地應(yīng)對(duì)風(fēng)險(xiǎn)，降低災(zāi)害損失。

在應(yīng)急供應(yīng)鏈風(fēng)險(xiǎn)動(dòng)態(tài)變化基礎(chǔ)上，現(xiàn)有研究對(duì)應(yīng)急物資的多目標(biāo)規(guī)劃與決策優(yōu)化的綜合性研究尚有不足，且在算法求解上仍有改進(jìn)空間。因此，本文提出離散時(shí)間馬爾科夫鏈—多目標(biāo)規(guī)劃模型(DTMC-MOP)，綜合考慮供應(yīng)率、時(shí)間、成本3個(gè)因素，動(dòng)態(tài)識(shí)別、分析、應(yīng)對(duì)應(yīng)急供應(yīng)鏈風(fēng)險(xiǎn)，并采用改進(jìn)自適應(yīng)NSGA-Ⅱ算法優(yōu)化求解模型，以實(shí)現(xiàn)應(yīng)急物資的最大滿足率供應(yīng)、最短時(shí)間供應(yīng)和最低成本供應(yīng)目標(biāo)，研究結(jié)果可為及時(shí)控制災(zāi)情和最大限度地降低災(zāi)害損失提供參考。

1 DTMC-MOP應(yīng)急供應(yīng)鏈模型的構(gòu)建

突發(fā)災(zāi)難事件發(fā)生時(shí)，應(yīng)急供應(yīng)鏈易受內(nèi)外部不同程度的風(fēng)險(xiǎn)擾動(dòng)。為保障應(yīng)急救援工作的順利開展，構(gòu)建由應(yīng)急管理部門、供應(yīng)商、制造商、需求方組成的4級(jí)動(dòng)態(tài)應(yīng)急供應(yīng)鏈網(wǎng)絡(luò)，基于離散時(shí)間馬爾科夫鏈(DTMC)模型，將狀態(tài)轉(zhuǎn)移概率及平穩(wěn)概率分布引入多目標(biāo)規(guī)劃(MOP)模型，構(gòu)建供應(yīng)物資滿足率最大、供應(yīng)時(shí)間最短、供應(yīng)成本最低的DTMC-MOP模型。

1.1 基于離散時(shí)間馬爾科夫鏈的應(yīng)急供應(yīng)鏈模型

1)定義

①離散時(shí)間馬爾科夫鏈

若隨機(jī)序列{Xn,n=0,1,2,…}為狀態(tài)離散的隨機(jī)過程，其中狀態(tài)空間I={i0,i1,…,in,j}，時(shí)間集合T={0,1,2,…}，如果對(duì)于任意狀態(tài)滿足式(1)：

(1)

則稱{Xn,n=0,1,2,…}是離散時(shí)間馬爾科夫鏈[13]。

②k步轉(zhuǎn)移概率及矩陣如式(2)所示:

(2)

(3)

③平穩(wěn)概率分布

若分布概率π=(π1,π2,…,πj)滿足式(4)：

(4)

則稱分布概率π是{Xn}的唯一平穩(wěn)概率分布。

2)DTMC應(yīng)急供應(yīng)鏈模型

突發(fā)災(zāi)難事件是1個(gè)隨機(jī)過程，且應(yīng)急供應(yīng)鏈的下一種狀態(tài)僅和當(dāng)前狀態(tài)有關(guān)，具有馬爾科夫性。通過DTMC模型動(dòng)態(tài)描述應(yīng)急供應(yīng)鏈所處狀態(tài)，直觀反映應(yīng)急供應(yīng)鏈的運(yùn)行狀況。

考慮到供應(yīng)鏈的韌性特征，即受到風(fēng)險(xiǎn)擾動(dòng)導(dǎo)致供應(yīng)鏈部分失效時(shí)仍能維持供應(yīng)狀態(tài)，并以最快的速度恢復(fù)到正常供應(yīng)狀態(tài)的能力。假設(shè)應(yīng)急供應(yīng)鏈在t時(shí)刻處于受不同程度干擾的有限個(gè)離散狀態(tài)，由于供應(yīng)鏈具有韌性，在t+1時(shí)刻可以向其他狀態(tài)轉(zhuǎn)移或保持該狀態(tài)不變。當(dāng)發(fā)生風(fēng)險(xiǎn)擾動(dòng)時(shí)，假設(shè)應(yīng)急供應(yīng)鏈存在以下4種狀態(tài)：狀態(tài)0即完全吸收干擾的正常供應(yīng)狀態(tài)；狀態(tài)1～3分別表示受到30%，60%，90%干擾的供應(yīng)狀態(tài)。任意2個(gè)狀態(tài)之間可相互轉(zhuǎn)移，用α1～α6和β1～β6(αi,βi≤1,i=1,…,6)表示狀態(tài)轉(zhuǎn)移概率。應(yīng)急供應(yīng)鏈DTMC模型如圖1所示。

圖1 應(yīng)急供應(yīng)鏈的DTMC模型

由圖1可知，若t時(shí)刻供應(yīng)鏈處于狀態(tài)1，則在t+1時(shí)刻分別有β1，α4，α5的概率向狀態(tài)0，2，3轉(zhuǎn)移，有1-β1-α4-α5概率保持狀態(tài)1；在t+2時(shí)刻，可從t+1時(shí)刻的狀態(tài)轉(zhuǎn)移至其他狀態(tài)。根據(jù)式(2)～(3)，該模型的狀態(tài)轉(zhuǎn)移概率矩陣Pij如式(5)所示：

(5)

通過了解應(yīng)急供應(yīng)鏈中每個(gè)狀態(tài)的轉(zhuǎn)移概率，能夠動(dòng)態(tài)識(shí)別、分析、應(yīng)對(duì)風(fēng)險(xiǎn)。當(dāng)供應(yīng)鏈狀態(tài)逐漸向中斷狀態(tài)(狀態(tài)3)轉(zhuǎn)移時(shí)，能夠及時(shí)提醒供應(yīng)鏈節(jié)點(diǎn)企業(yè)制定或修改應(yīng)急措施，以減輕干擾風(fēng)險(xiǎn)；當(dāng)供應(yīng)鏈狀態(tài)逐漸向正常供應(yīng)狀態(tài)(狀態(tài)0)轉(zhuǎn)移時(shí)，能夠體現(xiàn)出供應(yīng)鏈的韌性特征，并驗(yàn)證應(yīng)急措施的有效性。

1.2 基于DTMC-MOP的應(yīng)急供應(yīng)鏈模型

突發(fā)災(zāi)難事件的發(fā)生具有不確定性和突發(fā)性。當(dāng)發(fā)生突發(fā)事件時(shí)，災(zāi)區(qū)的應(yīng)急物資需求量會(huì)急劇增加，在此之前相關(guān)部門應(yīng)急物資的儲(chǔ)備量遠(yuǎn)小于需求量，因此，需要將應(yīng)急生產(chǎn)納入考慮范圍。應(yīng)急供應(yīng)鏈具有不確定性、強(qiáng)時(shí)效性和弱經(jīng)濟(jì)性等特征，為制定高效的應(yīng)急物資供應(yīng)方案，本文提出DTMC-MOP模型，研究供應(yīng)率、時(shí)間和成本3方面因素對(duì)應(yīng)急供應(yīng)鏈的影響。

1)模型假設(shè)

由于應(yīng)急供應(yīng)鏈涉及多個(gè)供應(yīng)環(huán)節(jié)及多個(gè)節(jié)點(diǎn)企業(yè)，為簡(jiǎn)化模型且不失一般性，本文提出以下3點(diǎn)假設(shè)：

假設(shè)1：應(yīng)急供應(yīng)鏈在0時(shí)刻處于狀態(tài)0，在t時(shí)刻進(jìn)行第1步狀態(tài)轉(zhuǎn)移，在t+1時(shí)刻進(jìn)行第2步狀態(tài)轉(zhuǎn)移，直到得到平穩(wěn)概率分布，結(jié)束狀態(tài)轉(zhuǎn)移。

假設(shè)2：制造商沒有原材料和產(chǎn)品庫(kù)存，應(yīng)急物資是收到供應(yīng)商提供的原材料后進(jìn)行生產(chǎn)制造，而后通過應(yīng)急物流送往需求地。

假設(shè)3：模型僅考慮原材料、應(yīng)急物資的生產(chǎn)時(shí)間和運(yùn)輸時(shí)間，其他耗時(shí)忽略不計(jì)。

2)DTMC-MOP模型

I為供應(yīng)鏈狀態(tài)i的集合，i∈I(i=0,1,2,3)；T為時(shí)刻t的集合，t∈T(t=0,t+1,…)；D為需求地d的集合，d∈D(d=1,2,…)；S為供應(yīng)商s的集合，s∈S(s=1,2,…)；M為制造商m的集合，m∈M(m=1,2,…)。

(6)

(7)

(8)

(9)

(10)

(11)

(12)

在DTMC-MOP模型中，式(6)描述最大供應(yīng)物資滿足率目標(biāo)，即最小化物資需求未滿足率；式(7)描述最短時(shí)間供應(yīng)目標(biāo)，即最小化原材料運(yùn)輸時(shí)間、應(yīng)急物資生產(chǎn)、運(yùn)輸時(shí)間的總和；式(8)描述最低成本供應(yīng)目標(biāo)，即最小化供應(yīng)商原材料成本和運(yùn)輸成本以及制造商物資生產(chǎn)和運(yùn)輸成本；約束(9)表示制造商M向需求地D供應(yīng)物資的數(shù)量不超過需求地的物資需求數(shù)量；約束(10)是需求地D的物資需求滿足率公式，它等于物資的實(shí)際供應(yīng)數(shù)量除以需求數(shù)量；約束(11)是離散時(shí)間馬爾科夫鏈中平穩(wěn)概率分布公式；約束(12)表示決策變量為正整數(shù)。

2 改進(jìn)NSGA-Ⅱ算法模型及評(píng)價(jià)

2.1 改進(jìn)NSGA-Ⅱ算法模型

對(duì)于多目標(biāo)規(guī)劃問題的求解，國(guó)內(nèi)外學(xué)者進(jìn)行大量研究，比較常見的求解方法包括將多目標(biāo)轉(zhuǎn)化為單目標(biāo)的精確算法[14-15]以及粒子群算法、模擬退火算法、遺傳算法等智能算法[16-18]。與前人研究相比，本文將采用收斂速度快、計(jì)算精確度高、計(jì)算復(fù)雜度低的NSGA-Ⅱ算法求解模型，但該算法在搜素性能、種群多樣性方面仍存在不足，為避免陷入局部最優(yōu)，本文通過改進(jìn)傳統(tǒng)NSGA-Ⅱ算法，進(jìn)一步適應(yīng)應(yīng)急供應(yīng)鏈風(fēng)險(xiǎn)快速、準(zhǔn)確響應(yīng)并決策的特性，提升算法的收斂速度和搜索能力。改進(jìn)自適應(yīng)NSGA-Ⅱ算法流程示意如圖2所示。

圖2 改進(jìn)自適應(yīng)NSGA-Ⅱ算法流程

1)初始化種群改進(jìn)

在應(yīng)急供應(yīng)鏈實(shí)際應(yīng)用中，初始數(shù)據(jù)較大且關(guān)系構(gòu)造復(fù)雜，本文通過以下3個(gè)步驟對(duì)數(shù)據(jù)及編碼方式進(jìn)行改進(jìn)處理：

①采用實(shí)數(shù)編碼方式對(duì)初始自變量進(jìn)行處理。

②運(yùn)用反編譯將數(shù)據(jù)編碼范圍限定為0～1的隨機(jī)數(shù)乘以基數(shù)的形式，可以大幅提高算法的運(yùn)行速度。

③對(duì)反編譯后的數(shù)據(jù)進(jìn)行歸一化處理，構(gòu)造種群規(guī)模N的初始種群E0。

2)自適應(yīng)遺傳算子改進(jìn)

①自適應(yīng)多點(diǎn)交叉

在錦標(biāo)賽選擇中得到N個(gè)新個(gè)體中隨機(jī)選擇2個(gè)個(gè)體進(jìn)行多點(diǎn)交叉，算法中交叉點(diǎn)的位置和數(shù)量是隨機(jī)產(chǎn)生的。為避免陷入局部最優(yōu)，本文采用自適應(yīng)交叉法進(jìn)行調(diào)節(jié)，從而大幅度提升算法的全局搜索能力。交叉概率Pc如式(13)所示：

(13)

式中：pfmax是該種群中適應(yīng)度最大值；pf′為2個(gè)個(gè)體中適應(yīng)度較大的值；pfmean為種群中平均適應(yīng)度。

②自適應(yīng)變異

隨機(jī)選擇交叉后的個(gè)體及基因位置以Pm的概率進(jìn)行變異，如式(14)所示：

(14)

自適應(yīng)變異算子0

3)自適應(yīng)精英保留策略改進(jìn)

將產(chǎn)生的子代種群Gn按式(15)隨機(jī)挑選出新子代精英種群NGn，同時(shí)為防止父代精英個(gè)體基因遺失，選擇父代種群En中Pareto等級(jí)低、擁擠度大的個(gè)體合并生成新種群Ln，直到Ln的規(guī)模達(dá)到N為止。該方法與傳統(tǒng)的父、子代種群合并生成新種群相比，前期精英規(guī)模較小，可以豐富種群的多樣性；后期隨精英規(guī)模增加，不斷提高種群的收斂性，可避免產(chǎn)生極端解，實(shí)現(xiàn)全局最優(yōu)，如式(15)所示：

(15)

式中：an為第n代精英保留規(guī)模的影響因子。

2.2 模型評(píng)價(jià)指標(biāo)

應(yīng)用多目標(biāo)化算法解決實(shí)際問題時(shí)，主要從分布性及收斂性2方面進(jìn)行分析。因此，本文將從基于歐氏距離的收斂性[19]和分布性指標(biāo)[20]進(jìn)行評(píng)估。

1)世代距離的收斂性指標(biāo)

收斂指標(biāo)ca可以評(píng)估算法理想的Pareto最優(yōu)前沿程度，定義如式(16)所示：

(16)

式中：Nd為算法求出的非支配向量個(gè)數(shù)；di為算法已知個(gè)體邊界與真正Pareto的歐幾里得距離。

2)空間評(píng)價(jià)的分布性指標(biāo)

指標(biāo)sa用于評(píng)價(jià)多目標(biāo)優(yōu)化算法所求的解集的分布性，如式(17)所示：

(17)

2.3 模型測(cè)試結(jié)果與分析

本文通過標(biāo)準(zhǔn)測(cè)試函數(shù)ZDT3和DTLZ2對(duì)該模型進(jìn)行測(cè)試與評(píng)估，并與傳統(tǒng)NSGA-Ⅱ算法對(duì)比，2種測(cè)試函數(shù)最優(yōu)前沿對(duì)比如圖3～4所示。

圖3 ZDT3測(cè)試函數(shù)最優(yōu)前沿比較

由圖3可知，改進(jìn)個(gè)體與參考解集基本重合，而傳統(tǒng)算法解集個(gè)體則與參考值存在偏差。由圖4可知，改進(jìn)算法求得的個(gè)體緊密附著于Pareto前沿，而未改進(jìn)算法求解的個(gè)體大多漂浮于Pareto前沿，甚至有少數(shù)個(gè)體脫離。由此，改進(jìn)的NSGA-Ⅱ算法更接近最優(yōu)Pareto前沿，算法精度相對(duì)更高。

圖4 DTLZ2測(cè)試函數(shù)最優(yōu)前沿比較

3 算例分析

假設(shè)我國(guó)某地區(qū)發(fā)生突發(fā)災(zāi)難事件，應(yīng)急管理部門需要制定應(yīng)急預(yù)案，并根據(jù)實(shí)際情況及時(shí)對(duì)應(yīng)急方案進(jìn)行調(diào)整。根據(jù)該地區(qū)目前災(zāi)情和以往發(fā)生災(zāi)害情況，該地區(qū)有2個(gè)需求地D1,D2，需要應(yīng)急生產(chǎn)和運(yùn)輸某應(yīng)急物資共35 000件；由2個(gè)供應(yīng)商S1,S2提供原材料，1 kg原材料可生產(chǎn)1件應(yīng)急物資，其中S1最多可供應(yīng)19 000 kg，S2最多可供應(yīng)17 000 kg；3個(gè)制造商M1,M2,M3負(fù)責(zé)應(yīng)急物資的生產(chǎn)和運(yùn)輸。本文運(yùn)用DTMC-MOP模型模擬應(yīng)急生產(chǎn)與調(diào)度情況，制定最佳應(yīng)急方案。該模型具體參數(shù)設(shè)置如下：

應(yīng)急供應(yīng)鏈初始狀態(tài)概率分布P0及狀態(tài)轉(zhuǎn)移概率矩陣Pij如式(18)所示：

(18)

本文運(yùn)用MATLAB R2020b軟件編程求解模型。首先依據(jù)式(4)、初始概率分布P0和狀態(tài)轉(zhuǎn)移概率矩陣Pij經(jīng)過200次迭代實(shí)現(xiàn)應(yīng)急供應(yīng)鏈狀態(tài)的動(dòng)態(tài)轉(zhuǎn)移，求得唯一平穩(wěn)概率分布π=(0.461 0 0.294 8 0.169 8 0.074 4)，將表1～5的參數(shù)依據(jù)平穩(wěn)概率分布進(jìn)行加權(quán)處理，并輸入如下參數(shù)：Cs1=6.3，Cs2=6.5，元；Cm1=18，Cm2=19.5，Cm3=19，元；k1=89，k2=92，k3=90，%；Q1=20 000，Q2=15 000，件。其次，對(duì)于改進(jìn)的NSGA-Ⅱ算法基礎(chǔ)參數(shù)設(shè)置如下：種群規(guī)模N=200，自適應(yīng)交叉概率Pc在0.2～0.8之間，變異概率Pm在0～0.1之間，最大迭代次數(shù)MG=200，運(yùn)行2種算法程序得到Pareto最優(yōu)前沿，如圖5所示。

表參數(shù)設(shè)置

表參數(shù)設(shè)置

目標(biāo)函數(shù)Z1供應(yīng)物資需求未滿足率大小、Z2供應(yīng)時(shí)間長(zhǎng)短以及Z3供應(yīng)成本高低3者之間的關(guān)系。由圖5可知，最優(yōu)解在三維空間中形成1個(gè)分布均勻的曲面，能夠較好地收斂于Pareto最優(yōu)前沿，改進(jìn)自適應(yīng)NSGA-Ⅱ算法曲線分布在傳統(tǒng)算法曲線下方，獲得更貼近理想值的Pareto前沿，目標(biāo)結(jié)果更優(yōu)。

表參數(shù)設(shè)置

表參數(shù)設(shè)置

表5 ti參數(shù)設(shè)置

圖5 Pareto最優(yōu)前沿

從改進(jìn)自適應(yīng)NSGA-Ⅱ算法的Pareto最優(yōu)前沿中分別選取單個(gè)目標(biāo)函數(shù)值最小的3組具有代表性的解，見表6。根據(jù)決策者的偏好差異可以選擇不同的解決方案，若決策者傾向于物資未滿足率最低，最佳方案是方案1；若決策者追求物資供應(yīng)時(shí)間最短，應(yīng)該選擇方案2；若決策者希望實(shí)現(xiàn)最低的物資供應(yīng)成本，方案3是最好的選擇；若決策者想綜合考慮3方面的因素，可以在其他Pareto最優(yōu)前沿中選擇相適應(yīng)的應(yīng)急方案。

表6 Pareto典型解

4 結(jié)論

1)提出DTMC-MOP模型，制定有效的應(yīng)急物資生產(chǎn)與調(diào)度方案?；贒TMC模型中的狀態(tài)轉(zhuǎn)移概率矩陣刻畫應(yīng)急供應(yīng)鏈的動(dòng)態(tài)性特征，可用于識(shí)別、分析、應(yīng)對(duì)應(yīng)急供應(yīng)鏈風(fēng)險(xiǎn)。在此基礎(chǔ)上結(jié)合MOP模型，構(gòu)建供應(yīng)物資滿足率最大、供應(yīng)時(shí)間最短、供應(yīng)成本最低的DTMC-MOP模型，運(yùn)用該模型可以很好地解決因風(fēng)險(xiǎn)擾動(dòng)引起的生產(chǎn)、運(yùn)輸能力變化，得到3個(gè)目標(biāo)的總體最優(yōu)策略。

2)采用改進(jìn)自適應(yīng)的NSGA-Ⅱ算法求解模型，優(yōu)化算法的收斂性和分布性，得到精度更高、更加貼近理想Pareto最優(yōu)前沿的結(jié)果。通過改進(jìn)初始化種群、自適應(yīng)交叉和變異算子、自適應(yīng)精英保留策略，提高算法的全局、局部搜索能力，決策者可以依據(jù)應(yīng)急管理核心目標(biāo)或是不同的偏好選擇相適應(yīng)的應(yīng)急方案。