混合粒子群算法求解帶軟時間窗的VRPSPD問題

范厚明，劉文琪，徐振林，耿 靜

1.大連海事大學(xué) 交通運(yùn)輸工程學(xué)院，遼寧 大連 116026

2.大連海事大學(xué) 戰(zhàn)略與系統(tǒng)規(guī)劃研究所，遼寧 大連 116026

1 引言

同時集配貨車輛路徑問題（Vehicle Routing Problem with Simultaneous Pickup and Delivery，VRPSPD）是VRP中應(yīng)用廣泛的的變體之一，其擴(kuò)展為每個客戶處擁有集貨和配貨需求且要求車輛服務(wù)一次同時滿足這兩種需求。這種情況在實際生活中也廣泛的存在，如：飲料行業(yè)空瓶的回收、包裝物以及裝載設(shè)備的循環(huán)再利用，工業(yè)有毒物品的回收凈化再處理等[1]，在減少資源浪費的同時，又能夠產(chǎn)生經(jīng)濟(jì)效益?，F(xiàn)有研究表明，這種配送模式與單向配送相比能夠有效地降低成本，提高效益[1-2]。因此隨著環(huán)保意識的提高和逆向物流的發(fā)展，該類問題日益受到研究者的關(guān)注。為了更加貼近現(xiàn)實，考慮客戶對于車輛到達(dá)時間有要求的實際情況，本文求解了更一般化的帶軟時間窗VRPSPD問題，即允許車輛到達(dá)時間違背時間窗約束，并產(chǎn)生相應(yīng)等待或耽擱成本，這樣處理不僅可以更好地從配送企業(yè)角度衡量工作效率，同時也能更好地保證一定的服務(wù)水平。

VRPSPD最早由Min[3]提出，解決了中心圖書館與地方圖書館之間圖書分發(fā)與回庫的問題。Dethloff[4]首次從逆向物流的角度對問題進(jìn)行了定義，并采用一種考慮車輛剩余裝載能力的插入啟發(fā)式算法進(jìn)行求解。由于VRPSPD是NP-hard問題[5]，因此更多的此類問題的求解方法主要集中于能夠在較短時間內(nèi)獲得滿意解的啟發(fā)式算法、模擬退火算法[2]、遺傳算法[6-7]、蟻群算法[8-9]、禁忌搜索[10]、量子進(jìn)化算法[11]，以及混合算法[12]等，獲得了一定的研究成果。然而，由于VRPSPDTW（VRPSPD with Time Window）的決策目標(biāo)、約束條件更為復(fù)雜，使求解難度進(jìn)一步加深，因此相關(guān)研究相對較少，仍有進(jìn)一步研究的需要。曹二保等[5]首次采用差分進(jìn)化算法進(jìn)行求解，并提出一種根據(jù)解滿足約束的情況，進(jìn)行不斷調(diào)節(jié)得出最小所需車輛數(shù)的方法，解決了客戶數(shù)為8和40的帶時間窗的同時集配貨車輛路徑問題（VRPSPDTW）；Wang和Chen等[13]擴(kuò)展了Solomon[14]關(guān)于VRPTW（Vehicle Routing Problem with Time Window）的測試算例，通過CIM（Cheapest Insertion Method）的多種變形產(chǎn)生初始解，并提出一種共同演化遺傳算法對VRPSPDTW進(jìn)行了求解。王超等[15]采用模擬退火算法求解該問題，其中通過RCRS（Residual Capacity and Radial Surcharge）插入準(zhǔn)則獲得初始解，在此基礎(chǔ)上Wang等[16]又進(jìn)一步提出并行模擬退火算法對其求解，均與Wang和Chen等[13]的遺傳算法進(jìn)行比較。黃務(wù)蘭等[17]通過改進(jìn)的全局人工魚群算法與平行模擬退火[16]在VRPSPDTW求解結(jié)果上的比較，驗證了其改進(jìn)算法的有效性。王超等[18]采用離散布谷鳥算法對VRPSPDTW問題進(jìn)行了求解，通過Rank值的非參數(shù)假設(shè)檢驗對算法進(jìn)行了驗證。Belgin等[19]對二級同時取送貨車輛路徑問題進(jìn)行了研究，并采用變鄰域下降搜索結(jié)合局部搜索算法對該問題進(jìn)行了求解。

雖然現(xiàn)有研究對于求解VRPSPDTW已獲得一定的成果，但是對更貼近現(xiàn)實配送業(yè)務(wù)要求的帶軟時間窗VRPSPD研究卻寥寥無幾。鄧愛民等[20]建立了帶軟時間窗VPRSPD的數(shù)學(xué)模型，采用增加了記憶功能的改進(jìn)模擬退火算法對其進(jìn)行了求解，并通過算例求解驗證了算法的可行性，但在求解結(jié)果上仍有改進(jìn)空間。通過上述軟硬時間窗模型在整體成本的比較來看，即使車輛在某客戶處違背時間窗約束，存在總配送成本相對較低的情況，即求解軟時間窗約束的問題有可能產(chǎn)生在硬時間窗約束下難以獲得的可行解，并且該類問題的通用性很強(qiáng)，可以經(jīng)過一些參數(shù)或條件的改變轉(zhuǎn)化為其他VRP擴(kuò)展問題，因此與硬時間窗約束相比，軟時間窗約束具有不可替代的優(yōu)勢。從查閱的文獻(xiàn)中可知，目前將粒子群算法應(yīng)用求解于帶軟時間窗的VRPSPD問題的研究還較少，并且結(jié)合不同的啟發(fā)式算法進(jìn)行優(yōu)勢互補(bǔ)，已然成為學(xué)者研究的重要方向之一。周蓉等[21]提出混合的離散粒子算法對VRPSPDTW進(jìn)行了求解，結(jié)合了變鄰域下降（Variable Neigborhood Descent，VND）和模擬退火算法的選擇機(jī)制，其中是對每個粒子進(jìn)行一定概率地選擇VND。然而在后續(xù)的實驗結(jié)果中得出，不同結(jié)合方式在一定程度上影響算法的性能。為了更好地探索和提供不同算法結(jié)合的思路，本文提出結(jié)合以VND為適應(yīng)性擾動機(jī)制主體的混合粒子群優(yōu)化算法（Particle Swarm Optimization，PSO），對帶軟時間窗的VRPSPD進(jìn)行求解。

PSO是一種基于群體智能搜索算法，模擬鳥類飛行的仿生算法，有著概念簡單清晰，參數(shù)較少，魯棒性好等優(yōu)點，且在各類多維連續(xù)空間優(yōu)化問題上表現(xiàn)出良好的特性[22]，在VRP及相關(guān)領(lǐng)域也有廣泛應(yīng)用，并取得較好的結(jié)果[22-27]。在PSO的參數(shù)中，慣性權(quán)重具有平衡全局和局部搜索的作用，因此本文所提算法在PSO中加入了適應(yīng)性權(quán)重。但是由于其進(jìn)化只是根據(jù)粒子個體自身經(jīng)驗信息和種群中的信息共享，所以容易陷入局部最優(yōu)，出現(xiàn)早熟的現(xiàn)象，且爬山能力弱，不易跳出。本文提出結(jié)合以VND為主體的適應(yīng)性擾動機(jī)制（Adaptive Perturbation Mechanism，APM）的混合粒子群搜索算法，能夠更好地根據(jù)種群陷入局部極值的情況進(jìn)行適應(yīng)性擾動，其與以往采用擾動機(jī)制的文獻(xiàn)[28-30]的不同之處：

（1）算法在種群迭代過程中，很有可能不只一次陷入局優(yōu)，且每次陷入局優(yōu)的情況不同，本文提出線性變化的調(diào)取APM準(zhǔn)則，這樣可以在充分發(fā)揮PSO全局搜索性能的前提下，提高搜索效率。

（2）在Shaking階段中，不是所有的鄰域結(jié)構(gòu)都被采用，也不是隨機(jī)或依已有概率來進(jìn)行鄰域結(jié)構(gòu)的選取[5，16，28，30]，而是根據(jù)適應(yīng)權(quán)重采取輪盤賭的方式進(jìn)行選取。

（3）對于擾動后產(chǎn)生的劣解不是一味地舍去，而是根據(jù)一定準(zhǔn)則接受（詳見3.4.3小節(jié)），以便跳出局優(yōu)。

本文主要研究帶軟時間窗的VRPSPD問題，設(shè)計了結(jié)合以VND為主體的適應(yīng)性擾動機(jī)制的混合PSO算法進(jìn)行求解，以車輛派遣成本、行駛成本和時間窗懲罰成本之和為總優(yōu)化目標(biāo)。實驗結(jié)果分析表明，算法具有良好的尋優(yōu)性能。

2 問題描述及數(shù)學(xué)模型

2.1 問題描述

VRPSPDSTW問題可以描述如下：多臺同質(zhì)車輛從物流中心出發(fā)，為一系列客戶點進(jìn)行服務(wù)，且每個客戶處均設(shè)有獨立的服務(wù)時間窗。軟時間窗約束允許車輛到達(dá)時間早于或晚于所設(shè)時間窗，并根據(jù)偏離時間窗的長短進(jìn)行懲罰。每輛車在訪問客戶點時，先卸貨后取貨，一次性滿足客戶的集配貨需求，并且在任何服務(wù)節(jié)點均不得超過車輛容量限制，在服務(wù)完所分配的客戶后，返回物流中心。要求在滿足所有客戶需求的前提下，使得車輛派遣成本、行駛成本和時間懲罰成本之和最小?？勺鋈缦录僭O(shè)：

（1）已知各客戶和物流中心的位置，且客戶處的送取貨量已知。

（2）每個客戶只能一輛車訪問一次，并且要求車輛同時滿足集配貨兩種需求。

（3）物流中心擁有多臺同質(zhì)車輛，即容量一致，且能夠滿足所有客戶的需求。

（4）車輛平均行駛速度已知且不變。

（5）客戶的服務(wù)時間窗已知。

2.2 問題描述模型建立

為了方便建模描述VRPSPDSTW問題，本文采用如下的符號體系。

Q：車輛最大容量；

Di：客戶i處的配貨量，非負(fù)值且小于Q；

Pi：客戶i處的集貨量，非負(fù)值且小于Q；

ETi：客戶i的時間窗下限；

ELi：客戶i的時間窗上限；

Tik：車輛k在客戶i的到達(dá)時間；

WTi：客戶i的等待時間，即WTi=max{ETi-Tik,0}；

si：客戶i的服務(wù)時間；

dij：車輛從客戶i到客戶 j的距離（采用歐幾里得距離得到）；

tij：車輛從客戶i到客戶 j的行駛時間；

n：總的客戶數(shù)，其中采用i、j表示客戶編號；

m：總的車輛數(shù)，其中采用k表示車輛編號；

cv：每輛車的派送成本；

cd：單位距離運(yùn)輸成本；

CTi：車輛早于時間窗到達(dá)時，單位時間機(jī)會損失成本；

CLi：車輛晚于時間窗到達(dá)時，單位時間懲罰成本；

yijk：表示經(jīng)由客戶i到客戶 j的車輛k，離開客戶i的車輛載重量；

xijk：為決策變量。若車輛k從客戶i直接到達(dá)客戶 j，xijk=1，否則為0。

建立相應(yīng)的數(shù)學(xué)模型如下：

式（1）為目標(biāo)函數(shù)，其中第一部分表示車輛派遣費用；第二部分表示運(yùn)輸費用；最后一部分是由早到或晚到時間懲罰成本構(gòu)成；式（2）確保派出的車輛數(shù)不超過最大車輛數(shù)；式（3）各節(jié)點車輛進(jìn)出平衡，當(dāng)i為0時，表示車輛從中心出發(fā)后又返回中心；式（4）確保每個客戶只被一輛車訪問，且只能訪問一次；式（5）為消除子回路約束；式（6）表示路徑上任意節(jié)點處的車輛載重均不超過最大載重量；式（7）表示車輛從中心出發(fā)時的載重量與所要配送的客戶的配貨量相等，且不超過Q；式（8）表示車輛返回物流中心，載重量為所服務(wù)客戶總的取貨量，且不超過Q；式（9）車輛在經(jīng)由前后客戶節(jié)點后的載重平衡，即車輛離開節(jié)點j的載重量等于車輛到達(dá)時的載重量去掉 j的配貨量，增加 j的取貨量；式（10）表示時間約束，其中M 為一個任意大的數(shù)；式（11）、（12）為變量屬性約束。

3 混合PSO算法設(shè)計

3.1 PSO基本原理

PSO算法最初是由Kennedy和Eberhart于1995年提出一種新的仿生進(jìn)化算法[31]，模擬鳥群覓食行為，通過鳥之間的協(xié)作和個體經(jīng)驗來引導(dǎo)鳥群的整體運(yùn)動。在PSO算法中，每個備選解抽象成一個無質(zhì)量的粒子，若干個粒子形成一個種群，搜索過程為每個粒子追隨兩個極值，一個是粒子本身歷史搜索的個體極值 pbest，另一個是全種群到目前為止搜索到的群體極值gbest。所有粒子都有一個根據(jù)目標(biāo)函數(shù)所確定的適應(yīng)值進(jìn)行評價，有相同維度的位置和速度向量，速度來決定各個粒子飛行的方向和距離，其數(shù)學(xué)描述如下：

設(shè)種群中每個粒子代表D維空間中的一個解，在第t代第i個粒子位置表示為，速度表示，該粒子第t代的個體最優(yōu)位置為 pbestt，種群第t代的種群最優(yōu)位置為gbestt，在第(t+1)代狀態(tài)粒子速度和位置更新公式如下：

其中，c1、c2表示加速因子，r表示[0，1]之間的隨機(jī)數(shù)，w表示慣性因子，其大小的設(shè)置影響著算法的收斂速度和結(jié)果，較大的w有利于跳出局部最優(yōu)解，有利于提高算法的全局搜索能力，較小的w有利于算法收斂，會增強(qiáng)算法的局部搜索能力。

3.2 粒子編碼

為了充分發(fā)揮PSO在連續(xù)空間優(yōu)化問題良好的搜索性能，本文采用實數(shù)編碼方式[25]，即對于擁有n個客戶、m輛車的問題來說，粒子的速度和位置均為n維實數(shù)向量，速度每一維的整數(shù)部分表示所用車輛編號，在同一車輛編號（即整數(shù)部分相同）的小數(shù)部分，大小順序?qū)?yīng)該客戶在此路徑中的訪問順序。故每個粒子的位置取值范圍為1～m.99，速度取值范圍為-(m.99-1)～(m.99-1)。為了保證粒子在可行范圍內(nèi)搜索，限定位置變化范圍為[xmin,xmax]，速度變化范圍為[vmin,vmax]。

3.3 初始化種群

首先根據(jù)3.2節(jié)編碼方式和設(shè)定范圍隨機(jī)生成粒子的速度和位置，然后其中的一半粒子根據(jù)最鄰近方法并轉(zhuǎn)換成實數(shù)編碼，產(chǎn)生新的位置，由于粒子的速度和位置是一一對應(yīng)關(guān)系，所以這一半粒子的速度是根據(jù)式（14）將前后兩個位置相減得到。

由于算法采用的是實數(shù)編碼，所以本文提出如下映射方法，將最鄰近方法產(chǎn)生的解映射為實數(shù)形式，即：在同一車輛訪問的節(jié)點集合中，假設(shè)共有n個節(jié)點，車輛編號為k，則對于該路徑上的第i個節(jié)點來說，其位置大小為i×0.99/(n+1)+k。

3.4 適應(yīng)性擾動機(jī)制APM（Adptive Perturbation Mechanism）

在PSO算法中，gbest記錄著當(dāng)前種群經(jīng)歷的最優(yōu)狀態(tài)，影響所有個體位置狀態(tài)的變化過程。若一味地接受更優(yōu)解，很有可能導(dǎo)致算法的搜索逐漸聚集在局優(yōu)附近，探索能力減弱。因此，為了改善上述問題，本文通過以VND為主體根據(jù)算法搜索進(jìn)程情況進(jìn)行適應(yīng)性擾動，不僅可以擴(kuò)展搜索的空間，而且可以根據(jù)設(shè)定的閾值接受劣解，以便跳出局優(yōu)。同時為了防止因擾動而丟失搜索過的最好解，將Gbest作為歷史搜索到的最優(yōu)解，gbest作為種群極值。

對于擾動機(jī)制來講，一個是擾動的時機(jī)，它可以有效地平衡PSO和VND兩者之間搜索優(yōu)勢的發(fā)揮，另一個是構(gòu)成擾動部分采用變鄰域結(jié)構(gòu)，這兩者都對擾動部分產(chǎn)生的效果起著關(guān)鍵性的作用。

3.4.1 線性啟動準(zhǔn)則

為了更好地發(fā)揮PSO的全局搜索能力，同時也為APM過程提供一個較好質(zhì)量的擾動解，故本文提出如下線性啟動準(zhǔn)則：

其中iter表示PSO中的迭代次數(shù)，noimp_max和noimp_min分別表示最大、最小允許歷史最優(yōu)解Gbest未改進(jìn)次數(shù)，若Gbest在連續(xù)多代未得到改進(jìn)則認(rèn)為算法陷入局部極值，所以只有當(dāng)Gbest的未改進(jìn)次數(shù)noimp_times大于viacriterion時，啟動APM。

3.4.2鄰域結(jié)構(gòu)

在進(jìn)行VND搜索之前，需要先定義一組鄰域結(jié)構(gòu)，本文采用了路徑間、路徑內(nèi)等常見的鄰域結(jié)構(gòu)來進(jìn)行擾動。

為了提高執(zhí)行操作的有效性，在路徑間鄰域變換中，首先隨機(jī)選擇一條路徑，通過計算此路徑與解中其他路徑重心之間的距離，并將各個距離與所有距離之和的比值作為選擇概率，根據(jù)輪盤賭規(guī)則選擇另一條路徑r2。

（1）Cross1：由于所選兩個路徑所含客戶數(shù)不同，所以分別隨機(jī)選擇交叉點。交叉點將路徑劃分為兩部分，路徑r1的第一部分與r2的第二部分結(jié)合形成新的路徑，同理r1、r2的剩余兩個部分結(jié)合成另一條新路徑（見圖1）。

圖1 Cross1鄰域結(jié)構(gòu)示意圖

（2）Cross2：此交叉操作與Cross1類似，不同之處在于路徑r1的第一部分與r2的第一部分結(jié)合形成新的路徑，同理r1、r2的剩余兩個部分結(jié)合成另一條新路徑，其中r1節(jié)點的順序不變，r2節(jié)點的順序要進(jìn)行逆轉(zhuǎn)（見圖2）。

圖2 Cross2鄰域結(jié)構(gòu)示意

（3）Shift（0，n）：在r1中任選一點或兩點插入到r2中，即n取1或2。

（4）Exchange(m,n)：在r1、r2中分別選取連續(xù)的m、n個節(jié)點，將r1中的m個節(jié)點和r2的n個節(jié)點進(jìn)行交換，在本文中，m取1或2，n取1。

（5）采用的路徑內(nèi)鄰域結(jié)構(gòu)如下，設(shè)解s中任選一條路徑r。

（6）Swap：在r中交換任意兩點的位置。

（7）插入：在r中任選一點或兩點插入到路徑中另外位置。

（8）2-opt：在r中任選兩點，將兩點間的節(jié)點順序逆轉(zhuǎn)。

3.4.3 適應(yīng)性選擇鄰域策略

本文采用適應(yīng)性選擇策略[32]來選擇路徑間某一鄰域結(jié)構(gòu)進(jìn)行搜索。與隨機(jī)選擇某一鄰域進(jìn)行變化相比，可以使搜索過較好解的鄰域結(jié)構(gòu)有更大機(jī)會被選中，減少搜索的盲目性，其內(nèi)容如下。

起初，每個路徑間鄰域結(jié)構(gòu)的權(quán)重賦予為1，擾動次數(shù)每隔50代進(jìn)行更新權(quán)重，并重置使用次數(shù)和所得分?jǐn)?shù)。根據(jù)如下公式進(jìn)行更新[32]：

σj,t計算準(zhǔn)則：在同一周期中，經(jīng)過路徑間鄰域結(jié)構(gòu)變化后，若其適應(yīng)度值優(yōu)于歷史搜索到的最優(yōu)解Gbest，則此鄰域結(jié)構(gòu)加30分，并設(shè)noimp_times為0，同時更新Gbest和gbest；若其所得解未進(jìn)一步改善Gbest，但卻優(yōu)于Gbest+Gbest×δ，δ為偏差系數(shù)，δ∈(0,1)，則加10分，并累加noimp_times，Gbest不變，但更新gbest，起到擾動PSO的作用；如果均不滿足上述兩種情況，而是獲得了滿足約束解，則加6分，并替換種群中最差的粒子；否則分?jǐn)?shù)不變。在本文所提算法中，δ取0.1。

3.4.4擾動搜索過程

本文設(shè)計的算法主要由路徑間VND和路徑內(nèi)VND組成，為了方便描述分別表示為inter-VND、intra-VND。擾動部分其過程如下。

擾動搜索步驟如下：

步驟1輸入初始解s，定義一組路徑間鄰域結(jié)構(gòu)Nm，其中m=1,2,…,mmax，路徑內(nèi)鄰域結(jié)構(gòu)Nn，其中n=1,2,…,nmax。

步驟2 Shaking：根據(jù)3.4.3小節(jié)中的策略選擇鄰域為N，在該鄰域多次循環(huán)，獲得其中的最小解即s1←N(s)，與 f()s進(jìn)行比較，若優(yōu)則進(jìn)行替換s。

步驟3 inter-VND其過程的偽代碼如下[9]：

1.m←1；

2.whilem＜mmax

3.inter←1；

4.whileinter＜mmax

5.s2←Nm(s)；/找到在鄰域Nm中的最好解/

6.if f(s2)＜f(s1)

7.s1←s2；

8.m←1；

9.else

10.m←m+1；

11.end

12.inter←inter+1；

13.end

14.end

步驟4 intra-VND。路徑內(nèi)的鄰域變化與上述步驟相同，其中不同的是輸入的初始解為s1，采用的是路徑內(nèi)鄰域變化結(jié)構(gòu)，以及所包含鄰域的個數(shù)為nmax。

步驟5根據(jù)3.4.3小節(jié)進(jìn)行選擇鄰域，結(jié)束擾動。

值得注意的是，搜尋解s在某種鄰域結(jié)構(gòu)的局優(yōu)時，較易忽略循環(huán)的設(shè)定次數(shù)[9]。所以為了更好地提高搜索效率，本文在此過程（即在上述步驟5）中加入自適應(yīng)性循環(huán)次數(shù)變化，可以根據(jù)在該種鄰域結(jié)構(gòu)搜索到解的質(zhì)量來調(diào)整循環(huán)次數(shù)。

假設(shè)當(dāng)前鄰域結(jié)構(gòu)為N，輸入的初始解為s，初設(shè)循環(huán)次數(shù) pmax，首先在該鄰域隨機(jī)獲得一個擾動解s′，由于鄰域變化本身帶有隨機(jī)性，所以在搜索過程中，獲得的鄰域解的質(zhì)量有所不同，其質(zhì)量的優(yōu)劣決定了在當(dāng)前鄰域搜索的次數(shù)多少。如果搜索到解s″優(yōu)于當(dāng)前s′，則代替當(dāng)前解，在此基礎(chǔ)上若滿足設(shè)定的條件，則增加循環(huán)次數(shù) pmax+1。如果不優(yōu)于解s′，則累計循環(huán)計數(shù)器 p，若在此基礎(chǔ)上不滿足設(shè)定的條件，則減少循環(huán)次數(shù) pmax-1，其中設(shè)定的條件為滿足約束。同時，為了進(jìn)一步控制循環(huán)的進(jìn)程，計算在同一鄰域中搜索到可行或不可行解連續(xù)出現(xiàn)的次數(shù)，若不小于ε，則變化循環(huán)次數(shù)。根據(jù)先前的實驗，ε設(shè)為3。主要步驟如下：

1.s'←N(s);p←1;mm=0；

2.while p＜pmax

3.s''←N(s)；

4.iff(s'')＜f(s')

5.s'←s''；

6.ifs'滿足條件 /可以設(shè)定為滿足可行約束/

7.pmax←pmax+1；

8.mm←0；

9.end

10.else

11.p←p+1

12.ifs'不滿足設(shè)定的條件

13.mm=mm+1；

14.ifmm＞=ε]

15.pmax←pmax-1；

16.end

17.end

18.end

19.end

3.5 適應(yīng)性權(quán)重

在PSO部分，本文采用如下適應(yīng)權(quán)重，其更新公式如下，設(shè)權(quán)重最小值wmin，最大值為wmax，第i個粒子的適應(yīng)值為 fi，種群最小值為 fmin，平均值為 favg，對于粒子i的權(quán)重值計算如下：

3.6 算法求解過程

算法求解過程如下，其流程圖見圖3。

圖3 混合粒子群算法流程圖

步驟1采用上述3.3節(jié)方法產(chǎn)生初始粒子的位置和速度，計算適應(yīng)度值并初始化種群的個體 pbest和種群的gbest、Gbest，以及最大迭代次數(shù)、加速因子等相關(guān)參數(shù)。

步驟2計算粒子的適應(yīng)性權(quán)重（公式（17）），并根據(jù)公式（13）和（14）更新粒子位置和速度，計算適應(yīng)值。

步驟3若優(yōu)于粒子個體最優(yōu)適應(yīng)值 pbest或種群最優(yōu)適應(yīng)值gbest，則相應(yīng)地更新 pbest或gbest。

步驟4種群所有粒子全部更新完以后，若gbest優(yōu)于歷史最優(yōu)Gbest，則更新Gbest，設(shè)noimp_times為0，否則累加未改進(jìn)的次數(shù)。

步驟5判斷是否滿足擾動準(zhǔn)則，是則依3.4.4小節(jié)方法進(jìn)行擾動搜索，進(jìn)行下一步，否則轉(zhuǎn)至步驟8。

步驟6對擾動后的解依3.4.3小節(jié)進(jìn)行評分并記錄鄰域使用次數(shù)，以及相應(yīng)地進(jìn)行更新。

步驟7判斷是否達(dá)到下一評分周期，若達(dá)到則更新權(quán)重，并重置分?jǐn)?shù)和使用次數(shù)。

步驟8若迭代次數(shù)達(dá)到最大允許迭代次數(shù)且未改進(jìn)次數(shù)達(dá)到終止的最大次數(shù)，則終止算法，否則返回步驟2。

4 算例驗證與結(jié)果分析

為了驗證所提算法的可行性和有效性，本文進(jìn)行了3組實驗。第1、2組實驗分別選取了文獻(xiàn)[13]和文獻(xiàn)[21]中的算例進(jìn)行硬時間窗問題的測試，第3組實驗選取文獻(xiàn)[20]中求解帶軟時間窗VRPSPD問題的算例測試。本文所設(shè)計算法在Matlab 2015a平臺上編譯，并在PC機(jī)上運(yùn)行實驗（測試環(huán)境CPU：雙核奔騰2.80 GHz；內(nèi)存：2 GB；操作系統(tǒng)：Windows 7）。

實驗1實驗選取由Wang和Chen[13]生成求解VRPSPDTW的算例，分別與遺傳算法[13]、模擬退火[15]、平行模擬退火[16]以及全局人工魚群算法[17]進(jìn)行比較。

在這組實驗中不同規(guī)模算例所取參數(shù)如下：W_min和W_max分別表示慣性權(quán)重的最大、最小值，Popsize表示種群規(guī)模，Itermax表示最大迭代次數(shù)（見表1）。加速因子c1=c2=1.496 18，在擾動部分中，Shaking階段在同一鄰域循環(huán)次數(shù)為3倍的Popsize，pmax在路徑間、路徑內(nèi)分別設(shè)為2倍的Popsize、Popsize。值得注意的是，50個客戶規(guī)模的算例中，為了充分地進(jìn)行擾動，shaking階段采用的鄰域有所不同，采用的是doublereplace（1，1），triple shift（1，0），triple swap（1，1），double cross[30]以及exchange（m，n）（其中 m 和 n 的組合為（3，2）、（3，1）和（2，2））。

表1 不同規(guī)模算例所用參數(shù)

表2列出了本文所提算法運(yùn)行10次獲得的最好解，并與其他算法在總運(yùn)輸距離的比較。其中將Wang和Chen[13]采用GA計算的結(jié)果作為標(biāo)準(zhǔn)1，其他算法的結(jié)果進(jìn)行“標(biāo)準(zhǔn)化”，即第二至五列均為與GA的比值，最后一列給出了所提算法計算得到的結(jié)果，最后一行列出了各算例的比值平均值。由表2可知，除算例Rcdp5004外，其他算例均能找到當(dāng)前已知最好解或更優(yōu)解，從比值的平均值來看，本文所提算法為0.937 5，分別與SAA[15]、p-SA[16]、IGAFSA[17]相比都較小，其中 Rcdp1007改進(jìn)最大。在9個算例中，更新了6個算例的最好解。在求解50個客戶規(guī)模的算例運(yùn)行時間均在820 s以內(nèi)獲得，雖然運(yùn)行時間較長（采用了多鄰域循環(huán)所致），但是仍在可以接受的范圍之內(nèi)。綜上，可以驗證所提算法的可行性和有效性。

表2 各算法運(yùn)行結(jié)果比較

Rcdp1001軌跡圖如圖4；Rcdp2504軌跡圖如圖5。

圖4 Rcdp1001軌跡圖

圖5 Rcdp2504軌跡圖

實驗2根據(jù)文獻(xiàn)[21]擴(kuò)展的Solomon[14]算例的方法，生成配送和取貨需求，在RC類算例中形成客戶規(guī)模分別為10、25和50共9個算例。該組實驗采用的相關(guān)參數(shù)與實驗1相同。表3為本文算法與文獻(xiàn)[21]運(yùn)行結(jié)果的比較，從表中的結(jié)果可以得到，本文算法的節(jié)省車輛數(shù)（NV）和總行駛距離（TD）結(jié)果都優(yōu)于離散粒子群算法DPSO[21]（Discrete Particle Swarm Optimization），在算例RCdp25107的比較結(jié)果可知，雖然所用車輛數(shù)比DPSO多1輛，但是總行駛距離改進(jìn)了29.83%。在各算例運(yùn)行時間均在475 s之內(nèi)，雖然相比較而言，耗時較長，但仍在可接受范圍內(nèi)。綜上表明，本文算法優(yōu)勢明顯，同時也說明不同結(jié)合方式很大程度上影響著算法的性能。

表3 與離散粒子群算法結(jié)果比較

實驗3本文根據(jù)文獻(xiàn)[20]中，包含1個物流中心和20個客戶點的算例，用來求解帶有軟時間窗的VRPSPD問題，通過計算總需求量與車輛載重的比值可知，最少需要車輛數(shù)為3，為了保證產(chǎn)生的粒子的多樣性且有效地減少搜索的盲目性，本文設(shè)定車輛數(shù)為4。該實驗所選用的參數(shù)與實驗1中25個客戶規(guī)模的算例相同，在模型中的參數(shù)與文獻(xiàn)[20]中相同。

表4列出的是加入適應(yīng)性權(quán)重的基本PSO（用BPSO表示）、SA[20]和本文所提算法運(yùn)行10次的結(jié)果，從中可以看出本文算法在10次運(yùn)行中均能找到較高質(zhì)量的解。從表中得知，與BPSO相比，在平均值、最好值和最差值上分別改進(jìn)了12.84% 、17.91%、14.69%，說明在PSO的基礎(chǔ)上，所提算法可以有效地跳出局部最優(yōu)；與文獻(xiàn)[20]的SA比較，分別改進(jìn)了5.08%、12.95%、4.05%。本文算法所獲的最好解其總配送成本為175.2元，總運(yùn)輸距離為113.1 km，所用總的車輛數(shù)為3，與前兩個算法所用車輛數(shù)均少1。具體的配送路線為：第1輛車：0-14-13-5-4-19-6-12-0；第2輛車：0-8-9-18-10-11-3-15-0；第3輛車：0-7-20-1-17-16-2-0。算法在10次運(yùn)行的平均耗時為160.9 s，雖然算法耗時較多，但仍在可接受的范圍內(nèi)。

表4 各算法運(yùn)行結(jié)果比較

如圖6為本文算法運(yùn)行最好結(jié)果（即第10次實驗）總配送成本和運(yùn)輸成本的迭代變化圖；如圖7是其中擾動種群全局極值的變化圖。從兩者比較來看，在1 000代左右以及1 600代到1 900代這兩個階段的擾動，均能有效地引導(dǎo)算法跳出局優(yōu)；在2 000代以后算法最終收斂。

圖6 迭代變化圖

圖7 擾動圖

綜上，擾動部分可以更進(jìn)一步優(yōu)化PSO的搜索進(jìn)程，起到改進(jìn)的作用，同時也驗證了算法的可行性和有效性。

5 結(jié)論語

本文針對更貼近現(xiàn)實物流配送網(wǎng)絡(luò)的帶軟時間窗VRPSPD問題，建立了以車輛派遣成本、行駛成本和時間窗懲罰成本之和最小為目標(biāo)的車輛路徑優(yōu)化模型，設(shè)計了混合PSO算法進(jìn)行求解。為了克服PSO容易陷入局部最優(yōu)，爬山能力弱等問題，在Gbest未改進(jìn)次數(shù)達(dá)到線性擾動啟動準(zhǔn)則計算的值時，啟動適應(yīng)性擾動機(jī)制，即進(jìn)行VND擾動，其中在Shaking階段采取了適應(yīng)性選擇鄰域策略，并在每個鄰域搜索中應(yīng)用可變的循環(huán)次數(shù)，以便提高對解空間的探測能力和搜索效率。

通過與現(xiàn)有算法求解結(jié)果比較，本文所設(shè)計的混合PSO能夠有效地求解該問題，同時也表明了擾動機(jī)制可以很好地跳出局部最優(yōu)，并在可接受的時間范圍內(nèi)獲得滿意解，驗證了算法的有效性。今后，將針對不同算法之間的結(jié)合方式，以及求解更一般化的車輛路徑問題進(jìn)行進(jìn)一步研究。