，，，,2

(1.浙江工業(yè)大學 特種裝備制造與先進加工技術教育部重點實驗室，浙江 杭州 310014；2.嘉興職業(yè)技術學院 機電與汽車分院，浙江 嘉興314036)

隨著全球氣候變暖，低碳經(jīng)濟逐漸成為世界各國關注的焦點，而物流業(yè)是碳排放的主要來源之一，因此低碳環(huán)境下的選址路徑問題(Location routing problem，LRP)受到學者的關注.在模型上，曹劍東等[1-2]以總成本最小化為優(yōu)化目標研究了考慮碳排放因素的車輛路徑問題(Vehicle routing problem，VRP).在實際的物流配送系統(tǒng)中，車輛路徑的優(yōu)化目標往往不止1個，可能同時考慮包括最短配送路徑、最少總成本、最少碳排放量和最少車輛數(shù)等因素.Jemai等[3-4]研究了最小化運輸路徑和碳排放的多目標VRP模型；魯建廈等[5-6]研究了最小化成本和車輛數(shù)的LRP模型.然而在現(xiàn)實生活中，物流企業(yè)一方面在配送時要盡量降低運營成本，另一方面還要減少碳排放、增加客戶滿意度.客戶的滿意度主要體現(xiàn)在車輛的準時到達，車輛在客戶要求的時間內(nèi)到達，則客戶的滿意度最高，提前或者推遲到達都會引起客戶滿意度的下降.在求解算法上，求解多目標VRP或LRP模型的算法基本是線性權重處理法或者是NSGA-II算法等傳統(tǒng)的啟發(fā)式算法，不同LRP變種問題的特殊性限制了這類算法的通用性，對于不同LRP問題都必須設計相對應的搜索策略，所以設計一種通用性算法使其較好地運用于一類相似的問題是目前研究的熱點，超啟發(fā)式算法就是這樣的方法之一.

基于此，在傳統(tǒng)LRP基礎上，研究考慮碳排放帶時間窗的選址路徑新模型，在優(yōu)化配送中心位置和數(shù)量的同時確定最佳的配送路線，降低配送過程中的碳排放，建立了包含碳排放、配送成本和客戶滿意度的多目標LCLRPTW(Low-carbon location-routing problem with time windows)模型，并結合禁忌搜索設計超啟發(fā)式算法進行求解.

1 選址-路徑問題的描述及數(shù)學規(guī)劃模型

1.1 問題描述

建立多目標LCLRPTW模型描述如下：有M個配送中心，每個配送中心擁有相同容量的K輛車，車輛最大載重量為Qk，有n個客戶的運輸任務需要完成，第i個客戶節(jié)點的需求量為di，且有maxdi

1.2 碳排放量的計算

由于要滿足同一條路線上各個客戶的需求，車輛在運輸過程中的裝載量會發(fā)生變化，所以CO2的排放量不是一成不變的，除了行駛距離外還需要考慮每一段運輸過程中實際的貨物裝載量.根據(jù)文獻[8]，CO2排放量和燃油消耗量成正比例關系，CO2排放量與行駛距離成正比例關系，與車輛裝載量成正線性相關.在計算車輛燃油消耗和CO2排放量的時候，只考慮車輛行駛距離和裝載量對兩者的影響，暫不考慮行駛速度、交通路況和天氣等因素的影響.燃油消耗量和CO2排放量的計算公式為

F=GD(aL+b)

(1)

式中：F為運輸過程的燃油消耗量；G為地形坡度因子；D為車輛的行駛距離；L為載貨重量；a,b分別為燃油消耗參數(shù).

計算出燃油消耗量，需要轉化為CO2排放量，同樣根據(jù)文獻[8]有

ECO2=Fη

(2)

式中：ECO2為CO2排放量；η為燃油轉換系數(shù).

為了驗證模型的有效性，取地形坡度因子G為單位1，根據(jù)實際問題可以得到

Eij=ηDij[aQij+b]

(3)

式中:Eij為車輛從客戶i后到客戶j的CO2排放量；Dij為客戶i到客戶j的行駛距離；Qij為離開客戶i后前往客戶j的車輛裝載的將要派送的客戶的貨物總量.

1.3 數(shù)學模型

對模型中的變量進行如下定義：M{m|m=1,2,…,M}為一系列配送中心；C{i|i=1,2,…,n}為一系列需要服務的客戶；V{k|k=1,2,…,K}為屬于各個配送中心的車輛；S{M∪C}為配送中心和客戶總和；Km為屬于配送中心m且同一車型的車輛；di為客戶i的需求；Cc為單位CO2排放成本；Cv為派車成本；Cr為單位車輛運輸成本；Cd為單位車輛折舊成本；Cf為單位燃油成本；Cm為配送中心開放成本；Qk為車輛的容量；Qm為配送中心的容量；[Ei,Li]為時間窗；T0為配送中心開始送貨時間；Tik為車輛k到達客戶i的時間；tij為客戶i～j的行駛時間；sik為車輛k服務客戶i的時間；α為單位時間車輛提前到達客戶點的懲罰系數(shù)；β為單位時間車輛延誤到達客戶點的懲罰系數(shù)；決策變量1：Xijk，當車輛k從客戶i～j為1，否則為0；決策變量2：Zm，配送中心開放為1，否則為0.

LCLRPTW數(shù)學模型如下：

(4)

(5)

(6)

(7)

(8)

(9)

(10)

(11)

(12)

(13)

(14)

(15)

T0=0

(16)

Tik+sik+tij=Tjk?i,j∈C,k∈Km

(17)

ai≤sik≤bi?i∈S,k∈Km

(18)

Xijk∈{0,1} ?i,j∈S,k∈Km

(19)

Zm∈{0,1} ?m∈M

(20)

其中：式(4～6)為3個目標函數(shù)；式(4)代表總碳排放量最少；式(5)代表總成本最少，第1項為配送中心開放成本，第2項為派車成本，第3項為車輛運輸成本和汽車折舊成本，第4項為燃油成本，第5項為CO2排放成本，第6項為早到等待成本，第7項為超時成本；式(6)代表等待和遲到的總時間最少；式(7)保證每個客戶均被訪問1次；式(8，9)表明車輛從配送中心出發(fā)，必須回到原配送中心；式(10)保證每1個運輸車輛的路徑最多從1個配送中心駛出；式(11)保證任何2個配送中心的車輛不會在同1條路徑上；式(12)保證訪問完客戶后必須離開；式(13)保證每個配送中心訪問的顧客總需求小于配送中心的容量；式(14)保證車的載重不大于它的載重能力；式(15)保證每個客戶的需求均被滿足；式(16)表示規(guī)定開始的送貨時間為0；式(17)表示車輛k到達客戶j處的時刻等于車輛k到達客戶i處的時刻與其在客戶i處服務時間，以及從客戶i處到客戶j處的行駛時間之和；式(18)表示車輛k到達客戶i處的時間需滿足客戶規(guī)定的時間窗；式(19，20)是對決策變量的描述.

2 求解算法設計

2.1 超啟發(fā)式算法

超啟發(fā)式算法(Hyper-heuristicalgorithm，HHA)是近年來發(fā)展起來的一種新型啟發(fā)式算法，旨在提高算法的通用性[9-10].超啟發(fā)式算法提供一種高層啟發(fā)式策略，通過管理或操縱一系列底層啟發(fā)式算子，用于求解各種組合優(yōu)化問題.

圖1給出了超啟發(fā)式算法的框架，分為2個部分[9]：在控制域部分，由智能計算專家進行設計高層啟發(fā)式策略，包括如何利用底層啟發(fā)式算子構造可行解或改進解的質量；在問題域部分，由應用領域專家提供一系列底層啟發(fā)式算子(Low-level heuristic,LLH)、問題定義和目標函數(shù)等信息.由于2個部分之間實現(xiàn)了領域屏蔽，只要修改問題域的LLH、問題定義和目標函數(shù)等信息，一種超啟發(fā)式算法可以方便地移植到新的問題上.超啟發(fā)式算法已經(jīng)被廣泛運用在多種領域，大多用于求解排課問題[11]，流水車間調(diào)度問題[12]，裝箱問題[13]等約束較少的組合優(yōu)化問題，對于更為復雜的問題實際應用潛力巨大.因此提出一種禁忌搜索的超啟發(fā)式算法對多目標LCLRPTW模型進行求解，為超啟發(fā)式算法進一步在LRP領域的應用提供參考.

圖1 超啟發(fā)式算法框架Fig.1 Framework of hyper-heuristic algorithm

2.2 基于禁忌搜索的超啟發(fā)式算法

根據(jù)超啟發(fā)式算法的特點，算法設計需要考慮問題域中底層啟發(fā)式算子的構建和控制域中如何設計禁忌搜索作為高層啟發(fā)式策略.

2.2.1 底層啟發(fā)式算子

在研究的LCLRPTW模型中，底層啟發(fā)式算子根據(jù)文獻[14]的分類標準，構建變異算子、破壞與重構算子、局部搜索算子和交叉算子.算子具體詳情如下：

1) 變異算子.變異算子是一種對解產(chǎn)生微小擾動的突變算子.所有這些算子都將返回任何滿足約束的路徑，無論目標函數(shù)值是改進或惡化.LLH1～ LLH4為路徑變異算子；LLH5，LLH6為配送中心選址變異算子.

LLH1：2-opt.選擇1條路徑，交換相鄰兩客戶點的位置.

LLH2：Or-opt.選擇1條路徑，把相鄰2客戶點插入到其他位置.

LLH3：Interchange.選擇2條路徑，交換任意2個客戶點位置.

LLH4：Shift.選擇2條路徑，1條路徑上的任意1個客戶點插入到另1條路徑合適的位置.

LLH5：Shift.選擇1條路徑，使用1個新的配送中心替換此路徑的配送中心.

LLH6：Interchange.選擇2條路徑，交換2個配送中心.

2) 破壞與重構算子.破壞與重構算子共有3類，分別為徑向破壞與重構算子、隨機破壞與重構算子和序列破壞與重構算子[15].使用第1類算子，即根據(jù)與“基準客戶”的接近程度從路徑中刪除一些客戶.

LLH7：Location-based radial ruin.隨機選擇1個客戶作為基準客戶，根據(jù)客戶位置接近基準客戶位置的原則，以[1%,10%]的概率將客戶從路徑中刪除.

LLH8：Time-based radial ruin，在整個調(diào)度期間內(nèi)隨機選擇1個時間，根據(jù)其時間窗口的接近時間以[25%,75%]的概率從解決方案中刪除客戶.

3) 局部搜索算子.與變異算子一樣，經(jīng)常是以交換或插入的形式產(chǎn)生新路徑.但與變異算子不同的是局部搜索算子只會返回1個改進解.

LLH9：Interchange.與LLH3一樣，不同的是此算子只接受改進解.

LLH10：Shift.與LLH4一樣，不同的是此算子根據(jù)衡量準則來選擇客戶，被選客戶插入到另一條路徑的任意位置，且只接受改進解.

LLH11：2-opt*.選擇2條路徑，根據(jù)衡量準則選擇交換位置點，交換此位置后的所有客戶點，且只接受改進解.

LLH12：GENI.計算不同路線上的任意2個客戶的距離，采用最短的距離作為基準距離，選擇移除距離最接近基準距離的客戶，且只接受改進解.

4) 交叉算子.選擇2個父代路徑作為輸入，生成1條子代路徑的方法，通?？梢酝ㄟ^1點、2點和均勻交叉等方法完成.

LLH13：Combine，選擇2個父代，以[25%,75%]的概率將1個父代復制得到1條子路徑， 添加來自另一個父代中非沖突的路徑，并隨機插入剩余的客戶.

LLH14：Longest Combine，選擇2個父代，所有路徑以服務客戶點數(shù)量的降序來考慮，任何不重復客戶的路徑都會被添加生成1條子路徑，并隨機插入剩余的客戶.

2.2.2 高層啟發(fā)式策略

高層啟發(fā)式策略主要是考慮對底層啟發(fā)式算子的選擇策略，即決定在1次迭代過程中選擇哪些底層啟發(fā)式算子用于改進當前解.

基于禁忌搜索的高層啟發(fā)式策略使用得分的方法來評價底層啟發(fā)式算子的表現(xiàn)，再依據(jù)分數(shù)決定在迭代過程中使用哪些底層啟發(fā)式算子改進當前解.每個底層啟發(fā)式算子均有1個相同的初始分數(shù)，采用強化學習的方法進行分數(shù)更新，當算子改進當前解時，給算子加分，反之則減分.最后根據(jù)分數(shù)，通過貪心法選擇底層啟發(fā)式算子，即總是選擇分數(shù)最高的算子.具體過程：在搜索開始時，每個底層啟發(fā)式算子k均有1個分數(shù)，設rk=0(假設初始分數(shù)均為0分).每個算子的分數(shù)允許在間隔[rmin,rmax]中減少和增加，rmin和rmax分別為最低和最高的分數(shù).每次迭代選擇1個不在禁忌表內(nèi)且分數(shù)最高的算子k，計算目標函數(shù)值f1并求得與上一個目標函數(shù)值f2之間的變化量，即Δ=f2-f1.如果結果有改進(若目標為最小化問題，即Δ>0)，則增加算子k的分數(shù)，例如rk=rk+α，其中α為1個正數(shù).否則減少分數(shù)，例如rk=rk-α，并把算子k放入固定禁忌長度的禁忌表內(nèi)，根據(jù)先進先出原則赦免禁忌表內(nèi)的算子.

對于多目標問題，對高層啟發(fā)式策略作如下設計：

1) 由于底層啟發(fā)式算子的性能不再只針對1個目標進行評估，而是針對每個目標進行評估，因此目標函數(shù)變化量Δ變?yōu)棣，算子分數(shù)rk變?yōu)閞k(u)，其中u=1,2,3(3為目標個數(shù)).

2) 由于上述改變，因此在每次迭代過程中以均等的概率隨機選取1個目標函數(shù)進行求解.

2.3 求解LCLRPTW過程

LCLRPTW模型采用首先進行配送中心選址分配，再進行路徑優(yōu)化的方法進行求解.即首先設計一種重心法選址方法進行配送中心初始定位，以確定配送中心位置和其服務的客戶群；然后求解不同客戶群的配送路徑.具體流程如下：

步驟1初始化種群.每條染色體均采用客戶直接編碼的方式：假設有9個客戶點和有3個候選配送中心，其中1～9表示客戶點，0(1)-0(3)表示候選配送中心.對這9個客戶進行隨機排列，假設生成的1條染色體2-5-6-4-3-7-8-1-9，隨機生成100條染色體構成初始種群P.

步驟2將種群中每個個體轉換成可行路徑.按照約束條件(此處為車輛容量和配送中心庫存)，依次將解的每個客戶納入到車輛路徑中，當1個客戶違背約束時，就重新開辟1條新的路徑并將這個客戶納入該路徑.例如0-2-5-6-0-4-3-0-7-8-1-9-0，代表著共有3條配送路徑0-2-5-6-0，0-4-3-0，0-7-8-1-9-0.類似地，若每1條染色體的每1段路徑均可行，則該染色體為有效染色體.

步驟3進行配送中心初始化選址.對于多個配送中心選址問題，為了使模型簡單化，很多時候可以將它變成多個單一配送中心選址問題來處理[16].結合重心法的思想，方案如下：1) 根據(jù)上述生成的3條配送路徑，由客戶坐標和需求計算出各條路徑的重心；2) 根據(jù)求得重心，根據(jù)距離最小化原則，給每條路徑分配一個配送中心；3) 形成初始化選址方案.可以得到3條路線：0(2)-2-5-6-0(2)，0(1)-4-3-0(1)，0(3)-7-8-1-9-0(3).

步驟4通過2.2.2中設計的高層啟發(fā)式策略對初始種群P中每條染色體的路徑進行優(yōu)化.

1) 設置相關算法參數(shù).包括底層啟發(fā)式算子的初始分數(shù)，加減分值α，禁忌長度，并令禁忌表為空.

2) 選擇目標函數(shù).在3個目標函數(shù)中(u,v,w)，以均等的概率隨機選擇1個目標函數(shù)u.

3) 選擇底層啟發(fā)式算子.根據(jù)底層啟發(fā)式算子的分數(shù)，選擇1個不在禁忌表內(nèi)且分數(shù)最高的算子k(若選中LLH5和LLH6這2個算子，即代表對配送中心進行操作).

4) 計算目標函數(shù)值.使用3)選擇的算子k計算2)所選目標u的函數(shù)值.

5) 計算目標函數(shù)值的變化量Δu.如果Δu>0，則算子k分數(shù)增加，即rk(u)=rk(u)+α,其中α=1，否則分數(shù)減少，即rk(u)=rk(u)-α，并把算子k放入禁忌表中，根據(jù)先進先出原則赦免禁忌表內(nèi)的算子.

6) 計算目標函數(shù)v，w，如果Δv>0，則rk(v)=rk(v)+α，否則rk(v)=rk(v)-α；如果Δw>0，則rk(w)=rk(w)+α，否則rk(w)=rk(w)-α.

7) 進行非支配排序，更新非支配解集.

8) 若滿足算法終止準則結束算法，否則轉到2).

9) 產(chǎn)生非支配解集.

3 應用實例及結果分析

3.1 應用實例

為了體現(xiàn)算法的實用性，采用某地區(qū)的一家物流公司作為實際案例進行求解[17]，該案例包括3個配送中心和30個客戶，其中每個客戶的坐標、需求量、時間窗以及服務時間均已知.客戶點詳細信息見文獻[17]表5.2，配送中心相關信息如表1所示.基于禁忌搜索的多目標超啟發(fā)式算法(Multi-objective hyper heuristic algorithm based on tabu search，MOHH_TS)以及進行對比的NSGA-II算法[18]均采用MATLAB編程，運行在CPU為Intel Core i5，2.6 GHz，RAM為4 G內(nèi)存的計算機平臺上，其他參數(shù)設置如下[8,17]：初始種群數(shù)為100，燃油消耗參數(shù)a=6.208×10-3，b=0.212 5，燃油轉換系數(shù)η=2.68，車輛行駛速度v=50 km/h，單位時間車輛提前到達或延誤到達客戶點的懲罰系數(shù)分別為α=0.5，β=1.

表1 配送中心相關信息Table 1 Related information of depot centers

3.2 結果分析

首先對MOHH_TS算法和NSGA-II算法在解決多目標LCLRPTW問題的性能進行比較，通過收斂代數(shù)、Pareto面和目標函數(shù)的收斂情況3個方面進行分析.

3.2.1 收斂代數(shù)分析

求解多目標LCLRPTW問題的收斂代數(shù)見表2，取運行10次的平均收斂代數(shù).收斂準則為連續(xù)10代不再出現(xiàn)更優(yōu)個體.

表2 MOHH_TS和NSGA-II的收斂代數(shù)Table 2 Convergence times of MOHH_TS and NSGA-II 次

由表2可以看出：MOHH_TS和NSGA-II算法均可以找到最優(yōu)解，但MOHH_TS算法的平均迭代數(shù)為183.2次少于NSGA-II算法的232.6次，說明MOHH_TS算法相比NSGA-II算法收斂速度更快，效率更高，這在求解更大規(guī)模問題時將會更加明顯；同時，MOHH_TS算法相比較于NSGA-II算法的收斂代數(shù)曲線更加平穩(wěn)，說明MOHH_TS算法更加穩(wěn)定.

3.2.2 Pareto面分析

為了能夠在相同的平臺進行比較，設置2個算法均運行300代，得到的Pareto面如圖2所示(圓圈代

表MOHH_TS算法，星號代表NSGA-II算法).

由圖2可以看到：MOHH_TS算法得到的Pareto數(shù)量更多，且更容易達到每個目標函數(shù)的極值點；而NSGA-II算法得到的Pareto解分布得較為分散，對于3個目標函數(shù)也不能均等對待.同時，從圖3中各個Pareto解的位置分布情況來看MOHH_TS算法的解幾乎落在同1個Pareto面上[19]，說明MOHH_TS算法收斂得較好.

圖2 MOHH_TS和NSGA-II的Pareto最優(yōu)解Fig.2 Pareto optimal solution of MOHH_TS and NSGA-II

3.2.3 目標函數(shù)分析

為了清楚地比較2個算法的各個目標函數(shù)的優(yōu)劣情況，把所有代中各個目標函數(shù)的平均值進行比較，見圖3(a～c)(粗線代表MOHH_TS算法，細線代表NSGA-II算法).

圖3 目標函數(shù)收斂情況比較Fig.3 Comparison of objective function convergence

由圖3可知：MOHH_TS算法的3個目標函數(shù)值好于NSGA-II算法，且曲線初期下降速度明顯，說明算法收斂速度更快；MOHH_TS算法的末端值更小，說明算法更有效，能更好地避免早熟現(xiàn)象，能找到更接近最優(yōu)解的Pareto解.

其次為了了解MOHH_TS算法中底層啟發(fā)式算子對解質量的影響情況，運行300代，運行10次，統(tǒng)計14個底層啟發(fā)式算子的平均接受率，如圖4所示.

圖4 算子的平均接受率Fig.4 Average acceptance rate of LLH

由圖4可以看出：14個底層啟發(fā)式算子的平均接受率均超過0.6，其中算子LLH2，LLH3，LLH8，LLH13的平均接受率相對較高，超過了0.8，對解質量的改進作用較大；而算子LLH5和LLH6的平均接受率相對較低，不超過0.7.4個局部搜索算子(LLH9～LLH12)的接受率為1，這是由局部搜索算子本身的性能所決定的，因為產(chǎn)生的所有候選解為改進解均會被接受.總體來說，14個底層啟發(fā)式算子均能對解有不錯的改進作用，證明了算子構建的有效性.

最后對于多目標優(yōu)化問題，決策者往往從實際需要出發(fā)選擇最為合適的解，從所得的Pareto解集中分別選擇碳排放最少、總成本最少和客戶滿意度最高(等待和遲到的總時間最少)的3組解作為配送方案，3組最優(yōu)解如表3所示.

表3 最優(yōu)解Table 3 Optimal solutions

對應的最優(yōu)配送路線如表4所示.

表4 最優(yōu)路線Table 4 Optimal routes

由此可見，無論是從算法的收斂代數(shù)、收斂速度和穩(wěn)定性，還是從最優(yōu)解來看，MOHH_TS算法可以有效地解決多目標LCLRPTW模型.

4 結 論

針對物流配送中選址路徑問題，在車輛路徑安排時考慮碳排放的因素，首先構建了同時考慮碳排放、總成本和客戶滿意度的多目標LCLRPTW模型；其次構建了和LCLRPTW問題相關的底層啟發(fā)式算子，提出了基于禁忌搜索的超啟發(fā)式算法；最后進行仿真實驗，并和NSGA-II算法進行比較.從底層啟發(fā)式算子的平均接受率來看，證明了算子構建的有效性.從算法的收斂代數(shù)、Pareto面和目標函數(shù)的收斂情況來看，證明了MOHH_TS算法可以有效地解決所提出的多目標LCLRPTW模型，根據(jù)得到的Pareto最優(yōu)解也有利于決策者從實際需要出發(fā)選擇最為合適的配送方案.由于超啟發(fā)式算法具有較好的通用性，可以將該算法應用到其他組合優(yōu)化問題上，而不需改變框架，只需提供1組與問題相關的底層啟發(fā)式算子.