求解不確定型車輛路徑問題的弱魯棒優(yōu)化方法

孫 亮，王 冰，郭 棟，徐 藝

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

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

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

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

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

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

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

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

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

1.1 基本定義

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

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

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

定義4極點：?x∈UT,對于?δ>0，x+δ?UT，則稱x是UT的極點。

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

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

1.2 問題描述

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

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

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

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

3)客戶點需求必須且只能由一輛車來服務(wù)完成；

4)每一條車輛路線上客戶點的需求量之和不超過Q。

1.3 問題建模

1.3.1 決策變量

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

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的特點在于通過控制違約車輛數(shù)，間接實現(xiàn)對總懲罰成本的控制。 不難看出，LRO1是NLRO在ε=0時的特殊情形。

1.4 NLRO的性質(zhì)分析

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

由于

所以

則ZNLRO≤ZLRO1。

□

2 求解方法

自設(shè)計遺傳算法(Automatic design of Genetic Algorithms， AGA)分為兩部分：①根據(jù)預(yù)期適應(yīng)度函數(shù)值最小的原則，篩選出有利于問題求解的遺傳算法的要素；根據(jù)適應(yīng)度函數(shù)值最小的原則，篩選出包含當(dāng)前算法優(yōu)化解的種群。②將兩者組合為一個新的遺傳算法，繼續(xù)對問題進(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)染色體解碼的具體過程如下：隨機生成新序列{i1,…,in}，依次將該序列中的客戶分配給第1輛車，直到第1輛車的剩余負(fù)載無法再為其他客戶服務(wù)。重復(fù)上述過程，直到所有客戶都能夠被車輛服務(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維表示選擇操作選中的個體數(shù)ik，第7維表示交叉概率pc，第8維表示變異概率pm。遺傳算子集合及其構(gòu)成詳見表1。

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

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

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

令

(17)

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

令

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

則a,b可采用下式計算：

(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)之間的隨機數(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òu)造的遺傳算法能夠在有限時間內(nèi)發(fā)現(xiàn)模型最優(yōu)解的概率，P(Anew)表示新生成遺傳算法在有限時間內(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}

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

又

ri=1-(1-ri)

顯然

所以

從而，命題成立。

□

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

□

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

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

□

3 仿真分析

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

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

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

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

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

(24)

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

3.1.2 測試環(huán)境

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

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

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

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

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

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

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

3.3 AGA的求解性能分析

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

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

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

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

4 結(jié)論

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

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