帶時間窗的多車型需求可拆分攬收配送問題

張 凱，靳 鵬，崔 勇

合肥工業大學 管理學院，合肥230009

傳統物流企業利用自有車輛為客戶提供攬收/配送服務。由于自有車輛具有使用成本高、管理難度大等缺點，在共享經濟大背景下，物流企業開始采用一種基于外部車輛的新型攬收配送模式，利用外部車輛為客戶提供攬收/配送服務，外部車輛主要是社會閑置車輛。新型攬收配送模式可以減少企業自有車輛的使用，降低物流成本，并對提高社會閑置車輛利用率和環境保護方面具有重要意義。

本文研究了基于新型攬收配送模式下的一種帶時間窗的多車型需求可拆分攬收配送問題（Multi-Vehicle Split Pickup and Delivery Problem with Time Windows，MVSPDPTW），使用社會閑置車輛作為攬收/配送任務的資源，社會閑置車輛擁有自己的起點、終點以及最早出發時間和最晚返回時間，車輛必須在最早出發時間后從起點出發為客戶提供攬收/配送服務，并在最晚返回時間前返回終點。在服務客戶的過程中，車輛需要在客戶時間窗內提供服務，并且允許需求拆分，同一客戶需求可由多輛車完成，以便充分利用車輛的裝載能力。

車輛路徑問題（Vehicle Routing Problems，VRP）[1]自1959 年被提出以來，一直是物流研究領域的熱點問題。目前，在VRP 的研究中一般限制客戶的需求只能由一輛車來完成，屬于需求不可拆分的VRP。文獻[2]在1989年首次研究了需求可拆分的VRP（Split Delivery Vehicle Routing Problems，SDVRP），指出對客戶需求進行合理的拆分有利于提高車輛裝載率和節約總成本，并在1994 年證明SDVRP 是NP-hard[3]問題。目前已有不少學者從不同角度對SDVRP展開了研究。由于客戶通常要求在給定的時間內得到服務，因此一種考慮時間窗約束的SDVRP[4-7]得到了研究，文獻[4]首次將時間限制引入到SDVRP，構造了考慮時間窗約束的SDVRP模型，提出了一種兩階段啟發式算法進行求解。文獻[5]提出一種禁忌搜索算法求解考慮時間窗約束的SDVRP，設計了多種新型鄰域搜索算子。文獻[6-7]研究了考慮軟時間窗約束的SDVRP，軟時間窗約束通過添加懲罰函數允許客戶在給定的時間窗外被服務。在SDVRP中只考慮了客戶的配送需求，近些年隨著逆向物流的發展，客戶往往同時含有攬收/配送需求，因此一種考慮攬收需求的SDVRP[8-9]得到了研究，文獻[8]首次將需求可拆分引入攬收配送問題，揭示了需求可拆分帶來的好處，證明對于給定的一組攬收和配送位置，當客戶需求量略大于車輛載重量一半時可以帶來最大的收益。文獻[9]提出了一個改進蟻群算法求解需求可拆分的攬收配送問題，實驗結果表明該算法可以提高配送中心的運作效率。

以往對SDVRP 的研究基本建立在單一車型上，但是在實際配送中，車輛往往擁有多種車型[10-12]。文獻[10]研究了多車型的SDVRP，車輛的裝載能力與最大行駛里程各不相同，提出了一種模擬退火算法進行求解。文獻[11]研究需求可拆分的多車型開放式車輛路徑問題，車輛從不同的配送中心出發為客戶服務，在服務完最后一個客戶后結束服務而不需要返回配送中心，提出了一種禁忌搜索算法進行求解。文獻[12]分別研究了單車型和多車型的SDVRP，通過實驗結果說明了不同裝載能力的車輛更符合實際需求并且可以有效提高提高車輛的裝載率。

上述多數文獻大都將攬收需求、配送需求、多車型、客戶時間窗等單獨研究，但是隨著物流技術的快速發展，用戶需求更加復雜，這種單一模式研究已經難以滿足實際應用需求。因此，本文綜合考慮上述多種要素，研究了一個帶時間窗的多車型需求可拆分攬收配送問題，針對這個問題本文以執行任務車輛行駛路徑總長度最小為目標函數，建立了一個混合整數線性規劃模型，提出了一種高效禁忌模擬退火（Tabu Simulated Annealing，TSA）算法進行求解。

1 問題描述

在MVSPDPTW中存在給定的一個配送中心，一組地理位置不同的客戶和一組多車型車輛，MVSPDPTW的目標是選擇合適的車輛執行完所有客戶任務，使得執行任務車輛行駛路徑總長度最小。

0表示配送中心，T={1,2,…,n}表示客戶點集合，每個客戶點對應一個攬收/配送任務，將對應攬收任務的客戶點稱為pickup 點，對應配送任務的客戶點稱為delivery點。對于客戶點i∈T，[ei,li] 表示i的時間窗，ei表示最早開始服務時間，li表示最遲開始服務時間。qi表示i需要攬收/配送的貨物量。每個客戶點可以被多個車輛訪問，但只允許被同一車輛至多訪問一次。K={1,2,…,m}表示車輛集合，對于k∈K，sk、nk、Qk、ok、fk分別表示車輛k的起點、終點、最大容量、最早出發時間以及最晚返回時間。如圖1所示，車輛需要在最早出發時間后從起點出發，在pickup點攬收貨物送回配送中心，在配送中心裝載貨物配送至delivery點，最終在最晚返回時間前返回終點。車輛執行任務的過程中允許多次訪問配送中心。

圖1 車輛執行任務示意圖

2 數學模型

2.1 虛擬化處理

由于車輛可能多次訪問配送中心，如果模型中以一個變量去區分車輛第幾次訪問配送中心，則建模難度和求解難度都會顯著增加。為了降低難度，采用了以下處理方式：T={1,2,…,n}表示客戶點集合，n為實際客戶點的數量，對于編號為i的pickup 點，構建一個編號為n+i的虛擬delivery點與之對應；對于編號為j的delivery點，則重新編號為n+j，并構建一個編號為j的虛擬pickup 點與之對應。經過上述處理后，T=TP?TD，TP={1,2,…,n} 表示pickup 點集合，TD={n+1,n+2,…,2n}表示delivery點集合，車輛需要從點i(i∈TP)攬收貨物配送至點n+i。該種處理方式將配送中心替換成了n個虛擬客戶點，所有虛擬客戶點的位置與配送中心位置相同且無服務時間窗限制；每個虛擬客戶點攬收/配送的貨物量與對應實際客戶點相等。如圖2 所示，圖1中的配送中心經過虛擬化處理后替換成了方框內的虛擬客戶點，其中i*=n+i。

圖2 虛擬化示意圖

2.2 模型符號

經過上述處理后，將MVSPDPTW在全有向圖上圖上定義G=(V,A),V={S?N?T?0}表示頂點集，A={(i,j):i,j∈V,i≠j}表示弧集。

S：車輛的起點集合，S={s1,s2,…,sm}；

N：車輛的終點集合，N={n1,n2,…,nm}；

T：客戶點集合，T={1,2,…,2n}；

dij：點i和點j之間的距離，i,j∈V；

tij：點i和點j之間的時間，i,j∈V；

bik：車輛k到達點i的時間，k∈K,i∈T?0；

aik：車輛k在點i的開始服務時間，k∈K,i∈T?0；

rik：車輛k在點i的服務時長，k∈K,i∈T?0；

qik：車輛k在點i的裝載/卸載的貨物量，k∈K,i∈T?0；

wik：車輛k經過點i后承載貨物容量，k∈K,i∈T?0；

xijk：表示車輛k是否直接從點i到點j，如果是則為1，否則為0；

M：表示非常大的整數。

2.3 線性模型

其中，公式（1）為目標函數，使得執行任務車輛行駛路徑總長度最小；公式（2）、（3）表示車輛從起點出執行任務并最終返回終點；公式（4）、（5）表示每個客戶點可以被多個車輛訪問，但只允許被同一車輛至多訪問一次；公式（6）、（7）、（8）表示車輛需要從pickup點攬收貨物配送至對應的delivery點；公式（9）表示車輛從起點出發以及返回終點的承載貨物容量均為0；公式（10）表示容量約束；公式（11）、（12）表示車輛在客戶點承載貨物容量變化；公式（13）表示每個客戶點的分割貨物量之和等于該客戶點的總貨物量；公式（14）、（15）表示車輛需要在最早出發時間后從起點出發為客戶提供服務，并在最晚返回時間前返回終點；公式（16）、（17）表示車輛到達客戶點時間和在客戶點開始服務時間之間的關系；公式（18）表示車輛只能在客戶時間窗內提供服務。

3 算法設計

為了解決MVSPDPTW，提出了一種禁忌模擬退火算法，設計了兩種新的鄰域搜索算子，分別用于修復違反車輛容量約束以及換車操作。在算法內加入禁忌機制以及違反約束懲罰機制，實現了搜索空間的有效裁剪，提高了算法的全局尋優能力。首先利用動態貪婪算法生成初始可行解，并將初始可行解作為禁忌模擬退火算法的輸入，最終輸出全局最好解。

3.1 初始解生成

模擬退火算法本質上是一種迭代求解算法，初始解決定了迭代的起點，因此一個好的初始解有助于算法找到較好的最終解，本文參考文獻[5]中初始解的生成方法設計了一種簡單有效的動態貪婪算法生成初始解。

公式（20）、（21）分別表示車輛以及客戶點的選擇準則，車輛的選擇以服務時間最大化為原則。客戶點的選擇以當前規劃路徑中最后的訪問對象i和待選擇客戶點j之間總的旅行時間和等待時間最小化為原則。利用公式（21）為當前車輛依次選擇服務的客戶點，檢驗所選客戶點是否滿足插入約束，包括時間窗約束、容量約束、客戶點訪問次數約束以及需求約束，其中需求約束是指所有配送的貨物來源于配送中心，所有攬收的貨物送回配送中心。如果滿足約束則將客戶點插入到當前車輛規劃路徑中，重復這個過程直到所有的客戶需求被滿足。動態貪婪算法的具體流程如下：

輸入：車輛集合K、客戶點集合T、配送中心0。

輸出：初始可行解Sinitial。

步驟1 如果集合T不為空集，轉步驟2，否則，算法結束輸出結果。

步驟2 利用公式（20）從K中選出當前車輛k。

步驟3 構造待選擇客戶點集合C=T。

步驟4 如果C不為空集，轉步驟5，否則，轉步驟9。

步驟5 利用公式（21）從C中選出客戶j。

步驟6 檢驗j是否滿足插入約束，如果滿足則將j插入到當前車輛規劃路徑Rk中。

步驟7 檢驗j的需求是否全部被滿足，如果滿足則從T中刪除j。

步驟8 將j從C中刪除，轉步驟4。

步驟9 檢驗Rk最后一個訪問對象是否為配送中心0，如果是，轉步驟11，否則，轉步驟10。

步驟10 檢驗0是否滿足時間窗約束，如果滿足將0插入到Rk中，轉步驟3，否則，轉步驟11。

步驟11 從K中刪除k，轉步驟1。

3.2 違反約束懲罰機制

在迭代過程中接受不可行解有助于增大鄰域搜索范圍以及提高算法跳出局部最優的能力，幫助算法在迭代步驟[13]更好的可行解，因此本文設計了一種違反約束懲罰機制。整體約束懲罰函數設計如公式（22），B表示路徑規劃解的適應度值，β表示參數懲罰權重，η(i,k)、φ(i,k)、λ(i,k)分別表示Rk中點i容量約束評價值、時間窗約束評價值以及需求約束評價值，它們的計算方式由公式（23）、（24）、（25）給出。Rk={ok,…,nk}表示車輛k的執行任務路徑，它按照車輛的訪問順序排列而成。

3.3 禁忌機制

禁忌搜索是一種模擬人類思維過程的元啟發式算法，它的基本思想是避免對已搜索過的空間進行重復搜索，迫使算法去探索其他新的解決方案。

將每次迭代的當前解[14]作為禁忌對象存儲在禁忌表中，在以后的l（禁忌表的長度）次迭代中，該對象被設置成禁忌狀態，這與在禁忌結構中常用的邊移動禁忌[15-16]不同。禁忌表的長度指禁忌表儲存對象的個數，當儲存對象個數達到上限，再有對象加入時，則采用先進先出的原則將前面的對象從禁忌表中移除，再進行加入。

3.4 多鄰域結構

多鄰域結構可以增加解的多樣性，提高算法找到全局最好解的可能性。本文采用了五種鄰域搜索算子，其中遷移算子、遷移分割算子、交換算子來源于經典鄰域搜索算子[17]的引用，修復算子和換車算子針對本文在迭代過程中接受不可行解以及多車型的特點而設計，分別用于修復違反容量約束和換車操作。

為了更好地解釋算子操作，定義了以下符號，Rk(k∈K)和Rh(h∈K)分別表示車輛k和車輛h的執行任務路徑，對于客戶點分別表示在算子操作前后車輛k在客戶i裝載/卸載的貨物量。公式（26）給出了B(Rk)的計算方式，B(Rk)表示路徑Rk的適應度值。

任務集：在已規劃好的車輛路徑中，以車輛的起點、終點以及配送中心等非客戶點為界，兩個相鄰的非客戶點之間連續的客戶點組成一個任務集。如果兩個相鄰的非客戶點之間不包含客戶點，則不構成任務集。如圖3所示，客戶點1和2組成一個任務集，終點與配送中心之間不含客戶點則不構成任務集。

圖3 任務集示意圖

遷移算子：遷移算子用于將客戶點從原路徑中刪除插入到其他路徑中，如圖4所示，對于客戶點i∈Rk，將i從Rk中刪除，以B(Rh)最小化為原則將i插入到Rh中。

圖4 遷移算子操作示意圖

遷移分割算子：如圖5所示，客戶i由車輛k和h共同服務，客戶j由車輛h服務，i和j屬于路徑Rh中同一個任務集。將i從路徑Rk中刪除，然后以B(Rk)最小化原則將j插入到Rk中，此時i由車輛h服務，j由車輛k和h共同服務，k和h在i和j裝載/卸載的貨物量調整為。

圖5 遷移分割算子操作示意圖

交換算子：交換算子用于交換兩條路徑中的客戶點，如圖6 所示，客戶點i∈Rk，客戶點j∈Rh，將i和j進行交換，交換后，i∈Rh,j∈Rk。

圖6 交換算子操作示意圖

修復算子：修復算子用于修復違反車輛容量約束，如圖7 所示，假定車輛k在服務客戶i時承載貨物容量超過了車輛的最大容量Qk，此時在i所屬任務集中任選一個客戶點j,j可以為i，任選一條非空路徑Rh，將j以B(Rh)最小化原則插入到Rh中，k和h在j裝載/卸載的貨物量調整為：的計算方式由公式（23）給出。

圖7 修復算子操作示意圖

換車算子：對于非空路徑Rk以及空路徑Rh，將Rk中的所有客戶點構成客戶點集合C，然后按照初始解中路徑的生成方式，將C中的客戶點插入到Rh中。如果C中所有的客戶點都成功插入到Rh中，那么路徑Rk清空，否則的話將插入到Rh的客戶點從Rk中刪除。

圖8 換車算子操作示意圖

由于單個客戶點只允許被同一車輛至多訪問一次，所以在算子操作后，如果出現如圖9所示情況，Rk中含有兩個客戶點i，則以B(Rk)最小化原則選擇一個i從Rk中刪除，將k兩次服務i裝載/卸載的貨物量合并。

圖9 合并客戶點示意圖

3.5 禁忌模擬退火算法描述

首先利用動態貪婪算法生成初始可行解，并將其設置為第0代的當前解Snow和全局最好解Sbest，然后開始迭代，在每次迭代中，利用設計的準則更新Snow和Sbest。最后當當前溫度小于設置的最低溫度時，算法終止輸出Sbest，具體算法流程如下：

輸入：初始溫度T0、終止溫度Tend、冷卻速率RT、迭代步長L、鄰域解數量Nu、初始可行解Sinitial。

輸出：全局最好解Sbest。

步驟1 將Sinitial設置為Snow以及Sbest。

步驟2 初始化當前溫度T'=T0，初始化當前溫度下迭代次數L'=0。

步驟3 初始化禁忌表，將Snow加入到禁忌表中。

步驟4 如果T'>Tend，轉步驟5，否則，算法結束輸出結果。

步驟5 如果L'

步驟6 從五個鄰域算子中隨機選擇一個鄰域算子轉換當前解生成鄰域解。

步驟7 重復步驟6直到鄰域解的數量達到Nu。

步驟8 從Nu個鄰域解中選出最優的沒有被禁忌的鄰域解Scan。

步驟9 如果Scan

步驟10 如果>Random(0,1),Scan將成為新的Snow,Random(0,1)用于生成（0，1）上均勻分布的隨機數據，否則保留Snow。

步驟11 如果Scan為可行解并且滿足Scan

步驟12 如果Snow發生變化，則將Snow加入禁忌表。

步驟13 轉步驟5。

4 仿真算例

本文利用文獻[18]使用的考慮時間窗約束的VRP（Vehicle Routing Problem with Time Windows，VRPTW）數據集、文獻[19]使用的考慮時間窗約束的VPRSD（Split Delivery Vehicle Routing Problem with Time Windows，SDVRPTW）數據集以及構造的MVSPDPTW 仿真數據集對TSA 算法進行了測試并對實驗結果進行對比分析，每個算例都運行20 次。所有實驗均在i5-8500 CPU 3.0 GHz、內存16 GB 的臺式計算機上、myeclipse2017的環境下進行。

4.1 參數設置方法

通過設置合理的參數可以有效提高算法的性能[20]，首先參考文獻[21-22]中的參數設置確定本文各個參數的測試區間，然后選取不同規模算例，通過更改參數組合進行多組測試，最后確定參數的最佳組合。

各參數取值范圍以及步長具體參考表1，為簡化計算復雜度，本文采用控制變量法依次確定各參數的最佳值。TSA 參數設置的最終結果如下所示：初始溫度T0=100、終止溫度Tend=0.1、迭代長度L=30、降溫速率RT=0.93、懲罰權重β=4、鄰域解數量Nu=2n+150（n為客戶數量）、禁忌表的長度l=8。

表1 算法參數組合

4.2 基于VRPTW、SDVRPTW數據集的數值實驗

由于目前在MVSPDPTW 領域沒有公認的數據集可供使用，導致沒有辦法對本文算法直接進行測試，MVSPDPTW 是經典VRPTW 以及SDVRPTW 的拓展，問題中很多約束條件相似，因此如果本文算法可以很好地解決VRPTW 以及SDVRPTW，那么有理由相信本文算法也可以有效地求解MVSPDPTW。

4.2.1 基于VRPTW數據集的數值實驗

從文獻[18]使用的VRPTW數據集中有針對性地選擇了27 個VRPTW 算例進行測試，VRPTW 算例中車輛容量皆為200。表2記錄了TSA算法求解每一個算例結果Rbest和算例已知最好解Kbest對比，Gap表示Rbest與Kbest之間差距，Gap=(Rbest-Kbest)/Kbest。Kbest數據來源于http：//web.cba.neu.edu/～msolomon/problems.htm。

表2 TSA算法求解結果與已知最好解對比

根據表2中數據可知，TSA算法更新了4個VRPTW算例（算例編號C101，客戶數量25；算例編號C102，客戶數量50；算例編號C107，客戶數量50；算例編號C109，客戶數量50）的已知最好解。在剩余23個VRPTW算例中，TSA 算法的求解結果雖然劣于已知最好解，但是它們之間的差距均保持在1%之內。實驗結果說明了TSA算法求解VRPTW效果好。

4.2.2 基于SDVRPTW數據集的數值實驗

從文獻[19]使用的SDVRPTW 數據集中有針對性地選擇了12 個SDVRPTW 算例進行測試，SDVRPTW算例中車輛容量皆為100。表3記錄了TSA算法和分支與價格和削減[19（]Branch and Price and Cut algorithm，BPC）算法求解12 個SDVRPTW 算例結果對比。Gap表示TSA 算法求解每一個算例結果(Rbest)和BPC 算法求解每一個算例結果(Cbest)之間的差距，Gap=(Rbest-Cbest)/Cbest。

從表3中數據可知，TSA算法在客戶數量為25的小規模算例下的求解結果與BPC算法相同。隨著客戶數量的增大，雖然TSA算法的求解結果劣于BPC算法，但是它們之間的差距均保持1%以內。在求解時間方面，在小規模算例下，BPC算法求解時間小于TSA算法，但是隨著客戶數量的增大，BPC算法求解時間遠超TSA算法，并且在相同客戶數量下不同算例的求解時間上，TSA 算法基本保持不變而BPC 算法則差異明顯。實驗結果充分說明了TSA 算法求解SDVRPTW 效果接近BPC算法。

表3 TSA算法與BPC算法求解結果對比

4.3 基于仿真數據集的數值實驗

基于4.2.2 節的12 個SDVRPTW 算例構建12 個MVSPSPTW 仿真算例，不改變SDVRPTW 算例中配送中心位置、客戶位置、客戶貨物量以及客戶時間窗，僅隨機選擇半數客戶將其配送任務改為攬收任務，對于車輛的起點和終點不再設置為配送中心，而是在客戶所屬100×100的正方形區域內隨機生成；車輛的開始服務時間為0～800間隨機生成的整數，結束服務時間等于500～900間隨機生成的整數加上車輛的開始服務時間；車輛最大容量為70～120 間隨機生成的整數。此外，當客戶數量為25時隨機生成15輛車，當客戶數量為50時隨機生成30輛車，保證車輛資源充足。

表4 記錄了TSA 算法、傳統模擬退火（Simulated Annealing，SA）算法、粒子群算法（Particles Swarm Optimization，PSO）求解MVSPDPTW 仿真算例結果對比，PSO算法的設計參考了文獻[23-24]中的處理方式，傳統SA 算法與PSO 算法的參數通過4.1 節參數設置方法來確定。Gap1、Gap2 分別表示TSA算法求解每一個算例結果(Rbest)與SA算法求解每一個算例結果(Tbest)、PSO算法求解每一個算例結果(Pbest)之間的差距，Gap1=(Tbest-Rbest)/Rbest，Gap2=(Pbest-Rbest)/Rbest。

根據表4中數據可知，本文所提出的TSA算法對比傳統SA 算法和PSO 算法求解性能更優。主要體現在：（1）TSA算法的求解結果最優，相比傳統SA算法最高提升了2.40%，平均提升了1.14%；相比于PSO算法最高提升了5.19%，平均提升了3.52%。（2）在求解時間方面，SA 算法最長，TSA 算法次之，PSO 算法最短，但TSA 算法與PSO算法耗時總體差距不大。此外，TSA算法求解相同客戶數量算例的時間基本保持不變，具有較好的穩定性。TSA 算法求解性能更優的原因主要在于：（1）由于本文多車型的特點，因此選車結果會對最終解的優劣產生至關重要的影響，換車算子有利于尋找到最合適的車輛執行任務。除此之外多種算子配合的方式有利于擴大鄰域搜索范圍，提高算法的全局尋優能力。（2）TSA算法內加入了禁忌機制以及違反約束懲罰機制，實現了搜索空間的有效裁剪，避免算法陷入局部最優的同時節約了求解時間，相比之下傳統SA 算法和PSO 算法易陷入局部最優。

表4 TSA、SA、PSO算法求解結果對比

5 總結展望

本文研究了基于新型攬收配送模式下的一種帶時間窗的多車型需求可拆分攬收配送問題。針對這個問題建立一個混合整數線性規劃模型。在構造模型的過程中，通過為每個客戶點構建一個對應虛擬攬收/配送點的方式減少模型中部分變量的增加，并使用大M法對非線性約束做了線性化處理。設計了一種高效TSA 算法，TSA 算法在傳統SA 算法框架內加入了禁忌機制以及違反約束懲罰機制，實現了搜索空間的有效裁剪，提高了算法的全局尋優能力。TSA 算法更新了27 個VRPTW算例中4個算例的已知最好解；12個SDVRPTW算例的求解結果中，TSA算法與BPC算法的差距均不超過1%，且求解時間更短；就12個MVSPDPTW仿真算例而言，TSA 求解結果均不差于傳統SA 算法以及PSO 算法，且分別最高提升了2.40%和5.19%。實驗結果充分證明了本文所提出的TSA 算法求解MVSPDPTW 的可行性和有效性。

本文所研究的內容依舊有許多亟待解決的問題，將來的研究工作可以從以下兩個方面展開：（1）考慮需求依訂單拆分而不是任意連續拆分；（2）考慮使用外部車輛和自有車輛共同配送方式。