陳治亞，高輝，徐光明，劉吉華

(1.中南大學(xué) 交通運(yùn)輸工程學(xué)院，湖南 長(zhǎng)沙 410075；2.五邑大學(xué) 軌道交通學(xué)院，廣東 江門(mén) 529020）

帶時(shí)間窗的車輛路徑問(wèn)題(Vehicle Routing Problem with Time Windows,VRPTW)是在車輛路徑問(wèn)題(Vehicle Routing Problem,VRP)基礎(chǔ)上衍生出的組合優(yōu)化問(wèn)題。可描述為：一列車隊(duì)從配送中心裝載產(chǎn)品后，向指定客戶分送貨物，在滿足客戶時(shí)間窗、限載量等約束條件下，規(guī)劃一條最優(yōu)路徑，達(dá)到諸如成本最小、路程最短等目標(biāo)。時(shí)間窗約束的VRP問(wèn)題[1-2]，可分為軟時(shí)間窗VRP(VRP with Soft Time Windows,VRPSTW)、硬時(shí)間窗VRP(VRP with Hard Time Windows,VRPHTW)、模糊時(shí)間窗VRP(VRP with Fuzzy Time Windows,VRPFTW)。例如，CAO等[3]提出基于時(shí)變交通流的多模糊時(shí)間窗車輛路徑問(wèn)題。多目標(biāo)VRP(Multi-objective Vehicle Routing Problem,MOVRP)問(wèn)題[4]根據(jù)現(xiàn)實(shí)情況構(gòu)建多個(gè)目標(biāo)，能更好滿足配送中心或客戶的需求。例如，僅以路徑最短為目標(biāo)，在遇到交通擁堵時(shí)，最短路徑會(huì)成為耗時(shí)最長(zhǎng)的路徑[5]。因此，除了考慮最短路徑，還需考慮運(yùn)輸時(shí)間、客戶滿意度、路徑均衡度、環(huán)境保護(hù)等因素。學(xué)者們考慮較多的是雙目標(biāo)車輛路徑問(wèn)題，如同時(shí)考慮最小配送成本和最大客戶滿意度，運(yùn)輸成本和路徑間的平衡[6-7]。在雙目標(biāo)的基礎(chǔ)上，WANG等[8]研究行駛路程最短、行駛時(shí)間最小、等待時(shí)間最短、車輛數(shù)最少和客戶滿意度最大5個(gè)目標(biāo)問(wèn)題。隨著生活水平的不斷提高，多目標(biāo)VRP會(huì)更加豐富。VRP屬于NP-hard問(wèn)題，目前較多采用智能算法進(jìn)行求解[9]，如禁忌搜索算法和遺傳算法。多目標(biāo)VRP由于目標(biāo)之間的不可比較性，較多學(xué)者采用化多目標(biāo)為單目標(biāo)進(jìn)行解的優(yōu)劣性判斷。蟻群算法是一種模擬進(jìn)化算法，具有許多優(yōu)良的性質(zhì)，在車輛路徑問(wèn)題中得到廣泛的應(yīng)用[10-11]?；谝陨涎芯?，本文以線路里程最短和線路均衡度最好為雙目標(biāo)，考慮隨機(jī)需求的影響性，建立帶硬時(shí)間窗約束的雙目標(biāo)車輛路徑規(guī)劃模型(Multi-objective Vehicle Routing Problem with Hard Time Windows,MOVRPHTW)。采用非支配排序策略選擇迭代后的精英蟻群釋放信息素，獲得貨物配送方案的Pareto最優(yōu)車輛路徑解集，以供配送中心參考。

1 問(wèn)題描述與模型建立

1.1 問(wèn)題描述

客戶需求量受到多種因素的影響，會(huì)出現(xiàn)一定的波動(dòng)性。當(dāng)以期望固定客戶需求，規(guī)劃配送路徑方案時(shí)，若客戶需求量小，則浪費(fèi)車輛運(yùn)力；當(dāng)客戶需求量大，則部分客戶需求無(wú)法得到滿足，需再派車進(jìn)行二次服務(wù)，進(jìn)而增加車輛配送成本。考慮需求的隨機(jī)波動(dòng)性，本文研究的MOVRPHTW問(wèn)題，可以描述為：在一個(gè)產(chǎn)品配送中心，有一列車輛數(shù)為m的車隊(duì)裝載貨物從配送中心出發(fā)，向一組位置、服務(wù)時(shí)間、需求量期望和方差已知的n個(gè)客戶分發(fā)貨物，分發(fā)完貨物后返回配送中心，尋找一個(gè)配送方案，使配送方案線路里程最短和線路均衡度最好。

1.2 模型假設(shè)及符號(hào)說(shuō)明

基本假設(shè)：

1)客戶點(diǎn)位置、服務(wù)時(shí)間、時(shí)間窗已知，由一輛車進(jìn)行配送；

2)客戶需求不確定，滿足相互獨(dú)立的隨機(jī)分布，其中客戶需求期望uj和方差σj已知；

3)配送中心位置，出發(fā)時(shí)間、返回時(shí)間已知；

4)客戶在規(guī)定時(shí)間窗內(nèi)接受服務(wù)，車輛早到等待到ei時(shí)接受服務(wù)，不允許晚到；

5)客戶之間及配送中心和客戶之間的行駛距離固定，距離滿足三角不等式；

6)車輛從配送中心出發(fā)，完成配送任務(wù)后，在規(guī)定時(shí)間前返回配送中心。

假設(shè)運(yùn)輸網(wǎng)絡(luò)[12]G=(V,A)，模型參數(shù)及相關(guān)變量如表1所示。

表1 模型參數(shù)及變量Table 1 Model parameters and variables

1.3 考慮隨機(jī)需求的運(yùn)輸可靠性分析

客戶點(diǎn)的需求是隨機(jī)滿足一定分布，車輛ξk所服務(wù)的運(yùn)輸量如式(2)。

配送中心根據(jù)實(shí)際運(yùn)營(yíng)情況給定一個(gè)可靠性概率θ，車輛k的運(yùn)載能力滿足可靠性約束式(3)，即車輛所服務(wù)客戶的貨物需求總量的概率大于可靠性所要求的概率。

對(duì)于式(2)，客戶點(diǎn)j的貨物需求量服從期望μj，方差σ2j已知的分布。根據(jù)大數(shù)定律[13-14]，車輛運(yùn)輸需求總量服從ξk~N(μk,σk)的正態(tài)分布，其中車輛運(yùn)載期望μk和方差σ2k，如式(4)～(5)所示。

(ξk-μk)/σk服從標(biāo)準(zhǔn)正態(tài)分布，如式(6)所示，其中?k通過(guò)給定的θ可查。當(dāng)車輛限載量滿足式(7)，則滿足式(3)，進(jìn)而可把式(3)轉(zhuǎn)化為式(8)，即車輛的限載量W大于期望和方差的額外部分之和。

1.4 多目標(biāo)車輛路徑優(yōu)化模型

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

約束條件：

式(9)和式(10)是MOVRPHTW的目標(biāo)函數(shù)，目標(biāo)(9)是各線路長(zhǎng)度之和，即線路里程，目標(biāo)(10)是各線路長(zhǎng)度方差，即線路均衡度。式(11)表示1個(gè)客戶只能由1輛車服務(wù)。式(12)～(14)表示車輛在線路上的流量限制。式(16)～(17)是硬時(shí)間窗限制，式(18)是隨機(jī)需求下的運(yùn)輸量限制約束。

2 快速非支配排序的精英蟻群算法

2.1 配送路徑編碼和解碼

根據(jù)MOVRPHTW特點(diǎn)，設(shè)計(jì)以車輛為單位，依次對(duì)需要服務(wù)的客戶進(jìn)行編碼，構(gòu)造行車路徑。假設(shè)n=10，客戶編號(hào)為{1,2,3,4,5,6,7,8,9,10}，需要服務(wù)的客戶產(chǎn)品需求量為{2,3,1,1,2,2,1,1,1,2}，配送中心的車輛編號(hào)為{11,12,13}，1個(gè)可行車輛路徑的編碼及解碼，如圖1所示。

圖1 編碼與解碼Fig.1 Encoding and decoding

路徑編碼序列為：以車輛編號(hào)開(kāi)頭并包含客戶編號(hào)的一組數(shù)據(jù)序列；解碼序列為：車輛11依次服務(wù)客戶1，6，9，運(yùn)載量為 5；車輛12依次服務(wù)客戶3，7，2，運(yùn)載量為 5；車輛13依次服務(wù)客戶8，10，5，4，運(yùn)載量為6。

2.2 螞蟻狀態(tài)轉(zhuǎn)移概率

在螞蟻狀態(tài)概率轉(zhuǎn)移公式中，考慮信息素、路徑距離啟發(fā)函數(shù)、時(shí)間窗寬度及螞蟻等待時(shí)間對(duì)螞蟻轉(zhuǎn)移概率的影響，如式(19)：

其中：τij是信息素啟發(fā)因子，ηij=是距離啟發(fā)因子，widthj=lj-ej是客戶的時(shí)間窗寬度因子，r0=0.5轉(zhuǎn)移規(guī)則參數(shù)，waitj=max{(ej-tjk),0}是客戶的等待時(shí)間，α,β,δ,γ是其參數(shù)重要程度因子。

2.3 基于貪心算法的種群初始化策略

蟻群算法采用隨機(jī)方式生成初始蟻群對(duì)求解結(jié)果影響較小，但相關(guān)研究表明：雖然非隨機(jī)方式生成的初始可行解會(huì)缺乏多樣性，但其有可能產(chǎn)生高質(zhì)量的最優(yōu)近似解。本文對(duì)雙目標(biāo)函數(shù)進(jìn)行歸一化處理，即貪心目標(biāo)：

采用貪心算法和隨機(jī)方式構(gòu)建螞蟻初始蟻群,初始化策略步驟如下所示。

步驟1：車輛隨機(jī)訪問(wèn)第1個(gè)服務(wù)的客戶點(diǎn)。

步驟2：再訪問(wèn)下一點(diǎn)時(shí)，分別把未訪問(wèn)客戶點(diǎn)加入路徑中。如果車輛不滿足運(yùn)載可靠性約束，則車輛返回配送中心，開(kāi)始下一條配送；如果車輛滿足運(yùn)載可靠性約束，若車輛早到達(dá)客戶，則等待后繼續(xù)服務(wù)，若車輛晚到，開(kāi)始下一條配送，即運(yùn)輸路徑必須滿足客戶硬時(shí)間窗約束。車輛能繼續(xù)訪問(wèn)的客戶點(diǎn)即為待訪問(wèn)客戶集合，從待訪問(wèn)客戶集合中，采用如下偽代碼訪問(wèn)客戶點(diǎn)：

注：其pi用于增加初始解的多樣性。

2.4 鄰域搜索操作

VRP隨著客戶點(diǎn)的增多，解空間規(guī)模呈指數(shù)級(jí)增加，單一的蟻群迭代很難搜索到較滿意解，為增加局部搜索能力，本文從配送方案的線路設(shè)計(jì)2種鄰域搜索(Neighborhood Search,NS)操作[15]:

1)路徑內(nèi)破壞與修復(fù)操作：線路內(nèi)隨機(jī)選擇一個(gè)客戶點(diǎn)進(jìn)行破壞與修復(fù)操作，若修復(fù)后的線路不滿足約束條件，則從不滿足的客戶點(diǎn)進(jìn)行拆分，直到滿足約束條件為止，客戶3的破壞與修復(fù)，如圖2所示。

圖2 路徑破壞與修復(fù)Fig.2 Route damage and repair

2)路徑間的交換操作：從不同的2條線路隨機(jī)選2個(gè)客戶點(diǎn)進(jìn)行線路間的交換操作。改變線路的不可行性時(shí)，進(jìn)行線路拆分，直到滿足約束條件為止，交換操作如圖3所示。

圖3 路徑交換Fig.3 Route switching

為避免鄰域操作中，線路拆分帶來(lái)的車輛數(shù)增加，把拆分后最短線路客戶依次從前往后加入到配送方案中其他線路最短客戶中，直到拆分最短線路客戶點(diǎn)全部分配完為止。若存在拆分最短線路客戶點(diǎn)都不滿足約束條件，則此次鄰域搜索操作失敗。這樣既能避免車輛數(shù)增加也能提高車輛運(yùn)載效率。

2.5 非支配排序的精英螞蟻選擇策略

在線路里程和均衡度的解空間{f(x),σ(x)}中，解集之間存在著2種情況[16]：

當(dāng)f1(x)

當(dāng)f1(x)σ2(x)或f1(x)>f2(x),σ1(x)<σ2(x)時(shí)，解1和解2無(wú)法比較優(yōu)劣。

在線路里程和線路均衡度雙目標(biāo)車輛路徑問(wèn)題中，通常存在第2種情況的解集，其對(duì)應(yīng)目標(biāo)空間中的目標(biāo)矢量所構(gòu)成的曲線稱作Pareto最優(yōu)前沿[17-18]。本文設(shè)計(jì)的精英螞蟻信息素釋放策略如下所示：

1)快速非支配排序算子

在Pareto最優(yōu)解集中依據(jù)個(gè)體非劣解水平對(duì)蟻群分層，指引螞蟻向最優(yōu)解集方向搜索[19]。這是一個(gè)循環(huán)適應(yīng)分級(jí)過(guò)程：1)找出群體中非支配解集，記為第1個(gè)非支配層，賦予非支配序irank=1，并從蟻群中去除；2)繼續(xù)尋找余下蟻群中非支配解集，記為第2個(gè)非支配層，賦予非支配序irank= 2；3)繼續(xù)循環(huán)，直到整個(gè)蟻群被分層，同一層內(nèi)非支配序irank相等。

2)螞蟻擁擠距離算子

為了對(duì)同一非支配層內(nèi)的螞蟻個(gè)體進(jìn)行選擇性排序，引入擁擠度概念。計(jì)算步驟如下：

步驟1：螞蟻初始化距離，L[i]d=0。

步驟2：非支配層排序。同一非支配層個(gè)體按目標(biāo)函數(shù)值升序排列。

步驟3：螞蟻個(gè)體邊緣選擇優(yōu)勢(shì)。給定一個(gè)大數(shù)M，令L[1]d=L[end]d=M。

步驟4：根據(jù)式(21)計(jì)算擁擠距離。

步驟5：對(duì)線路里程和均衡度的擁擠距離求和，得到個(gè)體擁擠度進(jìn)行降序排列。

其中L[i+1]g表示第i+1個(gè)螞蟻的第g個(gè)目標(biāo)函數(shù)值。fmaxg,fming分別為蟻群中第g個(gè)目標(biāo)函數(shù)的最大值和最小值。

分別對(duì)線路里程和均衡度重復(fù)步驟2～4，得到螞蟻i的擁擠度L[i]，優(yōu)先選擇非支配層靠前，擁擠度較大的螞蟻，以維持螞蟻種群多樣性。

3)精英螞蟻選擇算子

精英螞蟻選擇策略即保留蟻群中的優(yōu)良個(gè)體進(jìn)行信息素釋放，以防止Pareto最優(yōu)解的丟失。按照非支配排序irank從低到高排序，擁擠度從大到小排序，選擇一定數(shù)量的螞蟻種群進(jìn)行信息素釋放[20]。螞蟻k的信息素釋策略如式(22)～(23)：

其中：ρ是信息素蒸發(fā)系數(shù)；Q是信息素總量。

2.6 算法步驟

設(shè)計(jì)基于非支配排序的精英蟻群算法求解MOVRPHTW問(wèn)題，核心思想有4個(gè)部分：種群初始化、概率轉(zhuǎn)移公式、局部鄰域搜索、信息素更新。在每次蟻群迭代結(jié)束后，依據(jù)2.5節(jié)選擇前25只螞蟻釋放信息素，算法流程如圖4所示。

圖4 非支配排序的精英蟻群算法Fig.4 Elite ant colony algorithm for non-dominated sorting

3 算例仿真實(shí)驗(yàn)

實(shí)驗(yàn)在同一環(huán)境中進(jìn)行，其中CPU主頻為2.30 GHz，內(nèi)存32 G，操作系統(tǒng)為64位Windows7，編程語(yǔ)言為MATLAB R2018，以Solomon中C1類部分客戶進(jìn)行仿真實(shí)驗(yàn)。

3.1 隨機(jī)需求和時(shí)間窗影響分析

為更好體現(xiàn)不同需求和時(shí)間窗帶來(lái)的路徑差異性，以式(9)線路里程為目標(biāo)，選取C101類前35個(gè)客戶，以固定需求為隨機(jī)需求期望，根據(jù)需求量隨機(jī)產(chǎn)生不同范圍方差，此算例中客戶點(diǎn)期望和方差如表2所示，算法的相關(guān)參數(shù)如表3所示，θ對(duì)應(yīng)的?k=1.28。在考慮不同需求和時(shí)間窗情況下建立3種情形，如表4所示，其中軟時(shí)間窗無(wú)論何時(shí)到達(dá)都接受服務(wù)，采用式(24)增加時(shí)間懲罰里程，車輛k從i到j(luò)的懲罰里程為：

表2 客戶點(diǎn)期望和方差Table 2 Customer point expectations and variance

表3 算法參數(shù)水平Table 3 Combination of ACO parameter values

表4 不同情形下的約束條件Table 4 Constraints in different situations

3種情形下的配送路徑和運(yùn)載量如表5所示，配送方案路線圖如圖5所示。鑒于文章篇幅所限，只列出情形2下線路5的服務(wù)時(shí)間及里程懲罰值，如表6所示。情形3下線路5的到達(dá)時(shí)間和等待時(shí)間，如表7所示。

表6 情形2下線路5服務(wù)時(shí)間及懲罰里程Table 6 Line 5 service time and penalty mileage of scenario 2

表7 情形3下線路5到達(dá)時(shí)間及等待時(shí)間Table 7 Line 5 arrival time and waiting time of scenario 3

圖5 3種情形下配送方案Fig.5 Route map of the distribution plan in three scenarios

表5 不同情形下配送路徑Table 5 Distribution route under fixed demand and stochastic demand

在硬時(shí)間窗下，當(dāng)考慮固定需求和隨機(jī)需求時(shí)，即情形1和情形3。如情形1中，4條配送路徑都滿足固定需求時(shí)的車輛運(yùn)載能力，但在隨機(jī)需求下，路徑1，路徑2和路徑3都不滿足運(yùn)載可靠性約束，即配送路徑方案不可行；情形3中，5條路徑都滿足運(yùn)載可靠性約束。情形3和情形1相比，線路里程增加127.49，因此，配送企業(yè)應(yīng)考慮需求的隨機(jī)性對(duì)線路里程的增加程度。

在隨機(jī)需求下，當(dāng)考慮時(shí)間懲罰里程和硬時(shí)間窗時(shí)，即情形2和情形3。如情形2中的線路5，車輛到達(dá)客戶點(diǎn)后即開(kāi)始服務(wù)，但由于服務(wù)客戶點(diǎn)34時(shí)，在接受時(shí)間外到達(dá)產(chǎn)生值為4.4的線路懲罰里程；情形3中，線路5服務(wù)客戶點(diǎn)30時(shí)，車輛早到，那么等待時(shí)間為19.9，在到客戶點(diǎn)30左時(shí)間窗時(shí)開(kāi)始接受服務(wù)。

3.2 參數(shù)靈敏度分析

為分析設(shè)計(jì)算法參數(shù)對(duì)Pareto解的分布產(chǎn)生影響，取C101類前50個(gè)客戶，設(shè)置螞蟻數(shù)量為100，迭代200次，進(jìn)行多次實(shí)驗(yàn)，分析參數(shù)對(duì)Pareto解質(zhì)量的影響。

由圖6可知，總體上隨著權(quán)重因子的增大，解的質(zhì)量變好，但因子權(quán)重過(guò)大，解的質(zhì)量和多樣性變差。

圖6 蟻群參數(shù)α,β,γ,δ靈敏度分析Fig.6 Sensitivity analysis α,β,γ,δ of ant colony parameters

1)參數(shù)α。隨著參數(shù)α的增大，解的質(zhì)量較好,原因是由于精英螞蟻對(duì)解的指引性較強(qiáng)。

2)參數(shù)β。當(dāng)β=3時(shí)，解的質(zhì)量變差，可能是由于里程啟發(fā)因子權(quán)重過(guò)大，導(dǎo)致在路徑選擇中忽略了線路均衡度，具體原因有待進(jìn)一步研究。

3)對(duì)于參數(shù)δ。當(dāng)δ=3時(shí)，解的質(zhì)量、多樣性、分布性變差，可能是由于時(shí)間窗寬度因子權(quán)重過(guò)大，螞蟻優(yōu)先服務(wù)時(shí)間緊迫性高的客戶，這類分析有助于研究客戶時(shí)間緊迫性對(duì)配送路程的影響。

4)對(duì)于參數(shù)γ。當(dāng)γ=3時(shí)，解的質(zhì)量最差，可能是由于等待時(shí)間因子權(quán)重過(guò)大，導(dǎo)致螞蟻過(guò)度等待。

由圖7可知，對(duì)于信息素?fù)]發(fā)因子ρ，當(dāng)ρ較小時(shí)解的質(zhì)量較好，可能是由于此類客戶地理位置分布較為集中，路徑選擇方案較多。

3.3 雙目標(biāo)算例分析

根據(jù)3.1節(jié)和3.2節(jié)，設(shè)置非支配排序蟻群算法相關(guān)參數(shù)，求解最后一次迭代的第1層非支配Pareto解集，如表8所示Pareto解曲線如圖8所示，Pareto解8的配送路徑方案如圖9所示，以供配送中心根據(jù)實(shí)際情況合理選擇配送路徑方案。

表8 C101類客戶Pareto解值Table 8 Class C101 customer Pareto solutions values

圖8 C101類客戶Pareto解Fig.8 Class C101 customer Pareto solutions

4 結(jié)論

1)從配送中心的角度出發(fā)，配送費(fèi)用不僅與行駛里程有關(guān)，還可能與線路均衡度有關(guān)，構(gòu)建雙目標(biāo)車輛路徑問(wèn)題：以線路里程最短和線路均衡度最好為目標(biāo)函數(shù)，考慮車輛運(yùn)輸可靠性約束，建立帶硬時(shí)間窗的車輛路徑模型，并設(shè)計(jì)基于非支配排序的精英蟻群算法進(jìn)行求解。

2)用Solomon中的C101算例，以運(yùn)輸可靠性考慮不同隨機(jī)需求、硬時(shí)間窗約束的影響，可以發(fā)現(xiàn)隨機(jī)需求和硬時(shí)間窗雖然帶來(lái)線路里程的增加，但運(yùn)輸更加可靠。因此配送企業(yè)應(yīng)結(jié)合客戶的配送特點(diǎn)，如服務(wù)時(shí)間的優(yōu)先級(jí)等，合理安排配送路徑。

3)根據(jù)蟻群算法參數(shù)的靈敏度進(jìn)行分析，合理設(shè)置蟻群參數(shù)，求C101類客戶配送方案的Pareto解，可以發(fā)現(xiàn)不同非支配解的線路均衡度變化不大，但線路里程和車輛數(shù)差別較大，配送企業(yè)應(yīng)結(jié)合路徑均衡度對(duì)司機(jī)薪酬和車輛狀況的影響程度以及市場(chǎng)對(duì)產(chǎn)品需求的波動(dòng)情況合理安排配送路徑。