基于拉格朗日松弛算法的列车运行图优化

Train Working Diagram Optimization Based on Lagrange Relaxation Algorithm

摘要: 本文主要研究了列车运行图优化问题。在我们所研究的问题当中，根据单线轨道区间的特点，为了使列车保证经济效益同时也要合理地在规定时刻内运行，这就需要我们寻找一个优化的列车运行图。由于在铁路线上运行的列车种类很多导致速度要求也不尽相同，我们主要针对各种速度不同的列车运行图(即非平行运行图)进行展开研究。本文的主要工作是：一是为进一步深入研究该问题，提出了一类基于时–空图构建的列车非平行运行图优化模型；二是针对该模型我们设计相应的启发式拉格朗日松弛算法进行求解；三是通过数值实验对比了标准次梯度法和修正次梯度法的收敛情况。

Abstract: This paper mainly studies the optimization problem of train working diagram. In the problem we studied, according to the characteristics of the single track segment, in order to ensure the economic benefits of the train while running reasonably within the specified time, we need to find an optimized train working diagram. Because there are many types of trains running on the railway line, which lead to different speed requirements, we mainly study the train working diagrams with different speeds (i.e., non-parallel operation diagrams). The main contributions of this paper are as follows: firstly, in order to further study this problem, a kind of optimization model of non-parallel train operation graph based on time-space graph is proposed; secondly, we design the corresponding heuristic Lagrange relaxation algorithm to solve the model; thirdly, the convergence of the standard sub-gradient method and the modified sub-gradient method is compared by numerical experiments.

1. 引言

近年来，使用运筹学方法对铁路运输组织过程中的各种问题进行计算与优化，得到了广泛的关注与应用 [1] 。本文主要研究铁路单线区间上的列车运行图优化问题(TWDP)，即在一段单线区间上，任意相邻车站之间均为单轨铁路，在规划时间段内有若干列车等待通过。如图1所示，一趟列车通过区间所对应的收益是一个关于时刻的函数。使得函数值最大的时刻就是这趟列车通过区间的最优时间。一般情况下列车通过区间的实际时间与最优时间相差越大收益就越小。目标是通过规划列车运行图使得所有列车通过本区间的总收益达到最大。由于在铁路线上运行的列车很多，无论是在客运专线、货运专线或客货列车混合运行有线，各种客货列车的速度和要求也不尽相同。因此，按照列车速度的特点，可将列车运行图分为平行运行图或非平行运行图。采用平行运行图可以达到最大通过能力，但只能在铁路区间通过能力比较紧张的情况下使用。在通常情况下，采用的是非平行运行图。所以，本文主要按照列车运行速度的不同研究非平行运行图。如图2所示的非平行列车运行图：纵坐标间的1、2表示车站，横坐标表示时间，纵坐标表示里程，每条斜线代表一趟列车。1、2号列车由1站出发驶向2站，3、4号列车由2站出发驶向1站。

TWDP的一般数学模型可以归结为一类带软时间窗约束的车辆路径问题或者基于时–空网络的优化问题 [2] ，而我们所考虑的问题更偏向于后者。有些模型虽然包括了单线轨道区间的例子，但在优化调度计划以确定应该解决列车技术动作时，由于无法量化的限制，该方案可能无法实现 [3] [4] 。在客运方面，一些研究将旅客与列车运行图编制系统有机地结合并应用到实际中 [5] [6] 。此外，廖正文等 [7] 针对累积流变量模型的特点提出基于双线运行图的0~1整数规划模型。最近几年，由于中国高速铁路迅速发展，因而高质量的运行图优化应以最好的效益满足运输需求。高如虎等 [8] 研究了增开列车条件下高速铁路列车运行图调整，并综合考虑列车车站进路的影响。孟学雷等 [9] 以时间要素为研究对象，采用概率统计理等理论与技术手段，计算分析时间偏差与偏离差等指标，提出了基于列车运行图稳定性与区间缓冲时分和停站缓冲时分的相互关系的运行图稳定性研究。列车运行图问题在后续的混合多轨、单轨铁路网和适应动态需求环境下铁路快速运输线路的铁路节点中被提出 [10] [11] [12] 。在列车运行图的优化方法基础上，王凯等 [13] 针对国内外计算机编制列车运行图存在的不足方面有进一步的研究方向。对于铁路列车运行图绘制的需求和现状，乐建炜等 [14] 阐述了实现列车运行图在线可视化的技术要点、实现步骤和生成结果，为后续参与运输调度领域信息系统规划、设计和研发提供技术支持与参考，提高我国铁路运行图在线绘制的技术水。综上所述，目前关于铁路区间上优化列车运行图的研究主要集中于如何压缩列车通过时间，

Figure 1. A profit function diagram of the train

图1. 某列车的收益函数示意图

Figure 2. Non-parallel train working diagram

图2. 非平行列车运行图

提升区间通过能力。在我们所研究的问题当中通过设置收益函数可以转化为上述问题，因此本文所研究的问题更加一般化。

针对列车运行图问题规模太大或者结构太复杂的情况，我们采取分解方法的策略来解决整个数学模型。主要的分解方法包括分支–定界法和拉格朗日松弛算法。一些列车运行图问题通过采用分支–定界法，从而有效获得具有保证最优性的可行调度 [14] [15] 。但对于大规模离散优化模型来说，拉格朗日松弛算法并不像分支–定界法一样采用直接搜索最优解的策略，而是希望通过分解模型来计算其最优目标函数值的上界和下界 [12] 。众所周知，由于组合优化问题的复杂性和模型的特点，20世纪70年代引入的拉格朗日松弛技术得到实现。自此，拉格朗日松弛算法的应用已经发展到为组合优化问题提供了最好的现有算法，可能更好地解决相关的实际规模问题。在Song等 [16] 的研究中，他们不仅采用了标准的拉格朗日松弛方法，而且还使用了增广的拉格朗日松弛技术来处理车辆路径的对称性问题。Hsnyb等 [17] 提出了一种可以分解为每个路径的若干子问题拉格朗日松弛方法，其中车辆循环约束被松弛。Caprara等 [18] 推导出一个以拉格朗日方式松弛的整数线性规划模型，松弛被嵌入在一个启发式算法中，该算法广泛地使用了与拉格朗日乘子相关的对偶信息。在其他情况下，Barrena等 [11] 和Fischer等 [19] 可以分别使用该算法来提高分支–切割算法和切割平面方法的效率。此外，许多启发式算法被应用于列车运行图问题 [20] [21] 。我们试图强调目前存在的各种方法，并且Fisher等 [22] 提到关于这类问题优化方法的更全面的调查。

本文的结构做如下安排。在第2节中，提出了一个单线轨道区间上的列车运行图优化模型。在第3节中，利用包括拉格朗日松弛法在内的分解技术构造了其上界。接下来，我们将某些约束对偶化为目标函数，而松弛问题可以分解为两组具有附加相同约束的子问题。其次，我们应用启发式方法，通过解决一组拉格朗日松弛解，进一步得到了一个下界。考虑到优化问题的大小，计算时间适中。在第4节中，我们为其生成了实际的实例，并进行了数值实验。该方法可以提供可行解，生成下界来评估所提方法的有效性。在第5节中，本文展示了研究的结论和未来潜在的研究成果。

2. 数学模型

我们考虑的列车运行图问题中，如图3所示，有两个车站——车站1和车站2，一条连接车站1与车站2的单线轨道区间铁路。记E表示从车站1出发且向东行驶的列车集合；记W表示从车站2出发且向西行驶的列车集合。设时间集合为 $T_m=\left\{1,2,\cdots ,m\right\}$ 。为了方便讨论，我们将上述所有列车分为慢、快两种，分别记 $S_1$ 为慢车集合，对应通过该区间时间为 $h_1$ ； $S_2$ 为快车集合，对应通过该区间时间为 $h_2$ 。定义决策变量为

$x_t^n\in \left\{0,1\right\}$ (1)

Figure 3. A stretch of single track railway

图3. 一段单线轨道铁路

其中， $i=1,2,\cdots ,t$ 表示时刻，n表示列车。车如果列车n在时刻t从车站i出发，则 $x_t^n=1$ ；否则 $x_t^n=0$ 。如图4所示，图中的斜线表示列车。其中，斜线1与斜线2和斜线3方向相反，它们相交表示在某时间的同一区间内冲突；斜线4与斜线5表示方向相同且在某时间的同一区间内冲突；斜线6与斜线7重合表示它们始终占用同一段区间。我们的模型应该能够保证出发列车不会造成上述冲突。

Figure 4. Conflicts that may occur

图4. 可能发的冲突

首先，我们考虑在某个时刻从车站1和车站2离开的列车不能超过一列。否则，在未来的某个时间，一个区间将会有不止一列列车。换句话说，变量需要满足

${\displaystyle \underset{r\in E}{\sum }{x}_t^r}+{\displaystyle \underset{l\in W}{\sum }{x}_t^l}\le 1,\forall t\in T_m$ (2)

如图5所示，是当h₁ = 3、h₂ = 2时的列车调度示意图。当调度开往车站2方向的某辆慢车在时刻t进入轨道，则在后续间隔的1~2个单位时间内，不允许由反方向列车进入；往前1个单位时间内，不允许反方向进入；往前2单位时间内不允许反方慢车向列车进入。如图6所示，当调度开往车站2方向的某辆快车在时刻t进入轨道，则在后续间隔的1个单位时间内，不允许由反方向列车进入；往前1个单位时间内，不允许反方向进入；往前2单位时间内不允许反方慢车向列车进入。对于同向列车来说，实际上无需考虑不同车速，只要满足式(1)的约束，就不会有冲突。在我们的模型中，可以允许多辆列车同时到达终点，因为终点站通常有很多股道。开往车站1的列车也类似讨论。

Figure 5. Conflicts that may occur with a slow train

图5. 某辆慢车可能发生的冲突

Figure 6. Conflicts that may occur with a express train

图6. 某辆快车可能发生的冲突

其次，根据上述分析，我们考虑如下解决冲突的主要约束

${\displaystyle \underset{l\in S_1}{\sum }{\displaystyle \underset{k=1}{\overset{h_1-1}{\sum }}{x}_{t-k}^l}}+{\displaystyle \underset{r\in S_1}{\sum }{x}_t^r}+{\displaystyle \underset{l\in W}{\sum }{\displaystyle \underset{k=1}{\overset{h_1-1}{\sum }}{x}_{t+k}^l}}\le 1,\forall t\in T_m$ (3)

${\displaystyle \underset{r\in S_1}{\sum }{\displaystyle \underset{k=1}{\overset{h_1-1}{\sum }}{x}_{t-k}^r}}+{\displaystyle \underset{l\in S_1}{\sum }{x}_t^l}+{\displaystyle \underset{r\in E}{\sum }{\displaystyle \underset{k=1}{\overset{h_1-1}{\sum }}{x}_{t+k}^r}}\le 1,\forall t\in T_m$ (4)

${\displaystyle \underset{l\in S_1}{\sum }{\displaystyle \underset{k=1}{\overset{h_1-1}{\sum }}{x}_{t-k}^l}}+{\displaystyle \underset{l\in W}{\sum }{\displaystyle \underset{k=1}{\overset{h_2-1}{\sum }}{x}_{t-k}^l}}+{\displaystyle \underset{r\in S_1}{\sum }{x}_t^r}+{\displaystyle \underset{l\in W}{\sum }{\displaystyle \underset{k=1}{\overset{h_1-1}{\sum }}{x}_{t+k}^l}}\le 1,\forall t\in T_m$ (5)

${\displaystyle \underset{l\in S_1}{\sum }{\displaystyle \underset{k=1}{\overset{h_1-1}{\sum }}{x}_{t-k}^l}}+{\displaystyle \underset{l\in W}{\sum }{\displaystyle \underset{k=1}{\overset{h_2-1}{\sum }}{x}_{t-k}^l}}+{\displaystyle \underset{r\in S_2}{\sum }{x}_t^r}+{\displaystyle \underset{l\in W}{\sum }{\displaystyle \underset{k=1}{\overset{h_2-1}{\sum }}{x}_{t+k}^l}}\le 1,\forall t\in T_m$ (6)

${\displaystyle \underset{r\in S_1}{\sum }{\displaystyle \underset{k=1}{\overset{h_1-1}{\sum }}{x}_{t-k}^r}}+{\displaystyle \underset{r\in E}{\sum }{\displaystyle \underset{k=1}{\overset{h_2-1}{\sum }}{x}_{t-k}^r}}+{\displaystyle \underset{l\in S_1}{\sum }{x}_t^l}+{\displaystyle \underset{r\in E}{\sum }{\displaystyle \underset{k=1}{\overset{h_1-1}{\sum }}{x}_{t+k}^r}}\le 1,\forall t\in T_m$ (7)

${\displaystyle \underset{r\in S_1}{\sum }{\displaystyle \underset{k=1}{\overset{h_1-1}{\sum }}{x}_{t-k}^r}}+{\displaystyle \underset{r\in E}{\sum }{\displaystyle \underset{k=1}{\overset{h_2-1}{\sum }}{x}_{t-k}^r}}+{\displaystyle \underset{l\in S_2}{\sum }{x}_t^l}+{\displaystyle \underset{r\in E}{\sum }{\displaystyle \underset{k=1}{\overset{h_2-1}{\sum }}{x}_{t+k}^r}}\le 1,\forall t\in T_m$ (8)

其中，当 $h_2=1<h_1$ 时，约束(3)~(4)起作用；当 $1<h_2<h_1$ 时，约束(5)~(8)起作用。

最后，我们的目的是找到一个能使利润最大化的调度。在这种情况下，每组列车都有一个最佳时间，当列车在最佳时间出发时获得最大的利润，出发时间早于或晚于最佳时间将会减少利润。一般来说，实际出发时间与最优时间的差异越大，利润就越小。然而，这种差异和利润之间的关系通常不是线性的。设 $D\left(r\right)$ 表示在车站1发车的时刻对所有 $r\in E$ ， $D\left(l\right)$ 表示在车站2发车的时刻对所有 $l\in W$ ， $V(*)$ 表示列车发车所获利润。基于这样的假设所构建的模型，规划的目标是总利润最大，表达式为

$\max {\displaystyle \underset{r\in E}{\sum }V\left[D\left(r\right)\right]}+{\displaystyle \underset{l\in W}{\sum }V\left[D\left(l\right)\right]}$ (9)

综上所述，式(1)~(9)建立了关于单线区间上列车运行图优化的0~1数学规划模型，记为P。其中式(9)为优化目标，式(1)~(8)为约束条件。

3. 拉格朗日松弛算法

本节将介绍用于解决第2节所提到优化模型P的方法——拉格朗日松弛算法。在3.1中，我们叙述该方法的整体思路和框架；在3.2中，我们将详细介绍该算法的细节，其中包括松弛问题的构造、原问题P的启发式算法和求解拉格朗日界的次梯度优化。

3.1. 框架概述

拉格朗日松弛算法的关键组成部分是松弛约束的选择、可行解生成和乘子更新方法。为简单起见，设所要研究的数学规划P的目标函数值 $\psi$ 为，其对应拉格朗日松弛目标函数值为 $\phi$ 。给定最大化问题中的乘子，松弛问题对原问题有一个上界，即对所有可行解有 $\psi \le \phi$ 成立。拉格朗日松弛算法设计思想如下：松弛约束的选择是实现拉格朗日松弛技术的一个重要步骤。基本上，该方法通过引入拉格朗日乘子来放宽耦合约束，将一个大规模的问题分解为子问题。我们首先求解拉格朗日松弛问。然而，如果放宽了不同的约束条件，所得到的模型可能难以解决。因此，我们稍后将研究约束放松的策略(见3.2.1节)。其次，我们将讨论原始问题的拉格朗日启发式计算，以找到一个可行的解，例如通过调整松弛解来获取。此后，我们更新了乘子，以更接近解决方案。目前，第k次迭代为止的最佳解决方案的目标值 ${\varphi }_{best}^k=\min \left\{{\phi }^s|s=1,\cdots ,k\right\}$ 提供了一个原问题上限的最优值，和当前最好的可行的解决方案目标值 ${\psi }_{best}^k=\min \left\{{\psi }^s|s=1,\cdots ,k\right\}$ 提供了一个下界。因此，如果 ${\phi }_{best}^k-{\psi }_{best}^k$ 足够小或达到最大迭代次数k_max，我们可以终止迭代过程。由于原始启发式可能无法找到最优解，当该算法终止时，存在对偶性间隙。我们根据式(10)来计算一个相对间隙。为了提供了一个完整的拉格朗日松弛算法以求解铁路路段通行能力计算的整数规划模型。我们将上述思想概括成如图7的程序框图。

$\epsilon =\frac{{\varphi }_{best}^k-{\psi }_{best}^k}{{\psi }_{best}^k}\times 100\%$ (10)

Figure 7. Lagrange relaxation algorithm framework

图7. 拉格朗日松弛算法框架

3.2. 算法细节

3.2.1. 松弛问题的构造

对于一个给定的整数规划问题，线性松驰是去掉模型中的整数约束，而拉格朗日松驰则是放松模型中的部分线性约束，保留整数约束和其他的线性约束。这些被松驰的约束并不是被完全去掉，而是在利用拉格明日乘子在目标函数上增加相应的惩罚项，对不满足这些约束条件的解进行较大的惩罚。本文模型P的块结构启发了我们用拉格朗日松弛法来解决这个问题。我们将该优化问题中约束条件(2)被对偶化为具有拉格朗日乘子的序列的目标函数，得到了式(11)的拉格朗日松弛问题，记为RP。

$\begin{array}{l}\text{min}{\displaystyle \underset{r\in E}{\sum }V\left[D\left(r\right)\right]}+{\displaystyle \underset{l\in W}{\sum }V\left[D\left(l\right)\right]}+\lambda \left(1-{\displaystyle \underset{r\in E}{\sum }{x}_t^r}-{\displaystyle \underset{l\in W}{\sum }{x}_t^l}\right)\\ \text{s}\text{.t}\text{.}\lambda \ge \text{0},式\text{(1)},式\text{(3)~(8)}\end{array}$ (11)

一般来说，不可能保证找到的使原问题目标函数值等于拉格朗日松弛问题目标函数值，但这种情况经常发生在特定的问题实例中。实际上，原问题目标函数值 ≤ 拉格朗日松弛问题目标函数值允许使用拉格朗日松弛问题代替线性松弛问题来提供一个原问题的下界在分支–定界算法中。虽然这是拉格朗日松弛最明显的用途，但它还有许多其他用途。我们将在3.2.2节中讨论确定的方法。它可以作为选择分支变量和选择下一个分支要探索的媒介。通过扰动其相近可行解，通常可以得到原问题的良好可行解。

3.2.2. 原问题P的启发式算法

本小节将讨论如何使用松弛解来获得原问题的可行解。现实情况是，基于现有的电脑效率情况，我们能够接受通过使用一些启发式的形式去寻求大型离散优化间题的解决方案，如果剩余约束包含一些不等式。P的解可以是最优的，但不可行。我们希望修改可行解，以满足这些约束条件。因此，拉格朗日松弛算法需要一个启发式算法来寻找良好的可行解。许多启发式算法选择沿用精确解法，不过这些算法会在到达规定时间后停止搜索，从发现的可行方案中选取最佳的那个。随着分支定界、分支剪切以及大规模计算分析技术的发展，越来越多的大型的实例可以通过这种方式成功解决。尽管如此，这些问题的复杂度的指数增长最终会导致这些问题难以被如上的启发式算解决。

我们针对本文P模型提出的启发式算法，不再追求和精确方法表面上相似。相反，这些方法充分探求了有效的问题结构，寻找能够得到最优解的最有机会、最直观的驱动，从而有希望获得好的可行的解。利用这些方法获得的结果可能确实很好，但是通常它们的优劣需要经过大量的实例评估。这种严格的启发式方法很少会根据离数学上能够证明的离最优解的远近而决定是否停止运行。在过去研究报中，根据使用轨道的优先级列表，依次调度列车。本文利用一个贪婪启发式的方法构造了原问题的可行解。在P中，为了完成一个可行的解决方案的构建，只需要将冲突列车分配给空闲时间。现在，我们很快就想到了一些需要解决的两个问题。即：1) 我们如何解决相互冲突的列车？2) 如何调整剩余的列车？首先，获得最多收益的火车将首先被派遣。在当前可用时刻，其余列车将以高值优先于低值排列。我们设置了两条时间线，以帮助记录当前可用的时刻。我们的方法背后的一般思想现在可以总结为以下算法1的形式。

3.2.3. 求解拉格朗日界的次梯度优化

目前，我们介绍了拉格朗日松弛问题的主要思想及其原理，下面讨论如何选择合适的拉格朗日乘子，提高拉格朗日松弛算法求解P的效率，并分析其给出的P目标函数值的上下界。在RP中，记

$\phi \left(\lambda \right)={\displaystyle \underset{r\in E}{\sum }V\left[D\left(r\right)\right]}+{\displaystyle \underset{l\in W}{\sum }V\left[D\left(l\right)\right]}+\lambda \left(1-{\displaystyle \underset{r\in E}{\sum }{x}_t^r}-{\displaystyle \underset{l\in W}{\sum }{x}_t^l}\right)$ (12)

得到了对偶目标函数值，从而得到了相应的拉格朗日对偶问题式(13)。

$\begin{array}{l}\text{min}\phi \left(\lambda \right)\\ \text{s}\text{.t}\text{.}\lambda \ge \text{0}\end{array}$ (13)

现在，我们对式(13)中的拉格朗日对偶目标函数做一个简单解释。如果拉格朗日松弛问题的可行域是有界的，那么松弛问题的最优目标函数值将是关于拉格朗日乘子的分段线性凸函数，该函数对应着最大的以所有可行解为参数的关于乘子的线性函数。该函数之所以是分段的，原因在于对于一段连续的不同乘子，对应的最优解可能是相同的；而当乘子取到其他值时，可能得到和之前完全不同的最优解，此时最优目标函数就会发生变化，从而其是不可微的。

若一个优化问题，其目标函数可导，梯度方向(目标函数关于决策变量的偏导向量)对应着解的改进方向。例如在某一个最优解处，拉格朗日对偶问题目标函数由对偶面唯一决定，因为对应的拉格朗日松弛问题存在唯一的最优解。在此处，梯度存在并且得到改进方向。但是，对于不可微的目标函数，拉格朗日对偶问题搜索最优乘子的难点出现在对偶函数分段线性的分段交界处，此时，拉格朗日松弛问题将会存在两个最优解。为了解决以上问题，可以对梯度的概念进行推广得到次梯度，它是一个用所有可能的梯度方向的凸组合来表达的方向。

对于一个给定的整数规划模型和拉格朗日松弛模型，存在许多改进搜索算法努力尝试找到最优或近似最优的的原问题目标函数的界。本小节介绍两种次梯度搜索优化算法，即标准次梯度法和修正次梯度法，这些算法是最速下降梯度算法的推广。拉格朗日对偶问题的次梯度搜索算法的思想如下：给定一个整数规划问题，拉格朗日分解策略选择最佳松驰问題时的目标是使得该松弛问题的最优目标函数值是原问题的最好的界。因此，在选择最优乘子的过程中，该策略不断求解乘子对应的拉格朗日松弛问题，并根据得到的结果判断是否能一步改进乘子的选择。如果可以，则更新乘子，并重新计算对应的松弛问题；否则，则终止搜索过程。

首先介绍一种基于迭代公式来优化非光滑目标函数的经典有效方法——标准次梯度法，该算法的具体细节如下。

第一，拉格朗日乘子的初始值 ${\lambda }^0$ 被设置为0。基于这个初始值，将使用式(14)确定下一个 $\lambda$ 的值。其中， ${\gamma }_k$ 和 $g_k$ 分别表示第k次迭代时的步长和次梯度。

${\lambda }^{k+1}={\lambda }^k-{\gamma }_k\cdot g_k$ (14)

第二，在当前拉格朗日对偶问题的解 ${\lambda }^k$ 处，该算法虽然可以取到所有次梯，即松弛问题最优解不满足的被松弛约束。例如，在式(15)中， $g_k$ 这个向量是 $\phi$ 在 ${\lambda }_k$ 处的次梯度。因为该方法易于编程，且效果良好在许多实际问题上，它已成为最多的问题的流行方法。

$g_k=\left[1-{\displaystyle \underset{r\in E}{\sum }{\hat{x}}_t^r\left(\lambda \right)}-{\displaystyle \underset{l\in W}{\sum }{\hat{x}}_t^l\left(\lambda \right)}\right]$ (15)

第三，如果步长收敛到0，但是步长之和并不收敛到0，那么沿着次梯度方向的搜索步数也将收敛。在我们的例子中，步长参数将由式(16)更新。其中， $\phi \left({\lambda }_k\right)$ 是 $\phi$ 在k次迭代的值。在我们的实施过程中，由启发式方法生成一个可行的解来代替最优解 ${\bar{\phi }}_k$ 。 $\Vert *\Vert$ 表示普通的欧几里得范数。

${\gamma }_k=\frac{\phi \left({\lambda }^k\right)-{\bar{\phi }}_k}{{\Vert g_k\Vert }^2}$ (16)

定理 设任意步长序列 $\left\{{\lambda }_k\right\}$ ， $\left\{{\gamma }_k\right\}$ 是由该规则 ${\lambda }_{k+1}={\lambda }_k-{\gamma }_k\cdot g_k$ 生成的。对于无约束最优化问题，当目标函数 $\phi \left(\lambda \right)$ 满足：

• $\phi \left(\lambda \right)$ 是凸函数，

• $\phi \left(\lambda \right)$ 至少存在一极小点 ${\mu }^{*}$ ，即 ${\phi }_D=\phi \left({\mu }^{*}\right)>-\infty$ ，

• $\phi \left(\lambda \right)$ 满足利普希茨条件，即 $\left|\phi \left(x\right)-\phi \left(y\right)\right|\le G{\Vert x-y\Vert }_2，\forall x，y\in R^n$ ，G是一个常数，

• ${\gamma }_k\to 0，{\displaystyle \underset{i=0}{\overset{k}{\sum }}{\gamma }_i}=\infty 和{\displaystyle \underset{i=0}{\overset{k}{\sum }}{\gamma }_i^2}=\infty$ ，

则对任意 $k\ge 0$ 我们有 ${\phi }_{best}^k\to {\phi }_D$ 成立，其中 ${\phi }_D=\min \left\{\phi \left(\lambda \right)|\lambda \ge 0\right\}$ 。

证明 $\forall k\ge 0$ ， ${\lambda }_{k+1}={\lambda }_k-{\gamma }_k\cdot g_k\left({\gamma }_k\ge 0\right)$ 。现在我们考虑等式(17)

${\Vert {\lambda }_{k+1}-{\lambda }^{*}\Vert }_2^2={\Vert {\lambda }_k-{\lambda }^{*}\Vert }_2^2-2{\gamma }_k\left({\lambda }_k-{\lambda }^{*}\right)+{\gamma }_k^2{\Vert g_k\Vert }_2^2$ (17)

由凸函数的定义我们有式(18)成立

${\Vert {\lambda }_{k+1}-{\lambda }^{*}\Vert }_2^2-{\Vert {\lambda }_k-{\lambda }^{*}\Vert }_2^2\le -2{\gamma }_k\left[\phi \left({\lambda }_k\right)-\phi \left({\lambda }^{*}\right)\right]+{\gamma }_k^2{\Vert g_k\Vert }_2^2$ (18)

对任意k，显然将式(18)求和可以得到式(19)

${\displaystyle \underset{i=1}{\overset{k}{\sum }}\left({\Vert {\lambda }_{k+1}-{\lambda }^{*}\Vert }_2^2-{\Vert {\lambda }_k-{\lambda }^{*}\Vert }_2^2\right)}={\Vert {\lambda }_{k+1}-{\lambda }^{*}\Vert }_2^2-{\Vert {\lambda }_1-{\lambda }^{*}\Vert }_2^2\le -2{\displaystyle \underset{i=1}{\overset{k}{\sum }}{\gamma }_i\left[\phi \left({\lambda }_i\right)-\phi \left({\lambda }^{*}\right)\right]}+{\displaystyle \underset{i=1}{\overset{k}{\sum }}{\gamma }_i^2{\Vert g_i\Vert }_2^2}$ (19)

因为对任意k， ${\Vert {\lambda }_{k+1}-{\lambda }^{*}\Vert }_2^2$ 是大于零的，故重新改写式(19)，显然有式(20)成立

$2{\displaystyle \underset{i=1}{\overset{k}{\sum }}{\gamma }_i\left[\phi \left({\lambda }_i\right)-\phi \left({\lambda }^{*}\right)\right]}\le {\Vert {\lambda }_1-{\lambda }^{*}\Vert }_2^2+{\displaystyle \underset{i=1}{\overset{k}{\sum }}{\gamma }_i^2{\Vert g_i\Vert }_2^2}$ (20)

我们考虑最后一部分式(21)，其中第一个不等式是由于在第k，次迭代时的最优解，即 ${\phi }_{best}^k\le \phi \left({\lambda }_i\right),$ $i=1,\cdots ,k;$ 第二个不等式是因为目标函数的次梯度满足利普希茨条件

${\phi }_{best}^k-{\phi }_D\le \frac{{\Vert {\lambda }_1-{\lambda }^{*}\Vert }_2^2}{2{\displaystyle \underset{i=1}{\overset{k}{\sum }}{\gamma }_i}}+\frac{{\displaystyle \underset{i=1}{\overset{k}{\sum }}{\gamma }_i^2{\Vert g_i\Vert }_2^2}}{2{\displaystyle \underset{i=1}{\overset{k}{\sum }}{\gamma }_i}}\le \frac{R^2}{2{\displaystyle \underset{i=1}{\overset{k}{\sum }}{\gamma }_i}}+\frac{G^2{\displaystyle \underset{i=1}{\overset{k}{\sum }}{\gamma }_i^2}}{2{\displaystyle \underset{i=1}{\overset{k}{\sum }}{\gamma }_i}}$ (21)

当 ${\gamma }_k\to 0,{\displaystyle \underset{i=0}{\overset{k}{\sum }}{\gamma }_i}\to \infty$ 和 ${\displaystyle \underset{i=0}{\overset{k}{\sum }}{\gamma }_i^2}=\infty$ 成立时，一定有式(22)成立。故原命题成立。

$\frac{R^2}{2{\displaystyle \underset{i=1}{\overset{k}{\sum }}{\gamma }_i}}\to 0,\frac{G^2{\displaystyle \underset{i=1}{\overset{k}{\sum }}{\gamma }_i^2}}{2{\displaystyle \underset{i=1}{\overset{k}{\sum }}{\gamma }_i}}\to 0$ (22)

£

在定理中，以迭代格式进行更新的次梯度在理论上是收敛的。其次，本文给出了一种式(23)的有效算法——修正次梯度法。Kasimbeyli等 [23] 证明，如果 ${\bar{\phi }}_k$ 等于最优对偶值，则通过该修正方法得到的 ${\mu }_{k+1}$ 比通过式(14)得到的 ${\mu }_{k+1}$ 更接近最优集。在上述方法中，像这样的次梯度算法是发展曲折的。式(23)中的目标值 ${\bar{\phi }}_k$ 根据 ${\bar{\phi }}_k=\left(1-{\theta }_k\right){\phi }_{best}^k+{\theta }_k{\psi }_{best}^k$ ， $0\le {\theta }_k\le 1$ 进行更新。通过测试，我们实现了算法的收敛性和精度之间的权衡。

$\begin{array}{l}\text{min}{\Vert \lambda -{\lambda }_k\Vert }_2^2\\ \text{s}\text{.t}\text{.}{\bar{\phi }}_k\ge \phi \left({\lambda }_k\right)+{\text{g}}_k\left(\lambda -{\lambda }_k\right)\text{,}\lambda \ge 0\end{array}$ (23)

4. 数值实验与分析

本节中，我们将对P模型进行数值实验。其中，约束为避免对向列车冲突，目标是与列车发车时间有关的利润函数。对于该模型，我们利用拉格朗日贪婪启发式算法进行求解采用贪婪的启发式拉格朗日算法对其进行数值实验。

4.1. 参数设置和数据集生成

在拉格朗日松弛算法中，停止准则是将对偶性间隙设置为1%或者最大迭代次数k_max被设置为30。基于计算技术和信息技术研制的各类软件求解工具在科研领域得到越来越广泛的应用，目前优化求解器包括Gurobi、LINGO、GAMS和MATLAB……等等。MATLAB中的优化工具箱对目标函数梯度的相关参考可通过求解器估计或由用户向程序提供。其支持CFOR-TRAN、C++、JAVA等语言的调用。通过启用内置的并行算法或通过自定义的并行算法实现，从而缩短求解时间。我们应用MATLAB通过调用非线性规划求解器包来实现算法过程。

我们从太原铁路局车务段的日班计划中收集了一组货运列车数据，如表1所示。被测试的例子将安排30列火车，其中1号列车~17号列从车站1开往车站2，13号列车~30号列车从车站2开往车站1。其中：t_OPTI表示在该时刻发车或到站可获得最大利润的最优时间；t_min和t_max表示允发车的最大时间差；v_dep表示在t_OPTI时刻获得的利润。如果所有列车都以最高利润时刻发车，我们可以得到当前的最高总利润51。然而，它实际上是不可能获得的。因为它会造成一些冲突，从而导致带来大量的经济损失。由第3节中的描述我们可知，本文的NLP问题主要体现在优化目标函数上。

Table 1. Test train No.1~30

表1. 测试数据列车1号~30号

如图8所示，在非平行运行图的数据中，相反方向列车有13对发生冲突，而同一方向的列车存在共同竞争同一轨道的情况，如较粗的黑线表示不止一辆列车占用该轨道。

Figure 8. Conflicts that may occur

图8. 无约束列车运行图

4.2. 计算结果与分析

对于本文输入的非线性0~1数学数规划例子，我们用第4节中的次梯度对偶迭代算法进行了检验。图9、图10分别采用标准次梯度法和修正次梯度法对第3节的优化模型求解得到的调度列车运行图。经过实验调度与图8中列车无约束最优调度列车运行图有明显的对比，即在尽可能减小损失的情况下解决了所有列车的冲突。图9和图10分别获得利润为26.05和26.75。

Figure 9. The result of standard sub-gradient dual iterative

图9. 标准次梯度对偶迭代结果

Figure 10. The result of modified sub-gradient dual iteration

图10. 修正次梯度对偶迭代结果

Figure 11. The best NS-TWDP value and dual value for two dual iterations

图11. 两个对偶迭代的最佳 NS-TWDP 值和对偶值

上述我们给出了列车运行图优化模型的直观调度结果，但并不能很清晰地得到一些算法结论。现在，我们从数值方面的结果观察一些结论。首先，在图11中，我们以P为例展示了在测试示例中给出的最佳原始目标函数值和对偶值。横轴表示迭代次数，纵轴表示当前迭代次数下所有列车的利润值。上曲线对应于拉格朗日松弛结果的值，这是原问题的上界；下曲线表示原问题的可行解的值，是原问题的下界。由该图可知，修正次梯度法的收敛速度比标准次梯度法快，这可能是由于标准次梯度迭代格式和步长的选择所致。在P模型的计算结果中，虽然对偶间隙能反应一定精度问题, 但不能完全参考。标准次梯度的对偶间隙为(33.61 − 26.05)/26.05 = 29%，修正次梯度的对偶间隙为(26.99 − 26.75)/26.75 = 0.8%，两者的精度虽然相差接近28%，但最终实际获得的利润26.75 − 26.05 = 0.7缺相差不大。综上所述的实验结果，标准次梯度在收敛速度和精度方面优于校正后的次梯度，与之前研究的文献结果一致 [4] [24] 。

5. 结论

本文通过分析单线轨道区间上如何找到最佳调度的列车运行图问题，建立了该问题的数学规划模型。 在实际情景中，一方面由于列车数量较多使得模型中决策变量与约束条件的规模较大，另一方面为了不影响车流通畅，可用于计算的时间有限。综合上述考虑，本文给出了一种计算模型的拉格朗日启发式算法。数值实验的结果表明，使用拉格朗日松弛算法，对一个日班计划中需要被规划的列车可在合理的计算范围时间内较快得到结果；并通过启发式方法找到最优目标函数在最坏情况下的上界或下界。即使对偶间隙为29%，但实际利用修正次梯度迭代得到的上下界差33.61 − 26.05 = 7.65符合实际预设效果。未来我们着眼于可能将这些问题扩展到单轨线路的网络情况，合并到我们的模型中。而且可以通过某种形式的拉格朗日分解将网络分解为单轨线，然后使用本文的技术来实现。

致谢

对提供指导和帮助的梁东岳老师、给予转载和引用权的资料、文献所有者，表示感谢。