張 坤,張惠珍,馬 良,張 博

(上海理工大學(xué) 管理學(xué)院,上海 200093)

選址路徑問題(Location-Routing Problem,LRP)是一類經(jīng)典的組合優(yōu)化問題,需要同時(shí)確定候選設(shè)施建立的位置以及車輛行駛的路徑。隨著Cooper[1]首次提出LRP,將設(shè)施選址問題(Facility Location Problem,FLP)與車輛路徑問題(Vehicle Routing Problem,VRP)結(jié)合起來后,學(xué)界與業(yè)內(nèi)逐漸發(fā)現(xiàn)選址與路徑?jīng)Q策的耦合性很高,因此,LRP也成為了供應(yīng)鏈、物流以及運(yùn)籌優(yōu)化領(lǐng)域的熱點(diǎn)問題。截至目前,將標(biāo)準(zhǔn)LRP與實(shí)際應(yīng)用相結(jié)合,多種擴(kuò)展LRP 已被提出,如有容量限制的LRP(Capacitated LRP,CLRP)[2]、帶時(shí)間窗的CLRP(CLRP with Time Windows,CLRPTW)[3]和多車型 CLRP (CLRP with Heterogeneous Fleet,CLRPHF)[4]等。

近年來,隨著生態(tài)環(huán)保意識(shí)的普遍提升,低碳物流已成為運(yùn)籌優(yōu)化和物流領(lǐng)域的另一個(gè)研究熱點(diǎn),通過在物流優(yōu)化問題(LRP,VRP等)中加入低碳物流要素,旨在提高物流環(huán)節(jié)效率,減少碳排放量。為此,國(guó)內(nèi)外學(xué)者深入研究了車輛行駛過程中的碳排放計(jì)算方法。Score等[5]和Barth等[6]基于具體車型的詳細(xì)參數(shù)對(duì)車輛的瞬時(shí)油耗進(jìn)行估計(jì)。Charis等[7]提出了一種數(shù)據(jù)驅(qū)動(dòng)模型,通過對(duì)發(fā)動(dòng)機(jī)類型以及車輛屬性等一系列歷史數(shù)據(jù)估算車輛在行駛中的碳排放。在提高估算精度的同時(shí),許多學(xué)者亦對(duì)估算方法做了適當(dāng)?shù)暮?jiǎn)化以適應(yīng)優(yōu)化計(jì)算。Zhang等[8]設(shè)計(jì)了一種考慮行駛距離、車輛平均速度和車輛負(fù)載等參數(shù)的平均消耗模型以解決低碳背景下的VRP。Xiao等[9]在VRP 中考慮了車輛負(fù)載對(duì)油耗、碳排放的影響,其統(tǒng)計(jì)結(jié)果顯示油耗量與車輛負(fù)載線性相關(guān)。張春苗等[10]和王舜等[11]在LRP中加入碳排放計(jì)算方法,提出了LCLRP(Low-Carbon LRP),并討論了物流網(wǎng)絡(luò)中運(yùn)輸成本與碳排放的關(guān)系。

此外,宏觀層面的政策研究也帶來了一定啟示。已有研究表明,合理的碳交易定價(jià)機(jī)制和碳足跡測(cè)算方法是實(shí)現(xiàn)節(jié)能減排的關(guān)鍵之一[12-13]。碳定價(jià)機(jī)制是當(dāng)前國(guó)際上用于減少碳排放的主要政策工具,其以市場(chǎng)化手段解決了碳排放負(fù)外部性的內(nèi)生化問題[14],有望成為中國(guó)溫室氣體減排的主要政策[15]。因此,基于碳定價(jià)背景下的物流優(yōu)化,可以準(zhǔn)確計(jì)算運(yùn)輸過程中的碳排放,同時(shí)在決策層面統(tǒng)一碳排放量與物流成本的量綱,降低企業(yè)的決策難度。為此,本文在CLRP模型的基礎(chǔ)上加入油耗計(jì)算公式[9],并在碳定價(jià)背景下構(gòu)建了計(jì)算總成本的方法,建立了碳定價(jià)背景下的LCLRP 模型,該模型更貼近現(xiàn)實(shí)情況,能為決策者提供較好的決策依據(jù)和支持。

基本LRP 及其擴(kuò)展問題均屬于NP-Hard 問題。對(duì)于大規(guī)模LRP,精確算法往往無法在可以接受的時(shí)間內(nèi)對(duì)其求解。因此,在有限時(shí)間內(nèi)可以求得滿意解的啟發(fā)式算法成為了求解該問題的主要方法。Wang等[16]和胡大偉等[17]在求解LCLRP時(shí)加入了變鄰域搜索算法(Variable Neighborhood Search,VNS)。冷龍龍[18]提出了一種新的超啟發(fā)式方案求解異構(gòu)車隊(duì)的同時(shí)送取貨選址路徑(LRP with Simultaneously Pickup and Delivery。LRPSPD)。Cao等[19]利用文化基因算法求解了低碳背景下的電動(dòng)車路徑問題。這些方法在實(shí)際求解中均取得了不錯(cuò)的效果。

由無免費(fèi)午餐定理可知,沒有任何一種算法對(duì)所有同類問題均能表現(xiàn)出良好的性能,因此將不同優(yōu)化算法的搜索機(jī)制及策略相融合,構(gòu)建求解性能更高的混合算法已成為目前智能優(yōu)化領(lǐng)域的一大研究方向。灰狼算法(Grey Wolf Optimizer,GWO)是一種模擬狼群狩獵與圍攻行為來搜索最優(yōu)解的新興智能算法[20],其能夠較好地跳出局部最優(yōu),具有較強(qiáng)的全局搜索能力且所需調(diào)整的參數(shù)少。分布估計(jì)算法(Estimation of Distribution Algorithms,EDA)是一類基于種群迭代產(chǎn)生的新興啟發(fā)式算法,通過學(xué)習(xí)優(yōu)良解的概率分布以實(shí)現(xiàn)種群進(jìn)化,其全局搜索能力強(qiáng)且易于實(shí)現(xiàn),并在生產(chǎn)調(diào)度領(lǐng)域表現(xiàn)出了良好的效果[21-23]。為了擴(kuò)展GWO 及EDA的應(yīng)用領(lǐng)域,并為L(zhǎng)RP及其擴(kuò)展問題提供一種有效的求解方法,本文將GWO 及EDA 相結(jié)合,使其相互取長(zhǎng)補(bǔ)短,設(shè)計(jì)了一種分布估計(jì)灰狼算法(Grey Wolf Optimizer with Estimation of Distribution Algorithms,GWOEDA),并將其用于所構(gòu)建的LCLRP模型的求解。

1 基于碳定價(jià)機(jī)制的低碳選址-路徑模型(LCLRP)

本文考慮的LCLRP可被描述為:給定M個(gè)候選設(shè)施點(diǎn)和N個(gè)需求點(diǎn)(任意需求點(diǎn)間均聯(lián)通),從M個(gè)候選設(shè)施點(diǎn)中選擇開放若干設(shè)施并規(guī)劃配送路徑,以完成對(duì)N個(gè)需求點(diǎn)的配送任務(wù),并使總成本(包括設(shè)施啟用成本、運(yùn)輸成本和碳排放成本)最小。

1.1 碳排放成本

構(gòu)建優(yōu)化模型時(shí),由于車輛詳細(xì)參數(shù)以及行駛過程的瞬時(shí)速度難以預(yù)測(cè),故通常將車輛行駛時(shí)的速度設(shè)為勻速[9-11],在模型中僅考慮車輛的負(fù)載與行駛距離。車輛行駛過程中的油耗與負(fù)載呈線性關(guān)系,可近似表示為[9]

式中:ρ(h)為車輛載貨時(shí)單位行駛距離的油耗量,h為其負(fù)載量;ρ*與ρ0分別為滿載與空載時(shí)的油耗量;m為車輛最大容量。

以往研究中,減少油耗量通常被視作減碳的等價(jià)手段[10-11],在此基礎(chǔ)上,本文模型合理地計(jì)算了總運(yùn)輸成本,并添加了碳稅項(xiàng)以表現(xiàn)決策方在碳定價(jià)背景下的碳排放成本。載重為h的車輛從點(diǎn)i行駛至點(diǎn)j的碳排放成本為

式中,碳排放量由單位油耗量ρ(h)、i與j之間的距離oij和燃油轉(zhuǎn)換系數(shù)η相乘得到。碳排放量與碳排放成本(碳稅或碳交易價(jià)格)可以通過固定價(jià)格φ進(jìn)行轉(zhuǎn)化[13]。

1.2 模型假設(shè)與符號(hào)說明

為便于構(gòu)建模型,給出如下假設(shè):

(1) 候選設(shè)施點(diǎn)的位置、最大容量以及建設(shè)成本已知。

(2) 需求點(diǎn)的位置和需求量已知。

(3) 車輛單位負(fù)載單位行駛距離的運(yùn)輸成本、車輛最大負(fù)載以及啟用成本已知。

(4) 車輛必須由某設(shè)施點(diǎn)出發(fā)并回到該設(shè)施點(diǎn)。

(5) 設(shè)施點(diǎn)間無車輛往來。

(6) 每輛車只被啟用一次。

模型所用符號(hào)說明:

集合

I——需求點(diǎn)集合,I={1,2,…,n}

J——候選設(shè)施點(diǎn)集合,J={n+1,n+2,…,n+m}

G——G=I∪J

K——車輛

集合

參數(shù)

di——需求點(diǎn)i的需求量,i∈I

ej——設(shè)施j的開放成本,j∈J

cj——設(shè)施j的最大容量,j∈J

sk——車輛k的啟用成本,k∈K

mk——車輛k的最大容量,k∈K

rk——車輛k的燃油價(jià)格,k∈K

oij——點(diǎn)i至點(diǎn)j的距離,i,j∈G

φ——交易市場(chǎng)中單位指標(biāo)碳排放價(jià)格

ηk——車輛k的油耗轉(zhuǎn)換碳排放 量系數(shù),k∈K

——車輛k滿載時(shí)單位行駛距離的油耗量,k∈K

決策變量

xijk——0-1變量(若車輛k從節(jié)點(diǎn)i訪問j則為1,否則為0),?i,j∈G,k∈K

yjk——0-1變量(若車輛k分配給設(shè)施點(diǎn)j則為1,否則為0),?j∈J,k∈K

zj——0-1變量(若設(shè)施點(diǎn)j啟用則為1,否則為0),?j∈J

hjk——實(shí)數(shù)變量,表示車輛k離開節(jié)點(diǎn)j時(shí)的負(fù)載,?j∈G,k∈K

1.3 LCLRP的非線性混合整數(shù)規(guī)劃模型

在上述假設(shè)條件與參數(shù)設(shè)定的基礎(chǔ)上,本文所構(gòu)建的碳定價(jià)背景下的LCLRP模型描述如下:

目標(biāo)函數(shù)式(3)最小化總成本,包括:設(shè)施開放成本、車輛啟用成本、油耗成本和碳排放成本。式(4)為碳定價(jià)情景下的碳排放計(jì)算公式。約束條件式(5)～(7)表示只有開放的設(shè)施才能服務(wù)需求點(diǎn),且每個(gè)需求點(diǎn)必須有且僅有一輛車為其提供服務(wù);式(8)表示車輛駛?cè)胄枨簏c(diǎn)后也必須駛出;式(9)表示每輛車至多被啟用一次;式(10)表示每個(gè)需求點(diǎn)必須被服務(wù)一次;式(11)表示設(shè)施點(diǎn)間無車輛往來;式(12)表示每輛車不能超過其最大負(fù)載;式(13)表示設(shè)施點(diǎn)不能超過其最大容量;式(14)表示車輛負(fù)載的連續(xù)性;式(15)為消除子回路約束;式(16)表示車輛離開設(shè)施點(diǎn)時(shí)負(fù)載為0;式(17)～(19)為決策變量取值范圍。

混凝土澆筑主要質(zhì)量控制點(diǎn)之一“澆筑前”，嚴(yán)格控制“人、材、機(jī)”。人員方面，配置有經(jīng)驗(yàn)的工人進(jìn)行混凝土振搗；材料方面，在現(xiàn)場(chǎng)進(jìn)行混凝土坍落度試驗(yàn)，保證混凝土和易性適合肋拱施工；在機(jī)械方面，對(duì)臨時(shí)發(fā)電機(jī)、拌合站、罐車、吊車進(jìn)行檢查保養(yǎng)；倉面使用2個(gè)50 mm軟軸振搗棒，備用2個(gè)50 mm軟軸振搗棒。

2 分布估計(jì)灰狼算法

原始GWO 通過模仿自然界中灰狼的分級(jí)和狩獵制度對(duì)解空間進(jìn)行搜索[20],種群中適應(yīng)度最優(yōu)的個(gè)體被稱作α狼,排名第2及第3的個(gè)體為β狼和δ狼,其余狼均為ω狼。在狩獵(搜索)階段,ω狼同時(shí)利用3個(gè)領(lǐng)導(dǎo)狼的信息來更新自身位置,逐步逼近獵物(最優(yōu)解),進(jìn)行捕獵[20]。

GWOEDA 在GWO 的基礎(chǔ)上加入EDA,以提高算法對(duì)于路徑問題的求解性能,在整個(gè)灰狼種群的進(jìn)化過程中,每個(gè)個(gè)體中的路徑編碼依據(jù)EDA估計(jì)得到的概率模型來生成。GWOEDA 的基本思路為:首先,初始化概率模型生成初始種群;然后,利用灰狼算法的進(jìn)化規(guī)則更新ω狼,并對(duì)領(lǐng)導(dǎo)狼執(zhí)行鄰域搜索進(jìn)一步改進(jìn)優(yōu)良解的質(zhì)量;最后,根據(jù)種群中的精英個(gè)體更新概率模型。上述過程循環(huán)迭代直至滿足終止條件。

2.1 編碼方式

使用整數(shù)編碼表示個(gè)體(可行解),每個(gè)個(gè)體包含一條或多條配送路徑;每條路徑由一組需求點(diǎn)和一個(gè)設(shè)施點(diǎn)組成,表示某一配送車輛從該設(shè)施點(diǎn)出發(fā),按照需求點(diǎn)的排列順序完成收集任務(wù)。例如:給定具有3個(gè)候選設(shè)施點(diǎn)和10個(gè)需求點(diǎn)的LCLRP,圖1給出了該問題的一個(gè)可行解編碼。該可行解由3條收集路徑組成,路徑1車輛由11號(hào)設(shè)施出發(fā),完成對(duì)需求點(diǎn)6和5的收集任務(wù)后返回11號(hào)設(shè)施;路徑2和3均從12號(hào)設(shè)施出發(fā),路徑2完成需求點(diǎn)7,1,3的收集任務(wù),路徑3完成需求點(diǎn)2,10,8,4,9的收集任務(wù)。候選設(shè)施點(diǎn)13未出現(xiàn)在該可行解的編碼中,表示其不開放。

圖1 解的表示

2.2 概率模型

利用EDA 算法中的概率模型表示LCLRP 中節(jié)點(diǎn)間的相關(guān)性,本文中的概率模型設(shè)置為二維矩陣,矩陣元素的值越大,表示連接對(duì)應(yīng)兩個(gè)節(jié)點(diǎn)的路徑在精英灰狼中出現(xiàn)的次數(shù)越多。

為了保留迭代中的歷史信息,對(duì)每一代種群采樣得到的矩陣At進(jìn)行累加得到相關(guān)性矩陣,即

式中,ρ為相關(guān)性矩陣的衰減比例,目的是防止算法過早地陷入局部最優(yōu)。

利用距離矩陣D和相關(guān)性矩陣Ct構(gòu)建概率矩陣,即

Ct中的最大值在αt對(duì)應(yīng)位置轉(zhuǎn)化為最小值,以保證其符合距離矩陣D的價(jià)值排序。矩陣DTB中的數(shù)據(jù)元素值越小,表示對(duì)應(yīng)路徑在最優(yōu)解中出現(xiàn)的概率越大。生成初始種群時(shí)僅利用距離矩陣(即h=1)生成初始解,并且由于迭代初期種群的適應(yīng)度值較差,在算法迭代初期,令相關(guān)性矩陣Ct所占比重較小。而隨著迭代次數(shù)增加,相關(guān)性矩陣所占比重逐漸增加且距離矩陣所占比重由h逐漸衰減至0。

2.3 解的生成

利用概率矩陣構(gòu)建可行解,具體流程如下:

步驟1隨機(jī)選擇一個(gè)未被訪問的需求點(diǎn)作為初始點(diǎn);

步驟2構(gòu)造當(dāng)前需求點(diǎn)的允許訪問節(jié)點(diǎn)集合allowk;

步驟3根據(jù)式(24)計(jì)算allowk中節(jié)點(diǎn)的選擇概率,并依概率選擇下個(gè)需求點(diǎn)

步驟4若滿足2.4節(jié)中的路徑終止條件,則依據(jù)2.5節(jié)中的設(shè)施分配原則,為該配送路徑分配設(shè)施,并返回步驟1;否則,轉(zhuǎn)步驟2。

2.4 路徑終止條件

LCLRP中一條路徑的運(yùn)輸成本與車輛負(fù)載呈正相關(guān),因此,在車輛接近滿載時(shí)繼續(xù)訪問下一節(jié)點(diǎn)所帶來的成本可能會(huì)超過啟用新車輛。為此,需要控制車輛的平均載重率,平衡每輛車的行駛距離和負(fù)載。

通過路徑終止條件控制車輛的平均載重率。路徑終止概率p由下式計(jì)算得到,

若訪問下一需求點(diǎn)將違反車輛容量約束,則路徑終止概率p=1,終止路徑;若尚未違反車輛容量約束,則生成[0,1]的隨機(jī)數(shù)rand,當(dāng)rand

式(25)中,dtemp為下一個(gè)將被加入路徑的節(jié)點(diǎn)的需求量,通過調(diào)整式(25)中的參數(shù)a,可以控制車輛的平均載重率,且不會(huì)固定所有解中路徑的數(shù)量。

2.5 倉庫分配原則

對(duì)所有可用設(shè)施點(diǎn),一方面,按其已被使用容量進(jìn)行升序排列,得到每個(gè)可用設(shè)施點(diǎn)的排序值rankCT;另一方面,按其與當(dāng)前配送路徑中所有需求點(diǎn)的距離之和降序排列,得到每個(gè)可用設(shè)施點(diǎn)的排序值rankDT。綜合考慮每個(gè)可用設(shè)施點(diǎn)的排序值rankCT和rankDT,由下式可計(jì)算得到所有可用設(shè)施點(diǎn)對(duì)當(dāng)前路徑的得分GradeT,選擇GradeT值最高的候選設(shè)施為當(dāng)前路徑提供服務(wù),

式中,v為調(diào)整容量與距離在評(píng)分時(shí)所占的權(quán)重。

2.6 多父代交叉

采用多父代交叉[24](Multi-parent Sequential Constructive Crossover,MPSCX)的方法模擬原始GWO 算法中ω狼利用α、β和δ狼的位置信息對(duì)其更新。具體方法為:首先,剔除ω狼與領(lǐng)導(dǎo)狼中的設(shè)施節(jié)點(diǎn),保留需求點(diǎn)的排序不變;其次,隨機(jī)選擇某一待訪問需求點(diǎn)作為一條配送路徑的第1 個(gè)節(jié)點(diǎn),并將ω狼與領(lǐng)導(dǎo)狼中需求點(diǎn)排列順序中位于當(dāng)前需求點(diǎn)之后的待訪問需求點(diǎn)組成允許訪問集合allowk;再次,利用概率矩陣判斷allowk中各需求點(diǎn)的重要程度,依概率選擇所訪問的需求點(diǎn);最后,通過路徑終止條件決定是否終止路徑并分配設(shè)施點(diǎn)。

給定如圖2所示的ω、α、β和δ狼中需求點(diǎn)的排列順序,假設(shè)隨機(jī)選擇8號(hào)需求點(diǎn)為某一配送路徑的初始訪問節(jié)點(diǎn),則8 號(hào)需求點(diǎn)的允許訪問集合allowk為{4,2,6,5}。根據(jù)式(24)得到訪問集合中待訪問需求點(diǎn)4,2,6,5 的選擇概率分別為0.8、0.1、0.04和0.06,若依概率選擇4號(hào)需求點(diǎn)作為下個(gè)訪問的需求點(diǎn),則4號(hào)需求點(diǎn)的allowk為{9,5,2,10},再依據(jù)概率從該允許訪問集合中選擇所訪問的需求點(diǎn)。依此類推,每次進(jìn)行節(jié)點(diǎn)選擇時(shí)判斷是否需要終止路徑,插入候選設(shè)施點(diǎn)。

圖2 多父代交叉

2.7 兩點(diǎn)交換

兩點(diǎn)交換算子首先采用輪盤賭方法選擇一條路徑,并從中選擇一個(gè)需求點(diǎn),將其與另一條路徑中的節(jié)點(diǎn)交換,步驟如下:

步驟1計(jì)算出個(gè)體中所有車輛在各自路徑上的花費(fèi),并依據(jù)輪盤賭法分配選擇概率。

步驟2依概率選擇其中一條路徑,選擇該路徑中距離設(shè)施點(diǎn)最遠(yuǎn)的節(jié)點(diǎn)記作p1,該路徑所使用的設(shè)施點(diǎn)記作F1。找到距離p1最近的設(shè)施點(diǎn)記為F2。

步驟3計(jì)算由F2服務(wù)的所有需求點(diǎn)對(duì)應(yīng)的scores,其值由與F2及F1距離的差值決定,依概率選擇scores高的節(jié)點(diǎn)作為p2,即

步驟4交換p1、p2的位置,若滿足約束條件,則輸出交換后的解;若不滿足約束條件,則輸出初始解。

圖3是一個(gè)兩點(diǎn)交換的例子。假設(shè)對(duì)于3條路徑的運(yùn)輸成本(路徑1 >路徑2 >路徑3),用輪盤賭法從路徑1中選擇距離其設(shè)施點(diǎn)11較遠(yuǎn)的節(jié)點(diǎn)7。接著利用式(27)計(jì)算分?jǐn)?shù)得到節(jié)點(diǎn)8作為被交換節(jié)點(diǎn),交換7與8并判斷編碼的可行性。

圖3 兩點(diǎn)交換算子

2.8 單點(diǎn)插入

單點(diǎn)插入會(huì)依概率選擇成本高的路徑,從中選擇距離其設(shè)施點(diǎn)最遠(yuǎn)的需求點(diǎn)插入其他路徑。

步驟1計(jì)算出個(gè)體中所有車輛在各自路徑上的成本,并依據(jù)輪盤賭法選擇花費(fèi)較大的路徑。

步驟2依概率選擇其中一條路徑,選擇該路徑中距離設(shè)施點(diǎn)最遠(yuǎn)的節(jié)點(diǎn)記作p1,該路徑所使用的設(shè)施點(diǎn)記作F1。找到距離p1最近的設(shè)施點(diǎn)記為F2。

步驟3隨機(jī)選擇F2中的一條路徑插入p1。

圖4中存在3條路徑,依據(jù)各自的路徑成本用輪盤賭法選出路徑1,找到其中距離設(shè)施點(diǎn)最遠(yuǎn)的節(jié)點(diǎn)7,距離節(jié)點(diǎn)7最近的設(shè)施點(diǎn)為12,隨機(jī)將節(jié)點(diǎn)7插入12所在的路徑中。

圖4 單點(diǎn)插入算子

2.9 算法流程

GWOEDA 的具體流程如下:

步驟1根據(jù)式(22)初始化概率模型,此時(shí)h=1,即等價(jià)于距離矩陣。

步驟2利用貪心算法,根據(jù)2.3節(jié)中的個(gè)體生成規(guī)則生成初始種群。

步驟3對(duì)種群中的個(gè)體進(jìn)行排序,確定領(lǐng)導(dǎo)狼與被領(lǐng)導(dǎo)狼,并采樣種群,更新概率圖和系數(shù)h。

步驟4更新ω狼的位置,生成隨機(jī)數(shù)rand(0,1),若小于系數(shù)p,則使用多父代交叉算子進(jìn)行引導(dǎo),否則使用生成初始種群的方法替換該狼,以增加種群多樣性。

步驟5使用兩種鄰域搜索算子對(duì)領(lǐng)導(dǎo)狼改進(jìn)次。改進(jìn)后的結(jié)果若比之前更優(yōu),則改進(jìn)被接受;若劣于原個(gè)體,則以概率b接受。

步驟6滿足終止條件,則輸出結(jié)果,否則返回步驟3。

2.10 時(shí)間復(fù)雜度分析

從種群規(guī)模popsize,迭代次數(shù)iteration,父代數(shù)量parent,路徑截?cái)嘞禂?shù)k,以及鄰域搜索系數(shù)c和算法流程評(píng)估最壞時(shí)間復(fù)雜度。

在算法主程序調(diào)用的所有函數(shù)中:初始化距離矩陣的時(shí)間復(fù)雜度為O(n2);初始化個(gè)體編碼的時(shí)間復(fù)雜度為O(n+f(k)),其中f(k)為不同的系數(shù)k所決定的路徑數(shù)量,由此可得生成初始種群的時(shí)間復(fù)雜度為O(popsize(n+f(k)));對(duì)每個(gè)個(gè)體進(jìn)行解碼的時(shí)間復(fù)雜度為O(popsize(n+f(k)));利用快速排序?qū)ΨN群個(gè)體進(jìn)行排序的時(shí)間復(fù)雜度為O(popsiez log popsize);對(duì)非父代個(gè)體一次MPSCX 算子的時(shí)間復(fù)雜度為O((n+f(k))(n+f(k)+1)/2);兩點(diǎn)交換與單點(diǎn)插入算子的時(shí)間復(fù)雜度均為O(2n)。記popsize=S,iteration=I,parent=P,計(jì)算得到GWO 的時(shí)間復(fù)雜度為

式中,g(c)為算法中由系數(shù)c所決定的鄰域搜索次數(shù)。

GWOEDA 在GWO 的基礎(chǔ)上每次迭代更新一次分布矩陣,其時(shí)間復(fù)雜度為O(n2),則GWOEDA的時(shí)間復(fù)雜度為

3 算例分析

為了驗(yàn)證分布估計(jì)策略對(duì)算法性能的提升效果,使用結(jié)構(gòu)相同的兩種灰狼算法進(jìn)行對(duì)比,一種僅利用距離矩陣作為決策信息(固定式(22)中的h=1),另一種利用分布估計(jì)迭代更新的概率矩陣作為決策信息。所有程序均由Matlab R2019b編寫,并運(yùn)行在配置為16.0 GB RAM,AMD Ryzen 3 700X 3.6 GHz,操作系統(tǒng)為64位的Windows 10主機(jī)。

3.1 測(cè)試算例

為了驗(yàn)證基于分布估計(jì)策略的混合灰狼算法在求解LCLRP 模型上的有效性,實(shí)驗(yàn)部分求解了CLRP 標(biāo)準(zhǔn)算例 (http://prodhonc.free.fr/Instances/instances_us.htm)中的Prins數(shù)據(jù),并在其基礎(chǔ)上進(jìn)行相應(yīng)改造來進(jìn)行實(shí)驗(yàn)。

對(duì)于Prins標(biāo)準(zhǔn)算例,其運(yùn)輸成本計(jì)算公式為

將其改造后的運(yùn)輸成本計(jì)算公式為

根據(jù)對(duì)現(xiàn)階段國(guó)內(nèi)市場(chǎng)情況的調(diào)研,式(31)中:單位燃油價(jià)格rk=250;燃油轉(zhuǎn)換碳排放系數(shù)ηk=2.63;單位碳排放定價(jià)φ=20。式(4)中:空載單位距離油耗ρ0=0.77;滿載單位距離油耗ρ*=1.54。

表1 給出了標(biāo)準(zhǔn)算例庫中測(cè)試所用的8 個(gè)算例,為了區(qū)別原始CLRP 算例與經(jīng)改造所得的LCLRP算例,將LCLRP算例編號(hào)用*加以標(biāo)記。

表1 算例數(shù)據(jù)

3.2 參數(shù)設(shè)置

為了得到最優(yōu)的參數(shù)組合,對(duì)GWO 與GWOEDA分別采用L16(43)和L25(45)正交表進(jìn)行參數(shù)選取,選用算例3作為測(cè)試數(shù)據(jù)。由于MPSCX算子對(duì)于父代個(gè)數(shù)有相當(dāng)程度的敏感性[24],故算法中對(duì)原始灰狼算法中的領(lǐng)導(dǎo)狼個(gè)數(shù)進(jìn)行調(diào)整,使其作為參數(shù)進(jìn)行選擇。表2 中,popsize 為種群規(guī)模,parent為領(lǐng)導(dǎo)狼個(gè)數(shù),h為分布初始化時(shí)距離矩陣所占比重,p為選擇全局搜索算子的系數(shù),c為影響鄰域搜索次數(shù)的系數(shù)(系數(shù)越小,每次迭代中鄰域搜索的比重越大)。參數(shù)水平如表2所示。

表2 參數(shù)水平

本文主要列出GWO 與GWOEDA 在不同parent水平上的參數(shù)分析。圖5所示為兩種算法在不同parent水平下的正交實(shí)驗(yàn)結(jié)果。圖中曲線顯示了算法求解性能隨著parent數(shù)量增加而提高的過程,同時(shí),隨著parent數(shù)量的增加,算法運(yùn)行時(shí)間也有一定增長(zhǎng)。

圖5 GWO 與GWOEDA 的parent參數(shù)水平分析

為了更進(jìn)一步分析parent對(duì)算法性能與求解時(shí)間的影響,將在其余參數(shù)確定的情況下擴(kuò)展parent的取值范圍進(jìn)行算法效率分析。實(shí)驗(yàn)方法如下:固定一組相同的初始種群,設(shè)置迭代次數(shù)為1 000次。記錄不同parent水平下每次迭代的運(yùn)行時(shí)間tn和當(dāng)前迭代搜索到的最優(yōu)解bestn。記錄每次迭代中最優(yōu)解的變化impn=bestn-bestn-1,并得到最優(yōu)解在單位時(shí)間內(nèi)的變化率efficiencyn=impn/tn。

圖6所示為GWO 與GWOEDA 在parent水平6、9、12、15、18、21下算法運(yùn)行的累計(jì)時(shí)間曲線。不同的parent取值在很大程度上影響算法的求解時(shí)間,parent=21的計(jì)算時(shí)間相比parent=12增加了50%～90%,相比parent=6 增加了150%～180%。

圖6 GWO 和GWOEDA 不同parent下的求解時(shí)間

圖7、8分別刻畫了GWO 和GWOEDA 在不同parent水平下搜索到的最優(yōu)解bestn以及對(duì)應(yīng)的最優(yōu)解變化率。不同parent參數(shù)在兩種算法中的表現(xiàn)較為一致:較小的取值(parent=6)在時(shí)間成本上表現(xiàn)良好,但在最優(yōu)解bestn和最優(yōu)解變化率efficientn上表現(xiàn)不佳;而較大的取值(parent=21)通?？梢允顾惴ㄔ诘跗诳焖偈諗?但相對(duì)的是高額的時(shí)間成本,并且隨著算法迭代次數(shù)的上升,在最優(yōu)解的表現(xiàn)上與其余取值相差不大。當(dāng)parent=12時(shí),算法在求解能力bestn上與更高數(shù)量的父代表現(xiàn)相近,并且由于其較小的時(shí)間成本,在迭代次數(shù)接近1 000時(shí)具有最佳的求解效率。

圖8 GWO 和GWOEDA 不同parent下的求解效率

綜上可知,取parent=12 作為最終取值,完整的參數(shù)設(shè)置如表3所示。

表3 參數(shù)設(shè)置

3.3 分布估計(jì)策略性能分析

分別使用GWO 與GWOEDA 對(duì)8個(gè)CLRP算例和8個(gè)LCLRP算例進(jìn)行求解。將每個(gè)算例使用兩種算法獨(dú)立運(yùn)行20次,最大迭代次數(shù)為800,求解結(jié)果如表4、5所示。其中,Avg為20次的平均結(jié)果,Max和Min表示最好與最差結(jié)果。

表5 GWO與GWOEDA在LCLRP數(shù)據(jù)集上的結(jié)果

表4、5 的結(jié)果表明,在標(biāo)準(zhǔn)CLRP 與LCLRP中GWOEDA 都表現(xiàn)出了更好的效果,對(duì)于8個(gè)標(biāo)準(zhǔn)CLRP算例以及與其相應(yīng)的8個(gè)LCLRP 算例,GWOEDA 對(duì)絕大多數(shù)算例的求解結(jié)果均優(yōu)于GWO 算法的求解結(jié)果。這說明,使用概率模型的GWOEDA 算法不僅在尋優(yōu)能力上更有優(yōu)勢(shì),而且算法的自適應(yīng)性更強(qiáng)。由于加入了分布估計(jì)的步驟,GWOEDA 相比GWO 在求解時(shí)花費(fèi)了更多時(shí)間,但仍在可接受范圍內(nèi)。

3.4 低碳模型的效用分析

為了更好地說明本文所構(gòu)建的LCLRP模型的可行性及其對(duì)物流配送減少碳排放量的實(shí)際應(yīng)用性,將利用GWOEDA 算法求解標(biāo)準(zhǔn)CLRP算例求得的最優(yōu)配送路徑代入所構(gòu)建的LCLRP 模型中,計(jì)算得到的目標(biāo)函數(shù)值如表6所示。表6給出了標(biāo)準(zhǔn)CLRP 和LCLRP 的最優(yōu)目標(biāo)函數(shù)值,并將CLRP所得到的最優(yōu)路徑代入LCLRP 的目標(biāo)函數(shù)中得到結(jié)果RS。比較表6給出的3組目標(biāo)函數(shù)值,不難發(fā)現(xiàn),LCLRP 算例的最優(yōu)目標(biāo)函數(shù)值遠(yuǎn)小于CLRP 算例得到的路徑。這表明,本文所構(gòu)建的LCLRP模型可以有效降低物流配送中的碳排放成本,從而降低物流配送總成本。

表6 CLRP的最優(yōu)解在LCLRP上的成本

3.5 GWOEDA的有效性評(píng)估

為了說明GWOEDA 在求解CLRP上的性能,本節(jié)用GWOEDA 求解了10組Barreto算例,并將其求解結(jié)果與GRASP[25]、LRGTS[26]、2-phaseHGTS[27]、HGA[28]、HPSO[29]和GWO 的求解結(jié)果對(duì)比分析,如表7 所示。其中,BKOV 是算例已知的最優(yōu)解,BOV 是算法求解到的最好結(jié)果,T為運(yùn)算時(shí)間。GRASP、LRGTS和2-phaseHGTS由C++編寫,運(yùn)行環(huán)境參照文獻(xiàn)[25-27],HGA、HPSO、GWO 和GWOEDA 均 由Matlab 編寫,并運(yùn)行在配置為16.0 GB RAM,AMD Ryzen 3700X 3.6 GHz的主機(jī)上。由于后4種算法為進(jìn)化算法,為了保證公平的比較,本節(jié)中將終止條件統(tǒng)一設(shè)置為:最好結(jié)果經(jīng)過n次迭代未有改進(jìn),其中n為相應(yīng)算例規(guī)模,此時(shí)表明算法已找到全局或局部最優(yōu)解。

表7 Barreto算例的求解結(jié)果

由表7 可知,GWOEDA 與HPSO 均能得到9個(gè)最好結(jié)果,在性能上領(lǐng)先其余算法。在時(shí)間效率層面,GWOEDA 雖然得到了規(guī)模在50以下問題的所有最優(yōu)解,但在小規(guī)模問題的求解時(shí)間上表現(xiàn)并不突出。與GWO 對(duì)比時(shí),GWOEDA 對(duì)于更大規(guī)模的問題在性能和效率上會(huì)表現(xiàn)更好?？赡艿脑蚴?由于分布估計(jì)策略需要對(duì)每一代的種群進(jìn)行漸進(jìn)學(xué)習(xí),同時(shí)按照概率選擇生成個(gè)體,故搜索能力雖然更強(qiáng),但收斂效率較慢。

4 結(jié)語

針對(duì)碳定價(jià)背景下的低碳選址-路徑問題進(jìn)行建模,并設(shè)計(jì)了一種基于分布估計(jì)策略的混合灰狼算法對(duì)其進(jìn)行求解。算法通過構(gòu)建概率分布模型引導(dǎo)狼群搜索,有效地提高了算法的搜索性能。算例分析表明,本文所構(gòu)建的LCLRP 模型能夠有效降低物流配送中的碳排放量,所設(shè)計(jì)的GWOEDA 算法具有較好的優(yōu)化性能。

低碳物流配送優(yōu)化問題本身屬于多目標(biāo)優(yōu)化問題,并在現(xiàn)實(shí)配送中涉及很多不確定性因素。因此,在后續(xù)研究中將進(jìn)一步構(gòu)建不確定性多目標(biāo)LCLRP 模型,并設(shè)計(jì)其相應(yīng)的GWOEDA 求解算法。