外文翻译

利用混合整数线性规划求解固定与移动闭塞条件下的多列车操纵优化问题

王义惠

摘要：

本文主要解决基于固定闭塞与移动闭塞条件下的多列车操纵优化问题，利用婪式法以及同步法进行求解。对于两种求解方法，文中都将列车操纵优化问题转化为混合整数规划问题。特别地，文中将能量消耗目标函数以及列车模型利用分段仿射模型进行近似处理。在模型建立过程中，考虑了线路基本阻力、分段速度限制和最大牵引力等影响因素。并且，文中将固定闭塞与移动闭塞条件下的前车运行进行离散化，利用分段仿射处理转化为混合整数规划问题。利用同步法与婪式法进行求解，发现同步法的仿真结果的控制更加优化，但是需要花费更多的计算时间。并且将两种方法的仿真结果与伪谱法的仿真结果进行比较，发现伪谱法在求解优化控制中表现更佳。主要表现为婪式混合整数规划法的求解结果比婪式伪谱法的求解结果在节能与准时方面表现更差，但是求解时间呈现一至两个数量级的差异，婪式伪谱法需要花费更多的时间进行求解。同样的，同步混合整数规划法与同步伪谱法存在相同的差异性。

关键词：运行曲线，列车操纵，信号系统，混合整数规划，伪谱法

1.简介

随着能源价格的增长与环境问题的严重化，交通运输系统的节能变得越来越重要。轨道交通因其安全，快捷，准时与舒适的服务，在交通运输系统可持续发展中扮演着至关重要的作用。列车运行能耗占到中国铁路年运营与维护能耗的13%到16%，因此将其作为主要优化目标，哪怕是小小的一个提升，也会为铁路系统节约一大笔资金。为了达到这个目的，一些辅助驾驶系统已经被开发，例如FreightMiser、Metromiser以及driving style manager。随着铁路系统的现代化，自动驾驶系统在确保准确停车，准点，节能和驾驶舒适方面起到越来越重要的作用。铁路控制中心或者自动化控制系统主要依据列车监控系统提供的线路基本阻力、速度限制、最大牵引力和制动力等数据来决定列车运行曲线。

在参考文献中，对于列车操纵优化的研究开始于1960年代，主要致力于解决列车从一个站到相邻站的运行曲线。自从发现该研究对于列车运行节能有帮助，多种方法

被应用到这个研究中，主要分成两种类别：解析法与数值优化。对于解析法，Howlett 等人将极大值原理运用到列车运行操纵优化中，并将操纵工况划分为四种情况：全力牵引，恒速运行，惰行以及全力制动。但是在该方法中，如果将更多现场状况的约束引入到模型中，形成非线性约束条件，很难得到解。随着当今计算机计算能力的提升，越来越多数值优化方法被应用到这个研究中，比如婪式法以及遗传算法等，但是通过这些方法不一定能得到优化方案。另一方面，二次多项式优化被应用到求解该问题，主要特点是通过分段仿射处理列车模型。该理论主要由Va?ak提出，在先前的文章中，我们希望将运行曲线优化问题转化为混合整数规划问题，并通过现成的求解工具进行全局求解。

但是先前提到的方法都没有考虑铁路信号系统的影响，例如固定闭塞信号系统与移动闭塞信号系统。固定闭塞信号系统是通过轨旁信号机对固定闭塞区间进行保护，列车只有在信号机显示通过状态下，才能进入对应的闭塞区间。移动闭塞信号系统需要通过列车与区间控制者之间的沟通，得知列车位置以及速度的实时数据，对于每一趟列车定义移动闭塞。与固定闭塞相比，移动闭塞情况下，运行列车运行间隔减小，因此增加了区间通过能力。Lu和Feng考虑在移动闭塞区间下，同一方向，同一轨道，前车对后车造成的约束下，对后车运行曲线进行优化，对前后车使用相同的遗传算法进行运行曲线优化，使总的能源消耗降低。Gu,Lu和Tang利用非线性优化对后车运行曲线进行优化，在前车运行或者停车的两种情况都被考虑，并且提出对应的后车运行策略。Ding等人将移动闭塞区间下多车运行的约束条件考虑到列车运行节能模型中，并且将全力牵引，惰行，全力制动三种工况结合到算法中，通过既定规则安排工况顺序，最终得到优化的运行曲线。

在本文中，移动闭塞区间以及固定闭塞区间条件下，由于前车造成的对于后车的约束进行建模，先将其离散化，再通过分段仿射处理转化为线性约束。因此，这些约束可以融入到混合整数规划模型中，并且通过比其他方法更快地进行求解。文中，运用婪式法以及同步法对问题进行求解，并且与伪谱法进行对比。在过去地十年中，伪谱法在数值优化领域被大量使用，例如轨道转移，轨迹引导，引力控制，因此文中选择其作为仿真比较。

接下来的论文内容如下：第二部分，列车模型以及混合整数规划在单车运行控制的应用基于先前论文进行说明；第三部分，介绍固定闭塞区间以及移动闭塞区间的原则；第四部分，将固定闭塞条件下，后车跟驰前车受到的约束转换到混合整数规划模型中；第五部分，介绍了移动闭塞条件下的混合整数规划模型；第六部分，利用北京亦庄的数据进行仿真计算；第七部分，进行总结与展望。

2.列车运行模型与混合整数规划模型

本部分对于单车运行优化模型建立与混合整数规划模型的运用进行总结。 2.1优化控制模型 连续质点模型经常被运用到列车运行优化控制，如下：

mρdE ?ds =u (s )?R b (v )?R l (s,v ),

dt ?ds =√?, （1）

其中，m 代表的列车单位质量，ρ代表动态质量参数，E

?代表单位质量的动能，等价于0.5v 2，v 代表列车运行速度，s 代表列车的位置，u 代表控制变量即牵引力或者制动力，上界为最大牵引力u max ，下界为最大制动力u min ，因此u min

单位质量的动能E

?=0.5v 2和时间t 是相对固定的，位置参数s 在模型中是固定参数。运行曲线优化问题可以建立如下：

J =∫max?(0，u (s ))s

end s start ds （2） s.t.

u min ≤u (s )≤u max ，

0

?(s )≤E ?max (s )， E

?(s start )=E ?start ，E ?(s end )=E ?end ， t (s start )=0，t (s end )=T ，

其中，目标函数J 是不考虑再生制动的列车运行能量消耗，E ?(s )等价于0.5V m ax 2(s )，

V m ax (s )代表最大运行速度，主要取决于列车参数以及线路状况，是关于位置s 的分段

常数函数，s start ，E

?(s start )，t (s start )分别代表行程的起始位置，起始单位动能，起始时刻，s end ，E

?(s end )，t (s end )分别代表行程的结束位置，结束单位动能，结束时刻，T 代表计划运行时间又时刻表或者重订时刻表给定。文中假设E ?(s )>0，即列车速度大

于零，在路途中不停车，这个假设在现实中可以通过定义起车与停车为极小值来实现，而且列车运行管理系统也不希望列车在路途中途停车。

2.2转化规则

本小节将介绍通过Williams 提出的转化规则将非线性优化控制问题转化混合整数规划问题。首先将非线性的列车模型转化为混合逻辑模型，然后将优化控制问题的约束条件利用分段仿射处理转化为混合整数线性约束，最终形成混合整数规划问题。

考虑声明f ?(x ?)≤0，f ?：R n →R 表示仿射，x ?∈X ，X ∈R n ，并且使

M ?=max x ?∈X f ?(x ?)， m ?=min x

?∈X f ?(x ?)

如果我们引入逻辑变量δ∈{0，1}，那么可以得到如下等式：

当且仅当{f ?(x ?)≤M ?(1?δ)f ?(x ?)≥ε+(m ??ε)δ

成立时，[f ?(x ?)≤0]?[δ=1]成立 其中，ε是一个小正整数（典型的机械精度），引入其将严格等式转化为非严格等式，适合混合整数规划问题的框架。

两个逻辑变量的乘积δ1δ2，可用另一个辅助逻辑变量δ3代替，也就是说

[δ3=1]?[δ1=1]∧[δ2=1]

等价于

{?δ1+δ3≤0,

?δ2+δ3≤0,δ1+δ2?δ3≤1.

乘积δf ?(x ?)可用辅助变量z =δf ?(x ?)代替，并且满足[δ=0]?[z =0]以及[δ=

1]?[z =f ?(x ?)]，因此z =δf ?(x ?)等价于 {

z ≤M ?δ,z ≥m ?δ,z ≤f ?(x ?)?m ?(1?δ),z ≥f ?(x ?)?M ?(1?δ).

不等式集（5）-（7）均基于f ?为仿射。 2.3混合整数规划方法

在Wang 的论文中将列车运行的空间连续模型离散化，将起点站与终点站之间的位置区段[S start ,S end ]不包含中间站离散化成N 个间隔，并且假设线路，列车的特性 以及牵引力或者制动力为恒量在每一个区段[s k ,s k+1]，区段长度?s k =s k+1?s k ，k =1,2,….,N 。然后将空间离散模型通过PWA 公式转化为分段仿射模型，并且利用

2.2中描述的转化规则将分段仿射模型建立成为混合逻辑动态模型，如下：

x (k +1)=A k x (k )+B k u (k )+C 1,k δ(k )+C 2,k δ(k +1)

+D 1,k z (k )+D 2,k z (k +1)+e k , （8）

R 1,k δ(k )+R 2,k δ(k +1)+R 3,k z (k )+R 4,k z (k +1)≤R 5,k u (k )+R 6,k x (k )+R 7,k （9） 其中，x (k )=[E (k ),t (k )]T ，δ(·)和z (·)包含二元变量以及辅助变量，是由于引入转化规则。同时，等式（9）存在上下界约束变量E (k )，t (k )以及u (k )。为了方便起见，我

们用E (k )和t (k )分别标记E ?(s k )和t (s k )，混合逻辑动态模型的系数矩阵主要由列车模型，

PWA 近似值以及上下界等因素决定。

列车操作优化模型的目标为不考虑再生制动情况下的列车运行能耗，可以描述如下：

J =∑max(0,u (k ))N k=1?s k

在Wangl的论文中，引入新的辅助变量w(k)将优化控制模型重新如下建立，辅助变量w(k)的引入主要是为了转化等式（10）中的max(0,u(k))。

min

V?

C J T V?

约束于：

F1V?≤F2x(1)+f3

F4V?=F5x(1)+f6

其中，C J=[0…0??s1…?s N]T，V?=[u?T?δ?T?z?T?w?T]T，

u?=[u(1) u(2)?

u(N)]，δ?=[

δ(1)

δ(2)

?

δ(N+1)

]，z?=[

z(1)

z(2)

?

z(N+1)

],?w?=[

w(1)

w(2)

?

w(N+1)

]

进行合适的定义矩阵和矢量F1，F2，f3，F4，F5和f6，混合整数规划问题可以通过现有的商业和免费软件进行求解，例如CPLEX，Xpress-MP，GLPK。

3.区间闭塞信号系统的原则

区间闭塞信号主要是用来保持同一线路上前后相邻列车的安全距离，固定闭塞信号系统以及移动闭塞信号系统将在接下来阐述。

3.1固定闭塞信号系统

固定闭塞信号系统（FBS）在当今的列车控制系统中被普遍使用。在固定闭塞信号系统中，线路被分成不同的区段，主要依据最高列车运行速度，最差列车制动率以及信号显示的制式，信号显示的制式决定了信号可见的颜色，例如绿灯，黄灯，红灯等。每一个区段只能被一辆列车占用，该规则通过轨道电路进行检验。此外，区段通过轨边信号或者车载信号进行保护，轨边信号在铁路上的经典信号系统，车载信号被使用的越来越多，尤其是在高速铁路中，主要由于在高速下，列车司机无法有效识别轨边信号。固定闭塞信号系统主要分成两大类型，单区段信号系统和多区段信号系统。在单区段信号系统中，信号的显示根据对应区段的状况，同时需要有远信号灯为接近列车提供必要的信息。在多区段信号系统中，信号的显示需要根据对应的两个或者多个信号决定，一个简单的例子就是双区段信号系统伴随三中显示制式的信号灯，也就是红，黄和绿被称为三显式信号系统。我们将阐述三显式信号系统在固定闭塞信号系统运用中的限制，但是该方法论可以被运用到其他制式的固定闭塞信号系统中。

图一 固定闭塞信号系统下的速度编码

三显示信号系统如图一所示，在装备列车自动防护系统情况下，每一个区段在轨道电路始端携带速度编码，速度编码主要包含两部分，本区段的批准速度编码以及后区段的指导速度编码。速度编码由轨道电路控制设备进行编译，由轨道进行传输，通过车载天线传输到列车，如果列车想要进入零速区段或已被占用的区段，以及速度超过批准速度，列车会进行紧急制动。图一中显示了速度编码如何工作，当某一列车进入区段4并且接近信号机3，信号机3将会收到v max /v yellow 编码，去确定该区段的运行速度v max 以及下一区段的目标速度v yellow 。当列车进入区段2，编码变成v yellow /v min ，由于下一区段被列车1占用，因此当列车接近区段2末端时，速度需要降到v min （通常情况下v min =0）。如果列车试图进入区段1，车载系统将会检验最小速度编码（v min /v min ），并且进行紧急制动。

为了确保列车操纵不受信号系统的影响，也就是说列车操纵不受前面列车的影响，最小间隔被引入，最小间隔也就是车站发车间隔。在不被干扰的固定闭塞信号系统中，最小间隔可以被定义为：

H min,FBS =L a V max [

2+INT { L r F +

(v max F )2(2a b F )

L a } ] +v max F 2a b F +t d L +√2(L t +L s )a acc L 其中，L a 代表区段长度，L r F 表示后车在反应时间t r F 内的运行距离，反应时间主要受司机以及列车参数影响，v max F 表示后车的最高速度，a b F 表示最大制动率，t d L 表示前车的

停留时间，L t L 表示前车的车体长度，L s 表示安全间隔的长度，a acc L 表示前车的加速度。 3.2移动闭塞信号系统

随着铁路系统列车运行密度的增加，固定闭塞信号系统下的运输能力无法符合显示需求，即使固定闭塞信号系统下的运输能力可以通过缩短区段长度来增加，但是信号设备以及线路设备的安装使用费用并不合理。因此，移动闭塞信号系统（MBS ）由于其良好的表现被采用。

移动闭塞信号系统依赖于列车于区域控制室之间连续双向的联系，区域控制室计算通过计算区域内列车的移动权限限制来确保列车之间保持安全距离，

更加详细地阐

述，移动权限限制包含前车的车尾以及安全距离，也就是后车允许运行的位置。在文献（Pearson ，1973）中阐述了4种移动闭塞信号系统：空间移动闭塞信号系统，时间移动闭塞信号系统，纯粹的移动闭塞信号系统，相关的移动闭塞信号系统。Takeuchi （2003）等人评估了前三种系统，并且将其与固定闭塞信号系统基于稳定性和扰动性进行比较，纯粹的移动闭塞信号系统表现最佳。同时，其表示相关的移动闭塞信号系统没有被铁路系统所使用，尽管其经常在公路系统中被使用。因此，本文将围绕纯粹的移动闭塞信号系统进行展开，相关的方法论可以在其他移动闭塞信号系统中使用。同时，纯粹的移动闭塞信号系统是现已使用的移动闭塞信号系统的基础，其原理如下图二所示。

图二 纯粹的移动闭塞信号系统的原则

在纯粹的移动闭塞信号系统中，列车最小追踪间隔是基于后车的瞬间制动距离加上一个安全距离，利用最小追踪间隔可以避免前车突然停车时，后车对其的冲撞。在现实中，最小追踪间隔由于司机反应时间或者列车反应时间会变得更大，并且需要考虑列车长度。因此，实际移动闭塞信号系统下，最小追踪间隔计算公式如下： s L (t )?s F (t )≥L r F +(v F (t ))

2(2a b F )+S SM +L t L

其中，s L (t )和s F (t )分别代表前后车在时间t 的位置，v F (t )代表后车的速度，a b F 代表最

大制动率，S SM 代表安全距离，L r F 代表在司机或列车反应时间t r F 内的后车运行距离，L t

L 代表前车的车长，反应时间t r F 可从历史数据中获得。最小追踪间隔可以通过相邻列车

最小追踪时间间隔近似求得：

t F (s )?t L (s )≥t r F +t b F (s )+t safe

F (s ) 其中，t F (s )和t L (s )分别代表前后车经过位置s 的时刻，后车的制动时间t b F (s )以及由于

安全距离和前车车长造成的时间间隔t safe

F (s )可以如下计算： t b F (s )=v F (s )/a b F

t safe F (s )=(S SM +L t L )/v F (s )

其中，v F (s )代表后车在位置s 的速度。

为了使列车操纵不受信号系统的影响，也就是不被前车的运行影响（原则上缓行或者停车），最小追踪时间间隔这个概念被引入。车站发车间隔如下定义：

H min,MBS =t d L +t in?out =t d L +t r F +t b,max F +t safe L

同时，t in?out =t r F +t b,max F +t safe L ，其中t b,max F 表示后车从最高速度制动到停止所需

的时间，也就是说t b,max F =v max F /a b F ，t safe L 表示前车完全驶出安全部分（特殊的部分，

用来保护前车），其中包含安全距离，t safe L =√2(S SM +L t L +L s )/a acc L ，前车的加速度

a acc L 在计算式被视为常量。

4.FBS 的约束条件以及转化为MILP 问题

在固定闭塞信号系统中，由于前车造成的约束条件首先被离散化到每一个网格点s k ，k =1,2,…,N +1，然后这些逻辑约束被转化线性约束，这些线性约束可以被MILP 方法利用两种求解方法被利用到固定闭塞信号系统下多车运行操纵优化问题的求解。

4.1离散化FBS 约束条件

我们假设对于每个m ∈{1,2,…,M }，存在一个指数l m ∈{1,2,…,N +1}，下式成立： S FB,m =s lm 并且定义一个分段常数等式，如下：

?(k )=m ，对于l m

由于前车运行造成的限制条件可以转化到每个网格点s k ，如下所诉：

如果前车与后车位于同一个闭塞区间，也就是说t F (k )∈(t L (l ?(k )),t L (l ?(k )+1)]，这种情况在现实中是不允许发生的，那么后车的速度必须降低到最低速度，也就是说： v F (k )=V min ②如果前车在后车前一个闭塞区间，也就是说t F (k )∈(t L (l ?(k )+1),t L (l ?(k )+2)]，那么后车在位置s l ?(k )的速度需小于或等于V yellow ，在位置s l ?(k )+1的速度需等于V min ，也就是说：

v l ?(k )

F ≤V yellow v l ?(k )+1

F =V min 在间隔[s l ?(k ),s l ?(k )+1]内，假设列车减速度为常量，根据位置，速度以及加速度之间的

关系，我们可以得到：

2a l ?(k )(s l ?(k )+1?s l ?(k ))=V min 2?V yellow 2

2a l ?(k )(s k ?s l ?(k ))=V ?yellow,k 2?V yellow 2

其中，a l ?(k )表示减加速度，V

?yellow,k 表示在位置s k 的最高速度，通过消除式子（23）与式子（24）中的a l ?(k )，得到：

V ?yellow,k =√V yellow 2+(V min 2?V yellow 2)s k ?s l ?(k )s l ?(k )+1?s l ?(k )

因此，V

?yellow,k 为常量，等式右边均为常量，可得到如下约束： V F (k )≤V ?yellow,k 如果前车在后车前两个闭塞区间，也就是说t F (k )∈(t L (l ?(k )+2),t L (l ?(k )+3)]，那么后车在位置s l ?(k )的速度需小于或等于V m ax ，

在位置s l ?(k )+1的速度需小于或等于V yellow ，也就是说：

v l ?(k )

F ≤V m ax v l ?(k )+1

F ≤V yellow 与求解V

?yellow,k 相似，我们可以得到： V ?max,k =√V max 2+(V yellow 2?V max 2)s k ?s l ?(k )s l ?(k )+1?s l ?(k )

其中，V

?max,k 表示列车在位置s k 的最高速度，V ?max,k 同样也是一个常量，因此，我们可以得到如下约束：

V F (k )≤V

?max,k 4.2将FBS 中的约束条件引入MILP 求法

为了将上诉的逻辑约束转化为线性约束，引入以下二进制变量：

[t F (k)≤t L (l ?(k )+1)]?[δ1(k)=1]，

[t F (k)≤t L (l ?(k )+2)]?[δ2(k)=1]，

[t F (k)≤t L (l ?(k )+3)]?[δ3(k)=1]??

以及必要的额外约束：

t F (k )>t L (k ) 上式意味着后车经过位置s k 晚于前车经过该位置，以避免列车冲突。同时，基于变量δ1(k)，δ2(k)，δ3(k)的定义，可以得到下面的逻辑式：

[δ1(k )=1]?[δ2(k )=δ3(k )=1]，

[δ2(k )=1]?[δ3(k )=1]

三显示固定闭塞信号系统环境下，由前车造成的约束可以如下重建：

δ1(k)v F (k)≤V min

(1?δ1(k)δ2(k)v F (k))≤V

?yellow,k (1?δ2(k)δ3(k)v F (k))≤V

?yellow,k 式

（33）中小于等于情况，当δ1(k )=0，δ1(k)v F (k)等于0，并不等于V min 。 通过定义

M ?i =T max F ?t L (l ?(k )+i )≥max?(t F (k )?t F (l ?(k )+i ))，

m ?i =min?(t F (k )?t F (l ?(k )+i ))≥T min F ?t L (l ?(k )+i )，

通过转化特性式子（5）逻辑约束式（30）等价于下式：

t F (k )?t L (l ?(k )+i )≤M

?i (1?δi (k)), t F (k )?t L (l ?(k )+i )≥ε+(m ?i ?ε)?δi (k)

对于i =1,2,3，其中，T max F 表示前车到达终点站的时刻，T min F 表示后侧滑离开始发站

的时刻，ε表示精度。同时，我们定义二进制变量δ4(k )=δ1(k )δ2(k )和δ5(k )=δ2(k )δ3(k )去处理非线性条件δ1(k )δ2(k )和δ2(k )δ3(k )。根据第2.2节的转化特性式子，δ4(k )和δ5(k )的定义与式子（6）的等价。同时，辅助变量z i F (k)被引入去解决非线性因素δi (k)v F (k),被定义如下：

z i F (k )=δi (k )v F (k )??????i =1?5

这个定义与式（7）等价，因此固定闭塞信号系统下由于前车造成的约束条件可以转化为混合整数规划的条件。

4.3基于FBS 系统下的多车优化控制问题

现在有两种发放可以用于求解基于FBS 系统下的多车优化控制问题，为了方便起见，我们将问题考虑为两车的运行曲线优化控制问题，但是方法论可以被应用到多车问题中。其中一种为婪式算法（时序法），前车的运行速度曲线首先被确定，然后后车的速度曲线优化基于前车的结果。另外一种方法是同时计算前后车的运行速度曲线。 使用婪式算法，前车的运行速度曲线转化为MILP 问题进行求解，然后后车的控制优化求解与Wang （2011）中的求解相同，只是加入FBS 系统下前车产生的约束条件。式（8）和式（9）中的协同矩阵由前车确定，前车的运行速度曲线确定之后，

t L (l ?(k )+1)变为已知值，因此M ?i 和m ?i 变为常量，那么式（36）变为线性约束。

当前车与后车的优化同时进行的时候，两车的模型可由形式（8）和（9）确定，优化控制问题可以被写成形式（11）和（12），但是其中包含由于列车以及FBS 系统产生的约束条件。但是当前车的运行速度曲线被优化后，t L (l ?(k )+1)也变为已知值，约束条件（36）可以被写成：

(T max F ?t L (l ?(k )+i ))δi (k )≤?t F (k )+T max F ，

(T min F ?t L (l ?(k )+i )?ε)δi (k )≤t F (k )?ε?t L (l ?(k )+i )

为了解决非线性问题（38），我们定义：

z i L (k )=t L (l ?(k )+i )δi (k )??????i =1?3 与z i F (k )相似，辅助变量z i L (k )等同于线性约束，依据第2.2节阐述的转化规则。

5. MBS 的约束条件以及转化为MILP 问题

5.1离散化MBS 约束条件

回忆第2.3节中，将列车操纵优化的混合逻辑动态模型离散化成N 个空间，伴随着网格点s k ，k =1,…,N +1。这里我们将MBS 系统的约束在网格点s k 离散化如下：

t F (k )≥t L (k )+t r F +t b F (k )+t safe

F (k )??????k =1,2,…,N t F (k )≥t L (k )+t r F +t b,max F +t d L +t safe L ??????k =N +1

同时，一些调解约束被引入，为了确保网格点减的位置能够满足由于MBS 系统造成的约束。根据式子（15），得到对于s ∈[s k ,s k+1]满足以下约束：

t F (s )?t r F ?t b F (s )?t safe F (s)≥t L (s)

如果我们假设式子（42）左边为仿射方程，在区间[s k ,s k+1]内，我们可以得到：

(1?α)(t F (k )?t r F ?t b F (k )?t safe

F (k ))+α(t F (1+k )) ?t r F ?t b F (k +1)?t safe

F (k +1)≥t L (s +α?s k ) 其中，α在有限子集S α∈[0,1)，举例来说S α={0.1,0.2,…,0.9}，如果前车是固定的，那么t L (s +α?s k )为已知量，对于α=0和α=1是检索的（除了如果k =N ?1）。但是如果前车的运行速度曲线是未知的，那么将同时优化前后车的运行速度曲线，在这种情况下，条件t L (s +α?s k )是未知的。如果我们假设式子（42）右边为仿射方程，也就是：

t L (s +α?s k )=(1?α)t L (k )+αt L (k +1)，

那么需要检验式子（42）在点k 和k +1（也就是α=0和α=1）是否成立，如果成立，那么线性式子（42）将在任何中间点成立。

5.2将MBS 中的约束条件引入MILP 求法

声明约束（40），（41）和（43）在t L (k)，t L (k +1)和t F (k)是线性的，但是因为安全时间间隔（17）对于列车速度v F (k)是非线性，所有约束在v F (k)和v F (k +1)是非线性的。更多的，单位质量的动能是E F (k)而非E F (k )=0.5(v F (k))2（具体参考第2章），

因此对于E F (k)，制动时间t b F (k)和安全时间间隔t safe F (k)是非线性方程，如下：

t b F (k )=1

a b F √2E F (k ) t safe F (k )=(S SM +L t L F

非线性方程f 1(?):E F →√2E F 和f 2(?):E F →1/√2E F 可以通过PWA 方程进行近似化处理，存在很多种近似化处理方法。在本文中，我们首先选择采用基于复杂性和精确性的多层次PWA 方法，然后用最小二乘法对f 1(?)对应的仿射方程的间隔长度和特征量进行

优化，然后对f 2(?)进行相同的处理。如果我们考虑二段式，那么方程f 1(?)和f 2(?)可以被PWA 近似处理如下：

f 1,PWA (E F (

k ))={α1E F (k )+β1?????E min ≤E F (k )≤E 1α2E F (k )+β2?????E 1≤E F (k )≤E max f 2,PWA (E F (

k ))={λ1E F (k )+μ1?????E min ≤E F (k )≤E 1λ2E F (k )+μ2?????E 1≤E F (k )≤E max 其中，α1，α2，β1，β2，λ1，λ2，μ1，μ2和E 1优化参数。现在约束（40）可以被近似化如下线性方程：

t F (k )≥t L (k )+t r F +1b F (α1E F (k )+β1)+(S SM +L t L )(λ1E F (k )+μ1)?? E min ≤E F (k )≤E 1 t F (

k )≥t L (k )+t r F +1b F (α2E F (k )+β2)+(S SM +L t L )(λ2E F (k )+μ2)?? E 1≤E F (k )≤E max 同样的，约束条件（41）和（42）可以同样被转化为线性方程，这些由于MBS 系统造成的近似约束可以转入到MILP 问题。

5.3基于MBS 系统下的多车优化控制问题

第4.3节中提到的婪式算法和同步算法可以被应用到MBS 系统中。运用婪式算法，当前车的运行速度曲线已知时，变量t L (k)和t L (s +α?s k )变成已知量对于后车，后车的运行速度曲线的求解和Wang （2011）的文章中阐述的相同，但是需要考虑第

5.2节中提到的关于MBS 系统形成的约束条件。当同时优化前后车的运行速度曲线时，采用形式（8）和形式（9），然后写成MILP 问题的格式（11）和（12）。与单车问题相比较，多车问题的参数，二进制变量，辅助变量以及约束条件呈现线性增加，因此问题的规模变得更大，需要更多的求解时间。因为同步算法为全局优化算法，其结果将比婪式算法更加优秀。

6案例仿真

为了体现婪式算法和同步算法在解决基于固定闭塞信号系统和移动闭塞信号系统的多车运行速度曲线优化的表现，以北京亦庄地铁线路为基础进行仿真。

6.1建立

我们将MILP 方法与伪谱法进行比较，同时考虑婪式算法和同步算法，为了简单起见，我们仅仅考虑北京亦庄铁路线路的两个站：宋家庄站和肖村站，两个站之间的线路长度伪1332m ，速度限制和线路坡度如下图三所示：

图三宋家庄站与肖村站间速度限制和线路坡度

列车与线路参数如下表一所示：

表一列车与线路参数

在文献（Hansen，2008）中，旋转质量参数取1.06，本文中也取相同的值。根据第2.1节中的假设单位质量动能稍大于正临界值E min。在本次仿真中，我们取最小单位质量的动能为0.1J，对于列车的速度，最大牵引力为非线性方程，方程的最大值我们取315kN。列车的最大加速度和最大减加速度，取0.8m/s2，为了确保低驾驶水平的实际可以按照速度曲线进行操作。当取列车质量为2.78×105kg时，最大牵引力和最大制动力分别为222.4kN和?222.4kN。本文中考虑的目标函数为，不存在再生制动情况下的列车能量消耗。

在本次仿真中，两趟列车从宋家庄站向肖村站发车，发车间隔为75s。我们考虑两种情况：FBS系统和MBS系统，假设前车在整个仿真过程中有故障，其最高速度为40km/h，也就是11.1m/s，并且假设前后车在肖村站进入不同的站台股道，后车将在该站超车。计算最小追踪间隔的参数如表二所示：

表二计算最小追踪间隔的参数

列车的长度取90m，司机的反应时间取1s。对于FBS系统，我们假设两站之间包含

四个等长的闭塞区间，也就是长度为333m，根据式子（13）计算最小追踪间隔为98.4s。对于MBS系统，根据表二中的数据，求得t in?out=44.6s，最小追踪间隔为69.6s。指出初始的最小追踪间隔为75s，大于MBS系统的最小追踪间隔，小于FBS系统的最小追踪间隔。

在MILP问题中，区间[s k,s k+1]对应的间隔长度?s k的取决于限制速度，线路坡度，固定闭塞区间长度等因素，如果分割的空间数量N越大，那么MILP算法的精度越高，但是需要更长的计算时间。对于这种情况，我们取N为20，40和60，那么间隔长度对应分别为66.6m，33.3m和22.2m。运行速度曲线优化问题中的非线性条件如列车模型中的微分方程，通过PWA方程进行近似处理，如果处理PWA方程时，取更多的子功能，那么近似处理的精度将更高，同时对两段化和三段化处理进行比。我们在Matlab 环境下，用CPLEX对问题进行求解。

伪谱法是解决优化控制问题表现最好的算法，利用基于正交配置点，如高斯点，正交多项式对状态方程和控制方程进行近似处理，微分方程可以用代数方程组进行近似处理。伪谱法近似处理的误差可以用过增加正交配置点的数量来降低，我们分别取高斯点数量为20，40，80和120。现存在很多工具箱可以求解伪谱法，其中一个是PROPT，支持代数微分方程，具有多个求解工具，如MINOS和SNOPT，来求解非线性问题。本次仿真中，我们在Matlab和SNOPT环境下，通过Tomlab进行求解。6.2基于FBS系统的仿真结果

表三展示了基于FBS系统的MILP和伪谱法用婪式算法和同步算法的结果表现对比，结果表现中包含的例如能量消耗和时间差异，均通过非线性方程进行求解。总能量消耗包含了前后车的能量消耗，时间差异表示前后车实际运行时间与计划运行时间的差值。每趟车的能量消耗与实际运行时间和计划运行时间相关，如果实际运行时间比计划运行时间长，那么列车的平均运行速度比较低，那么总能量消耗就比较少。CPU计算时间表示问题求解的总用时，对于婪式算法，是前后车分别计算时间，对于同步算法，是前后车同步计算完成所需的时间。

表三基于FBS系统的MILP和伪谱法用婪式算法和同步算法的结果表现对比

更多的，图4，5展示了，基于FBS系统，两种不同方法求得的优化控制输入量和速度位置曲线，其中，MILP方对应的空间数量N和伪谱法对应的正交配置点都为40，并且MILP方法对应的PWA法采用三段式进行近似处理。速度位置曲线由优化控制输入量求得，优化控制输入量通过求解非线性列车模型得到。可以从图中看出，在FBS系统中，后车的行驶受到前车的影响，如在666m位置处信号机对后车现实黄灯，主要是由于求解中采用的最小追踪间隔伪75s，小于基于FBS系统计算所得的最小追踪间隔98.4s，因此后车必须降速，来满足黄灯对应的速度限制，40km/h，也就是11.1m/s。

图四基于FBS系统的75s最小追踪间隔下采用婪式算法求解MILP和伪谱法所得的运行速度曲线

图五基于FBS系统的75s最小追踪间隔下采用同步算法求解MILP和伪谱法所得的运行速度曲线在图四中，利用婪式法求解MILP问题对应的能量消耗，时间误差，计算时间分别为112.59MJ，2.21s，96.56s，利用婪式法求解伪谱法问题对应的能量消耗，时间误差，计算时间分别为110.04MJ，1.43s，1381.71s，显而易见建立MILP模型求解所得能量消耗和时间误差比建立伪谱法模型求解所得的值要大，但是计算时间短的多。

在图五中，利用同步法求解MILP问题对应的能量消耗，时间误差，计算时间分别为106.58MJ，1.73s，184.76s，利用同步法求解伪谱法问题对应的能量消耗，时间误差，计算时间分别为106.44MJ，1.68s，1521.30s，当与婪式法求解MILP问题对应的能量消耗，时间误差相比，变得更小，主要是由于前后车同时进行优化，并且前车的优化过程中考虑后车的运行，但是计算时间变得更长，是因为优化问题的规模变成增长了一倍。当与同步法求解伪谱法问题作比较，呈现相同的差异。

可以从表三看出，在问题离散化程度相同的情况下，利用婪式法求解MILP问题对应的能量消耗，时间误差比利用婪式法求解伪谱法问题对应的要大，但是伪谱法所需的求解时间是MILP法求解时间的一到两个数量级。当增加MILP方对应的空间数量N和伪谱法对应的正交配置点时，求得的能量消耗，时间误差更小，但是两种方法之间仍然呈现相同的规律。还可以看出利用同步法求解MILP问题对应的能量消耗比用婪式法求解MILP问题对应的能量消耗要小，但是求解时间更长，对于伪谱法的求解结果呈现相同的规律。

6.3基于MBS系统的仿真结果

对于MBS系统，婪式算法和同步算法的计算结果如表四所示，图六和图七中展示了速度位置曲线和优化控制输入量，基于与图四和图五中相同的参数设置，同样的，速度位置曲线以及结束时间误差等参数是基于优化控制输入量求得的，从图中可以看出后车并没有受到前车的干扰，主要是由于设置的最小追踪间隔为75s大于基于MBS 求得的最小追踪间隔。

表四基于MBS系统的MILP和伪谱法用婪式算法和同步算法的结果表现对比

图六基于MBS系统75s最小追踪间隔下采用婪式算法求解MILP和伪谱法所得的运行速度曲线

图七基于MBS系统75s最小追踪间隔下采用同步算法求解MILP和伪谱法所得的运行速度曲线

在图六中，利用婪式法求解MILP问题对应的能量消耗，时间误差，计算时间分别为67.89MJ，1.68s，75.98s，利用婪式法求解伪谱法问题对应的能量消耗，时间误差，计算时间分别为67.42MJ，1.66s，729.24s，利用婪式法求解MILP问题对应的能量消耗，时间误差比利用婪式法求解伪谱法问题对应的稍大，但是利用婪式法求解伪谱法问题对应的计算时间是利用婪式法求解伪谱法问题对应的值大一个数量级。同样的，利用同步算法进行求解，两种方法呈现相同的规律。当增加MILP方对应的空间数量N和伪谱法对应的正交配置点时，求得的能量消耗，时间误差更小，但是计算时间迅速增加。对于婪式法或者同步算法求解MILP问题，当PWA近似处理式，使用多端式，结束时间误差也会变小。伪谱法相较于MILP方法，在优化控制方面表现的更加优秀，但是由于优化问题规模的增加，其计算负担更大。

6.4小结

上述仿真结果显示，与婪式伪谱法相比，婪式MILP法得到的能量消耗和时间误差更大，但是计算时间却要小一到两个数量级，同样地，同步法下两种方法呈现相同的差异性。更多的，婪式MILP法求得的能量消耗要比同步MILP法求得的要大，但是计算时间更短，对于婪式伪谱法和同步伪谱法呈现相同的差异性。

7结论和展望

在本文中，我们提出了婪式法和同步法对多车运行优化控制进行求解。对于婪式算法，前车的运行速度曲线先被求得，基于前车的数据，对后车的运行速度曲线进行求解。对于同步算法，对前后车的运行速度曲线同时进行求解。在固定闭塞信号系统和移动闭塞信号系统中，由于前车造成的约束条件被引入到运行速度曲线优化问题中，列车模型和约束条件中的非线性项目通过分段仿射公式进行近似处理为线性项目，然后可以将问题重新建立为混合整数规划模型。还对MILP法与伪谱法进行两种不同求解实际问题结果表现的比较，结果显示MILP法与伪谱法有差不多的表现，并且需要的计算时间更短。

将来需要做的工作是MILP法与伪谱法进行更加详细及广泛的比较（利用其他非线性求解方法，如MINOS和KNITRO），以及动态算法，其他文献中提出的方法和框架（例如AMPL，APMonitor和ASCEBND），并且进行其他实际实例和情况的仿真以及对比。