需求不确定的铁路车流径路优化模型

Optimization Model of Railway Wagon Flow Path with Uncertain Demand

摘要

引入不确定理论中的鲁棒优化描述铁路货运日常运输组织工作中的车流量波动，提出需求不确定的铁路车流径路优化模型。首先，考虑线路通过能力和网络流量守衡约束，以绕道率阈值参数控制车流径路里程的合理性，以走行成本最小化为优化目标，建立需求确定的车流径路优化模型；然后，引入盒式不确定集描述车流量的波动性，设计鲁棒代价调节模型的保守程度，基于鲁棒对等理论建立需求不确定的车流径路优化模型；最后，通过小型算例验证模型的正确性，并基于沈阳铁路局集团有限公司所辖路网的实际算例验证模型的实用性。结果表明：部分车流量发生波动时，对应的最优径路方案随之发生变化；鲁棒代价能够权衡最优方案的成本和违反通过能力约束的概率，小型算例和实际算例中鲁棒代价分别取5和11时，得到的车流径路方案处于系统成本与运输需求满足的平衡点。该模型及鲁棒优化理论可为管理者制定车流径路方案提供决策参考。

Abstract

The robust optimization in the uncertainty theory is introduced to describe the traffic flow fluctuation in the daily transportation organization work of railway freight transport, and an optimization model of railway wagon flow path with uncertain demand is proposed. First, considering the constraints of line carrying capacity and network flow balance, setting a detour rate threshold parameter to measure the rationality of the wagon flow path mileage, and regarding the minimization of the travel cost as the optimization target, the optimization model of wagon flow path with the fixed demand is established. Then, the traffic flow fluctuation is described by introducing the box uncertainty set, the conservatism of the model is adjusted by designing the price of robustness, and the optimization model of wagon flow path with uncertain demand is constructed based on robust counterpart theory. Finally, the validity of the model is verified by a small-scale example, and the practicability of the model is verified by a real freight transportation network in Shenyang Railway Bureau Group Co., Ltd. The results show that the corresponding optimal path plan changes when part of the traffic flow fluctuates. The price of robustness can balance the optimal plan cost and the probability of violating the carrying capacity constraint. When the price of robustness is 5 and 11 in the small-scale example and real example respectively, the obtained wagon flow path plan is at the balance point between the system cost and the transportation demand satisfaction. This model and robust optimization theory can provide decision-making reference for managers to make the appropriate plan of wagon flow path.

Graphical abstract

关键词

铁路货运 / 通过能力 / 车流径路 / 运输需求波动 / 不确定理论 / 鲁棒代价

Key words

Railway freight transport / Carrying capacity / Wagon flow path / Fluctuation of transportation demand / Uncertainty theory / Price of robustness

铁路车流径路是编制货物列车编组计划最主要的基础和依据之一。在径路优化问题中，比较典型的代表是车辆径路优化问题（Vehicle Routing Problem，VRP）。最早的VRP问题可以追溯到1954年由Dantzig等^［1］提出的较大规模的旅行商问题（Traveling-Salesman Problem，TSP）。此后，TSP作为VRP问题的一种特殊形式得到了广大学者的研究和优化，而TSP本身又属于理论与算法设计一直备受关注的典型NP-hard组合优化问题^［2-4］。在关于运输网络的研究中，径路优化问题通常需要满足相应约束如网络流量守恒约束、运输能力约束和运行时间约束等^［5-8］。铁路运输行业由于经营模式的特殊性，其径路优化理论与其他运输方式相比有自己的特点，如为了达到运输规模经济的目的，车站间产生的原始流量和后方相邻技术站产生的到站相同的车流，需要在某一支点站集结后编入同一编组去向后发出，这也就形成了铁路运输独有的树形径路。我国对于铁路车流径路优化问题已有较为系统化的研究，通常以运输成本、运输距离或走行费用为优化目标建立数学规划模型。如，刘畅^［9］以公铁两网综合运输成本最小化为目标，构建了公铁货流转移情形下的铁路车流径路联合优化模型。何世伟等^［10］以车流走行费用最小化和路网点线能力协调为目标，考虑技术站解编作业能力、路网线路能力和车流改编次数约束，建立了铁路车流径路优化模型。车流径路优化问题属于NP-完全问题，对于小规模的案例可以采用商业软件进行求解，如Lingo和CPLEX等；对于具有实际规模的案例，大部分学者采用了启发式算法。如，薛锋等^［11］基于传统的车流分配及径路优化模型，根据罚函数的形式改进模型，并采用惩罚项动态更新的策略来回避传统固定惩罚因子难以取合适值的缺陷。

总体来看，关于铁路车流径路优化和编组计划优化问题的研究大多是在参数确定的情况下展开的，没有考虑影响编组计划编制的关键参数可能会在计划实施过程中发生变化的情况，如车流量会随着时间发生波动、集结系数会受车组大小的不均衡性和集结中断次数与时间的影响，但这些参数在优化过程中常取为定值，使确定环境下得到的车流径路和编组计划无法在实际应用中到达最优效果。有部分学者采用随机规划或模糊规划的方法来描述和解决这种参数不确定的优化问题。如李宵寅等^［12］将车流量、技术站编组去向数、技术站改编车数、线路区间通过列车数等的波动性分别考虑到列车编组计划的优化问题中，利用模糊机会约束规划理论对参数的波动性进行数学处理。陈崇双等^［13］在讨论固定重量形式分组列车的总集结消耗时，假设到达批中的车辆数服从泊松分布及间隔时间服从指数分布。严余松等^［14］假设车流量波动服从一定的概率分布，设计了基于随机模拟的改进分支定界法对模型进行求解，结果发现随机参数的模拟规模越大，最优编组计划与车流径路方案稳定性越好。随机规划需要掌握不确定参数的具体分布规律，但由于数据采集过程失误或不可抗力因素干扰，会导致分布规律很难获得，进而导致随机规划在编组计划与车流径路优化中的应用受到限制；模糊规划需要依据决策者的个人经验给出不确定参数的模糊隶属度函数，但这往往会导致模糊规划存在较大的主观性，实际应用时可能需要经过反复多次调整。

鲁棒优化是针对传统不确定优化方法不足，由鲁棒控制理论发展而来的一种方法。相比于传统方法，鲁棒优化不需要确定不确定参数的分布模型或模糊隶属度函数，且鲁棒优化的约束条件严格成立，最优解对参数变化的敏感性较低。Gabrel等^［15］在不确定性需求分布未知的情况下，对这种不确定性采用盒式不确定集的表示方法，并引入鲁棒参数对发生波动的客户需求数量进行限制。Wouter等^［16］对鲁棒时间窗约束下的库存-路径问题展开研究，同样采用盒式不确定集对不确定性参数进行表示。

为缓解铁路网络运输能力与增长的货运量需求之间不匹配的局面，增强车流径路优化模型在实际生产中的应用性，本文针对铁路车流径路优化中涉及的车流量关键因素，考虑线路通过能力，对车流径路鲁棒优化问题展开研究。

1 需求确定的铁路车流径路优化模型

在车流量确定的情况下（即采用计划车流），铁路车流径路优化模型考虑线路通过能力和运输径路里程限制，以最小化运输成本为优化目标，以便给出符合能力约束和里程约束的最佳车流径路方案。

1.1　模型假设

（1）将线路最大通过能力定义为扣除旅客列车占用能力后的剩余通过能力。为了适应日常货流波动，线路需要保留一定的后备能力，模型中涉及的通过能力是指在最大通过能力的基础上，按照通过能力后备系数再扣除一部分能力。该系数可根据各铁路方向的具体情况进行取值，一般单线取20%，双线取15%。

（2）将货运需求定义为日均车流且不考虑货物品类。尽管实际中的车流量在不同时期是变化的，如我国供暖季一般是当年11月至来年3月，铁路煤炭运输也会随之出现阶段性运力偏紧的状态。运输需求的这种波动性正是本文重要的研究背景，但由于某一周期内每支车流的日车流量波动情况难以用数学形式表达，而假设其服从某一概率分布或对应某一模糊隶属度函数都会存在不完全合乎实际的情况。因此以年为周期，对日均车流量的波动性进行描述。

1.2　模型参数与变量定义

V为铁路网络中所有车站的集合；i,j,s和t均为集合V的索引，即铁路网络中的车站；n为该铁路网络上的车站数量，V=\{i\left|i=\mathrm{1,2},\cdots ,n\right.\}；E为铁路网络中所有线路区段的集合；（s，t）为端点为s站和t站的线路区段，且线路方向为s站至t站，（s，t）\in E；G为铁路实际网络对应的支点网络，G=（V，E）；N_{ij}为i站到j站的车流量，即计划车流（编组计划执行期间的未来期望值）；l_{st}为弧段（s，t）的长度；L_{ij}为i站到j站的最短径路里程；C_{st}为弧段（s，t）的通过能力；m为列车平均编成辆数；\alpha为单线铁路的线路通过能力后备系数；\phi为车流走行径路的绕道率阈值：x_ij，st 为0-1变量，若i站至j站的车流径路中包括区段（s，t）则取值为1，否则取值为0。

1.3　模型构建

在铁路网络中车流量确定的情况下，考虑线路的通过能力及车流径路的里程限制等约束，以最小化车流在网络中的走行费用为优化目标Z_{\mathrm{1}}，根据车流不可拆分原则，建立以下混合整数线性规划模型，并将其命名为模型MIP_Ⅰ

Z_{\mathrm{1}}=\mathrm{m}\mathrm{i}\mathrm{n}\left(\sum _{\left(s,t\right)\in E}\sum _{i,j\in V}\sum _{i\ne j}{l}_{st}{N}_{ij}{x}_{ij,st}\right)（1）

s.t.

\sum _{\left(i,t\right)\in E}\sum _{i,j\in V}\sum _{i\ne j}{x}_{ij,it}^{}=\mathrm{1}（2）

\sum _{\left(s,j\right)\in E}\sum _{i,j\in V}\sum _{i\ne j}{x}_{ij,sj}^{}=\mathrm{1}（3）

\sum _{\left(i,j\right)\in E}\sum _{i,j\in V}\sum _{i\ne j}{x}_{ij,ij}^{}=\mathrm{1}（4）

\sum _{\left(s,t\right)\in E\backslash \left\{\left(i,j\right)\right\}}\sum _{i\in V}\sum _{j\in V|\left(i,j\right)\in E}{x}_{ij,st}^{}=\mathrm{0}（5）

\sum _{\left(s,t\right)\in E}\sum _{s\in V}\sum _{s\ne j}{x}_{ij,st}^{}=\sum _{\left(t,q\right)\in E}\sum _{q\in V}\sum _{q\ne i}{x}_{ij,tq}^{}\hspace{1em}\forall i\in V,

j\in V\backslash \left\{k|\left(i,k\right)\in E\right\},t\in V,t\ne i,

t\ne j,i\ne j（6）

\begin{array}{l}\sum _{\left(s,i\right)\in E}\sum _{s\in V}{x}_{ij,si}^{}+\sum _{\left(j,t\right)\in E}\sum _{t\in V}{x}_{ij,jt}^{}=\mathrm{0}\hspace{1em}\forall i\in V,\\ j\in V\backslash \left\{k|\left(i,k\right)\in E\right\}\end{array}（7）

\sum _{\left(s,t\right)\in E}{l}_{st}{x}_{ij,st}^{}\le \left(\mathrm{1}+\phi \right)L_{ij}\hspace{1em}\forall i\in V,j\in V,i\ne j（8）

\sum _{i,j\in V}\sum _{i\ne j}{N}_{ij}{x}_{ij,st}^{}\le \alpha m{C}_{st}\hspace{1em}\forall \left(s,t\right)\in E（9）

\sum _{\left(k,t\right)\in E}{N}_{ij}{x}_{ij,kt}^{}=\sum _{\left(t,k\right)\in E}{N}_{ij}{x}_{ij,tk}^{}\hspace{1em}\forall i\in V,j\in V,

t\in V,t\ne i,t\ne j,i\ne j（10）

x_{ij,st}^{}\in \left\{\mathrm{0,1}\right\}\hspace{1em}\forall i\in V,j\in V,\left(s,t\right)\in E,i\ne j（11）

模型MIP_Ⅰ中，式（1）为目标函数，表示所有车辆的走行成本之和。式（2）—式（11）为约束条件，式（2）表示对于从i站到j站的车流径路，其首条弧段应以i站为起点且仅能选择某一相邻车站作为该弧段的终点；类似地，式（3）表示对于从i站到j站的车流，其径路最末弧段应以j站为终点且最末弧段有且仅有一条；式（4）和式（5）表示对于相邻车站之间的车流，其仅能通过该相邻车站间的弧段进行运输，而不能选择其他弧段绕道运输；式（6）确保了车流径路的连续性，表示对于从i站到j站的车流，若其走行径路中包括t站，则组成该径路的弧段中，以t站为起点的弧段之和应等于以t站为终点的弧段之和；式（7）可辅助剔除路网中某一车流对应径路的不合理弧段，如对于从i站到j站的车流，其车流径路中显然不应该包括以i站为终点或以j站为起点的弧段；式（8）为车流径路里程约束，引入绕道率阈值\phi的概念，即对某一车流而言，可供其选择的车流径路包括最短径路，同时也包括物理里程在最短径路里程\phi倍之内的径路；式（9）为线路通过能力约束；式（10）为流量平衡约束；式（11）为变量取值范围约束。

2 需求不确定的铁路车流径路优化模型

2.1　问题描述

鲁棒优化以集合的形式对不确定因素进行描述，相比于随机优化和模糊规划，无须考虑不确定因素的分布规律或隶属度函数。采用鲁棒优化不确定集中最为常见的一种形式，即盒式不确定集（Box Uncertainty Set），其一般可表示为^［17］

B=\left\{\left(p_{\mathrm{1}},\cdots ,p_n\right)\left|\left|p_w-p_w^{\mathrm{*}}\right|\le {\sigma }_w,\forall w\right.\right\}（12）

式中：B为盒式不确定集合；p_{\mathrm{1}}和p_n分别为第1个和第n个不确定因素的取值；p_w和p_w^{\mathrm{*}}分别为不确定因素的取值和均值；w为不确定因素的索引；{\sigma }_w为不确定因素的最大偏离值。

相比于椭球不确定集（Ellipsoidal Uncertainty Set）或多面体不确定集（Polyhedral Uncertainty Set），盒式不确定集更易出现所有不确定参数均在区间集的上界或下界进行取值的情况，导致求解结果过于保守。为此，采用BERTSIMAS等^［14］提出的鲁棒对等理论并引入鲁棒代价（Price of Robustness）参数\Gamma，来调节模型的保守程度。以式（13）和式（14）这一经典的Soyster线性规划模型^［18］（简称为Soyster模型）为例介绍该鲁棒对等理论。值得注意的是，Soyster模型与一般的线性划模型形式相同，但其中的参数与变量并无实际含义；其与一般的线性规划模型的区别仅在于约束条件中含有不确定性的参数{\tilde{K}}_{ab}。

Z=\underset{b\in B}{\mathrm{m}\mathrm{a}\mathrm{x}}({{\rm P}}_b{\nu }_b)（13）

\sum _b{\tilde{K}}_{ab}{v}_b\le {{\rm H}}_a\hspace{1em}\forall a\in A（14）

式中：Z为该线性规划模型的目标函数值；{{\rm P}}_b和{{\rm H}}_a均为确定性参数，角标a和b分别表示其为属于凸集合A和B的变量；v_b为决策变量且其取值为正数。

设Soyster模型中参数{\tilde{K}}_{ab}\in \left[K_{ab}-{\widehat{K}}_{ab},K_{ab}+{\widehat{K}}_{ab}\right]，其中K_{ab}和{\widehat{K}}_{ab}分别为不确定性参数{\tilde{K}}_{ab}的均值和绝对波动量，则约束条件式（14）可转化为

\left\{\begin{array}{l}\sum _b{\tilde{K}}_{ab}{v}_b+\varphi \left(u,\Gamma \right)\le {{\rm H}}_a\hspace{1em}\forall a\in A\\ \sum _b{\tilde{K}}_{ab}{v}_b+\varphi \left(u,\Gamma \right)=\sum _b{K}_{ab}{v}_b+\\ \mathrm{m}\mathrm{a}\mathrm{x}\left[\sum _b{\widehat{K}}_{ab}{u}_b+\left({\Gamma }_a-⌊{\Gamma }_a⌋\right){\widehat{K}}_{a{e}_a}{u}_e\right]\end{array}\right.（15）

式中：\varphi \left(u,\Gamma \right) 为变量u与\Gamma组成的多项式，该多项式中含有不确定性参数的波动量；u_b为辅助变量，用于分离不确定性模型中的确定性参数和变量与不确定性参数和变量；{\Gamma }_a为第a个约束条件对应的不确定性水平，人为设定取值；{\widehat{K}}_{a{e}_a}为不确定性参数，其中角标e_a表示第a个约束条件对应的第e个不确定性参数，第a个约束条件共包括{\Gamma }_a个不确定性参数。

对于Soyster模型，当鲁棒代价参数{\Gamma }_a取值不为整数时，则不确定元素个数不超过 ⌊{\Gamma }_a⌋，同时可取{\widehat{a}}_{ij}的波动范围为 \left({\Gamma }_a-⌊{\Gamma }_a⌋\right){\widehat{K}}_{ab}；当{\Gamma }_a取值为整数时，则有\varphi \left(u,\Gamma \right)=\mathrm{m}\mathrm{a}\mathrm{x}\left(\sum _b{\widehat{K}}_{ab}{u}_b\right)。为简化计算，将鲁棒代价{\Gamma }_a取为整数。

为了使约束条件式（15）等价于式（14），需要添加约束条件

-u_b\le v_b\le u_b\hspace{1em}\forall b\in B（16）

若上述Soyster模型存在最优解v_b^{\mathrm{*}}，则必然有u_b=\left|v_b^{\mathrm{*}}\right|，即辅助变量u_b的取值应等于决策变量v_b的取值。

观察发现，当{\Gamma }_a取值为0时，鲁棒问题回归为确定问题；当{\Gamma }_a取值为\left|A\right|（即集合A中的元素个数）时，所有不确定元素均取上界值，这即是Soyster鲁棒对等理论^［18］。尽管这种情况下得到的鲁棒解在不确定因素发生任何变动时仍为可行解，但该解会使得极大化问题的目标函数值过小、极小化问题的过大，即求解结果过于保守。而采用鲁棒对等理论可以通过控制\Gamma的取值来调节鲁棒水平，有效避免出现最优解过于保守的情况。

2.2　辅助参数与变量定义

以模型MIP_Ⅰ为基础，采用上述鲁棒对等理论对需求不确定情形下的铁路车流径路优化问题进行描述，并对需求不确定条件下的铁路车流径路优化模型进行关于集合、参数、决策变量及对应角标的补充定义。{\Omega }_{ij}为车流量{\tilde{N}}_{ij}的不确定集合；{\tilde{N}}_{ij}为不确定环境下i站到j站车流量的实际取值；{\overline{N}}_{ij}和{\underline{N}}_{ij}分别为不确定性参数{\tilde{N}}_{ij}的上限和下限；{\theta }_{ij}为不确定性参数{\tilde{N}}_{ij}的波动系数；{\Gamma }_{st}为鲁棒代价参数，用以限制弧段 \left(s,t\right) 上发生车流量波动的车流数量；\lambda为单位转换系数；M为一个足够大的正数；z_{ij}为0-1变量，若i站到j站的车流量发生了波动则取值为1，否则取值为0；{\Upsilon }_{st}为连续型变量，为s站至t站的线路区段通过能力的补偿能力；{\delta }_{st}和{\eta }_{ij}_，st 为连续型变量，均为对偶变换时产生的变量；y_{ij}_，st 为0-1变量，用于分离鲁棒约束条件中的确定性参数与变量和不确定性参数与变量，其取值必须与决策变量x_{ij}^{}_，st 保持一致。

2.3　鲁棒模型构建

在考虑上述确定环境下铁路车流径路优化模型的基础上，对需求不确定的情况进行描述，模型中其他参数均是确定的。现将车流量的不确定集{\Omega }_{ij}定义为

{\tilde{N}}_{ij}\in {\Omega }_{ij}\equiv \left\{{\tilde{N}}_{ij}:{\underline{N}}_{ij}\le {\tilde{N}}_{ij}\le {\overline{N}}_{ij}\right\}（17）

其中，

\left\{\begin{array}{l}{\underline{N}}_{ij}=N_{ij}\left(\mathrm{1}-{\theta }_{ij}\right)\\ {\overline{N}}_{ij}=N_{ij}\left(\mathrm{1}+{\theta }_{ij}\right)\end{array}\right.\hspace{1em}\mathrm{0}\le {\theta }_{ij}\le \mathrm{1}

则线路通过能力约束式（9）可转化为

\left\{\begin{array}{l}\sum _{i,j\in V}\sum _{i\ne j}{N}_{ij}{x}_{ij,st}^{}+\varphi \left(y,\Gamma \right)\le \\ \alpha m{C}_{st}+m{\Upsilon }_{st}\hspace{1em}\forall \left(s,t\right)\in E\\ \sum _{i,j\in V}\sum _{i\ne j}{N}_{ij}{x}_{ij,st}^{}+\varphi \left(y,\Gamma \right)=\\ \sum _{i,j\in V}\sum _{i\ne j}{N}_{ij}{x}_{ij,st}^{}+\underset{\left(s,t\right)\in E}{\mathrm{m}\mathrm{a}\mathrm{x}}\left(\sum _{i,j\in V}\sum _{i\ne j}{\theta }_{ij}{N}_{ij}{y}_{ij,st}^{}\right)\end{array}\right.（18）

-y_{ij,st}^{}\le x_{ij,st}^{}\le y_{ij,st}^{}\hspace{1em}\forall i\in V,j\in V,\left(s,t\right)\in E,i\ne j（19）

y_{ij,st}^{}\in \left\{\mathrm{0,1}\right\}\hspace{1em}\forall i\in V,j\in V,\left(s,t\right)\in E,i\ne j（20）

{\Upsilon }_{st}\ge \mathrm{0}\hspace{1em}\hspace{1em}\hspace{1em}\forall \left(s,t\right)\in E（21）

考虑到模型MIP_Ⅰ中的通过能力不仅扣除了旅客列车占用能力，还扣除了一定的后备能力（单线设为20%），因此若(C_{st}+{\Upsilon }_{st})/C_{st}不超过1.2，即便最优解中补偿通过能力{\Upsilon }_{st}取值大于零，都可以认为未超过线路的实际货运通过能力。此外，如果令x_{ij}^{\mathrm{*}}_，st 为鲁棒车流径路问题的最优解，那么有y_{ij}^{}_，st=\left|x_{ij}^{*}\right._，\left.{}_{st}^{}\right|^［18］。

引入辅助决策变量z_{ij}后，式（18）左侧\varphi（y,{\Gamma }_{st}）可以用线性规划目标函数Z_{\mathrm{L}\mathrm{P}}表达为

Z_{\mathrm{L}\mathrm{P}}=\underset{\left(s,t\right)\in E}{\mathrm{m}\mathrm{a}\mathrm{x}}\left(\sum _{i,j\in V}\sum _{i\ne j}{\theta }_{ij}{N}_{ij}{y}_{ij,st}^{}{z}_{ij}\right)（22）

s.t.

\sum _{i,j\in V}\sum _{i\ne j}{z}_{ij}\le {\Gamma }_{st}（23）

\mathrm{0}\le z_{ij}\le \mathrm{1}\hspace{1em}\forall i\in V,j\in V,i\ne j（24）

每个线性规划问题都有对应的对偶问题。原问题（LP）可行有界，其对偶问题的目标函数与原问题一致，由此将LP模型转化为如下以目标函数Z_{\mathrm{D}\mathrm{P}}表达的对偶模型。

Z_{\mathrm{D}\mathrm{P}}=\underset{\left(s,t\right)\in E}{\mathrm{m}\mathrm{i}\mathrm{n}}\left(\sum _{i,j\in V}\sum _{i\ne j}{\eta }_{ij,st}+{\delta }_{st}{\Gamma }_{st}\right)（25）

s.t.

{\delta }_{st}+{\eta }_{ij,st}\ge {\theta }_{ij}{N}_{ij}{y}_{ij,st}\hspace{1em}\forall i\in V,j\in V,i\ne j（26）

{\eta }_{ij,st}\in \left\{\mathrm{0,1}\right\}\hspace{1em}\forall i\in V,j\in V,i\ne j（27）

{\delta }_{st}\ge \mathrm{0}（28）

将对偶问题DP带入至约束条件式（18），可得到如下基于鲁棒对等理论的铁路车流径路优化模型MIP_Ⅱ，其目标函数Z_{\mathrm{2}}为总车流走行费用与总补偿通过能力相关成本之和。

\begin{array}{l}{Z}_{\mathrm{2}}=\mathrm{m}\mathrm{i}\mathrm{n}\left\{\sum _{\left(s,t\right)\in E}\sum _{i,j\in V}\sum _{i\ne j}\left[l_{st}\left(N_{ij}{x}_{ij,st}+{\eta }_{ij,st}\right)\right]+\right.\\ \left.\sum _{\left(s,t\right)\in E}{l}_{st}{\delta }_{st}{\Gamma }_{st}+\lambda \sum _{\left(s,t\right)\in E}{\Upsilon }_{st}\right\}\end{array}（29）

s.t.

式（2）—式（8）、式（11）、式（19）—式（21）、式（26）—式（28）

\sum _{i,j\in V}\sum _{i\ne j}\left(N_{ij}{x}_{ij,st}+{\eta }_{ij,st}\right)+{\delta }_{st}{\Gamma }_{st}\le \alpha m{C}_{st}+m{\Upsilon }_{st}

\forall \left(s,t\right)\in E（30）

{\eta }_{ij,st}\le M{y}_{ij,st}\hspace{1em}\forall i\in V,j\in V,\left(s,t\right)\in E,i\ne j（31）

\sum _{\left(k,t\right)\in E}\sum _{k\in V}{\eta }_{ij,kt}^{}{y}_{ij,kt}^{}=\sum _{\left(t,k\right)\in E}\sum _{k\in V}{\eta }_{ij,tk}^{}{y}_{ij,tk}^{}\hspace{1em}\forall i\in V,

j\in V,t\in V,t\ne i,t\ne j,i\ne j（32）

模型MIP_Ⅱ中，式（29）右侧的最后一项为总补偿通过能力相关成本，其中参数\lambda具有2个作用，一是将由于车流量波动产生的补偿通过能力换算为吨公里，二是调整总补偿通过能力这一项成本的数量级；式（30）为考虑车流量波动的带有补偿的线路通过能力约束条件；约束条件式（31）用来辅助判断决策变量的取值；式（32）为流量平衡约束条件。

模型MIP_Ⅱ的目标函数式（29）考虑了2类成本，第1类是车流的总走行成本Z_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{v}\mathrm{e}\mathrm{l}}，通过计算每条弧段上车流量与该条弧段里程的乘积并求和得到；第2类是总的线路通过能力相关成本Z_{\Upsilon }。为便于说明其含义，将式（29）拆分为

\begin{array}{l}{Z}_{\mathrm{2}}=Z_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{v}\mathrm{e}\mathrm{l}}+Z_{\Upsilon }=\\ \sum _{\left(s,t\right)\in E}\sum _{i,j\in V}\sum _{i\ne j}{l}_{st}\left(N_{ij}{x}_{ij,st}+{\eta }_{ij,st}\right)+\\ \sum _{\left(s,t\right)\in E}{l}_{st}{\delta }_{st}{\Gamma }_{st}+\lambda \sum _{\left(s,t\right)\in E}{\Upsilon }_{st}\end{array}（33）

其中，

Z_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{v}\mathrm{e}\mathrm{l}}=\sum _{\left(s,t\right)\in E}\sum _{i,j\in V}\sum _{i\ne j}{l}_{st}\left(N_{ij}{x}_{ij,st}+{\eta }_{ij,st}+{\delta }_{st}{\Gamma }_{st}\right)（34）

之所以在模型MIP_Ⅱ的目标函数中考虑线路因车流量波动而增加的补偿通过能力，是因为最小化总补偿通过能力在理论层面上有较为重要的意义。尽管这部分增加的补偿能力可以通过线路预留的后备能力提供，而不必对铁路线路进行改造，但考虑到运营安全等因素，应尽可能控制后备能力的使用程度。

若采用d_{st}^{}表示车流量波动时弧段（s，t）上总车流量大小，则可以发现，d_{st}^{}与式（34）中的N_{ij}{x}_{ij}_，st+{\eta }_{ij}_，st+{\delta }_{st}{\Gamma }_{st}有如下关系

d_{st}^{}=\sum _{i,j\in V}\sum _{i\ne j}\left(N_{ij}{x}_{ij,st}+{\eta }_{ij,st}\right)+{\delta }_{st}{\Gamma }_{st}（35）

式（35）中，N_{ij}{x}_{ij}_，st 为车流量发生波动时，分配到弧段（s，t）上不确定参数{\tilde{N}}_{ij}中确定的部分，即无论车流量{\tilde{N}}_{ij}是否发生波动，计划车流量N_{ij}都需要分配到铁路运输网络的弧段上。而车流{\tilde{N}}_{ij}是否发生波动，波动后弧段上的车流总量又如何计算，这就需要进一步分析变量{\delta }_{st}和{\eta }_{ij}^{}_，st。

变量{\delta }_{st}和{\eta }_{ij}^{}_，st 来自文中的对偶模型（DP），根据强对偶理论，若对偶问题与原问题都有可行解，则它们都有最优解，且目标函数的最优值必然相等。线性规划模型LP的目标函数表示弧段（s，t）上的最大车流量波动值。相应地，其对偶模型DP的目标函数也具有同样的含义。

为更好地利用鲁棒代价\Gamma说明铁路货物运输系统的鲁棒性，引入运输需求满足度S_{st}的概念，定义为不违反线路区段（s，t）通过能力约束的概率，其数学表达形式为

S_{st}\left(x,\Gamma \right)=\mathrm{P}\mathrm{r}\left(\sum _{i,j\in V}\sum _{i\ne j}{\tilde{N}}_{ij}{x}_{ij,st}^{\mathrm{*}}\le \alpha m{C}_{st}+m{\Upsilon }_{st}\right)（36）

式中：Pr（·）为概率函数；x_{ij}^{\mathrm{*}}_，st为模型MIP_Ⅱ的最优解。

则最优解满足约束条件（18）的概率S_{st}可由下式计算^［19］

\begin{array}{l}{S}_{st}\left(x,\Gamma \right)=\mathrm{P}\mathrm{r}\left(\sum _{i,j\in V}\sum _{i\ne j}{\tilde{N}}_{ij}{x}_{ij,st}^{\mathrm{*}}\le \alpha m{C}_{st}+m{\Upsilon }_{st}\right)\ge \\ \mathrm{1}-\left(\mathrm{1}-\varpi \right)\psi \left(\left|E\right|,⌊g⌋\right)-{\displaystyle \sum _{l=⌊v⌋+\mathrm{1}}^n}\psi \left(\left|E\right|,h\right)\end{array}（37）

\psi \left(\left|E\right|,h\right)=\left\{\begin{array}{l}\frac{\mathrm{1}}{{\mathrm{2}}^{\left|E\right|}}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}h=\mathrm{0}\mathrm{或}\mathrm{h}=\left|E\right|\\ \frac{\mathrm{1}}{\sqrt[]{\mathrm{2}\mathrm{\pi }}}\left(\sqrt[]{\frac{\left|E\right|}{(\left|E\right|-h)h}}\right)\times \\ \mathrm{e}\mathrm{x}\mathrm{p}\hspace{1em}\left(\left|E\right|\mathrm{l}\mathrm{n}\left(\frac{\left|E\right|}{\mathrm{2}\left|E\right|-\mathrm{2}h}\right)+\right.\\ \hspace{1em}\hspace{1em}\hspace{1em}\left.h\mathrm{l}\mathrm{n}\left(\frac{\left|E\right|-h}{h}\right)\right)\hspace{1em}\mathrm{其他}\end{array}\right.（38）

其中，

g=\left({\Gamma }_{st}+\left|E\right|\right)/\mathrm{2}

\varpi =g-⌊g⌋

式中：\left|E\right| 为铁路网络上区段的总数量；h为可能的鲁棒代价参数，其取值为正整数，且h\in \left[\mathrm{0},\left|E\right|\right]；g和\varpi均为参数，可通过鲁棒代价参数{\Gamma }_{st}和铁路网络区段数量\left|E\right|求得。

3 算例验证

3.1　车流量确定与波动时的车流径路优化对比

为对比需求确定即车流量确定（对应模型MIP_Ⅰ）及需求不确定即车流量波动（对应模型MIP_Ⅱ）下铁路车流径路优化结果，以2019 INFORMS Rail Application Section（RAS）竞赛^［20］的公开数据集为例进行计算及分析。该虚拟铁路网中共包括8个支点及56支OD需求量，对应拓扑图如图1所示。图中：Y1—Y8为支点站。分析公开数据集可发现，该路网中线路上下行的通过能力均相同，但实际上受到线路平纵断面的影响，即便是单线铁路，其上下行通过能力一般也会有所差别；本文提出的优化模型无论上下行能力是否相同，都同样能够适用。

设置参数取值：车流量不确定集中的\theta为10%，列车平均编成辆数m为50，绕道率阈值\phi为1.25，线路通过能力后备系数为20%，鲁棒代价为2，单位换算系数\lambda为10⁴。

考虑到车流量确定和波动时，2种车流径路优化模型均为线性模型，采用大规模优化器Gurobi 9.0.2进行求解，硬件运行环境为Intel（R）Core（TM）i5-7267U CPU（3.10 GHz）。基于该配置条件，车流量确定时的优化方案走行成本为2 119 046 tkm，车流量波动（{\Gamma }_{st}=\mathrm{2}）时的优化方案走行成本为2 220 692 tkm，由于车流量波动而产生的补偿通过能力为18.9车；车流量的波动共造成了11支车流的径路发生了变化，对应的2种车流径路方案见表1和表2。表中：（i∶j）表示该行对应的车流径路起点为i站且终点为j站；{\zeta }_{ij}^{\mathrm{*}}和{\zeta }_{ij}^{}分别为车流量确定和波动（{\Gamma }_{st}=\mathrm{2}）时的最优车流径路里程；{\varphi }_{ij}^{\mathrm{*}}和{\varphi }_{ij}^{}为表示车流量确定和波动时最优车流径路的绕道率。相邻支点站间的车流必然通过两站间的弧段运输，因此表1和表2中均未列出相邻支点站间的车流径路方案；此外，2种径路方案均仅列出了车流量波动时径路发生变化的方案。

表3给出了发生波动后的车流量和初始车流量。将表3与表1和表2中的最优车流径路方案进行对比，可以分析出车流量波动后径路发生变化的原因。根据表1和表2给出的车流径路方案，以车流N_{\mathrm{16}}为例，该车流在需求确定及波动2种情形下的具体走行径路及对该径路产生影响的其他车流径路分别如图2和图3所示。图中：蓝色线条表示路网结构；红色、黑色箭头线分别表示车流N_{\mathrm{16}}的走行径路和其他车流的走行径路。综合表3以及图2和图3可得到如下结论。

（1）在车流量确定的情形下，Y1—Y6的车流量为109车，通过径路Y1→Y2→Y3→Y4→Y5→Y6进行运输；其他占用该线路区段的车流包括N_{\mathrm{13}}=\mathrm{127}，N_{\mathrm{14}}=\mathrm{114}，N_{\mathrm{15}}=\mathrm{120}，N_{\mathrm{18}}=\mathrm{98}，N_{\mathrm{23}}=\mathrm{102}，N_{\mathrm{24}}=\mathrm{124}，N_{\mathrm{28}}=\mathrm{47}，N_{\mathrm{73}}=\mathrm{53}，N_{\mathrm{78}}=\mathrm{73}。

（2）在车流量波动的情形下，Y1—Y6的车流量虽仍为109车，但通过径路Y1→Y2→Y7→Y4→Y5→Y6运输；其他占用该线路区段的车流包括N_{\mathrm{13}}=\mathrm{139.7}，N_{\mathrm{14}}=\mathrm{114}，N_{\mathrm{15}}=\mathrm{120}，N_{\mathrm{18}}=\mathrm{98}，N_{\mathrm{23}}=\mathrm{102}，N_{\mathrm{24}}=\mathrm{136.4}，N_{\mathrm{25}}=\mathrm{85}，N_{\mathrm{26}}=\mathrm{85}，N_{\mathrm{28}}=\mathrm{47}，N_{\mathrm{73}}=\mathrm{53}，N_{\mathrm{78}}=\mathrm{73}。

（3）车流量确定时，通过弧段（\mathrm{Y}\mathrm{2},\mathrm{Y}\mathrm{3}）的车流量为967车；而当车流量波动（{\Gamma }_{st}=\mathrm{2}）时，不仅通过（\mathrm{Y}\mathrm{2},\mathrm{Y}\mathrm{3}）的车流量变为1 000.1车，而且因途经弧段（\mathrm{Y}\mathrm{2},\mathrm{Y}\mathrm{3}）的车流量N_{\mathrm{13}}和N_{\mathrm{24}}分别取其不确定集的上界值，导致车流量确定时不经过弧段（\mathrm{Y}\mathrm{2},\mathrm{Y}\mathrm{3}）的车流量N_{\mathrm{24}}和N_{\mathrm{25}}也会改为经由该弧段通过，而原本途经（\mathrm{Y}\mathrm{2},\mathrm{Y}\mathrm{3}）的车流量N_{\mathrm{16}}则会改为经由弧段（\mathrm{Y}\mathrm{2},\mathrm{Y}\mathrm{7}）通过。这是由于在以车流走行费用最小化作为优化目标的背景下，该弧段会发生如下车流变化：部分车流量确定时经过此处的车流会在车流量发生波动后改道其他弧段，而原本未从此处通过的车流则会改道于此处。

车流量确定和波动时的线路通过能力占用率、车流量波动时对通过能力的补偿分别见表4和表5。表中：d_{st}^{\mathrm{*}}为车流量确定时弧段（s，t）承担的总车流量大小，按式（39）计算；{\lambda }_{st}^{\mathrm{*}}和{\lambda }_{st}分别为车流量确定和波动时的弧段（s，t）通过能力占用率，按式（40）和式（41）计算。

d_{st}^{\mathrm{*}}=\sum _{i,j\in V}\sum _{i\ne j}{N}_{ij}{x}_{ij}^{st}\hspace{1em}\hspace{1em}\forall \left(s,t\right)\in E（39）

{\lambda }_{st}^{\mathrm{*}}=d_{st}^{\mathrm{*}}/\alpha C_{st}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\forall \left(s,t\right)\in E（40）

{\lambda }_{st}^{}=d_{st}^{}/\left(\alpha C_{st}+{\Upsilon }_{st}\right)\hspace{1em}\forall \left(s,t\right)\in E（41）

综合表4和表5可以发现：车流量确定情况下最优解的线路通过能力占用率最高已达到99.80%，对应弧段（\mathrm{Y}\mathrm{8},\mathrm{Y}\mathrm{3}），因此在车流量波动（{\Gamma }_{st}=\mathrm{2}）的情形下，部分线路区段如（\mathrm{Y}\mathrm{4},\mathrm{Y}\mathrm{5}）和（\mathrm{Y}\mathrm{7},\mathrm{Y}\mathrm{4}），必然要在现有通过能力的基础上补偿更多的通过能力，才能满足货物运输的需求；由于模型中参与计算的线路通过能力均指去除旅客列车占用能力后的剩余能力，并在此基础上扣除了一部分能力作为后备能力，保留这部分后备能力的作用之一是为了适应日常车流量的波动，因此，若弧段上总车流量的波动不超过该弧段现有通过能力的20%，那么可以认为该弧段的后备通过能力足以满足发生波动的这部分车流量的运输需求，如弧段（\mathrm{Y}\mathrm{4},\mathrm{Y}\mathrm{5}）在车流量波动（{\Gamma }_{st}=\mathrm{2}）时需要补偿13.6车的通过能力，而该弧段最大通过能力为750车，后备能力为150车，需要补偿的通过能力在其后备能力的范围之内，因此无须对铁路线路进行改造，即可满足车流量波动情形下的货物运输需求。

3.2　不同鲁棒水平下的车流径路优化结果

为了挑选经济合理的鲁棒车流径路方案，需要衡量系统的成本与运输需求满足度。基于上述案例，假设车流量偏移值仍为日均车流的10%，鲁棒代价{\Gamma }_{st}的取值从0增加到9，得到不同鲁棒代价下的路网系统成本与运输需求满足度变化曲线如图4所示。从图4可以看出：随着鲁棒代价取值的增加，最优解的走行成本从2 119 046 tkm增加到2 391 459 tkm，系统鲁棒水平最高（{\Gamma }_{st}=\mathrm{9}）时总成本相比于需求确定（{\Gamma }_{st}=\mathrm{0}）时的最优解增加了102.14%，即当铁路网中所有车流量均取其上界时，系统总成本会增加一倍。

该铁路货运系统中，运输需求满足度的变化及不同鲁棒代价下需要的线路补偿通过能力见表6。从表8可以发现：当铁路网络中车流量确定时，货运系统的运输需求满足度仅为35.6%；当鲁棒代价{\Gamma }_{st}=\mathrm{5}时，系统的运输需求满足度达到90.5%，此时走行成本仅增加9.7%，总成本仅增加68.6%，即1 453 600 tkm，相当于每增加1%的运输需求满足度需要对应增加26 477.23 tkm的系统总成本。结合图4可以发现，{\Gamma }_{st}\le \mathrm{5}时运输需求满足度的增幅处于较高水平，但当{\Gamma }_{st}>5后增幅逐渐放缓，说明系统成本的增加无法换来快速增长的运输需求满足度，也就是系统总成本增速将远远大于运输需求满足度的增速，例如当{\Gamma }_{st}=\mathrm{9}时，系统的运输需求满足度达到99.8%，但此时走行成本增加12.9%，总成本增加102.1%，即2 164 413 tkm，意味着每增加1%的运输需求满足度须对应增加33 713.6 tkm的系统总成本。综上分析，可认为鲁棒代价{\Gamma }_{st}=\mathrm{5}时的最优鲁棒车流径路方案最为经济合理，且能够作为可行方案应用于超过90%以上的运输需求发生波动的不确定情形。

不同鲁棒水平下铁路货运系统的走行成本不同，其最优运输方案也不同。随着系统鲁棒性的增强，走行成本及补偿通过能力也有所增加，对应的运输需求满足度也随之增加。管理者可以通过对比不同鲁棒代价下的走行成本及补偿通过能力，选择合适的运输需求满足度；并通过控制系统的鲁棒性，得到能够满足大部分或全部车流量波动的最佳车流径路方案。

3.3　实际算例

以中国铁路沈阳局集团有限公司（简称“沈阳局”）管辖范围内的部分主要车站组成的货物运输网络为例，对上文中提出的模型进行应用与验证。沈阳局的管辖范围以沈阳为中心，辐射辽宁省、吉林省、黑龙江省南部、河北省东北部的部分地区等，电气化铁路总里程达7 843.9 km，实际路网结构如图5所示，简化路网由25个车站组成，如图6所示。依据简化路网随机生成552组OD数据，其中非零车流量496组。假设车流量的波动可用盒式不确定集表示，波动上界取10%。基于以上数据，分别对车流量确定和车流量波动时的车流径路优化模型进行求解。

求解发现当{\Gamma }_{st}=\mathrm{6}时，不确定模型的求解时间已达到6 521 s。为合理地减少求解时间，本例设定优化间隙为2.00%，即当目标下界和上界之间的差距小于2.00%乘以现有目标函数值时，求解器将终止并获得当前最优方案。此外，图6的简化路网中共包括39个线路区间，相应的鲁棒代价值应从0增加到39；但当鲁棒代价增加到23时，系统的运输需求满足度已达到了99.99%，此时走行成本增加了12.64%，总成本增加了52.74%，因此，仅列出鲁棒代价由0增加到23时的最优成本，见表7。

从表9可以看出：车流量确定时的总走行成本最小，为3 997 925 tkm；车流量波动且当鲁棒代价{\Gamma }_{st}=\mathrm{11}时，系统的运输需求满足度为94.55%，走行成本增加12.10%，总成本增加43.11%；当{\Gamma }_{st}值继续增加时，运输需求满足度上升程度明显放缓，这表明即使增加线路通过能力，也难再大幅提升运输需求满足度。图7和图8分别给出车流量确定和波动（{\Gamma }_{st}=\mathrm{11}）时的部分车流径路方案示意图。结合图6，由图7和图8可以发现：车流量发生波动后，车流的走行径路更多地集中在通辽—义县、义县—盘锦、盘锦—锦州区段，这是由于车流量确定条件下，途经沈阳西站的车流密度相较于路网上其他节点更密集；当途经沈阳西站的车流量发生波动时，需要先满足流量较大的车流运输需求。致使原本经过沈阳西站的车流量相对较小的车流需要转移到其他区段上。

我国铁路货物列车编组计划的编制为逐级负责制，在国铁集团的集中领导下编制全路跨局列车编组计划；各路局负责编制跨局区段及管内列车编组计划。本案例在规模上属于中等水平，说明本模型可以为各铁路局制定管内列车编组计划时提供必要的适应日常车流量波动的车流径路方案。

4 结 语

车流径路是编制货物列车编组计划最主要的依据之一。通常情况下，我国铁路货物列车编组计划是在车流径路已知的情况下进行编制的。车流径路方案的制定基于固定的日均车流量，但在现场执行中可能会出现实际车流量超过计划日均车流量的情况，因此，在确定径路方案时考虑车流量的波动性，能够为后续制定编组计划提供更符合实际应用的编制依据。本文选择鲁棒优化的理论和方法来描述车流量的波动性，考虑线路通过能力约束及网络中流量平衡约束，设定绕道率阈值参数来控制车流径路里程的合理性，建立了需求确定的铁路车流径路优化模型。在此基础上，建立了基于盒式不确定集的需求不确定的径路优化模型，通过计算能力约束条件满足的概率，并引入了运输需求满足度指标，用以反映最优径路方案的稳定程度。通过小型算例和实际算例对模型进行验证，证实了本模型的正确性与实用性，表明将鲁棒优化理论与方法应用到需求不确定下的车流径路与编组计划的一体化编制问题中具有重要的理论应用价值。

在小型算例中，对比车流量确定和波动时的径路方案发现，Y1—Y6的车流量均为109车，车流量确定时通过弧段（\mathrm{Y}\mathrm{2},\mathrm{Y}\mathrm{3}）的车流量为967车，而当车流量波动（{\Gamma }_{st}=\mathrm{2}）时则变为1 000.1车，且当部分车流量发生波动时，对应的最优车流径路方案随之发生变化；分析不同的鲁棒代价对目标函数值和运输需求满足度的影响发现，鲁棒代价可以控制最优方案的成本和违反能力约束的概率之间的平衡关系，具体到本算例，{\Gamma }_{st}=\mathrm{5}时处于系统成本与运输需求满足的平衡点，此时的车流径路方案最为合理，能够作为可行方案应用于超过90%以上的运输需求发生波动的不确定情形。

在沈阳局管内部分主要车站组成的货物运输网络中，分别求解得到车流量确定和波动时的2种车流径路方案。求解结果表明{\Gamma }_{st}=\mathrm{11}时的车流径路方案最优，此时系统的运输需求满足度为94.55%，走行成本增加12.10%，总成本增加43.11%，这意味着在总成本较低的情况下系统运输需求满足度达到较高的水平，而此后难以通过增加成本来大幅提升运输需求满足度。

通过扩大案例的规模，说明了基于需求不确定的铁路车流径路鲁棒优化模型的应用性，管理者可以通过权衡运输需求满足度和成本来选择合适的车流径路方案。但案例规模的扩大也带来了求解时间的快速增长，在未来的研究中，可以对求解该类模型的启发式算法进行研究设计，以求能够在短时间内得到合理的全铁路网车流径路方案。

参考文献

[1] DANTZIG G， FULKERSON R， JOHNSON S. Solution of a Large-Scale Traveling-Salesman Problem ［J］. Journal of the Operations Research Society of America， 1954， 2 （4）： 393-410.

[2] EL-SHAMIR EZUGWU A， ADEWUMI A O， FRÎNCU M E. Simulated Annealing Based Symbiotic Organisms Search Optimization Algorithm for Traveling Salesman Problem ［J］. Expert Systems with Applications， 2017， 77 （7）： 189-210.

[3] 陈科胜，鲜思东，郭鹏.求解旅行商问题的自适应升温模拟退火算法［J］.控制理论与应用，2021，38（2）：245-254.

[4] CHEN Kesheng， XIAN Sidong， GUO Peng. Adaptive Temperature Rising Simulated Annealing Algorithm for Traveling Salesman Problem ［J］. Control Theory & Applications， 2021， 38 （2）： 245-254. in Chinese

[5] TINIÇ G O， KARASAN O E， KARA B Y， et al. Exact Solution Approaches for the Minimum Total Cost Traveling Salesman Problem with Multiple Drones ［J］. Transportation Research Part B： Methodological， 2023， 168： 81-123.

[6] MALANDRAKI C， DASKIN M S. Time Dependent Vehicle Routing Problems： Formulations， Properties and Heuristic Algorithms ［J］. Transportation Science， 1992， 26 （3）： 185-200.

[7] MOGHDANI R， SALIMIFARD K， DEMIR E， et al. The Green Vehicle Routing Problem： a Systematic Literature Review ［J］. Journal of Cleaner Production， 2021， 279： 123691.

[8] KUCUKOGLU I， DEWIL R， CATTRYSSE D. The Electric Vehicle Routing Problem and Its Variations： a Literature Review ［J］. Computers & Industrial Engineering， 2021， 161： 107650.

[9] KONSTANTAKOPOULOS G D， GAYIALIS S P， KECHAGIAS E P. Vehicle Routing Problem and Related Algorithms for Logistics Distribution： a Literature Review and Classification ［J］. Operational Research， 2022， 22 （3）： 2033-2062.

[10] 刘畅.货物运输“公转铁”与铁路车流径路联合优化模型［J］.铁道学报，2020，42（12）：18-27.

[11] LIU Chang. Joint Optimization Model for Diverting Freight Flow from Road to Railway and Routing Railway Car Flow ［J］. Journal of the China Railway Society， 2020， 42 （12）： 18-27. in Chinese

[12] 何世伟，刘明玮，冯骁 等.考虑路网点线能力协调的铁路车流径路优化模型［J］.北京交通大学学报，2021，45（1）：1-7.

[13] HE Shiwei， LIU Mingwei， FENG Xiao， et al. Optimizing Model for Car Routing Considering Point-Line Capacity Coordination in Railway Network ［J］. Journal of Beijing Jiaotong University， 2021， 45 （1）： 1-7. in Chinese

[14] 薛锋，刘泳博，户佐安 等.基于动态罚函数的铁路车流分配与径路优化模型［J］.西南交通大学学报，2022，57（5）：941-948，959.

[15] XUE Feng， LIU Yongbo， HU Zuoan， et al. Railcar Traffic Distribution and Route Optimization Model Based on Dynamic Penalty Function ［J］. Journal of Southwest Jiaotong University， 2022， 57 （5）： 941-948， 959. in Chinese

[16] 李宵寅，魏玉晓.基于不确定参数的列车编组计划优化模型研究［J］.铁道运输与经济，2015，37（11）：26-31.

[17] LI Xiaoyin， WEI Yuxiao. Study on Optimization Model of Train Formation Plan Based on Uncertain Parameter ［J］. Railway Transport and Economy， 2015， 37 （11）： 26-31. in Chinese

[18] 陈崇双，王慈光，杨运贵，等.不确定条件下分组列车在编成站的集结参数研究［J］.铁道学报，2011， 33（5）：1-7.

[19] CHEN Chongshuang， WANG Ciguang， YANG Yungui， et al. Research on the Accumulation Parameter of Multi-Section Train in Formation Station under the Conditions of Uncertainty ［J］. Journal of the China Railway Society， 2011， 33 （5）： 1-7. in Chinese

[20] 严余松，户佐安，李宵寅.基于车流量波动的列车编组计划与车流径路综合优化［J］.交通运输系统工程与信息，2017，17（4）：124-131.

[21] YAN Yusong， HU Zuoan， LI Xiaoyin. Comprehensive Optimization of Train Formation Plan and Wagon-Flow Path Based on Fluctuating Wagon-Flow ［J］. Journal of Transportation Systems Engineering and Information Technology， 2017， 17 （4）： 124-131. in Chinese

[22] GABREL V， LACROIX M， MURAT C， et al. Robust Location Transportation Problems under Uncertain Demands ［J］. Discrete Applied Mathematics， 2014， 164： 100-111.

[23] LEFEVER W， TOUZOUT F A， HADJ-HAMOU K， et al. Benders’ Decomposition for Robust Travel Time-Constrained Inventory Routing Problem ［J］. International Journal of Production Research， 2021， 59 （2）： 1-25.

[24] BEN-TAL A， NEMIROVSKI A. Robust Solutions of Uncertain Linear Programs ［J］. Operations Research Letters， 1999， 25 （1）： 1-13.

[25] SOYSTER A L. Technical Note-Convex Programming with Set-Inclusive Constraints and Applications to Inexact Linear Programming ［J］. Operations Research， 1973， 21 （5）： 1154-1157.

[26] BERTSIMAS D， SIM M. The Price of Robustness ［J］. Operations Research， 2004， 52 （1）： 35-53.

[27] INFORMS. RAS Problem Solving Competition ［EB/OL］. ［2023-10-13］.

基金资助

国家重点研发计划先进轨道交通专项(2018YFB1201402)