申曉寧，游 璇，陳慶洲，潘紅麗，黃 遙

(1.南京信息工程大學(xué)自動化學(xué)院,江蘇 南京 210044；2.江蘇省大氣環(huán)境與裝備技術(shù)協(xié)同創(chuàng)新中心,江蘇 南京 210044；3.江蘇省大數(shù)據(jù)分析技術(shù)重點實驗室，江蘇 南京 210044)

1 引言

隨著計劃免疫和常態(tài)化防疫工作的深入開展，疫苗的需求量不斷攀升，疫苗冷鏈物流成為現(xiàn)代物流中的一個重要分支[1]。因此，低碳疫苗冷鏈優(yōu)化配送已成為亟需解決的一類實際物流應(yīng)用問題。該問題對固定數(shù)量的配送車輛進行路徑規(guī)劃，它的成功求解，能使疫苗的冷鏈配送服務(wù)更有秩序、效率更高，及時滿足各個地區(qū)的疫苗需求，提高當(dāng)?shù)氐娜嗣窠】蹬c醫(yī)療衛(wèi)生服務(wù)水平，同時有效減少二氧化碳的排放量，改善周圍大氣環(huán)境。

目前，國內(nèi)外研究人員已經(jīng)分別對疫苗的運輸策略和配送方式等進行了研究。de Carvalho等[2]提出了一個多目標(biāo)、多周期的混合整數(shù)線性規(guī)劃模型，用于可持續(xù)疫苗供應(yīng)鏈的設(shè)計和規(guī)劃，并以一家歐洲公司為典型案例進行了詳細研究?；诎ǚ咒N商和零售商的疫苗供應(yīng)鏈，Qi等[3]提出了一個基本模型研究分銷商通過冷鏈或非冷鏈運輸疫苗的條件和冷鏈運輸?shù)牟呗?。鄭浩[4]根據(jù)疫苗的保存條件、流通特征等，對長途和短途的疫苗冷鏈配送問題分別進行了詳細的研究，提出了能夠規(guī)避風(fēng)險且降低成本的運輸策略。陳勇等[5]建立了多品種疫苗的共同配送模型，有利于提升疫苗運輸?shù)臅r效性，并減少運輸成本。鄭春姣[6]以泰豐生物制品有限公司為例，分析現(xiàn)有疫苗冷鏈配送方案存在的不足，并從組織結(jié)構(gòu)、配送方式和冷藏包裝等方面對方案進行了改進，促進了公司在疫苗冷鏈業(yè)務(wù)上的發(fā)展。

煙花算法是由譚營等[7]在2010年首次提出的一種元啟發(fā)式算法。它模擬煙花在夜空中爆炸的過程，將煙花爆炸產(chǎn)生的火花看作問題的一個解，將煙花爆炸產(chǎn)生多個火花的過程看作其搜索鄰域的過程。目前，煙花算法已成功應(yīng)用于復(fù)雜網(wǎng)絡(luò)關(guān)鍵點集確定[8]、多區(qū)域電力系統(tǒng)調(diào)度[9]、電能交易[10]等多個實際問題的求解。本文提出一種雙檔案協(xié)同進化離散多目標(biāo)煙花算法，簡稱DTAMOFWA(Discrete Two-Archive-based Multi-Objective FireWorks Algorithm)。通過仿真實驗驗證了改進策略的有效性，并將DTAMOFWA與4種已有算法進行對比，實驗結(jié)果表明，所提算法在低碳疫苗冷鏈配送問題上具有更好的求解性能。

2 低碳疫苗冷鏈配送問題的數(shù)學(xué)模型

為了使疫苗在規(guī)定的客戶(包括醫(yī)院、各級疾控中心、診所等)時間窗內(nèi)盡快送達并產(chǎn)生較少的碳排放量，本文將低碳疫苗冷鏈配送問題建模為約束多目標(biāo)優(yōu)化模型，最小化運輸成本和平均客戶不滿意度，同時考慮可用車數(shù)量約束、車輛容量約束和模糊時間窗約束[11]。

2.1 問題描述

低碳疫苗冷鏈配送問題的描述如下：存在1個配送中心(編號為0)和n個客戶(編號為1～n)，由一定數(shù)量的同類型配送車輛對各個客戶進行配送服務(wù)，目標(biāo)是在滿足約束的條件下使運輸?shù)目偝杀竞涂蛻舨粷M意度盡可能地小。配送方案需滿足客戶對疫苗需求量和模糊時間窗的要求，每個子回路中所有客戶的需求量之和不超出配送車輛的容量限制，且實際配送車輛數(shù)不超過配送中心可用車輛數(shù)。每輛車從配送中心出發(fā)和回到配送中心時均需進行一定時間的消毒。配送車輛在疫苗配送過程中每間隔一定時間需對制冷設(shè)備添加一次制冷劑。

2.2 符號和變量

低碳疫苗冷鏈配送問題數(shù)學(xué)模型中的符號和變量含義如下所示：

n:客戶數(shù)量；

K:可用車輛數(shù)；

Q:運輸車輛的最大容量(t)；

w:運輸車輛的自身重量(t)；

qi:客戶i的需求量(t)；

dij:從客戶i處到客戶j處的車輛行駛距離(km)；

lij:從客戶i處行駛到客戶j處時的車輛載重(t)；

vij:從客戶i處到客戶j處的車輛行駛速度(km/h)；

Ce:碳稅(元/噸)；

FE:燃油排放參數(shù)(t/L)；

a:車輛行駛加速度(m/s2)；

g:重力加速度常量(m/s2)；

θij:從客戶i處到客戶j處這一路段的路面坡度；

Cr:滾動阻力系數(shù)；

Cd:牽引力系數(shù)；

A:車輛正面表面積(m2)；

ρ:空氣密度(kg/m3)；

P1:單位貨物在行駛單位距離時制冷產(chǎn)生的碳排放量(t/km)；

P2:每小時司機的薪酬(元/小時)；

tdis:車輛在配送中心的消毒時間(h)；

tsi:車輛在客戶i處的服務(wù)時間(h)；

P3:制冷劑單價(元)；

T:添加制冷劑的間隔時間(h)；

[ETi,LTi]:客戶i可接受的時間窗；

[ETi′,LTi′]:客戶i期望的時間窗；

ti:車輛到達客戶i的時刻；

t0:車輛的出發(fā)時刻；

U0:疫苗在出發(fā)時的品質(zhì)；

σ：疫苗品質(zhì)的衰減指數(shù)。

2.3 數(shù)學(xué)模型

低碳疫苗冷鏈配送問題的數(shù)學(xué)模型中，定義決策變量如式(1)和式(2)所示:

(1)

(2)

目標(biāo)函數(shù)和約束條件如式(3)～式(12)所示：

minZ1=C11+C12+C2+C3=

(3)

(4)

(5)

其中，式(3)表示運輸成本目標(biāo)函數(shù)，包括油耗碳排放成本C11、制冷碳排放成本C12、司機薪酬C2和制冷劑成本C3。采用一種綜合的碳排放計量模型，同時考慮行駛距離、速度等多種因素對車輛耗油量的影響，通過燃油排放參數(shù)與耗油量的乘積獲得車輛在該路段行駛產(chǎn)生的碳排放量[12]。式(4)表示平均客戶不滿意度目標(biāo)函數(shù)，包括平均時間不滿意度和平均品質(zhì)損失程度。式(5)為車輛到達客戶i處時的時間滿意度。式(6)～(12)為約束條件，式(6)表示每個客戶只能由一輛車服務(wù)；式(7)和式(8)表示對于每個被服務(wù)的客戶點，一定會有一輛車從某個地點行駛到該客戶點，并從該客戶點離開；式(9)為經(jīng)典的子回路消除約束，保證每輛車的行駛路線中沒有子回路；式(10)為車輛容量約束；式(11)為車輛數(shù)量約束；式(12)為模糊時間窗約束。

3 雙檔案協(xié)同進化的離散多目標(biāo)煙花算法DTAMOFWA

針對低碳疫苗冷鏈配送問題數(shù)學(xué)模型的特點，本文提出雙檔案協(xié)同進化的離散多目標(biāo)煙花算法DTAMOFWA，設(shè)計了消除車輛數(shù)量和容量約束的解碼方式、部分映射爆炸算子和雙外部檔案協(xié)同進化機制，以提高算法的求解性能。

3.1 個體編碼和消除車輛約束的解碼方式

DTAMOFWA采用整數(shù)編碼，對于包含n個客戶(編號為1～n)的問題，每個個體的編碼均為一串整數(shù)序列，其中每個整數(shù)的取值在1～n。由于低碳疫苗冷鏈配送問題帶有較多強約束條件，導(dǎo)致搜索過程中容易產(chǎn)生大量不可行解，不利于最優(yōu)配送方案的求解。因此，在DTAMOFWA中設(shè)計了一種消除車輛數(shù)量和容量約束的解碼方式。首先，對整數(shù)編碼序列按車輛容量約束和模糊時間窗約束進行解碼，即從配送中心出發(fā)，依次訪問編碼中的各個客戶。若某客戶違反上述2個約束中的任一個，則該子回路結(jié)束，并將違反約束的客戶作為下一個新子回路的第1個訪問對象，直到所有的客戶均被訪問，解碼結(jié)束。此后，判斷子回路數(shù)量，即使用的車輛數(shù)是否超過問題可用車輛數(shù)，若不超過，則所得解為可行解；若超過，則對原整數(shù)編碼個體僅按車輛容量約束重新解碼，并采用一種新型修復(fù)算子，使得已經(jīng)滿足車輛容量約束的序列進一步滿足車輛數(shù)量約束，從而得到僅違反模糊時間窗約束的解。

該修復(fù)算子對客戶訪問序列按車輛容量約束解碼后，當(dāng)子回路數(shù)量達到可用車輛數(shù)時，將剩余未訪問的所有客戶放入一個集合中。對該集合依次進行插入操作和置換操作，刪除滿足特定條件的客戶編號。若集合變?yōu)榭占?，則該修復(fù)操作終止；若集合無法變?yōu)榭占?，即?jīng)插入和置換操作后，仍不能同時滿足車輛數(shù)量和容量約束，則隨機初始化一串訪問序列替代原整數(shù)編碼個體，并重新解碼。插入操作是指:依次判斷集合中每個編號是否能在不違反車輛容量約束的前提下插入某個已有子回路。若能，則將該編號插入已有子回路中并將其從集合中刪除；若不能，則保留在集合中等待進行置換操作。置換操作是指：對集合中的每個編號，在不違反車輛容量約束的前提下將其與每個子回路最小客戶需求量對應(yīng)的編號進行置換。若置換成功，則對當(dāng)前集合進行插入操作；若置換失敗，則將該編號繼續(xù)留在集合中。以客戶數(shù)量為8，可用車輛數(shù)為3，車輛最大容量為12為例，2種操作如圖1所示。圖1中，qi表示客戶i的需求量，Y1為執(zhí)行操作前的解，0為配送中心的編號，子回路分別為k1,k2,k3和k4，Y′1為執(zhí)行操作后所得解。

Figure 1 Repair operator for eliminating vehicle quantity constraint and vehicle capacity constraint

本文所提解碼方式能有效消除車輛數(shù)量和容量約束，生成解大部分為可行解，少數(shù)不可行解也可通過修復(fù)算子將其轉(zhuǎn)化成只違反模糊時間窗約束的解。該操作大大增加了尋優(yōu)過程中可行解的獲取數(shù)量，降低了因不可行解數(shù)量大而為搜索帶來的難度。同時，算法僅需對模糊時間窗約束進行處理，提高了算法的搜索效率。

3.2 部分映射爆炸算子

由于低碳疫苗配送問題屬于組合優(yōu)化問題，原來求解連續(xù)優(yōu)化問題的爆炸算子在本文中并不適用。為了使產(chǎn)生的爆炸火花更均勻地分布在以煙花為圓心、一定幅度為半徑的圓形區(qū)域內(nèi)，充分利用有限的計算資源進行搜索，本文設(shè)計了部分映射爆炸算子，考慮通過個體編碼中發(fā)生變化的編號個數(shù)，來表示決策空間中在煙花周圍產(chǎn)生的位移大小。該算子采用新的爆炸半徑計算公式，并融合部分映射交叉操作生成固定數(shù)量的爆炸火花。部分映射交叉所選取的編碼片段長度由所計算的爆炸半徑?jīng)Q定。

假設(shè)煙花種群大小為N，問題決策變量維數(shù)為n，目標(biāo)個數(shù)為m。所提爆炸算子計算第i個煙花Xi的爆炸半徑的方法如下：首先將Xi的m個目標(biāo)值分別進行歸一化并相乘，得到乘積M(Xi)，如式(13)所示;其后，再將M(Xi)歸一化的結(jié)果與決策變量維數(shù)n相乘，以確定Xi的爆炸半徑Ai，即Xi需要發(fā)生爆炸的決策變量個數(shù)即客戶點個數(shù)，如式(14)所示：

(13)

i=1,2,…,N

(14)

其中，fs(Xi)、fsmax和fsmin分別表示煙花Xi的第s個目標(biāo)值、當(dāng)前種群中所有N個煙花在第s個目標(biāo)上的最大值和最小值；M(Xi)為各目標(biāo)歸一化后結(jié)果的乘積；Mmax和Mmin為當(dāng)前種群中所有N個煙花對應(yīng)乘積中的最大值和最小值;[·]表示四舍五入取整操作。根據(jù)式(14)計算煙花的爆炸半徑，當(dāng)某個煙花靠近一個極端解，即某一目標(biāo)值很小時，歸一化后的乘積較小，因此得到的爆炸半徑也較小；當(dāng)2個目標(biāo)值都較大時，歸一化后的乘積較大，因此得到的爆炸半徑也較大；當(dāng)2個目標(biāo)值均比較折衷時，此時煙花在目標(biāo)空間上處于當(dāng)前Pareto前沿的中間位置。若該煙花更靠近坐標(biāo)原點，說明其更接近真實Pareto最優(yōu)前沿(最小化問題)，歸一化后的乘積較小，得到的爆炸半徑較?。蝗粼摕熁x最優(yōu)前沿較遠，則歸一化后的乘積與前一種情況相比較大，得到的爆炸半徑相對較大。根據(jù)上述原理，所提爆炸半徑計算公式能夠合理地控制不同適應(yīng)度煙花的爆炸半徑，使得適應(yīng)度較好的煙花實現(xiàn)較小范圍的精細搜索，適應(yīng)度較差的煙花實現(xiàn)較大范圍的搜索。

令每個煙花產(chǎn)生4種爆炸半徑的爆炸火花，其爆炸半徑如式(15)所示：

(15)

每個爆炸火花的爆炸半徑確定后，對每個煙花，將編碼中的配送中心0刪除，得到一串客戶的訪問序列。隨機選取4個其它煙花的客戶訪問序列，分別與當(dāng)前煙花進行部分映射交叉操作。ri1,ri2,ri3和ri4分別為每次交叉所選取的片段長度，交叉后按各自原配送中心位置解碼可得到2個爆炸火花。以客戶數(shù)量為8的序列為例，設(shè)爆炸半徑即交叉選取的片段長度為4，則部分映射爆炸算子的操作過程如圖2所示。X1和X2為2個煙花個體;子回路為k1,k2和k3;B1和B2為X1和X2的客戶訪問序列。B′1和B′2為將部分片段交叉后得到的客戶訪問序列。此時，B′1未交叉的部分與已交叉部分有重復(fù)的客戶點，B′2也是如此。根據(jù)圖2中的映射關(guān)系，將B′1和B′2未交叉部分的重復(fù)客戶點分別進行替換，得到2個沒有重復(fù)客戶點的客戶訪問序列S1和S2，最后按X1和X2的配送中心位置解碼得到2個爆炸火花E1和E2。

Figure 2 Partial mapping explosion operator

所提部分映射爆炸算子通過控制煙花的爆炸半徑，交換煙花個體中不同長度的片段，使得個體之間能夠進行不同程度的信息交流，讓煙花根據(jù)自身特點在范圍不同的區(qū)域內(nèi)進行搜索，從而實現(xiàn)了算法在局部尋優(yōu)和全局尋優(yōu)之間的平衡。但是，由于該算子只對客戶訪問序列進行操作，配送中心的位置保持不變，限制了種群的多樣性。

3.3 雙外部檔案協(xié)同進化機制

一些已有的多目標(biāo)優(yōu)化算法中處理約束的方法是直接將搜索過程產(chǎn)生的不可行解淘汰，只對可行解進行操作。但是，該方法不僅會丟失不可行解中的有用信息，且還可能導(dǎo)致種群同化現(xiàn)象嚴(yán)重。另外，在約束優(yōu)化問題中，在可行域邊界的不可行解約束違反程度較小，若在其附近進行局部搜索，則易于得到大量可行解。為解決上述問題，本文采用可行解檔案和不可行解檔案協(xié)同進化機制，如圖3所示。DTAMOFWA設(shè)置了可行解檔案和不可行解檔案2個不重疊的外部檔案，分別保存可行解和不可行解中的精英個體。

Figure 3 Schematic diagram of the coevolutionary mechanism of two external archives

對不可行解檔案實施可行性搜索，可能生成新的可行解和不可行解。為了避免對同一個不可行解進行重復(fù)搜索，若不可行解檔案中的個體數(shù)量大于煙花種群規(guī)模N，則每次從不可行解檔案中隨機選擇N個個體進行可行性搜索；否則，對不可行解檔案中的所有個體進行可行性搜索。按照3.1節(jié)消除車輛數(shù)量和容量約束的解碼方式進行解碼，種群中的不可行解均只違反模糊時間窗約束。因此，可行性搜索策略對某個配送方案中違反時間窗約束程度最大的子回路，按客戶可接受到達時間進行排序，所得的新配送方案稱為可行性搜索火花。為了防止產(chǎn)生的新解相似度過高，導(dǎo)致出現(xiàn)早熟收斂的情況，所提可行性搜索策略每次僅對違反約束程度最大的一個子回路進行調(diào)整。以客戶數(shù)量為8，車輛數(shù)為3的序列為例，可行性搜索如圖4所示。其中，k1為違反時間窗約束程度最大的子回路。編號下方的時間為客戶的可接受到達時間，優(yōu)先選擇可接受到達時間早的客戶，以降低不可行解的約束違反程度。

Figure 4 Feasibility search strategy

對于不可行解檔案經(jīng)可行性搜索后產(chǎn)生的可行解，將其與子代中的非支配個體共同進行目標(biāo)驅(qū)動的啟發(fā)式擴展搜索。該搜索操作利用問題目標(biāo)中包含的啟發(fā)信息——碳排放量和客戶不滿意度，將個體中碳排放和客戶不滿意度最大的子回路的客戶編號重新排序。排序時依次選擇碳排放量或客戶不滿意度最小的路段終點編號作為下一個訪問客戶。經(jīng)過該操作后，用所得解更新可行解檔案，使其在可行空間不斷向Pareto最優(yōu)前沿進化。此外，利用可行性搜索產(chǎn)生的不可行解，根據(jù)時間窗約束違反程度更新不可行解檔案，使其在約束空間向可行區(qū)域不斷移動，達到2個外部檔案協(xié)同進化的效果。若不可行解檔案規(guī)模超過預(yù)設(shè)外部檔案的最大規(guī)模，依次選取其中約束違反程度最小的解進入不可行解檔案；否則，不可行解檔案保持不變。

采用雙外部檔案協(xié)同進化機制，能夠充分利用不可行解中的有用信息不斷推進搜索進程，提高搜索效率。此外，由于可行性搜索沒有對子回路中的客戶進行增刪操作，因此得到的新配送方案依然不違反車輛的容量和數(shù)量約束。

3.4 算法的求解步驟

步驟1初始化。對問題模型及算法中的變量和相關(guān)參數(shù)進行初始化設(shè)置。均勻隨機生成規(guī)模為N的煙花種群POP(POPulation)，根據(jù)消除車輛約束的解碼方式進行解碼，并計算其目標(biāo)向量。將POP中的可行解加入可行檔案Feasible，不可行解加入不可行解檔案Infeasible。設(shè)置當(dāng)前目標(biāo)評價次數(shù)Eva=N，最大目標(biāo)評價次數(shù)為Evamax。

步驟2部分映射爆炸算子。每個煙花根據(jù)4種不同的爆炸半徑通過部分映射交叉操作產(chǎn)生8個爆炸火花。更新Eva=Eva+ 8*N。

步驟3變異算子。對每個煙花編碼進行隨機的兩點交換操作。更新Eva=Eva+N。

步驟4將爆炸火花和變異火花中的非支配解加入NDS(Non-Dominated Set)，并將其中的不可行解加入Infeasible。

步驟5對Infeasible進行可行性搜索，得到可行性搜索火花，其中可行解集為F，不可行解集為I，將Infeasible更新為I。更新Eva=Eva+N。

步驟6對NDS和F的并集實施目標(biāo)驅(qū)動的啟發(fā)式擴展搜索，得到擴展搜索火花SS(Search Sparks)。更新Eva=Eva+ 2*|NDS|。

步驟7更新種群。分別根據(jù)Pareto支配概念和ε支配概念[13],用NDS和SS的并集更新POP和Feasible。若規(guī)模超出最大預(yù)設(shè)值，則根據(jù)擁擠距離對POP和Feasible進行裁剪。使用I更新Infeasible。

步驟8選擇策略。從POP和Feasible中各隨機選擇N/2個個體，若Feasible中個體數(shù)量小于N/2，則用POP中的個體補充。在選出的煙花種群中隨機選擇一個個體，計算該個體與其他煙花的相似度，即個體編碼中對應(yīng)位置編號相同的位置個數(shù)。將相似度高于80%的個體進行隨機長度的循環(huán)移位操作。

步驟9終止條件。若Eva≤Evamax，則轉(zhuǎn)移到步驟2；否則循環(huán)結(jié)束，輸出Feasible。

4 仿真實驗與結(jié)果分析

本文所有實驗均在Intel(R)Core(TM)i5-10210U CPU@1.60 GHz，運行內(nèi)存為12 GB的計算機上采用Matlab R2019a軟件運行。根據(jù)對實際疫苗冷鏈配送過程的調(diào)研，生成不同規(guī)模的4個仿真實例，以驗證3種改進策略的有效性和所提算法DTAMOFWA的整體性能。所有對比算法均在仿真實例上分別進行30次獨立實驗，每次實驗的最大目標(biāo)函數(shù)評價次數(shù)Evamax為50 000。實驗中采用反世代距離IGD(Inverted Generational Distance)[14]和超體積HV(HyperVolume)[15]作為綜合性能評價指標(biāo)，對于計算HV時所需的參考點獲取方式如下：將所有算法求得的非支配解在各目標(biāo)上的最差值(最小化問題為最大值，最大化問題為最小值)分別放大1.5倍。將所有算法運行30次所得的全部非支配解合并成一個新的解集，再從該解集中確定出所有的Pareto非支配解作為該問題的近似Pareto最優(yōu)解集。DTAMOFWA的種群規(guī)模N設(shè)置為10，可行解檔案和不可行解檔案的最大規(guī)模設(shè)置為100。

4.1 仿真實例的生成

本文低碳疫苗冷鏈配送問題是一個實際場景中的應(yīng)用問題，由于無法獲取到實際工程數(shù)據(jù)，因此根據(jù)對相關(guān)領(lǐng)域的調(diào)研，在一定參數(shù)范圍內(nèi)隨機生成了4個不同規(guī)模的仿真實例，分別記為C-n15-k6、C-n20-k8、C-n35-k15和C-n55-k23。其中可用車輛數(shù)的估算方式如式(16)所示[16]。

(16)

其中，μ表示配送任務(wù)難度，實驗中設(shè)置μ=0.5；車輛的最大容量Q=15 t。各仿真實例中客戶位置的橫縱坐標(biāo)均在[0,200]內(nèi)隨機產(chǎn)生，配送中心的坐標(biāo)為[100,100]。在[1,5]的整數(shù)范圍內(nèi)隨機生成每個客戶的需求量，并設(shè)置相應(yīng)的時間窗限制和服務(wù)時間。構(gòu)造一個(n+1)維的方陣表示編號為0～n的地點中任意兩點間路段的行駛速度，速度大小均在[30,50]均勻隨機生成，單位為km/h。低碳疫苗冷鏈配送問題模型中的參數(shù)設(shè)置如表1所示[17,18]。

4.2 改進策略的有效性驗證

針對應(yīng)用問題的特點，本文共提出了3種改進策略。為了驗證改進策略的有效性，分別改變DTAMOFWA的相應(yīng)策略，得到3種對比算法。首先，將相鄰配送中心之間的編碼作為一個子回路，用這種解碼方式替換本文所提消除車輛約束的解碼方式，得到對比算法DTAMOFWA-D。其次，將DTAMOFWA的部分映射爆炸算子替換為已有離散煙花算法DMOFWA[19]的爆炸算子，得到對比算法DTAMOFWA-E。最后，在算法搜索過程中直接淘汰所有不可行解，只對可行解進行操作，得到對比算法DTAMOFWA-V。所有算法的種群規(guī)模N均設(shè)為10。實驗結(jié)果如表2所示，表中mean和std分別表示評價指標(biāo)的平均值和標(biāo)準(zhǔn)差，每組數(shù)據(jù)的最佳值加粗表示；“+”“-”“=”為顯著性水平為0.05的Wilcoxon秩和檢驗結(jié)果，分別表示本文算法顯著優(yōu)于、劣于對比算法以及與對比算法沒有顯著差別。

Table 1 Parameter settings in the model of low-carbon cold chain distribution of vaccines

由表2可知，DTAMOFWA在4個仿真實例上的IGD和HV平均值均為最佳；在C-n15-k6和C-n20-k8上的IGD標(biāo)準(zhǔn)差最優(yōu);同時在C-n20-k8上的HV標(biāo)準(zhǔn)差最優(yōu);在剩余實例上IGD和HV標(biāo)準(zhǔn)差與最優(yōu)值也十分接近。Wilcoxon秩和檢驗結(jié)果也顯示DTAMOFWA在IGD和HV上均顯著優(yōu)于其它對比算法。由此表明3種改進策略均是可行且有效的，都能在一定程度上提高DTAMOFWA的收斂性和多樣性。DTAMOFWA的解碼方式能夠消除車輛約束，使得解碼后的解大部分為可行解或僅違反時間窗約束的解，避免了大量不可行解對搜索進程的干擾。部分映射爆炸算子能夠增強對局部區(qū)域的搜索，更好地實現(xiàn)全局搜索和局部搜索之間的平衡，提高了算法的求解精度。DTAMOFWA采用雙外部檔案協(xié)同進化機制，保留尋優(yōu)過程中產(chǎn)生的約束違反程度較小的不可行解。利用可行性搜索產(chǎn)生的可行解輔助可行解檔案進行尋優(yōu)，同時利用可行性搜索產(chǎn)生的不可行解更新不可行解檔案，引導(dǎo)其逐漸向可行域移動，不可行解中的有用信息因此得以利用，提高了算法的求解效率。DTAMOFWA對不可行解檔案進行可行性搜索的策略，消除或降低了個體對模糊時間窗約束的違反程度，有利于加速生成質(zhì)量更高的可行解。

Table 2 Comparison results of verifying the effectiveness of three improved strategies on IGD and HV

4.3 所提算法的性能驗證

本節(jié)實驗使用的對比算法包括3種文獻中已有的求解帶時間窗約束車輛路徑問題的算法HMOEA(Hybrid Multi-Objective Evolutionary Algorithm)[20]、hybrid DE algorithm(hybrid Differential Evolution algorithm)[21]和MOEA-FCS(Multi-Objective Evolutionary Algorithm with Fast Convergence Strategy)[22]。此外，將DTAMOFWA的多目標(biāo)處理框架改為與已有多目標(biāo)煙花算法S-MOFWA[23]相同的框架，編碼方式、個體生成方式及約束處理方法保留不變，由此得到第4種對比算法，記為S-DMOFWA。由于對比算法對可用車輛數(shù)沒有限制，因此統(tǒng)一在對比算法中使用本文所提消除車輛約束的解碼方式。實驗中DTAMOFWA和S-DMOFWA的種群規(guī)模設(shè)置為10，其余算法的種群規(guī)模為100。實驗結(jié)果如表3所示。

由表3可知，針對不同規(guī)模的仿真實例，DTAMOFWA在IGD和HV上均獲得了最佳平均值和標(biāo)準(zhǔn)差。根據(jù)Wilcoxon秩和檢驗結(jié)果可知，DTAMOFWA在所有仿真實例上的IGD和HV的平均值和標(biāo)準(zhǔn)差均顯著優(yōu)于4種對比算法的。可見，DTAMOFWA能夠搜索到收斂精度更高、分布更加寬廣、支配的目標(biāo)空間體積更大的Pareto最優(yōu)前沿，且算法的求解穩(wěn)定性更佳。這是因為：(1)DTAMOFWA通過消除車輛約束的解碼方式，將問題轉(zhuǎn)化為僅帶有模糊時間窗約束的問題，降低了其求解難度;(2)部分映射爆炸算子生成分布較為均勻的爆炸火花，實現(xiàn)了對煙花周圍局部區(qū)域的精細搜索，提高了算法的收斂性;(3)利用可行性搜索產(chǎn)生的不可行解更新不可行解檔案，保留約束違反程度較小的個體，提高了不可行解轉(zhuǎn)化為可行解的概率；(4)對不可行解檔案進行可行性搜索后，使產(chǎn)生的可行解參與可行解檔案中目標(biāo)驅(qū)動的啟發(fā)式擴展搜索，既能增加解的多樣性，又能強化算法對于當(dāng)前Pareto前沿附近區(qū)域的局部搜索，使得算法的求解精度得到提高。因此，本文所提算法在每次求解過程中都能對整個可行域進行充分的探索和利用，同時挖掘不可行解的有用信息，使算法盡可能找到更多的可行非支配解。

Table 3 Comparison results of DTAMOFWA and other algorithms on IGD and HV

為了直觀地比較各算法的收斂速度及精度，以C-n20-k8為例，所有算法的IGD和HV隨目標(biāo)評價次數(shù)變化的收斂曲線如圖5所示。在每種算法運行過程中，每500次目標(biāo)評價次數(shù)對IGD和HV值進行一次采樣，因此在一個收斂曲線中共有100個數(shù)據(jù)點。由圖5可知，在C-n20-k8上，所提算法DTAMOFWA的IGD和HV收斂精度均高于對比算法的，收斂速度也較快。在迭代初期，算法HMOEA的IGD和HV取值均最小，下降速度也超過了DTAMOFWA的，但不久后它的取值基本保持在一個較差水平。這是由于其種群規(guī)模大于DTAMOFWA的，初始解的多樣性比較好，但經(jīng)過一定的迭代次數(shù)后，種群開始陷入局部最優(yōu)，導(dǎo)致所得Pareto前沿曲線在此之后基本不變。由圖5還可觀察到，HMOEA和hybrid DE algorithm的IGD值均在迭代前期出現(xiàn)了退化現(xiàn)象。原因可能是這2種算法缺乏精英保留機制，導(dǎo)致某些適應(yīng)度較好的解被淘汰。所提算法DTAMOFWA由于采用有效的約束處理機制，充分利用可行解和不可行解中的重要進化信息，使得可行解檔案和不可行解檔案不斷更新和進化。同時，變異算子使其具備一定的跳出局部最優(yōu)的能力。因此，DTAMOFWA能夠在前期和中期以較快的速度收斂到精度較高的近似Pareto最優(yōu)前沿，并在后期保持穩(wěn)定，具有更好的收斂性和多樣性。在其它仿真實例上獲取的IGD和HV收斂曲線也能得到類似的結(jié)論。

綜上，本文所提算法是一種能有效處理約束的多目標(biāo)元啟發(fā)式算法，在不同規(guī)模的低碳疫苗冷鏈配送問題中均表現(xiàn)出了較好的收斂性和多樣性，求解準(zhǔn)確且穩(wěn)定，并具有良好的可擴展性。

Figure 5 IGD and HV convergence curves of five algorithms on C-n20-k8

5 結(jié)束語

本文建立了低碳疫苗冷鏈配送問題的約束多目標(biāo)優(yōu)化模型。該模型以最小化企業(yè)運輸成本和平均客戶不滿意度作為優(yōu)化目標(biāo)，同時滿足可用車數(shù)量約束、車輛容量約束和模糊時間窗約束，更貼合實際應(yīng)用，但同時也增加了問題求解的難度。針對該模型強約束條件多的特點，提出了一種雙檔案協(xié)同進化的離散多目標(biāo)煙花算法。該算法設(shè)計了3種改進策略：采用消除車輛約束的解碼方式，部分映射爆炸算子，以及雙外部檔案協(xié)同進化的機制。在規(guī)模遞增的4個仿真實例中，分別驗證了各改進策略的有效性，并將所提算法與4種已有算法進行對比實驗。結(jié)果表明，所提算法在收斂性和多樣性方面均優(yōu)于對比算法，且對問題規(guī)模具有可擴展性，適用于求解低碳疫苗冷鏈配送問題這類約束多目標(biāo)優(yōu)化問題。