湯希峰 ，何 杰 ，張 浩

（1.河海大學(xué)土木與交通學(xué)院，江蘇 南京 210098；2.東南大學(xué)交通學(xué)院，江蘇 南京 211189）

城市是生產(chǎn)、生活的主要集中地，也是物流活動(dòng)最為集中的地方.兩層的物流網(wǎng)絡(luò)作為一種可行的城市物流解決方案正被國內(nèi)外許多大中城市所采用[1-4]，其核心問題之一是兩階段選址-路徑問題 （2ELRP）.經(jīng)典的2E-LRP 可以描述為：重型或中型貨車負(fù)責(zé)將貨物從一個(gè)位于城郊的物流園區(qū)運(yùn)送至若干位于城市中心區(qū)周邊的配送中心，貨物在這些配送中心經(jīng)過重新配裝后再由輕型或微型貨車送達(dá)終端客戶，其目標(biāo)是尋求最優(yōu)的配送中心選址以及第一階段和第二階段貨車的最優(yōu)路徑以實(shí)現(xiàn)總的物流成本最小化[5].

眾所周知，貨車相較于客車容易導(dǎo)致交通事故和交通擁堵，更是污染物和碳排放的主要貢獻(xiàn)者.以我國為例，2019 年全國貨車保有量雖然只占汽車總量的10.71%，貨車CO、HC、NOx、PM 排放量卻占到了汽車排放總量的29.7%、26.3%、83.5%、90.1%[6].隨著生態(tài)環(huán)境保護(hù)問題日漸受到重視，如何優(yōu)化2 層網(wǎng)絡(luò)結(jié)構(gòu)的物流運(yùn)作以減少貨車的碳排放引起了國內(nèi)外學(xué)者越來越多的研究興趣.Govindan 等[7]圍繞生鮮配送，建立了以物流成本最小化和環(huán)境影響最小化為目標(biāo)的2E-LRP 多目標(biāo)優(yōu)化模型，并設(shè)計(jì)了多目標(biāo)粒子群優(yōu)化算法（MOPSO）和改進(jìn)的多目標(biāo)可變鄰域搜索算法（AMOVNS）2 種求解算法；胡大偉等[8]針對燃油車和電動(dòng)車不同的行駛與排放特性，分別建立了燃油車和電動(dòng)車的2 階段開放式選址-路徑問題模型，并提出了一種改進(jìn)的模擬退火算法；劉成清等[9]以經(jīng)典的2 階段開放式選址-路徑問題模型為基礎(chǔ)，通過考慮車輛的燃油消耗和CO2排放以及車輛路徑選擇的靈活性，建立了基于路徑靈活性的兩階段開放式低碳選址-路徑問題模型，并運(yùn)用Cplex進(jìn)行了求解；Pitakaso 等[10]以燃油消耗最小化為目標(biāo)，建立了考慮碳排放的綠色兩階段選址路徑問題（G2ELRP）優(yōu)化模型，并設(shè)計(jì)了相應(yīng)的可變鄰域策略的自適應(yīng)搜索（VANSAS）算法.

一方面，鑒于車輛的碳排放受到車型、車速、載重、行駛距離以及駕駛行為等諸多因素的影響[11-12]，現(xiàn)有文獻(xiàn)大多基于車輛的燃油消耗進(jìn)行碳排放的計(jì)算，所需參數(shù)較多，實(shí)用性不強(qiáng)；另一方面，現(xiàn)有的用于求解各種2E-LRP 的算法多適用于求解小規(guī)模問題 （客戶數(shù)量n≤100, 候選配送中心數(shù)量m<10），求解大規(guī)模問題 （n>100,m≥10） 時(shí)計(jì)算耗時(shí)較長，很難實(shí)際應(yīng)用.為此，本文基于一種簡單實(shí)用的以排放因子為主要參數(shù)的碳排放計(jì)算方法，建立了考慮碳排放的2E-LRP 優(yōu)化模型，并針對所研究問題的計(jì)算復(fù)雜性，設(shè)計(jì)了一種可快速求解大規(guī)模問題的兩階段混合算法.

1 模型構(gòu)建

1.1 問題描述

基于圖論，考慮碳排放的2E-LRP 可以描述為G={V,E}，V=V0∪Vs∪Vc，V0為物流園區(qū)，Vs= （1,2,…,m）為配送中心候選集，Vc= （m+1,m+2,…,m+n）為客戶集；E=E1∪E2，E1= {（i,j）:i,j∈V0∪Vs;i≠j}為物流園區(qū)與配送中心之間的路線集，E2= {（i,j）:i,j∈Vs∪Vc;i,j?Vs×Vs;i≠j} 為配送中心與客戶以及客戶與客戶之間的路線集；H為可供物流園區(qū)使用的重型或中型貨車的集合，每臺(tái)重型或中型貨車具有相同的載重能力，記為QH，滿載和空載時(shí)的碳排放系數(shù)分別為εvfH、εveH（單位： kg/km, 下同）；K為可供配送中心共用的輕型或微型貨車的集合，每臺(tái)輕型或微型貨車也具有相同的載重能力，記為QK（QK

圖1 2E-LRP 示意Fig.1 An illustration of the 2E-LRP

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

若用qijh和qijk分別為貨車h（h∈H）在路線 （i,j）∈E1、貨車k（k∈K）在路線（i,j）∈E2上行駛時(shí)的載重量；決策變量xijh=1 表示貨車h（h∈H） 經(jīng)過路線 （i,j）∈E1，否則xijh=0；yijk=1 表示貨車k（k∈K）經(jīng)過路線 （i,j）∈E2，否則yijk=0；zi=1 表示配送中心i（i∈Vs） 被選用，否則zi=0；uij=1 表示客戶j由配送中心i提供服務(wù)，否則uij=0；則以碳排放最小化為目標(biāo)的2E-LRP 數(shù)學(xué)模型可表示為

式 （1） 為目標(biāo)函數(shù)，表示從物流園區(qū)到配送中心以及配送中心到客戶2 個(gè)階段的貨車碳排放量之和最小，與經(jīng)典2E-LRP 模型的目標(biāo)函數(shù) （關(guān)于車輛通行成本或車輛行駛距離的函數(shù)） 不同，是一個(gè)關(guān)于車輛載重和車輛行駛距離的函數(shù)；式（2）、（3） 表示2 個(gè)階段中某一貨車要么不被使用，要么至多被使用一次；式（4）、（5） 保證2 個(gè)階段中貨車路線的連續(xù)性，同時(shí)保證貨車完成任務(wù)后返回原先出發(fā)的物流園區(qū)或配送中心；式（6） 表示客戶只能由一個(gè)配送中心提供服務(wù)；式（7） 表示客戶若由某一配送中心提供服務(wù)，則必被從該配送中心出發(fā)的貨車訪問一次；式（8） 表示某一配送中心要么不被選用，若被選用，所服務(wù)的客戶總需求量不能超過其作業(yè)能力；式（9）、（10） 表示2 個(gè)階段中貨車服務(wù)的客戶需求均不能超過其載重能力；式（11）、（12） 為2 個(gè)階段中路線的完整性約束；式（13）、（14）表示2 個(gè)階段中貨車載重量在相鄰路段上的數(shù)量關(guān)系；式（15） ～（18）為決策變量的取值約束.

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

由于設(shè)施選址問題和車輛路徑問題都屬于NPhard 問題，所以本文研究的考慮碳排放的2E-LRP是更為復(fù)雜的NP-hard 問題，求解難度很大，尤其是隨著問題規(guī)模的增大，計(jì)算耗時(shí)呈指數(shù)增加.鑒于此，本文設(shè)計(jì)了一種可快速求解大規(guī)模問題的2 階段混合算法（TSHA ），算法流程如圖2 所示.第1 階段，通過將2E-LRP 轉(zhuǎn)化成不考慮車輛路徑的2 階段設(shè)施選址問題 （2E-FLP），直接求解配送中心選址方案和客戶分配方案；第2 階段，基于得到的配送中心選址和客戶分配方案，物流園區(qū)到被選用的配送中心以及每個(gè)配送中心到所分配客戶的車輛路徑問題就變成了獨(dú)立且較易求解的VRP 問題，進(jìn)而加以解決.

圖2 兩階段混合算法的算法流程Fig.2 Flowchart of the TSHA

2.1 第1 階段算法

由于迭代算法在求解過程中非常低效，為提高所建模型的求解效率，TSHA 在第1 階段先將考慮碳排放的2E-LRP 轉(zhuǎn)化成不考慮車輛路徑的2EFLP，模型如下：

式 （19） 為目標(biāo)函數(shù)，表示客戶到配送中心以及配送中心到物流園區(qū)的以需求量為權(quán)重的距離之和最小，約束條件對應(yīng)原問題模型中的式 （6）、（8）、式（15）、（16）.不難看出，上述2E-FLP 模型為線性規(guī)劃模型，可直接運(yùn)用Cplex 求解器進(jìn)行求解，從而得到配送中心選址和客戶分配方案.

2.2 第2 階段算法

基于得到的配送中心選址和客戶分配方案，物流園區(qū)到被選用的配送中心以及每個(gè)被選用配送中心到所分配客戶的車輛路徑問題都變成了獨(dú)立的VRP 問題.為求解以碳排放最小化為目標(biāo)的VRP 問題，TSHA 在第2 階段采用了一種改進(jìn)的蟻群算法，如圖3 所示.

圖3 改進(jìn)蟻群算法的偽代碼Fig.3 Pseudocode of the improved ant colony algorithm

與經(jīng)典的蟻群算法相同[13]，改進(jìn)的蟻群算法也包含2 個(gè)關(guān)鍵步驟：1） 主要負(fù)責(zé)路線的構(gòu)建，即基于當(dāng)前節(jié)點(diǎn)逐步選擇下一個(gè)將要訪問的節(jié)點(diǎn)；2） 主要負(fù)責(zé)信息素濃度的更新，具體包括局部信息素更新和全局信息素更新.

從式 （1） 碳排放的計(jì)算式可以看出，貨車的碳排放與其載重量和行駛距離密切相關(guān).為此，改進(jìn)的蟻群算法采用了一種新的路線構(gòu)建規(guī)則，如式 （20） 、式（21） 所示.

式中：Ω為尚未被訪問的節(jié)點(diǎn)集，Ω∈Vs或Vc；τij為當(dāng)前節(jié)點(diǎn)i與未到訪節(jié)點(diǎn)j之間路線上的信息素濃度；ηij為節(jié)點(diǎn)i與節(jié)點(diǎn)j之間距離的倒數(shù)；α、β為重要度系數(shù)；q0為給定的常數(shù)；q為 [0, 1] 上的隨機(jī)數(shù).

當(dāng)隨機(jī)產(chǎn)生的q≤q0時(shí)，選擇使取值最大的j作為下一個(gè)將要訪問的節(jié)點(diǎn)；當(dāng)q>q0時(shí)，則根據(jù)的值所占的比例選擇下一個(gè)將要訪問的節(jié)點(diǎn)j，在TSHA 算法中采用輪盤賭的方式進(jìn)行選擇.

每只螞蟻都按照上述的偽比例選擇規(guī)則構(gòu)建路線，當(dāng)螞蟻構(gòu)建出一條路線（即達(dá)到貨車的載重限制時(shí)，返回原先出發(fā)的物流園區(qū)或配送中心），然后重新出發(fā)構(gòu)建下一條路線，如此反復(fù)，直至遍訪所有客戶.當(dāng)螞蟻構(gòu)建出一個(gè)完整的路線方案后，對該路線方案中的每條路線進(jìn)行局部信息素更新，第t+ 1 次信息素濃度為

式中：ρ（0 ≤ρ≤1） 為信息素的揮發(fā)系數(shù)；τ0為初始的信息素濃度，取值為1/（mEnns），Enns為基于最近鄰搜索算法得到的貨車碳排放量.

當(dāng)蟻群中每只螞蟻均構(gòu)建出一個(gè)完整的路線方案后，從中選擇一個(gè)碳排放量最小的路線方案作為當(dāng)次最優(yōu)的路線方案，為在當(dāng)次最優(yōu)路線方案和當(dāng)前最優(yōu)路線之間取得平衡，亦即鼓勵(lì)隨后的螞蟻在當(dāng)次最優(yōu)路線和當(dāng)前最優(yōu)路線附近的解空間繼續(xù)搜索，改進(jìn)的蟻群算法采用Ting 等[14]推薦的策略對當(dāng)次最優(yōu)和當(dāng)前最優(yōu)路線方案中的每條路線進(jìn)行全局信息素更新，如式（23）、（24）所示.

式中：EIW為當(dāng)次最差路線方案產(chǎn)生的碳排放.

如此，通過若干次循環(huán)后即可得到所有以碳排放最小化為目標(biāo)的VRP 的最優(yōu)路線方案.考慮到由物流園區(qū)提供服務(wù)的配送中心數(shù)量較少而由每個(gè)配送中心提供服務(wù)的客戶數(shù)量相對較多，對得到的配送中心到客戶的VRP 最優(yōu)路線方案，本文運(yùn)用3 種局部搜索算法對其再進(jìn)行改進(jìn)：cross_relocation 算法將從某個(gè)配送中心出發(fā)路線上的客戶逐一插入從其他配送中心出發(fā)的路線中，若得到的路線可行且碳排放量小于原路線，則用改進(jìn)后的路線替代原路線；cross_swap 算法將從不同配送中心出發(fā)的兩條路線上的客戶兩兩交換并重新生成新的路線，若新路線可行且碳排放量較小，則替代原路線；最后，再對每一條路線運(yùn)用3-opt 算法進(jìn)行改進(jìn).至此，改進(jìn)后的路線方案即可作為最終的最優(yōu)路線方案.

3 算例測試

由于目前還沒有可用于測試2E-LRP 碳排放的標(biāo)準(zhǔn)算例集，所以本文選取了Prodhon 標(biāo)準(zhǔn)算例集[15]中全部6 個(gè)最大規(guī)模的算例稍加修改，用以測試本文所提出的模型和算法.Prodhon 算例集中6 個(gè)最大規(guī)模的算例均包括200 個(gè)客戶和10 個(gè)候選配送中心.算例的修改包括去掉原算例中的配送中心選址費(fèi)用、增加不同類型貨車的碳排放系數(shù)，其中重型或中型貨車滿載和空載時(shí)的碳排放系數(shù)分別為2.392、1.638 kg/km，輕型或微型貨車滿載和空載時(shí)的碳排放系數(shù)分別為1.096、0.772 kg/km，其他數(shù)據(jù)保持不變.

設(shè)計(jì)的TSHA 算法采用MATLAB 2018b 編程并嵌入Cplex V12.7 求解器，改進(jìn)的蟻群算法控制參數(shù)包括：求解物流園區(qū)到被選用的配送中心VRP 問題時(shí)，循環(huán)次數(shù)設(shè)定為選用的配送中心個(gè)數(shù)的1/2，蟻群中螞蟻的數(shù)量與選用配送中心個(gè)數(shù)相同；求解配送中心到被分配客戶的VRP 問題時(shí)循環(huán)次數(shù)均設(shè)定為被分配客戶數(shù)量的1/2，蟻群中螞蟻的數(shù)量均等于被分配的客戶的數(shù)量；重要度系數(shù)α= 2，β= 1；給定常數(shù)q0= 0.5；揮發(fā)系數(shù)ρ= 0.2.針對每個(gè)算例，獨(dú)立運(yùn)行程序20 次，并分別記下計(jì)算結(jié)果和時(shí)間.

3.1 對TSHA 算法思想的驗(yàn)證

為驗(yàn)證TSHA 的算法思想，本文先編寫了基于同樣算法思想的TSHA-Ⅱ算法用以求解Prodhon 算例集中6 個(gè)最大規(guī)模的原算例.TSHA-Ⅱ與TSHA略微不同的是：1） 用物流成本代替碳排放作為目標(biāo)函數(shù)；2） 改進(jìn)的蟻群算法中，路線構(gòu)建不考慮客戶的貨運(yùn)需求量，即式 （20） 、（21） 中去掉啟發(fā)因子dj.TSHA-Ⅱ算法的求解結(jié)果與目前文獻(xiàn)中兩種最新算法的求解結(jié)果對比如表1 所示.表中：B為目前文獻(xiàn)所給出的最優(yōu)解；Cmin為基于某種算法得到的最優(yōu)解；Tmin為某種算法得到最優(yōu)解所需的時(shí)間；G為某個(gè)算法求得的最優(yōu)解與B之間的差異，即G=(Cmin-B)/B.

表1 TSHA-II 算法與LNS-2E、Two-stage 算法求解2E-LRP 原算例的結(jié)果對比Tab.1 Comparison of the results obtained by LNS-2E, two-stage algorithm, and TSHA-II in solving original 2E-LRP instances

通過對比可以看出：TSHA-Ⅱ算法求解大規(guī)模的2E-LRP 算例時(shí)，求解質(zhì)量稍遜于Breunig 等[16]提出的LNS-2E 算法，但明顯優(yōu)于Dai 等[17]提出的Twostage 算法；求解效率雖稍遜于后者，但大大優(yōu)于前者.這表明，TSHA-Ⅱ算法能夠在求解質(zhì)量與求解效率之間取得較好的平衡，可以在較短的時(shí)間內(nèi)得到較高質(zhì)量的解，TSHA 的算法思想可行有效.

3.2 對考慮碳排放的2E-LRP 算例的求解

算法思想通過驗(yàn)證后，TSHA 再被用于求解考慮碳排放的2E-LRP 大規(guī)模算例，計(jì)算結(jié)果如表2所示.表中：Eavg為20 次運(yùn)算得到的平均碳排放量；Emax、Emin分別為20 次運(yùn)算得到的最大和最小碳排放量；Tavg為20 次運(yùn)算的平均耗時(shí).

表2 TSHA 算法求解考慮碳排放的2E-LRP 算例的結(jié)果Tab.2 Results obtained by TSHA in solving 2E-LRP instances considering carbon emissions

從表2 可以看出，TSHA 在求解考慮碳排放的2E-LRP 大規(guī)模算例時(shí)，求解結(jié)果的最大值、最小值和平均值以及計(jì)算時(shí)長的最小值和平均值之間相差很小，表現(xiàn)穩(wěn)定.可見，TSHA 可以作為求解考慮碳排放的2E-LRP 大規(guī)模算例的一種有效算法.

4 結(jié) 論

1） 相較于經(jīng)典的以物流成本最小化為目標(biāo)的2E-LRP 模型，本文基于一種簡單實(shí)用的以排放因子為主要參數(shù)的碳排放計(jì)算方法，建立了以碳排放最小化為目標(biāo)的2E-LRP 優(yōu)化模型.

2） 對Prodhon 標(biāo)準(zhǔn)算例集中全部6 個(gè)最大規(guī)模算例的測試，表明TSHA-Ⅱ算法比現(xiàn)有算法更具兼顧求解結(jié)果和求解效率的優(yōu)勢，也證明了TSHA算法思想的可行性和有效性.

3） TSHA 算法在求解考慮碳排放的2E-LRP 大規(guī)模算例時(shí)表現(xiàn)得高效而且穩(wěn)定，能夠滿足實(shí)際應(yīng)用場景的計(jì)算要求.

4） 本文提出的模型和算法有助于物流企業(yè)通過優(yōu)化城市物流運(yùn)作實(shí)現(xiàn)車輛的減排目標(biāo)，服務(wù)國家“雙碳”戰(zhàn)略.