林明錦，王建新，王 超

(1.重慶大學(xué) 機(jī)械與運(yùn)載工程學(xué)院，重慶 400044；2.太原理工大學(xué) 經(jīng)濟(jì)管理學(xué)院，山西 太原 030002；3.北京工業(yè)大學(xué) 經(jīng)濟(jì)與管理學(xué)院，北京 100124)

0 引言

隨著經(jīng)濟(jì)的全球化、移動(dòng)互聯(lián)網(wǎng)技術(shù)的日趨成熟以及物流基礎(chǔ)設(shè)施的完善，以集群式供應(yīng)鏈為支撐的集成制造模式得以快速發(fā)展。在企業(yè)層面，由于受不同制造中心(以下簡(jiǎn)稱(chēng)客戶(hù))主生產(chǎn)計(jì)劃的制約、市場(chǎng)需求波動(dòng)導(dǎo)致產(chǎn)量變更等不確定因素的影響，各級(jí)供應(yīng)商不僅需要按照主生產(chǎn)計(jì)劃將客戶(hù)所需的零部件按時(shí)配送給客戶(hù)，還要盡可能地動(dòng)態(tài)調(diào)整配送計(jì)劃以滿(mǎn)足客戶(hù)的動(dòng)態(tài)不確定需求，以提升市場(chǎng)響應(yīng)能力，這使得各大供應(yīng)商求解面向客戶(hù)動(dòng)態(tài)需求的帶時(shí)間窗的車(chē)輛路徑問(wèn)題(Vehicle Routing Problem with Time Windows, VRPTW)變得愈加復(fù)雜[1]。傳統(tǒng)制造業(yè)的供應(yīng)商更多關(guān)注如何在滿(mǎn)足時(shí)間要求并遍歷有限客戶(hù)的約束下尋找起止于配送中心的最短車(chē)輛路徑，現(xiàn)代供應(yīng)商為快速響應(yīng)客戶(hù)動(dòng)態(tài)需求，開(kāi)始傾向于在能夠輻射客戶(hù)的區(qū)域建立快速響應(yīng)客戶(hù)的中轉(zhuǎn)站。在政策方面，隨著北京、上海、廣州、重慶等地中重型貨車(chē)24/12小時(shí)限行相關(guān)政策試點(diǎn)的實(shí)施，跨省市長(zhǎng)途運(yùn)輸?shù)拇笮拓涇?chē)直接為客戶(hù)提供送貨服務(wù)受到了嚴(yán)格限制，迫使供應(yīng)商選擇以中轉(zhuǎn)站為分界點(diǎn)，采用不同額定運(yùn)力的車(chē)輛滿(mǎn)足干線(xiàn)與支線(xiàn)零部件供應(yīng)的需求，進(jìn)一步催生了對(duì)考慮動(dòng)態(tài)度和時(shí)間窗的兩級(jí)車(chē)輛路徑問(wèn)題(Two-Echelon Vehicle Routing Problem with Time Window considering Dynamic Degree, 2E-VRPTWDD)的研究。

2E-VRPTWDD是兩級(jí)車(chē)輛路徑問(wèn)題(Two-Echelon Vehicle Routing Problem, 2E-VRP)的進(jìn)一步擴(kuò)展，其中一級(jí)路徑為車(chē)輛從配送中心出發(fā)遍歷所有中轉(zhuǎn)站時(shí)形成的路徑，二級(jí)路徑為車(chē)輛由中轉(zhuǎn)站出發(fā)遍歷所有客戶(hù)時(shí)形成的路徑。在服務(wù)集成制造系統(tǒng)的物流運(yùn)輸中，2E-VRPTWDD指在制造單元時(shí)間和動(dòng)態(tài)需求約束下，將零部件從配送中心到客戶(hù)的配送過(guò)程劃分為兩級(jí)，一級(jí)配送網(wǎng)絡(luò)的節(jié)點(diǎn)要素包括配送中心和中轉(zhuǎn)站，二級(jí)配送網(wǎng)絡(luò)的節(jié)點(diǎn)要素包括中轉(zhuǎn)站和客戶(hù)。供應(yīng)商則通過(guò)兩級(jí)協(xié)調(diào)運(yùn)作，為客戶(hù)提供響應(yīng)市場(chǎng)波動(dòng)需求的零部件供應(yīng)服務(wù)。

針對(duì)2E-VRP，葛顯龍等[2]以?xún)杉?jí)配送中車(chē)輛的總配送距離為目標(biāo)函數(shù)，構(gòu)建了基于場(chǎng)景動(dòng)態(tài)度的兩級(jí)配送網(wǎng)絡(luò)優(yōu)化模型，并采用禁忌搜索算法求解約束優(yōu)化模型，然而該模型并未考慮客戶(hù)對(duì)個(gè)性化時(shí)間窗的現(xiàn)實(shí)需求；YAN等[3]研究了冷鏈物流配送中兩級(jí)配送的中心選址問(wèn)題，但在選址過(guò)程中未考客戶(hù)的動(dòng)態(tài)度；LIU等[4]研究了具有客戶(hù)需求更新的兩級(jí)物流供應(yīng)商批量訂購(gòu)策略，所提需求更新策略對(duì)動(dòng)態(tài)車(chē)輛路徑問(wèn)題的研究具有借鑒意義；MU等[5]構(gòu)建了帶有時(shí)間約束的階段性車(chē)輛路徑優(yōu)化模型。路徑優(yōu)化約束模型的求解方法分為以分支界定法[6]、動(dòng)態(tài)規(guī)劃法[7]為主的精確式求解方法，以及以模擬退火算法[8]、遺傳算法[9]等為代表的啟發(fā)式算法。為了在短時(shí)間內(nèi)獲得路徑優(yōu)化問(wèn)題的近似最優(yōu)解，較多文獻(xiàn)采用啟發(fā)式算法求解車(chē)輛路徑問(wèn)題約束優(yōu)化模型。例如，WEI等[8]采用模擬退火算法，具體求解過(guò)程為隨機(jī)生成初始解，判斷是否滿(mǎn)足約束條件，隨后隨機(jī)生成滿(mǎn)足約束條件的新解，并對(duì)比新解與舊解，根據(jù)模擬退火準(zhǔn)則判斷是否接受新解，直到達(dá)到設(shè)計(jì)好的終止條件，輸出算法的最好解；LI等[9]采用遺傳算法，與模擬退火算法的不同在于新解的生成規(guī)則以及解的接受規(guī)則，遺傳算法的新解主要通過(guò)交叉和變異產(chǎn)生，是否接受解則由適應(yīng)度函數(shù)決定。其他求解車(chē)輛路徑優(yōu)化約束模型的方法還包括離散布谷鳥(niǎo)算法[10]、蟻群算法[11]、粒子群優(yōu)化算法[12]、迭代變鄰域下降算法[13]等。針對(duì)多目標(biāo)優(yōu)化的建模問(wèn)題，KAROONSOONTAWONG等[14]構(gòu)建了用于求解VRPTW的層次多目標(biāo)模型(hierarchical multi-objective formulation)，ZHANG等[15]構(gòu)建了兩級(jí)優(yōu)化模型，VINCENT等[16]構(gòu)建了字典序優(yōu)化模型。這3種多目標(biāo)優(yōu)化模型均在數(shù)學(xué)模型中設(shè)置了多個(gè)目標(biāo)函數(shù)，不同在于求解規(guī)則，層次多目標(biāo)模型是通過(guò)算法求解多目標(biāo)問(wèn)題，即先對(duì)主目標(biāo)進(jìn)行優(yōu)化，再在主目標(biāo)非劣化的情況下對(duì)次目標(biāo)進(jìn)行優(yōu)化[14]；兩級(jí)優(yōu)化模型是將最小化的目標(biāo)函數(shù)作為約束嵌入約束條件[15]；字典序優(yōu)化模型則是通過(guò)在滿(mǎn)足第1個(gè)目標(biāo)的可行解范圍內(nèi)尋找滿(mǎn)足第2個(gè)目標(biāo)的最佳可行解，在滿(mǎn)足第2個(gè)目標(biāo)的可行解范圍內(nèi)尋找滿(mǎn)足第3個(gè)目標(biāo)的最佳可行解，依次類(lèi)推，最終找到滿(mǎn)足多目標(biāo)的可行解[16]。

綜上可知，現(xiàn)有文獻(xiàn)針對(duì)2E-VRP，已從不同角度展開(kāi)了深入研究并取得一定成果，然而鮮有文獻(xiàn)從考慮客戶(hù)動(dòng)態(tài)需求與時(shí)間要求的角度對(duì)2E-VRP進(jìn)行研究。另外，客戶(hù)的動(dòng)態(tài)度越大，對(duì)路徑優(yōu)化的挑戰(zhàn)性越高，然而少有文獻(xiàn)以動(dòng)態(tài)度為指標(biāo)對(duì)2E-VRP進(jìn)行差異優(yōu)化，以盡可能滿(mǎn)足客戶(hù)需求動(dòng)態(tài)化與時(shí)間個(gè)性化的要求。因此，本文在借鑒現(xiàn)有2E-VRP研究成果基礎(chǔ)上，設(shè)計(jì)了基于動(dòng)態(tài)需求獲取準(zhǔn)則及滿(mǎn)足客戶(hù)動(dòng)態(tài)需求的差異策略，并將其融入所構(gòu)建的2E-VRPTWDD優(yōu)化模型，旨在為提升供應(yīng)商對(duì)客戶(hù)動(dòng)態(tài)需求的快速響應(yīng)力提供理論支撐。為提升求解2E-VRPTWDD的效率，對(duì)串行模擬退火算法的初始方案生成規(guī)則和鄰域搜索策略進(jìn)行了改進(jìn)，并借鑒現(xiàn)有并行算法思想設(shè)計(jì)了改進(jìn)的并行模擬退火算法，最后以H公司為案例對(duì)所設(shè)計(jì)的方法進(jìn)行了驗(yàn)證。

1 問(wèn)題描述及數(shù)學(xué)模型

1.1 問(wèn)題描述

所求解的2E-VRPTWDD(如圖1)中，在任意靜態(tài)時(shí)間段內(nèi)均包括2E-VRP網(wǎng)絡(luò)優(yōu)化。一級(jí)路徑中主要優(yōu)化由配送中心和中轉(zhuǎn)站組成的帶時(shí)間約束的有向連通子網(wǎng)絡(luò)，二級(jí)路徑中主要優(yōu)化由中轉(zhuǎn)站和客戶(hù)組成的帶時(shí)間窗約束的有向連通子網(wǎng)絡(luò)。其中，動(dòng)態(tài)度指在車(chē)輛執(zhí)行配送任務(wù)過(guò)程中動(dòng)態(tài)客戶(hù)與所有客戶(hù)的比例。由于車(chē)輛在執(zhí)行任務(wù)過(guò)程中可能出現(xiàn)動(dòng)態(tài)客戶(hù)，而動(dòng)態(tài)度會(huì)隨時(shí)間的增加而增加，車(chē)輛須做出是否響應(yīng)動(dòng)態(tài)客戶(hù)的決策，在較小動(dòng)態(tài)度時(shí)響應(yīng)動(dòng)態(tài)客戶(hù)雖然能夠快速滿(mǎn)足客戶(hù)需求，但是配送成本亦隨之增加；在較大動(dòng)態(tài)度時(shí)響應(yīng)動(dòng)態(tài)客戶(hù)，雖然能夠降低配送成本，但是不能快速滿(mǎn)足客戶(hù)需求。

2E-VRPTWDD可被描述為兩個(gè)相關(guān)的有向連通網(wǎng)絡(luò)圖中邊的選擇問(wèn)題。基本假設(shè)為：由于中重型車(chē)輛在城市中行駛受限，在不同級(jí)配送網(wǎng)絡(luò)中采用不同類(lèi)型車(chē)輛，在同一級(jí)配送網(wǎng)絡(luò)中采用相同類(lèi)型車(chē)輛；如果可供調(diào)度的車(chē)輛到達(dá)中轉(zhuǎn)站/客戶(hù)的時(shí)間早于中轉(zhuǎn)站/客戶(hù)的最早時(shí)間窗，則需在中轉(zhuǎn)站/客戶(hù)處等待，但車(chē)輛的到達(dá)時(shí)間不能晚于中轉(zhuǎn)站/客戶(hù)要求的最晚時(shí)間窗；一級(jí)配送網(wǎng)絡(luò)的送達(dá)時(shí)間窗應(yīng)滿(mǎn)足二級(jí)配送網(wǎng)絡(luò)中客戶(hù)的時(shí)間窗要求；二級(jí)配送網(wǎng)絡(luò)中存在動(dòng)態(tài)客戶(hù)與靜態(tài)客戶(hù)。其中，靜態(tài)客戶(hù)指在[Ti,Ti+1]周期內(nèi)需求量保持不變的客戶(hù)，動(dòng)態(tài)客戶(hù)包括在配送過(guò)程中新增的客戶(hù)和原有客戶(hù)新增需求衍生的新客戶(hù)。因此，求解問(wèn)題為在考慮動(dòng)態(tài)度的情形下，使不同級(jí)車(chē)輛行駛和使用成本最小化。約束條件為：每個(gè)中轉(zhuǎn)站/客戶(hù)均能在一個(gè)配送周期內(nèi)被唯一車(chē)輛服務(wù)一次，每輛車(chē)只服務(wù)起訖于配送中心/中轉(zhuǎn)站的一條路徑，且滿(mǎn)足車(chē)輛裝載能力、中轉(zhuǎn)站/客戶(hù)時(shí)間窗約束。

本文用到的數(shù)學(xué)符號(hào)及含義如下：

G1,G2為一、二級(jí)有向連通網(wǎng)絡(luò)圖,G1=(V1,E1),G2=(V2,E2)；

V1,V2為一、二級(jí)網(wǎng)絡(luò)圖的節(jié)點(diǎn),V1=V1_d∪V1_c,V2=V2_d∪V2_c；

V1_d，V2_d為一級(jí)網(wǎng)絡(luò)中的配送中心和二級(jí)網(wǎng)絡(luò)中的中轉(zhuǎn)站；

V1_c，V2_c為一、二級(jí)網(wǎng)絡(luò)中服務(wù)的節(jié)點(diǎn)，其中一級(jí)網(wǎng)絡(luò)中服務(wù)的節(jié)點(diǎn)類(lèi)型為中轉(zhuǎn)站，二級(jí)網(wǎng)絡(luò)中服務(wù)的節(jié)點(diǎn)類(lèi)型為客戶(hù)；

N1,N2為一、二級(jí)網(wǎng)絡(luò)中所服務(wù)節(jié)點(diǎn)的數(shù)量；

E1，E2為一、二級(jí)網(wǎng)絡(luò)圖的邊，E1={(i,j)|i,j∈V1,i≠j}，E2={(i,j)|i,j∈V2,i≠j}；

c1,c2為一、二級(jí)車(chē)輛的單位行駛距離成本；

K1，K2為一、二級(jí)車(chē)輛的集合，K1={1,2,3,…,K1}，K2={1,2,3,…,K2}；

k為車(chē)輛編碼，k∈K1∪K2；

Q1_k(k∈K1)，Q2_k(k∈K2)為一、二級(jí)車(chē)輛k的額定裝載量；

q1_i(i∈V1_c)，q2_i(i∈V2_c)為一級(jí)網(wǎng)絡(luò)中的中轉(zhuǎn)站及二級(jí)網(wǎng)絡(luò)中客戶(hù)的需求量；

[ET1_i,LT1_i]，[ET2_i,LT2_i]為服務(wù)一級(jí)網(wǎng)絡(luò)中的中轉(zhuǎn)站及二級(jí)網(wǎng)絡(luò)中客戶(hù)的時(shí)間窗；

ΔT為緩沖時(shí)間；

TS1_i，TS2_i為服務(wù)一級(jí)網(wǎng)絡(luò)中的中轉(zhuǎn)站及二級(jí)網(wǎng)絡(luò)中客戶(hù)的時(shí)間；

dij為節(jié)點(diǎn)i,j之間的歐式距離；

tij為車(chē)輛通過(guò)節(jié)點(diǎn)i,j的行駛時(shí)間；

S為車(chē)輛路徑中所服務(wù)的不同類(lèi)型的節(jié)點(diǎn)集合；

|S|為節(jié)點(diǎn)集合中節(jié)點(diǎn)數(shù)量的絕對(duì)值；

p2_i為客戶(hù)i的動(dòng)態(tài)需求增量概率；

σ2_i-min為客戶(hù)i在評(píng)估區(qū)間內(nèi)歷史需求量方差的最小值；

σ2_i-max為客戶(hù)i在評(píng)估區(qū)間內(nèi)歷史需求量方差的最大值；

Δd2_i為滿(mǎn)足客戶(hù)i動(dòng)態(tài)增量需求的補(bǔ)充量；

η為滿(mǎn)足動(dòng)態(tài)增量需求的調(diào)節(jié)參數(shù)；

paccept為啟用動(dòng)態(tài)需求配送概率閾值；

Dyn為動(dòng)態(tài)度；

Ti,Ti+1為周期時(shí)間刻度值；

Arrays.Sort為升序排序函數(shù)；

A為時(shí)間窗跨度相對(duì)距離權(quán)重值；

d1_oi，d2_oi為第一、二級(jí)節(jié)點(diǎn)i到配送中心之間的距離；

y1_ik(i∈V1_c,k∈K1)，y2_ik(i∈V2_c,k∈K2)為一、二級(jí)車(chē)輛開(kāi)始服務(wù)節(jié)點(diǎn)的時(shí)間；

x1_ijk，x2_ijk為一、二級(jí)中的0-1決策變量，車(chē)輛從節(jié)點(diǎn)i直接行駛至節(jié)點(diǎn)j取值為1，否則為0；

1.2 模型構(gòu)建

(1)動(dòng)態(tài)需求獲取準(zhǔn)則

由于原有客戶(hù)的動(dòng)態(tài)需求增量與歷史波動(dòng)需求量之間存在一定的馬爾科夫相關(guān)性，葛顯龍等[2]根據(jù)這一特性，在解決客戶(hù)動(dòng)態(tài)需求配送問(wèn)題時(shí)提出需求配額準(zhǔn)則。該準(zhǔn)則根據(jù)客戶(hù)的歷史需求水平推算當(dāng)前客戶(hù)的缺貨概率，并通過(guò)與既定的概率閾值比較確定是否為動(dòng)態(tài)需求增量補(bǔ)貨。然而，在確定補(bǔ)貨概率時(shí)，該文獻(xiàn)以歷史需求量方差和評(píng)估區(qū)間內(nèi)最小歷史方差的差值，與歷史方差最大波動(dòng)差值的比作為衡量指標(biāo)，并未指出歷史具體時(shí)刻的需求方差。為此，采用歷史均值對(duì)該準(zhǔn)則進(jìn)行改進(jìn)與完善，所設(shè)計(jì)的原有客戶(hù)動(dòng)態(tài)需求增量概率

(1)

p2_i越大表明該客戶(hù)在評(píng)估的歷史區(qū)間內(nèi)需求量波動(dòng)越大，潛在動(dòng)態(tài)需求增量概率亦越大。為避免根據(jù)概率函數(shù)無(wú)差別為所有客戶(hù)提供增量配送服務(wù)時(shí)對(duì)無(wú)增量需求客戶(hù)的干擾，引入客戶(hù)動(dòng)態(tài)增量需求滿(mǎn)足函數(shù)

(2)

基于客戶(hù)動(dòng)態(tài)增量需求滿(mǎn)足函數(shù)，將對(duì)有動(dòng)態(tài)增量的客戶(hù)提供歷史最大需求量與實(shí)際配送量的差額配送量；對(duì)無(wú)動(dòng)態(tài)增量需求的，則通過(guò)概率閾值直接取消(Ti+Ti+1)/2時(shí)刻的動(dòng)態(tài)配送優(yōu)化，以保障初始2E-VRPTWDD優(yōu)化車(chē)輛配送網(wǎng)絡(luò)的穩(wěn)定性。

(2)基于動(dòng)態(tài)度的配送策略

2E-VRPTWDD中的動(dòng)態(tài)度指在某一時(shí)間區(qū)間內(nèi)動(dòng)態(tài)客戶(hù)數(shù)量與所有客戶(hù)數(shù)量的比值，是制定滿(mǎn)足動(dòng)態(tài)客戶(hù)需求策略的重要依據(jù)。動(dòng)態(tài)度

(3)

由于客戶(hù)的動(dòng)態(tài)需求對(duì)配送網(wǎng)絡(luò)的敏捷性響應(yīng)提出了較高要求，為提升配送網(wǎng)絡(luò)優(yōu)化的穩(wěn)定性，基于動(dòng)態(tài)度高低制定兩種滿(mǎn)足動(dòng)態(tài)需求的配送策略。當(dāng)動(dòng)態(tài)度水平較低時(shí)，在[Ti,(Ti+Ti+1)/2]配送周期內(nèi)，采用配送網(wǎng)絡(luò)局部修復(fù)策略，以盡可能降低對(duì)現(xiàn)有配送網(wǎng)絡(luò)的破壞；當(dāng)動(dòng)態(tài)度較高時(shí)，則在(Ti+Ti+1)/2時(shí)刻啟動(dòng)全局更新策略，以滿(mǎn)足新增客戶(hù)和原有客戶(hù)新增的需求。其中，動(dòng)態(tài)度閾值決定是否啟動(dòng)新的車(chē)輛執(zhí)行配送任務(wù)。當(dāng)原有車(chē)輛裝載量能夠滿(mǎn)足新增動(dòng)態(tài)需求時(shí)，定義動(dòng)態(tài)度為較低水平，采用局部修復(fù)策略；當(dāng)原有車(chē)輛的裝載量不能滿(mǎn)足新增動(dòng)態(tài)需求時(shí)，定義動(dòng)態(tài)度為較高水平，采用全局更新策略。

(3)一級(jí)配送網(wǎng)絡(luò)時(shí)間窗約束準(zhǔn)則

一級(jí)配送網(wǎng)絡(luò)中所服務(wù)中轉(zhuǎn)站的時(shí)間與二級(jí)配送網(wǎng)絡(luò)中所服務(wù)客戶(hù)的時(shí)間之間存在相關(guān)性。因此，需基于中轉(zhuǎn)站所服務(wù)客戶(hù)群的最小最早開(kāi)始時(shí)間與最大最晚開(kāi)始時(shí)間確定中轉(zhuǎn)站的時(shí)間窗約束。各中轉(zhuǎn)站的時(shí)間窗約束

(ET1_i,LT1_i)=

(4)

式中：(ET1_i,LT1_i)為一級(jí)配送網(wǎng)絡(luò)中中轉(zhuǎn)站i的時(shí)間窗；ET2_iJ為二級(jí)配送網(wǎng)絡(luò)中中轉(zhuǎn)站i所服務(wù)的客戶(hù)群的最早開(kāi)始時(shí)間窗，LT2_iJ為二級(jí)配送網(wǎng)絡(luò)中中轉(zhuǎn)站i所服務(wù)客戶(hù)群的最晚開(kāi)始時(shí)間窗；ΔT為一級(jí)配送網(wǎng)絡(luò)時(shí)間窗與二級(jí)配送網(wǎng)絡(luò)時(shí)間窗約束存在實(shí)際意義的可連續(xù)的過(guò)渡性緩沖時(shí)間，用于保障配送車(chē)輛有足夠的時(shí)間在中轉(zhuǎn)站中轉(zhuǎn)零部件。

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

層次多目標(biāo)模型由一系列約束條件和多個(gè)目標(biāo)函數(shù)組成，其目標(biāo)是在滿(mǎn)足約束條件的情況下，盡可能使多個(gè)目標(biāo)均實(shí)現(xiàn)最小化或最大化[14]。本文參考文獻(xiàn)[17-18]構(gòu)建的時(shí)間依賴(lài)型車(chē)輛路徑問(wèn)題(Time Dependent Vehicle Routing Problem, TDVRP)模型，構(gòu)建以車(chē)輛使用數(shù)量成本最小化為主要目標(biāo)函數(shù)，以車(chē)輛行駛距離成本最小化為次要目標(biāo)函數(shù)的層次多目標(biāo)模型。在求解過(guò)程中，優(yōu)先以主目標(biāo)進(jìn)行迭代更新求解，在主目標(biāo)函數(shù)值未得到優(yōu)化的情況下，再以次目標(biāo)進(jìn)行迭代更新求解，直至達(dá)到指定的終止條件后輸出最好解。本文構(gòu)建的層次多目標(biāo)模型如下：

1)目標(biāo)函數(shù)

(5)

(6)

2)約束條件

(7)

(8)

?i∈V1_c,?k∈K1；

(9)

?i∈V2_c,?k∈K2；

(10)

?i∈V1_c,?k∈K1；

(11)

?i∈V2_c,?k∈K2；

(12)

x1_ijk(y1_ik+TS1_i+tij)≤y1_jk,

?i∈V1_c,?j∈V1_c,?k∈K1；

(13)

x2_ijk(y2_ik+TS2_i+tij)≤y2_jk,

?i∈V2_c,?j∈V2_c,?k∈K2；

(14)

?S∈V1_c,?S≠?；

(15)

?S∈V2_c,?S≠?；

(16)

(17)

(18)

(19)

(20)

?h∈V1_c,?k∈K1；

(21)

?h∈V2_c,?k∈K2；

(22)

(23)

(24)

(25)

(26)

3)決策變量取值范圍

y1_ik∈R+,?i∈V1_c,?k∈K1；

(27)

y2_ik∈R+,?i∈V2_c,?k∈K2；

(28)

x1_ijk∈{0,1},

?i,j∈V1_c,i≠j,?k∈K1；

(29)

x2_ijk∈{0,1},

?i,j∈V2_c,i≠j,?k∈K2。

(30)

其中：式(5)和式(6)分別為主目標(biāo)函數(shù)和次目標(biāo)函數(shù)，主目標(biāo)函數(shù)為車(chē)輛使用成本最小化，次目標(biāo)函數(shù)為車(chē)輛行駛距離成本最小化；式(7)和式(8)限制車(chē)輛在兩級(jí)路徑中的實(shí)際裝載量不能超過(guò)額定裝載量；式(9)～式(12)約束每級(jí)車(chē)輛為客戶(hù)提供服務(wù)的時(shí)間必須在客戶(hù)要求的時(shí)間窗內(nèi)；式(13)和式(14)為車(chē)輛服務(wù)客戶(hù)的順序約束；式(15)和式(16)約束每級(jí)車(chē)輛行駛過(guò)程中不能形成子回路；式(17)～式(20)約束每級(jí)的每個(gè)客戶(hù)均得到唯一車(chē)輛的一次服務(wù)；式(21)和式(22)約束每級(jí)車(chē)輛服務(wù)完任意客戶(hù)后必須從該客戶(hù)所在位置離開(kāi)；式(23)～式(26)約束每級(jí)所有車(chē)輛從送配送中心出發(fā)后必須返回配送中心；式(27)～式(30)為決策變量的取值范圍。

2 算法設(shè)計(jì)

因?yàn)橐肟蛻?hù)動(dòng)態(tài)需求增加了求解2E-VRPTWDD的復(fù)雜度，所以在串行模擬退火算法的基礎(chǔ)上設(shè)計(jì)并行模擬退火算法以滿(mǎn)足動(dòng)態(tài)車(chē)輛路徑問(wèn)題。當(dāng)動(dòng)態(tài)度較低時(shí)，直接用并行算法滿(mǎn)足修復(fù)性策略；當(dāng)動(dòng)態(tài)度較高時(shí)，在[Ti,(Ti+Ti+1)/2]時(shí)刻啟用全局更新策略對(duì)車(chē)輛路徑進(jìn)行全局更新。

2.1 基于時(shí)間窗精致度的初始方案生成規(guī)則

Solomon提出推進(jìn)插入啟發(fā)式算法求解經(jīng)典VRPTW[19]，并建立了著名的Solomon測(cè)試算例庫(kù)。該算法生成初始解的步驟為：選取插入成本最小的未分配路徑的客戶(hù)到當(dāng)前的可行路徑序列中，并判定插入后的路徑序列是否滿(mǎn)足車(chē)輛裝載能力和時(shí)間窗約束，滿(mǎn)足則繼續(xù)插入未分配路徑的客戶(hù)，直至超過(guò)車(chē)輛額定裝載約束為止，從而生成一條飽和可行路徑序列；依次重復(fù)上述步驟，直至所有客戶(hù)均被插入可行路徑序列為止，生成初始方案。在選擇潛在待插入客戶(hù)時(shí)，有基于距離準(zhǔn)則和基于時(shí)間準(zhǔn)則兩種選擇策略。SOLOMON[20]提出基于距離準(zhǔn)則策略，即優(yōu)先選擇距離物流配送中心較遠(yuǎn)的客戶(hù)作為潛在備選點(diǎn)；CZECH等[21]則提出基于時(shí)間準(zhǔn)則策略，即優(yōu)先選擇最早開(kāi)始時(shí)間最小的客戶(hù)作為潛在備選點(diǎn)。本文結(jié)合文獻(xiàn)[17-18]提出的基于時(shí)間窗精致度策略確定潛在待插入客戶(hù)，潛在客戶(hù)排序函數(shù)

R(c1)=Arrays.Sort

(A×(LT1_i-ET1_i)-d1_oi)。

(31)

R(c2)=Arrays.Sort(A×

(LT2_i-ET2_i)-d2_oi)。

(32)

該式說(shuō)明時(shí)間窗跨度越短、距離配送中心越遠(yuǎn)的客戶(hù)排名越靠前，越會(huì)被優(yōu)先考慮插入當(dāng)前路徑。

2.2 基于路徑內(nèi)外鄰域搜索策略

鄰域搜索指基于某種策略對(duì)原始路徑進(jìn)行變換產(chǎn)生新解的過(guò)程，按照涉及原始路徑條數(shù)可分為單條路徑內(nèi)部鄰域搜索(Or-opt和2-opt)和兩條路徑間的鄰域搜索(2-opt*和Swap/shift)。搜索策略規(guī)則如下：

(1)Or-opt 隨機(jī)選定一條可行路徑上連續(xù)的若干個(gè)客戶(hù)，對(duì)其在可行路徑上的位置進(jìn)行整體調(diào)整，從而產(chǎn)生新的可行路徑[22]。

(2)2-opt 隨機(jī)選取一條可行路徑上的兩個(gè)點(diǎn)，將第1個(gè)點(diǎn)之前的路徑不變生成第1條新的可行路徑；將第1個(gè)與第2個(gè)點(diǎn)之間的路徑倒序后，添加到第1條新的可行路徑中生成第2條新的可行路徑；將第2個(gè)點(diǎn)之后的路徑不變，添加到第2條新的可行路徑中生成第3條新的可行路徑；最終以目標(biāo)函數(shù)值為評(píng)價(jià)準(zhǔn)則，保留4條可行路徑中最優(yōu)的一條路徑作為當(dāng)前的最優(yōu)可行解[23]。

(3)2-opt*將兩條可行路徑上分別被某一節(jié)點(diǎn)切斷的滯后連續(xù)路徑段互換位置，產(chǎn)生兩條新的可行路徑[24]。

(4)Swap/shift 將兩條可行路徑的兩個(gè)節(jié)點(diǎn)互換或單向插入，產(chǎn)生兩條新的可行路徑[25]。

基于路徑內(nèi)外鄰域的4種搜索策略如圖2所示。在設(shè)計(jì)搜索策略過(guò)程中，這4種策略被選擇的概率均設(shè)定為1/4。具體實(shí)現(xiàn)過(guò)程為隨機(jī)產(chǎn)生一個(gè)[0，1]內(nèi)的隨機(jī)數(shù)，若該隨機(jī)數(shù)的取值落在[0，0.25)，則選擇Or-opt進(jìn)行鄰域搜索；若隨機(jī)數(shù)取值落在[0.75，1)，則選擇Swap/shift進(jìn)行鄰域搜索，更新當(dāng)前可行路徑。

2.3 并行模擬退火優(yōu)化設(shè)計(jì)

1983年，KIRKPATRICK等[26]為解決局部最優(yōu)解的問(wèn)題提出模擬退火算法，算法核心思想是在1953年Metropolis提出Metropolis準(zhǔn)則的基礎(chǔ)上融合了退火過(guò)程。物體溫度越高，其內(nèi)部的分子和原子狀態(tài)越不穩(wěn)定，溫度越低，其內(nèi)部的分子和原子狀態(tài)越穩(wěn)定，模擬退火算法正是通過(guò)模擬這一過(guò)程來(lái)尋找原子狀態(tài)相對(duì)穩(wěn)定的局部最優(yōu)解。在退火過(guò)程中，通過(guò)Metropolis概率接受準(zhǔn)則來(lái)迫使算法跳出局部最優(yōu)解，進(jìn)而尋找全局最優(yōu)解，其策略為：在迭代過(guò)程中，如果系統(tǒng)整體能量梯度下降，則將狀態(tài)的轉(zhuǎn)移概率設(shè)定為1；如果迭代過(guò)程中系統(tǒng)整體能量梯度上升，則狀態(tài)轉(zhuǎn)移能否被接受取決于在[0，1]內(nèi)產(chǎn)生的隨機(jī)數(shù)，以及與能量和溫度相關(guān)的動(dòng)態(tài)概率大小的關(guān)系，如果小于該動(dòng)態(tài)概率，則這種狀態(tài)轉(zhuǎn)移被接受，否則被拒絕。這一過(guò)程有效地使算法迭代過(guò)程跳出局部最優(yōu)解，繼續(xù)迭代尋找下一個(gè)新的局部最優(yōu)解。隨著溫度的不斷下降以及上述迭代過(guò)程的不斷重復(fù)，會(huì)產(chǎn)生若干局部最優(yōu)解，最終通過(guò)適應(yīng)度函數(shù)值挑選出所有局部最優(yōu)解中的全局最好解。串行模擬退火算法的偽代碼如下：

Pseudo code of serial simulated annealing algorithm

1.Start

2.Set x=0 and T=m

3.Public static void main(create new(y(x)))

4.i=random.randint(0,len(y(x)-1))

5.j=random.randint(0,len(y(x)-1))

6.y(x)[i],y(x)[j]=y(x)[j],y(x)[i]

7.Return y(x)

8.While(T>T-min)

9.{ Cost savings=f(y(x+1))-f(y(x))

10. if(Cost savings≥0)

11. y(x+1)=y(x)

12. else

13.{ If(random(0,1)

14. y(x+1)=y(x) }

15.x=x+1

16.T=r*T }

17.End

18.Output optimal solution and related parameters

其中：f(x)為系統(tǒng)在x狀態(tài)下的目標(biāo)函數(shù)值；y(x)為系統(tǒng)在x時(shí)所處的狀態(tài)；y(x+1)為系統(tǒng)在x之后所處的狀態(tài)；r為溫度下降速度調(diào)節(jié)參數(shù)；T為系統(tǒng)整體的溫度值；T-min為算法終止的溫度下限值。模擬退火算法的求解質(zhì)量受初始溫度和降溫速率影響較大，高的初始溫度與緩慢的降溫速率有利于尋找到全局最好解，然而需要花費(fèi)大量的計(jì)算時(shí)間。為此，有學(xué)者將并行移動(dòng)、多馬爾科夫鏈融入模擬退火過(guò)程以實(shí)現(xiàn)算法的并行化計(jì)算[5]，其核心思想是將一條馬爾科夫鏈分裂為在一定時(shí)間、空間內(nèi)既能相互獨(dú)立生長(zhǎng)又能彼此交互信息的多條馬爾科夫鏈。馬爾科夫的決策策略可分為同步和異步兩種策略，基于馬爾科夫異步策略的并行模擬退火算法，實(shí)際上是將計(jì)算均勻分布到不同的線(xiàn)程上，各個(gè)線(xiàn)程獨(dú)立運(yùn)行模擬退火算法，當(dāng)各線(xiàn)程運(yùn)行結(jié)束時(shí)，相互對(duì)比并選擇較優(yōu)的局部最優(yōu)解更替當(dāng)前階段的全局最好解?？紤]到求解效率與動(dòng)態(tài)度的關(guān)系，基于馬爾科夫同步策略所設(shè)計(jì)的用于求解2E-VRPTWDD的并行模擬退火算法流程如圖3所示。

在基于馬爾科夫同步策略的并行模擬退火算法中，嵌入向前插入啟發(fā)式(Push Forward Insertion Heuristic，PFIH)算法，在初始階段，由主線(xiàn)程向各個(gè)分線(xiàn)程傳輸一個(gè)基于時(shí)間窗精致度的初始解{x1,x2，…,xn}，隨后各個(gè)分線(xiàn)程開(kāi)始獨(dú)立運(yùn)行一個(gè)以Metropolis接受準(zhǔn)則為主的過(guò)程，當(dāng)經(jīng)過(guò)一個(gè)固定的路徑尋優(yōu)周期后，各個(gè)分線(xiàn)程單向傳遞并比較彼此Metropolis接受準(zhǔn)則下得到的局部最好解{f(x1),f(x2),…,f(xn)}，選擇最小值作為分線(xiàn)程下一次迭代的初始解。各分線(xiàn)程運(yùn)行結(jié)束后，將彼此得到的局部最好解傳遞到主線(xiàn)程，按照接近最優(yōu)解的程度重置各初始解{x1,x2,…,xn}，并將其循環(huán)輸入各個(gè)分線(xiàn)程，以實(shí)現(xiàn)并行運(yùn)算。其中，所設(shè)計(jì)的并行模擬退火算法偽代碼如下：

Pseudo code of parallel simulated annealing algorithm

1.Set T=m and r=r0

2.Generate initial solution{x1,x2,…,xi} based on PFIH

3.Send {x1,x2,…,xi} to thread {p1,p2,…,pi}

4.For 1:pi

5.Public static void main (create new(y(x)))

6.Char a=random. randint (0,1)

7.Switch (a)

8. Case '0.00≤a<0.25'

9. System.out.printIn (initial solution xi+1based on Or-opt strategy)

10. Break;

11.Case '0.25≤a<0.50'

12. System.out.printIn (initial solution xi+1based on 2-opt strategy)

13. Break;…;

14. Return y(xi+1)

15.While(T>T-min)

16.{ Cost savings=f(y(xi+1))-f(y(xi))

17. If (Cost savings≥0)

18. y(xi+1)=y(xi)

19. else

20.{ If (random (0,1)

21. y(xi+1)=y(xi) }

22.xi=xi+1

23.T=r*T }

24.End

25.System.out.printIn (Update initial solution{x1,x2,…,xi})

26.Choose min{f(y(x1),f(y(x2),…,f(y(xi)} as the current optimal solution

27.If go=false and stop, else continue;

28.Send update {x1,x2,…,xi} to thread{p1,p2,…,pi}…

29.End loop

30.Output optimal solution and related parameters

3 算法性能測(cè)試

PERBOLI等[27]研究了2E-VRP，并公布用于對(duì)該問(wèn)題進(jìn)行測(cè)試的set2～set5不同系列的基準(zhǔn)數(shù)據(jù)(https://www.univie.ac.at/prolog/research/TwoEVRP)。

然而，在PERBOLI研究的2E-VRP中，每個(gè)中轉(zhuǎn)站可以跨區(qū)域?yàn)樗锌蛻?hù)提供配送服務(wù)且沒(méi)有時(shí)間窗約束。在實(shí)際運(yùn)營(yíng)中，考慮到配送員對(duì)區(qū)域的熟悉度，以及對(duì)不確定事件的處理能力，企業(yè)主管通常會(huì)安排最佳配送員負(fù)責(zé)固定區(qū)域內(nèi)的配送。因此，不采用PERBOLI給出的2E-VRP數(shù)據(jù)作為測(cè)試基準(zhǔn)數(shù)據(jù)。另外，本文研究的2E-VRPTWDD屬于兩個(gè)帶時(shí)間窗的車(chē)輛路徑聯(lián)合優(yōu)化問(wèn)題，但目前還沒(méi)有相關(guān)的測(cè)試算例。

因此，選取SOLOMON[20]在1987年給出的帶時(shí)間窗的車(chē)輛路徑R，C，RC系列測(cè)試數(shù)據(jù)，作為進(jìn)一步驗(yàn)證所設(shè)計(jì)模型和算法的基準(zhǔn)數(shù)據(jù)，數(shù)據(jù)下載網(wǎng)址為http://web.cba.neu.edu/～msolomon/。R系列包括R101～R112，且每個(gè)測(cè)試集的100個(gè)客戶(hù)與一個(gè)配送中心呈現(xiàn)隨機(jī)特征，分布在100×100的坐標(biāo)系內(nèi)；C系列包括C101～C109，且每個(gè)測(cè)試集的100個(gè)客戶(hù)與一個(gè)配送中心呈現(xiàn)聚類(lèi)特征，隨機(jī)分布在100×100的坐標(biāo)系內(nèi)；RC系列包括RC101～RC108，且每個(gè)測(cè)試集的100個(gè)客戶(hù)與一個(gè)配送中心呈現(xiàn)聚類(lèi)與隨機(jī)混合特征，分布在100×100的坐標(biāo)系內(nèi)。測(cè)試環(huán)境參數(shù)為MacBook Air 13.3 Core i5，1.8 GHz CPU雙核，8 G內(nèi)存，128 G SSD；Windows 10 64 bit；Java JDK-8u251編程環(huán)境。算法基本參數(shù)包括并行線(xiàn)程數(shù)量為6，連續(xù)未能找到相對(duì)局部改進(jìn)解的算法終止次數(shù)為10，溫度下降系數(shù)為0.8，溫度與成本比例系數(shù)為1；選擇R101，R105，C101，C105，RC101，C105作為對(duì)比實(shí)例。計(jì)算結(jié)果與當(dāng)前已知最好解及文獻(xiàn)[28]的計(jì)算結(jié)果對(duì)比如表1所示。

表1 與已知最好解及文獻(xiàn)[28]優(yōu)化結(jié)果的對(duì)比

由表1可知，本文設(shè)計(jì)的并行模擬退火算法在尋優(yōu)能力上，明顯整體優(yōu)于文獻(xiàn)[28]設(shè)計(jì)的分散搜索算法。所求的R105最好解相比已知最好解的優(yōu)化率達(dá)到1.31%，使用的車(chē)輛數(shù)由Solomon給出的已知最好解的14輛降低到13輛。所求的R105的最優(yōu)路徑所需的13輛車(chē)的行駛路徑分別為：1 48 37 20 9 85 18 61 90 1；1 60 93 99 100 88 58 44 97 1；1 34 66 72 10 82 4 69 55 25 81 1；1 64 63 12 65 50 47 49 1；1 22 74 76 42 23 75 59 1；1 73 40 24 68 57 5 56 26 1；1 53 83 8 91 11 51 2 1；1 96 15 45 39 87 92 101 94 1；1 32 89 19 7 95 1；1 29 13 30 80 79 35 36 78 1；1 43 16 3 41 54 27 1；1 28 70 77 31 52 21 67 33 71 1；1 6 84 46 62 17 86 38 98 14 1。其中，采用并行模擬退火算法優(yōu)化過(guò)程中，主目標(biāo)函數(shù)和次目標(biāo)函數(shù)的迭代收斂圖如圖4所示。結(jié)果表明，所設(shè)計(jì)的算法能夠較快地收斂到相對(duì)較優(yōu)解，較好地滿(mǎn)足VRPTW優(yōu)化的需要。

為進(jìn)一步測(cè)試算法的穩(wěn)定性，在相同的參數(shù)配置下，對(duì)R101,R105,C101,C105,RC101,RC105分別測(cè)試10次，結(jié)果如圖5所示。

圖5為運(yùn)行10次的結(jié)果，在10次測(cè)試過(guò)程中，R101的最大值為1 746.37，最小值為1 646.73，平均值為1 689.69，而已知最好解為1 645.79，相對(duì)已知最好解的平均誤差率為2.67%。R105的最大值為1 469.01，最小值為1 359.32，平均值為1 407.76，而已知最好解為1 377.11，相對(duì)已知最好解的平均誤差率為2.23%。6個(gè)測(cè)試集中相對(duì)已知最好解的平均誤差率最大為C105(4.85%)，最小為RC105(0.95%)?？偠灾?，本文算法均能找到相比已知最好解較好的解，由于啟發(fā)式算法的局限性導(dǎo)致每次求解結(jié)果均有一定波動(dòng)，而10次測(cè)試的相對(duì)平均誤差率均控制在5%的可接受范圍內(nèi)，且在10次內(nèi)基本能夠求得與已知最好解接近甚至更優(yōu)的解，說(shuō)明所設(shè)計(jì)的算法具有一定的穩(wěn)定性，能夠用于求解2E-DVRPTWDD。

為分析Or-opt，2-opt，2-opt*，Swap/shift 4種鄰域操作算子對(duì)算法改進(jìn)的效果，在不同組合下，分別對(duì)算法運(yùn)行10次，結(jié)果表明，融合Or-opt，2-opt，2-opt*，Swap/shift 4種鄰域操作算子能有效提升算法性能，主要原因在于不同的鄰域操作策略代表不同的新解生成規(guī)則，通過(guò)生成隨機(jī)數(shù)的方式隨機(jī)選擇4種鄰域策略時(shí)，能夠有效擴(kuò)大可行解的搜索范圍，有利于在搜索空間內(nèi)尋找到相對(duì)更好的解。相反，當(dāng)采用單一鄰域策略時(shí)，由于鄰域搜索操作規(guī)則的一致性，容易導(dǎo)致算法陷入局部最優(yōu)解，不利于在更大的空間范圍內(nèi)搜索到相對(duì)更好的解。為進(jìn)一步探索并行數(shù)量對(duì)優(yōu)化結(jié)果的影響，通過(guò)采用控制變量法的方式設(shè)置不同的實(shí)驗(yàn)組，分別對(duì)其進(jìn)行測(cè)試實(shí)驗(yàn)，具體實(shí)驗(yàn)組的設(shè)置規(guī)則為：在其他參數(shù)保持不變的情況下(τ=10，σ=1.0，γ=1.0，β=0.7，δ=1.8)，將并行線(xiàn)程數(shù)量分別設(shè)置為1，8，10，15，分別對(duì)應(yīng)A，B，C，D 4個(gè)相互對(duì)照的實(shí)驗(yàn)組，每個(gè)實(shí)驗(yàn)組單獨(dú)運(yùn)行10次，統(tǒng)計(jì)每次運(yùn)行所得的車(chē)輛行駛距離成本和車(chē)輛使用數(shù)量成本。

就實(shí)例R101的求解結(jié)果而言，10次測(cè)試均能找到優(yōu)于或與已知最好解接近的解，R101已知最好解中的車(chē)輛行駛距離成本為1 645.79，車(chē)輛使用成本為19；所找到的最好解中的車(chē)輛行駛距離成本為1 624.15，相比已知最好解的優(yōu)化率為1.31%；車(chē)輛使用數(shù)量成本為18，相比已知最好解的優(yōu)化率為5.26%。在4組實(shí)驗(yàn)組中，性能表現(xiàn)最好的是D組，10次測(cè)試的車(chē)輛行駛距離成本平均值為1 679.72，車(chē)輛使用數(shù)量成本為18。表明隨著并行線(xiàn)程數(shù)量的增加，算法的表現(xiàn)性能趨于良好。不同線(xiàn)程數(shù)對(duì)優(yōu)化迭代過(guò)程的影響如圖6所示，圖中橫坐標(biāo)表示算法外層迭代步長(zhǎng)，縱坐標(biāo)表示車(chē)輛的行駛距離成本。

圖6反映了在不同線(xiàn)程參數(shù)設(shè)置下，不同測(cè)試實(shí)驗(yàn)中的迭代收斂過(guò)程。在不同的并行線(xiàn)程設(shè)置下，起始位置均相同，表明初始解的生成只與初始解的構(gòu)造規(guī)則有關(guān)，與并行線(xiàn)程數(shù)量多少無(wú)關(guān)。在求解R101時(shí)，初始解中的車(chē)輛行駛距離成本均為3 458.79，車(chē)輛使用數(shù)量成本均為48。在10次鄰域內(nèi)搜索不到更優(yōu)解便終止算法的條件設(shè)置下，4組實(shí)驗(yàn)的解趨于穩(wěn)定的平均迭代步長(zhǎng)分別為11.6，10.6，9.2，9.8，表明并行數(shù)量對(duì)算法達(dá)到穩(wěn)定狀態(tài)的速度有一定影響，當(dāng)并行數(shù)為10時(shí)，平均迭代步長(zhǎng)最小，為9.2步。然而，隨著并行數(shù)量的增加，算法達(dá)到穩(wěn)定狀態(tài)的耗時(shí)也會(huì)增加，從而增加了尋得相對(duì)最優(yōu)解的時(shí)間成本，這是因?yàn)椴⑿芯€(xiàn)程數(shù)量會(huì)占用更多內(nèi)存，并進(jìn)行更多的內(nèi)部循環(huán)迭代。當(dāng)并行線(xiàn)程設(shè)置為1時(shí)，平均耗時(shí)為90 s左右，而當(dāng)并行線(xiàn)程數(shù)量設(shè)置為15時(shí)，平均耗時(shí)為300 s左右。接下來(lái)，以Perboli給出的2E-VRP中的2eVRP_200-10-3實(shí)例為基準(zhǔn)數(shù)據(jù)，對(duì)比分析連續(xù)兩級(jí)網(wǎng)絡(luò)優(yōu)化和兩個(gè)子網(wǎng)絡(luò)優(yōu)化的結(jié)果，如圖7所示，該實(shí)例包括200個(gè)客戶(hù)和10個(gè)中轉(zhuǎn)站。

其中，圖7a為采用連續(xù)2E-VRP優(yōu)化方法得到的優(yōu)化結(jié)果，車(chē)輛行駛距離成本為2 142.36，在一級(jí)配送網(wǎng)絡(luò)中共需4輛車(chē)執(zhí)行配送任務(wù)，在二級(jí)配送網(wǎng)絡(luò)中共需52輛車(chē)執(zhí)行配送任務(wù)。圖7b為兩個(gè)子網(wǎng)絡(luò)的優(yōu)化結(jié)果，即首先根據(jù)配送員對(duì)區(qū)域的熟悉度對(duì)200個(gè)節(jié)點(diǎn)進(jìn)行聚類(lèi)，以確定中轉(zhuǎn)站的位置，構(gòu)建由配送中心和中轉(zhuǎn)站節(jié)點(diǎn)組成的一級(jí)配送網(wǎng)絡(luò)；其次，構(gòu)建由各中轉(zhuǎn)站與客戶(hù)組成的二級(jí)配送網(wǎng)絡(luò)。優(yōu)化結(jié)果顯示，雖然采用連續(xù)2E-VRP優(yōu)化方法在車(chē)輛成本方面具有優(yōu)勢(shì)，但是由于在優(yōu)化過(guò)程中忽略了配送員對(duì)配送區(qū)域的熟悉度，導(dǎo)致出現(xiàn)跨區(qū)域配送，進(jìn)一步降低了配送效率；當(dāng)采用兩個(gè)子網(wǎng)絡(luò)的優(yōu)化方法時(shí)，雖然車(chē)輛行駛距離成本略高于連續(xù)2E-VRP結(jié)果，但是能夠滿(mǎn)足配送員對(duì)配送區(qū)域熟悉度的要求。因此，將連續(xù)2E-VRP優(yōu)化問(wèn)題轉(zhuǎn)換為兩個(gè)VRP的優(yōu)化問(wèn)題，在提升配送效率方面具有現(xiàn)實(shí)意義。

4 案例研究

H公司作為一家服務(wù)于集成制造產(chǎn)業(yè)群的供應(yīng)商，為上百個(gè)不同的客戶(hù)提供各種零部件的供應(yīng)服務(wù)，各大客戶(hù)由于市場(chǎng)追加訂單、原有零部件供應(yīng)中途故障等不確定因素的影響，會(huì)不定時(shí)地向公司發(fā)出緊急供應(yīng)訂單需求以確保主生產(chǎn)線(xiàn)穩(wěn)定運(yùn)行。

現(xiàn)選取某一周期內(nèi)公司服務(wù)的100個(gè)靜態(tài)客戶(hù)和32個(gè)動(dòng)態(tài)客戶(hù)為案例進(jìn)行研究。鑒于可比較性，使用Solomon給出的C101中的100個(gè)客戶(hù)坐標(biāo)、需求、時(shí)間窗和服務(wù)時(shí)間數(shù)據(jù)表示公司服務(wù)的客戶(hù)，對(duì)比分析公司采用的新方法和原有方法的效果，10個(gè)二級(jí)配送中心和100個(gè)靜態(tài)客戶(hù)的地理分布與需求如表2所示。

表2 靜態(tài)客戶(hù)與所屬中轉(zhuǎn)站及配送中心信息

公司在確定中轉(zhuǎn)時(shí)，考慮到配送員對(duì)區(qū)域的熟悉度，會(huì)采用聚類(lèi)方法確定一個(gè)能夠輻射固定區(qū)域客戶(hù)的中轉(zhuǎn)站。受中、重型車(chē)輛通行的限制，在一級(jí)配送網(wǎng)絡(luò)中使用核定裝載量為10 t的送貨車(chē)，在二級(jí)配送網(wǎng)絡(luò)中使用核定裝載量為2 t的卡車(chē)。在制定滿(mǎn)足動(dòng)態(tài)客戶(hù)的車(chē)輛行駛路徑優(yōu)化方案時(shí)，采用文獻(xiàn)[29]給出的動(dòng)態(tài)客戶(hù)滿(mǎn)足方法，即將完整的一個(gè)物流配送周期劃分為兩部分，在1/2周期時(shí)刻啟動(dòng)滿(mǎn)足客戶(hù)動(dòng)態(tài)需求的配送策略，以滿(mǎn)足前半周期內(nèi)的新增客戶(hù)和新增需求量。執(zhí)行物流配送任務(wù)過(guò)程中的基本參數(shù)如表3所示。

表3 執(zhí)行物流配送任務(wù)過(guò)程中的基本參數(shù)

根據(jù)新增概率和接受動(dòng)態(tài)新增概率閾值規(guī)則確定在[Ti,(Ti+Ti+1)/2]區(qū)間內(nèi)服務(wù)的5個(gè)動(dòng)態(tài)客戶(hù)和未被服務(wù)的靜態(tài)客戶(hù)。在[Ti,(Ti+Ti+1)/2]時(shí)期內(nèi)新增客戶(hù)的數(shù)據(jù)如表4所示。

表4 在[Ti,(Ti+Ti+1)/2]時(shí)期新增的客戶(hù)數(shù)據(jù)

4.1 考慮動(dòng)態(tài)度和時(shí)間窗的兩級(jí)車(chē)輛路徑優(yōu)化

在初始化帶時(shí)間約束的2E-VRP優(yōu)化過(guò)程中，僅考慮對(duì)應(yīng)一個(gè)配送周期[Ti,Ti+1]的靜態(tài)客戶(hù)。采用所設(shè)計(jì)的模型和算法，得到的初始周期的最好2E-VRP路徑如圖8所示。

一級(jí)物流配送網(wǎng)絡(luò)中，共調(diào)用了兩輛額定裝載量為10 t的貨車(chē)為市區(qū)物流配送中心配送貨物，車(chē)輛1的配送路線(xiàn)為0→2→3→4→5→7→0，車(chē)輛2的配送路線(xiàn)為0→6→8→9→10→1→0。兩輛車(chē)的裝載率分別達(dá)到92%，89%，說(shuō)明車(chē)輛的額定裝載量得到了較充分的利用，為車(chē)輛預(yù)留一定空間也符合實(shí)際情況。根據(jù)接受滿(mǎn)足動(dòng)態(tài)新增概率閾值為0.5可知1，2，3，7，8存在動(dòng)態(tài)客戶(hù)。假定新增客戶(hù)動(dòng)態(tài)度為0.5，原有客戶(hù)新增需求動(dòng)態(tài)度為0.2，則根據(jù)式(3)生成[Ti,(Ti+Ti+1)/2]內(nèi)的新增動(dòng)態(tài)客戶(hù)數(shù)量如表5所示。

表5 新增動(dòng)態(tài)客戶(hù)數(shù)量

在滿(mǎn)足新增動(dòng)態(tài)客戶(hù)需求階段，共服務(wù)5個(gè)中轉(zhuǎn)站的24個(gè)新增客戶(hù)和8個(gè)原有客戶(hù)追加需求變成的新客戶(hù)，新增動(dòng)態(tài)客戶(hù)用空心圓表示，在(Ti+Ti+1)/2時(shí)刻優(yōu)化得到的2E-VRPTWDD車(chē)輛路徑如圖9所示。

表6所示為2E-VRPTWDD配送模式下的裝載率、服務(wù)客戶(hù)數(shù)和路徑成本。在設(shè)定2E-VRP網(wǎng)絡(luò)中車(chē)輛額定裝載量不變的情況下，由于動(dòng)態(tài)度的約束關(guān)系，使得(Ti+Ti+1)/2時(shí)刻滿(mǎn)足動(dòng)態(tài)客戶(hù)的車(chē)輛裝載率明顯低于Ti時(shí)刻滿(mǎn)足靜態(tài)客戶(hù)車(chē)輛的裝載率，在一定程度上造成車(chē)載浪費(fèi)。因此，實(shí)際應(yīng)用場(chǎng)景中，在滿(mǎn)足動(dòng)態(tài)客戶(hù)需求時(shí)可考慮采用多車(chē)型配送模式，以減少車(chē)輛裝載的浪費(fèi)。

表6 2E-VRPTWDD配送模式下的裝載率及路徑對(duì)比

4.2 新舊方案對(duì)比及敏感性分析

受道路交通的限制，公司在原有送貨方案中使用額定裝載量為2 t的卡車(chē)為100個(gè)靜態(tài)客戶(hù)和32個(gè)動(dòng)態(tài)客戶(hù)提供零部件供應(yīng)服務(wù)。在一個(gè)配送周期[Ti,Ti+1]的起點(diǎn)，公司以行駛成本和車(chē)輛使用成本最小化為目標(biāo)，調(diào)用10輛車(chē)分別為100個(gè)客戶(hù)提供供貨服務(wù)，其車(chē)輛行駛成本為828.94 km。在(Ti+Ti+1)/2時(shí)刻調(diào)用5輛車(chē)為32個(gè)動(dòng)態(tài)客戶(hù)提供零部件供貨服務(wù)，其車(chē)輛行駛成本為445.47 km。公司原有的兩級(jí)車(chē)輛行駛路徑優(yōu)化如圖10所示。

表7所示為采用原有方案時(shí)的車(chē)輛裝載率和路徑成本。就總成本而言，所設(shè)計(jì)的2E-VRP路徑成本為1 222.79，略?xún)?yōu)于H公司原有的1 274.41。原有配送方案中各路徑上的車(chē)輛裝載率與所設(shè)計(jì)的2E-VRPTWDD模式下各路徑車(chē)輛裝載率保持一致，這是因?yàn)檫@一過(guò)程中均使用裝載量為2 t的車(chē)輛為相同的客戶(hù)群提供服務(wù)，但就各路徑的車(chē)輛行駛成本而言，所設(shè)計(jì)的2E-VRPTWDD模式下的路徑成本均值要優(yōu)于原始方案。

表7 原有2E-VRP網(wǎng)絡(luò)裝載率及路徑對(duì)比

響應(yīng)并滿(mǎn)足客戶(hù)的需求時(shí)間是衡量供應(yīng)商運(yùn)營(yíng)效率的重要指標(biāo)之一，下面從響應(yīng)時(shí)間的角度對(duì)新舊方案進(jìn)行對(duì)比分析。假設(shè)新舊方案中車(chē)輛的行駛速度均相同且均為單位行駛速度，車(chē)輛從配送中心駛出，服務(wù)完對(duì)應(yīng)的客戶(hù)再返回配送中心記為一個(gè)響應(yīng)滿(mǎn)足客戶(hù)時(shí)間，當(dāng)一級(jí)配送網(wǎng)絡(luò)中的10 t大貨車(chē)服務(wù)多個(gè)客戶(hù)群時(shí)，取其平均值作為響應(yīng)滿(mǎn)足各客戶(hù)的時(shí)間。新舊方案的響應(yīng)客戶(hù)時(shí)間如圖11所示。

由圖11可知，新方案的平均響應(yīng)時(shí)間略?xún)?yōu)于舊方案，說(shuō)明所設(shè)計(jì)的方法能夠較好地滿(mǎn)足動(dòng)態(tài)客戶(hù)的快速需求，進(jìn)一步表明所設(shè)計(jì)的方法能夠滿(mǎn)足供應(yīng)商對(duì)客戶(hù)動(dòng)態(tài)需求的響應(yīng)。另外，新方案響應(yīng)時(shí)間的波動(dòng)程度優(yōu)于舊方案，說(shuō)明所設(shè)計(jì)的方法具有較強(qiáng)的穩(wěn)定性。本節(jié)將分析動(dòng)態(tài)度對(duì)車(chē)輛行駛距離成本和車(chē)輛使用數(shù)量成本的影響，將動(dòng)態(tài)度設(shè)置為由0開(kāi)始，每次增加0.05個(gè)單位動(dòng)態(tài)度，直到1結(jié)束的20個(gè)梯度，每個(gè)梯度下運(yùn)行10次取最優(yōu)的距離和車(chē)輛成本作為最終成本，得到的距離和車(chē)輛成本隨動(dòng)態(tài)度變化的規(guī)律如圖12所示。

圖12顯示，距離和車(chē)輛成本總體隨動(dòng)態(tài)度的增加而增加，主要原因是當(dāng)客戶(hù)數(shù)量和需求量增加時(shí)，公司需要安排更多車(chē)輛為增量客戶(hù)提供配送服務(wù)，這一現(xiàn)象與實(shí)際情況相符。車(chē)輛使用成本隨動(dòng)態(tài)度的增加而呈現(xiàn)出階段性增加的特征，例如：動(dòng)態(tài)度在0.55～0.65之間的車(chē)輛使用成本均為6，在0.65～0.70之間的車(chē)輛使用成本均為7，即動(dòng)態(tài)度在一定范圍內(nèi)變化并不會(huì)增加車(chē)輛使用數(shù)量成本，進(jìn)一步表明所設(shè)計(jì)的方法具有良好的魯棒性，其原因是當(dāng)動(dòng)態(tài)客戶(hù)新增需求量在未超過(guò)現(xiàn)有車(chē)輛運(yùn)力時(shí)，現(xiàn)有車(chē)輛有能力為動(dòng)態(tài)客戶(hù)提供配送服務(wù)，而不需要重新啟動(dòng)新的車(chē)輛。因此，在局部動(dòng)態(tài)度變化過(guò)程中，動(dòng)態(tài)度與車(chē)輛數(shù)量使用成本無(wú)關(guān)。

5 結(jié)束語(yǔ)

在各大城市陸續(xù)實(shí)施中、重型貨車(chē)限行政策的背景下，2E-VRP受到重視，而在車(chē)輛執(zhí)行任務(wù)過(guò)程中，客戶(hù)存在動(dòng)態(tài)需求的不確定性，企業(yè)是否及時(shí)響應(yīng)客戶(hù)的動(dòng)態(tài)需求對(duì)于能否提升客戶(hù)滿(mǎn)意度和運(yùn)營(yíng)效率非常重要。當(dāng)企業(yè)接受客戶(hù)的動(dòng)態(tài)需求變更后，原有車(chē)輛路徑方案可能難以滿(mǎn)足新客戶(hù)的訂單需求，因此制定合適的響應(yīng)客戶(hù)動(dòng)態(tài)需求的策略并及時(shí)更新車(chē)輛路徑，對(duì)企業(yè)的重要性不言而喻。為此，本文針對(duì)配送環(huán)節(jié)中客戶(hù)動(dòng)態(tài)需求的不確定性，以及配送時(shí)間個(gè)性化導(dǎo)致的動(dòng)態(tài)配送網(wǎng)絡(luò)響應(yīng)不及時(shí)等現(xiàn)實(shí)問(wèn)題，設(shè)計(jì)基于客戶(hù)動(dòng)態(tài)增量概率閾值和動(dòng)態(tài)度響應(yīng)增量需求的優(yōu)化路徑更新策略，構(gòu)建了以車(chē)輛使用成本和距離成本最小化為主、次優(yōu)化目標(biāo)的2E-VRPTWDD數(shù)學(xué)模型，并設(shè)計(jì)了改進(jìn)的并行模擬退火算法，旨在為考慮動(dòng)態(tài)客戶(hù)需求的2E-VRP優(yōu)化提供方法與理論支撐。通過(guò)本文研究得到的結(jié)論如下：

(1)引入客戶(hù)動(dòng)態(tài)增量需求滿(mǎn)足函數(shù)和動(dòng)態(tài)度衡量指標(biāo)，能夠有效降低動(dòng)態(tài)配送網(wǎng)絡(luò)在空間與時(shí)間維度的求解復(fù)雜度。

(2)相比SOLOMON算例中R105的已知最好解，優(yōu)化率達(dá)到1.31%，說(shuō)明所設(shè)計(jì)的并行模擬退火算法在求解質(zhì)量方面具有一定優(yōu)勢(shì)。

(3)通過(guò)對(duì)比分析新舊方案可知，2E-VRPTWDD模式下的車(chē)輛調(diào)度方案對(duì)提升供應(yīng)商響應(yīng)客戶(hù)需求速度有積極的作用。

(4)通過(guò)分析動(dòng)態(tài)度大小對(duì)優(yōu)化結(jié)果的敏感性可知，車(chē)輛調(diào)度方案在局部動(dòng)態(tài)度變化范圍內(nèi)具有較強(qiáng)的魯棒性，總之，車(chē)輛使用成本與距離成本隨著動(dòng)態(tài)度的增加而增加。

然而，車(chē)輛的裝載率在滿(mǎn)足客戶(hù)動(dòng)態(tài)增量需求配送階段普遍偏低，說(shuō)明在無(wú)差異車(chē)輛配置的情況下存在車(chē)輛裝載浪費(fèi)，而現(xiàn)實(shí)中的物流配送企業(yè)通常配備無(wú)差異的車(chē)輛。因此，下一階段將對(duì)多種車(chē)型靈活組合的動(dòng)態(tài)配送問(wèn)題展開(kāi)研究，以進(jìn)一步探索多車(chē)型的最優(yōu)組合配送策略，為物流企業(yè)減少裝載浪費(fèi)提供科學(xué)的指導(dǎo)。