考慮碳排放的兩階段選址-路徑問題及其算法

湯希峰 ，何 杰 ，張 浩

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

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

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

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

1 模型構建

1.1 問題描述

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

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

1.2 數學模型

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

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

2 算法設計

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

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

2.1 第1 階段算法

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

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

2.2 第2 階段算法

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

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

與經典的蟻群算法相同[13]，改進的蟻群算法也包含2 個關鍵步驟：1） 主要負責路線的構建，即基于當前節點逐步選擇下一個將要訪問的節點；2） 主要負責信息素濃度的更新，具體包括局部信息素更新和全局信息素更新.

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

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

當隨機產生的q≤q0時，選擇使取值最大的j作為下一個將要訪問的節點；當q>q0時，則根據的值所占的比例選擇下一個將要訪問的節點j，在TSHA 算法中采用輪盤賭的方式進行選擇.

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

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

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

式中：EIW為當次最差路線方案產生的碳排放.

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

3 算例測試

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

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

3.1 對TSHA 算法思想的驗證

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

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

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

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

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

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

從表2 可以看出，TSHA 在求解考慮碳排放的2E-LRP 大規模算例時，求解結果的最大值、最小值和平均值以及計算時長的最小值和平均值之間相差很小，表現穩定.可見，TSHA 可以作為求解考慮碳排放的2E-LRP 大規模算例的一種有效算法.

4 結 論

1） 相較于經典的以物流成本最小化為目標的2E-LRP 模型，本文基于一種簡單實用的以排放因子為主要參數的碳排放計算方法，建立了以碳排放最小化為目標的2E-LRP 優化模型.

2） 對Prodhon 標準算例集中全部6 個最大規模算例的測試，表明TSHA-Ⅱ算法比現有算法更具兼顧求解結果和求解效率的優勢，也證明了TSHA算法思想的可行性和有效性.

3） TSHA 算法在求解考慮碳排放的2E-LRP 大規模算例時表現得高效而且穩定，能夠滿足實際應用場景的計算要求.

4） 本文提出的模型和算法有助于物流企業通過優化城市物流運作實現車輛的減排目標，服務國家“雙碳”戰略.