孫 亮，王 冰，郭 棟，徐 藝

(1. 山東理工大學(xué) 交通與車(chē)輛工程學(xué)院， 山東 淄博 255049；2. 上海大學(xué) 機(jī)電工程與自動(dòng)化學(xué)院， 上海 200072)

開(kāi)放式車(chē)輛路徑問(wèn)題(Open Vehicle Routing Problem，OVRP)指企業(yè)租用車(chē)輛來(lái)完成針對(duì)客戶的配送任務(wù)，在滿足一定約束條件下確定相應(yīng)的車(chē)輛行駛路線以有序服務(wù)客戶，實(shí)現(xiàn)決策目標(biāo)最優(yōu)化。企業(yè)所租用車(chē)輛從企業(yè)出發(fā)，完成配送任務(wù)后，不必返回企業(yè)。

從物流配送的實(shí)際營(yíng)運(yùn)過(guò)程看，影響第三方物流模式營(yíng)運(yùn)效果的因素主要有兩個(gè)：①旅行時(shí)間的不確定性；②客戶期望的服務(wù)時(shí)間段。因此，針對(duì)不確定型OVRP優(yōu)化方法的研究，對(duì)于提升不確定環(huán)境下第三方物流模式營(yíng)運(yùn)效果具有重要的理論和現(xiàn)實(shí)意義。

魯棒優(yōu)化方法[1]采用不確定數(shù)據(jù)的邊界特性描述模型參數(shù)的不確定性，有效避免了隨機(jī)優(yōu)化方法在闡述參數(shù)不確定性上過(guò)度依賴先驗(yàn)知識(shí)及服從概率分布假定的弊端。

目前，求解不確定型車(chē)輛路徑問(wèn)題(Vehicle Routing Problem, VRP)的魯棒優(yōu)化方法主要包含兩類：

1)最壞場(chǎng)景魯棒優(yōu)化方法：該方法利用不確定集的邊界值，將不確定的優(yōu)化模型轉(zhuǎn)化為確定的線性規(guī)劃模型，發(fā)現(xiàn)一個(gè)對(duì)所有觀測(cè)值可行的最優(yōu)解，并確保最壞實(shí)現(xiàn)下的最優(yōu)目標(biāo)函數(shù)值達(dá)到最優(yōu)。Sungur等提出了需求不確定的VRP的最壞場(chǎng)景魯棒優(yōu)化模型，給出了三種需求限制約束的魯棒對(duì)應(yīng)式[2]。Hu等提出了描述需求和旅行時(shí)間不確定VRP特征的魯棒優(yōu)化模型，并分別使用變鄰域搜索方法和分支定界法對(duì)各自提出的模型進(jìn)行求解[3]。Agra等針對(duì)旅行時(shí)間不確定的VRP，提出此類問(wèn)題的魯棒優(yōu)化模型[4]。劉洋等以總服務(wù)成本最小化為優(yōu)化目標(biāo)，建立了路徑規(guī)劃問(wèn)題的魯棒優(yōu)化模型[5]。Cao等提出針對(duì)需求不確定的OVRP的魯棒優(yōu)化模型，使用禁忌算法對(duì)該模型進(jìn)行求解[6]。

2) 弱魯棒優(yōu)化方法：該方法是以最小的約束違背來(lái)保證目標(biāo)函數(shù)值始終不超過(guò)合理范圍，從而改善最壞值與預(yù)期值之間的偏差[7]。Han等提出了基于情景描述旅行時(shí)間不確定性的VRP的弱魯棒優(yōu)化模型[8-9]。Wu等針對(duì)旅行時(shí)間不確定的VRP, 提出了一個(gè)抑制目標(biāo)函數(shù)惡化程度的弱魯棒優(yōu)化模型[10]。

本文以旅行時(shí)間不確定的開(kāi)放式車(chē)輛路徑問(wèn)題(Open Vehicle Routing Problems with Uncertain Travel time,OVRP-UT)為研究對(duì)象，提出一個(gè)描述該問(wèn)題特征的弱魯棒優(yōu)化模型。為了進(jìn)一步提高啟發(fā)式算法獲取最優(yōu)解的概率，在超啟發(fā)式算法的框架下，通過(guò)引入預(yù)測(cè)適應(yīng)度函數(shù)，提出一種自設(shè)計(jì)遺傳算法對(duì)弱魯棒優(yōu)化模型進(jìn)行求解。

1 OVRP-UT的弱魯棒優(yōu)化模型

1.1 基本定義

定義1懲罰成本：如果配送車(chē)輛完成時(shí)間超過(guò)客戶允許的最晚完成時(shí)間，根據(jù)合同向客戶支付的罰金稱為懲罰成本。

定義2違約車(chē)輛：產(chǎn)生懲罰成本的車(chē)輛稱為違約車(chē)輛。

定義3關(guān)于旅行時(shí)間的不確定集

定義4極點(diǎn)：?x∈UT,對(duì)于?δ>0，x+δ?UT，則稱x是UT的極點(diǎn)。

定義5極點(diǎn)集:UT中由極點(diǎn)構(gòu)成的子集，用ext(UT)表示。

定義6最壞實(shí)現(xiàn)：ext(UT)中極點(diǎn)對(duì)應(yīng)的模型參數(shù)的觀測(cè)值稱為最壞實(shí)現(xiàn)。

1.2 問(wèn)題描述

一個(gè)完全的賦權(quán)圖G=(V,E)，這里V={0,1,2,…,n}為點(diǎn)集，0表示車(chē)場(chǎng)，C={1,2,…,n}表示客戶集，客戶的數(shù)量為n；E={(i,j)|i,j∈V,i≠j}為邊集。每一條邊賦有一個(gè)旅行時(shí)間，旅行時(shí)間在UT內(nèi)任意取值。所有車(chē)輛具有相同的最大負(fù)載Q。每一個(gè)客戶i賦有一個(gè)需求di，di≤Q。企業(yè)共租用m輛車(chē)來(lái)完成配送任務(wù)，令M={1,…,m}。Bf是配送任務(wù)最晚完成時(shí)間，如果配送任務(wù)完成時(shí)間超過(guò)Bf，則會(huì)產(chǎn)生懲罰成本。ω1表示單位懲罰成本。

對(duì)于UT內(nèi)的任意實(shí)現(xiàn)，該問(wèn)題的優(yōu)化目標(biāo)是為每輛車(chē)確定合理的運(yùn)輸路線，使違約車(chē)輛數(shù)和總懲罰成本的加權(quán)和達(dá)到最小。調(diào)度方案必須滿足以下約束條件：

1)調(diào)度方案對(duì)UT內(nèi)的任意實(shí)現(xiàn)均保持可行；

2)車(chē)輛路線始于企業(yè)，終止于某個(gè)客戶，服務(wù)完成后，車(chē)輛不再返回企業(yè)；

3)客戶點(diǎn)需求必須且只能由一輛車(chē)來(lái)服務(wù)完成；

4)每一條車(chē)輛路線上客戶點(diǎn)的需求量之和不超過(guò)Q。

1.3 問(wèn)題建模

1.3.1 決策變量

sik:第k輛車(chē)到達(dá)客戶i的時(shí)刻。

1.3.2 弱魯棒優(yōu)化模型

(1)

s.t.

(2)

(3)

(4)

(5)

(6)

(7)

(8)

(9)

(10)

(11)

(12)

sik≥0,?i∈N,?k∈M

(13)

(14)

(15)

sk∈{0,1},?k∈M

(16)

根據(jù)文獻(xiàn)[8]的定義給出另一類弱魯棒優(yōu)化模型LRO1如下：

約束條件與NLRO相同。

LRO1的特點(diǎn)在于通過(guò)控制違約車(chē)輛數(shù)，間接實(shí)現(xiàn)對(duì)總懲罰成本的控制。 不難看出，LRO1是NLRO在ε=0時(shí)的特殊情形。

1.4 NLRO的性質(zhì)分析

性質(zhì)1ZNLRO表示NLRO的最優(yōu)解對(duì)應(yīng)的總懲罰成本，ZLRO1表示LRO1的最優(yōu)解對(duì)應(yīng)的總懲罰成本，ZNLRO≤ZLRO1。

由于

所以

則ZNLRO≤ZLRO1。

□

2 求解方法

自設(shè)計(jì)遺傳算法(Automatic design of Genetic Algorithms， AGA)分為兩部分：①根據(jù)預(yù)期適應(yīng)度函數(shù)值最小的原則，篩選出有利于問(wèn)題求解的遺傳算法的要素；根據(jù)適應(yīng)度函數(shù)值最小的原則，篩選出包含當(dāng)前算法優(yōu)化解的種群。②將兩者組合為一個(gè)新的遺傳算法，繼續(xù)對(duì)問(wèn)題進(jìn)行求解，輸出最終算法優(yōu)化解。

AGA的算法流程如圖1所示。

圖1 AGA流程圖Fig.1 Flow chart of AGA

2.1 遺傳算法的關(guān)鍵要素

1)染色體編碼方式：染色體共由兩部分組成，第一部分采用基于客戶的編碼方式表示具體的調(diào)度方案，第二部分記錄該調(diào)度方案的適應(yīng)度函數(shù)值。染色體編碼結(jié)構(gòu)如圖2所示。

圖2 遺傳算法染色體編碼結(jié)構(gòu)示意Fig.2 Chromosome representation of genetic algorithm

2)染色體解碼的具體過(guò)程如下：隨機(jī)生成新序列{i1,…,in}，依次將該序列中的客戶分配給第1輛車(chē)，直到第1輛車(chē)的剩余負(fù)載無(wú)法再為其他客戶服務(wù)。重復(fù)上述過(guò)程，直到所有客戶都能夠被車(chē)輛服務(wù)。

3)適應(yīng)度函數(shù)的選取：目標(biāo)函數(shù)(1)作為染色體的適應(yīng)度函數(shù)。

2.2 粒子群算法的關(guān)鍵要素

2.2.1 粒子的編碼與解碼方法

粒子向量由位置向量、速度向量、適應(yīng)度函數(shù)值三部分組成。

1)位置向量：構(gòu)造1×8維向量X，第1維表示選擇因子集合S中的操作類型，第2維表示交叉因子集合Cr中的操作類型，第3維表示變異因子集合Mu中的操作類型，第4維表示遺傳算法的種群規(guī)模Mga，第5維表示遺傳算法的迭代次數(shù)Tga，第6維表示選擇操作選中的個(gè)體數(shù)ik，第7維表示交叉概率pc，第8維表示變異概率pm。遺傳算子集合及其構(gòu)成詳見(jiàn)表1。

表1 遺傳算子集合及其構(gòu)成

注：不同算子的操作步驟，詳見(jiàn)文獻(xiàn)[12]。

2)速度向量：速度向量也是一個(gè)1×8維向量，用V表示，微粒的運(yùn)行速度限定在[Vmin,Vmax]。

令

(17)

式中，x0(i)表示該粒子構(gòu)造的遺傳算法的第i次迭代得到的最優(yōu)適應(yīng)度函數(shù)值。 令x0={x0(1),…,x0(q)}，對(duì)其做一次累加生成運(yùn)算，即

令

Y=(x0(2)x0(3) …x0(q))T

則a,b可采用下式計(jì)算：

(18)

粒子的編碼結(jié)構(gòu)如圖3所示。

圖3 粒子編碼結(jié)構(gòu)示意Fig.3 Particle representation

2.2.2 更新規(guī)則

速度和位置更新公式為：

i=1,2,3,4,5,6,7,8

(19)

i=1,2,3

(20)

Xi(t)=ceil(Xi(t-1)+Vi(t)),i=4,5,6

(21)

(22)

(23)

其中，ω0為慣性因子的初始值，ωe為慣性因子的最終值，T表示粒子群算法迭代的最大次數(shù)；r1和r2為(0，1)之間的隨機(jī)數(shù)；a1和a2表示加速因子。

2.3 AGA的算法步驟

AGA的偽代碼如算法1所示。

算法1 AGA的偽代碼

注：①以粒子的位置向量表示的遺傳算法參數(shù)作為輸入，采用標(biāo)準(zhǔn)遺傳算法[12]作為該粒子的解碼算法。

2.4 AGA的性質(zhì)分析

性質(zhì)2令P(Am)表示第m個(gè)粒子所構(gòu)造的遺傳算法能夠在有限時(shí)間內(nèi)發(fā)現(xiàn)模型最優(yōu)解的概率，P(Anew)表示新生成遺傳算法在有限時(shí)間內(nèi)能夠發(fā)現(xiàn)模型最優(yōu)解的概率，

表示AGA發(fā)現(xiàn)最優(yōu)解的概率事件，則AGA發(fā)現(xiàn)NLRO的最優(yōu)解的概率滿足如下性質(zhì)：

P(B)>max{A1,…,AM·T,Anew}

證明：粒子群算法的搜索過(guò)程具有隱含并行性，不同粒子的解碼過(guò)程彼此之間不互相影響，因此，Am和Aj是相互獨(dú)立的。令P(Ai)=ri,i=1,…,M·T+1，AGA發(fā)現(xiàn)最優(yōu)解的概率為：

又

ri=1-(1-ri)

顯然

所以

從而，命題成立。

□

性質(zhì)3AGA最壞情況下的時(shí)間復(fù)雜度與遺傳算法的時(shí)間復(fù)雜度相同。

□

性質(zhì)4AGA最壞情況下的空間復(fù)雜度為O(18M+L)。

證明：就空間復(fù)雜度而言，每個(gè)粒子需要18個(gè)變量空間進(jìn)行存儲(chǔ)，另計(jì)算適應(yīng)度函數(shù)值需要L個(gè)輔助變量空間，從算法1看，AGA存儲(chǔ)粒子和生成適應(yīng)度函數(shù)值是重復(fù)利用的，因此AGA的空間復(fù)雜度為O(18M+L)。

□

3 仿真分析

3.1 仿真環(huán)境與參數(shù)設(shè)置

3.1.1 算例選取與基本參數(shù)設(shè)置

在對(duì)NLRO和LRO1的優(yōu)化解進(jìn)行分析時(shí)，通?？梢越柚墨I(xiàn)[13]的標(biāo)準(zhǔn)測(cè)試數(shù)據(jù)(C1～C14,O1～O8,vrpnc1～vrpnc14,F11,F12)。此外，還可利用文獻(xiàn)[14]的標(biāo)準(zhǔn)測(cè)試數(shù)據(jù) (A-n34-k5，B-n34-k5)分析不確定參數(shù)對(duì)NLRO最優(yōu)解的影響。NLRO和LRO1中確定型參數(shù)的取值與標(biāo)準(zhǔn)測(cè)試數(shù)據(jù)中的對(duì)應(yīng)取值相同。令ω1=1,ε=1,θij=3,ρij=1,Λ=225,算法參數(shù)信息詳見(jiàn)表2。

表2 算法參數(shù)設(shè)置

遺傳算法初始種群的生成方式：首先，隨機(jī)生成一組滿足約束(2)～(8)的調(diào)度方案，并記錄該組可行解的每一輛車(chē)對(duì)應(yīng)的最大配送完成時(shí)間Maxroutetimek。其中，Bf采用如下方式生成：

(24)

然后，根據(jù)Bf的值，隨機(jī)生成滿足約束(2)～(16)的調(diào)度方案，并作為遺傳算法的初始種群。

3.1.2 測(cè)試環(huán)境

用MATLAB R2016a 64位實(shí)現(xiàn)了求解NLRO和LRO1的AGA，在 Intel(R) Core(TM) i5-6200U CPU @ 2.3 GHz，8 GB內(nèi)存的計(jì)算機(jī)上進(jìn)行了仿真實(shí)驗(yàn)。

3.2 算法優(yōu)化解的性能分析

針對(duì)表2中的模型與算法參數(shù)信息，根據(jù)表3和表4中的Bf，使用AGA對(duì)模型NLRO和LRO1分別獨(dú)立求解30次，分別取其中的最好解作為各自的算法優(yōu)化解，該解對(duì)應(yīng)的違約車(chē)輛數(shù)記錄在表3和表4中的veh_n和veh_nl，總懲罰成本分別記錄在fc_n和fc_nl這兩列中。每一次的運(yùn)算時(shí)間上限設(shè)定為3600 s。表3和表4中，cum列表示測(cè)試數(shù)據(jù)的客戶數(shù)，problem列表示測(cè)試數(shù)據(jù)的名稱。

表3 C類和O類問(wèn)題優(yōu)化性能分析

表4 vrpnc問(wèn)題的優(yōu)化性能分析

針對(duì)表3和表4中veh_n和veh_nl兩列數(shù)據(jù)，利用Wilcoxon秩和檢驗(yàn)進(jìn)行分析，在違約車(chē)輛數(shù)的優(yōu)化能力方面，兩者無(wú)顯著差異(p>0.05)。對(duì)于表3和表4中fc_n和fc_nl這兩列數(shù)據(jù)，利用Wilcoxon秩和檢驗(yàn)進(jìn)行分析，分析結(jié)果表明，NLRO產(chǎn)生的懲罰成本顯著小于LRO1(p<0.05)。LRO1單純以違約車(chē)輛數(shù)最小為優(yōu)化目標(biāo)，可能存在多個(gè)違約車(chē)輛數(shù)相同的算法優(yōu)化解，當(dāng)?shù)鷹l件終止時(shí)，輸出的不一定是總懲罰成本最小的算法優(yōu)化解。與LRO1不同，NLRO兼顧違約車(chē)輛數(shù)和總懲罰成本兩個(gè)因素。因此，NLRO算法優(yōu)化解對(duì)應(yīng)的總懲罰成本低于LRO1。

表3和表4數(shù)據(jù)表明：對(duì)中小規(guī)模和大規(guī)模測(cè)試數(shù)據(jù)而言，AGA均具有在有限時(shí)間內(nèi)發(fā)現(xiàn)最優(yōu)解的能力。根據(jù)性質(zhì)2，使用AGA求解表2和表3中的算例對(duì)應(yīng)的未發(fā)現(xiàn)最優(yōu)解的概率的范圍在[0.000 1,0.000 7]之間，因此，此時(shí)的算法優(yōu)化解表現(xiàn)出了性質(zhì)1所描述的最優(yōu)解相類似的性質(zhì)。

3.3 AGA的求解性能分析

以文獻(xiàn)[15]提到的超啟發(fā)式遺傳算法作為AGA的對(duì)比算法，記為S-HPSO。使用兩種算法對(duì)NLRO獨(dú)立求解30次。S-HPSO與AGA采用相同的初始染色體種群。S-HPSO的相關(guān)參數(shù)設(shè)置詳見(jiàn)文獻(xiàn)[15]。S-HPSO與AGA遍歷的可行解總數(shù)相同。其他參數(shù)設(shè)置詳見(jiàn)第3.1小節(jié)。分別將兩種算法對(duì)應(yīng)的30次運(yùn)算結(jié)果的平均值、標(biāo)準(zhǔn)差、最差適應(yīng)度函數(shù)值、最好適應(yīng)度函數(shù)值、運(yùn)算時(shí)間的均值記錄于表5中。表5中標(biāo)有“*”的測(cè)試數(shù)據(jù)表示算法搜索到最優(yōu)目標(biāo)函數(shù)值為0的最優(yōu)解。

在表5中，Avg_a和Avg_s分別表示兩種算法在30次運(yùn)算結(jié)果中所有最優(yōu)適應(yīng)度函數(shù)值的平均值，這兩列中的括號(hào)內(nèi)數(shù)字表示該測(cè)試數(shù)據(jù)搜索到模型最優(yōu)解的次數(shù)。Std_a和Std_s分別表示兩種算法在30次運(yùn)算結(jié)果當(dāng)中所有最優(yōu)適應(yīng)度函數(shù)值的標(biāo)準(zhǔn)差。Tavg_a和Tavg_s分別表示兩種算法對(duì)應(yīng)的平均運(yùn)行時(shí)間。

對(duì)表5中Avg_a和Avg_s這兩列數(shù)據(jù)進(jìn)行Wilcoxon秩和檢驗(yàn)， 結(jié)果表明AGA的平均求解性能顯著優(yōu)于S-HPSO的假設(shè)成立(p<0.05)。對(duì)表5中Tavg_a和Tavg_s這兩列數(shù)據(jù)進(jìn)行Wilcoxon秩和檢驗(yàn)，結(jié)果表明兩種算法的平均運(yùn)行時(shí)間不存在統(tǒng)計(jì)學(xué)差異(p>0.05)。雖然AGA增加了利用預(yù)期適應(yīng)度函數(shù)篩選更新規(guī)則環(huán)節(jié)，但這部分運(yùn)算時(shí)間的增加不足以引起統(tǒng)計(jì)量秩次的變化，因此，兩者在平均運(yùn)算時(shí)間上不存在統(tǒng)計(jì)學(xué)差異。

表5 算法統(tǒng)計(jì)特性的比較與分析

4 結(jié)論

對(duì)OVRP-UT建立弱魯棒優(yōu)化模型并求解，核心內(nèi)容如下： ①模型方面，基于現(xiàn)有的弱魯棒優(yōu)化的框架構(gòu)建了針對(duì)OVRP-UT的弱魯棒優(yōu)化模型；②算法方面，設(shè)計(jì)AGA求解NLRO，提高了啟發(fā)式算法發(fā)現(xiàn)最優(yōu)解的概率，且所得算法優(yōu)化解更接近于最優(yōu)解。

目前求解不確定型車(chē)輛路徑問(wèn)題的魯棒優(yōu)化模型的啟發(fā)式算法主要是元啟發(fā)式算法，下一步研究應(yīng)關(guān)注如何使用超啟發(fā)式算法求解不確定型車(chē)輛路徑問(wèn)題的魯棒優(yōu)化模型。