改進(jìn)混合粒子群算法求解帶時間窗的無人機與車輛協(xié)同路徑調(diào)度問題

摘 要：為提高物流配送效率，考慮時間窗、無人機換電以及無人機多點連續(xù)配送等因素，提出了一種帶時間窗的車輛與無人機協(xié)同配送問題，并設(shè)計一種帶局部搜索的混合粒子群算法進(jìn)行求解。該算法以混合粒子群算法為核心，通過構(gòu)建高效的編解碼策略實現(xiàn)了問題解空間到算法搜索空間的轉(zhuǎn)換。進(jìn)一步，該算法融合單點插入策略、車輛更換策略、無人機更換策略組成局部搜索策略，以此提高算法尋優(yōu)能力。實驗結(jié)果表明：所提模型比純車輛配送的模型效率更高，節(jié)省了31.51%的成本；所提算法優(yōu)于四種對比算法，優(yōu)化率最高達(dá)到82.08%。

關(guān)鍵詞：無人機； 車輛調(diào)度； 粒子群； 時間窗； 車輛路徑問題

中圖分類號：TP301 文獻(xiàn)標(biāo)志碼：A

文章編號：1001-3695（2024）08-013-2336-07

doi：10.19734/j.issn.1001-3695.2023.12.0608

Improved hybrid particle swarm optimization algorithm for vehiclerouting problem with drone and time window

Ye Liwei1， Wu Junhao1，2， Qi Yuanhang1，2， Luo Haoyu1，2， Huang Gewen3， Wang Fujie4

（1.School of Computer Science， University of Electronic Science & Technology of China， Zhongshan Institute， Zhongshan Guangdong 528402， China; 2.School of Automation， Guangdong University of Technology， Guangzhou 510006， China; 3.Information & Network Center， Jiaying University， Meizhou Guangdong 514015， China; 4.Elite Engineers College of Dongguan University of Technology， Dongguan Guangdong 523808， China）

Abstract：In order to improve the delivery efficiency of logistics， this paper proposed a vehicle routing problem with drone and time window considered by the time window， UAV（unmanned aerial vehicle） power change and multi-point continuous delivery of UAV. Then， this paper proposed a hybrid particle swarm optimization algorithm with local search to solve it. Based on the hybrid particle swarm optimization algorithm， the proposed algorithm transformed the problem solution space into the algorithm search space by constructing an efficient encoding and decoding strategy. Further， the proposed algorithm combined a single point insertion strategy， a vehicle replacement strategy and a UAV replacement strategy to form a local search strategy to improve the optimization ability. The experimental results show that the proposed model is more efficient than the model of pure vehicle distribution and saves 31.51% of the cost. The proposed algorithm is better than the four comparison algorithms， and its highest optimization rate is 82.08%.

Key words：drone; vehicle scheduling; particle swarm; time window; vehicle routing problem

0 引言

隨著無人機技術(shù)的不斷發(fā)展，無人機開始被應(yīng)用于物流配送領(lǐng)域。亞馬遜、谷歌、DHL、順豐等企業(yè)相繼開展了無人機配送項目的研究及實際應(yīng)用［1］。但是無人機存在載重能力小、續(xù)航時間短、配送距離短等缺點［2］。如何協(xié)調(diào)無人機與車輛以實現(xiàn)更加高效的配送就成為了當(dāng)前物流研究的一個熱點問題。因此，無人機與車輛聯(lián)合的車輛路徑問題（vehicle routing problem with drone，VRPD）［3，4］應(yīng)運而生。Murray等人［5］以配送時間最短為目標(biāo)，建立了兩種無人機車輛配送模型。進(jìn)一步，Ha等人［6］提出車輛無人機協(xié)同模型可以降低運營成本。王新等人［7］設(shè)計了一種自適應(yīng)大規(guī)模鄰域搜索方法，證明了卡車與無人機聯(lián)合配送的優(yōu)越性。李妍峰等人［8］則建立了混合整數(shù)規(guī)劃模型，大大降低了運輸成本和人力成本。在上述研究中，無人機具有固定的電量，在電量用盡后無法繼續(xù)配送。因此，部分學(xué)者提出無人機可以返回車上進(jìn)行換電服務(wù)的方案，以延長無人機的飛行時間和配送距離。

考慮無人機可換電因素，Raeesi等人［9］建立了時空同步的混合整數(shù)線性規(guī)劃模型，為電動車更換能量耗盡的電池。進(jìn)一步，顏瑞等人［10］構(gòu)建了卡車與無人機聯(lián)合的車輛路徑問題的混合整數(shù)線性規(guī)劃模型。范厚明等人［11］則進(jìn)一步考慮載重、續(xù)航、區(qū)域路網(wǎng)交通信息等約束建立模型。而在現(xiàn)實配送的過程中，客戶會要求貨物在指定的時間區(qū)間內(nèi)到達(dá)（時間窗），這種情況下的無人機換電時間變得不可忽略。因此，Huang等人［12］提出無人機在返回車輛后需要進(jìn)行一段時間的恢復(fù)才可重新出發(fā)。而Di Puglia等人［13］考慮到帶時間窗車輛路徑問題在物流配送中的廣泛應(yīng)用，建立了以最小化運輸成本為目標(biāo)的數(shù)學(xué)模型。值得注意的是，文獻(xiàn)［12，13］的無人機每次飛行只允許配送1個客戶，這并不符合現(xiàn)實的無人機現(xiàn)狀。因為隨著無人機技術(shù)的發(fā)展，現(xiàn)有的無人機容量得到了一定提高，并允許在1次飛行中為多個客戶服務(wù)。

因此，本文考慮了時間窗、無人機換電以及無人機多點連續(xù)配送等因素，提出了帶時間窗的無人機與車輛協(xié)同調(diào)度問題（vehicle routing problem with drone and time window，VRPDTW）。不難發(fā)現(xiàn)，VRPDTW是一個組合優(yōu)化問題。針對這類問題，智能算法能在合理的時間內(nèi)得到一個可行解［14］。因此，本文以混合粒子群算法［15］為核心，構(gòu)建了快速有效的編解碼策略以及局部搜索策略，提出了帶局部搜索的混合粒子群算法（hybrid particle swarm optimization with local search strategy，HPSO-LSS）求解VRPDTW。最后，通過實驗證明了本文模型及算法的有效性。

1 問題建模

1.1 問題描述

VRPDTW的配送流程為：面對一群有不同貨物需求的客戶，配送中心向客戶提供貨物，由一個車隊（每臺車攜帶單臺無人機）負(fù)責(zé)在客戶規(guī)定的時間窗內(nèi)將貨物送達(dá)并進(jìn)行卸貨服務(wù)。該過程包括：

a）車輛為無人機提供存儲貨物以及自動換電服務(wù)（無須人工操作）。無人機在完成某些客戶的配送任務(wù)后，回到車輛上自動換電（無人機恢復(fù)為其最大續(xù)航），可以再從該車輛上起飛并繼續(xù)為其他客戶送貨。

b）車輛可在某個節(jié)點（配送中心或者客戶點）發(fā)射或回收無人機。當(dāng)車輛的無人機起飛后，執(zhí)行完配送任務(wù)將降落到原車輛。

c）車輛和無人機可以同步進(jìn)行配送服務(wù)，即無人機外出執(zhí)行配送任務(wù)時，車輛也可以為其他客戶配送貨物。

d）如果車輛與無人機將在某一客戶點匯合，則默認(rèn)該客戶點使用車輛進(jìn)行配送，卸貨服務(wù)時間將從車輛抵達(dá)該客戶點的時刻開始計算，而無人機的自動換電服務(wù)時間則從車輛和無人機匯合的時刻開始計算。在這個過程中，無人機的換電服務(wù)以及車輛的卸貨服務(wù)可以同時進(jìn)行。此外，車輛需要完成上述兩個服務(wù)后才能繼續(xù)前往下一個客戶點，而無人機完成換電服務(wù)便可以重新起飛。

VRPDTW目標(biāo)在于確定車輛以及無人機的混合配送路徑，在滿足客戶需求的前提下，實現(xiàn)配送成本最小。其配送示意圖如圖1所示，包含了1個配送中心、2臺車輛（每臺車攜帶1臺無人機）以及12個客戶點。在圖1的路徑調(diào)度方案中，車輛1及其無人機負(fù)責(zé)配送客戶1～7，車輛2及其無人機負(fù)責(zé)配送客戶8～12。

1.2 數(shù)學(xué)模型

1.2.1 數(shù)學(xué)變量

1）集合與變量

N為客戶集合，N={1，2，…，P}，P為客戶的數(shù)量；N+為客戶集合與{P+1}的并集，N+=N∪{P+1}；N－為客戶集合與{0}的并集，N－=N∪{0}；NAll為所有節(jié)點集合，NAll=N∪{0，P+1}，其中0表示作為起點的配送中心，P+1表示作為終點的配送中心；K為車的數(shù)量；k為車的索引，k=1，2，…，K；D為無人機的最大發(fā)射次數(shù)；d為無人機發(fā)射次數(shù)索引，d=1，2，…，D；P為客戶的數(shù)量；p為客戶的索引，p=1，2，…，P;w、w′為節(jié)點索引，w=0，1，…，P+1，w′=0，1，…，P+1；θ、σ為節(jié)點索引，θ=0，1，…，P，σ=1，2，…，P+1;Wk為車輛的最大載重量;Vk為車輛的行駛速度；Ck為車輛配送貨物的每公里配送成本；Lk為車輛的最大行駛距離；Wd為無人機的最大載重量；Vd為無人機飛行速度；Cd為無人機配送貨物的每公里配送成本；Ld為無人機的最大續(xù)航；T為無人機自動換電服務(wù)時間；dw，w′為節(jié)點w到節(jié)點w′的距離；Wp為客戶p的貨物需求；［tTWSp，tTWEp］：客戶p的時間窗，tTWSp和tTWEp為時間窗下限和上限；tSp為客戶p的卸貨服務(wù)時間。

2）決策變量

ek，w，w′：若車輛k經(jīng)過路徑（w，w′），ek，w，w′=1，否則ek，w，w′=0；yk，d，w，w′：若車輛k的無人機第d次航程經(jīng)過路徑（w，w′），yk，d，w，w′=1，否則yk，d，w，w′=0;sp，k：若車輛k服務(wù)客戶p，sp，k=1，否則sp，k=0

；up，k，d：若車輛k的無人機第d次航程服務(wù)客戶p，up，k，d=1，否則up，k，d=0；s0k：若車輛k駐留在配送中心，s0k=1，否則s0k=0;uonw，k，d：若節(jié)點w為車輛k的無人機第d次航程起飛節(jié)點，uonw，k，d=1，否則uonw，k，d=0;uoffw，k，d：若節(jié)點w為車輛k的無人機第d次航程降落節(jié)點，uoffw，k，d=1，否則uoffw，k，d=0；tAw，k為車輛k到達(dá)節(jié)點w的時刻，tAw，k≥0；tAw，k，d為車輛k的無人機第d次航程到達(dá)節(jié)點w的時刻，tAw，k，d≥0；tSp，k，d為車輛k的無人機第d次航程在客戶p的服務(wù)實際開始時間，tSp，k，d≥0；tAp為客戶p處車輛/無人機抵達(dá)時間。

1.2.2 目標(biāo)函數(shù)

以最小化配送總成本為目標(biāo)，本文的目標(biāo)函數(shù)如式（1）所示。

min F=∑k∑w∑w′Ckek，w，w′dw，w′+∑k∑d∑w∑w′Cdyk，d，w，w′dw，w′（1）

其中：第一部分表示車輛總配送成本，第二部分表示無人機總配送成本。

1.2.3 約束條件

首先構(gòu)建車輛及無人機的載重、行駛距離以及連續(xù)性配送的約束，如式（2）～（9）所示。

∑pWpsp，k+∑p∑d（Wpup，k，d）≤Wk k=1，2，…，K（2）

∑w∑w′（ek，w，w′dw，w′）≤Lk k=1，2，…，K（3）

∑p（Wpup，k，d）≤Wd k=1，2，…，K；d=1，2，…，D（4）

∑w∑w′（yk，d，w，w′dw，w′）≤Ld k=1，2，…，K；d=1，2，…，D（5）

∑ksp，k+∑k∑dup，k，d=1 p∈N（6）

s0k+∑pek，p，P+1=s0k+∑pek，0，p=1 k=1，2，…，K；p∈N（7）

∑σek，p，σ=∑θek，θ，p=sp，k k=1，2，…，K；p∈N（8）

∑γ∈S∑γ′∈Sek，γ，γ′≤|S|-1 SN；k=1，2，…，K（9）

其中：式（2）～（5）分別為車輛的載重約束、車輛的最大行駛距離約束、無人機的載重約束與無人機的最大航程約束;式（6）表示每個客戶只能由一臺車輛或者一臺無人機提供服務(wù);式（7）（8）表示每個節(jié)點車輛到達(dá)的數(shù)量需要等于車輛離開的數(shù)量；式（9）表示車輛路徑不能出現(xiàn)子回路。進(jìn)一步，定義無人機起飛、飛行、降落的規(guī)則，如式（10）～（19）所示。

uoff0，k，d=0 k=1，2，…，K；d=1，2，…，D（10）

uonP+1，k，d=0 k=1，2，…，D；d=1，2，…，D（11）

uon0，k，d=∑pyk，d，0，p k=1，2，…，K；d=1，2，…，D（12）

uoffP+1，k，d=∑pyk，d，p，P+1 k=1，2，…，K；d=1，2，…，D（13）

yk，d，p，p′+sp，k+sp′，k≤2p∈N；p′∈N；k=1，2，…，K；d=1，2，…，D（14）

|（∑σek，w，σ）∑σyk，d，w，σ-（∑θek，θ，w′）∑θyk，d，θ，w′|-Ψ（1-ek，w，w′+s0k）≤0

k=1，2，…，K；d=1，2，…，D；w∈NAll；w′∈NAll（15）

∑σyk，d，p，σ=∑θyk，d，θ，p=up，k，d

p∈N；k=1，2，…，K；d=1，2，…，D（16）

∑w∑σyk，d，w，σ≥δ∑wuonw，k，d

k=1，2，…，K；d=1，2，…，D；δ∈R+（17）

uonp，k，d+∑θyk，d，θ，p=∑σyk，d，p，σ+uoffp，k，dp∈N；k=1，2，…，K；d=1，2，…，D（18）

∑σyk，d，p，σ（∑d′∈D＼d∑σyk，d′，p，σ）=0

p∈N；k=1，2，…，K；d=1，2，…，D（19）

其中：式（10）表示配送中心0不能作為無人機的降落點；式（11）表示配送中心P+1只能作為無人機的降落點；式（12）表示車輛攜帶的無人機如若在配送中心0起飛，則配送中心0為無人機的起飛點；式（13）表示車輛攜帶的無人機如若在配送中心P+1降落，則配送中心P+1為無人機的降落點；式（14）表示無人機起飛后至少要服務(wù)一個客戶才能降落回原車輛，式（15）表示無人機從車上起飛后，需要降落在車輛服務(wù)的下一個節(jié)點；式（16）～（18）表示無人機到達(dá)與離開的服務(wù)節(jié)點數(shù)量需要相同；式（19）表示同一無人機不能多次經(jīng)過同一節(jié)點。最后，進(jìn)行車輛與無人機和時間相關(guān)的約束，如式（20）～（29）所示。

tA0，k=0 k=1，2，…，K（20）

tAσ，k≥（tSp，k+tSp+dp，σVk－Ψ（1－sp，kek，p，σ））

σ∈N+；p∈N；k=1，2，…，K（21）

tSp，k≥max{tAp，k，tTWSp}－Ψ（1－sp，k）

p∈N；k=1，2，…，K（22）

tAσ，k≥tA0，k+d0，σVk-Ψ（1-ek，0，σuon0，k，d）σ∈N+；k=1，2，…，K；d=1，2，…，D（23）

tAσ，k≥max{tAθ，k，d，tAθ，k}+T+dθ，σVk-Ψ（1-yk，d，θ，σuoffθ，k，d）

θ∈N－；σ∈N+；k=1，2，…，K；d=1，2，…，D（24）

tA0，k，d=0；k=1，2，…，K;d=1，2，…，D（25）

tAσ，k，d≥tSp，k，d+tSp+dp，σVd－Ψ（1－up，k，dyk，d，p，σ）

σ∈N+;p∈N;k=1，2，…，K;d=1，2，…，D（26）

tSp，k，d≥max{tAp，k，d，tTWSp}－Ψ（1－up，k，d）

p∈N;k=1，2，…，K;d=1，2，…，D（27）

tAσ，k，d≥tAθ，k + dθ，σ Vd -Ψ（1-yk，d，θ，σ uonθ，k，d）

θ∈N-;σ∈N+;k=1，2，…，K;d=1，2，…，D（28）

tAσ，k，d≥max（tAθ，k，tAθ，k，d）+T+dθ，σVd-Ψ（1-yk，d，θ，σuoffθ，k，d）

θ∈N-;σ∈N+;k=1，2，…，K;d=1，2，…，D（29）

其中：式（20）表示車輛從配送中心0出發(fā)的時間定義為0時刻；式（21）表示車輛如若在某節(jié)點不需要與無人機匯合，則服務(wù)完該節(jié)點后即可前往下一節(jié)點；式（22）表示車輛開始服務(wù)客戶節(jié)點的時刻不可早于客戶時間窗的下限；式（23）表示車輛可直接從配送中心0前往下一節(jié)點；式（24）表示如若無人機需要與車輛匯合時，車輛需等待無人機換電完成才可前往下一節(jié)點；式（25）表示無人機從配送中心0出發(fā)的時間定義為0時刻；式（26）表示無人機需要在節(jié)點完成服務(wù)任務(wù)后才可前往下一節(jié)點；式（27）表示無人機開始服務(wù)客戶節(jié)點的時間不可早于客戶的時間窗下限；式（28）表示無人機需要等待車輛到達(dá)節(jié)點后才可起飛出發(fā)；式（29）表示無人機與車輛匯合后需進(jìn)行換電服務(wù)才可重新出發(fā)。

2 算法設(shè)計

2.1 混合粒子群算法

在粒子群算法（particle swarm optimization，PSO）中，種群里的每個粒子都有自己的位置、速度以及由優(yōu)化目標(biāo)函數(shù)決定的適應(yīng)度值。但是，傳統(tǒng)的PSO在粒子搜索的過程中，粒子僅僅學(xué)習(xí)了個體歷史最優(yōu)值和全局最優(yōu)值，這容易導(dǎo)致粒子群陷入局部最優(yōu)，難以求解復(fù)雜的組合優(yōu)化問題［16］。為解決這個問題，Liang等人［15］于2021年提出了一種粒子群算法的改進(jìn)算法——混合粒子群算法（hybrid particle swarm optimization，HPSO）。HPSO把種群劃分為精英子群和跟隨子群，以便不同類型的子群能夠分別負(fù)責(zé)不同的搜索任務(wù)。精英子群采用一種交叉學(xué)習(xí)策略，以增強全局搜索能力；而跟隨子群則引入了隨機粒子學(xué)習(xí)策略，以提高算法的局部搜索能力。

2.1.1 子種群的劃分

根據(jù)適應(yīng)度，升序排列所有的粒子，并根據(jù)預(yù)設(shè)的群體比率λγ來確定兩種子群的大小，以此來劃分精英子群和跟隨子群。HPSO可以通過調(diào)節(jié)群體比例λγ，進(jìn)而調(diào)節(jié)算法的全局搜索和局部搜索能力。

2.1.2 學(xué)習(xí)策略

1）基于交叉的綜合學(xué)習(xí)策略

在HPSO算法中，每個粒子可能在不同的維度取到最優(yōu)的值。對于精英種群，本文設(shè)計了一種水平交叉算子，可以使兩個不同維度的不同粒子之間進(jìn)行相同維度的交叉操作，進(jìn)而增強粒子間的相互學(xué)習(xí)，增加粒子的多樣性，提高尋優(yōu)能力。同時，為了控制搜索空間的大小，使得粒子具有跳脫局部最優(yōu)的能力，本文設(shè)計了一種垂直交叉算子，可以對同一個體最優(yōu)位置的不同維度進(jìn)行交叉操作。本文中速度更新公式與傳統(tǒng)PSO速度更新公式保持一致，如式（30）所示。水平交叉算子和垂直交叉算子可以在位置更新公式中體現(xiàn)，由γ1取的概率值和λc的比較，決定采用哪種交叉算子來更新粒子的位置，如式（31）所示。

vi，a=ωvi，a+c1r1（pbesti，a-xi，a）+c2r2（gbesta-xi，a）（30）

xi，a=r1pbesti，a+（1-r1）pbestj，a+c（pbesti，a-pbestj，a） γ1≤λcr2pbesti，a+（1-r2）pbesti，bγ1>λc（31）

其中：粒子i的速度vt={vt，1，vt，2，…，vt，a，…，vt，A}；粒子i的位置xi={xi，1，xi，2，…，xi，a，…，xi，A}；xi和vi的最大維度為A；a、b是xi和vi的任意維度，且a≠b；gbest代表全局最優(yōu)個體位置；pbesti代表粒子xi的個體最優(yōu)位置；pbestj代表粒子xj的個體最優(yōu)位置；ω是慣性權(quán)值；c1和c2是加速度系數(shù)；c是在［-1，1］均勻分布的隨機值；r1、r2、γ1、λc是在［0，1］的隨機值。

2）隨機粒子學(xué)習(xí)策略

對于跟隨子群，將粒子根據(jù)適應(yīng)度進(jìn)行排序，適應(yīng)度較差的粒子排在后面，為引導(dǎo)適應(yīng)度較差的粒子尋找更優(yōu)空間，本文設(shè)計了一種隨機粒子學(xué)習(xí)策略。假設(shè)排序后一粒子xi，前面有i-1個粒子{x1，…，xi-1}，隨機在這i-1個粒子中抽取一個粒子xe作為適應(yīng)度較差粒子xi的學(xué)習(xí)榜樣，進(jìn)行如式（32）所示的交叉運算，生成一個速度矢量?；谠撍俣仁噶?，讓適應(yīng)度較差的粒子學(xué)習(xí)適應(yīng)度較好粒子的經(jīng)驗，去探索適應(yīng)度較好粒子的空間，可以得到一個新的粒子位置，如式（33）所示。

vi，a=ωvi，a+c1r1（（r2pbesti，a+（1-r2）pbeste，a）-xi，a）（32）

xi，a=pbesti，a+c2r1（gbesta－pbesti，a） γ2≤λm

xi，a+vi，aγ2>λm（33）

其中：粒子i的速度vt={vt，1，vt，2，…，vt，a，…，vt，A}；粒子i的位置xi={xi，1，xi，2，…，xi，a，…，xi，A}；xi和vi的最大維度為A；a是xi和vi的任意一維度；e是從{x1，…，xi-1}中隨機選取的粒子xe的序號；gbest代表全局最優(yōu)個體位置；pbestj代表粒子xj的個體最優(yōu)位置；pbeste是較好粒子xe的個體最優(yōu)位置；ω是慣性權(quán)值；c1和c2為加速度系數(shù)；r1、r2、γ2、λm是在［0，1］的隨機值。

2.2 編解碼策略

2.2.1 編碼策略

本文設(shè)計的粒子位置包含客戶序列信息、運載工具序列信息、路徑總數(shù)和路徑分段序列信息四部分信息。假設(shè)客戶數(shù)為P，車的數(shù)量為K，每臺車上的無人機數(shù)量為1，則編碼粒子i的位置xi=（xi，1，xi，2，…，xi，2P+K+1）如圖2所示。

在圖2中，xi，1， … ，xi，P∈［-100，100］代表客戶索引； xi，P+1， … ，xi，2P∈{0，1}代表客戶對應(yīng)的配送工具（車輛或無人機）； xi，2P+1∈{1，2，…，K}代表實際路徑總數(shù)（實際使用的車輛數(shù)）； xi，2P+2， … ，xi，2P+K+1∈［0，100］代表路徑分段信息。本文將使用一個具體的例子進(jìn)行闡述，更好地描述本文的編碼策略。假設(shè)客戶的數(shù)量為10、車的數(shù)量為4，則粒子的編碼如圖3所示。圖3中粒子i位置的維度為10+10+1+4=25，進(jìn)而根據(jù)編碼的取值規(guī)則隨機便可得到粒子每個維度上的值。接下來將詳細(xì)介紹粒子的解碼方法。

2.2.2 解碼策略

為對應(yīng)上述的編碼策略，本文所設(shè)計的解碼策略包括以下四個步驟：

a）獲取客戶的配送順序。xi，1，…，xi，P的值分別對應(yīng)著客戶1～P，將xi，1，…，xi，P的值進(jìn)行升序排列，便能得到對應(yīng)的客戶配送序列xSi。本文將繼續(xù)沿用圖3的例子進(jìn)行闡述，將圖3中粒子前10維度的數(shù)據(jù)進(jìn)行升序排列后，得到客戶配送序列xSi=（8，1，5，10，3，6，4，2，7，9）。

b）獲取客戶的配送方法。xi，P+1，…，xi，2P的值分別對應(yīng)著客戶1～P的配送工具選擇，xi，P+p=0表示客戶p由汽車負(fù)責(zé)配送，xi，P+p=1表示客戶p由無人機進(jìn)行配送。將客戶的對應(yīng)配送方法與上述所獲得的配送序列相匹配結(jié)合，便可以得到客戶的配送序列與其配送工具的匹配信息。繼續(xù)沿用圖3中的例子，可以得到客戶1～P的配送工具序號為（0，1，0，1，1，1，0，0，1，1），接下來將其與配送序列進(jìn)行匹配結(jié)合，可以得到匹配的信息{（8，0），（1，0），（5，1），（10，1），（3，0），（6，1），（4，1），（2，1），（7，0），（9，1）}，為直觀表明匹配方法，匹配過程如圖4所示。從圖4可以看出，客戶1、3、7和8采用車輛配送；客戶2、4、5、6、9、10采用無人機進(jìn)行配送。

c）獲取實際的路徑總數(shù)以及每段路徑的客戶數(shù)。xi，2P+1的值表示實際的路徑總數(shù)，同時令λ=xi，2P+1。當(dāng)λ=1，則表示只有1條路徑，所有客戶均在這條路徑上配送。當(dāng)λ>1時，需要使用（xi，2P+2，…，xi，2P+M+1）中前λ維計算每段路徑的客戶數(shù)。令路徑1，2，…，λ′，…，λ分別對應(yīng)xi，2P+2，xi，2P+3，…，xi，2P+1+λ′，…，xi，2P+1+λ。對于路徑λ′（λ′=1，2，…，λ），其客戶數(shù)為

cosλ′= xi，2P+1+λ′∑λτ=1xi，2P+1+τ（P-λ）」+1 λ′≠λP-∑λ-1τ=1cOSτλ′≡λ（34）

其中： 」是向下取整數(shù)。

由圖3可知，xi，21=3表示路徑總數(shù)為3，cOS1、cOS2和cOS3計算過程如下所示。

cos1= 27.827.8+46.5+23.3×（10-3）」+1=2

cos2= 46.527.8+46.5+23.3×（10－3）」+1=4

cos3=10－2－4=4

在本例子中，由于路徑總數(shù)為3，因而作為路徑4值的xi，25不需要被使用。

d）獲取分段路徑信息及具體的配送路徑。根據(jù)cOS1，cOS2，…，cOSλ的值，從左到右劃分客戶配送序列，得到臨時分段信息。由于所有路徑均是從配送中心出發(fā)并回到配送中心，進(jìn)一步在每段臨時分段信息的始末加上配送中心“0”，再結(jié)合每個節(jié)點對應(yīng)的配送方法（默認(rèn)配送中心的配送工具為車輛），從而得到最終的分段路徑信息，過程如圖5所示。

在圖5中，分段路徑1信息為{（0，0），（8，0），（1，0），（0，0）}，分段路徑2信息為{（0，0），（5，1），（10，1），（3，0），（6，1），（0，0）}，分段路徑3信息為{（0，0），（4，1），（2，1），（7，0），（9，1），（0，0）}。根據(jù)這些分段信息，易得各個車輛包含其無人機的具體路徑。例如，車輛1的配送路徑為{（0，0），（8，0），（1，0），（0，0）}，其沒有應(yīng)用無人機；車輛2的配送路徑為{（0，0），（3，0），（0，0）}，其無人機的配送路徑包含2條，{（0，0），（5，1），（10，1），（3，0）}、{（3，0），（6，1），（0，0）}；車輛3的配送路徑為{（0，0），（7，0），（0，0）}，其無人機的配送路徑包含2條，{（0，0），（4，1），（2，1），（7，0）}，{（7，0），（9，1），（0，0）}。獲得具體路徑后，可根據(jù)目標(biāo)式（1）計算各個粒子的適應(yīng)度。

2.3 局部搜索策略

局部搜索策略可有效提高算法的搜索能力［17］。因此，本文設(shè)計了三種不同的局部搜索策略加強算法的求解能力，分別為隨機取點插入策略、隨機取點更換運輸工具策略和遍歷更換無人機策略。

a）單點插入策略：隨機取一個客戶點及其對應(yīng)的運輸工具插入到另一個客戶點前面，從而獲得一個新的配送方案。如果配送方案變得更優(yōu)，則更新配送方案；否則，保持原方案。其中，插入情況包括不同路徑間的單點插入、同段路徑的單點插入。

b）車輛更換策略：隨機取一個客戶，將它的運輸工具更換為車輛，從而獲得一個新的配送方案。如果配送方案變得更優(yōu)，則更新配送方案；否則，保持原方案。

c）無人機更換策略：隨機取一組路徑，從第一個客戶點開始，從左到右遍0+rKdtGu1+uX+d9J+IPM5ccXV/H6bVWUpaXnU/02E7Y=歷地依次更換客戶點的運輸工具為無人機（若原為使用無人機配送，則無須更換）。如果配送方案變得更優(yōu)，則更新配送方案；否則，繼續(xù)更新下一客戶點的運輸工具，直到該組路徑的客戶點遍歷完畢。

2.4 算法流程

求解VRPDTW問題的HPSO-LSSD算法流程如圖6所示。

3 實驗與分析

3.1 實驗算例及實驗環(huán)境

實驗數(shù)據(jù)集采用solomon標(biāo)準(zhǔn)數(shù)據(jù)集中的C109作為基礎(chǔ)算例。為了更符合現(xiàn)實場景數(shù)據(jù)，本文按比例適當(dāng)調(diào)整了C109的部分?jǐn)?shù)據(jù)，并引入車輛以及無人機實際的工作參數(shù)，如表1所示。此外，實驗的仿真環(huán)境是Windows 10，Intel Core i5-11400 CPU@2.60 GHz，16 GB RAM。同時，為了保證算例分析的公平性，每個算法的種群大小都設(shè)置為100，迭代次數(shù)為200次。每個算法獨立運行10次得到實驗結(jié)果。其中HPSO-LSSD和PSO-LSSD算法的相關(guān)參數(shù)設(shè)置如下：w=0.5，c1=2.0，c2=2.0，λγ=0.5，λc=0.5，λm=0.5。

3.2 算法性能實驗

為了驗證模型的有效性，本實驗假設(shè)VRPDTW的所有配送只允許由車輛來完成，其他條件不變，得到一個新模型VRPDTW_N。接著，使用HPSO-LSSD分別求解VRPDTW和VRPDTW_N，實驗結(jié)果如表2所示。

通過表2可以看出，在平均值方面，雖然VRPDTW比VRPDTW_N多使用了9.65元的無人機成本，但VRPDTW的總成本比VRPDTW_N的減少了51.58元，節(jié)省了31.51%的成本。由此證明，無人機的加入可以有效提高車輛配送的效率以及降低配送成本。VRPDTW在第6次實驗求解得到最優(yōu)方案，其配送方案如圖7所示，而具體的路徑信息如表3所示。表3的“具體配送信息”中，每一個節(jié)點以Z（Y，X）來表示，Z表示節(jié)點序號，Y表示配送工具（0為車輛，1為無人機），X表示車輛/無人機的到達(dá)時間（單位：s），節(jié)點0的到達(dá)時間以0計算。

從表3可以看出，該方案每段路徑的車輛/無人機都滿足續(xù)航約束、載重約束和時間約束。接下來以路徑4為例具體分析。路徑4的車輛/無人機的行駛/飛行里程分別為3.81 km、0.36 km、1.40 km、0.77 km和0.62 km，符合式（3）和（5）的約束。路徑4的車輛實際載重為60+5+10+15+15=105 kg，小于車輛最大載重量。同時，無人機的單次航程載重分別為5 kg、10 kg、15 kg和15 kg，也小于無人機的最大載重量。其中，路徑4的無人機在1次飛行中連續(xù)服務(wù)2個客戶，證明模型成功應(yīng)用到了方案中。此外，路徑4中每個客戶有且只有一臺車輛/無人機提供服務(wù)，也滿足配送過程中的時間窗約束。由此可見，HPSO-LSSD能有效地求解VRPDTW。

3.3 算法對比實驗

為了驗證本文算法的尋優(yōu)性能，本實驗采用四個算法來求解VRPDTW配送模型，所得的實驗結(jié)果與HPSO-LSSD作對比。其中，第1個算法是傳統(tǒng)PSO；第2個對比算法是HPSO；第3個對比算法是PSO-LSSD，即傳統(tǒng)PSO加上本文所提出的LSSD；第四個算法是變鄰域搜索算法（variable neighborhood search，VNS）［18］。5個算法獲得最優(yōu)解的迭代過程如圖8所示。

由圖8可以看出，HPSO-LSSD和PSO-LSSD、VNS在迭代20次前便收斂到最優(yōu)值，其收斂速度和效果明顯優(yōu)于HSPO和PSO。10次獨立實驗的結(jié)果如圖9所示。

在圖9中，HPSO-LSSD和PSO-LSSD、VNS所獲得的方案成本在75～175元，而HPSO和PSO所獲取的方案成本在375～500元。由此可見，HPSO-LSSD和PSO-LSSD、VNS的求解穩(wěn)定性明顯高于HPSO和PSO。具體的實驗數(shù)據(jù)如表FKlkjDq++VuZh+GEca0Efghs3k6NLBlRWiGiLed6Iqk=4所示。

從表4可以看出，無論是最優(yōu)解、平均解或者最差解，HPSO-LSSD均優(yōu)于其他四種算法。HPSO-LSSD所得最優(yōu)解的總成本分別比PSO-LSSD、HPSO、PSO、VNS節(jié)約了27.54元、304.91元、385.80、27.08元，優(yōu)化率分別達(dá)到25.30%、78.94%、82.59%、24.98%。此外，HPSO-LSSD的平均解比HPSO的平均解節(jié)約了316.76元，優(yōu)化率為73.72%；PSO-LSSD的平均解比PSO的平均解節(jié)約了395.77元，優(yōu)化率為77.80%。這證明了本文所提的局部搜索算法能有效提高算法的局部尋優(yōu)能力。

4 結(jié)束語

針對車輛與無人機協(xié)同配送問題，本文進(jìn)一步考慮時間窗、無人機換電以及無人機多點連續(xù)配送等因素，提出了VRPDTW，并設(shè)計了HPSO-LSSD進(jìn)行求解。為驗證VRPDTW以及HPSO-LSS的有效性，本文對solomon標(biāo)準(zhǔn)數(shù)據(jù)集中的C109算例進(jìn)行改造，設(shè)計了實驗算例。實驗結(jié)果表明：本文提出的VRPDTW的配送成本遠(yuǎn)遠(yuǎn)低于只采用車輛進(jìn)行配送的VRPDTW_N，節(jié)省了31.51%的成本；本文算法優(yōu)化策略高效可靠，HPSO-LSSD尋優(yōu)能力遠(yuǎn)優(yōu)于PSO-LSSD、HPSO、PSO、VNS，優(yōu)化率分別為25.71%、70.61%、82.08%、24.98%。

參考文獻(xiàn)：

［1］劉正元， 王清華. 無人機和車輛協(xié)同配送映射模式綜述與展望［J］. 系統(tǒng)工程與電子技術(shù)， 2023， 45（3）： 785-796. （Liu Zhengyuan， Wang Qinghua. Review and prospect of UAV and vehicle collaborative distribution mapping mode［J］. System Engineering and Electronic Technology， 2023， 45（3）： 785-796.）

［2］任璇， 黃輝， 于少偉， 等. 車輛與無人機組合配送研究綜述［J］. 控制與決策， 2021， 36（10）： 2313-2327. （Ren Xuan， Huang Hui， Yu Shaowei， et al. Summary of research on vehicle and UAV combined distribution［J］. Control and Decision-Making， 2021， 36（10）： 2313-2327.）

［3］雷定猷， 宋文杰， 張英貴. 平衡裝載約束下的車輛路徑問題研究［J］. 計算機應(yīng)用研究， 2020， 37（6）： 1622-1625，1641. （Lei Dingyou， Song Wenjie， Zhang Yinggui. Research on vehicle routing problem under balanced loading constraints［J］. Application Research of Computers， 2020， 37（6）： 1622-1625，1641.）

［4］Wang Xingyin， Poikonen S， Golden B. The vehicle routing problem with drones： several worst-case results［J］. Optimization Letters， 2017， 11： 679-697.

［5］Murray C C， Chu A G. The flying sidekick traveling salesman problem： optimization of drone-assisted parcel delivery［J］. Transportation Research Part C： Emerging Technologies， 2015， 54： 86-109.

［6］Ha Q M， Deville Y， Pham Q D， et al. On the min-cost traveling salesman problem with drone［J］. Transportation Research Part C： Emerging Technologies， 2018， 86： 597-621.

［7］王新， 王征， 徐偉. 面向多個無人機站點的車輛與無人機聯(lián)合配送路徑問題研究［J］. 運籌與管理， 2021， 30（5）： 31-37. （Wang Xin， Wang Zheng， Xu Wei. Research on vehicle and UAV joint distribution routing problem for multiple UAV sites［J］. Operation and Management， 2021， 30（5）： 31-37.）

［8］李妍峰， 李佳， 向婷. 需求可拆分的無人機與卡車協(xié)同路徑優(yōu)化問題［J］. 工業(yè)工程， 2022， 25（1）： 54-63，143. （Li Yanfeng， Li Jia， Xiang Ting. Cooperative path optimization problem of UAV and truck with split demand［J］. Industrial Engineering， 2022， 25（1）： 54-63，143.）

［9］Raeesi R， Zografos K G. The electric vehicle routing problem with time windows and synchronised mobile battery swapping［J］. Transportation Research Part B： Methodological， 2020， 140： 101-129.

［10］顏瑞， 陳立雙， 朱曉寧， 等. 考慮區(qū)域限制的卡車搭載無人機車輛路徑問題研究［J］. 中國管理科學(xué)， 2022， 30（5）： 144-155. （Yan Rui， Chen Lishuang， Zhu Xiaoning， et al. Research on vehicle routing problem of truck with UAV considering regional constraints［J］. China Management Science， 2022， 30（5）： 144-155.）

［11］范厚明， 張躍光， 田攀俊. 時變路網(wǎng)下多中心電動車-無人機協(xié)同配送路徑優(yōu)化［J］. 管理工程學(xué)報， 2023， 37（2）： 131-142. （Fan Houming， Zhang Yueguang， Tian Panjun. Multi-center electric vehicle-UAV cooperative distribution path optimization under time-varying road network［J］. Chinese Journal of Industrial Enginee-ring and Engineering Management， 2023， 37（2）： 131-142.）

［12］Huang S H， Huang Y H， Blazquez C A， et al. Solving the vehicle routing problem with drone for delivery services using an ant colony optimization algorithm［J］. Advanced Engineering Informatics， 2022， 51： 101536.

［13］Di Puglia Pugliese L， Guerriero F. Last-mile deliveries by using drones and classical vehicles［C］//Proc of Optimization and Decision Science： Methodologies and Applications.Cham：Springer， 2017： 557-565.

［14］戚遠(yuǎn)航， 蔡延光， 蔡顥， 等. 帶容量約束的供應(yīng)鏈物流運輸調(diào)度問題的雙層變鄰域蝙蝠算法［J］. 電子學(xué)報， 2019，47（7）： 1434-1442. （Qi Yuanhang， Cai Yanguang， Cai Hao， et al. Two-layer variable neighborhood bat algorithm for supply chain logistics transportation scheduling problem with capacity constraints［J］. Journal of Electronics， 2019， 47（7）： 1434-1442.）

［15］Liang Baoxian， Zhao Yunlong， Li Yang. A hybrid particle swarm optimization with crisscross learning strategy［J］. Engineering Applications of Artificial Intelligence， 2021， 105： 104418.

［16］張碩， 錢曉明， 樓佩煌， 等. 基于改進(jìn)粒子群算法的大規(guī)模自動導(dǎo)引車系統(tǒng)路徑規(guī)劃優(yōu)化［J］. 計算機集成制造系統(tǒng)， 2020， 26（9）： 2484-2496. （Zhang Shuo， Qian Xiaoming， Lou Peihuang， et al. Path planning optimization of large-scale automatic guided vehicle system based on improved particle swarm optimization［J］. Computer Integrated Manufacturing System， 2020， 26（9）： 2484-2496.）

［17］馬俊， 張紀(jì)會， 郭乙運. 考慮客戶分類的隨機時間車輛路徑優(yōu)化模型與算法［J］. 計算機應(yīng)用研究， 2022， 39（7）： 1979-1984. （Ma Jun， Zhang Jihui， Guo Yiyun. Random time vehicle routing optimization model and algorithm considering customer classification［J］. Application Research of Computers， 2022， 39（7）： 1979-1984.）

［18］Qi Yuanhang， Cai Yanguang. Hybrid chaotic discrete bat algorithm with variable neighborhood search for vehicle routing problem in complex supply chain［J］. Applied Sciences-Basel， 2021， 11（21）： 10101.