考虑动态随机客流的城市轨道交通列车时刻表仿真优化

Simulation Optimization of Urban Rail Transit Train Timetable Based on Dynamic Random Passenger Flow

摘要

为更加精准刻画动态随机客流与列车时刻表编制间的复杂匹配过程，综合考虑乘客动态随机到达过程、列车容量限制、车内拥挤度对乘客上下车速率影响、停站时间动态调整和列车运行安全防护等重要运营特征，构建城市轨道交通列车运行离散仿真模型。在此基础上，建立以最小化乘客平均候车时间为目标的随机非线性优化模型，并设计一种基于仿真的有限差分随机逼近算法，用于优化城市轨道交通发车方案。以国内某条线路的实际运营数据为例验证仿真模型及优化方法，结果表明所构建模型及算法具有较好的优化效果和计算效率，优化后的列车实绩运行时刻能更好地适应客流需求的动态随机性，且在不增加发车频次的前提下有效降低乘客平均候车时间，提升城市轨道交通运营管理水平。

Abstract

To more accurately depict the complex matching process between dynamic random passenger flow and train capacity, this paper comprehensively considered important operational characteristics such as passenger dynamic random arrival process, train capacity constraints, carriage congestion effects on passenger boarding and alighting rates, dynamic station dwell time adjustments, and train operation safety protection, and constructed a discrete simulation model for the operation of urban rail transit trains. Based on this, a random nonlinear optimization model was built to minimize average passenger waiting time, and a simulation-based finite difference stochastic approximation algorithm was designed to optimize the departure plan of trains. Taking the actual operational data of an urban rail transit line in China as an example to verify the simulation model and optimization methods, the results show that the constructed model and algorithm have effective optimization and computational efficiency. The optimized actual train operation time can better adapt to the dynamic randomness of passenger demand, significantly reducing the average passenger waiting time without increasing the frequency of train services, which is of great value for further improving the quality of urban rail transit operation and management.

Graphical abstract

关键词

城市轨道交通 / 动态随机客流 / 列车时刻表 / 系统仿真 / 有限差分算法

Key words

Urban Rail Transit / Dynamic Random Passenger Flow / Train Timetable / System Simulation / Finite Difference Algorithm

0　引言

截至2022年底，中国内地已开通308条城市轨道交通线路，涵盖了55个城市，运营里程达10 287.45 km^[1]。凭借运量大、速度快、安全、便捷、高效和低耗能等优点，城市轨道交通分担公共交通客流比例逐年提升，已成为城市综合交通运输网络的骨干。随着城市轨道交通规模的不断扩大，其服务水平和运营效率也面临着一系列新的挑战。2020年3月发布的《中国城市轨道交通智慧城轨发展纲要》提出城市轨道交通运营公司应更加精准地匹配运力和客流需求，以提高社会效益和经济效益，而编制以乘客出行需求为导向的列车时刻表则是精准匹配供需关系的关键。

目前，国内大多数城市轨道交通运营公司根据需求量大小将运营时段划分为早高峰、平峰、晚高峰等多个时段，各时段内一般采用等间隔发车模式；此外，也有一些运营公司会以小时为单位划分运营时段，并根据最大断面客流量，确定发车频率。这种时刻表编制方法简单易行，但由于其时段划分粒度相对较粗，每个时段内的运力配置方案不够精细化，无法适应时段内出行需求的动态性，导致时段内运力紧张和虚靡的现象并存^[2]，较大程度影响了服务水平和运输资源利用效率。为此，运营企业和众多学者开始聚焦动态客流下的城市轨道交通列车时刻表优化问题。史峰等^[3]面向时变客流需求，以最小化列车对数和车底使用数量为目标，构建城市轨道交通列车运行计划优化模型，并设计双向关联序列化算法和匈牙利算法对其求解。Niu等^[4]在给定时变客流OD需求矩阵的前提下，构建了以最小化乘客等待时间为目标的非线性整数规划模型，并利用GAMS求解器求解。Yin等^[5]基于动态客流需求，提出了一种以最小化运行成本(即能源消耗)和乘客等待时间为目标的列车时刻表优化模型，并开发一种基于拉格朗日松弛(LR)的启发式算法求解。然而，上述文献仅考虑了城市轨道交通客流的动态性，却忽略了客流的随机性特征。近年来，有学者针对每日变化波动的随机客流需求进行研究^[6]。Hassannayebi等^[7]针对不确定性的客流需求，提出了一种基于两阶段遗传算法的仿真优化方法，以最小化乘客等待时间，提高服务质量；张春田等^[8]以极小化铁路运营公司总运营成本为目标，利用轻鲁棒技术，构建不确定旅客需求下高速铁路鲁棒列车开行方案优化模型。Han等^[9]针对旅客需求的随机性，构建基于场景的列车时刻表框架，将旅客需求可能性根据乘客需求画像和数量，划分为一系列场景，设计高级自适应算法(AANSGA-II)求解。

尽管上述针对城市轨道交通时刻表的研究已取得不错进展，但这些文献主要采用数学建模方法，难以准确地描述乘客行为与列车运行过程间的复杂交互关系。例如，为降低数学模型的构建难度和求解复杂度，既有研究通常将乘客候车时间描述为列车到站间隔的1/2^[10]，这种简单的处理方法难以精准计算每名乘客的候车时间，更是无法描述大客流下部分乘客无法上车而导致其候车时间进一步增加的现象。相较数学建模方法而言，仿真更能够精细化描述每个乘客的出行行为和各个列车的运行状态^[11-12]。然而，既有文献中尚未有研究能够综合考虑客流需求的动态随机性、车内拥挤度对上下车乘客速率影响、停站时间随需求动态调整、列车容量限制、列车运行安全条件等因素建立仿真模型。综上，应针对城市轨道交通时刻表优化问题，综合考虑现实中城市轨道交通运营的动态性特征，建立城市轨道交通列车运行仿真模型，并进一步设计列车时刻表随机非线性优化模型及有限差分随机逼近算法，从而提高列车运力与动态随机客流的匹配程度，降低乘客平均候车时间。

1　仿真场景描述

1.1　线路场景

假设某条城市轨道交通线路共设有K个车站，车站1、车站K分别表示线路始发站和终点站。由于城市轨道交通线路通常包括上、下行2个方向，令f=\mathrm{1}表示下行方向，f=\mathrm{2}表示上行方向。由于2个方向相反，可以认为上行和下行之间的客流需求相互独立。为直观描述2个方向的运营情况，将K个车站抽象化为\mathrm{2}K个站台，其中站台1和站台\mathrm{2}K、站台2和站台2K-1等，在实际车站示意图中对应同一个车站，某条城市轨道交通线路运行示意图如图1所示。列车从车辆段发出，依次经过站台1，2，…，K-1后到达该方向的终点站台K，进行折返操作后，继续沿上行方向运行，依次经过站台K+1，K+2，…，2K-1后到达上行方向的终点站台\mathrm{2}K，至此该列车结束一次运行任务。

1.2　列车运行场景

列车在线路上的运行速度与列车自动驾驶系统(Automatic Train Operation，ATO)所采取的控车曲线有关。常用的控车曲线通常包括4个阶段：加速、匀速、惰行和减速。这些阶段的具体控制方式和曲线形状会因运营企业、线路条件、列车型号等因素而有所差异。为减少模型复杂性并提高计算效率，将列车在区间的运行过程简化为加速、匀速、制动3个过程。在加速阶段，列车会根据区间运行距离和限速要求，计算并设计出加速曲线，并在保证列车安全运行前提下，逐渐加速到目标速度。在匀速阶段，列车会根据当前位置和目标速度，维持以恒定的速度行驶。当列车与相邻列车间不满足安全间隔或即将到达站台时，应根据所处位置及自身制动能力进行减速。

1.3　客流到达场景描述

为描述乘客到达的动态随机过程，令\omega =\{p_i=(o_i，d_i，t_i)|i=\mathrm{1}，\mathrm{2}，\dots ，n\}表示随机客流场景样本，其中，p_i=(o_i，d_i，t_i)表示乘客i的出行信息，n表示该场景样本中最大乘客数，o_i和d_i分别为乘客i的出发车站、到达车站，t_i为该乘客到达出发车站站台的时刻。样本\omega可以通过一定的随机过程采样获得：首先，将研究时段划分为多个细小的子时段，分别统计各个子时段的客流强度；然后，根据泊松到达过程生成各个子时段的样本序列；最后，将所有时段的样本序列合并，构成整个研究时段的随机客流样本。重复多次该采样过程，可以获得客流随机空间下的多个客流场景，将其存放至场景集合W=\{{\omega }_{\mathrm{1}},{\omega }_{\mathrm{2}},\dots ，{\omega }_z\}，其中z表示最大场景个数，用于后续仿真及优化研究。

2　城市轨道交通运行仿真模型

2.1　模型假设

为清晰简洁地刻画乘客出行和列车运行过程，对仿真模型作如下合理性假设。

（1）所有乘客进入站台后，均按顺序排队候车，且候车乘客在各屏蔽门前呈均匀分布。

（2）候车乘客在当前列车容量充足的情况下均选择上车，并遵循“先下后上”和“先到先上”原则。

（3）每个站台每次只能容纳1列列车，所有列车均采用站站停的开行模式，且运行中不发生越行和跳停的情况。

2.2　仿真模型构成主体

由于城市轨道交通列车运行系统的复杂性和关联性，将其分解为3个相互独立又相互关联的部分，即乘客主体、列车主体以及线路主体，3个主体间的交互关系如图2所示。乘客主体是系统中的核心，其行为包括乘客进站、到达站台、排队候车、上下车、出站等；列车主体是乘客完成出行的载体，其行为包括发车、区间运行、停站上下车作业、入车辆段；线路主体包括车站、站台等元素，构成乘客主体和列车主体的仿真环境。

2.3　仿真更新规则

在仿真系统中，对研究时段进行离散化处理，将其划分为若干个固定长度的等时间间隔，即步长\mathrm{\Delta }t(为提高仿真的精确度，设\mathrm{\Delta }t=\mathrm{1}\mathrm{s})。仿真从研究时段内初始时刻开始，在每个时间步长内(即t\to t+\mathrm{\Delta }t时间段)，依据更新规则完成对乘客和列车行为的动态模拟。

2.3.1　乘客出行行为规则设计

为准确把握乘客的出行过程，随着仿真时间的推移，应对乘客所属行为状态进行更新。根据AFC数据已知乘客i刷卡通过闸机时的时间为t{}_i^{\mathrm{i}\mathrm{n}}，当系统时间t=t_i^{\mathrm{i}\mathrm{n}}时乘客进入城市轨道交通系统。乘客i在车站内走行t_k^{\mathrm{w}\mathrm{a}\mathrm{l}\mathrm{k}}时间后到达其上车站台k，即当系统时间t=t_i^{\mathrm{i}\mathrm{n}}+t_k^{\mathrm{w}\mathrm{a}\mathrm{l}\mathrm{k}}时，该名乘客到达站台，并按顺序进入候车队列中。乘客能否上车需进行3个条件判断(即乘客到达站台时间是否早于列车到达时间；是否满足列车容量限制；是否满足线路规定的停站时间限制)，若同时满足则该名乘客可以上车，否则应继续候车。乘客上车后随车运行，当列车运行到乘客目的站台时，乘客下车并完成本次出行。

2.3.2　列车区间运行规则设计

为保证行车安全，列车在更新时需充分考虑到前方相邻列车的运行情况。在时刻t对列车l进行更新时，已知此时列车运行速度为v_{l,t}，加速度为a_l，减速度为d_l，所在区间m长度为s_m，区间m最大运行速度v_m^{\mathrm{m}\mathrm{a}\mathrm{x}}。列车区间运行示意图如图3所示。

图3中，令s_l^{\mathrm{t}}表示列车l车头与前方相邻列车车尾的距离，s_l^{\mathrm{p}}表示列车l车头与前方站台停车点的距离。计算列车可行驶的安全距离{\widehat{s}}_{l,t}=\mathrm{m}\mathrm{i}\mathrm{n}\{s_l^{\mathrm{t}}，s_l^{\mathrm{p}}\}和列车当前时刻最大容许速度{\widehat{v}}_{l,t}=\sqrt[]{\mathrm{2}{d}_l{\widehat{s}}_{l,t}}。判断列车当前速度v_{l,t}和最大容许速度{\widehat{v}}_{l,t}的关系，需分情况讨论。

情况1：如果v_{l,t}\ge {\widehat{v}}_{l,t}，则列车应制动，更新列车l下一时刻的运行速度v_{l,t+\mathrm{\Delta }t}=\mathrm{m}\mathrm{a}\mathrm{x}\{\mathrm{0}，{\widehat{v}}_{l,t}-d_l\mathrm{\Delta }t\}和列车位置x_{l,t+\mathrm{\Delta }t}=x_{l,t}+({{\widehat{v}}^{\mathrm{2}}}_{l,t}-v_{l,t+\mathrm{\Delta }t}^{\mathrm{2}})/(\mathrm{2}{d}_l)。

情况2：如果，则列车应加速，计算列车正常加速后在下一时刻的速度v\mathrm{\text{'}}=\mathrm{m}\mathrm{i}\mathrm{n}\{v_{l,t}+a_l\mathrm{\Delta }t，v_m^{\mathrm{m}\mathrm{a}\mathrm{x}}\}，位置x\mathrm{\text{'}}=x_{l,t}+(v_{l,t}+v\mathrm{\text{'}})\mathrm{\Delta }t/\mathrm{2}和安全距离\widehat{s}\mathrm{\text{'}}={\widehat{s}}_{l,t}-(v_{l,t}+v\mathrm{\text{'}})\mathrm{\Delta }t/\mathrm{2}，并进一步判断v\mathrm{\text{'}}是否满足安全要求。

当v\mathrm{\text{'}}\le \sqrt[]{\mathrm{2}{d}_l\widehat{s}\mathrm{\text{'}}}，表示按预期速度加速后的下一时刻速度不会超过容许速度，则更新下一时刻的速度v_{l,t+\mathrm{\Delta }t}=v\mathrm{\text{'}}和位置x_{l,t+\mathrm{\Delta }t}=x\mathrm{\text{'}}。

当v\mathrm{\text{'}}>\sqrt[]{\mathrm{2}{d}_l\widehat{s}\mathrm{\text{'}}}，表示按预期速度加速后会导致下一时刻速度超出容许速度，因此需要通过求解以下方程计算下一时刻最大容许速度{\widehat{v}}_{l,t+\mathrm{\Delta }t}。

\frac{v_{l,t}+{\widehat{v}}_{l,t+\mathrm{\Delta }t}}{\mathrm{2}}\mathrm{\Delta }t+\frac{{{\widehat{v}}^{\mathrm{2}}}_{l,t+\mathrm{\Delta }t}}{\mathrm{2}{d}_l}={\widehat{s}}_{l,t}⑴

公式⑴可进一步转化为公式⑵。

{{\widehat{v}}^{\mathrm{2}}}_{l,t+\mathrm{\Delta }t}+d_l\mathrm{\Delta }t\cdot {\widehat{v}}_{l,t+\mathrm{\Delta }}+d_l\mathrm{\Delta }t\cdot v_{l,t}-\mathrm{2}{d}_l{\widehat{s}}_{l,t}=\mathrm{0}⑵

公式⑵可通过一元二次方程求根公式计算{\widehat{v}}_{l,t+\mathrm{\Delta }t}

{\widehat{v}}_{l,t+\mathrm{\Delta }t}=\frac{-d_l\mathrm{\Delta }t+\sqrt[]{(d_l{\mathrm{\Delta }t)}^{\mathrm{2}}-\mathrm{4}(d_l\mathrm{\Delta }t\cdot v_{l,t}-\mathrm{2}{d}_l{\widehat{s}}_{l,t})}}{\mathrm{2}}⑶

更新下一时刻列车速度v_{l,t+\mathrm{\Delta }t}={\widehat{v}}_{l,t+\mathrm{\Delta }t}和位置x_{l,t+\mathrm{\Delta }t}=x_{l,t}+(v_{l,t}+v_{l,t+\mathrm{\Delta }t})\mathrm{\Delta }t/\mathrm{2}。

2.3.3　列车动态停站时间确定及乘客上下车规则设计

列车的停站时间会受列车开关门时间、乘客上下车需求人数、车内拥挤度、列车最大容量限制、线路最大最小停站时间设置等多重因素影响，为此设计一种考虑列车容量限制和车内拥挤度的动态停站时间确定规则。

实际运营中乘客上下车速率和列车满载率密切相关，文献^[13]选取上海市某条地铁7个典型车站对乘客平均上下车时间进行计算分析，得到乘客平均上下车时间与满载率的关系如表1所示。表1中满载率为50%情况下的上下车速率为0.64 s/人，但随着车内乘客数增加，平均上下车时间大幅增加。因此，为有效模拟不同拥挤度下乘客上下车速率，根据表1中的数据选取50%，80%两个具有代表性的满载率作为分界点，设计乘客平均上下车时间{\overline{t}}^{\mathrm{a}\mathrm{b}}关于满载率\vartheta的分段函数，如公式⑷所示。

⑷

基于上述分段函数，设计列车动态停站时间计算流程如下。

步骤1：计算列车所需停站时间。

①当容量为c_l的列车l到达某站台时，已知车内人数为n_l，车内待下车人数为n_l^{\mathrm{a}\mathrm{l}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{t}}，站台候车人数为n^{\mathrm{w}\mathrm{a}\mathrm{i}\mathrm{t}}，列车l的单侧车门数为r_l，考虑列车容量限制，初步判定待上车人数n_l^{\mathrm{b}\mathrm{o}\mathrm{a}\mathrm{r}\mathrm{d}}=\mathrm{m}\mathrm{i}\mathrm{n}\{c_l-n_l+n_l^{\mathrm{a}\mathrm{l}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{t}}，n^{\mathrm{w}\mathrm{a}\mathrm{i}\mathrm{t}}\}。

②令A=n_l^{\mathrm{a}\mathrm{l}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{t}}表示还未下车的人数，设置累计下车时间T_A=\mathrm{0}。逐一对待下车乘客进行下车操作，直至所有乘客下车(即A=\mathrm{0})。每更新一名下车乘客时，依次计算列车满载率\theta =n_l/c_l，根据公式⑷计算当前乘客下车速率{\mu }^{\mathrm{a}}={\overline{t}}^{a,b}/r_l，更新待下车人数A=A-\mathrm{1}、车内人数n_l=n_l-\mathrm{1}、累计下车时间T_A=T_A+{\mu }^{\mathrm{a}}。

③令B=n_l^{\mathrm{b}\mathrm{o}\mathrm{a}\mathrm{r}\mathrm{d}}表示还未上车的人数，设置累计上车时间T_B=\mathrm{0}。逐一对待上车乘客进行上车操作，直至站台候车乘客全部上车(即B=\mathrm{0})或总停站时间到达最大停站时间限制(即T_{OC}+T_A+T_B\ge T_{\mathrm{m}\mathrm{a}\mathrm{x}}^{\mathrm{d}\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}}，其中T_{OC}表示列车开关门时间)。每更新一名上车乘客时，依次计算列车满载率\theta =n_l/c_l，根据公式⑷计算当前乘客上车速率{\mu }^{\mathrm{b}}={\overline{t}}^{a,b}/r_l，更新待上车人数B=B-\mathrm{1}、车内人数n_l=n_l+\mathrm{1}、累计上车时间T_B=T_B+{\mu }^{\mathrm{b}}。

步骤2：确定列车最终停站时间。列车最终停站时间不能小于最小停站时间T_{\mathrm{m}\mathrm{a}\mathrm{x}}^{\mathrm{d}\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}}，因此令T^{\mathrm{d}\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}}=\mathrm{m}\mathrm{a}\mathrm{x}\{T_{\mathrm{m}\mathrm{i}\mathrm{n}}^{\mathrm{d}\mathrm{w}\mathrm{e}\mathrm{l}\mathrm{l}}，T_{OC}+T_A+T_B\}。

2.4　总体仿真流程

基于上述核心更新规则，构建城市轨道交通列车运行仿真模型。仿真的输入信息包括列车发车方案、列车相关参数(如列车型号、加(减)速度、车门数、容量等)、线路相关参数(车站和站台情况、区间里程、区间限速、行车安全间隔时间设置等)、乘客相关参数(乘客OD需求、到达出发车站时间等)、运行相关参数(区间运行时分、列车运行计划等)，城市轨道交通系统仿真流程如图4所示。

3　基于仿真的列车时刻表优化模型及算法

所构建的仿真模型能够精细化描述任意客流场景下城市轨道交通列车运行动态过程，并准确统计每一名乘客的候车时间。在该仿真模型基础上，设计以平均候车时间最小为目标的列车时刻表优化模型及基于仿真的随机逼近算法。

3.1　模型构建

3.1.1　优化目标

城市轨道交通发车方案的制定需要从乘客角度出发，因此将目标函数设置为动态随机客流下乘客平均候车等待时间最小

\mathrm{m}\mathrm{i}\mathrm{n}E[T^{\mathrm{w}\mathrm{a}\mathrm{i}\mathrm{t}}(X，\omega )]⑸

式中：X=[x_l|l=\mathrm{1}，\mathrm{2}，...，L]为研究时段内线路上列车发车时刻向量；x_l为从始发车站发出第l个车次的时刻；T^{\mathrm{w}\mathrm{a}\mathrm{i}\mathrm{t}}(X，\omega )为在发车方案\mathrm{x}、客流场景\omega下由仿真得到的乘客平均候车时间，由于\omega具有随机性，故目标函数中的符号[\cdot ]是关于随机变量\omega的平均候车时间期望值。

3.1.2　约束条件

为保证行车安全和满足服务水平基本要求，设定如下约束条件。

约束1：乘客的平均最大候车时间不应大于某一阈值{\delta }_{\mathrm{1}}，即E[T_{\mathrm{m}\mathrm{a}\mathrm{x}}^{\mathrm{w}\mathrm{a}\mathrm{i}\mathrm{t}}(X，\omega )]\le {\delta }_{\mathrm{1}}。其中，T_{\mathrm{m}\mathrm{a}\mathrm{x}}^{\mathrm{w}\mathrm{a}\mathrm{i}\mathrm{t}}(X，\omega )是通过仿真得到的所有乘客中的最大候车时间。该约束条件是为了避免出现个别乘客候车时间过长的情况，虽然少数候车时间很长的乘客并不会对目标函数带来较大变化，但会非常影响个别乘客对运营企业的感观。

约束2：全日列车的平均最大满载率不应大于某一阈值{\delta }_{\mathrm{2}}，即E[{\theta }_{\mathrm{m}\mathrm{a}\mathrm{x}}(X，\omega )]\le {\delta }_{\mathrm{2}}。其中，{\theta }_{\mathrm{m}\mathrm{a}\mathrm{x}}(X，\omega )是通过仿真得到的全日所有列车中的最大满载率。该约束条件是为了避免个别列车出现极端拥挤的情况。

约束3：各个时段的行车间隔均应满足最小、最大追踪间隔约束。

I_{\mathrm{m}\mathrm{i}\mathrm{n}}\le x_{l+\mathrm{1}}-x_l\le I_{\mathrm{m}\mathrm{a}\mathrm{x}}\forall l=\mathrm{1}，\mathrm{2}，...，L-\mathrm{1}⑹

式中：I_{\mathrm{m}\mathrm{i}\mathrm{n}}和I_{\mathrm{m}\mathrm{a}\mathrm{x}}分别表示列车最小、最大追踪间隔，min。结合《地铁设计规范》对于城市轨道交通列车高峰和平峰的运行间隔的规定，列车运行时车头距不能小于最小追踪间隔2 min，最大运行间隔为10 min。

3.1.3　模型转换

由于模型的约束1和约束2包括仿真指标，因此将这2个约束条件添加至目标函数中，从而消除随机非线性的约束式，转化后的模型如下。

\mathrm{m}\mathrm{i}\mathrm{n}E[f(X，\omega )]=\mathrm{m}\mathrm{i}\mathrm{n}E[T^{\mathrm{w}\mathrm{a}\mathrm{i}\mathrm{t}}(X，\omega )+M_{\mathrm{1}}{P}_{\mathrm{1}}(X，\omega )+M_{\mathrm{2}}{P}_{\mathrm{2}}(X，\omega )]⑺

其中，M_{\mathrm{1}}和M_{\mathrm{2}}为充分大的正数，即惩罚因子；P_{\mathrm{1}}和P_{\mathrm{2}}为指示性函数，计算公式具体如下。

P_{\mathrm{1}}(X,\omega )=\left\{\begin{array}{l}\mathrm{1}{T}_{\mathrm{m}\mathrm{a}\mathrm{x}}^{\mathrm{w}\mathrm{a}\mathrm{i}\mathrm{t}}(X,\omega )>{\delta }_{\mathrm{1}}\\ \mathrm{0}{T}_{\mathrm{m}\mathrm{a}\mathrm{x}}^{\mathrm{w}\mathrm{a}\mathrm{i}\mathrm{t}}(X,\omega )\le {\delta }_{\mathrm{1}}\end{array}\right.⑻

P_{\mathrm{2}}(X,\omega )=\left\{\begin{array}{l}\mathrm{1}{\theta }_{\mathrm{m}\mathrm{a}\mathrm{x}}(X,\omega )>{\delta }_{\mathrm{2}}\\ \mathrm{0}{\theta }_{\mathrm{m}\mathrm{a}\mathrm{x}}(X,\omega )\le {\delta }_{\mathrm{2}}\end{array}\right.⑼

3.2　算法设计

城市轨道交通列车时刻表优化问题已被证明属于NP-Hard问题^[14]，随着问题中列车数量的增加，模型计算复杂度呈指数级增长，无法在多项式时间内枚举出全部可行的发车方案。此外，优化模型的目标函数带有关于随机客流场景\omega的期望符号[\cdot ]，目标函数缺乏解析形式且梯度值无法直接测量^[15]。

针对此模型，设计基于仿真的有限差分随机逼近(Finite Difference Stochastic Approximation，FDSA)算法。该算法首先需要选取一个初始可行解作为优化起点；在每次迭代优化中，通过不断采样随机客流样本、仿真并统计各客流样本在各扰动方案下的目标函数值，从而估计目标函数的梯度，优化决策变量，进而获得优化方案，得到基于仿真的FDSA算法流程如图5所示。令X^{(h)}表示算法在第h次迭代后产生的方案，算法详细步骤如下。

步骤1：设置扰动参数\beta、步长参数c、A、\alpha、最大更新步长{\mathrm{\Delta }}_{\mathrm{m}\mathrm{a}\mathrm{x}}、单次迭代客流样本采样数Z、最大迭代次数为h_{\mathrm{m}\mathrm{a}\mathrm{x}}；设置初始方案X^{(\mathrm{0})}=[x_{\mathrm{1}}^{(\mathrm{0})}，x_{\mathrm{2}}^{(\mathrm{0})}，...，x_l^{(\mathrm{0})}](可根据区间最大断面客流量计算或参照目前实际运营情况)；令迭代次数h=\mathrm{1}。

步骤2：随机采样Z个客流样本\{{\omega }_{\mathrm{1}}，{\omega }_{\mathrm{2}}，\dots， {\omega }_Z\}。

步骤3：对各个列车的发车时刻进行扰动，并在扰动后估计目标函数梯度值，具体如下。

①针对列车l，根据伯努利\pm \mathrm{1}分布按50%概率随机生成方向值\xi \in \{-\mathrm{1}，\mathrm{1}\}，由此获得沿发车时刻x_l方向的扰动方案X_l\mathrm{\text{'}}=[x_{\mathrm{1}}^{(h-\mathrm{1})}，x_{\mathrm{2}}^{(h-\mathrm{1})}，\dots ，x_l^{(h-\mathrm{1})}+\beta \cdot \xi ，\dots ]。对扰动后的方案进行可行性检验，即判断调整后的发车方案是否满足该线路上最大最小发车间隔约束，若不满足则令\xi =-\xi后重新计算扰动方案X_l\mathrm{\text{'}}=[x_{\mathrm{1}}^{(h-\mathrm{1})}，x_{\mathrm{2}}^{(h-\mathrm{1})}，\dots ，x_l^{(h-\mathrm{1})}+\beta \cdot \xi ,\dots ]，并再次判断，若仍不满足则令\xi =\mathrm{0}，计算扰动方案为X_l\mathrm{\text{'}}=[x_{\mathrm{1}}^{(h-\mathrm{1})},x_{\mathrm{2}}^{(h-\mathrm{1})},\dots ,x_l^{(h-\mathrm{1})}+\beta \cdot \xi ,\dots ]。仿真得到扰动方案x_l\mathrm{\text{'}}在各随机客流样本\omega \in \{{\omega }_{\mathrm{1}}，{\omega }_{\mathrm{2}}，\dots ，{\omega }_Z\}下的目标函数值f(X_l\mathrm{\text{'}}，{\omega }_Z)，并进一步统计平均值\overline{f}(X_l\mathrm{\text{'}})={\displaystyle \sum _{z=\mathrm{1}}^Z}f(X_l\mathrm{\text{'}}，{\omega }_Z)/Z。

②用g_l表示x_l方向上的梯度估计值，其计算公式为

g_l=\frac{\overline{f}(X_l\mathrm{\text{'}})-\overline{f}(X^{(h)})}{\beta \cdot \xi }⑽

步骤4：根据估计梯度值更新决策变量，分以下4个子步骤进行。

①令更新步长a^{(h)}=c/{(A+h)}^{\alpha }，更新各决策变量

{\widehat{x}}_l^{(h+\mathrm{1})}=x_l^{(h)}-\mathrm{m}\mathrm{a}\mathrm{x}\{\mathrm{m}\mathrm{i}\mathrm{n}\{g_l\cdot a^{(h)}，{\mathrm{\Delta }}_{\mathrm{m}\mathrm{a}\mathrm{x}}\}，-{\mathrm{\Delta }}_{\mathrm{m}\mathrm{a}\mathrm{x}}\}⑾

②将更新后的解映射至可行域，即保证新解满足约束公式⑹。此外，应尽量使更新解在向可行域映射时的调整量最小，由此构建模型如下。

\mathrm{m}\mathrm{i}\mathrm{n}{\displaystyle \sum _{l=\mathrm{1}}^L}\left|{\widehat{x}}_l^{(h+\mathrm{1})}-x_l^{(h+\mathrm{1})}\right|⑿

s.t. I_{\mathrm{m}\mathrm{i}\mathrm{n}}\le x_{l+\mathrm{1}}^{(h+\mathrm{1})}-x_l^{(h+\mathrm{1})}\le I_{\mathrm{m}\mathrm{a}\mathrm{x}}\forall l=\mathrm{1}，\mathrm{2}，...，L-\mathrm{1} ⒀

其中x_l^{(h+\mathrm{1})}表示向可行域映射后得到的新解。

③将上述模型线性化从而方便求解，引入中间变量{\mathrm{y}}_l=\left|{\widehat{x}}_l^{(h+\mathrm{1})}-x_l^{(h+\mathrm{1})}\right|，将原模型转化成如下。

\mathrm{m}\mathrm{i}\mathrm{n}{\displaystyle \sum _{l=\mathrm{1}}^L}{\mathrm{y}}_l⒁

s.t. \left\{\begin{array}{l}{I}_{\mathrm{m}\mathrm{i}\mathrm{n}}\le x_{l+\mathrm{1}}^{(h+\mathrm{1})}-x_l^{(h+\mathrm{1})}\le I_{\mathrm{m}\mathrm{a}\mathrm{x}}\forall l=\mathrm{1}，\mathrm{2}，\dots ，L-\mathrm{1}\\ {\widehat{x}}_l^{(h+\mathrm{1})}-x_l^{(h+\mathrm{1})}\le y_l\forall l=\mathrm{1}，\mathrm{2}，\dots ，L\\ -{\widehat{x}}_l^{(h+\mathrm{1})}+x_l^{(h+\mathrm{1})}\le y_l\forall l=\mathrm{1}，\mathrm{2}，\dots ，L\end{array}\right. ⒂

④对新解仿真得到其在各随机客流样本\omega \in \{{\omega }_{\mathrm{1}}，{\omega }_{\mathrm{2}}，\dots ，{\omega }_Z\}下的目标函数值f(X^{(h+\mathrm{1})}，{\omega }_Z)，并进一步统计平均值\overline{f}(X^{(h+\mathrm{1})})={\displaystyle \sum _{z=\mathrm{1}}^Z}f(X^{(h+\mathrm{1})}，{\omega }_Z)/Z。

步骤5：判断迭代终止条件。若h\ge h_{\mathrm{m}\mathrm{a}\mathrm{x}}则算法结束；否则，令h=h+\mathrm{1}，返回步骤2。

4　算例分析

4.1　相关参数

为验证仿真模型及优化算法效果，以国内某城市的某条地铁线路为背景进行算例分析。该线共包含25个站点，仿真模型相关参数如表2所示。根据2021年10月工作日的历史客流数据，以15 min为间隔划分全日研究时段，统计各时段内各OD的平均客流量，并将其拟合为泊松分布模型。

优化算法参数如表3所示。选取该城市轨道交通运营公司在该线路上现行的列车发车时刻表作为初始方案。仿真模型及有限差分随机逼近算法均采用C#.NET编程实现，优化算法步骤3中②的梯度估计采用并行计算，步骤4中③的线性规划模型通过调用Gurobi进行求解。测试环境的CPU型号为Intel i7-1165G7，内存为40 GB，操作系统为64位Windows 10。

4.2　结果分析

经测试，算法平均一次迭代所需时间约30 s，每次仿真所需运行时间为0.4 s。算法迭代500次后基本收敛，目标函数迭代过程如图6所示。在不增加发车频次的情况下，相较于初始方案，乘客平均等待时间从175.4 s减少到163.6 s，降低6.75%。

对该线路全日单方向180个车次优化前后的发车间隔进行对比，得到优化前后发车间隔对比如图7所示。可以看出，在保持现有的发车能力条件下，优化后的发车方案会根据客流需求的随机性和动态性调整发车间隔，使其呈现出非均匀间隔的特征。此外，在过渡段和高峰时段的发车间隔调整量小于平峰时段。

分别仿真初始方案和优化后方案的列车运行过程，得到优化前后的仿真实际运行图对比如图8所示。其中，红色线表示优化后的城市轨道交通列车运行图，而黑色线表示初始方案下的列车运行图。从图8中可以清晰地看出，优化后的方案将更多的运力配置到高峰时段，故其他部分时段(如9：30—11：30)的行车间隔有一定程度增加，但整体上优化后的运行图和客流需求更加匹配。

5　结束语

针对城市轨道交通列车时刻表优化问题，考虑客流需求的随机性与动态性，结合列车容量限制、列车停站时间及区间运行时间实时调整等运营特征，构建了城市轨道交通系统仿真模型。在此基础上，提出列车时刻表优化模型，以最小化乘客候车等待时间；并设计基于仿真的有限差分随机逼近优化算法求解。通过算例验证表明，经过本模型优化后的列车时刻表，可更好地适应客流需求的动态随机性，有效降低乘客的平均候车等待时间，提高城市轨道交通服务水平。在下一步的研究工作中，将考虑从单条线路优化延伸到网络层面的城市轨道交通线网优化，并着重考虑基于随机性换乘客流需求下的列车时刻表协同优化问题。

参考文献

[1] 中国城市轨道交通协会.城市轨道交通2022年度统计和分析报告[R].北京：中国城市轨道交通协会，2023.

[2] 孙晟凯，李思杰，刘志钢，等. 基于客流与列车动态交互的多站协同客流管控策略评价[J]. 铁道运输与经济，2021，43(12)：116-125.

[3] SUN Shengkai， LI Sijie， LIU Zhigang，et al. Evaluation on Multi-Station Collaborative Passenger Flow Control Strategy Based on Dynamic Interaction between Passenger Flow and Trains[J]. Railway Transport and Economy，2021，43(12)：116-125.

[4] 史　峰，涂　纤，赵　烁. 面向O-D时变需求的城轨列车运行计划优化方法[J]. 中国铁道科学，2020，41(1)：115-125.

[5] SHI Feng， TU Xian， ZHAO Shuo. Optimization Method of Urban Rail Train Operational Plan Based on O-D Time-Varying Demand[J]. China Railway Science，2020，41(1)：115-125.

[6] NIU H M， ZHOU X S， GAO R H. Train Scheduling for Minimizing Passenger Waiting Time with Time-Dependent Demand and Skip-Stop Patterns：Nonlinear Integer Programming Models with Linear Constraints[J]. Transportation Research Part B：Methodological，2015(76)：117-135.

[7] YIN J T， YANG L X， TANG T，et al. Dynamic Passenger Demand Oriented Metro Train Scheduling with Energy-Efficiency and Waiting Time Minimization：Mixed-Integer Linear Programming Approaches[J]. Transportation Research Part B：Methodological，2017(97)：182-213.

[8] GONG C C， SHI J G， WANG Y H，et al. Train Timetabling with Dynamic and Random Passenger Demand：A Stochastic Optimization Method[J]. Transportation Research Part C：Emerging Technologies，2021(123)：102963.

[9] HASSANNAYEBI E， SAJEDINEJAD A， MARDANI S. Urban Rail Transit Planning Using a Two-Stage Simulation-Based Optimization Approach[J]. Simulation Modelling Practice and Theory，2014(49)：151-166.

[10] 张春田，戚建国，杨立兴，等. 基于不确定旅客需求的高速铁路鲁棒列车开行方案研究[J]. 交通运输系统工程与信息，2022，22(1)：115-123.

[11] ZHANG Chuntian， QI Jianguo， YANG Lixing，et al. Robust Train Operation Plan Based on Uncertain Passenger Demands for High Speed Railway Corridors[J]. Journal of Transportation Systems Engineering and Information Technology，2022，22(1)：115-123.

[12] HAN Z Y， HAN B M， LI D W，et al. Train Timetabling in Rail Transit Network under Uncertain and Dynamic Demand Using Advanced and Adaptive NSGA-II[J]. Transportation Research Part B：Methodological，2021(154)：65-99.

[13] HASSANNAYEBI E， ZEGORDI S H， YAGHINI M，et al. Timetable Optimization Models and Methods for Minimizing Passenger Waiting Time at Public Transit Terminals[J]. Transportation Planning and Technology，2017，40(3)：278-304.

[14] LEI Y Z， LU G Y， ZHANG H X，et al. Optimizing Total Passenger Waiting Time in an Urban Rail Network：A Passenger Flow Guidance Strategy Based on a Multi-Agent Simulation Approach[J]. Simulation Modelling Practice and Theory，2022，117：102510.

[15] ZHU Y T， MAO B H， LIU L，et al. Timetable Design for Urban Rail Line with Capacity Constraints[J]. Discrete Dynamics in Nature and Society，2015：429219.

[16] 周　明. 上海地铁停站时间影响因素分析及优化[J]. 交通世界，2018(9)：166-168.

[17] ABDOLMALEKI M， MASOUD N， YIN Y F. Transit Timetable Synchronization for Transfer Time Minimization[J]. Transportation Research Part B：Methodological，2020，131 ：143-159.

[18] ZHANG Y， YAN H， LUO Y，et al. Simulation-based Optimization of Timetables Coordination in an Urban Rail Transit Network[J]. Simulation Modelling Practice and Theory，2024，130： 102857.

基金资助

国家自然科学基金项目(72101218)

四川省科技计划资助项目(2023NSFSC1033)

四川省科技计划资助项目(2023NSFSC0906)