廣東工業(yè)大學(xué)自動(dòng)化學(xué)院 李旭陽(yáng) 蔡延光

針對(duì)帶時(shí)間窗的物流運(yùn)輸調(diào)度問題，設(shè)計(jì)一種改進(jìn)的和聲搜索算法。該算法利用類電磁機(jī)制算法改進(jìn)和聲搜索的隨機(jī)產(chǎn)生規(guī)則，并且使用了和聲記憶庫(kù)擾動(dòng)策略和2-Opt局部搜索策略提高算法性能。結(jié)果表明：相比基本和聲搜索算法及其他啟發(fā)式算法，所設(shè)計(jì)的算法具有更好的收斂速度和收斂精度。

物流運(yùn)輸調(diào)度問題國(guó)外一般稱為車輛路徑規(guī)劃問題（Vehicle Routing Problem，VRP），該問題自提出以來(lái)就一直是研究的熱點(diǎn)。戚遠(yuǎn)航等提出了一種雙層變鄰域蝙蝠算法求解帶容量約束的物流問題，該算法采用變鄰域局部搜索策略加強(qiáng)算法的尋優(yōu)能力，很好的解決了帶容量約束的物流運(yùn)輸調(diào)度問題；鄧麗娟等提出了一種混合蟻群算法求解帶時(shí)間窗的VRP，該算法探索兩個(gè)目標(biāo)函數(shù)，獲得了很好的實(shí)驗(yàn)效果；蔡延光等提出了一種帶遺傳算子的自適應(yīng)蟻群算法，求解帶軟時(shí)間窗的車輛路徑問題。

2001年Geem ZW等人源于音樂的創(chuàng)作，提出了和聲搜索（Harmony Search，HS）算法，一種新的元啟發(fā)算法；高立群等把粒子群算法與和聲搜索算法相結(jié)合，提出自適應(yīng)和聲粒子群搜索算法；歐陽(yáng)海濱等研究了在非對(duì)稱區(qū)間內(nèi)，證明了和聲搜索算法的參數(shù)與即興創(chuàng)作過程的探索之間的關(guān)系，從而證明了算法的迭代收斂性，提出了一種改進(jìn)的和聲搜索算法；王艷等采用改進(jìn)的和聲搜索算法，用一種離散編碼方式，解決了車間調(diào)度問題。本文針對(duì)帶時(shí)間窗的物流運(yùn)輸調(diào)度問題，設(shè)計(jì)了改進(jìn)的和聲搜索算法。該算法利用類電磁機(jī)制算法改進(jìn)和聲搜索算法的隨機(jī)生成的規(guī)則，并且使用了和聲記憶庫(kù)擾動(dòng)策略和2-Opt局部搜索策略提高算法性能。實(shí)驗(yàn)結(jié)果表明：該算法比基本和聲搜索算法有更好的收斂速度和精度。

1 問題描述

帶時(shí)間窗物流運(yùn)輸調(diào)度(Vehicle Routing Problem with Time Windows，VRPTW)問題可描述為：某車場(chǎng)需要為N個(gè)客戶提供貨物配送服務(wù)，并要求在[Ei,Li]所表示的時(shí)間窗內(nèi)將貨物送達(dá)，其中Ei表示客戶i要求的最早貨物到達(dá)時(shí)間，Li表示客戶i要求的最晚貨物到達(dá)時(shí)間。若未按時(shí)送達(dá)則增加相應(yīng)的時(shí)間懲罰成本，用C1（單位：元/min）表示早到的懲罰系數(shù)，C2（單位：元/ min）表示晚到的懲罰系數(shù)；tik表示車輛k到達(dá)客戶i的實(shí)際時(shí)間；UTik表示運(yùn)輸車輛k在客戶i卸貨所需的時(shí)間；RTijk表示運(yùn)輸車輛k從客戶i行駛到客戶j所需的時(shí)間；每條客戶的貨物需求量為qi(i=1,2,3,...N)，該倉(cāng)庫(kù)有K輛車且車型相同，車的最大載重為W，客戶i到客戶j的距離為dij，其中i=0表示倉(cāng)庫(kù)，定義決策變量：，當(dāng)k經(jīng)過客戶i駛向客戶j時(shí)為1，否則為0，；，當(dāng)客戶i由車輛滿足時(shí)為1，否則為0。并建立數(shù)學(xué)模型如下：

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

約束條件：

式(1)為VRPTW的目標(biāo)函數(shù)，式(2)(3)表示封閉式的車輛路徑規(guī)劃，式(4)(5)(6)表示每條客戶僅要一輛車完成任務(wù)，式(7)(8)表示車的最大容量和行駛限制，式(9)表示參與任務(wù)的車輛數(shù)不能超過車場(chǎng)所擁有的車輛總數(shù)，式(10)(11)表示保證運(yùn)輸車輛在配送貨物在時(shí)時(shí)間上的連續(xù)性，避免時(shí)間沖突。

2 算法設(shè)計(jì)

2.1 和聲搜索算法

2.1.1 基本變量

（1）和聲記憶庫(kù)大小HMS：指和聲記憶庫(kù)中和聲的數(shù)量。

（2）記憶庫(kù)取值概率HMCR：在每次迭代中，有一組和聲將從庫(kù)中被選擇為一個(gè)新和聲的特定概率。

（3）微調(diào)概率PAR：選擇一組和聲需要用微調(diào)概率判斷是否要微調(diào)。如果PAR值過小，算法的優(yōu)化范圍減小，而PAR值過大，搜索算法會(huì)變成隨機(jī)選擇。

（4）音調(diào)微調(diào)帶寬 bw：指微調(diào)時(shí)的調(diào)整幅度。

（5）創(chuàng)作的次數(shù) Tmax：即算法的迭代次數(shù)。

根據(jù)相關(guān)研究，算法參數(shù)為HMS值可取10～50，HMCR可取0.70～0.95，PAR可取0.20～0.50。

2.1.2 和聲搜索算法步驟

Step 1：初始化算法，指定和聲搜索算法各項(xiàng)基本參數(shù)：HMS、HMCR、PAR、bw、Tmax。

Step 2：根據(jù)HMS生成和聲記憶庫(kù)，并計(jì)算目標(biāo)函數(shù)的值，組成如下所示的和聲記憶庫(kù)HM：Step 3：產(chǎn)生新的和聲X'=[x'1,x'

2,...,x'N]。新的和聲通過三種方式生成：①由HM隨機(jī)產(chǎn)生。②音調(diào)微調(diào)。③隨機(jī)選擇，其生成過程可表示為：

IF rand(0,1)

ELSE

END IF

其中，rand(0,1)是在0和1之間的隨機(jī)變量值，xi,min,xi,max分別是決策變量xi的最小、最大值。若新的和聲來(lái)自和聲記憶庫(kù)HM，則用微調(diào)概率來(lái)判斷是否要微調(diào)規(guī)則：

Step 4：更新記憶庫(kù)。若Step3得到的一組新的和聲比HM中最差和聲要好，則用新和聲替換最差和聲。

Step 5：重復(fù)Step3和Step4，直到滿足最大迭代次數(shù)Tmax。

2.2 類電磁機(jī)制算法

在物理學(xué)中，帶電粒子之間存在著相互吸引或相互排斥的庫(kù)侖力。Birbil和Fang受到這一現(xiàn)象的啟發(fā)，于2003年提出了類電磁機(jī)制（Electromagnetism-like Mechanism，EM）算法，并且用該算法成功解決了無(wú)約束優(yōu)化問題。

EM算法將待優(yōu)化問題的每個(gè)解模擬為電磁場(chǎng)中的帶電粒子，并根據(jù)解的適合度來(lái)確定每個(gè)粒子的電荷。每個(gè)帶電粒子會(huì)被其他粒子吸引或排斥，這取決于它們所攜帶的電荷和它們之間的距離。根據(jù)粒子受到的合力的大小，粒子會(huì)沿著合力的方向移動(dòng)一定的距離。由于算法是對(duì)電磁原理的一種類比，兩者之間并不完全相同，因此稱為類電磁機(jī)制。

EM算法主要由以下四個(gè)步驟組成：

（1）初始化。將一定數(shù)量的粒子隨機(jī)、均勻地分布在待優(yōu)化問題的可行解域里，然后計(jì)算出每個(gè)粒子所在位置的適應(yīng)度值，并將適應(yīng)度最好的粒子的位置記為Xbest。

（2）局部搜索。EM算法中的局部搜索是指以粒子當(dāng)前的位置為中心，按照一定的步長(zhǎng)搜索解空間的每個(gè)維度。如果找到了更好的解決方案，搜索將被終止，并執(zhí)行下一步。

（3）計(jì)算合力。計(jì)算合力是為了將EM算法的局部搜索能力和全局搜索能力相結(jié)合，提高算法的求解能力，并為下一步更新粒子位置提供相應(yīng)的信息。

粒子i所帶電荷量的多少qi決定了粒子i所受的庫(kù)倫力的大小，電荷量qi的計(jì)算公式為：

式中，n為待優(yōu)化問題的維度，m為解空間中帶電粒子的個(gè)數(shù)。由式(16)可知，適應(yīng)度越好的粒子所帶的電荷數(shù)越大，對(duì)其他帶電粒子的吸引或排斥能力也就越強(qiáng)。一個(gè)粒子受另一個(gè)粒子的庫(kù)侖力的方向由兩個(gè)粒子的適應(yīng)度好壞決定。在計(jì)算粒子i受到解空間中其他每個(gè)粒子的力的大小和方向后，即可通過累加方法得到作用在粒子i上的合力Fi：

（4）移動(dòng)粒子

計(jì)算完合力Fi之后，粒子i將沿著Fi的方向移動(dòng)，移動(dòng)后粒子i的位置為：

式中，步長(zhǎng)λ在[0,1]上服從矩形分布；Range為一個(gè)向量，向量中的每一個(gè)分量表示粒子在解空間的對(duì)應(yīng)維度上能夠移動(dòng)的距離。同時(shí)，為了保證每個(gè)粒子移動(dòng)后的位置不會(huì)超出解空間的范圍，EM算法對(duì)粒子所受的合力做了無(wú)量綱化處理。

這樣，粒子的位置得到更新，也就完成了EM算法的一次迭代。

2.3 和聲記憶庫(kù)擾動(dòng)策略

當(dāng)和聲搜索算法運(yùn)算一定的次數(shù)后，和聲記憶庫(kù)中的和聲多樣性水平降低，甚至出現(xiàn)大部分和聲的適應(yīng)度都相等的情形。此時(shí)，算法將會(huì)限入局部最優(yōu)解。這個(gè)時(shí)候如果能夠?qū)吐曈洃泿?kù)中的一部分和聲進(jìn)行擾動(dòng)，以提高和聲記憶庫(kù)的和聲多樣性水平，將會(huì)有利于算法跳出局部最優(yōu)解繼續(xù)尋優(yōu)，從而增加得到全局最優(yōu)解的概率。

2.3.1 擾動(dòng)時(shí)機(jī)

和聲記憶庫(kù)的擾動(dòng)需要在算法迭代至一定次數(shù)，且和聲記憶庫(kù)中的和聲相似程度達(dá)到一定閾值時(shí)進(jìn)行。設(shè)和聲記憶庫(kù)大小為HMS，算法最大迭代次數(shù)為Tmax，當(dāng)前迭代次數(shù)為gn，則進(jìn)行和聲記憶庫(kù)擾動(dòng)的時(shí)機(jī)需要同時(shí)滿足如下兩個(gè)條件：

②和聲記憶庫(kù)連續(xù)t次迭代未更新：

2.3.2 擾動(dòng)步驟

當(dāng)滿足擾動(dòng)時(shí)機(jī)時(shí)，擾動(dòng)步驟如下：

Step1：對(duì)HM中的和聲按照適應(yīng)度值排序，隨機(jī)選取(HMS/2)組非最優(yōu)和聲；

Step2：對(duì)Step1中選擇的(HMS/2)組非最優(yōu)和聲進(jìn)行混沌擾動(dòng)；

Step3：計(jì)算生成的新和聲的適應(yīng)度值，若優(yōu)于和聲記憶庫(kù)中的最差和聲，則用新和聲替換最差和聲，更新和聲記憶庫(kù)；否則需要重新擾動(dòng)。

2.4 局部搜索策略

在求解物流運(yùn)輸調(diào)度問題時(shí)，常用的局部搜索算法主要有2-Opt搜索，0-1搜索和1-1搜索等算法。本文采用2-Opt方法進(jìn)行局部搜索。一次2-Opt操作可描述為：在一條路徑中隨機(jī)選擇兩個(gè)點(diǎn)i和j，將i-1與j連接，i與j+1連接，并將i和j之間的路徑反向，若2-Opt操作后的路徑長(zhǎng)度比操作之前的長(zhǎng)度小，則更新路徑。2-Opt操作示意如圖1所示。

圖1 2-Opt操作示意圖

2.5 編碼策略

為實(shí)現(xiàn)VRPTW路線的解空間和和聲搜索算法的搜索空間的解空間之間的轉(zhuǎn)換，本文采用序列號(hào)編碼方案。對(duì)于n個(gè)客戶，k輛車的車輛路徑規(guī)劃問題，采用個(gè)n+k-1個(gè)的整數(shù)來(lái)編碼，編碼中的n位表示客戶，k-1位為車輛號(hào)，0表示配送中心。例如對(duì)于10個(gè)客戶，3輛車的車輛路徑規(guī)劃問題，有如下編碼：3，10，9，11，5，8，7，2，12，6，4，1，其中11和12為車輛標(biāo)識(shí)編碼，該編碼對(duì)應(yīng)和聲庫(kù)里面的一個(gè)和聲。該組編碼表示第一輛車行車路徑為0－3－10－9－0，第二輛車行車路徑為0－5－8－7－2－0，第三輛車行車路徑為0－6－4－1－0。這些編碼就組成了和聲記憶庫(kù)，接下來(lái)就是對(duì)庫(kù)中編碼尋找最優(yōu)解的過程。

2.6 算法步驟

將改進(jìn)策略整合到基本的和聲搜索算法，改進(jìn)的算法的具體步驟如下：

Step1：初始化算法各項(xiàng)參數(shù)；生成和聲記憶庫(kù)，并計(jì)算每個(gè)和聲的適應(yīng)度值，采用車輛的行駛總路程與時(shí)間窗懲罰之和的倒數(shù)作為適應(yīng)度值，適應(yīng)度越大，表示和聲越優(yōu)。

Step2：生成新和聲；

Step2-1：若生成的隨機(jī)數(shù)rand(0,1)

Step2-2：以當(dāng)前和聲記憶庫(kù)中的所有和聲作為類電磁機(jī)制算法中的粒子，以HM中的最優(yōu)和聲作為最優(yōu)粒子Xbest，根據(jù)式(16)與式(17)計(jì)算所有粒子的電荷量以及每個(gè)粒子受到的合力，并結(jié)合式(18)計(jì)算每個(gè)粒子移動(dòng)后的位置，選擇最優(yōu)的粒子賦值給Xnew，轉(zhuǎn)Step3；

Step2-3：若生成的隨機(jī)數(shù)rand(0,1)

Step3：2-Opt操作，將迭代次數(shù)gn加1。對(duì)Xnew進(jìn)行一次2-Opt操作，若操作后路徑長(zhǎng)度變短則更新Xnew；

Step4：更新和聲記憶庫(kù)。計(jì)算Xnew的適應(yīng)度值，如果優(yōu)于庫(kù)中的最差和聲，則使用Xnew替換最差和聲，否則，庫(kù)不變；

Step5：判斷是否滿足擾動(dòng)條件，若滿足則根據(jù)擾動(dòng)策略對(duì)HM進(jìn)行擾動(dòng)，否則轉(zhuǎn)Step6；

Step6：判斷gn是否達(dá)到最大迭代次數(shù)Tmax，若沒有則返回Step2，否則轉(zhuǎn)Step7；

Step7：輸出適應(yīng)度最好的和聲并輸出其對(duì)應(yīng)的適應(yīng)度值，該和聲即為最優(yōu)解。

3 實(shí)驗(yàn)與分析

3.1 實(shí)驗(yàn)數(shù)據(jù)

為驗(yàn)證算法的有效性，假設(shè)一個(gè)配送中心需要為14個(gè)客戶配送貨物，客戶坐標(biāo)隨機(jī)產(chǎn)生，各客戶所需的貨物量、時(shí)間需求如表1所示。設(shè)定車場(chǎng)坐標(biāo)(0,0)，編號(hào)為0擁有4輛運(yùn)輸車輛，每次最大載重量為15t，配送中心于上午8:00開始為各客戶運(yùn)輸，運(yùn)輸車輛勻速行駛，車速為2020km/h，運(yùn)輸車輛早到的懲罰系數(shù)為2元/min，晚到的懲罰系數(shù)為5元/min。

表1 客戶坐標(biāo)、需求量與時(shí)間窗

本章和聲搜索算法參數(shù)設(shè)置為：HMS=20，HMCR=0.75，PAR=0.3，Tmax=200。

3.2 實(shí)驗(yàn)結(jié)果與分析

算法運(yùn)行50次取平均結(jié)果，仿真結(jié)果與基本和聲搜索算法及其他三種算法的對(duì)比如表2所示，最優(yōu)配送線路如表3所示，最優(yōu)配送軌跡圖如圖2所示。

表2 結(jié)果對(duì)比

表3 最優(yōu)配送線路

圖2 最優(yōu)配送軌跡圖

由表2中的數(shù)據(jù)可知，基于類電磁機(jī)制的和聲搜索算法因?yàn)橐肓祟愲姶艡C(jī)制算法加快了收斂速度，并使用了和聲記憶庫(kù)擾動(dòng)策略和2-Opt局部搜索策略能夠有效的跳出局部最優(yōu)解，提高了算法的局部搜索能力，所以在50次運(yùn)算后，所得到的平均里程、平均搜索時(shí)間以及平均時(shí)間窗懲罰都要優(yōu)于基本和聲搜索算法所得到的平均結(jié)果，并且優(yōu)于其他幾種算法。

算法在迭代完成后，得到最優(yōu)配送路徑圖，配送車輛行駛路徑總長(zhǎng)度為91.65km，總時(shí)間窗懲罰為143元。其中，編號(hào)為1、2和4的配送車輛所執(zhí)行的配送任務(wù)中，分別出現(xiàn)59元、32元和52元的時(shí)間窗懲罰，而編號(hào)為3的運(yùn)輸車輛所執(zhí)行的配送任務(wù)沒有時(shí)間窗懲罰，即不存在早到和遲到的情況。由此可見，本文提出的改進(jìn)的和聲搜索算法能夠有效地解決帶時(shí)間窗的物流運(yùn)輸調(diào)度問題。

總結(jié)：本文設(shè)計(jì)了改進(jìn)的和聲搜索算法來(lái)解決帶時(shí)間窗的路徑規(guī)劃問題，該算法通過類電磁機(jī)制算法改進(jìn)和聲搜索算法的隨機(jī)產(chǎn)生規(guī)則，并且使用了和聲記憶庫(kù)擾動(dòng)策略和2-Opt局部搜索策略提高算法的搜索能力，性能得到了提高，具有更快的收斂速度和更高的精度。本文所提的模型在很大程度上貼合了實(shí)際情況，但是車輛過少，下一步準(zhǔn)備研究多車場(chǎng)下的模型，提出一種適合求解多車場(chǎng)下物流運(yùn)輸調(diào)度問題的算法。

基金項(xiàng)目：國(guó)家自然科學(xué)基金（61074147）；廣東省自然科學(xué)基金（S2011010005059）；廣東省教育部產(chǎn)學(xué)研結(jié)合項(xiàng)目（2012B091000171，2011B090400460）；廣東省科技計(jì)劃項(xiàng)目（2012B050600028，2014B010118004，2016A050502060）；廣州市花都區(qū)科技計(jì)劃項(xiàng)目（HD14ZD001）；廣州市科技計(jì)劃項(xiàng)目（201604016055）；廣州市天河區(qū)科技計(jì)劃項(xiàng)目（2018CX005）。