熊興隆，劉 佳，李 猛，馬愈昭

1.中國民航大學(xué) 天津市智能信號與圖像處理重點(diǎn)實(shí)驗室，天津 300300

2.中國民航大學(xué) 電子信息與自動化學(xué)院，天津 300300

3.中國民航大學(xué) 民航空管研究院，天津 300300

在物流配送領(lǐng)域，傳統(tǒng)單配送車送貨存在配送成本高、效率低、延誤高、環(huán)境污染等問題。隨著無人機(jī)在物流領(lǐng)域受到越來越多的關(guān)注，提出一種新型交付概念——無人機(jī)聯(lián)合配送車協(xié)同配送包裹。目前亞馬遜、德國郵政、谷歌等企業(yè)已經(jīng)嘗試將無人機(jī)應(yīng)用在其業(yè)務(wù)中以降低配送成本縮短配送時間[1-2]。無人機(jī)由于其移動消耗能量少，不受道路網(wǎng)絡(luò)交通擁堵影響，速度快等特性，相比傳統(tǒng)配送車送貨而言可降低交付成本縮短交付時間，但考慮其載荷及續(xù)航能力問題，配送車又顯得更有優(yōu)勢[3-5]。由于二者的互補(bǔ)性，將無人機(jī)整合到現(xiàn)有的“最后一英里”交付模式中可能產(chǎn)生協(xié)同效益。因此引發(fā)了“無人機(jī)-旅行商”問題（Travelling Salesman Problem with Drone，TSP-D）[6]針對無人機(jī)與配送車聯(lián)合路徑的研究。

2015年，Murray等人[6]首次將無人機(jī)引入TSP研究，引入混合線性規(guī)劃公式求解兩類問題模型，提出了一種基于先分群再排路線的啟發(fā)式算法解決了10個節(jié)點(diǎn)的案例，該算法目前只能解決小規(guī)模節(jié)點(diǎn)案例。Agatz 等人[7]引入一個整數(shù)規(guī)劃公式，提出了基于局部搜索和動態(tài)規(guī)劃的先排路線再分群的兩種啟發(fā)式算法，同樣只適用于小于10個節(jié)點(diǎn)的案例。Bouman等人[8]基于動態(tài)規(guī)劃提出了一種求解TSP-D 的精確方法，能更快找到10個節(jié)點(diǎn)問題示例的最優(yōu)解決方案。Sacramento等人[9]研究了TSP-D 的一個新變體，其中有多輛配送車存在，并對Murray等人[6]的混合線性規(guī)劃公式進(jìn)行了調(diào)整，提出了自適應(yīng)大鄰域搜索算法。

綜上，由于TSP-D 是TSP 一個新變體，目前較少文獻(xiàn)研究。為進(jìn)一步改進(jìn)TSP-D 的數(shù)學(xué)模型[10-12]，本文提出了一種新型優(yōu)化迭代算法對該問題進(jìn)行求解，在均勻和聚集兩類示例下隨機(jī)生成10個和11個節(jié)點(diǎn)分析驗證實(shí)驗結(jié)果，并與傳統(tǒng)TSP 解進(jìn)行對比，對算法的有效性進(jìn)行驗證。

1 TSP-D模型構(gòu)建

1.1 TSP-D問題描述

無人機(jī)聯(lián)合配送車的交付方式具體過程為：配送開始前，工作人員為無人機(jī)和配送車裝載好貨物，無人機(jī)?？吭谂渌蛙囓図敗Ｅ渌烷_始，無人機(jī)從車頂發(fā)射，二者分別前往各自客戶節(jié)點(diǎn)配送貨物。由于無人機(jī)的裝載能力與電池壽命有限，其每訪問完一次客戶后需與配送車匯合，裝載下一個客戶的包裹并刷新電池。因此，客戶節(jié)點(diǎn)可分為三類：卡車節(jié)點(diǎn)、無人機(jī)節(jié)點(diǎn)以及匯合節(jié)點(diǎn)。

無人機(jī)和配送車是否能同時到達(dá)匯合節(jié)點(diǎn)取決于二者在各自路徑上的耗時。三種情況，當(dāng)配送車先到達(dá)匯合節(jié)點(diǎn)時，需在節(jié)點(diǎn)處等待無人機(jī)到來，等待時間為Wt；若無人機(jī)先到達(dá)匯合節(jié)點(diǎn)或者二者同時到達(dá)匯合節(jié)點(diǎn)，配送車只需沿其路徑繼續(xù)行駛無需等待即Wt為0。因此，整個交付時間可看作配送車配送時長Ttruck與配送車等待時長Wt之和。TSP-D問題要求無人機(jī)和配送車的路線必須在倉庫開始和結(jié)束，目標(biāo)是最小化交付時間。

如圖1 說明了這種交付方式。直線代表配送車路徑，虛線代表無人機(jī)路徑。0 點(diǎn)代表倉庫，節(jié)點(diǎn)2、5、3、4、1、10是配送車服務(wù)的客戶節(jié)點(diǎn)，無人機(jī)服務(wù)的客戶節(jié)點(diǎn)為7、8、9、6，節(jié)點(diǎn)5、4、10為二者匯合節(jié)點(diǎn)。其中配送車路徑為0-2-5-3-4-1-10-0，無人機(jī)路徑為0-7-5-8-4-9-10-6-0。如圖2為傳統(tǒng)單配送車交付方式的模型。

圖1 TSP-D問題模型Fig.1 TSP-D problem model

圖2 TSP問題模型Fig.2 TSP problem model

1.2 約束條件

（1）無人機(jī)載重有限。無人機(jī)每訪問完一個客戶節(jié)點(diǎn)后需返回配送車拿取下一件包裹并刷新電池，接著再訪問下一個客戶節(jié)點(diǎn)，即無人機(jī)配送完貨物的下一個節(jié)點(diǎn)總是配送車節(jié)點(diǎn)。如圖3所示，0點(diǎn)為倉庫，圓形代表配送車，方形代表無人機(jī)。

圖3 無人機(jī)路徑示例Fig.3 Drone path example

（2）無人機(jī)續(xù)航時間有限。無人機(jī)每次從起飛到降落至配送車這段路徑耗時需在無人機(jī)的電池續(xù)航時間之內(nèi)。如圖3 中的各段無人機(jī)路徑0-6-2，2-7-5……每段路徑耗時都需小于等于無人機(jī)的最大續(xù)航時間才能完成配送任務(wù)。

（3）無人機(jī)與配送車匯合節(jié)點(diǎn)不能是上一個匯合節(jié)點(diǎn)之前的節(jié)點(diǎn)，即不能返回至配送車路過的節(jié)點(diǎn)。如圖4 所示配送車路徑為0-1-2-3-0，無人機(jī)路徑的第一段為0-4-2，第二段為2-5-1，由于配送車與無人機(jī)是同步運(yùn)動，二者在時間上也是同步的，顯然無人機(jī)的第二段路徑不成立，配送車此時已路過節(jié)點(diǎn)1，無人機(jī)若返回節(jié)點(diǎn)1則無法降落在配送車上。

圖4 不滿足條件的無人機(jī)路徑示例Fig.4 Example of drone path that does not meet criteria

（4）無人機(jī)不允許重復(fù)訪問客戶。

（5）配送車不允許重復(fù)發(fā)射或接收無人機(jī)。

（6）總配送時間最短。

1.3 符號說明

（1）集合：

C表示客戶節(jié)點(diǎn)集合；

R表示配送車路徑集合（在第一階段確定）；

D表示分配給無人機(jī)的客戶節(jié)點(diǎn)集合（在第一階段確定）；

P表示無人機(jī)路徑集合（在第二階段確定）。

（2）參數(shù)：

N表示客戶節(jié)點(diǎn)個數(shù)；

T表示無人機(jī)與配送車聯(lián)合配送總耗時；

EN表示無人機(jī)續(xù)航時間；

a、b、c表示節(jié)點(diǎn)位置；

r表示路徑；

pk表示配送車路徑中的節(jié)點(diǎn)k；

ti表示配送車到達(dá)i處的耗時；

di表示無人機(jī)到達(dá)i處的耗時。

（3）決策變量：

nar表示路徑r中a節(jié)點(diǎn)若為無人機(jī)服務(wù)節(jié)點(diǎn)即為1，否則為0；

sar表示若路徑r從節(jié)點(diǎn)a出發(fā)即為1，否則為0；

ecr表示若路徑r在節(jié)點(diǎn)c結(jié)束即為1，否則為0；

xr表示若路徑r為無人機(jī)路徑即為1，否則為0；

yab若配送車從a行駛到b(a≠b)即為1，否則為0；

xabc表示若無人機(jī)路徑從a發(fā)射，在c降落，途中訪問b即為1，否則為0。

1.4 模型建立

針對TSP-D 問題建立以下數(shù)學(xué)模型：公式（1）為目標(biāo)函數(shù)，表示無人機(jī)與配送車聯(lián)合配送的總耗時最短；式（2）計算聯(lián)合配送的總耗時，分為兩部分：第一部分為配送車配送時長；第二部分為配送車等待時長。小括號中表達(dá)式表示配送車在發(fā)射和回收節(jié)點(diǎn)間的耗時，大括號中表達(dá)式表示發(fā)射和回收節(jié)點(diǎn)a、c間配送車與無人機(jī)的時間差，即配送車等待時長；式（3）表示無人機(jī)不能在返回配送車之前被重復(fù)發(fā)射，如無人機(jī)在a點(diǎn)被發(fā)射，在c點(diǎn)被回收，假設(shè)k點(diǎn)為a與c之間的點(diǎn)，那么無人機(jī)不允許在k點(diǎn)被再次發(fā)射或回收；式（4）表示若無人機(jī)從a點(diǎn)發(fā)射在c點(diǎn)被回收，則配送車在訪問c之前需先訪問a；式（5）表示若無人機(jī)從a點(diǎn)發(fā)射在c點(diǎn)被回收，則a與c為分配給配送車的客戶點(diǎn)；式（6）表示若無人機(jī)從a發(fā)射在c被回收，則無人機(jī)耗時需在無人機(jī)續(xù)航時間內(nèi)；式（7）表示每個客戶節(jié)點(diǎn)只能被無人機(jī)訪問一次；式（8）表示一個節(jié)點(diǎn)作為發(fā)射節(jié)點(diǎn)至多一次；式（9）表示一個節(jié)點(diǎn)作為回收節(jié)點(diǎn)至多一次；式（10）至式（15）表示二進(jìn)制決策變量。

2 TSP-D求解方法

2.1 求解算法

TSP-D屬于NP難題（Non-deterministic Polynomialhard，NP-hard），由于其非確定性，傳統(tǒng)方法難以快速找到最優(yōu)解，為縮小搜索范圍，在合理時間內(nèi)得到最優(yōu)解[13-15]，故采用一種新型優(yōu)化迭代方法求解問題。

本文將問題分為兩階段：第一階段確定配送車路徑；第二階段確定無人機(jī)路徑。由于配送車無需考慮重復(fù)裝載貨物等問題，首先確定配送車路徑。不在配送車路徑上的節(jié)點(diǎn)即為無人機(jī)配送節(jié)點(diǎn)，因此第一階段解決了配送車路徑和客戶節(jié)點(diǎn)分配兩個問題。第二階段，固定第一階段配送車路徑和無人機(jī)配送節(jié)點(diǎn)，從配送車節(jié)點(diǎn)中確定滿足約束的匯合節(jié)點(diǎn)，得到無人機(jī)配送路線。方法如圖5、圖6所示，具體如下：

圖5 第一階段固定配送車路線Fig.5 Fixed truck route in the first stage

圖6 第二階段確定無人機(jī)路線Fig.6 Determined drone route in the second stage

第一階段，生成所有可能的配送車路徑。為確定配送車路徑，首先需要確定可能的配送車節(jié)點(diǎn)。從所有節(jié)點(diǎn)中去除可能與配送車組合的最多個無人機(jī)節(jié)點(diǎn)，剩余節(jié)點(diǎn)即為最少的配送車節(jié)點(diǎn)。從最少的配送車節(jié)點(diǎn)開始直至遍歷C-1 個節(jié)點(diǎn)（只有一個無人機(jī)的情況），生成所有可能的配送車路徑。由于無人機(jī)每配送完一次貨物需返回配送車拿取下件貨物。例如有10個配送節(jié)點(diǎn)，為滿足以上條件，無人機(jī)節(jié)點(diǎn)最多5 個最少1 個，則配送車對應(yīng)路徑即為5個節(jié)點(diǎn)至9個節(jié)點(diǎn)的任意排列組合；若有11 個配送節(jié)點(diǎn)，無人機(jī)節(jié)點(diǎn)最多6 個最少1 個才能保證每個無人機(jī)節(jié)點(diǎn)都有“返回點(diǎn)”，則配送車路徑為5到10個節(jié)點(diǎn)的隨機(jī)組合。在開始求解前，首先計算得到N個節(jié)點(diǎn)由“傳統(tǒng)單貨車配送”的最優(yōu)路徑解即TSP 解，為得到可行的配送車路徑并縮小搜索范圍，只保留路徑距離小于TSP解的配送車路徑，將配送車路徑按其耗時升序存儲在集合R中。算法流程圖如圖7所示。

圖7 第一階段算法流程圖Fig.7 Flow chart of first-stage algorithm

第二階段，固定第一階段確定的配送車路徑及無人機(jī)節(jié)點(diǎn)確定無人機(jī)路線。算法開始，將TSP解作為全局上界ub，對于固定的配送車路徑將相應(yīng)剩余的無人機(jī)節(jié)點(diǎn)在滿足“每個無人機(jī)都有返回節(jié)點(diǎn)”的條件下，隨機(jī)插空進(jìn)配送車路徑中，得到所有可能的無人機(jī)路徑；然后考慮無人機(jī)續(xù)航時間的約束，計算各段無人機(jī)起飛至降落到配送車的路徑耗時，保留滿足無人機(jī)續(xù)航時間的無人機(jī)路徑；接著結(jié)合對應(yīng)的配送車路徑，針對其相應(yīng)的無人機(jī)路徑計算聯(lián)合配送情況下配送車的等待時間Wt，得到聯(lián)合配送時間T為配送車耗時Truck與Wt之和；最后，若T小于TSP 解則保留該條聯(lián)合配送路徑，并將T作為新的全局上界否則不保留，以此方式更新改進(jìn)ub。由于T為Ttruck與Wt的總和，而在第一階段配送車路徑按Ttruck升序排列，因此當(dāng)?shù)恋趐條配送車路徑時，Ttruck若大于最新的ub，則程序直接終止，此時的ub即為最優(yōu)解，其對應(yīng)的路徑即為最優(yōu)聯(lián)合路徑。算法流程圖如圖8所示。

圖8 第二階段算法流程圖Fig.8 Flow chart of second-stage algorithm

2.2 算法示例

為理解算法，現(xiàn)以4 個節(jié)點(diǎn)為例說明算法流程，其中倉庫坐標(biāo)（15.18，15.76），配送節(jié)點(diǎn)坐標(biāo)分別為（6.21，16.42）、（19.76，10.86）、（26.56，15.82），單位為km。算法開始前求得4 節(jié)點(diǎn)的TSP 解為65.16 min。4 節(jié)點(diǎn)情況下，配送車節(jié)點(diǎn)個數(shù)至少為1，至多為2才能使得無人機(jī)有返回點(diǎn)。因此第一階段生成滿足以下兩個條件的配送車路徑：（1）路徑耗時小于TSP解；（2）避免重復(fù)路徑，要求最后一個客戶節(jié)點(diǎn)序號大于第一個客戶節(jié)點(diǎn)序號。例如，路徑0-2-3-0 與0-3-2-0 是相反方向的同一條路徑，為縮小搜索范圍，只保留第一條。最后按照路徑耗時由小到大排列存儲。第一階段確定的配送車路徑如表1所示。

表1 配送車路徑舉例Table 1 Example of delivery vehicle routes

第二階段，固定第一階段確定的配送車路徑，將剩余節(jié)點(diǎn)合理插入配送車路徑中生成滿足約束條件的無人機(jī)路徑，具體迭代過程如表2所示。

表2 無人機(jī)路徑及算法迭代過程舉例Table 2 Example of drone path and algorithm iteration process

表2中第2、3、6次迭代的聯(lián)合路徑總耗時都大于全局上界，因此ub未更新。第4 次迭代的路徑耗時小于ub，因此ub更新為39.76 min。針對表1 中第4 條配送車路徑0-2-3-0，對應(yīng)有3條可能的無人機(jī)路徑如表2中的4 至6 次迭代所示。其中無人機(jī)路徑0-1-3-0 耗時29.35 min，超出了無人機(jī)自身續(xù)航時間25 min，因此不再后續(xù)計算同時刪除該條路徑。由于配送車路徑耗時已經(jīng)按由小到大排序，而第7 次迭代配送車路徑耗時45.52 min 大于最新的全局上界39.76 min，因此算法停止，得到最優(yōu)解為39.76 min。

3 實(shí)驗結(jié)果分析

3.1 示例數(shù)據(jù)

由于TSP-D被引入文獻(xiàn)不久，目前沒有被廣泛接受的基準(zhǔn)數(shù)據(jù)，因此為研究評估本文算法的性能，本節(jié)通過隨機(jī)生成的兩類測試示例來求解TSP-D問題模型，并從運(yùn)行時間、節(jié)點(diǎn)個數(shù)等方面分析了算法在兩種示例上的結(jié)果。最后，從求解質(zhì)量和求解時間等方面對該算法的性能進(jìn)行了評價。

如圖9所示均勻示例，所有客戶節(jié)點(diǎn)坐標(biāo)均勻地分布在0～30 km，形成一個900 km2的配送網(wǎng)絡(luò)。如圖10所示聚集示例，所有客戶節(jié)點(diǎn)集中分布在以倉庫為中心的配送網(wǎng)絡(luò)中。表3和表4分別列出了均勻與聚集示例下，10、11節(jié)點(diǎn)的詳細(xì)實(shí)驗數(shù)據(jù)。

圖9 均勻分布坐標(biāo)Fig.9 Uniformly distributed coordinates

圖10 聚集類坐標(biāo)Fig.10 Aggregation class coordinates

實(shí)驗假設(shè)配送車速度為40 km/h，無人機(jī)速度60 km/h，無人機(jī)電池續(xù)航25 km/h。實(shí)驗使用Matlab語言編程，使用設(shè)備處理器為Intel?Core?i7-8700k CPU，內(nèi)存為32 GB。

3.2 結(jié)果分析

由于無人機(jī)速度大于配送車速度，因此考慮是否分配越多的無人機(jī)節(jié)點(diǎn)越能得到最優(yōu)解。為考察無人機(jī)節(jié)點(diǎn)個數(shù)與求解質(zhì)量間的關(guān)系，觀察表3 和表4 最后一列發(fā)現(xiàn)，最優(yōu)解的無人機(jī)節(jié)點(diǎn)個數(shù)始終小于可能的最多無人機(jī)節(jié)點(diǎn)個數(shù)。如10 節(jié)點(diǎn)情況下，最大可能的無人機(jī)節(jié)點(diǎn)個數(shù)為5（11 節(jié)點(diǎn)時為6），無論在均勻還是聚集示例下，最優(yōu)解的無人機(jī)節(jié)點(diǎn)個數(shù)始終小于5。這表明得到最優(yōu)路徑并不等同分配更多無人機(jī)節(jié)點(diǎn)。由于無人機(jī)與配送車的“同步性”，將大部分客戶分配給無人機(jī)會導(dǎo)致配送車等待時間提升，配送效率反而下降。

表3 均勻示例Table 3 Uniform samples

表4 聚集示例Tablel 4 Aggregation samples

運(yùn)行時間上分析，如表3和表4第5列所示，總體運(yùn)行時間大約控制在半小時以內(nèi)。由于該過程需要枚舉，每增加一個節(jié)點(diǎn)，配送車路徑和無人機(jī)路徑的排列組合呈指數(shù)增加，因此11 節(jié)點(diǎn)相比10 節(jié)點(diǎn)同條件下運(yùn)行時間明顯提高。

表5根據(jù)表3和表4第2、3列數(shù)據(jù)，計算TSP-D交付方式較傳統(tǒng)TSP 節(jié)省效率，結(jié)果顯示10 節(jié)點(diǎn)情況下平均配送時間可節(jié)省30%左右，而11 節(jié)點(diǎn)情況下平均可節(jié)省17%以上，表明將無人機(jī)整合入傳統(tǒng)配送車運(yùn)送模式的有效性以及在聚集示例下算法優(yōu)勢更大。

表5 TSP-D解提高百分比Tablel 5 Percentage increases in TSP-D solution

圖11～13 為聚集示例相同坐標(biāo)下TSP 與TSP-D 兩種配送模式的仿真結(jié)果，圖14～16為均勻示例相同坐標(biāo)下兩種配送模式的仿真結(jié)果。其中圖12、15為TSP-D配送模式求解過程中的兩個隨機(jī)解。圖11得到單配送車耗時110.93 min，圖13相同坐標(biāo)下配送車聯(lián)合無人機(jī)最優(yōu)解耗時83.61 min，配送時間節(jié)省24.63%；圖14單配送車耗時181.16 min，圖16相同坐標(biāo)下配送車聯(lián)合無人機(jī)最優(yōu)解耗時124.37 min，配送時間節(jié)省31.35%。

圖11 單配送車配送路徑（聚集）Fig.11 Single distribution vehicle distribution path（aggregation）

圖12 無人機(jī)-配送車聯(lián)合配送隨機(jī)路徑（聚集）Fig.12 Random path of drone-delivery vehicle joint delivery（aggregation）

圖13 無人機(jī)-配送車聯(lián)合配送最優(yōu)路徑（聚集）Fig.13 Optimal path of drone-delivery vehicle joint delivery（aggregation）

圖14 單配送車配送路徑（均勻）Fig.14 Single distribution vehicle distribution path（uniform）

圖15 無人機(jī)-配送車聯(lián)合配送隨機(jī)路徑（均勻）Fig.15 Random path of drone-delivery vehicle joint delivery（uniform）

圖16 無人機(jī)-配送車聯(lián)合配送最優(yōu)路徑（均勻）Fig.16 Optimal path of drone-delivery vehicle joint delivery（uniform）

4 結(jié)論

無人機(jī)聯(lián)合配送車送貨是一種新型的物流配送模式，TSP-D 主要解決配送車與無人機(jī)的路由問題，目標(biāo)減少“最后一英里”的運(yùn)輸時間，是對傳統(tǒng)TSP問題的改進(jìn)。為優(yōu)化路由決策，本文建立聯(lián)合配送模型并提出一種新型優(yōu)化迭代算法對該問題進(jìn)行求解，結(jié)果表明無人機(jī)聯(lián)合配送車這種新型交貨方式與傳統(tǒng)單配送車交貨相比在時間效率上有所提高，該算法能在合理時間內(nèi)得到11節(jié)點(diǎn)內(nèi)的最優(yōu)解。

就實(shí)際應(yīng)用而言，由于復(fù)雜的道路狀況、法規(guī)等的約束以及問題的NP-hard性質(zhì)，使得解決問題的復(fù)雜度提升，因此在實(shí)際約束條件的限制以及問題規(guī)模的大小等方面有待進(jìn)一步拓展研究。