Performance improvement of energy consumption, passenger time and robustness in metro systems: A multi-objective timetable optimization approach

Computers & Industrial Engineering 137 (2019) 106076

Contents lists available at ScienceDirect

Computers & Industrial Engineering journal homepage: www.elsevier.com/locate/caie

Performance improvement of energy consumption, passenger time and robustness in metro systems: A multi-objective timetable optimization approach

T

Xin Yang, Jianjun Wu , Huijun Sun, Ziyou Gao, Haodong Yin, Yunchao Qu ⁎

State Key Laboratory of Rail Traffic Control and Safety, Beijing Jiaotong University, Beijing 100044, China

ARTICLE INFO

ABSTRACT

Keywords: Metro systems Timetable optimization Energy consumption Passenger time Robustness

As an effective transport mode to reduce the impact of urban mobility with great capacity and low pollution, metro systems have received rapid development all over the world in recent years. Since a well-designed timetable can improve the service level of passengers, increase the revenue of operating companies, and reduce the environmental pollution, the issue of timetable design/optimization is usually represented as a multi-objective decision problem. This paper develops a multi-objective integer programming model which integrates the energy consumption, the passenger waiting time, and the robustness by optimizing departure and arrival time of trains at each station. A non-dominated sorting genetic algorithm II (NSGA-II) is designed to find the optimal solution. Finally, a case study of Yizhuang line in Beijing metro system is conducted to illustrate the efficiency and effectiveness of the suggested timetable optimization approach.

1. Introduction Metro system has received fast development all over the world due to its great capacity, high reliability and low pollution. In some metropolises such as London, New York and Beijing, the metro system has become a large-scale complex network (Yang, Yin, Wu, Qu, & Gao, 2019). The network mileages and number of lines and stations for some big cities around the world until the end of 2017 are listed in Table 1 (ChinaDaily, 2017; Wikipedia, 2018; Yang, Chen, Wu, Gao, & Tang, 2019). From the operation perspective, the networking operation has taken the place of the traditional single-line operation mode. Under the background of networking operation, as a crucial element of a metro system, timetable bears much more significance in terms of transport efficiency, operation cost and service quality. Although per capita energy consumption in metro systems is low, the total energy consumption is considerable (Cao, Wang, Liu, Li, & Xie, 2019). For example in Beijng, the Beijing Metro is the first large industry in electricity consumption. Minimization of the energy consumption is important to reduce the operation cost of metro operating companies. In the other hand, the essential mission of metro systems is transporting passengers efficiently and conveniently. Reduction of passenger travel time should be also taken into consideration. Therefore, it is significant to study a well-designed timetable considering the benefits of both metro passengers and operating companies. A number of existing ⁎

literature on the timetable design/optimization problem considered both passenger time and energy consumption as research objectives (Ghoseiri, Szidarovszky, & Asgharpour, 2004; Li, Wang, Li, & Gao, 2013; Xu, Li, & Li, 2016; Yang, Li et al., 2015; Yang, Ning, Li, & Tang, 2014; Yin, Yang, Tang, Gao, & Ran, 2017). In real-world metro systems, trains are stipulated to complete the travel task strictly according to the published timetable. However, there are a number of uncertain perturbations during the normal operation period (Cao, Ma, Xiao, Zhang, & Xu, 2017; Corman, D’Ariano, Marra, Pacciarelli, & Samà, 2017; Zhang, Li, & Yang, 2018). For example, crowding passengers at stations can lead to the delay of departure time; and trains often have an early departure perturbation at off-peak hours due to the over abundant dwell time. In these cases, the actual train record in time level (i.e., recorded timetable) has a disordered deviation in comparison with the published timetable. Fig. 1 provides an illustration of deviations between recorded and published timetables according to the real-world data from the Beijing Metro Yizhuang line (Zhang, 2018). Both energy consumption and passenger time are sensitive to the spatially and temporally relative position of trains scheduled by timetable. Therefore, it is difficult to achieve the optimal expected objectives under the condition of these perturbations, although the published timetable is the optimal one. For this issue, we introduce the robustness in the published timetable optimization problem to deal with the

Corresponding author. E-mail address: [email protected] (J. Wu).

https://doi.org/10.1016/j.cie.2019.106076 Received 16 October 2018; Received in revised form 11 July 2019; Accepted 17 September 2019 Available online 18 September 2019 0360-8352/ © 2019 Elsevier Ltd. All rights reserved.

Computers & Industrial Engineering 137 (2019) 106076

X. Yang, et al.

Table 1 The scale of metro systems for some metropolises around the world. City

Shanghai

Beijing

London

New York

Moscow

Seoul

Tokyo

Paris

Mileage Line Station

617 14 303

574 19 288

402 13 277

368 24 468

327 12 171

314 9 376

302 13 277

293 17 379

#Note. Units: mileage (kilometers), line (number), station (number).

Fig. 1. The illustration of deviations between recorded and published timetables.

2. Literature review

uncertain perturbations. Note that the traditional research on robustness in timetable problems aims to absorb perturbations for improving the punctuality rate. It is easily see that the robustness in this paper is a different concept, which is used to eliminate or reduce the influence of perturbations to energy consumption and passenger time. Overall, this paper develops a multi-objective integer programming model and a non-dominated sorting genetic algorithm II (NSGA-II) to improve the performance of energy consumption, passenger waiting time and robustness for metro systems. The main contributions of this paper are summarized as follows:

The timetable optimization problem for metro/rail systems has been widely studied due to the good achievements in improving operational quality. The literature reviewed in this section includes single-objective optimization to improve the performance of energy consumption/passenger time/robustness, and multi-objective optimization approaches. 2.1. Energy consumption Timetable optimization is considered as a preferential approach to improve the benefits of regenerative braking in metro systems for energy conservation. For example, Ramos, Peña-Alcaraz, FernándezCardador, and Cucala (2007) proposed a timetable optimization method to increase the utilization of regenerative braking energy by maximizing the overlapping time between adjacent trains. Domínguez, Fernández-Cardador, Cucala, and Pecharromán (2012) further developed a power flow model combined with the timetable optimization model to minimize the net energy consumption. Ye and Liu (2016) established a multiphase energy-efficient control method incorporating complex train running conditions and real-world train operation constraints to determine the optimal timetable and speed profile. Zhou, Tong, Chen, Tang, and Zhou (2017) formulated a unified modeling framework using space-time-speed grid networks to jointly optimize train timetable and speed profile for minimizing the energy consumption. Canca and Zarzo (2017) developed a mixed integer optimization model to determine an energy-efficient timetable with consideration of the passenger loading distribution. Huang et al. (2018) investigated an integrated timetable and driving strategy optimization approach for reducing the energy consumption of multiple trains by considering regenerative braking energy. Based on the practical data from the Beijing Metro Yizhuang line in China, the research group from Beijing Jiaotong University (BJTU) did a series of studies on energy-efficient timetable optimization problem in recent years. For example, Yang, Li, Gao, Wang, and Tang (2013) proposed an energy-efficient cooperative scheduling model to synchronize the braking and accelerating processes of adjacent trains for

(1) We firstly use the concept of robustness to tackle the uncertain dwell time in real-world train operations for improving the performance of timetable optimization approaches applied to realworld metro systems. (2) We develop a multi-objective integer programming model to solve the timetable optimization problem with dwell time uncertainty, where the energy consumption, passenger waiting time, as well as robustness are optimized of compromise. Traditional studies assumed the passenger arrival rate as a uniform distribution, such that the average passenger waiting time equals to a half of headway. As we use the real-world passenger arrival rate (i.e., a random distribution) in the formulation, we can capture the wispy sensibility of dwell time and section running time to the passenger waiting time. (3) We design a detailed NSGA-II procedure to solve the complicated multi-objective integer programming model, and come up with some indicators to evaluate the performance of the Pareto-optimal solutions. The remainder of this research is organized as follows. In Section 2, we review the literature on the timetable optimization problem for metro/rail systems. In Section 3, we formulate the problem as a multiobjective integer programming model. In Section 4, we design a NSGAII to find the optimal solution. In Section 5, we present a case study based on the real-world data from the Beijing Metro Yizhuang Line. Conclusions are provided in Section 6. 2

Computers & Industrial Engineering 137 (2019) 106076

X. Yang, et al.

increasing the recovered energy. On this basis, Li and Yang (2013) extended the model with consideration of the stochastic departure time delay for trains at busy stations to reduce the total energy consumption. Taking account of the speed profile, Li and Lo (2014a) proposed an integrated speed profile and timetable model to achieve a better performance on the net energy consumption. Furthermore, Li and Lo (2014b) developed a dynamic integrated model with adaptive headway and cycle time to determine the energy-efficient speed profile and timetable based on the passenger demand. Yang, Chen et al. (2015) presented an energy-efficient scheduling approach to increase the utilization of recovered energy, where the number of trains and the cycle time were kept unchanged in the optimization model. Yang, Chen et al. (2016) developed a two-phase stochastic programming model to determine the timetable and speed profile with uncertain train mass for minimizing the total tractive energy consumption. The more train timetable optimization approaches to reduce energy consumption in metro/rail systems with the analysis of their advantages and dis advantages can be found in Yang, Li et al. (2016) and Scheepmaker, Goverde, and Kroon (2017).

stage optimization model for determining robust rolling stock circulations for passenger trains. Benders decomposition was designed to solve the model in a short time for a representative number of disruption scenarios. Sels, Dewilde, Cattrysse, and Vansteenwegen (2016) formulated an efficient and robust timetable optimization model to minimize the excepted passenger time, including running time, dwell time, transfer time, the first delay time and the second delay time. Lu, Tang, Zhou, Yue, and Huang (2017) developed an integrated two-stage approach to consider the recovery-to-optimality robustness into the optimized timetable design without predefined structure information. Jovanović, Kecman, Bojović, and Mandić (2017) studied the optimal allocation of buffer times for improving the timetable robustness to ensure the reliability and punctuality of rail systems. Lee, Lu, Wu, and Lin (2017) formulated a simulation-based approach to adjust the time supplement and buffer time in a given passenger railway timetable to reduce the average delay. The proposed approach can quantify the relationship between timetable efficiency and robustness.

2.2. Travel time

As different stakeholders with different benefits and interests are involved, it is reasonable to consider the timetable optimization process as a multi-objective decision problem. For example, Higgins, Kozan, and Ferreira (1996) proposed a bi-objective mixed-integer programming model to reduce both train operation cost and delay time, and designed a branch and bound algorithm to solve the integer program. Ghoseiri et al. (2004) formulated a multi-objective optimization model for lowering the total passenger time and the fuel consumption cost, and applied the ε-constraint method to determine the Pareto frontier. Corman, D’Ariano, Pacciarelli, and Pranzo (2012) developed a bi-objective conflict detection and resolution optimization model to maximize the total value of satisfied connections and minimize the consecutive delays between trains. Li et al. (2013) established a multiobjective scheduling model to minimize the energy consumption, the passenger time and the carbon emission cost applying a fuzzy programming algorithm for obtaining the optimal solution. Yang et al. (2014) constructed a bi-objective integer programming method with dwell time and headway control to improve the utilization of recovered energy and shorten the total passenger waiting time at all stations. On the other hand, Yang, Li et al. (2015) developed another bi-objective optimization method to determine both the speed profile and timetable, and designed an optimal train control algorithm combined with an adaptive genetic algorithm to solve the formulation. Xu et al. (2016) developed a bi-objective programming model to optimize the all-day timetable with consideration of both service quality and energy efficiency. Yin et al. (2017) developed two mixed-integer linear programming models to minimize the energy consumption and passenger waiting time with consideration of dynamic passenger demand in metro systems. To clarify the contributions of our paper, the detailed features of some closely related studies in the literature are listed in Table 2.

2.4. Multi-objective optimization

Travel time in metro systems mainly includes passenger waiting time at stations, passenger time on the train, and transfer time. Passenger waiting time is focused in this paper. Some early studies calculated the passenger waiting time by various analytical formulas in order to derive optimal scheduling strategies. Vansteenwegen and Van Oudheusden (2007) proposed a two-phase timetable optimization model to decrease the passenger waiting time, where ideal running time supplements and continuous linear programming are applied to find the optimal solution. Niu and Zhou (2013) constructed a binary integer programming model incorporated with departure events and passenger loading to optimize the timetable under time-dependent demand in oversaturated conditions. Barrena et al. (2014a) presented three timetable optimization models to reduce passenger waiting time and designed a branch-and-cut algorithm to obtain the optimal solutions. Furthermore, considering a dynamic passenger demand, Barrena et al. (2014b) developed a non-periodic train timetabling model to minimize the average passenger waiting time. Niu, Zhou, and Gao (2015) developed a nonlinear integer programming model to minimize passenger waiting time with consideration of time-dependent demand and skipstop patterns. Tian and Niu (2019) developed a nonlinear integer programming model and a dynamic programming algorithm to synchronize train timetables for reducing the passenger time. 2.3. Robustness Traditional robustness concept in metro/rail systems has been studied for a long time and achieved some good performance. For example, Fischetti, Salvagnin, and Zanette (2009) proposed four different approaches to improve the robustness of a given train timetabling solution for the noncyclic timetable problem. Salido, Barber, and Ingolotti (2012) presented analytical and simulation methods based on the infrastructure topology and buffer times to measure robustness in a single railway line from the point of view of railway operators. Cacchiani and Toth (2012) surveyed the main studies dealing with the train timetabling problem in its nominal and robust versions. Periodic schedule and non-cyclic schedule were the two main variants of the nominal problem. In robust version of the problem, the aim was to avoid delay propagation as much as possible. Cacchiani, Caprara et al. (2012) firstly proposed a simple modification of the Lagrangian optimization scheme capable to deal with the train timetable robustness. Based on the realworld test cases from the Italian Railways, it showed that the modified Lagrangian Heuristic could produce within much shorter computing time solutions whose quality was comparable with those produced by existing approaches. Cacchiani, Caprara et al. (2012) formulated a two-

3. Model formulation This section proposes a multi-objective integer programming model to improve the performance of energy consumption, passenger waiting time and robustness by optimizing departure and arrival time of trains at each station and headway time. The following presentations focus on detailing each part of the formulation, i.e., notions, model assumptions, objective functions, constraints, optimization model and convergence property. 3.1. Notations We define a one-way metro line G = (N, A), where N is a finite set of stations and A is a finite set of sections between adjacent stations. We 3

Computers & Industrial Engineering 137 (2019) 106076

X. Yang, et al.

Table 2 The detailed features of closely related studies. Publication

Dwell time

Passenger arrival rate

Main contributions

Ghoseiri et al. (2004)

Certain

Not considered

Yang et al. (2014) Xu et al. (2016)

Certain Certain

Uniform distribution Uniform distribution

This paper

Uncertain

Real-world passenger arrival rate

Firstly formulated a multi-objective timetable optimization model with consideration of both energy consumption and passenger time Firstly considered the regenerative braking energy in the multi-objective timetable optimization model Firstly discussed the effect of passenger arrival rate to passenger time in the multi-objective timetable optimization model Firstly use the concept of robustness to tackle the uncertain dwell time, and study the real-world passenger arrival rate in the multi-objective timetable optimization model

Fig. 2. Illustration of a one-way metro line.

use (n, n + 1) to denote the section between station n and station n + 1 as shown in Fig. 2. Notations used throughout the paper are listed as follows and all boldface letters denote the corresponding vectors, and some of them are indicated in Figs. 2 and 3.

(A1). We only consider a one-way metro line in this paper, as we study trains running on the same direction. In addition, metro timetable is considered as cyclic in operation, i.e., the same speed sequence is assumed for all trains on the same section, but it is different for different sections. (A2). The process of a train running between adjacent stations includes three phases as shown in Fig. 3, i.e., maximum accelerating [din, cin), coasting [cin, bin), and maximum braking [bin, ai(n+1)). (A3). Some parameters such as train mass m, maximal tractive force fa, maximal braking force fb, basic running resistance r and additional running resistance g are considered as constant values. (A4). The number of arrival passengers to station n at time t is randomly generated in advance based on the real-world historical data.

(1) Indices i index of a train, i = 1, 2, …, I n index of a station, n = 1, 2, …, N (2) Parameters m train mass fa maximal tractive force maximal braking force fb r basic running resistance g additional running resistance by metro line curve and gradient η1 conversion efficiency during the accelerating process η2 conversion efficiency during the braking process number of arrival passengers at station n and time t pn(t) qin number of passengers left by train i − 1 to wait train i at station n length of section (n, n + 1) l(n, n+1) t1(n, n+1) minimum running time of a train on section (n, n + 1) t2(n, n+1) maximum running time of a train on section (n, n + 1) T1 minimum value of total travel time from origin to terminal station T2 maximum value of total travel time from origin to terminal station current value of total service time T3 t1n minimum dwell time at station n t2n maximum dwell time at station n minimum headway time hl hu maximum headway time

3.3. Objective functions Three objectives with randomness of dwell time are considered in the formulation, i.e., energy consumption is denoted by the expected value of total net energy consumption; passenger waiting time is denoted by the expected value of total time of passengers waiting at all stations; and robustness is represented by the mean square error of total net energy consumption and total passenger waiting time. (1) Energy consumption

(3) Random variables ξin dwell time of train i at station n, i.e., t1n ≤ ξin ≤ t2n ξ dwell time set, i.e., ξ = {ξin|1 ≤ i ≤ I, 1 ≤ n ≤ N}

Considering a section (n, n + 1), the departure time stamp of train i at station n is din, and the arrival time stamp of train i at station n + 1 is ai(n+1). According to the equations of motion, for each 1 ≤ n ≤ N − 1 and 1 ≤ i ≤ I, the speed of train i at time stamp t is formulated as

(4) Decision variables ain arrival time stamp of train i at station n a arrival time stamp set, i.e., a = {ain|1 ≤ i ≤ I, 1 ≤ n ≤ N} departure time stamp of train i at station n, i.e., din = ain + ξn din d departure time stamp set, i.e., d = {din|1 ≤ i ≤ I, 1 ≤ n ≤ N} hi headway time of train i and train i + 1 h headway time set, i.e., h = {hi|1 ≤ i ≤ I} (5) Intermediate Variables cin switching time stamp from accelerating to coasting phase of train i at section (n, n + 1) bin switching time stamp from coasting to braking phase of train i at section (n, n + 1) vcin speed value at the switching point from accelerating to coasting phase at section (n, n + 1) vbin speed value at the switching point from coasting to braking phase at section (n, n + 1) vin speed sequence of train i at section (n, n + 1)

3.2. Model assumptions To simplify the model formulation and solution procedure, we make the following assumptions according to some real-world operation properties of metro systems.

Fig. 3. Illustration of some parameters and variables.

4

Computers & Industrial Engineering 137 (2019) 106076

X. Yang, et al.

vin (t ) =

(fa

r + g )(t

vcin vbin

(r g )(t cin)/ m , if cin t < bin , (fb + r g )(t bin)/m , if bin t < ai (n + 1) ,

din)/m ,

if din

0 min{cin

t < cin, O1 (a , d, h , ) =

(1)

r + g )(cin

din)/ m ,

vbin = (fa

r + g )(cin

din )/ m

(r

g )(bin

cin)/m .

(fb + r

2 vcin /(fa

g )(ai (n + 1)

2 r + g ) + (vcin

2 g ) + vbin /(fb + r

0, if din t < bin , fb vin (t ) 2 , if bin t < ai (n + 1).

Fin (t ),

I

t = din

Bin (t ). i=1 n=1

t = din

c(i + 1) n ai (n + 1) . d (i + 1) n < ai (n + 1) ai (n + 1) < d (i + 1) n

(O1 (a , d , h , ) + O2 (a , d, h , )). i=1 n=1

g ) = 2l (n - 1, n) / m ,

(10)

We take the ratio of overlapping time and braking time as the conversion efficiency of the braking process, i.e., 3 (a ,

d , h, ) =

O (a , d , h , ) I 1N 1 i=1 n=1

(ai (n + 1)

. bin)

(11)

As the difference between the total required tractive energy and regenerative energy, the expected value of total net energy consumption is formulated as

JE (a, d, h, ) = E [JF (a, d , )

3 (a ,

d , h, ) JB (a , d , )]

(12)

This objective function of energy consumption is nonlinear. (2) Passenger waiting time

(4)

Considering a station n, the number of arrival passengers pn(t) at time t can be obtained based on the real-world smart card data. Set a referential time stamp d(i−1)n, i.e., the departure time stamp of train i − 1 at station n. As the capacity limitation of train i − 1, the number of passengers left at station n to wait train i is defined as qin. Note that qin is equal to zero during off-peak hours due to the few passengers. When train i arrives at station n, i.e., at the time stamp ain, the accumulated number of passengers left at station n is formulated as

(5)

ain

Qin = qin +

pn (t )

(13)

t = d(i 1) n + 1

Further, for each 1 ≤ i ≤ I and 1 ≤ n ≤ N − 1, the total time of passengers to wait train i at station n can be formulated as

Pin = Qin (ain = qin +

d (i

1) n )

ain t = d (i 1) n + 1

pn (t ) (ain

d (i

1) n )

(14)

In traditional researches, passengers are assumed to arrive at stations evenly in Fig. 5(a), in which passenger waiting time is only related to headway time and dwell time. In this paper, the number of passengers arriving at each station is related to running time in Fig. 5(b). The arrival time and departure time will vary with the change of running time, which will lead to a change in the number of passengers arriving at stations during the period. The passenger waiting time will change as well, that is, it will be affected by running time. Therefore, the total passenger waiting time during the period of I trains for completing trips from station 1 to station N is formulated as

(6)

N 1 ai (n + 1) 1

JB (a, d , h, ) =

d (i + 1) n, c(i + 1) n bin } bin bin , ai (n + 1) d (i + 1) n } bin

I 1 N 1

N 1 ai (n + 1) 1

i=1 n=1

.

c(i + 1) n < bin

O (a , d , h , ) =

With the consideration of random variables ξ (i.e., for each 1 ≤ i ≤ I and 1 ≤ n ≤ N − 1, we have ξin = din − ain), the total required tractive energy and regenerative energy of I trains for completing trips from station 1 to station N are formulated as

JF (a , d , h , ) =

a(i + 1) n

The total overlapping time from the time when the first train leaves the first station to the last train arrives the last station is

Remark 1. The tractive energy is only consumed during the accelerating phase, i.e., conversion efficiency η1 denotes the percentage of electricity to mechanical energy. The regenerative energy is only generated during the braking phase, i.e., conversion efficiency η2 denotes the percentage of mechanical energy to regenerated energy.

I

din

1)

(9)

(2)

The generated regenerative energy of train i at section (n, n + 1) at time stamp t is

Bin (t ) =

b(i + 1)(n

a(i + 1) n < din

min{c(i + 1) n O 2 (a , d , h , ) = min{ai (n + 1) 0

in which the first equation means after braking phase the speed of train i will reduce to zero from vbin; and the second equation is a distance constraint corresponding to speed sequence. Through Eqs. (2) and (3), the four intermediate variables cin, bin, vcin and vbin could be expressed by four formulas related to decision variables and given parameters, respectively. Combined them with Eq. (1), it is seen that the speed sequence vin can also be expressed with the parameters and decision variables. For each 1 ≤ i ≤ I and 1 ≤ n ≤ N − 1, the required tractive energy of train i at section (n, n + 1) at time stamp t is

fa vin (t )/ 1 , if din t < cin, 0, if cin t < ai (n + 1) .

a (i + 1) n

cin < b(i + 1)(n 1) b(i + 1)(n 1) cin a (i + 1) n

0

(3)

Fin (t ) =

1) ,

din

(8)

bin)/m = 0,

2 vbin )/(r

b(i + 1)(n

1) }

0

Furthermore, the switching speeds vci(n−1) and vbi(n−1) should satisfy the following equations:

vbin

min{a (i + 1) n

b(i + 1)(n

din }

where the three rows represent the speed sequence during the maximum accelerating phase, the coasting phase, and the maximum braking phase, respectively. Therefore, we define the speed sequence of train i at section (n, n + 1) as vin = {vin(t)|din ≤ t < ai(n+1)}. The switching speeds vcin and vbin at the beginning and ending points of coasting phase are determined as

vcin = (fa

din, cin

(7)

As shown in Fig. 4, only the regenerative energy generated during the shadow parts will be utilized and the rest regenerative energy will be wasted. The overlapping time of Zone 1 is shown in Eq. (8), which denotes that the train i is accelerating and the train i + 1 is braking at the same time. The overlapping time of Zone 2 is shown in Eq. (9), which represents that the train i + 1 is accelerating and the train i is braking at the same time.

I

N 1

J p (a , d , h , ) =

Pin. i=1 n =1

5

(15)

Computers & Industrial Engineering 137 (2019) 106076

X. Yang, et al.

Fig. 4. Overlapping time between accelerating and braking trains.

With the consideration of random variables ξ, the objective function of passenger waiting time is represented by the expected value of total time of passengers waiting at all stations, i.e.,.

maximum running times on this section, i.e.,

t1(n, n + 1)

t2(n, n + 1),

1

i

I, 1

n

N

1.

(18)

For each 1 ≤ i ≤ I and 1 ≤ n ≤ N − 1, the difference between departure time stamp and arrival time stamp should be the dwell time, i.e.,

(3) Robustness In this formulation, improving the robustness is defined to improve the performance of energy consumption and passenger waiting time under the condition of random dwell times. Therefore, we use the sum of mean square errors of total net energy consumption and total passenger waiting time to denote the robustness value. The robustness is formulated as

JR (a, d, h, ) = V [JE (a, d, h, )] + V [JP (a, d, h, )]

din

(2) Dwell time constraints

(16)

JP (a , d , h , ) = E [Jp (a, d, h, )].

ai (n + 1)

din

ain =

n,

1

i

I,

1

n

N

(19)

1

(3) Headway constraints With the cyclic timetable, the time space for all adjacent trains should be valued between the maximum and minimum headway time, i.e.,

(17)

hl < hi < hu ,

3.4. Constraints

1

i

I

1,

1

n

N

1.

(20)

(4) Total travel time constraints

The formulation must satisfy a series of constraints to ensure that the optimal timetable can achieve a safe and efficient train operation. For clarity, we set the time stamp of first train arriving at first station is zero, i.e., a11 = 0 .

For the optimal timetable, the expected difference between the arrival time stamps at terminal point (i.e., station N) and starting point (i.e., station 1) should be valued between the minimum and maximum total travel time, i.e.,

(1) Running time constraints

T1

For each 1 ≤ i ≤ I and 1 ≤ n ≤ N − 1, the running time of train i on section (n, n + 1) should be valued between the minimum and

E (aiN

ai1)

T2

(5) Total service time constraints

Fig. 5. Sensibility of section running time to the passenger waiting time. 6

(21)

Computers & Industrial Engineering 137 (2019) 106076

X. Yang, et al.

Fig. 6. Flow chart of the designed NSGA- II in this paper.

7

Computers & Industrial Engineering 137 (2019) 106076

X. Yang, et al.

For the optimal timetable and current timetable, the expected value of total service time of the same number of trains should be consistent, i.e.,

E (aIN

4. Solution algorithm Target conflict is the key problem for multi-objective optimization, that is, there is no such solution, which enables all the objectives to be optimal at the same time. Searching a Pareto non-dominated front with a set of equally good solutions is the basic goal of the multi-objective optimization. A set of equally optimal solutions are called as a Pareto front, which are non-dominated solutions and could reduce the conflict among the various targets to the minimum. The rest ones are dominated solutions that dominated by optimal solutions in Pareto front. One of the goals of the ideal multi-objective optimization is to find the Pareto optimal solutions as many as possible. Srinivas and Deb (1994) proposed NSGA to obtain efficient Pareto optimal solution set. But the computation is complex, the parents are easy to be covered and it needs to make share parameters. Deb, Pratap, Agarwal, and Meyarivan (2002) proposed NSGA-II that effectively overcomes the defects of NSGA. A detail description of the NSGA-II in this paper is given as follows and the flow chart is shown in Fig. 6. Some other multi-objective algorithms (e.g., NSGA-III, SPEA2, PAES, PESA-II, OMOPSO, etc.) can be also used to solve the developed model (Deb & Jain, 2014; Durillo & Nebro, 2011; Durillo, Nebro, & Alba, 2010).

(22)

a11) = T3.

It is seen that all the above constraints are linear. 3.5. Optimization model In the metro timetable optimization problem with dwell time uncertainty, we want to minimize the expected values of energy consumption and passenger waiting time, as well as the mean square errors of them. Therefore, based on the above analysis of objective functions and constraints, the timetable optimization problem is formulated as the following multi-objective integer programming model:

min {JE (a, d, ), JP (a, d, ), JR (a, d, )} s . t . ain , din , t1(n, n + 1)

ai (n + 1)

for all 1 i I , 1 n N for all 1 i I , 1 n N 1 din

for all 1

i

I, 1

n

N

1

(1) Population initialization

t2(n, n + 1), din ain = n, hl < hi < h u , T1 E (aiN ai1) T2, E (aIN a11) = T3.

A solution x = {a11, a12 , ...,a1N , d11, d12, ...,d1N 1, h1, h2 , ...,hI 1} , which is composed of decision variables ain, din and hi, is denoted by a chromosome X. Initialize the population according to the constraints in equation (23).

for all 1 i I , 1 n N 1 for all 1 i I 1, 1 n N 1 for all 1 i I (23)

(2) Non-dominated sorting

3.6. Convergence property

According to non-domination, the population is sorted into several fronts. The detailed steps of fast sort algorithm are as follows.

For proving the convergence property of a multi-objective programming model, the sufficient condition is to prove each objective in the model is convergent to the decision variables.

For each p P ● Initialize Sp =

Proposition 1. The model (23) is convergent to the decision variables, i.e., arrival time a, departure time d and headway time h

*if JE , p JE, q &JP ,p < JP, q &JR, p < JR,q or JE , p < JE, q &JP, p JE , p < JE, q &JP, p < JP, q &JR, p JR,q , then

Proof. In essence, this proof aims to show the net energy consumption function JE, the passenger waiting time JP, and the robustness function JR are all convergent to the decision variables. Arrival time and departure time can be represented by running time and dwell time.

• For the energy consumption function J

•

•

and calculate the fitness values JE , p , JP ,p and JR, p for p.

● Initialize np = 0. ● For each q P , calculate the fitness values JE , q , JP ,q and JR, q .

● Sp = Sp

JP, q &JR, p < JR,q or

{q} , add q to the set of solutions dominated by p.

*else ● np = np + 1, increase the domination of counter p. ● if np = 0 , then

* prank = 1, F1 = F1 {p} , p belongs to the first front. Initialize the front counter, i = 1. While Fi ● Q = , Q is used to store the members of the next front. ● For each p Fi *For each q Sp

= JF − JB, the tractive energy consumption JF is negative correlation with the running time on each section. Constraints (18) ensure that the running time has both upper and lower bounds, such that the limitation value of JF is determinate. The upper bound of regenerative energy is the total braking energy, and the lower bound of regenerative energy is zero. The limitation value of JB is determinate. Therefore, the energy consumption JE is convergent. For the passenger waiting time JP, the best condition is that every passenger arrives to the station with the same time of a train arriving, and the waiting time is zero. The worst condition is that every passenger arrives to the station when the train departs from the station, and the waiting time is a headway time. Therefore, the both upper and lower bounds of JP is determinate, i.e., the passenger waiting time JP is convergent. For the robustness function JR, it is defined with the sum of mean square errors of total net energy consumption and total passenger waiting time. The convergence of JR focuses on the stochastic dwell time. It is all known that the mean square error is a convergent statistical indicator. E

● nq = n q 1 ● if np = 0 , then

qrank = i + 1, Q = Q {q} , q belongs to the next front. Update the counter of the next front, i = i + 1. Store the next front, i.e. Fi = Q .

(3) Crowding distance Crowding distance represents the density of individuals surrounding a particular individual in the population. The calculation of crowding distance ensures the diversity of the population. For the individuals in the same front, the comparison among crowding distance values is meaningful. The calculation for crowing distance of l individuals in front Fi is as follows. ● For each individual x j in the population, set the crowding distance d (x j ) = 0 . ● For each objective, assign infinite distance to boundary individuals in Fi i.e.

d (x 1) = and d (x l) = ● For j = 2 to l − 1

8

.

Computers & Industrial Engineering 137 (2019) 106076

X. Yang, et al. *d (x j ) =

|JE , j + 1 max JE

JE , j 1 | min JE

+

|JP , j + 1 JPmax

JP , j 1 | JPmin

JE , j + 1 represents the value of individual value of individual x j

1

|JR, j + 1 JR, j 1 | max J min JR R x j + 1 in objective

*Generate a random number *Calculate the value of s

+

JE , and JE , j

1

is the

in objective JE .

s

=

(2 s ) 1

(4) Selection

1 m +1

s

1,

[2(1

s )]

1 m+1 ,

[0, 1)

s

< 0.5,

s

0.5,

where m is a non-negative real number. *The offspring in next population i + 1 is calculated below.

After the operation of non-dominated sorting and the calculation of crowding distance, each solution in the population has two attributes. Select the individuals with the minimal rank into the next offspring population. For the individuals with the same rank, the individual with the maximal crowding distance is better. The detail process of selection is as follows:

xs1, i + 1 = x s1, i +

s.

(6) Recombination and selection Combine the current generation population and the offspring population into a new population. Sort the new population based on nondomination and calculate the crowding distance of each individual. Fill all individuals of each front into the next generation population in the ascending order until the population size exceeds the current population size pop_size. When adding all individuals in front Fi into the next generation population, the population exceeds pop_size. For all individuals in front Fi , compare the crowding distance and select the ones with smaller crowding distance values until the size of the population reaches pop_size (Deb et al., 2002). When adding all individuals in front Fi into the next generation population, the population exceeds P.

● For each individual p P and q P , *if prank < qrank Select individual p into the next offspring population. *else if prank > qrank Select individual q into the next offspring population. *else if prank = qrank

● if d (p) > d (q) Select individual p into the next offspring population. ● else Select individual q into the next offspring population.

(5) Simulated Binary Crossover and Polynomial Mutation

(7) Stop criterion

In this paper, the individuals in the population are encoded as real numbers.

The stop criterion determines when the system has reached the predetermined maximum number of iterations max_gen or predetermined precision.

(1) Simulated binary crossover is proposed by Ded and his students in 1995. The two parents in population i are encoded as 1 1 1 1 x 1, i = {a11 , a12 , ...,a11N , d11 , d12 , ...,d11N 1, h11, h 21, ...,hI1 1} and x 2, i = 2 2 2 2 {a11 , a12 , ...,a12N , d11 , d12 , ...,d12N 1, h12, h22, ...,hI2 1} . The crossover probability is Pc and the process of generating offspring is given below.

5. Case study This section presents some numerical results using the real-world operation and passenger data from the Beijing Metro Yizhuang Line to illustrate the efficiency and effectiveness of the developed multi-objective programming model and solution algorithm.

x i , which satisfies the crossover probability, xsi is the sth

● For each variable xsi

variable of individual x i in population i *Generate a random number µs [0, 1) *Calculate the value of

s

=

(2µs )

1 c +1 ,

µs 1 c +1

1 2(1

µs )

5.1. Basic data

s

Beijing Metro Yizhuang Line consists of 14 stations from Songjiazhuang station to Huochezhan station and the length of each section is shown in Fig. 7. The current real-world operation and passenger data of the Beijing Metro Yizhuang Line are provided by the Beijing Metro Operation Company. The original timetable is presented in Table 3. The passenger arrival rate for a whole day at all 14 stations is shown in Fig. 8. The passenger arrival rate is used the averagely daily value of the smart card data during Dec. 2016. Note that the Huochezhan station is not opened to passengers yet, the passenger arrival rate at this station is zero. The dwell time of each train at each station is a random variable. And we call the operation of all trains on the whole line as a scenario. This paper sets 50 scenarios and calculates the expected values of the three objectives. According to the current timetable, we set a low limit and a high limit for each random variable. For one scenario, Fig. 9

0.5,

, µs > 0.5,

where c is a non-negative real number. *The offspring in next population i + 1 is calculated as

xs1,i + 1 = 0.5 (1 + xs2, i + 1 = 0.5 (1

1, i s ) xs 1, i ) x s s

+ (1 + (1 +

2, i s ) xs 2, i ) x s s

, .

The symbol xs1, i represents the sth variable in parent x 1,i . 1 1 1 1 , a12 , ...,a11N , d11 , d12 , ...,d11N 1, h11, h21, ...,hI1 1} (2) Use the parent x 1, i ={a11 to make polynomial mutation. The mutation probability is Pm and we can get the offspring below.

● For each xsi

x i , which satisfies the mutation probability

Fig. 7. Illustration of the Beijing Metro Yizhuang Line. 9

Computers & Industrial Engineering 137 (2019) 106076

X. Yang, et al.

Table 3 Original timetable for the Beijing Metro Yizhuang Line. Station

SJZ

XC

XHM

JG

YZQ

WHY

WY

Arrival time (s) Departure time (s) Station Arrival time (s) Departure time (s)

0 30 RJ 1112 1142

220 250 RC 1246 1276

358 388 TJN 1440 1470

545 575 JH 1620 1650

710 745 CQN 1790 1825

835 865 CQ 1927 1972

979 1009 YZ 2077 –

shows ten possible values of dwell time of ten trains at each station, in which different lines represent dwell time of different trains. The value and unit of remaining parameters are listed in Table 4. 5.2. Optimal timetable In this section, the NSGA-II is coded in MATLAB programming language to verify the efficiency and effectiveness of the proposed model. The calculation is based on the data of Beijing Metro Yizhuang Line at early peak hours, which takes 18.83 min to solve the proposed model. The feasible solutions for Beijing Metro Yizhuang Line are shown in Fig. 10. With one of the Pareto-optimal solutions, one of the best found timetables is shown in Table 5. Fig. 11 shows the comparison between the current timetable and the optimal timetable. We can see that the arrival time, departure time and headway time of the optimal timetable are different from that in the current timetable, which clearly shows the difference before and after the optimization.

Fig. 9. Possible values of dwell time. Table 4 Value and unit of some parameters.

5.3. System performance in Pareto-optimal solutions

Parameter

N

I

m

fa

fb

r

g

T

δ

Value Unit Parameter Value

14 – ƞ1 0.7

20 – ƞ2 0.3

311,800 kg pop_size 200

315,000 N max_gen 200

258,000 N Pc 0.9

2000 N Pm 0.5

500 N ƞc 20

2077 s ƞm 20

5 s

the performance of scheduling system. Therefore, the total service time and number of trains should not change too much. In this paper, the number of trains is fixed. We randomly select a robust solution from the Pareto-optimal solutions and compare the service time with that of the current timetable. The service time of both current and optimal timetable, from the first train arrives at station 1 to the last train arrives at station N, is shown in Table 6. As shown in Table 6, there are no significant changes on service

In order to validate the effectiveness of the proposed model and algorithm, we propose the following five indicators. (1) Service time The service time is closely related to passenger satisfaction, operation cost, energy consumption and so on, which is important to evaluate

Fig. 8. Passenger arrival rate for a whole day at stations of the Beijing Metro Yizhuang Line (average data during Dec. 2016). 10

Computers & Industrial Engineering 137 (2019) 106076

X. Yang, et al.

Table 6 Comparison of service time between current and optimal timetable. Train number

1

1

2

3

4

5

6

Current timetable Optimal timetable Train number Current timetable Optimal timetable Train number Current timetable Optimal timetable

0 0 7 2857 2863 14 3767 3734

2077 2095 8 2987 2966 15 3897 3849

2207 2228 9 3117 3112 16 4027 3996

2337 2327 10 3247 3217 17 4157 4156

2467 2428 11 3377 3342 18 4287 4330

2597 2550 12 3507 3440 19 4417 4453

2727 2715 13 3637 3579 20 4547 4547

Table 7 Comparison between the current timetable and optimal timetable.

Fig. 10. Pareto-optimal solutions with Beijing Metro Yizhuang Line.

Solution number

Energy consumption

Passenger waiting time

Robustness

Optimal timetable

1 2 3 4 5 6 7 8 9 10 Average

3581.95 3581.92 3583.43 3581.02 3614.24 3624.01 3624.26 3658.90 3660.53 3693.87 3620.41

594.70 592.79 592.71 593.51 592.23 591.55 591.49 591.85 593.21 595.72 592.98

51.80 51.71 50.18 61.57 52.98 55.25 54.80 66.79 50.05 50.06 54.52

Current timetable

–

3698.10

704.25

72.51

Reduction

–

2.10%

15.80%

24.81%

Table 5 Optimal timetable for the Beijing Metro Yizhuang Line. Station

SJZ

XC

XHM

JG

YZQ

WHY

WY

Arrival time (s) Departure time (s) Station Arrival time (s) Departure time (s)

0 31 RJ 1117 1147

225 255 RC 1248 1278

362 392 TJN 1442 1472

549 579 JH 1628 1658

711 746 CQN 1796 1831

837 867 CQ 1935 1980

986 1016 YZ 2095 –

time, suggesting that the optimal timetable can make reductions/improvement in energy consumption, passenger waiting time and robustness while the change on service time is very slight.

better in consideration on the uncertain dwell time. In fact, different optimal timetables can lead to different expected values of energy consumption, passenger waiting time and robustness. The smaller value of robustness is with the stronger adaptability of the timetable. And for the operator preferences to energy consumption and passenger waiting time, this paper provides different Pareto-robust timetable. In order to further prove the validity of the proposed model, we propose three indicators to measure the discrete degree of energy consumption and passenger waiting time.

(2) Expected values of objectives According to fifty scenarios of dwell time and the current running time, the expected values of energy consumption, passenger waiting time and robustness are calculated. In comparison with the current timetable, several Pareto-optimal solutions are selected to verify the effectiveness of the proposed model in Table 7. As shown in Table 7, the average value of energy consumption of optimal timetable is 3620.41 kW·h, which reduces 2.10% compared with 3698.10 kW·h of the current timetable. Passenger waiting time drops from 704.25 to 592.98 persons each hour, which means passenger waiting time decreases 15.80%. Robustness value reduces 24.81%, that is, the stability and adaptability of the optimal timetable is

(3) Variance indicators Variance indicators denote the discrete degree of the two main objectives (i.e., energy consumption and passenger waiting time) under the fifty scenarios of dwell time, which are described as follows:

Fig. 11. Current timetable and optimal timetable during morning peak hours. 11

Computers & Industrial Engineering 137 (2019) 106076

X. Yang, et al.

Table 8 Values of the three indicators. Solution number

VIE

VIP

MADE

MADP

RE

RP

Optimal timetable

1 2 3 4 5 6 7 8 9 10 Average

31.759 32.323 30.627 41.884 33.112 37.250 34.875 39.054 29.909 29.326 34.012

20.045 19.391 19.554 19.687 19.868 29.162 19.930 27.733 20.143 20.739 21.625

12.547 12.676 11.944 14.033 12.998 14.033 13.957 14.700 10.722 12.967 13.058

11.733 11.345 12.802 12.760 11.480 11.408 12.760 11.690 11.585 11.886 11.945

22.840 23.530 23.700 26.130 24.110 26.240 26.070 26.640 20.080 22.250 24.159

21.686 20.948 22.405 23.274 21.269 23.368 23.307 20.657 23.075 22.927 22.292

Current timetable

–

43.998

28.517

13.407

13.526

26.710

23.677

Reduction

–

22.70%

24.17%

2.61%

11.69%

9.55%

5.85%

VI E = V [JE (a, d, h, )], VI P = V [JP (a, d, h, )].

developed approach can improve the performance of total energy consumption, total passenger waiting time and robustness value by 2.10%, 15.80% and 24.81%. We consider the problem on a single line without extending it on a network. Also the energy consumption is calculated by kinematic equations but not an electrical energy flow model. In addition, the passenger behavior during loading and unloading is not considered when calculating the passenger waiting time. These limitations can be deeply studied in our future research. Timetable/schedule optimization in other transportation systems (e.g., air traffic, maritime traffic, etc.) is also usually considered as a multi-objective decision problem. The developed multi-objective integer programming model and solution algorithm can be applied to these areas, where the formulations of objectives and constraints should be modified according to the specific problem.

(24)

(4) Maximum absolute differences If the expected values of the objectives are the center of a Paretooptimal solution, the maximum absolute differences reflect the farthest distance from the center of the solution. Maximum absolute differences are described as follows:

MAD E = max{Abe |Abe = |JE (a, d, h, )

E (JE (a, d, h, ))|},

MAD P = max{Abp |Abp = |JP (a, d , h, )

E (JP (a, d, h, ))|}.

(25)

(5) Ranges

Acknowledgements

In a Pareto-optimal solution, range means the distance between the maximum value and the minimum value of the objectives, indicating the discrete degree of the objectives:

RE = max{JE (a, d, h, )}

min{JE (a, d , h, )},

RP = max{JP (a, d, h, )}

min{JP (a, d , h, )}.

This work was supported by the National Natural Science Foundation of China (Nos. 71701013, 71890972/71890970, 71525002, 71621001), the Beijing Municipal Natural Science Foundation (No. L181008), the Young Elite Scientists Sponsorship Program by CAST (No. 2018QNRC001), and the State Key Laboratory of Rail Traffic Control and Safety (No. RCS2019ZZ001).

(26)

For energy consumption, the reduction of variance indicator is 22.70% in Table 8, from 43.998 to 34.012. Maximum absolute difference is 13.058, reduced by 2.61% than 13.407 of current timetable. Range of optimal timetable is 24.159, which is 9.55% less than current value 26.710. For passenger waiting time, the optimal variance is 21.625, which is 24.17% less than current variance. The average farthest distance from the center of the Pareto-optimal solution is 11.945, which is 11.10% shorter than that of the current timetable. The difference between the maximum value and minimum value decreases by 5.85%, from 23.677 to 22.292. The results of the three indicators, variance indicators, maximum absolute differences and ranges, show that the discrete degree of the optimal timetable is lower, which means the optimal timetable is more stable and reliable.

References Barrena, E., Canca, D., Coelho, L. C., & Laporte, G. (2014b). Single-line rail rapid transit timetabling under dynamic passenger demand. Transportation Research Part B: Methodological, 70, 134–150. Barrena, E., Canca, D., Coelho, L. C., & Laporte, G. (2014a). Exact formulations and algorithm for the train timetabling problem with dynamic demand. Computers & Operations Research, 44, 66–74. Cacchiani, V., Caprara, A., & Fischetti, M. (2012a). A Lagrangian heuristic for robustness, with an application to train timetabling. Transportation Science, 46(1), 124–133. Cacchiani, V., Caprara, A., Galli, L., Kroon, L., Maróti, G., & Toth, P. (2012b). Railway rolling stock planning: Robustness against large disruptions. Transportation Science, 46(2), 217–232. Cacchiani, V., & Toth, P. (2012). Nominal and robust train timetabling problems. European Journal of Operational Research, 219(3), 727–737. Canca, D., & Zarzo, A. (2017). Design of energy-efficient timetables in two-way railway rapid transit lines. Transportation Research Part B: Methodological, 102, 142–161. Cao, Y., Ma, L., Xiao, S., Zhang, X., & Xu, W. (2017). Standard analysis for transfer delay in CTCS-3. Chinese Journal of Electronics, 26(5), 1057–1063. Cao, Y., Wang, Z., Liu, F., Li, P., & Xie, G. (2019). Bio-inspired speed curve optimization and sliding mode tracking control for subway trains. IEEE Transactions on Vehicular Technology, 68(7), 6331–6342. ChinaDaily (2017). Urban rail growth on fast track. [Online]. Available: http://www. chinadaily.com.cn/a/201712/28/WS5a442b96a31008cf16da3ddb.html. Corman, F., D’Ariano, A., Marra, A. D., Pacciarelli, D., & Samà, M. (2017). Integrating train scheduling and delay management in real-time railway traffic control. Transportation Research Part E: Logistics and Transportation Review, 105, 213–239. Corman, F., D’Ariano, A., Pacciarelli, D., & Pranzo, M. (2012). Bi-objective conflict detection and resolution in railway traffic management. Transportation Research Part C: Emerging Technologies, 20(1), 79–94.

6. Conclusion The main contribution of this paper is to use the concept of robustness to tackle the uncertain dwell time in real-world train operations for improving the performance of timetable optimization approaches. Based on this, we develop a multi-objective integer programming model to solve the metro timetable optimization problem with dwell time uncertainty, where the energy consumption, passenger waiting time, as well as robustness are minimized of compromise. The case study based on the Beijing Metro Yizhuang line shows that the 12

Computers & Industrial Engineering 137 (2019) 106076

X. Yang, et al.

train control and timetabling. European Journal of Operational Research, 257(2), 355–376. Sels, P., Dewilde, T., Cattrysse, D., & Vansteenwegen, P. (2016). Reducing the passenger travel time in practice by the automated construction of a robust railway timetable. Transportation Research Part B: Methodological, 84, 124–156. Srinivas, N., & Deb, K. (1994). Multi-objective optimization using non-dominated sorting in genetic algorithms. Evolutionary Computation, 2(3), 221–248. Tian, X., & Niu, H. (2019). A bi-objective model with sequential search algorithm for optimizing network-wide train timetables. Computers & Industrial Engineering, 127, 1259–1272. Vansteenwegen, P., & Van Oudheusden, D. (2007). Decreasing the passenger waiting time for an intercity rail network. Transportation Research Part B: Methodological, 41(4), 478–492. Wikipedia (2018). Urban rail transit. [Online]. Available: http://en.wikipedia.org/wiki/ Urban_rail_transit. Xu, X., Li, K., & Li, X. (2016). A multi-objective subway timetable optimization approach with minimum passenger time and energy consumption. Journal of Advanced Transportation, 50(1), 69–95. Yang, X., Chen, A., Li, X., Ning, B., & Tang, T. (2015a). An energy-efficient scheduling approach to improve the utilization of regenerative energy for metro systems. Transportation Research Part C: Emerging Technologies, 57, 13–29. Yang, X., Chen, A., Ning, B., & Tang, T. (2016a). A stochastic model for the integrated optimization on metro timetable and speed profile with uncertain train mass. Transportation Research Part B: Methodological, 91, 424–445. Yang, X., Chen, A., Wu, J., Gao, Z., & Tang, T. (2019). An energy-efficient rescheduling approach under delay perturbations for metro systems. Transportmetrica B: Transport Dynamics, 7(1), 386–400. Yang, X., Li, X., Gao, Z., Wang, H., & Tang, T. (2013). A cooperative scheduling model for timetable optimization in subway systems. IEEE Transactions on Intelligent Transportation Systems, 14(1), 438–447. Yang, X., Li, X., Ning, B., & Tang, T. (2015b). An optimisation method for train scheduling with minimum energy consumption and travel time in metro rail systems. Transportmetrica B: Transport Dynamics, 3(2), 79–98. Yang, X., Li, X., Ning, B., & Tang, T. (2016b). A survey on energy-efficient train operation for urban rail transit. IEEE Transactions on Intelligent Transportation Systems, 17(1), 2–13. Yang, X., Ning, B., Li, X., & Tang, T. (2014). A two-objective timetable optimization model in subway systems. IEEE Transactions on Intelligent Transportation Systems, 15(5), 1913–1921. Yang, X., Yin, H., Wu, J., Qu, Y., & Gao, Z. (2019). Recognizing the critical stations in urban rail networks: An analysis method based on the smart-card data. IEEE Intelligent Transportation Systems Magazine, 11(1), 29–35. Ye, H., & Liu, R. (2016). A multiphase optimal control method for multi-train control and scheduling on railway lines. Transportation Research Part B: Methodological, 93, 377–393. Yin, J., Yang, L., Tang, T., Gao, Z., & Ran, B. (2017). Dynamic passenger demand oriented metro train scheduling with energy-efficiency and waiting time minimization: Mixedinteger linear programming approaches. Transportation Research Part B: Methodological, 97, 182–213. Zhang, H., Li, S., & Yang, L. (2018). Real-time optimal train regulation design for metro lines with energy-saving. Computers & Industrial Engineering. https://doi.org/10. 1016/j.cie.2018.02.019. Zhang, L., 2018. The operation data for the Beijing Metro Yizhuang Line. Beijing Mass Transit Railway Operation Corporation Limited Technical report (in Chinese). Zhou, L., Tong, L. C., Chen, J., Tang, J., & Zhou, X. (2017). Joint optimization of highspeed train timetables and speed profiles: A unified modeling approach using spacetime-speed grid networks. Transportation Research Part B: Methodological, 97, 157–181.

Deb, K., & Jain, H. (2014). An evolutionary many-objective optimization algorithm using reference-point-based nondominated sorting approach, part I: Solving problems with box constraints. IEEE Transactions on Evolutionary Computation, 18(4), 577–601. Deb, K., Pratap, A., Agarwal, S., & Meyarivan, T. A. M. T. (2002). A fast and elitist multiobjective genetic algorithm: NSGA-II. IEEE Transactions on Evolutionary Computation, 6(2), 182–197. Domínguez, M., Fernández-Cardador, A., Cucala, A. P., & Pecharromán, R. R. (2012). Energy savings in metropolitan railway substations through regenerative energy recovery and optimal design of ATO speed profiles. IEEE Transactions on Automation Science and Engineering, 9(3), 496–504. Durillo, J. J., & Nebro, A. J. (2011). jMetal: A Java framework for multi-objective optimization. Advances in Engineering Software, 42(10), 760–771. Durillo, J. J., Nebro, A. J., & Alba, E. (2010). The jMetal framework for multi-objective optimization: Design and architecture. IEEE congress on evolutionary computation, CEC 2010 (pp. 4138–4325). . Fischetti, M., Salvagnin, D., & Zanette, A. (2009). Fast approaches to improve the robustness of a railway timetable. Transportation Science, 43(3), 321–335. Ghoseiri, K., Szidarovszky, F., & Asgharpour, M. J. (2004). A multi-objective train scheduling model and solution. Transportation Research Part B: Methodological, 38(10), 927–952. Higgins, A., Kozan, E., & Ferreira, L. (1996). Optimal scheduling of trains on a single line track. Transportation Research Part B: Methodological, 30(2), 147–161. Huang, Y., Yu, H., Yin, J., Hu, H., Bai, S., Meng, X., & Wang, M. (2018). An integrated approach for the energy-efficient driving strategy optimization of multiple trains by considering regenerative braking. Computers & Industrial Engineering, 126, 399–409. Jovanović, P., Kecman, P., Bojović, N., & Mandić, D. (2017). Optimal allocation of buffer times to increase train schedule robustness. European Journal of Operational Research, 256(1), 44–54. Lee, Y., Lu, L. S., Wu, M. L., & Lin, D. Y. (2017). Balance of efficiency and robustness in passenger railway timetables. Transportation Research Part B: Methodological, 97, 142–156. Li, X., & Lo, H. K. (2014a). An energy-efficient scheduling and speed control approach for metro rail operations. Transportation Research Part B: Methodological, 64, 73–89. Li, X., & Lo, H. K. (2014b). Energy minimization in dynamic train scheduling and control for metro rail operations. Transportation Research Part B: Methodological, 70, 269–284. Li, X., Wang, D., Li, K., & Gao, Z. (2013). A green train scheduling model and fuzzy multiobjective optimization algorithm. Applied Mathematical Modelling, 37(4), 2063–2073. Li, X., & Yang, X. (2013). A stochastic timetable optimization model in subway systems. International Journal of Uncertainty, Fuzziness & Knowledge-Based Systems, 21(Supp. 1), 1–15. Lu, C., Tang, J., Zhou, L., Yue, Y., & Huang, Z. (2017). Improving recovery-to-optimality robustness through efficiency-balanced design of timetable structure. Transportation Research Part C: Emerging Technologies, 85, 184–210. Niu, H., & Zhou, X. (2013). Optimizing urban rail timetable under time-dependent demand and oversaturated conditions. Transportation Research Part C: Emerging Technologies, 36, 212–230. Niu, H., Zhou, X., & Gao, R. (2015). Train scheduling for minimizing passenger waiting time with time-dependent demand and skip-stop patterns: Nonlinear integer programming models with linear constraints. Transportation Research Part B: Methodological, 76, 117–135. Ramos, A., Peña-Alcaraz, M., Fernández-Cardador, A., & Cucala, A. P. (2007). Mathematical programming approach to underground timetabling problem for maximizing time synchronization. Proceedings of international conference on industrial engineering and industrial management (pp. 88–95). . Salido, M. A., Barber, F., & Ingolotti, L. (2012). Robustness for a single railway line: Analytical and simulation methods. Expert Systems with Applications, 39(18), 13305–13327. Scheepmaker, G. M., Goverde, R. M., & Kroon, L. G. (2017). Review of energy-efficient

13