基于拉格朗日松弛的铁路行包运输方案编制方法研究

王 泽，潭宇燕，魏玉光

(北京交通大学 交通运输学院， 北京 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 模型构建

行包运输方案编制问题本质上是大规模有限时空网络资源利用问题。基于已构建的时空网络，将行包运输方案编制问题转化为多商品流问题。通过对时空网络中各项弧权的合理设定，如费用、时间或两者的加权求和，可以根据不同的优化目标对目标方程进行改进。以行包方案编制问题中最具代表性的总时间最短作为优化目标。由于最大在途时间约束、始发时间约束和中转接续时间约束已嵌入时空网络中，模型仅需考虑行李车载运能力约束和行包中转次数约束。

目标方程为

(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++语言编程实现，且在性能更强的工作站中引入并行计算技术，计算效率将进一步提升。

在未来研究中，我们会进一步考虑行包办理站装卸作业能力和仓储能力对行包运输方案编制问题的影响。