大规模货物列车运行图优化编制的一种拉格朗日松弛启发式算法

2020-04-16倪少权

铁道学报 2020年3期

江峰，倪少权

(1. 西南交通大学交通运输与物流学院，四川成都 610031；2. 西南交通大学全国铁路列车运行图编制研发培训中心，四川成都 610031)

编制列车运行图需考虑各类约束条件限制，其实质是疏解各列车运行线对运输资源的潜在占用冲突，进而确定各列车运行线的时空位置，是NP-hard问题[1-3]。目前货物列车运行图编制依靠人工推线求解，结果依赖始发时刻[1-2]，自动化水平较低。

作为铁路运输组织的核心问题，列车运行图优化引起了国内外学者的广泛关注。文献[4-5]以列车运行图排序模型研究单线货物列车运行图的优化编制方法。文献[6]考虑相邻区间衔接关系及车站间隔时间，优化单线非追踪平行列车运行图。文献[7]根据列车运行轨迹的有限性和相邻列车的强约束性构造时空域窗口，滚动求解列车运行图。文献[3]给定列车优先级，以追踪间隔时间和连发间隔时间为约束，提出多目标综合函数，建立列车运行图优化编制模型并求解。文献[8]采用广度搜索法生成高峰小时周期列车运行图，以深度搜索法添加非周期列车运行线，研究计算机编制周期性列车运行图关键技术。文献[9]以总效益最大为目标建立整数规划模型，对其进行拉格朗日松弛，以动态规划方法求解了包含17个车站，18列旅客列车及8列货物列车的列车运行图。文献[10]考虑路网中总旅行时间最优，基于时空网络构建货物列车径路优化模型，采用基于拉格朗日松弛的启发式算法求解。

大规模条件下列车运行图优化编制问题无多项式时间有效算法，通常需对原问题进行分解以缩小规模、加快求解速度[4-5,7,9-10]。现有研究主要通过列车始发时刻及沿途总停时优化列车旅行速度[2,4-7,9-10]，但部分研究简化了约束条件(如忽略天窗约束)[9-10]，部分研究解决问题规模较小[3,6,8-9]，部分研究缺乏对所得解的优劣评价[4-5,7-8]。

货物列车运行图受旅客列车运行图限制，是典型的非周期列车运行图，现有基于拉格朗日松弛的列车运行图优化编制方法具有优化解质量可评估、模型拓展性好等优点，是非周期列车运行图优化编制的有效方法[9-10]。

本文根据我国货物列车运行图优化编制需要，以时空网络中的路径优化描述列车运行图优化编制问题，通过对时空网络时间轴的离散化处理，将出发、到达间隔时间等线路能力约束及停站需求、天窗限制等运输组织约束表示为时空网络中的节点选择限制，从而疏解不同列车运行线对运输资源的潜在占用冲突。上述问题的实质是一类最短路径求解问题，通过给定各列车运行线初始利润，并对始发时刻调整及总停时延长设置罚数，计算各列车运行线实际利润，以全图所有列车运行线总利润最大为目标构建整数规划模型，根据模型特点对其进行拉格朗日松弛求解，基于拉格朗日松弛问题所得时空网络中各节点罚数，以启发式算法迭代优化货物列车运行图。通过参数设置提出不同优化策略，对京九线北京西至阜阳区段内439列货物列车进行运行图优化编制实验并分析结果，验证了算法有效性。

1 问题描述

货物列车运行图优化编制是在旅客列车运行图基础上，满足各列车始发时间域及总停站时间限制，对各列车始发时刻及全程总停时进行优化的过程。

作为列车具体时刻时空位置的图解，列车运行图可用有向时空网络G=(N,A)描述，其中:点集N中元素代表时空域中某列车运行线所有经过的节点，在每个车站，可分为到达节点U及出发节点W；弧集A中元素代表列车运行线组成部分，由车站到达及出发节点确定。S为全体车站集合;T=(1,…,t)为全体列车集合;Sj=(fj,…,lj)⊆S为列车j∈T所经车站集合(共s站)，列车j的运行线可表示为一定约束下G中自fj起，途经fj+1,…,lj-1至lj止的一条路径，其包含的弧集合为Aj,Aj⊆A、节点集合为Nj,Nj⊆N。

对全图设置一个虚拟出发点σ和虚拟终到点ε标注各列车的始发、终到，则有

N={σ,ε}∪

(W1∪…∪Ws-1)∪(U2∪…∪Us)

列车j的运行线可用与其有关的弧子集Aj表示，其包含有：

(1) 由虚拟出发点σ至列车j始发站出发节点v∈Wfj的始发弧(σ,v)。

(2) 在中间站i,i∈Sj{fj,lj}由该站到达节点u∈Ui至该站出发节点v∈Wi的停留弧(u,v)，其中出发节点v∈Wi由该站到达节点及该列车运行线在站停留时间确定。

(3) 在区间(i,i+1)内由车站i出发节点v∈Wi至i+1站到达节点u∈Ui+1的运行弧(v,u)，在出发节点确定情况下，到达节点u∈Ui+1由列车在区间内的运行时间确定。

(4) 从列车j终点站到达节点u∈Ulj至虚拟终到点ε的终止弧(u,ε)。

其中各弧的时间跨度含义为：运行弧(v,u)时间跨度为列车j在区间(i,i+1)的运行时分，由列车区间运行标尺确定；停留弧(u,v)时间跨度为列车j在车站i的停站时间，当列车在该站通过时取值为0。

一个包含5个车站、5条列车运行线的简单示例，见图1。其中，K1、K2为旅客列车，H1、H2、H3为货物列车，A、B、C为技术站，其余为中间站，设货物列车在技术站作业时间、各站出发、到达间隔时间均为4 min，以时空网络中路径表示各列车运行线，由于旅客列车运行线的限制，引起了H2、H3在b站的技术停站及H3在B站停时的延长。

由于运输资源时空分配的有限性，各列车运行线可能产生占用运输资源的潜在冲突(如图1中H3、K1次在B站，H2、K1次在b站，H2、H3次在b站)，货物列车运行图优化编制的关键是疏解上述潜在冲突，即在各类约束条件及旅客列车运行线造成的时空网络节点选择限制下，受限时空网络中的最短路径求解问题[4,9-10]。

2 货物列车运行图优化编制模型

2.1 问题假设

列车运行图优化编制需考虑多种影响因素，为更好地确定研究范围，针对列车运行线铺画这一列车运行图优化编制本质问题，做出如下假设：

(1) 我国主要干线均已复线化，因此以上下行线路固定使用的复线区间为背景；列车在区间内可追踪运行，追踪间隔时间由相邻车站出发、到达间隔时间最大值决定。

(2) 车站间隔时间与列车运行状态有关，为简化分析，假设各站出发、到达间隔时间固定且已知。

(3) 车站进路安排与车站站型有关，暂不考虑进路安排问题，忽略进路冲突并假设各站能力充足。

(4) 不考虑列车运行图均衡性，假设各列车最优始发时刻给定，并可在一定范围内调整。

(5) 旅客列车运行线不可调整。

上述假设可简化问题建模及数据存储且不影响问题一般性。

2.2 约束条件

为确定时空网络中各节点位置关系，将1 d以分钟为单位离散为1到q，q=1 440，θ(v)为节点v∈N所处时刻，两节点间时间差在θ(u)>θ(v)的情况下，可表示为Δ(v,u)=θ(u)-θ(v)；在θ(u)<θ(v)的情况下，可表示为Δ(v,u)=θ(u)-θ(v)+q，以uv表示Δ(v,u)>Δ(u,v)，即节点v在时空网络中对应时间早于节点u。

且Nj与Aj间一一对应。

对列车j及其途径节点v∈∪j∈TNj：xa为0-1变量，xa=1表示弧a为任意列车运行线组成部分，其他为0；yv为0-1变量，yv=1表示节点v被某列车运行线经过，其他为0；zjv为0-1变量，zjv=1表示列车j的运行线经过节点v，其他为0。

对两相邻列车j、k，在车站i∈(Sj{lj})∩(Sk{lk})包含的运行弧(v1,u1)∈Aj，(v2,u2)∈Ak，v1,v2∈Wi，u1,u2∈Ui+1需满足行车组织约束、停站需求约束、天窗约束、列车运行线唯一性约束及相关连接约束。行车组织约束包括：

(1) 列车在i站出发时间差不小于该站出发间隔时间di，即Δ(v1,v2)≥di,v1v2，见图2。

设w为列车运行线可选节点(下同)，相关约束可表示为对任意i∈S{s},v∈Wi，Δ(v1,v2)

( 1 )

(2) 两列车在i站到达时间差不小于该站到达间隔时间ai，即Δ(u1,u2)≥ai,u1u2，见图3。

相关约束可表示为对任意i∈S{1},u∈Ui，Δ(u1,u2)

( 2 )

(3) 两列车运行线在区间内不相交，设v1v2,u2u1，见图4。

设区间内列车k速度高于列车j，相关约束可表示为对相邻车站i和i+1，列车运行线j、k相应节点(v1,u1)∈Aj，(v2,u2)∈Ak，v1,v2∈Wi，u1,u2∈Ui+1的选择限制，即图4中4个节点中至多3个节点被选择

zjv1+zju1+zkv2+zku2≤3

( 3 )

式(3)包括对出发节点及到达节点的选择限制，在区间运行标尺固定情况下，可将其表示为对出发节点的选择限制[10]。

θ(v3)=max[θ(v1)+di,θ(v1)+ai+1+tj-tk]

式中：tj、tk分别为列车j、k在区间(i,i+1) 内的运行时分。

由式( 1 )～式( 3 )可知，v1、v3及v2、v4至多有一个节点被列车j或k的运行线选择，见图5。

同时到达节点可行性可由式( 2 )保证，则区间不相交约束可表示为

( 4 )

基于2.1节所做假设，列车追踪间隔时间可由式( 1 )、式( 2 )保证。

各列车运行线需满足停站需求约束及天窗约束。

停站需求约束可表示为：

(4) 停站需求约束限制最短停时内出发节点选择

( 5 )

天窗约束保证天窗期内任何列车不得从该站出发，见图7。

天窗约束可表示为：

(5) 天窗约束限制v1与v2间出发节点的选择

( 6 )

此外，列车运行线还应满足唯一性约束：

(6) 各列车运行线至多选择1个始发弧

( 7 )

(7) 各列车运行线进入与离开节点v弧数量相等

( 8 )

(8) 各节点至多被一条列车运行线选择

( 9 )

2.3 目标函数

货物列车运行图优化目标通常为减少全程总停时、提高旅行速度[4-5]，其他优化目标可在上述目标下通过合理确定始发时刻实现。本文在各货物列车运行线理想始发时刻给定条件下，以始发时刻调整及列车全程停时总延长为主要考虑因素构建目标函数。

列车运行线由Gj=(Nj,Aj)中弧构成，需满足2.2节所提出的各类约束条件。全图列车运行线总利润最大为目标构建目标函数为

(12)

pa=p(σ,v)+p(u,v)

p(σ,v)=πj-αjυ(v)

p(u,v)=-βjμ(u,v)

υ(v)=|θ(v)-θ(fj)|

式中：pa为列车j运行线中所有构成弧a∈Aj利润和，由构成该运行线的始发弧和停留弧利润之和表示；p(σ,v)为始发弧(σ,v)，v∈Wfj的利润；p(u,v)为停留弧μ(u,v)，u∈Ui,v∈Wi,i∈Sj{fj,lj}的利润；πj为给定列车j的运行线的一个初始利润；αj、βj为相应参数；υ(v)为实际始发点v对应出发时刻与计划出发时刻θ(fj)的偏离；μ(u,v)为第i站实际停时与计划停时的差值。

约束条件为式( 1 )、式( 2 )、式( 4 )～式(11)。

3 模型求解

3.1 约束条件处理

原模型中式( 5 )、式( 6 )与列车特征及编图区段运输组织条件有关，求解过程中通过对列车信息及时空网络的预处理考虑此部分约束，具体做法为：

同时，类似对式( 5 )、式( 6 )的处理方式，将列车运行图中旅客列车运行线占用的节点及各旅客列车运行线引起的不满足式( 1 )、式( 2 )、式( 4 )、式( 9 )、式(11)的相关节点设为禁选。

3.2 原问题的拉格朗日松弛

式(12)及式( 1 )、式( 2 )、式( 4 )～式( 9 )构成了一个0-1整数规划模型，实际问题对应规模庞大，无法在有效时间内精确求解，通常以最优可行解为求解目标[9-11]。

拉格朗日松弛技术通过将影响问题求解的难约束添加拉格朗日乘子后吸收入目标函数，从而减小问题规模、加速求解，广泛应用于大规模整数规划的启发式算法设计[11-12]。基于文献[9]所述方法，对原问题进行拉格朗日松弛加速求解，根据松弛问题所得对偶信息设计启发式算法对原问题最优可行解进行迭代求解。

拉格朗日松弛的关键是确定松弛约束，即模型中影响求解速度的难约束[11-12]。原问题中行车组织约束式( 1 )、式( 2 )、式( 4 )反映列车运行线间相互影响，是一类簇约束。以约束式( 1 )为例，对每一个可能出发节点w，违反约束式( 1 )的所有节点均不可被其他列车运行线选择，即与节点w不兼容，因此约束式( 1 )实际上包含多个约束条件，其数量与di有关。此类簇约束会使模型规模急剧增加，是影响模型求解的主要因素。

由此确定对原模型中的行车组织约束进行松弛。对约束式( 1 )、式( 2 )、式( 4 )分别添加拉格朗日乘子后加入目标函数得原问题的拉格朗日松弛问题

(13)

式中：C为与节点v不相容弧的集合，含义为由相应约束决定的经过与节点v不兼容节点的弧的集合；λC1、λC2为对应的拉格朗日乘子。

对式(13)进行整

(14)

拉格朗日松弛问题中行车组织约束被解除，在拉格朗日乘子一定的条件下其求解是一类简单的最短路问题，可由线性规划方法快速求解[9-11]，对应各列车运行线的解即为其拉格朗日松弛解，其利润为拉格朗日利润。计算全图各列车运行线拉格朗日利润，若大于零则将其加入解集从而求得D(Z)。由于原问题进行拉格朗日松弛后目标函数增加、可行域扩大，因此有D(Z)>Z，即D(Z)为原问题的一个上界[11-12]。

当且仅当拉格朗日乘子全部为0时，D(Z)最优解与Z最优解相等并达到全局最优，此时各列车运行线均为满足约束条件的实际可行解，但该最优解通常无法求得[11]。拉格朗日乘子不为0时，D(Z)中各列车运行线的拉格朗日松弛解不完全符合行车组织约束，可能为实际非可行解。

3.3 基于拉格朗日松弛的启发式算法

D(Z)中包含各节点拉格朗日罚数这一对偶信息，其意义是时空网络中各节点的选择“费用”。当某节点的拉格朗日罚数为0时，表明各列车运行线在该节点不产生潜在占用冲突，即各列车运行线在该节点的选择“费用”为0；当某节点的拉格朗日罚数不为0时，表明多条列车运行线在该节点存在潜在占用冲突，即各列车运行线在该节点的选择“费用”不为0。列车运行线铺画对运输资源占用的潜在冲突就以各节点拉格朗日罚数的形式体现，由此可根据各节点拉格朗日罚数作为启发信息求解各列车运行线实际可行解。

求解列车运行线实际可行解需考虑行车组织约束并确定各列车运行线的铺画顺序。各节点选择“费用”会影响列车运行线的拉格朗日利润，列车运行线拉格朗日利润较高说明其与其他列车运行线的潜在冲突较少，由此根据各列车运行线拉格朗日利润降序确定铺画顺序，优先铺画潜在冲突较少的列车运行线。

在始发时间域内考虑始发时刻调整、总停时延长引起的罚数并选择拉格朗日罚数为0的节点规避运输资源占用冲突，以实际利润最大为目标,根据深度搜索原则遍历具体列车运行线的所有可行解[8]，若该列车运行线实际利润大于零，说明其存在实际可行解并存储所得结果，全图各列车运行线实际可行解之和即为原问题的一个可行解，记为F(Z)。上述求解过程考虑了拉格朗日罚数为0这一节点选择限制，可行域缩小，因此有F(Z)≤Z，即F(Z)是Z的一个下界[11-12]。

D(Z)目标函数值与拉格朗日乘子λh取值有关，可通过次梯度法更新λh优化D(Z)[11-12]，λh的更新可由Fisher提出的公式进行[12]

(15)

式中:UB、LB分别为算法当前最优上界和最优下界，分别对应D(Z)当前最优解及F(Z)当前最优可行解；ρ为给定的非负步长；τ为由当前D(Z)决定的全局次梯度向量；τh为λh对应约束的次梯度向量。

根据更新的拉格朗日乘子求解F(Z)即可对其进行迭代优化[9]。实际问题中，上述算法所得解的质量可由最优上界评估，满足一定条件时可认为已求得最优可行解[12]。

算法开始时，以全图所有列车初始利润之和作为全图总利润并存储为当前最优上界以开始迭代，上述启发式算法的终止条件设置为：(1)迭代次数达到设定上限；(2)算法运行时间达到设定上限；(3)当前最优上界与当前最优下界间差值小于设定值；(4)各拉格朗日乘子为0。算法流程见图8。

4 实例验证

以京九线北京西至阜阳区段实际数据为例进行验证。该区段包括55个车站，全长855 km，根据2015年5月编制的列车运行图数据，全线共开行列车711列，其中货物列车439列，列车构成见表1。

表1 列车构成

“*”表示客车标尺，归类为旅客列车。

根据实际数据设定各货物列车经停站及停时。以添加起停车附加时分的松弛运行标尺作为列车区间运行时分验证算法有效性，保证结果可行性的同时增加列车运行线冗余度。需说明的是上述松弛运行标尺可能会影响旅行速度，但实验结果显示原图货物列车平均旅行速度仍得到了不同程度改善，证明了算法有效性。

视各列车重要程度相同并设定初始利润为10 000、总停时最大延长为210 min，以列车实际始发时刻为初始布点、始发时间域为初始布点的±20 min，次梯度法求解拉格朗日乘子时起始步长ρ为1，若迭代10次后最优上界未发生改变则减半ρ值，当前最优上界与当前最优下界差值设为当前最优上下界差值的0.1%。

由2.3节知，算法计算与参数αj、βj有关。αj、βj除表示列车始发时刻偏离及总停时延长造成的罚数外，另一意义是列车始发时刻及全程停时总延长重要程度的界定，因此可根据优化需要调整参数αj、βj以不同策略对列车运行图进行优化编制。本文采取3种参数设置策略分6组进行测试并分析所得结果，其中1、2组为始发时刻优先策略，3、4组为旅行速度优先策略，5、6组为均衡策略。

列车运行图优化编制对时效性要求不高，为得到更好的结果，设置算法运行时间上限为12 h、迭代次数上限为500次，相关算法基于C语言编程实现，实验并在一台Intel Core2 2.33 GHz，2 GB内存的个人计算机上进行。

4.1 始发时刻优先策略

设置αj>βj优先选择列车运行线可行方案中始发时刻总调整值较小的解。统计最优上界、最优下界、全图货物列车始发时刻总调整值、全程停时总延长值、算法迭代次数、计算时间、优化后平均旅行速度、旅行速度改善比率等指标见表2，算法因迭代次数达到上限而终止。

不同参数下算法上界、下界及约束数量随迭代次数变化曲线见图9。

表2 始发时刻优先策略指标统计

列车种类参数设置 αjβj算法指标优化后平均旅行速度/(km·h-1)旅行速度改善比率/%技术直达始发直达煤炭、石油直达直通货物区段货物摘挂及小运转其他全图平均101最优下界7093799.039.590.03最优上界7107354.558.9821.91始发时刻总调整值/min39935.6315.98全程停时总延长值/min1221140.9412.28迭代次数∗50050.138.18计算时间/s3736141.7740.8340.440.0041.956.08技术直达始发直达煤炭、石油直达直通货物区段货物摘挂及小运转其他全图平均201最优下界7099557.039.56-0.05最优上界7102647.559.1422.24始发时刻总调整值/min38135.6315.98全程停时总延长值/min1282340.9112.21迭代次数∗50049.847.55计算时间/s3701141.6940.5640.440.0041.816.44

“*”表示算法因上下界差值达到设定标准而终止。

4.2 旅行速度优先策略

设置αj<βj优先选择列车运行线可行方案中总停时延长较小的解，结果见表3。

不同参数下算法上界、下界及约束数量随迭代次数变化曲线见图10。

表3 旅行速度优先策略指标统计

列车种类参数设置 αjβj算法指标优化后平均旅行速度/(km·h-1)旅行速度改善比率/%技术直达始发直达煤炭、石油直达直通货物区段货物摘挂及小运转其他全图平均110最优下界∗7021310.040.672.75最优上界7021308.559.8923.79始发时刻总调整值/min239037.3321.52全程停时总延长值/min963043.1318.29迭代次数5653.4215.28计算时间/s458743.2345.7540.650.4943.7211.30技术直达始发直达煤炭、石油直达直通货物区段货物摘挂及小运转其他全图平均120最优下界∗6922247.040.722.28最优上界6922243.060.5425.13始发时刻总调整值/min243337.4021.74全程停时总延长值/min926643.3218.82迭代次数2853.1014.59计算时间/s216043.2245.7240.650.5243.7211.30

“*”表示算法因上下界差值达到设定标准而终止。

4.3 均衡策略

设置αj=βj忽略优先级差异均衡选择列车运行线方案，结果见表4。

不同参数下算法上界、下界及约束数量随迭代次数变化曲线见图11。

表4 均衡策略全图指标统计

“*”表示算法因上下界差值达到设定标准而终止。

4.4 结果分析

实验结果显示，与原图相比拉格朗日松弛启发式算法所得各类货物列车旅行速度得到了不同程度的优化，全图货物列车平均旅行速度由39.28 km/h提升至41.81～43.72 km/h(视优化策略而定)。

算法性能受列车运行图编制策略的影响，具体体现在：

(1) 不同策略下全图列车始发时刻总调整值、全程停时总延长值不同。始发时刻优先策略下，全图始发时刻总调整值较小；旅行速度优先策略下全图停时总延长值较小；均衡策略介于二者之间，同时相同策略下算法对参数不敏感，见图12。

(2) 不同策略下全图总利润(最优下界)、列车平均旅行速度不同。全图总利润最大值为7 099 557.0，最小值为6 882 300.0，分别在旅行速度优先策略、均衡策略下达到；列车平均旅行速度最高为43.72 km/h，最低为41.81 km/h，分别在旅行速度优先策略、始发时刻优先策略下达到，见图13。