基于拉格朗日松弛的鐵路行包運輸方案編制方法研究

王 澤，潭宇燕，魏玉光

(北京交通大學 交通運輸學院， 北京 100044)

我國鐵路行包運輸少數采用行包專列車，多數仍依托旅客列車編掛的行李車，具有全天候、安全、快速、通達范圍廣、低碳環保的優點。旅客列車開行方案依據客流流量、流向以及變化規律確定，考慮行包流的因素較少，因此存在行包運輸能力與行包流量、流向不完全匹配的問題。解決該問題的關鍵是編制行包運輸方案，為行包指定明確的裝運和接續車次。目前，行包運輸方案的編制依舊采用傳統編制原則與人工經驗相結合的方式，自動化水平較低，且難以提高服務水平。

既有研究主要將行包運輸方案編制問題分解為行包運輸路徑生成和行包流分配兩個子問題來考慮。

對于行包運輸路徑生成問題，文獻[1-2]綜合考慮運輸成本、時間、能力和現場作業等因素，給出鐵路行包運輸路徑的基本形式和選擇策略，設計行包運輸路徑搜索算法，以實現對始發和中轉列車的合理選取。文獻[3]基于先直達、后1次中轉和2次中轉的路徑搜索策略，以準裝區段限制、作業接續時間等作為約束條件，設計基于發到站坐標位置網格圖的行包運輸路徑快速算法。文獻[4-5]基于旅客列車運行圖構建動態服務網絡，考慮行包辦理站裝卸作業時間、行包中轉次數等限制因素，以時效性為核心，為每個OD需求生成可行路徑集合。

對于行包流分配問題。文獻[6-7]基于行包合理路徑集，以利潤最大、成本最小或運輸時間最短為目標，考慮行李車載運能力、車站作業能力、裝卸作業時間等約束，建立線性規劃模型，使用求解器求解。文獻[8]針對行包到站分散的特點，以集中到發為目標建立模型，保證行包作業盡可能集中，降低運營管理成本。

構建行包路徑備選集可以簡化行包運輸網絡，但備選集的質量和合理性極大影響求解結果；當部分行李車運能緊張時，固化的路徑集合也無法適應變化條件下的中轉路徑調整。此外，既有研究普遍未能設計有效的求解算法，因此所能求解的模型規模極大依賴于求解器的計算能力，耗時長、精度差。

本文認為該問題實質上是行包流在既有旅客列車服務網絡上的配流問題，其網絡不能僅僅理解為行包運輸的物理網絡，而應當拓展為能體現旅客列車服務網絡特點的時空網絡。

受限于旅客列車運程，長程行包運輸往往需要由相互銜接的旅客列車配合完成，在各旅客列車的銜接點進行中轉作業。少量行包可以在中間站利用途經列車停站時間完成中轉裝卸作業，大批量的行包只能在前一個列車的終到站同時也是接續列車的始發站進行中轉作業。主要采用行包裝卸方式實現中轉，在一定條件下也可以采用行李車整車換掛方式。因此，這一類時空網絡配流問題又具有不同于其他時空網絡配流問題的特點和復雜性。

針對既有研究的不足，本文根據該問題的性質和特點將行包運輸物理網絡拓展為時空網絡，通過引入行包中轉弧，可實現在一定變化條件下行包中轉路徑的調整；為突破大規模網絡問題受限于求解器計算能力的瓶頸，提出基于拉格朗日松弛的求解算法，通過將原始問題分解為相互獨立的子問題，實現模型的高效精確求解；針對拉格朗日下界解不可行的缺點，設計相應的上界啟發式算法，進一步提高求解效率。

1 行包運輸時空網絡構建

行包運輸作業包括發送作業、在途運輸、到達作業3個主要過程，其中在途運輸包括途中運輸和中轉作業兩個子過程。行包運輸方案編制問題的核心是確定在途運輸過程中行包中轉作業的辦理車站以及行包流在各旅客列車間的合理分配。旅客列車按照列車運行圖開行，而列車運行圖是列車具體時空位置的圖解，因此行包的在途運輸過程可以用時空網絡來描述。

A為旅客列車車次集合，a∈A；r為規劃時段內車次a的開行列數序號，r=1,2,…,ra；tarr(a,r,i)、tdep(a,r,i)分別為車次a的第r列車在行包辦理站i的到達、發出時刻；K為行包集合，k∈K；O(k)、D(k)分別為行包k始發站、終到站；tEDT(k)為最早發出時刻，指行包k發送作業的最早完成時刻；tLAT(k)為最晚到達時刻，指行包k到達作業的最晚開始時刻；vk為該批行包的質量,t；ek為行包的最大在途時間,h，ek=tLAT(k)-tEDT(k)。行包作業過程及各項時間關系見圖1。

G′=(V′,E′)為行包物理網絡；V′為行包辦理站集，i，j∈V′；E′為相鄰站間的運行區間集，(i,j)∈E′。引入時間維度t∈T，將其拓展為時空網絡G=(V,E)，其中：V為時空點集，(i,t)∈V；E為時空弧集，(i,j,t,t′)∈E。行包物理網絡與對應的時空網絡見圖2。

時空網絡的構建步驟如下：

Step1將每個行包辦理站i拓展為列車到達層和列車出發層，按照時間順序分別排列列車到達時空點(i,t)∈Varr和出發時空點(i,t)∈Vdep，其中Varr為達到時空點集，Vdep為出發時空點集，以此來表示旅客列車的進站和出站操作。

Step2為每批行包k∈K添加起始時空點(i=O(k),t=tEDT(k))∈Vsource、終到時空點(j=D(k),t′=tLAT(k))∈Vsink及虛擬弧(i=O(k),j=D(k),t=tEDT(k),t′=tLAT(k))∈Evirtual，Vsource、Vsink、Evirtual分別為起始時空點集、終到時空點集、虛擬弧集。起訖時空點既表示其在途運輸過程的開始和結束，也體現了最大在途時間約束；虛擬弧則保證了起訖時空點間的連通性。

Step3添加區間運行弧和列車停站弧。區間運行弧(i,j,t,t′)∈Etrain表示行包隨列車a于時刻t=tdep(a,r,i)從行包辦理站i發出，然后于時刻t′=tarr(a,r,i)到達相鄰行包辦理站j的途中運輸過程；列車停站弧(i,i,t,t′)∈Etrain表示行包隨列車a于時刻t=tarr(a,r,i)到達行包辦理站i，在站?？?，然后于時刻t′=tdep(a,r,i)從該站發出的過程，Etrain為區間運行弧和列車停站弧并集。

Step5添加行包終到弧(i,i,t,t′)∈Esink。該弧表示行包k于時刻t=tarr(a,r,i)隨列車a抵達終到站i=D(k)，其中t′=tLAT(k)為行包的最晚到達時刻。

(1)

區間運行弧和列車停站弧的費用為相應的持續時間，行包終到弧的費用為0，行包虛擬弧的費用為行包的最大在途時間。為減少行包在站停留時間，引入在站懲罰系數βk=(|T|/ek)·(1+0.1gi)≥1，其中：gi為行包辦理站的車站等級(0、1、2、3、4、5分別對應特等站、一等站、二等站、三等站、四等站、五等站)；|T|為規劃時長。行包越緊急、行包辦理站的車站等級越低，其懲罰性越強。使用在站懲罰系數將行包始發弧和中轉弧的費用修正為βk·(t′-t)。各時空弧的費用為

(2)

2 行包運輸方案編制模型

2.1 問題假設和輸入

模型基于以下假設：

①行李車能力假設。僅以重量來衡量行李車的載運能力，不考慮行包體積對載運能力的影響；不考慮行包的混裝限裝規定。

②行包辦理站能力假設。假設行包辦理站的站存能力富裕，能夠滿足行包大量堆放的要求；不考慮裝卸工人、裝卸機械的數量和效率對行包裝卸作業的影響，假設行包在規定的中轉時間內均能完成中轉作業。

③運輸路徑唯一假設。每批行包不可拆分，僅能選擇唯一的運輸路徑來完成途中運輸。

④模型優化目標假設。僅以時間最短作為優化目標，不考慮行包作業各項收支對運輸方案的影響。

模型輸入為：

①各批行包k∈K的始發站O(k)、終到站D(k)、最早發出時間tEDT(k)、最晚到達時間tLAT(k)及重量vk。

②行包物理網絡G′=(V′,E′)。

③規劃時段范圍T內的列車時刻表，包括列車車次a∈A、開行列數r=1,2,…,ra，以及列車在各途經站的到達時刻tarr(a,r,i)和出發時刻tdep(a,r,i)。

2.2 模型構建

行包運輸方案編制問題本質上是大規模有限時空網絡資源利用問題?；谝褬嫿ǖ臅r空網絡，將行包運輸方案編制問題轉化為多商品流問題。通過對時空網絡中各項弧權的合理設定，如費用、時間或兩者的加權求和，可以根據不同的優化目標對目標方程進行改進。以行包方案編制問題中最具代表性的總時間最短作為優化目標。由于最大在途時間約束、始發時間約束和中轉接續時間約束已嵌入時空網絡中，模型僅需考慮行李車載運能力約束和行包中轉次數約束。

目標方程為

(3)

式中：X={xk(i,j,t,t′)}k∈K,(i,j,t,t′)∈E為0-1變量，若行包k選擇時空弧(i,j,t,t′)，則xk(i,j,t,t′)為1，否則為0。

約束條件：

(1)行包流守恒約束

(4)

式中：(i=O(k),t=tEDT(k))、(i=D(k),t=tLAT(k))分別為行包k的起訖時空點。該約束規定了各時空點的流量平衡。

(2)行李車載運能力約束

?(i,j,t,t)∈Etrain

(5)

該約束規定行李車載運的行包重量不可超過其最大載重量。

(3)行包中轉次數約束

(6)

(4)二元變量約束

(7)

由此得到原問題P為

(8)

s.t.

式(4)～式(7)

3 基于拉格朗日松弛的求解方法

在行包運輸網絡中，每批行包的中轉站點在一定范圍內都是可選擇的，中轉自由度的增加使可行徑路集增大，問題的規模也隨之擴大。為此，本文引入拉格朗日松弛技術。

拉格朗日松弛是選擇原問題P中的困難約束(式(5)、式(6))，添加拉格朗日乘子，將其乘積作為懲罰項帶入原目標方程中，從而將原問題分解為多個易于求解的子問題。通過求解拉格朗日松弛問題，可以獲得原問題P的最優邊界。經過拉格朗日乘子的不斷迭代更新，松弛解逐步逼近原問題的最優解[9]。

3.1 原問題的拉格朗日松弛

引入車載能力乘子λ={λ(i,j,t,t′)≥0}(i,j,t,t′)∈Etrain和行包中轉次數乘子μ={μk≥0}k∈K，將行李車載運能力約束和行包中轉次數約束松弛至原目標方程中，得到新的目標方程FLR(X,λ,μ)為

(9)

(10)

由此得到拉格朗日松弛問題LR為

(11)

s.t.

式(4)、式(7)、式(9)

給定拉格朗日乘子，松弛問題LR為|K|個相互獨立的最小費用路徑子問題，可通過標號設定算法求解[10]。求解該問題得到是原問題P的松弛域下界，為獲得最逼近上界可行解的下界，需要構造拉格朗日對偶問題LD為

(12)

s.t.

式(4)、式(7)、式(9)、式(11)

一般采用次梯度方法來求解該問題[11]，通過迭代更新拉格朗日乘子來逐步逼近原問題P的最優解。詳細求解步驟見3.3節。

3.2 上界啟發式算法

由于原問題P的可行域被擴大，拉格朗日對偶問題LD的下界解可能違背部分松弛約束[11]。因此，本節給出啟發式算法來獲得上界可行解。該算法結合下界解中行包流在時空網上的路徑信息，依照最晚到達時間的先后順序，檢索出不滿足中轉次數約束和行李車載運能力約束的行包，然后將其分配至運能充足且時間最短的可行路徑中，從而將不可行解調整為可行解。本算法先對每條時空弧進行行包的試分配，通過將試分配后流量大于能力的飽和時空弧排除，保證了新時空路徑的可行性。若新解優于已知最優可行解，則將其保留。

上界啟發代算法步驟為

Step0初始化

以“tLAT(k)降序”為主序、“vk升序”為輔序，對行包排序得到K′；初始化上界時空弧累計流量為

fUB(i,j,t,t′)←fLB(i,j,t,t′) ?(i,j,t,t′)∈Etrain

初始化上界可行解

Step1檢索超過最大中轉次數限制的行包：

O(|K|·|E|)

Step1.1更新已分配和待分配行包集合為

Step1.2更新時空弧累計流量為

Step1.3中轉弧費用為

Step1.4重置變量為

Step2檢索使行李車能力過載的行包：

O(|Ktrain|2·|E|)

Step2.2時空弧累計流量為

Step2.3重置變量為

Step2.4若fUB(i,j,t,t′)≤Ccap(i,j,t,t′) 則跳出Step2.1，繼續檢驗下一條列車時空弧。

否則返回Step2.1，檢索下一批行包。

Step3行包再分配：

O(|K|·(log2|V|+|Etrain|))

Step3.1識別能力飽和的時空弧：

對每條列車時空弧(i,j,t,t′)∈Etrain，若fUB(i,j,t,t′)+vk>Ccap(i,j,t,t′)，則ck(i,j,t,t′)←+∞。

Step3.4更新時空弧累計流量為

Step4計算上界并更新上界可行解：

O(|K|·|E|)

算法結束。

3.3 拉格朗日求解算法

在每步迭代中，拉格朗日求解算法基于當前各時空弧的懲罰費用(車載能力乘子λ和中轉次數乘子μ)更新弧權、分配各批行包至最小費用路徑中，并計算中轉次數、累計流量以及下界值。通過調用3.2節中的上界啟發式算法，下界解被調整為新的上界可行解。若上、下界值收斂至容許誤差范圍內，則算法結束。否則，基于當前行包流對行李車載運能力約束和行包中轉次數約束的違反程度，次梯度和拉格朗日乘子將得到更新，以作為各時空弧新的懲罰費用。若迭代次數達到設定的最大值，算法結束，否則將進入下一輪迭代。

拉格朗日求解算法步驟為

初始化：初始化迭代步數和最優下界值：n←0；zLB←-∞；初始化車載能力乘子和中轉次數乘子為

初始化上界可行解和上界值為

zUB←+∞

Step1對每批行包執行Step1.1 ～ Step1.3：O(|K|·(log2|V|+|E|))

Step1.1更新時空弧費用

區間運行弧和列車停站弧為

行包中轉弧為

?(i,i,t,t′)∈Etr

Step1.2搜尋最小費用路徑得到

Step1.3更新中轉次數為

Step2更新時空弧累計流量為

Step3更新下界值

Step4更新上界可行解為

O[|K|·|E|+|Ktrain|2·|E|+|K|·(log2|V|+|Etrain|)]

Step5計算誤差率：若(zUB-zLB)/zUB≤εgap，算法結束，否則執行Step6。

Step6更新次梯度及迭代步長：

O(|K|+|Etrain|)

車載能力次梯度為

中轉次數次梯度為

Step7更新乘子：O(|K|+|Etrain|)

車載能力乘子為

中轉次數乘子為

Step8更新迭代步數：n←n+1；若n>N，算法結束，否則返回Step1。

本文模型為多商品流問題，屬二元整數組合規劃，因此必存在有限最優解[9]。拉格朗日算法是針對該類問題的一種有效算法，已被廣泛應用。但在求解過程該算法松弛了模型的部分困難約束，擴大了相應可行域，導致拉格朗日求解算法收斂至平衡態時未必能夠得到原問題的最優可行解，僅為其下界[11]。盡管本文提出的上界啟發式算法能彌補解不可行的不足，但也依賴于下界解的質量，難以保證解的最優性。

4 案例驗證

相關算法基于Python編程語言實現，所有實驗均在一臺Intel Core i7-9750 H CPU @2.60 GHz, 16 GB RAM的個人計算機上進行。

4.1 小規模案例

以包含8個行包辦理站、單日10對列車的小規模網絡驗證算法的計算效率。線站示意圖見圖3。

列車時刻表見表1。每批行包的重量和初始行李車載運能力分別從均值為1.5，標準差為0.3以及均值為5，標準差為1.5的正態分布中隨機抽樣產生，單位為噸。運到期限按400 km內為3 d、每增加400 km遞增1 d的方法確定[12]。各辦理站間的運價里程和運到期限見圖4。假定行包的發送和到達作業共耗時1 d，則最大在途時間為其運到期限減去1 d。最早發出時間均為規劃時段首日的18:00。每列車編掛一輛行李車。最大在途時間為2～3、4～5、6 d及以上的行包最大中轉次數分別為1、2、3。最大迭代次數N為200，容許誤差率εgap為5%，中轉弧懲罰系數γ取1.5，系數α初值取2。求解器GUROBI 9.0保持默認設置。

表1 小規模案例列車時刻表(單日)

為詳細對比拉格朗日算法和求解器GUROBI的計算效率，設計4組不同規模的實驗。各組模型規模見表2,實驗信息見表3，時空網絡規模見表4。

表2 小規模案例模型規模

表3 小規模案例實驗信息

表4 小規模案例時空網絡規模

不同規模網絡下算法上界、下界隨迭代次數變化曲線見圖5，其中下界與上界分別可在第75、35步迭代時趨于穩定。誤差率隨迭代次數變化曲線見圖6。由圖6可知，不同規模網絡下的上下界誤差率均可在第75步迭代時達到6%左右，證明算法的收斂性較好。圖7為不同規模下兩者達到相同誤差率的耗時對比。由圖7可見，隨著規模的增大，拉格朗日求解算法較GUROBI的計算時間增長較為緩慢，由此證明了拉格朗日求解算法的高效性。若增大步長系數α的初值，預計算法的收斂效果會更好。

4.2 大規模實例

以哈爾濱鐵路局和沈陽鐵路局集團有限公司管轄范圍內的行包運輸網絡為背景，驗證模型和算法對真實案例的適用性。本實例包含行包辦理站205個，單日列車299對，行包500批。|T|=8 d。列車停站時間小于3 min則視作從該站通過。車站等級與各站單日經停列車數見圖8、圖9。最大在途時間為2～3、4～6 d的行包最大中轉次數分別為2、3、6 d及以上的行包不限制中轉次數。由于求解效率過低且內存占用較大，不使用GUROBI進行求解。為減少求解時間，設置拉格朗日步長系數α初值為5，并采取前100次迭代中每10步、后續每50步調用一次上界啟發式算法的求解策略。其他設定同4.1節。

時空網絡規模與模型規模見表5、表6。

表5 大規模實例時空網絡規模

表6 大規模實例模型規模

經300次迭代，耗時5 h 50 min，求得的上界、下界分別為16 734.375、15 607.310 t·h，誤差率為8.735%，上、下界及誤差率隨迭代次數的變化曲線見圖10。算法的收斂效果良好。

最大在途時間側面反映了行包的運距。運輸時間情況見圖11。由圖11可知，運距每增長400 km，運輸時間約增加4～5 h，在站停留時間約增加2 h。中轉情況見圖12,由圖12可知，隨著運距的增加，直達行包占比下降，中轉行包占比上升。

5 結束語

基于旅客列車行李車的時空特性，本文將行包運輸方案編制問題轉化為基于時空網絡的多商品流問題；并以時間最短為目標，考慮行李車載運能力、中轉次數及各項時間約束，建立了混合整數規劃模型。針對模型規模龐大、求解困難的問題，設計了基于拉格朗日松弛的求解算法，將原問題分解為一系列易于求解的子問題；并給出了上界可行解啟發式算法。

算例結果表明，與商用求解器相比，本文提出的模型與算法具備高效求解大規模行包運輸方案編制問題的能力。若算法采用更為高級的C++語言編程實現，且在性能更強的工作站中引入并行計算技術，計算效率將進一步提升。

在未來研究中，我們會進一步考慮行包辦理站裝卸作業能力和倉儲能力對行包運輸方案編制問題的影響。