闫 璐,张 琦,丁舒忻,王荣笙,
(1. 中国铁道科学研究院研究生部,北京 100081;2. 中国铁道科学研究院集团有限公司通信信号研究所,北京 100081;3. 中国铁道科学研究院集团有限公司国家铁路智能运输系统工程技术研究中心,北京 100081)
高速铁路列车运行存在“高密度”的特征[1],一旦受到突发事件影响,往往会引发列车晚点,并在线路和路网中传播,情况严重时会给运输组织和旅客正常出行造成极大影响。根据晚点的原因和程度,可将引发晚点的内外界扰动分为2 种类型:一种是一般干扰,如旅客乘降作业超时等;另一种是严重干扰,如因风、雨、雪等自然灾害或设备故障产生的长时间封锁等[2]。无论哪种干扰,均将导致列车运行在一定程度上偏离既定计划的后果,此时就需要进行列车运行调整。列车运行调整是实现“按图行车”的关键,也是高速铁路列车安全、高效、有序、正点运行的基础[3]。
国内外众多学者已对列车运行调整问题开展了大量研究[4-7]。现阶段的研究通常考虑多个优化目标,其中减少列车总晚点时间多为主要目标,其他优化目标包括提高旅客满意度、减少取消列车数量等。这些需要优化的目标之间往往互相冲突,需要权衡决策者对不同优化目标的偏好程度,从而得到不同的帕累托(Pareto)最优解。目前大部分研究采用前决策(a priori)的方式[8-14],考虑决策者已经提供了全局偏好信息,在此基础上搜寻该偏好下的一个Pareto 最优解。LI 等[8]针对车站股道封锁不确定恢复时间,考虑最小化晚点成本和股道调整成本,建立混合整数非线性规划模型,采用2 阶段法求解,第1阶段通过贪心算法求解列车调整后的到发时间和次序,并在该决策方案下采用遗传算法调整车站股道,完成第2 阶段求解。文献[9—13]均考虑最小化取消列车趟数和总晚点时间,通过商用求解器CPLEX 或GUROBI 求解。WANG 等[14]针对初始晚点,以最小化列车总晚点时间和严重晚点列车趟数为优化目标,建立混合整数规划模型,并提出遗传算法结合粒子群算法的求解方法。目前大部分文献都采用对不同优化目标线性加权的形式确定决策者偏好信息,这样虽然不需要大量计算,但决策者的全局偏好信息并不能准确获得。
还有部分研究采用后决策(a posteriori)方式,即先搜索整个Pareto前沿,再在该前沿中选择最偏好的解。该方法并不需要知道决策者的全局偏好信息,而且可找到整个Pareto 前沿的近似解集,但需要很大的计算代价。SHAKIBAYIFAR 等[15]针对区间封锁,以最小化列车终点站晚点时间和最小化列车发车运行偏差建立双目标优化模型,采用多目标邻域搜索算法求解。BINDER 等[16]针对车站存在多个股道封锁,考虑最小化旅客不满意度(换乘时间)、运营成本和调整成本,建立整数规划模型,通过CPLEX 结合ε-约束法进行求解。ALTAZIN 等[17]针对列车晚点,考虑最小化恢复时间、旅客不满意度(等待时间和车内时间)、调整操作(如取消车次、列车折返、增加停站、跳停等)总数、总晚点时间和晚点事件个数,采用启发式方法求解并通过仿真进行评价。然而上述研究并未计算整个Pareto前沿,缺失了部分非支配解。
针对高速铁路列车运行中的晚点情况,综合考虑区间运行干扰和车站停车作业干扰,通过调整列车到发时刻和列车次序,以最小化列车总晚点时间和列车到发时刻调整次数为双优化目标,建立高速铁路列车运行调整的混合整数非线性规划模型,并通过定义辅助矩阵将其转化为线性规划模型;提出改进的ε-约束法求解模型。通过算例以及加权法和先到先服务方法对比,验证模型和求解算法。
以最小化列车总晚点时间和最小化列车到发时刻调整次数为双优化目标,以车站作业、区间运行、追踪间隔等时间条件为约束条件,构建高速铁路列车运行调整的双目标优化模型。
结合实际情况首先做如下假设:
(1)列车运行调整的方式包括调整列车到发时刻和列车次序;
(2)突发事件最终表征为列车在区间运行时和在车站停车作业时受到的干扰;
(3)车站均满足接发车能力限制。
定义模型参数:i,l为列车编号;N为列车总数;j为车站编号;J为车站总数;k为区间编号;K为区间总数;τsij为列车i在车站j的图定到站时刻;τeij为列车i在车站j的图定发车时刻;αi为列车i的始发站车站编号;βi为列车i的终到或交出站车站编号;为列车i在车站j的最小作业时间(停站时间);为列车i在区间k的最小运行时间;Ik为区间k的最小追踪间隔时间;M为足够大的正整数;为列车i在区间k受到干扰时增加的区间运行时间(区间运行干扰时间);为列车i在车站j受到干扰时增加的车站停车作业时间(车站停车作业干扰时间);为含有区间运行干扰的列车及对应的区间集合;为含有车站停车作业干扰的列车及对应的车站集合。
定义决策变量:xsij为列车i在车站j的实际到站时刻;xeij为列车i在车站j的实际发车时刻;qilk为0-1 决策变量,表示列车i和列车l在区间k的实际次序,当列车i先于列车l进入区间k时取值为1,否则取值为0。
受到晚点影响后调度部门需要进行列车运行调整,使列车尽快恢复正点运行。对于受到影响的列车来说,在不同车站调整到发时间,会影响车站当前作业,增加相关作业人员负担。因此,以最小化列车总晚点时间和最小化列车到发时刻调整次数为目标函数。
1)列车总晚点时间最少
定义列车总晚点时间等于运行调整之后所有列车到站时刻与图定到站时刻的差值加上列车发车时刻与图定发车时刻的差值,则列车总晚点时间最少的目标函数f1(x)可表示为
其中,
式中:x为实际到站时刻和实际发车时刻。
2)列车到发时刻调整次数最少
定义列车到发时刻调整次数等于所有列车到站和发车的总晚点次数,即晚点1次就需要调整1次,则列车到发时刻调整次数最少的目标函数f2(x)可表示为
式中:sgn(·)为符号函数,返回对应参数的正负号(参数为0 则返回值为0),由于调整之后列车到站或发车时刻不允许早于图定时刻,返回值不为负数,若调整后列车到站或发车时刻晚于图定时刻,则返回值为1,记录晚点1 次;若等于图定时刻,则返回值为0,不记录晚点次数。
为了保证列车运行安全,合理利用接发列车能力、车站通过能力以及区间通过能力,建立的模型应满足如下约束条件。
1)车站停站时间约束
列车i在车站j存在停车作业干扰时对运行图的影响如图1所示。图中:黑色实线表示计划运行线,对应停站时间为图定停站时间;红色虚线表示最小停站时间下的运行线;蓝色虚线表示含有车站停车作业干扰时间下的运行线,此时列车停站时间在图定停站时间基础上增加干扰时间,该干扰时间为车站停站缓冲时间吸收之后的干扰时间。
图1 车站停车作业存在干扰时对运行图的影响
对于含有车站停车作业干扰的列车,调整之后的列车停站时间必须不小于该列车在该站的图定停站时间与车站停车作业干扰时间之和;对于不含车站停车作业干扰的列车,调整之后的列车停站时间必须不小于该列车在该站的最小停站时间,即
式中:(i,j)∈Tdwell,dis和(i,j)∉Tdwell,dis分别为含有、不含车站停车作业干扰的列车及对应的车站。
2)区间运行时间约束
列车i在区间k存在区间运行干扰时对运行图的影响如图2所示。图中:黑色实线表示计划运行线,对应区间运行时间为图定运行时间;红色虚线表示最小区间运行时间下的运行线;蓝色虚线表示含有区间运行干扰时间下的运行线,此时列车区间运行时间在图定运行时间基础上增加的干扰时间,即为区间运行缓冲时间吸收之后的干扰时间。
图2 区间运行时间存在干扰时对运行图的影响
对于含有区间运行干扰的列车,调整之后的列车区间运行时间必须不小于该列车在该区间的图定运行时间与区间运行干扰时间之和;对于不含区间运行干扰的列车,调整之后的列车区间运行时间必须不小于该区间的最小运行时间,即
式中:(i,k)∈Trun,dis和(i,k)∉Trun,dis分别指含有、不含区间运行干扰的列车及对应的区间。
3)列车到发间隔时间约束
对于在同一区间内运行的相邻列车,其到站和发车间隔时间必须不小于最小追踪间隔时间,即
式中:∨和∧分别表示“或”和“且”,用于计算相邻列车始发车站编号最大值和相邻列车终到或交出车站编号最小值;qilk与qlik均为列车在区间次序的决策变量,二者取值不同,即qilk=1 时qlik=0,qilk=0时qlik=1。
4)列车到发时刻约束
列车在车站实际发车时刻必须不小于图定发车时刻与停车作业干扰时间之和,实际到站时刻必须不小于图定到站时刻与区间运行干扰时间之和,即
5)决策变量约束
列车实际到发时刻决策变量必须为非负变量,列车在区间的次序决策变量必须为0-1变量,即
式(2)中存在sgn(·)符号函数,因此需要将其转换为线性模型再进行处理。定义辅助矩阵c1=[c1,ij]N×J和c2=[c2,ij]N×J,赋值为
通过式(11)将式(2)中的参数替换,使后者转换为式(12)所示的混合整数线性规划模型。由此建立的列车运行调整双目标优化模型P0为
该模型属于NP-hard问题[3],无法在多项式时间内找到整个Pareto前沿,需要将其转换为单目标优化模型(例如加权法或ε-约束法),通过商用求解器(如CPLEX,GUROBI 等)对一定规模下的问题在合理时间内得到Pareto前沿。
参考文献[18],针对本文基于双目标优化的高速铁路列车运行调整问题,引入如下定义。
定义1(Pareto 支配):目标向量u和v为2 个可行解集对应的优化目标,目标向量u支配v(记作u≻v),当且仅当ub≤vb,∀b∈{1,2,…,B}(B为目标个数),且u≠v。
定义2(Pareto 有效性):对于可行解x,不存在其他可行解y,其目标函数值构成的目标向量f(x)和f(y)满足Pareto支配关系f(y)≻f(x),则说明解x为Pareto有效性解,又称为非支配解。
定义3(Pareto 解集,Pareto Set):所有满足Pareto 有效性的解的集合称为Pareto 解集,又称为非支配解集,记为APS。
定义4(Pareto 前沿,Pareto Front):所有Pa⁃reto 最优目标向量的集合称为Pareto 前沿,记为APF={f(x)|x∈APS}。
定义5(理想点,Ideal Point):向量zIdeal=(zIdeal1,…,zIdealB)中每个元素都为每个目标函数fb(x)的最优解,可以通过zIdealb=minfb(x)求得理想点。
定义6(最差点,Nadir Point):向量zNadir=(zNadir1,…,zNadirB)中的每个元素都为每个目标函数fb(x)在Pareto 前沿上的最大值。其中每个元素zNadirb定义为zNadirb=max {fb(x)|x∈APF}。
求解多目标优化问题的方法包括ε-约束法和线性加权法,其原理都是通过将原问题转换为多个不同的单目标优化问题,从而得到整个Pareto 前沿。其中,ε-约束法是通过将原有多个优化目标分为主目标和其他的次目标,并将其他的次目标作为约束进行单目标优化求解。相比于线性加权法,ε-约束法不需要设置权重,无需对目标函数进行归一化处理,并能够得到非凸前沿。
对模型P0,定义ε2∈[zIdeal2,zNadir2]为次目标优化上界,采用ε-约束法将目标函数f2(x)转化为约束条件,则可将其转化为模型P1,即
通过修改约束右边的上限ε2,可以得到对应的Pareto前沿。
然而ε-约束法还存在以下3 个缺陷[19]:①在Pareto 前沿之外存在一些不必要的计算;②无法保证得到所有位于Pareto前沿上的解;③当优化目标超过2个时会显著增加求解时间。因此,需要对ε-约束法进行改进。目前主流的改进方法如增广ε-约束法(包括AUGMECON[19]和AUGMECON2[20]),虽然通过引入松弛变量、确定理想点和最差点等形式确定了约束范围,但求解时增加了额外变量,从求解效率来看优势不足。有必要进一步改进ε-约束法,在克服上述3个缺陷的基础上实现高效求解。
分析式(2)可知,优化目标f2(x)均为整数情况,可利用理想点和最差点估计出Pareto前沿的最大个数。利用每次求解约束优化问题得到对应的次优化目标值f2(x∗),便能够得到下一次优化的约束上界。相比于传统ε-约束法固定约束增量下的求解,这种方法能够减少不必要的重复非支配解的计算;相比于增广ε-约束法,这种方法不仅可减少对辅助变量的使用,还能够在无须为设置均匀分布的格点的基础上保证求解到整个Pareto前沿。
采用改进ε-约束法求解模型P1的步骤如下。
步骤1:计算理想点,求解不含式(12)的优化模型P0,得到zIdeal1;求解不含式(1)的优化模型P0,得到zIdeal2。
步骤2:计算最差点,求解不含式(12)的优化模型P0,求解时增加额外约束f2(x)=zIdeal2,得到zNadir1;求解不含式(1)的优化模型P0,求解时增加额外约束f1(x)=zIdeal1,得到zNadir2。
步骤3:确定次优化目标取值范围,设定约束为次优化目标最大值;根据理想点和最差点得到目标函数f2(x)的范围为[zIdeal2,zNadir2],求解最多(zNadir2-zIdeal2+1)个约束单目标模型P1 得到对应Pareto前沿,并令ε2=zNadir2,得到对应的模型P1。
步骤4:在设置的约束下计算调整方案,求解模型P1 得到当前最优解x∗,记录解[f1(x∗),f2(x∗)]。
步骤5:更新约束值,令ε2=f2(x∗)-1。
步骤6:判断约束值是否为次优化目标最小值,如果ε2≠zIdeal2,则转步骤4;如果ε2=zIdeal2,则继续步骤7。
步骤7:输出Pareto 前沿,将步骤4 中不同约束下求解模型P1 得到的解合并,去除其中的支配解,最终构成了模型P0对应的Pareto前沿。
此外,若决策者仅对部分Pareto 前沿感兴趣,可以在目标函数f2(x)的范围[zIdeal2,zNadir2]内选择优化求解部分模型P1,得到部分感兴趣的非支配解。
选取京沪高铁北京南—泰安区段为例,设置3 种干扰场景,采用本文模型和改进ε-约束求解算法进行列车运行优化调整;并与同样采用本文模型但模型求解时分别采用加权法和先到先服务(First-Come-First-Served,FCFS)的启发式策略(简称为FCFS 法)得到的结果进行比较,验证本文模型和改进ε-约束法的合理性、有效性。
北京南—泰安区段共有7 个车站,自北京南开始依次编号为1,2,…,7。6 个区间的最小运行时间见表1。车站最小停站时间为2 min,相邻列车最小追踪间隔为4 min,M取值为1 000。下行开行列车40趟。
表1 区间最小运行时间tmin,run ik
根据不同的干扰情况(区间运行干扰、车站停车作业干扰),设置以下3个场景。
场景1:仅由车站停车作业干扰组成。设置列车2在北京南停车作业干扰时间为20 min,列车20在廊坊停车作业干扰时间为20 min,列车30 在北京南停车作业干扰时间20 min。
场景2:仅由区间运行干扰组成。设置列车4在北京南—廊坊区间运行干扰时间为15 min,列车18在廊坊—天津南区间运行干扰时间为20 min,列车32在北京南—廊坊区间运行干扰时间为20 min。
场景3:由上述2种干扰共同组成。设置列车3在北京南停车作业干扰时间20 min,列车25在廊坊停车作业干扰时间10 min,列车33在北京南停车作业干扰时间15 min;列车6 在北京南—廊坊区间运行干扰时间为15 min,列车15在廊坊—天津南区间运行干扰时间为20 min,列车28在北京南—廊坊区间运行干扰时间为20 min。
对于上述场景中未提及的列车及其途经区间和车站,均为正常场景,即其对应的区间运行干扰时间和车站停车作业干扰时间均为0。
为了验证本文提出的改进ε-约束法的性能,采用了下面经典的多目标优化性能指标。
非支配解个数(Number of non-dominated so⁃lutions,NNS):该性能指标描述了算法得到的近似Pareto前沿,NNS值越大说明该算法求解能力越强。
逆代距(Inverted generational distance,IGD)[21]:该指标综合评价了非支配解集的收敛性和多样性,IGD值越小,说明非支配解集A的质量越好。假设Ψ∗为1个均匀分布在Pareto真实前沿的集合,Ψ为1个非支配解集,则Ψ的IGD值IIGD定义为
式中:v为解集Ψ∗中的1 个解;d(v,Ψ)为v和Ψ中所有点的最短欧氏距离;|Ψ∗|为集合Ψ∗的势。
超体积(Hypervolume,HV)[21]:该指标通过非支配解和参考点,计算出超体积。该性能指标同时评价了非支配解集的收敛性和多样性。HV 的值越大,对应算法的性能越好。对于双目标优化问题中的一个非支配解集A,其HV值IHV计算如下
式中:xi为非支配解集A中的1个解;Vi为解xi对应的目标值f(xi)和参考点为边界构成的矩形体积。
针对3.1 节中给出的3 种干扰场景,在Intel Core i5-8265U CPU 1.60GHz,8GB 内存,操作系统Windows 10,64 位主机上分别采用改进ε-约束法、加权法和FCFS 法求解本文模型。其中改进ε-约束法和加权法均采用商用求解器GUROBI 9.1.0,通过YALMIP 工具包[22]在Matlab R2018b进行仿真求解。GUROBI 各参数采用默认值。加权法以w1f1(x)+(1-w1)f2(x)为优化目标,其中w1取范围[0,1]之间的均匀间隔0.02 为权重,权重组合个数为51。
3.3.1 求解得到的Pareto前沿和解
3 种方法得到3 种场景下的Pareto 前沿和解如图3所示。由图3可以得到如下结论。
图3 不同场景下的非支配解的分布
(1)采用改进ε-约束法求解,得到了3 种场景下的全部Pareto前沿。
(2)采用加权法求解,在场景1 中的前沿为Pareto 前沿的一部分(共7 个解),但无法得到对应的非凸前沿(共12 个解),即调整次数为77,64,63,60,57,52,51,50,49,43,42 和40的解。而在场景2 和场景3 中,仅得到了左上角的部分Pareto前沿,而右下角得到的是支配解。这些实验结果说明处理多目标优化问题中,加权法存在一些缺点。分别为:无法得到非凸前沿;受不同优化目标函数值的范围影响,需要对目标函数进行归一化来避免加权后的单目标优化问题过于偏向优化某1 个目标函数,造成重复计算;加权法在边界值时由于仅优化单个目标,没有办法保证另1个优化目标同时最小,会出现支配解。而本文提出的改进ε-约束法能够克服这些缺点。
(3)采用FCFS 法求解,该策略以先到先服务的启发式规则完成运行计划调整,没有专门对不同的优化目标进行处理,因此仅能得到1个解。对于场景1 中干扰类型均为车站停车作业干扰的情况,FCFS 法得到的解并不理想,与Pareto 前沿的差距较大。而对于场景2 和场景3 中含有区间运行干扰的情况,FCFS法相对效果有所提升。
3.3.2 调整后的列车运行图
1)场景1
该场景下的列车计划运行图如图4所示,不同车次的运行线通过粗细不同区分。运用改进ε-约束法(优化目标为(629,67))和FCFS 法(优化目标为(927,76))进行列车运行调整后,得到的运行图分别如图5和图6所示,图中红色虚线表示运行调整后与原计划运行图不一致的列车运行线。由于加权法得到的结果大部分和改进ε-约束法相同,这里不做具体分析。由图4—图6可以看出:对于车站停车作业干扰的情况,改进ε-约束法会更多地通过调整受影响的列车次序完成调整,而FCFS 法依然按照干扰之后的列车次序调整到发时刻;针对第1 个优化目标,FCFS 法得到的结果是改进ε-约束法的1.47 倍,总晚点时间多298 min;在终点站时,改进ε-约束法下仅有1 趟列车存在晚点,其他受影响的列车在到达终点站之前均通过调整实现了恢复,而FCFS 法下还存在6 趟车到达终点站晚点情况。
图5 场景1下运用改进ε-约束法按优化目标(629,67)调整后的列车运行图
图6 场景1下运用FCFS法按优化目标(927,76)调整后的列车运行图
2)场景2
该场景下,运用改进ε-约束法(优化目标为(843,86))和FCFS 法(优化目标为(857,87))进行列车运行调整后,得到的运行图分别如图7和图8所示。由图7和图8可以看出:改进ε-约束法和FCFS 法对于场景2 前2 个干扰影响到的车次,均未调整发车次序,而是仍按原有次序依次调整到发时刻,仅在第3个干扰下分别对列车次序进行了调整,这是因为场景2 中的干扰均属于区间运行干扰的缘故;对于第1 个优化目标列车总晚点时间,FCFS法得到的结果是改进ε-约束法的1.02倍,总晚点时间比后者多14 min;对于运行至终点站仍存在晚点的列车数量,改进ε-约束法下为2 列,而FCFS法为3列。
图7 场景2下运用改进ε-约束法按优化目标(843,86)调整后的列车运行图
图8 场景2下运用FCFS法按优化目标(857,87)调整后的列车运行图
3)场景3
该场景下,运用改进ε-约束法(优化目标为(1 180,115))和FCFS 法(优化目标为(1 415,115))进行列车运行调整后,得到的运行图分别入图9和图10所示。从图9和图10可以看出:在存在多种类型干扰的场景3 下,尽管对于区间运行干扰下FCFS 法和改进ε-约束法的调整方案相似,但由于FCFS 法没有调整次序的能力,其无法得到全局最优解;在相同的第2 个优化目标下,对于第1 个优化目标列车总晚点时间,FCFS 法是改进ε-约束法的1.20 倍,总晚点时间比后者多235 min;对于运行至终点站仍存在晚点的列车数量,改进ε-约束法为4列,而FCFS法为7列。
图9 场景3下运用改进ε-约束法按优化目标(1 180,115)调整后的列车运行图
图10 场景3下运用FCFS法按优化目标(1 415,115)调整后的列车运行图
3.3.3 3种求解算法的性能指标比较
3 种场景下3 种求解算法的性能指标见表2。此外,表中还给出了算法的总计算时间和对应得到非支配解的平均计算时间,其中加粗数据表示该算法对应性能指标最优。从表2中可以看出:改进ε-约束法在NNS,IGD 和HV 这3 种性能指标中均优于加权法和FCFS 法,且得到了整个Pareto 前沿,但由于算出了所有的非支配解,所以其总计算时间较长,即在场景2 和场景3 中大于加权法和FCFS法;加权法并没有算得全部Pareto前沿,还存在一些不必要的计算,在场景1中总时间最长;对于非支配解的平均计算时间,加权法同样大于改进ε-约束法;FCFS 法作为一种启发式策略,其计算时间会明显优于其他基于求解器的精确算法,但其计算结果并不能保证最优。
表2 不同场景下的算法性能
综上,本文提出的改进ε-约束法,能够有效求解提出的基于双目标优化的高速铁路列车运行调整模型的整个Pareto前沿。对于其计算时间在部分场景中较长的问题,实际中调度员可以有选择的控制改进ε-约束法中的次优化目标对应的约束范围,控制得到的非支配解的数量和偏好,降低调度员对不感兴趣区域搜索所需要的时间。
(1)以列车运行调整的列车总晚点时间和列车到发时刻调整次数2 个优化目标,以车站作业、区间运行、追踪间隔等时间为约束条件,构建高速铁路列车运行调整的双目标优化模型;对模型中的非线性项进行线性化处理;提出的改进ε-约束法对模型求解。
(2)提出的改进ε-约束法,得到了包括非凸前沿的整个Pareto前沿,可以为调度部门提供不同的列车调整决策方案,并且在非支配解个数NNS,逆代距IGD 和超体积HV 这3 个优化性能评价指标上均优于加权法和FCFS 法。与加权法相比,大部分非支配解对应列车总晚点时间和列车到发时刻调整次数相同,但改进ε-约束法的计算时间更短;对于边界值,加权法可能会得到支配解,但改进ε-约束法对应调整后的总晚点时间更短。与FCFS法相比,改进ε-约束法能够得到多样化的结果,在列车到发时刻调整次数相近情况下,总晚点时间更短,终点站晚点列车个数更少。
(3)当更加复杂的干扰场景,求解问题规模增大后,计算整个Pareto前沿会消耗过多时间,可以在Pareto前沿中选择调度员感兴趣的搜索方向,搜索部分非支配解。此外,还可以通过设计一些启发式方法和多目标进化计算算法,在有限时间获得收敛性、多样性好的近似前沿。