基于蟻群選擇超啟發(fā)算法的低碳選址—路徑問(wèn)題

王 舜，趙燕偉,冷龍龍，張春苗，蔣海青

(浙江工業(yè)大學(xué) 特種裝備制造與先進(jìn)加工技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室，浙江 杭州 310014)

0 引言

選址—路徑問(wèn)題(Location-Routing Problem, LRP)是一個(gè)經(jīng)典的組合優(yōu)化問(wèn)題，可以看作是配送中心選址分配問(wèn)題(Location Allocation Problem, LAP)與車(chē)輛路徑問(wèn)題(Vehicle Routing Problem, VRP)[1-2]的組合問(wèn)題，是結(jié)合了兩個(gè)復(fù)雜決策問(wèn)題的NP-hard問(wèn)題。配送中心選址與客戶地理位置、車(chē)輛配送路線規(guī)劃之間的相互依賴(lài)關(guān)系有利于幫助企業(yè)全面協(xié)調(diào)物流系統(tǒng)的各個(gè)環(huán)節(jié)，以最少成本、最短配送距離等滿足客戶的需求。

通過(guò)考慮現(xiàn)實(shí)物流的附加要求和各種路線建設(shè)限制，衍生出許多類(lèi)型的LRP，其中包括帶容量約束的LRP(Capacitated LRP, CLRP)[3-7]、帶同時(shí)取送貨的LRP(LRP with Simultaneous Pickup and Delivery, LRPSPD)[8-11]、帶時(shí)間窗的LRP(LRP with Time Windows, LRPTW)[12-16]、多目標(biāo)LRP(Multi-Objective LRP, MOLRP)[17-22]等。然而，基本LRP以及衍生出不同類(lèi)型的LRP大多是以總運(yùn)行成本為優(yōu)化目標(biāo)，少數(shù)以最小化配送時(shí)間、最大化客戶滿意度為目標(biāo)。但是在全球低碳經(jīng)濟(jì)以及綠色經(jīng)濟(jì)的響應(yīng)下，大量二氧化碳等溫室氣體的排放造成了日趨嚴(yán)重的環(huán)境污染，因此在物流配送過(guò)程中減少二氧化碳?xì)怏w的排放已成為物流公司不得不解決的問(wèn)題?；诖?，本文提出了低碳選址—路徑問(wèn)題(Low Carbon LRP, LCLRP)，并設(shè)計(jì)了一種基于蟻群選擇機(jī)制的進(jìn)化式超啟發(fā)算法(Hyper-Heuristic algorithm of Ant colony Selection mechanism, HHAS)進(jìn)行問(wèn)題模型求解，為物流公司提供參考。據(jù)筆者所知，目前尚未有學(xué)者將超啟發(fā)算法用于求解LCLRP，本文采HHAS算法求解LCLRP，將大洪水、模擬退火、只接受好解作為高層自適應(yīng)接受準(zhǔn)則機(jī)制，以引導(dǎo)全局搜索并加快收斂。

考慮二氧化碳排放的LRP文獻(xiàn)還比較欠缺，但是在VRP中考慮二氧化碳排放以及燃油消耗的文獻(xiàn)比較豐富。車(chē)輛的碳排放量受到很多因素的影響，Demir等[23]指出影響碳排放量的因素主要有車(chē)輛參數(shù)、交通狀況、坡度、環(huán)境以及駕駛?cè)怂郊傲?xí)慣等。車(chē)輛參數(shù)包括車(chē)輛型號(hào)、尺寸大小、重量等等；交通狀況包括道路擁擠情況、紅綠燈等；環(huán)境主要指自然環(huán)境，如濕度、氣溫等。車(chē)輛運(yùn)行過(guò)程中碳排放的計(jì)算十分復(fù)雜，部分文獻(xiàn)將燃油消耗通過(guò)數(shù)學(xué)轉(zhuǎn)換間接獲得。例如，Scott等[24]從地形坡度和車(chē)輛裝載率測(cè)量車(chē)輛的燃油消耗；Zhang等[25]考慮了行駛距離、車(chē)輛裝載量?jī)煞矫鏈y(cè)量車(chē)輛的燃油消耗量；張春苗等[26]將車(chē)輛的重量和運(yùn)輸距離作為影響燃油消耗的主要因素；Palmer[27]建立了貨車(chē)碳排放模型，研究了速度對(duì)燃油消耗的影響。常見(jiàn)模型如下：

(1)基于速度與裝載量的燃油消耗模型

1)Barth等[28]考慮了車(chē)輛的速度與裝載量對(duì)車(chē)輛油耗的關(guān)系，建立了貨車(chē)燃油消耗模型：

Fn=λ(khNhVhd/v+Mhvαd+βhγdv2)。

式中：kh表示發(fā)動(dòng)機(jī)摩擦系數(shù)；Nh表示發(fā)動(dòng)機(jī)轉(zhuǎn)速；Vh表示發(fā)動(dòng)機(jī)排量；d表示行駛距離；v表示車(chē)輛瞬時(shí)速度；Mh表示車(chē)輛總重量；α、βh為車(chē)輛特定常數(shù)；λ、γ為常數(shù)?？梢缘贸觯簁hNhVhd/v為發(fā)動(dòng)機(jī)相關(guān)模塊，與運(yùn)行時(shí)間呈線性數(shù)學(xué)關(guān)系；Mhvαd為重量模塊，與車(chē)輛總重量呈線性關(guān)系；βhγdv2為速度模塊，與速度呈二次方關(guān)系。

2)Zhang等[25]考慮了行駛距離、汽車(chē)行駛速度、車(chē)輛裝載量等，假設(shè)增加車(chē)輛載重M會(huì)增加p百分比的燃油消耗，提出的燃油消耗模型為：

式中：FCij表示從節(jié)點(diǎn)i行駛到節(jié)點(diǎn)j所產(chǎn)生的燃油消耗，LPHij表示車(chē)輛從節(jié)點(diǎn)i到節(jié)點(diǎn)j單位時(shí)間所消耗的燃油消耗，dij表示節(jié)點(diǎn)i到節(jié)點(diǎn)j的行駛距離，vij表示節(jié)點(diǎn)i到節(jié)點(diǎn)j的平均速度。

(2)瞬時(shí)油耗模型

Akcelik等[29]考慮了車(chē)輛自身特性，如車(chē)輛重量、行駛速度、加速度等，建立了車(chē)輛瞬時(shí)燃油消耗模型：

在時(shí)間t內(nèi)的燃油消耗量為

其中：RT是車(chē)輛的牽引力(單位：kN)；M是車(chē)輛裝載量(單位：kg)；v是連續(xù)速度(單位：km/h)；a是加速度(單位：m/s2)；β1是車(chē)輛燃油消耗參數(shù)(單位：g/kJ)；β2是車(chē)輛加速時(shí)相關(guān)數(shù)值g/(kJ·m/s2)。

1 問(wèn)題描述

1.1 碳排放模型

上述兩個(gè)燃油消耗與二氧化碳排放計(jì)算模型考慮了車(chē)輛的自身參數(shù)，建立了與車(chē)輛速度以及加速度有關(guān)的油耗模型，但是在VRP以及LRP中車(chē)輛的速度與加速度很難進(jìn)行具體描述。根據(jù)日本政府網(wǎng)[30]在2010年的報(bào)告，同一類(lèi)車(chē)輛的油耗與載重量成正比例關(guān)系。鑒于此，Xiao等[31]將車(chē)輛考慮為勻速行駛，只考慮載重量與行駛距離對(duì)車(chē)輛燃油排放的影響，將車(chē)輛重量分成空載重量Q0和載重量Q1兩部分，則燃油消耗率(Fuel Consumption Rate, FCR)

φ(Q1)=α(Q1+Q0)+b。

(1)

將車(chē)輛可以承受的最大載重量定義為Q，則空載時(shí)的FCR以及滿載時(shí)的FCR就可以分別表示為φ0=αQ0+b和φ*=α(Q0+Q)+b，其中α=(φ*-φ0)/Q。因此，φ(Q1)=φ0+(α(φ*-φ0)/Q)*Q1，說(shuō)明FCR與車(chē)輛的載重呈現(xiàn)出線性關(guān)系，其中截距為φ0，斜率為φ*-φ0/Q。

車(chē)輛在路徑中的燃油消耗可以表示為[32]：

(2)

式中：xij為決策變量，當(dāng)車(chē)輛由客戶點(diǎn)i駛向客戶點(diǎn)j時(shí)，xij=1，否則xij=0；因?yàn)棣読j與載重量有關(guān)，所以總?cè)加拖腃fuel可以依靠調(diào)整裝貨卸貨順序來(lái)優(yōu)化。令yij為從點(diǎn)i運(yùn)送到點(diǎn)j的貨物重量，因此根據(jù)式(3)可以得到：

(3)

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

LCLRP可定義為一個(gè)完整的定向網(wǎng)絡(luò)G=(N,A)，N由ND與NC構(gòu)成，表示物流系統(tǒng)內(nèi)所有點(diǎn)的集合。其中ND為候選配送中心集合，NC為客戶點(diǎn)集合，每個(gè)候選倉(cāng)庫(kù)的容量與位置已知，每個(gè)客戶的需求量與地理位置已知。車(chē)輛從配送中心出發(fā)，按照安排的路徑依次服務(wù)客戶，滿足客戶的送貨要求，模型的目標(biāo)是選擇配送中心位置并分配車(chē)輛，設(shè)計(jì)車(chē)輛路徑以使車(chē)輛的總碳排放量最小化。各符號(hào)定義如表1所示。

表1 參數(shù)符號(hào)含義

選址—路徑問(wèn)題的數(shù)學(xué)模型可表示如下：

(4)

s.t.

(5)

(6)

(7)

xik≤zik,?i∈NC,k∈ND；

(8)

xki≤zik,?i∈NC,k∈ND；

(9)

?i,j∈NC,i≠j,?k∈ND；

(10)

(11)

(12)

xij∈{0,1},?i,j∈N；

(13)

Yk={0,1},?k∈ND；

(14)

zik={0,1},?i∈NC,k∈ND。

(15)

其中：式(4)為目標(biāo)函數(shù)，表示車(chē)輛碳排放量最??；式(5)～式(6)保證了每個(gè)客戶點(diǎn)都要被訪問(wèn)且只能被訪問(wèn)一次；式(7)表示每個(gè)客戶點(diǎn)只能分配給一個(gè)配送中心；式(8)～式(10)表示所服務(wù)車(chē)輛必須返回原配送中心；式(11)保證每個(gè)配送中心訪問(wèn)的顧客總需求小于配送中心的容量；式(12)保證車(chē)的載重量不超過(guò)載重能力,S為某車(chē)輛所服務(wù)的客戶集合；式(13)～式(15)為決策變量。

2 考慮燃油消耗的選址路徑問(wèn)題優(yōu)化

2.1 超啟發(fā)算法概述

常見(jiàn)的用于求解選址路徑問(wèn)題的算法主要有遺傳算法(Genetic Algorithm, GA)[32]、禁忌搜索(Tabu Search, TS)[33]、粒子群算法(Particle Swarm Optimization, PSO)[34]等，這些算法都是對(duì)解空間進(jìn)行直接操作的。超啟發(fā)式算法(Hyper-Heuristic Algorithms, HHA)是近年來(lái)發(fā)展起來(lái)的一種新型啟發(fā)式算法，可以簡(jiǎn)單闡述為“尋找啟發(fā)式算法的啟發(fā)式算法”，具有更高級(jí)智能系統(tǒng)的特征，是未來(lái)啟發(fā)式優(yōu)化領(lǐng)域研究的主流方向之一，其嚴(yán)格定義如下[35- 36]：超啟發(fā)式算法提供了一種高層次啟發(fā)式方法，通過(guò)管理或操縱一系列底層次啟發(fā)式算法(Low-Level Heuristic, LLH)，以產(chǎn)生新的啟發(fā)式算法。這些新啟發(fā)式算法被用于求解各類(lèi)組合優(yōu)化問(wèn)題。

如圖1所示為超啟發(fā)算法框架，可分為上層控制領(lǐng)域與底層問(wèn)題領(lǐng)域兩個(gè)部分。在控制領(lǐng)域，由智能計(jì)算機(jī)專(zhuān)家進(jìn)行高層啟發(fā)式策略的設(shè)計(jì)，可以分為選擇策略與接受準(zhǔn)則兩個(gè)部分，選擇策略用于搜索底層啟發(fā)式算子構(gòu)成的空間，并根據(jù)底層啟發(fā)式算子的歷史性能選擇底層啟發(fā)式算子；接受準(zhǔn)則根據(jù)子代解的質(zhì)量判斷是否代替父代解。圖1中：Scurrent表示子代解，Sprev表示父代解。在問(wèn)題領(lǐng)域，由領(lǐng)域?qū)＜姨峁┑讓訂l(fā)式算子以及問(wèn)題描述、目標(biāo)函數(shù)等信息。底層啟發(fā)式算子與待解問(wèn)題相關(guān)，用于直接搜索問(wèn)題空間。問(wèn)題領(lǐng)域與控制領(lǐng)域存在領(lǐng)域屏蔽，用于隔離上層控制領(lǐng)域與底層問(wèn)題領(lǐng)域。由于控制領(lǐng)域與問(wèn)題領(lǐng)域之間存在領(lǐng)域屏蔽，只要修改問(wèn)題領(lǐng)域的底層啟發(fā)式算子空間、問(wèn)題描述以及目標(biāo)函數(shù)，就可以將算法應(yīng)用到其他問(wèn)題中。由此可見(jiàn)，超啟發(fā)算法具有很強(qiáng)的通用性。因此，超啟發(fā)算法已經(jīng)應(yīng)用到諸多領(lǐng)域中，如考試時(shí)間表問(wèn)題[37]、圖形著色問(wèn)題[38]、排班問(wèn)題[39]等。

表2給出了不同選擇策略與接受準(zhǔn)則，選擇策略與接受準(zhǔn)則可以?xún)蓛蛇M(jìn)行組合，根據(jù)選擇策略提供的啟發(fā)式機(jī)制來(lái)選擇一個(gè)或一系列的底層啟發(fā)式算子，并由底層啟發(fā)式算子對(duì)父代解進(jìn)行操作獲得子代解，最后由接受準(zhǔn)則所提供機(jī)制判斷是否接受子代解。如SR+AM就表示選擇策略為簡(jiǎn)單隨機(jī)，接受準(zhǔn)則為接受所有解的高層啟發(fā)式策略。

表2 選擇策略和接受準(zhǔn)則

2.2 選擇策略描述

本文設(shè)計(jì)了蟻群選擇機(jī)制作為超啟發(fā)算法的高層進(jìn)化式選擇策略，提出一種基于蟻群選擇機(jī)制的超啟發(fā)算法[40],即將蟻群選擇機(jī)制作為高層啟發(fā)式策略，評(píng)估底層啟發(fā)式算子的性能，并管理和操縱一系列底層啟發(fā)式算子，蟻群選擇機(jī)制并非直接對(duì)問(wèn)題空間進(jìn)行搜索操作，而是運(yùn)用蟻群選擇機(jī)制對(duì)底層算子所組成的空間進(jìn)行搜索操作。每個(gè)頂點(diǎn)代表底層啟發(fā)式算子，如圖2所示，假設(shè)有m只螞蟻，均勻分布在底層啟發(fā)式算子之中，每只螞蟻根據(jù)每個(gè)路徑中的信息素來(lái)確定下一個(gè)頂點(diǎn)，一旦螞蟻到達(dá)一個(gè)新的頂點(diǎn)，便應(yīng)用頂點(diǎn)所對(duì)應(yīng)的底層啟發(fā)式算子。例如，1-1-2-3就表示其中一只螞蟻的路徑，路徑表示底層啟發(fā)式算子的選擇順序；如圖3所示，弧(i,j)表示螞蟻從底層啟發(fā)式算子i移動(dòng)到底層啟發(fā)式算子j。但是不同于TSP，在超啟發(fā)算法中，算法對(duì)底層啟發(fā)式算子進(jìn)行操作，因此每只螞蟻允許多次訪問(wèn)同一個(gè)點(diǎn)，意味著每個(gè)底層啟發(fā)式算子在一次迭代中可以進(jìn)行多次操作。

參考文獻(xiàn)[40]，采用可見(jiàn)度函數(shù)ηj來(lái)評(píng)判底層啟發(fā)式算子的優(yōu)化性能:

(19)

式中：Ikj(t)是在螞蟻k的路徑中，底層啟發(fā)式算子j在迭代次數(shù)為t時(shí)的改進(jìn)值(可以為負(fù))；Tkj(t)是在螞蟻k的路徑中，底層啟發(fā)式算子j在迭代次數(shù)為t時(shí)所占用的CPU運(yùn)行時(shí)間；γ是一個(gè)位于0到1之間的權(quán)重值。

每條弧上的信息素

(20)

式中：ρ(0<ρ<1)表示信息素的蒸發(fā)系數(shù)，1-ρ表示信息素的持久系數(shù)；Pk(t)表示螞蟻k周游過(guò)程所產(chǎn)生的選擇底層啟發(fā)式算子的順序；#ij(Pk(t))表示在螞蟻k的路徑中弧(i，j)所出現(xiàn)的次數(shù)；I(Pk(t))表示螞蟻k在路徑中改進(jìn)的量；T(Pk(t))表示螞蟻k完成路徑周游所占用的CPU時(shí)間。

蟻群選擇的決策過(guò)程需要算法結(jié)合信息素和可見(jiàn)度兩個(gè)量來(lái)進(jìn)行，因此，引入一個(gè)變量V來(lái)結(jié)合信息素和能見(jiàn)度：

Vij(t)=αηj(t)+βτij(t)。

(21)

這里，為了使表現(xiàn)不好的底層啟發(fā)式算子也有一定的概率被選中，引入變量PV：

PVij(t)=max{Vij(t),QσVij(t)}。

(22)

式中

(23)

式中：H為底層啟發(fā)式算子的集合;n=|H|表示底層啟發(fā)式算子的個(gè)數(shù)；ε和σ保證表現(xiàn)性能不佳的底層啟發(fā)式算子仍有一個(gè)很小但不為0的概率被選中。若弧(i,j)未出現(xiàn)在路徑中，則PVij(t)=0；若所有底層啟發(fā)式算子的性能表現(xiàn)都不好，則Q=0，并將所有PV的值都設(shè)為1，保證所有底層啟發(fā)式算子都以相等的概率被選取。最后，弧(i,j)被選擇的概率

(24)

2.3 接受準(zhǔn)則描述

接受準(zhǔn)則是超啟發(fā)算法的上層框架，用于判斷子代解是否能替換父代解，其表現(xiàn)性能的好壞直接影響超啟發(fā)算法的收斂速度與優(yōu)化精度。

(1)接受所有解 無(wú)論改進(jìn)與否，接受所有解。

(2)概率接受 接受所有的改進(jìn)解，對(duì)于非改進(jìn)解，有50%的概率接受。

(3)模擬退火 與模擬退火算法類(lèi)似。設(shè)置初始溫度T、降溫系數(shù)λ，每次迭代后根據(jù)溫度更新退火溫度：T=T×λ。每次迭代使用底層啟發(fā)式算子，若得到了改進(jìn)解，則接受；若得到了非改進(jìn)解，則以一定的概率接受，概率p=eΔ/T。

(4)大洪水 設(shè)置下界LB，LB為當(dāng)前迭代的最優(yōu)解。每次迭代使用底層啟發(fā)式算子，若得到了改進(jìn)解，且子代解ftemp小于一閾值，則接受該子代解，該閾值為L(zhǎng)B+(f0-LB)×(1-iter/Tmax)。其中：ftemp表示子代解，f0表示初始解，iter表示當(dāng)前迭代次數(shù)，Tmax表示最大迭代次數(shù)。

2.4 算法流程

HHAS算法用于求解LCLRP模型時(shí)，首先進(jìn)行配送中心選址分配，再采用路徑優(yōu)化的方法進(jìn)行求解，具體流程如下：

步驟1初始化。設(shè)螞蟻運(yùn)行路徑長(zhǎng)度為L(zhǎng)P。首先進(jìn)行種群初始化，求出初始種群的解S，并將S賦值于最佳解Sb。對(duì)于底層啟發(fā)式算子j∈H，初始化可見(jiàn)度ηj=0。對(duì)于路徑中每條弧(i,j),設(shè)置一個(gè)初始值τij(t)=0。將m只螞蟻均勻分布在n個(gè)頂點(diǎn)上，并賦予每只螞蟻初始種群與初始解S，并將螞蟻k的解賦值于矩陣Sk中。

步驟3蟻群選擇機(jī)制參數(shù)更新。根據(jù)式(19)更新底層啟發(fā)式算子的可見(jiàn)度ηj，根據(jù)式(20)更新各個(gè)底層啟發(fā)式算子之間的信息素τij，最后根據(jù)式(24)更新選擇底層啟發(fā)式算子概率probabilityijk(t)。

步驟4滿足停止條件(即達(dá)到迭代次數(shù))，則輸出結(jié)果，否則返回步驟2。

假設(shè)采用蟻群選擇機(jī)制作為選擇策略，模擬退火作為接受準(zhǔn)則，超啟發(fā)算法組合的偽代碼如算法1所示。

算法1超啟發(fā)算法偽代碼。

1:Input:Low-level heuristics(LLH1-LLH13)

2: Initial Solution:S

3: Number of iterations:C

4: Pheromones:τ

5: Evaporation coefficient:ρ

6: Initial temperature:T

7: Minimum temperature:Tmin

8: Cooling coefficient:r

9:Output:Best Solution:S′

10:Solution S←initialiseSolution()

11:for i←1:C

12: Solution S′←S

13: LLHi←SelectHeuristics(AntSelect(LLHs))

14: S′←ApplyHeuristic(LLHi)

15: Δ←Function(S)-Function(S′)

16: while T>Tmindo

17: if Δ>0 then

18: S=S′

19: else

20: if exp(Δ/T)>r then

21: S←S′

22: end if

23: T=T*r

24: end if

25: end while

26: τ←(1-ρ)×τ+α

27: Update probability

28:end for

2.5 底層啟發(fā)式算子描述

在本文研究的LCLRP模型中，底層啟發(fā)式算子根據(jù)文獻(xiàn)[41]的分類(lèi)標(biāo)準(zhǔn)，構(gòu)建以下算子：變異算子(mutation heuristics)、破壞重組算子(ruin-recreate heuristics)、局部搜索算子(local search heuristics)、交叉算子(crossover heuristics)。底層啟發(fā)式算子用于直接對(duì)解空間進(jìn)行操作，具體如下：

(1)變異算子 變異算子是一種對(duì)解微小擾動(dòng)的突變算子。所有這些算子都將返回任何滿足約束的路徑，無(wú)論目標(biāo)函數(shù)值是改進(jìn)或惡化。

1)路徑變異算子。

LLH1：2-opt。如圖4所示，任意選擇一條路徑，交換相鄰兩客戶點(diǎn)的位置；

LLH2：Or-opt。如圖5所示，任意選擇一條路徑，把相鄰兩客戶點(diǎn)插入到其他位置；

LLH3：Interchange。如圖6所示，任意選擇兩條路徑，交換任意兩個(gè)客戶點(diǎn)位置；

LLH4：Shift。如圖7所示，任意選擇兩條路徑，在一條路徑上找一個(gè)任意的客戶插入到另一條路徑中合適的位置。

2)配送中心變異算子。

LLH5：Interchange。如圖8所示，任意選擇兩條路徑，交換配送中心；

LLH6：Shift。如圖9所示，任意選擇一條路徑，開(kāi)啟一個(gè)新的配送中心，并替換這條路徑的配送中心。

(2)局部搜索算子 與變異算子一樣，通常是以交換或插入的形式產(chǎn)生新路徑。然而，與變異算子不同的是局部搜索算子只會(huì)返回一個(gè)等于或優(yōu)于初始的解。

1)LLH7：2-opt*。任意選擇兩條路徑，根據(jù)衡量準(zhǔn)則選擇一個(gè)位置，交換此位置后的所有客戶點(diǎn)，并且只接受好解。

2)LLH8：Shift。和LLH4一樣。不同的是LLH8根據(jù)衡量準(zhǔn)則來(lái)選擇客戶，備選客戶插入到另一路徑的任意位置，并且只接受好解。

3)LLH9：Interchange。和LLH3一樣，不同的是此算子只接受好解。

4)LLH10：GENI。計(jì)算不同路徑上的任意兩個(gè)客戶的距離，采用最短距離作為基準(zhǔn)距離，選擇移除距離最接近基準(zhǔn)距離的客戶，并且只接受好解。

(3)破壞重組算子 選擇一些客戶根據(jù)其與“基準(zhǔn)客戶”的接近程度從路徑中刪除。

LLH11:Location-based Radial Ruin。隨機(jī)選擇一個(gè)客戶作為“基準(zhǔn)客戶”，根據(jù)客戶位置接近“基準(zhǔn)客戶”的原則，以[1%，10%]的概率將客戶從路徑中刪除。

(4)交叉算子 選擇兩個(gè)父代路徑作為輸入，生成一條子代路徑，通?？梢酝ㄟ^(guò)一點(diǎn)、兩點(diǎn)和均勻交叉來(lái)完成。

1)LLH12：Combine。隨機(jī)選擇兩個(gè)父代，以[25%,75%]的概率將一個(gè)父代復(fù)制得到一條子路徑，添加來(lái)自另一個(gè)父代中非沖突的路徑，并隨機(jī)插入剩余客戶。

2)LLH13：Longest Combine。隨機(jī)選擇兩個(gè)父代，所有路徑以服務(wù)客戶點(diǎn)數(shù)量的降序來(lái)考慮，任何不重復(fù)客戶的路徑都會(huì)被添加生成一條子路徑，并隨機(jī)插入剩余客戶。

3 實(shí)驗(yàn)結(jié)果及分析

3.1 算法參數(shù)設(shè)定正交實(shí)驗(yàn)

關(guān)于燃油消耗的參數(shù)，車(chē)輛空載時(shí)單位距離燃油消耗系數(shù)φ0=1，車(chē)輛滿載時(shí)單位距離燃油消耗系數(shù)φ*=2[7]，燃油轉(zhuǎn)換系數(shù)C0=2.68。不考慮配送中心和車(chē)輛的固定成本。

HHAS算法參數(shù)設(shè)定如下：螞蟻路徑長(zhǎng)度LP=11；螞蟻數(shù)量m=(10，15，20)；權(quán)重參數(shù)α=(0.6，0.7，0.8),β=(0.6，0.7，0.8),γ=(0.6，0.7，0.8)；增強(qiáng)系數(shù)ε=0.001,σ=1.001[41]。

由于本節(jié)尚未確定接受準(zhǔn)則的性能，每次迭代隨機(jī)選擇一種接受準(zhǔn)則(GD、SA、OI)。

對(duì)于HHAS算法，單個(gè)螞蟻在一次循環(huán)中所經(jīng)過(guò)的路徑，即為底層啟發(fā)式算子集中的一個(gè)解，m只螞蟻在一次循環(huán)中所經(jīng)過(guò)的路徑，則表現(xiàn)為底層啟發(fā)式算子的一個(gè)子集，螞蟻數(shù)量的變化決定了算法的穩(wěn)定性；而α、β以及γ反映了蟻群在路徑搜索中隨機(jī)因素作用的強(qiáng)度。因此，為了使所提算法獲得更好的算法性能，針對(duì)關(guān)鍵參數(shù)螞蟻數(shù)量m、權(quán)重參數(shù)α、β、γ，設(shè)置3個(gè)參數(shù)水平(如表4)，選擇規(guī)模L(34)的正交實(shí)驗(yàn)，在各種參數(shù)組合下選擇中等規(guī)模的基準(zhǔn)實(shí)例Christ50x5ba進(jìn)行10次獨(dú)立實(shí)驗(yàn)，每次實(shí)驗(yàn)迭代1 000次，測(cè)試結(jié)果如表5所示。HHAS算法使用MATLAB 2015編程，計(jì)算機(jī)環(huán)境為Intel(R)Core(TM)i5-3230M CPU @2.60 GHz, 4 GB RAM, Windows 10系統(tǒng)。

表4 參數(shù)水平分布

表5 不同參數(shù)正交實(shí)驗(yàn)結(jié)果

由表5，可得如表6所示各參數(shù)響應(yīng)值。由表5與表6可知，m、α、β、γ均會(huì)對(duì)算法性能造成影響：螞蟻數(shù)量m增多，可以提高算法的全局搜索能力以及算法的穩(wěn)定性，但m增大后，會(huì)使大量的曾被搜索過(guò)的路徑上的信息素變化趨于平均，信息素正反饋?zhàn)饔貌幻黠@，螞蟻數(shù)量減小，搜索的隨機(jī)性減弱，使算法的全局性能降低，容易出現(xiàn)過(guò)早停滯；α、β、γ的大小反映了蟻群在路徑搜索中隨機(jī)性因素作用的強(qiáng)度，較大的值將造成隨機(jī)性減弱，而當(dāng)值過(guò)小時(shí)，則易使蟻群的搜索陷入局部最優(yōu)。如圖10所示折線圖，在當(dāng)m、α、β、γ處于參數(shù)水平2時(shí)所得到的解最優(yōu)，即m=15,α=0.7,β=0.7,γ=0.7。

表6 各參數(shù)響應(yīng)值

3.2 不同接受準(zhǔn)則實(shí)驗(yàn)結(jié)果對(duì)比

對(duì)于超啟發(fā)算法，不同的高層啟發(fā)式策略的性能效果是不同的，為了使HHAS算法達(dá)到更好的性能，通過(guò)實(shí)驗(yàn)對(duì)比來(lái)選擇出一個(gè)最佳的接受準(zhǔn)則：選擇策略均為蟻群選擇機(jī)制，接受準(zhǔn)則分別為GD、SA以及OI。每個(gè)基準(zhǔn)實(shí)例獨(dú)立運(yùn)行10次，迭代次數(shù)為1 000次。實(shí)驗(yàn)結(jié)果如表7所示，其中HHAS算法中的參數(shù)m=15,α=0.7,β=0.7,γ=0.7。

由表7可以看出，在18個(gè)基準(zhǔn)實(shí)例測(cè)試中，AS+OI取得最優(yōu)解的個(gè)數(shù)為13次，AS+SA取得最優(yōu)解5次，而AS+GD只獲得了0次最優(yōu)解，其中，在小規(guī)模實(shí)例中，AS+OI獲得最小解6次，AS+SA獲得最小解2次；在中等規(guī)模實(shí)例中，AS+OI獲得最小解6次，AS+SA獲得最小解2次；在大規(guī)模實(shí)例中，AS+OI獲得最小解1次，AS+SA獲得最小解1次，AS+GD獲得最小解0次。而且，在算法穩(wěn)定性上，OI接受準(zhǔn)則所得結(jié)果標(biāo)準(zhǔn)差比GD、SA接受準(zhǔn)則所得到的結(jié)果標(biāo)準(zhǔn)差小，算法穩(wěn)定性更高。

表7 不同接受準(zhǔn)則計(jì)算碳排放量實(shí)驗(yàn)結(jié)果

例如在圖11，Gaskell22x5的求解結(jié)果中：GD接受準(zhǔn)則得到的最佳解為2 685.0 kg，平均值2 688.1 kg，標(biāo)準(zhǔn)差為0.9；SA接受準(zhǔn)則得到的最佳解為2 407.7 kg，平均值為2 426.0kg，標(biāo)準(zhǔn)差為8.7；OI接受準(zhǔn)則得到的最佳解為2 288.5 kg,平均值為2 288.9 kg，標(biāo)準(zhǔn)差為0.2；Gaskell36x5的計(jì)算結(jié)果中：GD接受準(zhǔn)則得到的最佳解為3 385.0 kg，平均值為3 393.3 kg，標(biāo)準(zhǔn)差為1.9；SA接受準(zhǔn)則得到的最佳解為3 214.0 kg，平均值為3 243.1 kg，標(biāo)準(zhǔn)差為8.8；OI接受準(zhǔn)則得到的最佳解為2 576.5 kg，平均值為2 580.5 kg，標(biāo)準(zhǔn)差為2.4；Perl85x7的計(jì)算結(jié)果中，GD接受準(zhǔn)則獲得的最佳解為6 831.6 kg，平均值為6 927.5 kg，標(biāo)準(zhǔn)差為25.6；SA接受準(zhǔn)則獲得的最佳解為6 692.8 kg，平均值為6 750.4 kg，標(biāo)準(zhǔn)差為12.9；OI接受準(zhǔn)則獲得的最佳解為6 287.6 kg，平均值為6 310.2 kg，標(biāo)準(zhǔn)差為4.8。綜上所述，AS+OI在求解不同規(guī)模的實(shí)例時(shí)表現(xiàn)均優(yōu)于AS+GD與AS+SA，接受準(zhǔn)則為OI時(shí)超啟發(fā)算法性能最優(yōu)，在接下來(lái)的實(shí)驗(yàn)分析中將使用OI接受準(zhǔn)則進(jìn)行實(shí)驗(yàn)。

3.3 算法求解基本LRP有效性分析

為說(shuō)明HHAS算法對(duì)求解LRP的有效性，用HHAS算法求解基本LRP模型，并與LRGTS算法[5]、GRASP算法[42]的計(jì)算結(jié)果進(jìn)行對(duì)比。HHAS算法運(yùn)行環(huán)境與3.1節(jié)中的一致，GRASP算法與LRGTS算法均使用Visual C++編碼，運(yùn)行環(huán)境為DELL PC Optiplex GX260,2.4 GHz Pentium 4處理器，512 MB運(yùn)行內(nèi)存，Windows XP操作系統(tǒng)。HHAS算法求解Barreto算例的流程與2.4節(jié)中的算法流程一致，蟻群選擇機(jī)制參數(shù)為m=15,α=0.7，β=0.7，γ=0.7，接受準(zhǔn)則為OI，不同之處在于優(yōu)化的目標(biāo)函數(shù)為總成本而不是車(chē)輛碳排放?？偝杀綶42]：

式中：Oi表示配送中心開(kāi)放成本;ci表示配送過(guò)程中的路徑成本;F表示車(chē)輛成本;xijk為決策變量，當(dāng)車(chē)輛k由節(jié)點(diǎn)i行駛向節(jié)點(diǎn)j時(shí)xijk=1，否則xijk=0。其他約束條件與1.2節(jié)中的一致。得到的仿真結(jié)果計(jì)算比較結(jié)果見(jiàn)表8。其中：Q表示車(chē)輛的最大載重量，BKS表示已知最優(yōu)解，cost表示成本，gap表示成本相對(duì)LB的偏差率，cpu表示程序運(yùn)行時(shí)間，HHAS算法參數(shù)設(shè)置為m=15，α=0.7，β=0.7，γ=0.7。每個(gè)實(shí)例獨(dú)立運(yùn)行10次，迭代次數(shù)為5×(M+N+K)2,其中：M表示客戶數(shù)量，N表示配送中心數(shù)量，K表示車(chē)輛數(shù)量。

表8 算法求解LRP問(wèn)題實(shí)驗(yàn)結(jié)果

續(xù)表8

如表8所示，HHAS算法相對(duì)LRGTS算法與GRASP算法時(shí)間效率較差，但在3個(gè)算法對(duì)比過(guò)程中，HHAS算法在13個(gè)算例中總共獲得12次最優(yōu)，優(yōu)化效果明顯優(yōu)于LRGTS算法與GRASP算法，在有效時(shí)間內(nèi)，可以獲得更好的解，說(shuō)明HHAS算法對(duì)于求解LRP是有效的。

3.4 算法求解LCLRP有效性分析

為了進(jìn)一步說(shuō)明HHAS算法對(duì)求解低碳LRP的有效性，將求解結(jié)果與量子進(jìn)化算法(Quantum Evolutionary Algorithm, QEA)[44]對(duì)比。對(duì)比實(shí)驗(yàn)所使用的實(shí)例采用LRP基準(zhǔn)實(shí)例的Prodhon[5]算例，目標(biāo)函數(shù)為碳排放量，結(jié)果用CE表示，cost表示成本。QEA的實(shí)驗(yàn)數(shù)據(jù)從文獻(xiàn)[26]中獲得，算法使用MATLAB編程，計(jì)算機(jī)環(huán)境為intel(R)Xeon(R)CPU E3-1230 v5 @3.40 GHz,4 GB RAM, Windows 10系統(tǒng)。對(duì)比實(shí)驗(yàn)結(jié)果如表9和圖12所示。

表9 求解CECLRP對(duì)比實(shí)驗(yàn)結(jié)果

如表9所示，HHAS算法求解得到的碳排放量均優(yōu)于QEA算法，14個(gè)算例求解中最大的改進(jìn)率可以達(dá)到14.3%，平均改進(jìn)率為5.7%?；诖?，該實(shí)驗(yàn)驗(yàn)證了HHAS算法對(duì)于求解LCLRP數(shù)學(xué)模型的有效性。

3.5 LCLRP與CLRP對(duì)比實(shí)驗(yàn)分析

為證明所提LCLRP模型的有效性，與以車(chē)輛總行駛距離為目標(biāo)的CLRP模型所得到的數(shù)據(jù)進(jìn)行對(duì)比。本文選擇了算法參數(shù)為m=15,α=0.7,β=0.7,γ=0.7，接受準(zhǔn)則為OI的HHAS算法對(duì)于求解LCLRP和CLRP實(shí)例的計(jì)算結(jié)果。CLRP模型以車(chē)輛的總成本為目標(biāo)。每個(gè)基準(zhǔn)實(shí)例都運(yùn)行10次，每次迭代次數(shù)為5×(M+N+K)2，計(jì)算得出車(chē)輛行駛的總成本(Cost)與總碳排放量(CE)。LCLRP與CLRP運(yùn)行10次所取得的最小值(min)與平均值(mean)，結(jié)果如表10所示。其中，RelatedCost是在LCLRP以碳排放量為目標(biāo)所得到的相關(guān)成本，RelatedCE是在CLRP以總成本為目標(biāo)所得到的相關(guān)碳排放量。表格最后兩列描述的是差值百分比，計(jì)算公式如下：

表10 CECLRP與CLRP實(shí)驗(yàn)對(duì)比結(jié)果

CEdev.=((CE-RelatedCE)/

RelatedCE)×100，

(25)

Costdev.=((RelatedCost-Cost)/Cost)×100。

(26)

由表10可以看出，求解LCLRP模型得到的平均碳排放量比求解CLRP模型得到的平均碳排放量低5.7%，而總成本增加了4.0%；在比較運(yùn)行10代所獲得的最佳解時(shí)，求解LCLRP模型得到的碳排放量比求解CLRP模型得到的碳排放量低5.7%，距離卻只增加了3.7%。由圖13可以看出，減少的碳排放量明顯低于增加的總成本，其中基準(zhǔn)實(shí)例Gaskell22x5、Gaskell32x5a、Gaskell32x5b以及Christ100x10的優(yōu)化效果最為明顯。以Christ100x10為例，比較以碳排放量最少為優(yōu)化目標(biāo)與以總成本最小為優(yōu)化目標(biāo)的路徑，結(jié)果如圖14所示。

如圖14，當(dāng)以總成本最小為優(yōu)化目標(biāo)時(shí)，求解得到的碳排放量為2 819.3 kg，對(duì)應(yīng)的總成本為836.4元，而當(dāng)以碳排放量最少為優(yōu)化目標(biāo)時(shí)，求解得到的碳排放量值為2 516.9 kg，對(duì)應(yīng)的總成本為863.4元，碳排放節(jié)約了約10.7%，而總成本只增加了約3.1%，優(yōu)化效果明顯。由圖中可以看出，當(dāng)以碳排放量最少為優(yōu)化目標(biāo)得出的車(chē)輛行駛距離大于以總成本最小為優(yōu)化目標(biāo)的車(chē)輛行駛距離時(shí)，碳排放量反而減少了。因此，簡(jiǎn)單的行駛路徑減少并不一定意味著碳排放量減少，合理安排車(chē)輛路徑，優(yōu)化卸貨順序即減少車(chē)輛的載重量可以有效減少物流過(guò)程中的碳排放量。

4 結(jié)束語(yǔ)

本文考慮了環(huán)境污染這一現(xiàn)狀，將碳排放因素考慮到車(chē)輛路徑中，建立了以碳排放最少為優(yōu)化目標(biāo)的LCLRP模型，并提出一種基于蟻群選擇機(jī)制的超啟發(fā)算法，構(gòu)造了基于LCLRP模型的底層啟發(fā)式算子。通過(guò)對(duì)基準(zhǔn)實(shí)例的計(jì)算測(cè)試出算法最佳的參數(shù)組合以及高層策略組合；通過(guò)對(duì)比以碳排放量最少為優(yōu)化目標(biāo)函數(shù)的LCLRP與以總成本最小為優(yōu)化目標(biāo)函數(shù)的CLRP發(fā)現(xiàn)，以碳排放量最少為優(yōu)化目標(biāo)的LCLRP模型可以有效減少車(chē)輛行駛過(guò)程中的碳排放量，且總成本增加量少，模型的有效性得到了驗(yàn)證，可為物流公司提供一定的參考。

LRP以及VRP對(duì)于速度以及加速度的改變難以進(jìn)行預(yù)測(cè)，因此本文構(gòu)建的碳排放量模型僅考慮了車(chē)輛的裝載以及行駛距離。在未來(lái)的研究過(guò)程中，可將其他因素如路況、車(chē)輛速度等因素參數(shù)加入到碳排放模型中；另外，本文所設(shè)計(jì)的基于蟻群選擇機(jī)制的超啟發(fā)式算法效率相對(duì)已有算法的效率表現(xiàn)較差，如何提高算法效率也是下階段研究的問(wèn)題之一。