A stochastic model for the integrated optimization on metro timetable and speed profile with uncertain train mass

Transportation Research Part B 91 (2016) 424–445

Contents lists available at ScienceDirect

Transportation Research Part B journal homepage: www.elsevier.com/locate/trb

A stochastic model for the integrated optimization on metro timetable and speed profile with uncertain train mass Xin Yang a, Anthony Chen b,c,∗, Bin Ning a, Tao Tang a a

State Key Laboratory of Rail Traffic Control and Safety, Beijing Jiaotong University, Beijing 100044, China Department of Civil and Environmental Engineering, The Hong Kong Polytechnic University, Kowloon, Hong Kong c Key Laboratory of Road and Traffic Engineering, Tongji University, Shanghai 201804, China b

a r t i c l e

i n f o

Article history: Received 2 June 2015 Revised 14 June 2016 Accepted 15 June 2016 Available online 22 June 2016 Keywords: Integrated optimization Timetable Speed profile Uncertain train mass Metro systems

a b s t r a c t The integrated timetable and speed profile optimization model has recently attracted more attention because of its good achievements on energy conservation in metro systems. However, most previous studies often ignore the spatial and temporal uncertainties of train mass, and the variabilities of tractive force, braking force and basic running resistance on energy consumption in order to simplify the model formulation and solution algorithm. In this paper, we develop an integrated metro timetable and speed profile optimization model to minimize the total tractive energy consumption, where these real-world operating conditions are explicitly considered in the model formulation and solution algorithm. Firstly, we formulate a two-phase stochastic programming model to determine the timetable and speed profile. Given the speed profile, the first phase determines the timetable by scheduling the arrival and departure times for each station, and the second phase determines the speed profile for each inter-station with the scheduled arrival and departure times. Secondly, we design a simulation-based genetic algorithm procedure incorporated with the optimal train control algorithm to find the optimal solution. Finally, we present a simple example and a real-world example based on the operation data from the Beijing Metro Yizhuang Line in Beijing, China. The results of the real-world example show that, during peak hours, off-peak hours and night hours, the total tractive energy consumptions can be reduced by: (1) 10.66%, 9.94% and 9.13% in comparison with the current timetable and speed profile; and (2) 3.35%, 3.12% and 3.04% in comparison with the deterministic model. © 2016 Elsevier Ltd. All rights reserved.

1. Introduction Metro system has received rapid development around the world because of its high reliability and large transport capacity. For example, only 3 cities in mainland China (i.e., Beijing, Shanghai and Guangzhou) have metro lines in 20 0 0. By the end of 2014, the number of cities has increased to 22 with 88 metro lines and a total distance of more than 30 0 0 km (Online News, 2015). Along with this rapid expansion is the sharp increase in total energy consumption in metro systems. Almost half of the energy is used by the tractive system for accelerating trains. Take the Beijing Metro system as an example, the total energy consumption and the tractive energy consumption of some main metro lines in 2012 are shown in Fig. 1. ∗

Corresponding author. E-mail addresses: [email protected] (X. Yang), [email protected], [email protected] (A. Chen), [email protected] (B. Ning), [email protected] (T. Tang). http://dx.doi.org/10.1016/j.trb.2016.06.006 0191-2615/© 2016 Elsevier Ltd. All rights reserved.

X. Yang et al. / Transportation Research Part B 91 (2016) 424–445

15

x 10

425

7

Energy consumption (kW.h)

Annual total energy consumption Annual tractive energy consumption

10 49.93%

58.81% 74.64% 60.29%

5

53.98% 56.27% 56.04%

49.90%

0

1 Line

2 Line

5 Line

8 Line

10 Line

13 Line

Line

43.91%

15 Line Line shan Yizhuang Fang

Fig. 1. Energy consumption of some main metro lines in Beijing in 2012 (Source: Zhang (2013)).

Due to the rising environmental concerns and energy prices, energy conservation has attracted much attention in recent years. Energy-efficient train operations can make a significant contribution to reduce the tractive energy consumption. However, traditional studies on energy-efficient train operations, including timetable optimization (Ghoseiri et al., 2004; Yang et al., 2014), speed profile optimization (Ke et al., 2009; 2012) and integrated optimization (Li and Lo, 2014a; 2014b), usually assume the train mass, tractive force, braking force and basic running resistance as constant parameters. In realworld metro systems, the tractive force, braking force and basic running resistance are variable according to the train speed variation. Train mass consists of the empty vehicle mass and the variable mass of passengers on the train. The former is considered fix (or deterministic), while the latter is highly variable (or stochastic) at different inter-stations and at different operational periods (e.g., peak hours, off-peak hours and night hours). The purpose of this paper is to develop a two-phase stochastic programming model for the integrated timetable and speed profile optimization problem with the aim of minimizing the total tractive energy consumption. Compared to previous studies, this paper has the following contributions and differences. (1) Compared to Yang et al. (2013) and Yang et al. (2015a), which only optimize the timetable to reduce the energy consumption, this paper develops a two-phase stochastic programming model to optimize both timetable and speed profile for achieving better energy performance. (2) Compared to Li and Lo (2014a) and Li and Lo (2014b), which assume deterministic train mass and fixed tractive force, braking force and basic running resistance as constant parameters for simplifying the model formulation and solution algorithm, this paper explicitly considers the uncertainty of train mass and the variability of tractive force, braking force and basic running resistance in the integrated timetable and speed profile optimization problem. In addition, this paper analyzes the convergence property and formulation complexity of the developed model. (3) Compared to Yang et al. (2015b), which also assume the train mass, tractive force and braking force as constant parameters, and a predetermined maximum train speed in the optimal train control algorithm (Howlett and Pudney, 1995), this paper designs a more general optimal train control algorithm that can generate an energy-efficient speed profile without presupposing a maximum train speed. Furthermore, this paper provides an existence proof of an optimal speed profile for a given inter-station with a reasonable running time. The remainder of this paper is organized as follows. In Section 2, we review the literature on the energy-efficient train operations. In Section 3, we describe the problem statement and motivation of this study based on real-world metro systems, such as the Beijing Metro System. In Section 4, we develop a two-phase stochastic programming model to determine the optimal timetable and speed profile. In Section 5, we design a simulation-based genetic algorithm procedure incorporated with the optimal train control algorithm to find the optimal solution. In Section 6, we present a simple example and a realworld example based on the operation data from the Beijing Metro Yizhuang Line. Conclusions are provided in Section 7. 2. Literature review In metro systems, energy-efficient train operations have gained attention in recent decades because of the good achievements in reducing energy consumption. The literature on energy-efficient train operations include timetable optimization, speed profile optimization, and integrated timetable and speed profile optimization.

426

X. Yang et al. / Transportation Research Part B 91 (2016) 424–445

The train timetable optimization problem is one of the widely studied problems in rail transportation systems (Bešinovic´ et al., 2016; Cacchiani et al., 2010; Meng and Zhou, 2014; Niu et al., 2015; Xu et al., 2015; Zhou and Zhong, 2007). Many recent studies focused on improving the utilization of regenerative braking energy by synchronizing accelerating and braking trains to reduce energy consumption. For example, Ramos et al. (2007) proposed a timetable optimization model to maximize the overlapping time between adjacent trains so that the utilization of regenerative braking energy can be improved. Nasri et al. (2010) also developed a timetable optimization model by making use of the regenerative energy from braking trains to study the effect of headways and reserve times on the amount of energy consumption. Peña Alcaraz et al. (2012) extended the timetable optimization model to synchronize the braking and acceleration processes of all trains in the same electrical section for reducing energy consumption. They also developed a power flow model to calculate the regeneration saving factor of each synchronization within the electrical section. Yang et al. (2013) proposed a cooperative scheduling model to coordinate the accelerating and braking processes of adjacent trains so that the regenerative energy from the braking trains can be immediately used by the accelerating trains. Furthermore, Li and Yang (2013) extended the timetable optimization problem to consider the randomness of departure delay for trains at busy stations. Yang et al. (2015a) developed an energy-efficient scheduling model based on the real-world speed profile from the Beijing Metro Yizhuang Line to improve the utilization regenerative braking energy, where the cycle time and the number of trains were kept unchanged in the metro operations. Research on the speed profile optimization problem began as early as 1960s. Ishikawa (1968) firstly proposed an optimal control model to determine the speed profile. Under the optimal control approach, research studies are mainly grouped into two categories: exact solution algorithm and heuristic algorithm (Yang et al., 2016). For the exact solution algorithm, Howlett (1988) proved the problem can be formulated in an appropriate functional space to determine an optimal speed profile according to the Pontryagin type criterion. Howlett (1990) further formulated the speed profile optimization problem as a finite dimensional constrained optimization model, and used the Pontryagin principle to find the optimal strategy. To apply the theory into practice, speed limits, gradients, traction efficiency and track curvature were considered (Khmelnitsky, 20 0 0), which made the speed profile optimization problem more difficult. Howlett et al. (2009) provided an analytical solution method to prove that the cruising phase must be interrupted by the accelerating phase for steep uphill sections and the coasting phase for steep downhill sections. For the heuristic algorithm, Ke et al. (2012) presented a new method based on the MAX-MIN ant system algorithm to optimize train speed profiles for online optimization purposes. Lu et al. (2013) applied the ant colony optimization algorithm, genetic algorithm and dynamic programming to find the optimal speed profile and compared the performance of the three algorithms. Fernández-Rodríguez et al. (2015) developed a robust and energyefficient optimization method with considerations of train load variations, tractive and braking forces variations, and delays uncertainty. However, the model only optimized the speed profiles of inter-stations without accounting for the timetables at stations for a more global consideration of the energy consumption of the whole metro line. Due to the good achievements on energy conservation, the integrated timetable and speed profile optimization problem gained more attention recently. For example, Bocharnikov et al. (2010) developed a multi-train optimization model to consider both tractive energy consumption and utilization of regenerative energy. Ding et al. (2011) formulated a two-level optimization model to determine the running time and speed profile for each inter-station and designed a genetic algorithm to find a good solution. Su et al. (2013) considered the timetable and speed profile optimization problem for one train within one cycle operation to minimize the tractive energy consumption. Tuyttens et al. (2013) designed a genetic algorithm to determine the speed profiles for all trains to reduce the energy consumption by improving the utilization of regenerative energy. Li and Lo (2014a) proposed an integer programming model to determine the timetable and speed profile with the minimum net energy consumption based on the difference between the tractive energy consumption and the utilization of regenerative energy. Later, Li and Lo (2014b) developed a dynamic integrated optimization model with adaptive cycle time based on the passenger demand. The results showed that the proposed dynamic approach can reduce the net energy consumption by 7.86% compared to the static integrated optimization model (Li and Lo, 2014a). Gong et al. (2014) developed an integrated energy-efficient operation method to determine both the timetable and speed profile for the Shanghai Metro Line One. Yang et al. (2015b) developed an integrated optimization method with dwell time and maximum speed control to reduce the energy consumption and travel time for the Beijing Metro Yizhuang Line. Based on the review above, a summary of the integrated timetable and speed profile optimization problem is provided in Table 1. As can be seen, most of them ignored the uncertainty of train mass, and assumed the tractive force, the braking force and the basic running resistance as constants. 3. Problem statement In some real-world metro systems (e.g., the Beijing Metro system), metro managers divide the daily operation into three periods: peak hours, off-peak hours, and night hours based on the daily passenger demand profile. Trains operated during each period are typically assigned with a period-specific timetable and speed profile. However, the three assigned timetables and speed profiles are closely related. To facilitate operations and avoid faults during switching, dwell time at the same station, running time at the same inter-station, and turnaround time are all kept the same during the three operational periods. Further, the speed profile used by all trains at the same inter-station is also kept the same. Only the headways during the three operational periods are altered to accommodate the different passenger demand profiles. That is to say, the operations of metro trains are considered as cyclic and periodic (i.e., similar operational considerations are also adopted

X. Yang et al. / Transportation Research Part B 91 (2016) 424–445

427

Table 1 Recent publications on the integrated timetable and speed profile optimization problem and our work. Publication

Train mass

Tractive/braking force

Basic running resistance

Solution algorithm

Bocharnikov et al. (2010) Ding et al. (2011) Su et al. (2013) Tuyttens et al. (2013) Li and Lo (2014a) Li and Lo (2014b) Gong et al. (2014) Yang et al. (2015b) This paper

Constant Constant Constant Constant Constant Constant Constant Constant Uncertain

Constant Constant Constant Constant Constant Constant Constant Constant Variable

Constant Constant Constant Zero Constant Zero Variable Variable Variable

Genetic algorithm Genetic algorithm Iterative algorithm Genetic algorithm Genetic algorithm Lagrange multiplier Genetic algorithm Adaptive genetic algorithm Simulation-based genetic algorithm

Dwell time at station 2

Keep the same Keep the same Dwell time Dwell time at station 2 at station 2

Running time at inter-station (2, 3) Turnaround time

Running time at inter-station (2, 3)

Running time at inter-station (2, 3)

Turnaround time

Turnaround time

Station 3

Station 4

Station 2

Station 5

Headway

Increase

Headway

Increase

Headway

Station 1

Station 6 Timetable at peak hours Timetable at off-peak hours

Timetable at night hours

Time

Fig. 2. Illustration of timetable with fixed dwell time, fixed running time, fixed turnaround time, and variable headways during peak hours, off-peak hours and night hours.

in Li and Lo (2014a, 2014b); Niu and Zhou (2013). We only need to generate a single set of dwell time, running time, and speed profile and three different headways to determine three sets of timetable and speed profile for peak hours, off-peak hours, and night hours, respectively. For illustration, we provide Fig. 2 to illustrate the timetable with fixed dwell time, fixed running time, fixed turnaround time, and variable headways during the three operational periods. Remark 1. The above paragraph explains that the timetable in metro systems is periodic, which is different from the nonperiodic timetable in general railway systems. The periodicity of metro timetable at different operational periods can be observed in Fig. 3. From the view point of one train operating at different periods, the total travel time of one cycle trip (i.e., departing from station 1 arriving at station 6) is the same, but the headway is different (i.e., using a higher frequency and a correspondingly lesser headway in the peak hours to accommodate the higher passenger demand, and a lower frequency and greater headway in the non-peak hours and night hours). From the timetable perspective, one needs to add some additional trains to the line from off-peak hours to peak hours, and remove some trains from peak hours to night hours. This explains how the timetable and speed profile (to be optimized in Section 4) are used in different operational periods by adding or removing trains according to the variable headways to accommodate the different passenger demand patterns. As train mass is uncertain at different inter-stations and at different operational periods, this paper develops a stochastic optimization model considering the train mass at each inter-station as a stochastic variable based on distinct probabilities for three operational periods to determine the optimal timetable and speed profile. Note that this is not the same as three different deterministic optimization models for the three operational periods. Otherwise, we cannot ensure using the same dwell time, running time, and speed profile for the three operational periods with different headways to accommodate different passenger demand patterns.

428

X. Yang et al. / Transportation Research Part B 91 (2016) 424–445

Station 3

Station 4

Station 2

Station 5

Station 6

Station 1 The 1st cycle of The 2nd cycle of The 3rd cycle of train i train i train i Peak hours Night hours Off-peak hours Add some additional Remove some trains trains to the line from the line

Time

Fig. 3. Periodicity of metro timetable at different operational periods.

Fig. 4. Illustration of a metro line.

In the integrated model, the timetable and speed profile of one train completing one cycle trip are optimized to minimize the total tractive energy consumption. Once we obtain the optimal timetable and speed profile of one train completing one cycle trip, we can use the obtained dwell time, running time and speed profile to determine the optimal timetable and speed profile of all trains throughout the overall operational period with assigned headways. 4. Model formulation This section formulates the stochastic optimization problem as a two-phase expected value programming model to determine the optimal metro timetable and speed profile. The following discussions focus on detailing each part of the model, i.e., notation systems, uncertainty of train mass, system constraints, objective function, stochastic optimization model, convergence property and formulation complexity. 4.1. Notation systems The parameters and variables are listed below for the convenience of formulating the problem under consideration, and some of them are indicated in Fig. 4. Indices and parameters i Train index, i = 1, 2, · · · , I. n Station index, n = 1, 2, · · · , 2N. acin Current arrival time for train i at station n. acin Lower bound of arrival time for train i at station n. acin Upper bound of arrival time for train i at station n. dcin Current departure time for train i at station n. dcin Lower bound of departure time for train i at station n. dcin

Upper bound of departure time for train i at station n.

X. Yang et al. / Transportation Research Part B 91 (2016) 424–445

tcn tcn tt C Xn x v gn (x) Vn (x) F(v) B(v) r(v) m

τ 1n τ 2n τ 3n

429

Lower bound of running time at inter-station (n, n + 1 ). Upper bound of running time at inter-station (n, n + 1 ). Turnaround time at the first station and the terminal station. Total travel time, i.e., the period that a train from departure at station 1 to arrival time at station 2N. Length of inter-station (n, n + 1 ). Length unit. Train speed index. Gradient force per length unit x at inter-station (n, n + 1 ). Speed restriction per length unit x at inter-station (n, n + 1 ). Maximum tractive force. Maximum braking force. Basic running resistance. Empty vehicle mass. Mass of passengers on the train at inter-station (n, n + 1 ) during peak hours. Mass of passengers on the train at inter-station (n, n + 1 ) during off-peak hours. Mass of passengers on the train at inter-station (n, n + 1 ) during night hours.

Stochastic variables ξn Train mass at inter-station (n, n + 1 ). Decision variables din Optimal departure time for train i at station n. ain Optimal arrival time for train i at station n. kF (x, ξ n ) Output rate of tractive force per length unit x, kF (x, ξ n ) ∈ [0, 1]. kB (x, ξ n ) Output rate of braking force per length unit x, kB (x, ξ n ) ∈ [0, 1]. vn (x, ξ n ) Train speed per length unit x at inter-station (n, n + 1 ). tn (x, ξ n ) Running time per length unit x at inter-station (n, n + 1 ). Note that din − ain denotes the dwell time at station n; ai(n+1 ) − din denotes the running time at inter-station (n, n + 1 ); and kF (x, ξ n ), kB (x, ξ n ), vn (x, ξ n ), tn (x, ξ n ) and ai(n+1 ) − din can determine the optimal speed profile at inter-station (n, n + 1 ).

4.2. Uncertainty of train mass Train mass consists of the empty vehicle mass and mass of passengers on the train. The empty vehicle mass is generally fixed, but the mass of passengers on the train is spatially and temporally uncertain. That is, the passenger volumes are not only spatially uncertain across different inter-stations, but also temporally uncertain across different operational periods as indicated in Fig. 5. Take the Beijing Metro Yizhuang Line as an example. The daily operational period of this line is from 5:00 am to 11:00 pm. The peak hours are from 6:30 am to 9:30 am and from 4:30 pm to 7:30 pm, the night hours are from 9:30 pm to 11:00 pm, and the remaining operation time is off-peak hours. In order to obtain the accurate passenger mass during different operational periods at different inter-stations, the Beijing Jiaotong University (BJTU) group has counted the number of passengers on the train during peak hours, off-peak hours and night hours on March 23, 2014 (Sunday) and March 24, 2014 (Monday). The average numbers of passengers on the train at each inter-station during different operational periods are shown in Fig. 5. Note that the Songjiazhuang station (last station of the Yizhuang Line) is located in the city center of Beijing that connects to the Beijing Metro network (also see Fig. 13(a)); hence the passenger volume increases along the up direction to the Songjiazhuang station. In addition, the Yizhuang station has not yet opened to the public on March 23, 2014; hence there are only hired staffs (for counting the number of passengers) on the train between the Yizhuang and Ciqu stations. As described in the problem statement, we want to generate a set of dwell time, running time and speed profile suitable for the overall operational period. We use p1n , p2n and p3n to denote the probabilities for peak hours, off-peak hours and night hours. Thus, the train mass at inter-station (n, n + 1 ) can be defined as a discrete stochastic variable, i.e.,



m + τ1n , with probability p1n , ξn = m + τ2n , with probability p2n , m + τ3n , with probability p3n ,

(4.1)

where m is the empty vehicle mass; τ 1n , τ 2n and τ 3n are the masses of passengers on the train at inter-station (n, n + 1 ) during peak hours, off-peak hours and night hours, respectively; p1n , p2n and p3n are positive real numbers between (0, 1) satisfying p1n + p2n + p3n = 1. As the train mass at each inter-station is a stochastic variable, there are 2N-2 different stochastic variables in the formulation.

430

X. Yang et al. / Transportation Research Part B 91 (2016) 424–445

Peak hours Off-peak hours Night hours

1000 800

Up direction

600 400

Number of passengers

Songjiazhuang

Xiaocun

Xiaohongmen

Jiugong

Yizhuangqiao

0

Wenhuayuan

Wanyuan

Tongjinan

Jinghai

Ciqunan

Ciqu

Station

Yizhuang

0

Rongjing

200 Rongchang

Number of passengers

1200

200 400 600 800 1000

Peak hours Off-peak hours Night hours

Down direction

1200 Fig. 5. Average numbers of passengers on the train at each inter-station (Source: Xun et al. (2014)).

4.3. System constraints The intention of the integrated optimization model is to determine the optimal timetable and speed profile. The optimality of the timetable and speed profile is interdependent (i.e., both sets of decision variables are driven by the objective function given in Section 4.5). Therefore, we divide the modeling process into two phases: Given the speed profile, the first phase determines the timetable by scheduling the arrival and departure times for each station to obtain the dwell time at each station and the running time at each inter-station, and estimates the total tractive energy consumption for one cycle trip (i.e., departing from station 1 arriving at station 2N); the second phase determines the speed profile for each interstation with the scheduled arrival and departure times from the timetable, and calculates the tractive energy consumption for each inter-station. Remark 2. Train mass at each inter-station is a stochastic variable based on distinct probabilities for peak hours, off-peak hours and night hours (see Eq. 4.1). Therefore, the calculation of the tractive energy consumption for each inter-station is a stochastic evaluation. 4.3.1. Schedule phase For the schedule phase, the model should satisfy the following constraints: (1) Departure time constraints For each n ∈ [1, 2N − 1], the optimal departure time for train i at station n should be an integer and valued between the lower and upper bounds of departure time for train i at station n, i.e.,

din ∈ {0, 1, 2, · · · , C }, dcin ≤ din ≤ dcin ,

∀ n ∈ [1, N − 1] ∪ [N + 1, 2N − 1].

(4.2)

(2) Arrival time constraints For each n ∈ [2, 2N], the optimal arrival time for train i at station n should be an integer and valued between the lower and upper bounds of arrival time for train i at station n, i.e.,

ain ∈ {0, 1, 2, · · · , C }, acin ≤ ain ≤ acin ,

∀ n ∈ [2, N] ∪ [N + 2, 2N].

(4.3)

(3) Running time constraints For each n ∈ [1, N − 1] ∪ [N + 1, 2N − 1], the optimal running time for train i at inter-station (n, n + 1 ) should be valued between the lower and upper bounds of the current running time for train i at inter-station (n, n + 1 ), i.e.,

tcn ≤ |ai(n+1) − din | ≤ tcn ,

∀ n ∈ [1, N − 1] ∪ [N + 1, 2N − 1].

(4.4)

(4) Dwell time constraints For each n ∈ [2, N − 1] ∪ [N + 2, 2N − 1], the optimal dwell time for train i at station n should be equal to the current dwell time for train i at station n, i.e.,

din − ain = dcin − acin ,

∀ n ∈ [2, N − 1] ∪ [N + 2, 2N − 1].

(4.5)

X. Yang et al. / Transportation Research Part B 91 (2016) 424–445

431

(5) Turnaround time constraints The turnaround time in the optimal timetable should be equal to the current turnaround time, i.e.,

di(N+1) − aiN = tt ,

(4.6)

where the dwell time at stations N and N + 1 is included in the turnaround time. (6) Total travel time constraints The total travel time for one cycle in the optimal timetable should be equal to the current total travel time for one cycle, i.e.,

ai(2N ) − di1 = C.

(4.7)

4.3.2. Speed profile phase Given an inter-station (n, n + 1 ) with length Xn and scheduled running time Tn = ai(n+1 ) − din , denote tn (x, ξ n ), vn (x, ξ n ), gn (x) and Vn (x) as the running time, train speed, gradient force and speed restriction per length unit x, respectively. For each n ∈ [1, N − 1] ∪ [N + 1, 2N − 1], the model in the speed profile phase should satisfy the following constraints: (7) Speed restriction constraints The train speed should be no more than the speed restriction per length unit x and should be zero at the starting and ending points, i.e.,

E[vn (0, ξn )] = 0, E[vn (Xn , ξn )] = 0, E[vn (x, ξn )] ≤ Vn (x ),

∀ n ∈ [1, N − 1] ∪ [N + 1, 2N − 1],

(4.8)

∀ x ∈ (0, Xn ), n ∈ [1, N − 1] ∪ [N + 1, 2N − 1].

(4.9)

(8) Complete running time constraints The accumulated running time per length unit x at inter-station (n, n + 1 ) should be equal to the scheduled running time at inter-station (n, n + 1 ), i.e.,



Xn

E 0



tn (x, ξn )dx = Tn ,

∀ n ∈ [1, N − 1] ∪ [N + 1, 2N − 1].

(4.10)

(9) State variable constraints Based on the optimal train control theory, the general equations of motion are



E

dtn (x, ξn ) − dx

1 vn (x, ξn )



= 0,

∀ n ∈ [1, N − 1] ∪ [N + 1, 2N − 1],



(4.11)



dvn (x, ξn ) kF (x, ξn )F [vn (x, ξn )] + kB (x, ξn )B[vn (x, ξn )] + r[vn (x, ξn )] + gn (x ) E − = 0, dx ξn vn (x, ξn )

∀ n ∈ [1, N − 1] ∪ [N + 1, 2N − 1],

(4.12)

where the maximum tractive force F(v) is non-negative; the maximum braking force B(v) is non-positive; the basic running resistance r(v) is negative; the gradient force gn (x) can be positive, zero or negative. Eqs. (4.11) and (4.12) denote the time state and speed state per length unit. The complexity of evaluating these equations is much increased since the maximum tractive force, maximum braking force and basic running resistance are all variable according to the train speed variation. Take the trains operated in the Beijing Metro Yizhuang Line as an example, these real-world operating conditions are graphically depicted in Fig. 6 along with the formulas. 4.4. Objective function The operations of metro trains are considered as periodic by the following two points: First, trips from one train at different cycles are distributed with the same dwell time at stations and the same speed profile at inter-stations; second, trips from different trains are distributed with the same dwell time at stations and the same speed profile at inter-stations. For simplicity, we consider the total tractive energy consumption of one train completing one cycle trip (i.e., departing from station 1 arriving at station 2N) for calculating the objective function. Therefore, the total tractive energy consumption of all trains throughout the overall operational period can be simply summed. We first describe the tractive energy consumption of one train at one inter-station. For each n ∈ [1, N − 1] ∪ [N + 1, 2N − 1], the tractive energy consumption of train i at inter-station (n, n + 1 ) is formulated as

J (din , ai(n+1) , ξn ) =



Xn 0

kF (x, ξn )F [vn (x, ξn )] dx,

(4.13)

where x is per length unit; Xn is the length of inter-station n; F[vn (x, ξ n )] is the maximum tractive force per length unit x (see Fig. 6(b)); kF (x, ξ n ) is the output rate of tractive force per length unit x; kF (x, ξ n )F[vn (x, ξ n )] is the final output

432

X. Yang et al. / Transportation Research Part B 91 (2016) 424–445

Fig. 6. Basic running resistance and maximum tractive and braking forces for trains operated in the Beijing Metro Yizhuang Line (Xiao, 2014).

tractive force on train i per length unit x. The tractive energy consumption is obtained using the final output tractive force to integrate over the length of the inter-station. Remark 3. Equation (4.13) shows that the tractive energy consumption of train i at each inter-station contains an independently distributed discrete random variable. The total tractive energy consumption of train i completing one cycle trip is the sum of the tractive energy consumption of train i at each inter-station. Therefore, it contains 2N − 2 independently distributed discrete random variables. In the following proposition, we want to illustrate that the expected value of the total tractive energy consumption completing one cycle trip is the sum of the expected value of the tractive energy consumption at each inter-station. Proposition 1. (Linearity of the Expected Value) Let ξ 1 , ξ 2 , , ξ N be N independently distributed discrete random variables with finite expected values, we have E (ξ1 + ξ2 + · · · + ξN ) = E (ξ1 ) + E (ξ2 ) + · · · + E (ξN ).

X. Yang et al. / Transportation Research Part B 91 (2016) 424–445

Proof. We first prove E (ξ1 + ξ2 ) = E (ξ1 ) + E (ξ2 ). {y j with probability p j | j = 1, 2, · · · , J}. We have

E ( ξ1 + ξ2 )

L

Let

ξ1 = { z l

with probability

433

pl

| l = 1, 2, · · · , L} and ξ2 =

J

( zl + y j ) pl p j   z p p + Ll=1 Jj=1 y j pl p j j=1 l l j l=1  J L y p l=1 zl pl + j=1 j j E ( ξ1 ) + E ( ξ2 ) .

=

l=1 Jj=1 L

= = =

(4.14)

Without loss of generality, we have E (ξ1 + ξ2 + · · · + ξN ) = E (ξ1 ) + E (ξ2 ) + · · · + E (ξN ). Hence the expected value is linear, and the proof is complete.  For simplicity, we denote d = {din |n = 1, 2, · · · , 2N − 1}, a = {ain |n = 2, 3, · · · , 2N}, and ξ = {ξn |n = 1, 2, · · · , N − 1, N + 1, · · · , 2N − 1}. Based on the linearity of the expected value, the expected value of the total tractive energy consumption of train i completing one cycle trip (i.e., departing from station 1 arriving at station 2N) is the sum of the expected value of tractive energy consumption of train i at each inter-station. It is formulated as

E[JTotal (d, a, ξ )]

= =

 2N−1 E[ N−1 n=1 J (din , ai(n+1 ) , ξn ) +  n=N+1 J (din , ai(n+1 ) , ξn )] N−1 2N−1 E[ J ( d , a , ξ ) ] + n in i(n+1 ) n=1 n=N+1 E[J (din , ai(n+1 ) , ξn )].

(4.15)

This final objective function explicitly considers the uncertain train mass and the variable tractive force. 4.5. Stochastic optimization model Based on the above analysis, the integrated timetable and speed profile optimization problem with uncertain train mass is formulated as a two-phase stochastic programming model. Since the objective function contains stochastic variables, we take the expected value criterion to minimize the average total tractive energy consumption value. The two-phase stochastic programming model is formulated as:

⎧ min ⎪ ⎪ ⎪ ⎪ s.t. ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩

E[JTotal (d, a, ξ )] din ∈ {0, 1, 2, · · · , C }, dcin ≤ din ≤ dcin , ∀ n ∈ [1, N − 1] ∪ [N + 1, 2N − 1] ain ∈ {0, 1, 2, · · · , C }, acin ≤ ain ≤ acin , ∀ n ∈ [2, N] ∪ [N + 2, 2N] tcn ≤ |ai(n+1) − din | ≤ tcn , ∀ n ∈ [1, N − 1] ∪ [N + 1, 2N − 1] din − ain = dcin − acin , ∀ n ∈ [2, N − 1] ∪ [N + 2, 2N − 1] di(N+1) − aiN = tt ai(2N ) − di1 = aci(2N ) − dci1 where the expected value J (din , ai(n+1) , ξn ) of the JTotal (d, a, ξ ) is computed by  Xn E[J (din , ai(n+1) , ξn )] = min E { kF (x, ξn )F [vn (x, ξn )] dx} 0

s.t. E[vn (0, ξn )] = 0, E[vn (Xn , ξn )] = 0, ∀ n ∈ [1, N − 1] ∪ [N + 1, 2N − 1] E[vn (x, ξn )] ≤ Vn (x ), ∀ x ∈ (0, Xn ), n ∈ [1, N − 1] ∪ [N + 1, 2N − 1]  Xn E[ tn (x, ξn )dx] = ai(n+1) − din , ∀ n ∈ [1, N − 1] ∪ [N + 1, 2N − 1]

(4.16)

0

1 dtn (x, ξn ) E[ − ] = 0, ∀ n ∈ [1, N − 1] ∪ [N + 1, 2N − 1] dx vn (x, ξn ) dvn (x, ξn ) kF (x, ξn )F [vn (x, ξn )] + kB (x, ξn )B[vn (x, ξn )] + r[vn (x, ξn )] + gn (x ) E[ − ] = 0, dx ξn vn (x, ξn ) ∀ n ∈ [1, N − 1] ∪ [N + 1, 2N − 1].

We can see that the first phase determines the timetable by scheduling the arrival and departure times (i.e., a and d) constrained by the speed profile (i.e., kF (x, ξ n ), kB (x, ξ n ), vn (x, ξ n ) and tn (x, ξ n )), and the second phase determines the speed profile constrained by the scheduled arrival and departure times. 4.6. Convergence property The stochastic optimization model is usually converted into a deterministic one by replacing the objective function to its empirical mean for illustrating the convergence property (Dai et al., 20 0 0). We proof the convergence property of the developed model (4.16) as follows. Proposition 2. We denote c = {d, a} as the decision variable vector and X as the feasible domain. The objective function in model (4.16) can be rewritten as

H (c ) = E[JTotal (c, ξ )],

(4.17)

434

X. Yang et al. / Transportation Research Part B 91 (2016) 424–445

and we would like to find an c∗ ∈ X that solves the optimization model (4.16), i.e., c∗ should satisfy

H (c∗ ) = min H (c ).

(4.18)

c∈X

We want to show that c∗ and H(c∗ ) are convergent. Proof. We replace H(c) by its empirical mean based on a total of Z samples. Hence, the problem is to find an cZ∗ by solving

min HZ (c ) 

c∈X

Z 1 JTotal (c, ξk ), Z

(4.19)

k=1

where {ξk , k = 1, 2, · · · , Z } are samples of ξ . As E[|JTotal (c, ξ )|] < ∞ for all c ∈ X , the strong law of large numbers guarantees that

HZ (c ) → H (c ),

(4.20) cZ∗

c∗

for each fixed c ∈ X (Shiryayev, 1984). Then it is proved → in a certain sense as Z → +∞ (King and Rockafellar, 1993), and H (cZ∗ ) → H (c∗ ) when the samples are independently and identically distributed (Kaniovski et al., 1995). Since minimizing HZ (c) is a deterministic problem, the convergence of cZ∗ and H (cZ∗ ) can be proved (Birge and Louveaux, 2011), and we have obtained cZ∗ → c∗ and H (cZ∗ ) → H (c∗ ). Therefore, c∗ and H(c∗ ) are also convergent. The proof is complete.  4.7. Formulation complexity Solving the two-phase stochastic programs usually has a number of difficulties, and we analyze the complexity of the formulated two-phase stochastic programming model. For the first phase, we need to determine the timetable by scheduling the arrival and departure times for each station, and the number of all possible arrival and departure times combinations is

 = [(acin − acin + 1 )(dcin − dcin + 1 )]2N−2 .

(4.21)

In the numerical examples (to be shown in Section 5), the number  is 8126 . Clearly, it is a huge number. Hence, it is challenging to find the optimal arrival and departure times from all these possible combinations. For the second phase, we need to determine an energy-efficient speed profile for each inter-station with the scheduled arrival and departure times from the timetable. It has two levels of difficulties: (1) the tractive energy consumption and energy-efficient speed profile are nonlinear relative to the arrival and departure times due to the variable tractive force, braking force and basic running resistance shown in Fig. 6, and (2) we need to repeat a number of times to solve the stochastic program for obtaining the expected value considering the uncertain train mass. Obviously, the formulated two-phase stochastic programming model is very complex, and it is very difficult to find the optimal solution using the classical optimization methods such as branch-and-bound algorithm and Newton algorithm. 5. Solution method Evolutionary algorithms, particularly the genetic algorithm, have shown to be effective in solving complex problems. Due to its extensive generality, strong robustness, high efficiency and practical applicability, genetic algorithm has been widely adopted as a numerical method for solving many problems, including transportation network design (Cacchiani et al., 2010; Chen et al., 2010; Chen and Subprasom, 2007), transit network design (Bagloee and Ceder, 2011; Kim and Schonfeld, 2014), pedestrian crossing design (Yu et al., 2015), and train timetable optimization (Niu and Zhou (2013); Xu et al. (2015)). In this section, we first use Fig. 7 to show the whole structure of the solution method to the two-phase stochastic program. Then we design a simulation-based genetic algorithm procedure shown in Fig. 8 to find the optimal solution. The simulation-based genetic algorithm procedure mainly includes two modules, i.e., simulation module and genetic algorithm module. Simulation module is designed based on the optimal train control algorithm to determine the energyefficient speed profile and to calculate the expected value of tractive energy consumption for each inter-station. Genetic algorithm module (including the representation, selection, crossover and mutation operations) is adopted to determine the timetable by scheduling the arrival and departure times for each station. 5.1. Simulation module Simulation module present the detailed procedure to obtain the energy-efficient speed profile for an metro inter-station with a given length and running time. Proposition 3. (Existence of an optimal speed profile) Let Tn ≥ Tnmin (Tnmin denotes the minimum running time at inter-station (n, n + 1 )). There exists a control scheme kF = kF (x, ξn ) and kB = kB (x, ξn ), a corresponding speed sequence v = vn (x, ξn ), and a time sequence t = tn (x, ξn ) satisfying the differential constraints vv = kF /v + kB + r (v ) + g and t = 1/v, the initial conditions v(0 ) = 0 and v(Xn ) = 0, and time constraints t (0 ) = 0 and t (Xn ) = 0 in such a way that the cost J = 0Xn kF /v dx is minimized.

X. Yang et al. / Transportation Research Part B 91 (2016) 424–445

435

Fig. 7. Structure of the solution method.

1. Define maximum number of generations max_generation (index: j1) 2. Define population size pop_size (index: j2) 3. Define number of inter-stations of 2N-2 (index: n) 4. Generate initial population

GA module 1. Representation 2. Selection 3. Crossover 4. Mutation

j1 = 1 j1 = j1 + 1 j2 = 1 j2 = j2 + 1 n=1 n=n+1 Simulation module

No

No

No

n > 2N-2

Yes

j2 > pop_size

Yes

j1 > max_generation Yes

Report final solutions Fig. 8. Flowchart of the simulation-based genetic algorithm procedure.

Outline of proof. The Hamiltonian function (Howlett et al., 2009) is defined as

H=

β − kF α + 2 [kF + kB v + r (v )v + gv], v v

where the adjoint variables α = α (x ) and β = β (x ) are solutions to the differential equations

α =

β − kF α + 3 [2kF + kB v + r (v )v − r (v )v2 + gv] and β = 0. v2 v

436

X. Yang et al. / Transportation Research Part B 91 (2016) 424–445

We wish to maximize the Hamiltonian subject to the constraints kF ∈ [0, 1] and kB ∈ [0, 1]. Thus we define a Lagrangian function H = H + ρ kF + σ ( 1 − kF ) + τ kB + ω ( 1 − kB ), where ρ , σ , τ and ω are the Lagrange multipliers. Maximizing H subject to the given constraints, there exists five possible phases of the optimal control (Howlett and Pudney, 1995; Liu and Golovitcher, 2003): (1) (2) (3) (4) (5)

maximum acceleration: α > v, kF = 1 and kB = 0; speed-holding: α = v, 0 ≤ kF ≤ 1 and kB = 0; coasting: 0 < α < v, kF = 0 and kB = 0; partial braking: α = 0, kF = 0 and 0 ≤ kB ≤ 1; maximum braking: α < 0, kF = 0 and kB = 1.

At the same time, Howlett and Pudney (1995) also proved the speed-holding phase only occurs in long inter-stations and the partial braking phase cannot occur. The inter-stations in metro systems are often with short travel distance and approximatively constant gradient track (i.e., gradients only have small changes within an inter-station as the inter-station length is short). Albrecht et al. (2015a); Howlett and Pudney (1995), and Albrecht et al. (2015b) have shown that the most energy-efficient control strategy for a metro inter-station only contains the maximum acceleration, coasting and maximum braking. Remark 4. For the optimal train control strategy, the initial results by Milroy (1981) suggested that an energy-efficient speed profile should consist of three phases: maximum acceleration, coasting and maximum braking. At this early stage, it was not understood that the speed-holding phase should be used in some conditions. Then Lee et al. (1982) concluded the speed-holding phase should be included after the maximum acceleration phase when the running time is sufficient. This optimal control strategy is widely used in general railway systems (Albrecht et al., 2013; Howlett, 20 0 0). Later, the Transport Control Group at the South Australian Institute of Technology developed the Metromiser system (Howlett and Pudney, 1995) and showed that the speed-holding phase only occurs in long inter-stations. That is to say, the most energy-efficient control strategy only contains the maximum acceleration, coasting and maximum braking for a typical metro journey interstation with short travel distance. For the key principles of optimal train control strategies, readers can refer to the recent comprehensive review by Albrecht et al. (2015a, 2015b). Given an inter-station (n, n + 1 ) with length Xn and running time Tn , we divide the inter-station into many units (e.g., 1 m). For simplicity, we denote vn (x, ξn ) = {vn (x, ξn )|x = 1, 2, · · · , Xn }. We define X1 as the switching point from the maximum acceleration phase to the coasting phase, and X2 as the switching point from the coasting phase to the maximum braking phase. The procedure for obtaining the energy-efficient speed sequences of the inter-station (n, n + 1 ) is described as follows: Step 1. Set a matrix U = 0 and a real number W = 0. Step 2. Generate ξ n based on the probabilities. Step 3. Set X1 = 0 and X2 = Xn . Step 4. Set X1 = X1 + 1, X2 = X2 − 1, x = 0 and vn (x, ξ n ). Step 5. For 0 ≤ x < X1 , generate v2n (x + 1, ξn ) − v2n (x, ξn ) = 2(F [vn (x, ξn )] + r[vn (x, ξn )] + gn (x )). For X1 ≤ x < X2 , generate v2n (x + 1, ξn ) − v2n (x, ξn ) = 2(r[vn (x, ξn )] + gn (x )). For X2 ≤ x < Xn , generate v2n (x + 1, ξn ) − v2n (x, ξn ) = 2(B[vn (x, ξn )] + r[vn (x, ξn )] + gn (x )). Step 6. For 0 < x ≤ Xn , calculate tn (x, ξn ) = 1/vn (x, ξn ).  Xn Step 7. If vn (Xn , ξn ) = 0 and Tn = xx= t (x, ξn ), obtain vn (x, ξ n ) and J (din , ai(n+1 ) , ξn ). Otherwise, go to step 4. =1 n Step 8. Calculate U = U + vn (x, ξn ) and W = W + J (din , ai(n+1 ) , ξn ). Step 9. Repeat the step 2 to step 8 for Y times, where Y is a sufficiently large positive integer. Step 10. Return vn (x, ξn ) = U /Y and E[J (din , ai(n+1 ) , ξn )] = W/Y . 5.2. Genetic algorithm module Genetic algorithm usually starts with an initial set of randomly generated feasible solutions, which are encoded as chromosomes called a population. The decision variables, i.e., departure time d = {din |n = 1, 2, · · · , 2N − 1} and arrival time a = {ain |n = 2, 3, · · · , 2N}, constitute a chromosome L = {ln |n = 1, 2, · · · , 2N − 2} (see Fig. 9). A new population of chromosomes is generated following the evaluation, selection, crossover and mutation operations. Genetic algorithm will terminate after a given number of iterations of the above steps. The process can be summarized as follows: Step 1. Initialize parameters: population size pop_size, crossover probability Pc , mutation probability Pm , and max_generation. Set generation index j1 = 1. Step 2. Initialize pop_size feasible chromosomes as the initial population. Step 3. Calculate the evaluation function values for all chromosomes. Step 4. Select the chromosomes by spinning the roulette wheel. Step 5. Produce the next generation through crossover and mutation operations. Step 6. If j1 = max_generation, return the best found solution. Otherwise, set j1 = j1 + 1, and go to step 3.

X. Yang et al. / Transportation Research Part B 91 (2016) 424–445

437

Fig. 9. Structure of a chromosome.

Moving direction Station 1

Length: 1334 m Original running time: 110s

Station 2

Station 3

Length: 1286 m Original running time: 100s

Train mass: 200704 kg, with probability 0.33 Train mass: 205630 kg, with probability 0.33 200680 kg, with probability 0.59 204010 kg, with probability 0.59 200580 kg, with probability 0.08 201940 kg, with probability 0.08 Expected value of total tractive energy consumption is 12.41 k·Wh with the current timetable and speed profile Fig. 10. A simple one-way line.

20

Speed (m/s)

(1) Optimal speed profile chooses larger tractive and braking force

Current speed profile Optimal speed profile

15 10 5

Moving direction

0

Station 1

Station 2

Station 3

Time (s)

Current timetable Optimal timetable

(2) Optimal timetable reallocates running times as 106 s and 104 s Fig. 11. Obtained optimal timetable and speed profile in comparison with the current timetable and speed profile.

6. Numerical example In this section, we present two numerical examples to illustrate the developed model and solution method. Example 1 is a simple one-way line used to examine the effectiveness of the designed simulation-based genetic algorithm procedure. Example 2 is based on the operation data (e.g., speed restrictions, gradients, vehicle tractive and braking forces, basic running resistance, passengers demand, etc.) from the Beijing Metro Yizhuang Line to examine the practicability of the developed two-phase stochastic programming model to a real-world metro line. 6.1. Example 1: a simple line In example 1, a simple one-way line shown in Fig. 10 consists of 3 stations and 2 inter-stations. As the solution domain is limited, we use the enumeration method to accurately calculate the minimum expected value of total tractive energy consumption. The optimal value is 11.25 kW · h, which can reduce the total tractive energy consumption by (12.41 − 11.25 )/12.41 = 9.35%. The obtained optimal timetable and speed profile is shown in Fig. 11 in comparison with the current timetable and speed profile. For examining the effectiveness of the simulation-based genetic algorithm procedure, we perform it with pop_size = 6, max_generation = 15, Pc = 0.6 and Pm = 0.15 for 10 times. The computation results are presented in Fig. 12 and recorded

438

X. Yang et al. / Transportation Research Part B 91 (2016) 424–445

(b) Calculation results for the 4th, 5th and 6th time 1st 2nd 3rd

12.2 12

Best found value (kW.h)

Best found value (kW.h)

(a) Calculation results for the 1st, 2nd and 3rd time 12.4

11.8 11.6 11.4 11.2 0

5

10

12.4

4th 5th 6th

12.2 12 11.8 11.6 11.4 11.2 0

15

5

Iteration index

7th 8th 9th

12.2 12 11.8 11.6 11.4 5

15

(d) Calculation results for the 10th time

Best found value (kW.h)

Best found value (kW.h)

(c) Calculation results for the 7th, 8th and 9th time 12.4

11.2 0

10 Iteration index

10

12.4

10th

12.2 12 11.8 11.6 11.4 11.2 0

15

5

Iteration index

10

15

Iteration index

Fig. 12. Detailed computation results by the simulation-based genetic algorithm procedure for 10 times. Table 2 Best found minimum expected value of total tractive energy consumption by the simulation-based genetic algorithm procedure. Time

1st

2nd

3rd

4th

5th

6th

7th

8th

9th

10th

Computation result (kW · h) Error

11.25 0

11.25 0

11.31 0.53%

11.25 0

11.25 0

11.25 0

11.25 0

11.25 0

11.25 0

11.25 0

in Table 2. The results show that the average error by the simulation-based genetic algorithm procedure is 0.53%/10 = 0.053% in comparison with the optimal value 11.25 kW · h, which demonstrates that the designed simulation-based genetic algorithm procedure is effective.

6.2. Example 2: Beijing Metro Yizhuang Line In example 2, we firstly provide the real-world operation data of the Beijing Metro Yizhuang Line, and then present two cases to illustrate the practicability of the developed model. Case 1 shows the optimal timetable and speed profile of the developed model using the provided real-world operation data and makes a comparison with the current timetable and speed profile. Case 2 makes a comparison between the previous deterministic integrated optimization model (Li and Lo, 2014a) and the developed stochastic integrated optimization model (considering the uncertainty of train mass and the variability of tractive force, braking force and basic running resistance).

6.2.1. Basic data The Yizhuang Line links the downtown of Beijing and the Yizhuang Economic Development Zone (see lower-right corner of Fig. 13(a)), which consists of 14 stations covering a length of 22.73 km (see Fig. 13(b)). We obtained the current operation data of the Beijing Metro Yizhuang Line from the Beijing Mass Transit Railway Operation Corporation Limited (Zhang, 2014). The current length, running time, the bounds of running time for each inter-station and the dwell time for each station are provided in Table 3. The speed restriction and gradient in the down direction are shown in Fig. 14. The empty vehicle mass is 199,0 0 0 kg, and the masses of passengers on the train at each inter-station during different periods are provided in Table 4. The parameters in the basic running resistance, maximum tractive force and maximum braking force formulas are: λ1 = 3779.90 N, λ2 = 1.57, λ3 = 8.52 ∗ 10−4 , λ4 = 310, 0 0 0 N, λ5 = 86, 228 N, λ6 = 960.00, λ7 = 260, 0 0 0 N, λ8 = 6265.30, vα = 10.00 m/s, vmax = 22.22 m/s and vβ = 16.67 m/s. The remaining parameters are listed as follows: N = 14, tt = 240 s, and C = 4572 s.

X. Yang et al. / Transportation Research Part B 91 (2016) 424–445

Changping Line

Line 5

Line 15

Songjiazhuang

Line 8

Line 4

Airport Line

Line 13

Xiaocun Xiaohongmen Jiugong

Line 10

Line 6

Line 1

439

Wenhuayuan

Line 6 Line 2

Line 14

Down direction

Line 2 Line 10

Line 5

Yizhuangqiao Wanyuan

Line 1

Line 9

Up direction

Rongjing

Fangshan Line

Rongchang

Yizhuang Line

Jinghai

Tongjinan

Ciqunan

Ciqu Yizhuang

Line 4 (a) Beijing Metro Network

(b) Yizhuang Line

Fig. 13. Beijing Metro network and the Yizhuang Line.

Table 3 Current length, current running time, the bounds of running time and dwell time for the Beijing Metro Yizhuang Line. Inter-station# Departure - Arrival

Length (m)

Running time (s)

Lower bound of the running time (s)

Upper bound of the running time (s)

Dwell time at the arrival station (s)

Up direction Yizhuang–Ciqu Ciqu–Ciqunan Ciqunan–Jinghai Jinghai–Tongjinan Tongjinan–Rongchang Rongchang–Rongjing Rongjing–Wanyuan Wanyuan–Wenhuayuan Wenhuayuan–Yizhuangqiao Yizhuangqiao–Jiugong Jiugong–Xiaohongmen Xiaohongmen–Xiaocun Xiaocun–Songjiazhuang

1334 1286 2086 2265 2338 1354 1280 1538 993 1982 2366 1275 2631

110 100 141 150 162 103 101 111 90 135 157 105 195

106 96 137 146 158 99 97 107 86 131 153 101 191

114 104 145 154 166 107 105 115 94 139 161 109 199

45 35 30 30 30 30 30 30 35 30 30 30 –

Down direction Songjiazhuang–Xiaocun Xiaocun–Xiaohongmen Xiaohongmen–Jiugong Jiugong–Yizhuangqiao Yizhuangqiao–Wenhuayuan Wenhuayuan–Wanyuan Wanyuan–Rongjing Rongjing–Rongchang Rongchang–Tongjinan Tongjinan–Jinghai Jinghai–Ciqunan Ciqunan–Ciqu Ciqu–Yizhuang

2631 1275 2366 1982 993 1538 1280 1354 2338 2265 2086 1286 1334

190 108 157 135 90 114 103 104 164 150 140 102 105

186 104 153 131 86 110 99 100 160 146 136 98 101

194 112 161 139 94 118 107 108 168 154 144 106 109

30 30 30 35 30 30 30 30 30 30 35 45 –

# Note. Abbreviation of stations name: Songjiazhuang (SJZ), Xiaocun (XC), Xiaohongmen(XHM), Jiugong (JG), Yizhuangqiao (YZQ), Wenhuayuan (WHY), Wanyuan (WY), Rongjing (RJ), Rongchang (RC), Tongjinan (TJ), Jinghai (JH), Ciqunan (CQN), Ciqu (CQ), Yizhuang (YZ).

6.2.2. Case 1 This case provides the results of the developed two-phase stochastic programming model and makes a comparison with the current operating timetable and speed profile to illustrate the effectiveness of the developed model on energy conservation. Similar to example 1, we test the behavior of the simulation-based genetic algorithm procedure on a real case study by running for 3 times with pop_size = 30, max_generation = 120, Pc = 0.6 and Pm = 0.15 on a Windows 8.1 platform of personal computer with processor frequency of 2.4 GHz and memory size of 8 GB. The results presented in Fig. 15 show

440

X. Yang et al. / Transportation Research Part B 91 (2016) 424–445

25

Speed (m/s)

20 15 10 5 Speed restriction 0 0

0.5

1 Position in down direction (m) (a) Speed restriction

1.5

2 x 10

4

Gradient force (N ‰)

20 Gradient 10 0 -10 -20 -30 0

0.5

1 Position in down direction (m) (b) Gradient

1.5

2 x 10

4

Fig. 14. Speed restriction and gradient in the down direction. Table 4 Masses of passengers on the train at each inter-station during different periods. Inter-station (n, n + 1 ) At peak hours (kg) At off-peak hours (kg) At night hours (kg)

YZ–CQ 1680 1680 1680

CQ–CQN 6630 5010 2940

CQN–JH 7620 5640 30 0 0

JH–TJ 14,460 8820 4620

TJ-RC 28,140 18,330 8160

RC–RJ 33,990 22,590 10,260

RJ–WY 35,550 25,890 11,880

WY–WHY 38,070 28,590 12,840

WHY–YZQ 42,750 30,120 14,280

Inter-station (n, n + 1 ) At peak hours (kg) At off-peak hours (kg) At night hours (kg)

YZQ–JG 47,910 33,210 15,600

JG–XHM 57,960 36,360 16,320

XHM–XC 59,100 37,290 16,800

XC–SJZ 59,910 39,810 17,280

SJZ–XC 69,150 46,170 34,680

XC–XHM 65,790 43,740 33,120

XHM–JG 64,290 42,090 31,260

JG–YZQ 56,670 35,040 24,720

YZQ–WHY 51,900 31,890 21,540

Inter-station (n, n + 1 ) At peak hours (kg) At off-peak hours (kg) At night hours (kg)

WHY–WY 47,850 30,120 19,680

WY–RJ 38,940 25,890 18,120

RJ–RC 32,310 23,610 16,620

RC–TJ 20,640 17,460 13,680

TJ–JH 9840 8790 6300

JH–CQN 6210 5940 4800

CQN–CQ 5760 5790 4560

CQ-YZ 1680 1680 1680

that the algorithm is convergent after the 60th generation for all three experiments. Therefore, we perform the solution procedure with pop_size = 30, max_generation = 60, Pc = 0.6 and Pm = 0.15 in the following cases. We obtain the optimal timetable presented in Table 5 and also shown in Fig. 16 in comparison with the current timetable, and the optimal speed profile shown in Fig. 17 in comparison with the current speed profile. Based on the optimal timetable and speed profile, the total tractive energy consumptions of one train for completing one cycle trip during peak hours, offpeak hours and night hours for the optimal timetable and speed profile are 210.84 kW · h, 200.28 kW · h and 191.62 kW · h, respectively. They are also compared with those of the current timetable and speed profile in Table 6. The results show that the total tractive energy consumptions are reduced by 10.66%, 9.94% and 9.13% during peak hours, off-peak hours and night hours, respectively. As shown in Figs. 16 and 17, the maximum train speed, acceleration rate and deceleration rate of the optimal speed profile at each inter-station are larger than those of the current speed profile. It is because the optimal speed profile chooses larger tractive and braking forces to reduce the tractive energy consumption for each inter-station. Moreover, the optimal timetable increases the running time at inter-stations (e.g., XHM-XC, XC-SJZ, SJZ-XC, XC-XHM) with a larger train mass while decreases the running time at inter-stations (e.g., YZ-CQ, CQ-CQN, CQN-CQ, CQ-YZ) with a smaller train mass, such that the total tractive energy consumption can be further reduced. These changes in the timetable and speed profile, albeit minor in Figs. 16 and 17, can achieve an average of around 10% saving on the total tractive energy consumption. The results show that, in comparison with the current timetable and speed profile, the optimal timetable and speed profile: (1) reduces the maximum train speed at each inter-station by choosing larger tractive and braking forces, such that the tractive energy consumption for each inter-station can be reduced; and (2) increases the running time at inter-stations

X. Yang et al. / Transportation Research Part B 91 (2016) 424–445

441

218 1st experiment 2nd experiment 3rd experiment

Objective function value

216 214 212 210 208 206

Convergent

204 202

0

20

40

60 Iteration

80

100

120

Fig. 15. Convergence of the simulation-based genetic algorithm procedure. Table 5 Optimal timetable for one cycle for the Beijing Metro Yizhuang Line. Up station Arrival time (s) Departure time (s)

Yizhuang – 0

Ciqu 107 152

Ciqunan 250 285

Jinghai 423 453

Tongjinan 600 630

Rongchang 794 824

Rongjing 926 956

Up station Arrival time (s) Departure time (s)

Wanyuan 1055 1085

Wenhuayuan 1194 1224

Yizhuangqiao 1316 1351

Jiugong 1489 1519

Xiaohongmen 1679 1709

Xiaocun 1817 1847

Songjiazhuang 2046 –

Down station Arrival time (s) Departure time (s)

Songjiazhuang – 2286

Xiaocun 2480 2510

Xiaohongmen 2621 2651

Jiugong 2811 2841

Yizhuangqiao 2978 3013

Wenhuayuan 3105 3135

Wanyuan 3250 3280

Down station Arrival time (s) Departure time (s)

Rongjing 3384 3414

Rongchang 3517 3547

Tongjinan 3709 3739

Jinghai 3886 3916

Ciqunan 4053 4088

Ciqu 4186 4231

Yizhuang 4332 4572

Table 6 Total tractive energy consumptions of one train completing one cycle trip during different periods.

Current timetable and speed profile Optimal timetable and speed profile Energy saving rate

Peak hours

Off-peak hours

Night hours

236.01 kW · h 210.84 kW · h 10.66%

222.38 kW · h 200.28 kW · h 9.94%

210.87 kW · h 191.62 kW · h 9.13%

with a larger train mass while decreases the running time at inter-stations with a smaller train mass, such that the total tractive energy consumption can be further reduced. 6.2.3. Case 2 This case makes a comparison between the integrated optimization model (Li and Lo, 2014a) and the proposed stochastic integrated optimization model. For convenience, the two models are respectively named as the deterministic model and stochastic model for the comparison. For the deterministic model, we solve the optimal timetable and speed profile with the constant train mass (287,080 kg), tractive force (235,406 N), braking force (109,090 N) and basic running resistance (5742 N) (i.e., these variables take the same values in all inter-stations and during all periods). The calculation is based on the method provided in Li and Lo (2014a) by removing the calculation of regenerative braking energy. For the stochastic model, we solve the optimal timetable and speed profile with the uncertain train mass and variable tractive force, braking force and basic running resistance from the real-world operation data. The obtained total tractive energy consumptions during different operational periods for the two models (see the last row in Table 7) show that the developed stochastic model reduces the total tractive energy consumptions by 3.35%, 3.12% and 3.04% during peak hours, off-peak hours and night hours in comparison with the deterministic model.

442

X. Yang et al. / Transportation Research Part B 91 (2016) 424–445

15

15 Current timetable Optimal timetable

SJZ 13

Current timetable Optimal timetable

SJZ 14

Station index

Down direction

Up direction

11 9 7 5

13

12

3 1 0

YZ

1000

2000 3000 Time (s)

4000

YZ

5000

11 2250

2300

2350

2400 2450 Time (s)

2500

2550

2600

Fig. 16. Comparison of the current and optimal timetables for one cycle.

25

Speed (m/s)

Larger braking force

Larger tractive force

20

Current speed profile Optimal speed profile

15 Lower maximum speed

10 5 0 YZ

CQ

CQN

JH

TJ

RC

RJ

WY

WHY YZQ

JG

XHM

XC

SJZ

(a) Speed profile in up direction

25

Speed (m/s)

Larger braking force

Larger tractive force

20 15 10

Lower maximum speed

5 0 SJZ

XC

XHM

JG

YZQ WHY

WY

RJ

RC

TJ

JH

CQN

CQ

YZ

(b) Speed profile in down direction Fig. 17. Comparison of the current and optimal speed profiles for one cycle.

The obtained tractive energy consumption for each inter-station during different operational periods for the two models are presented in Table 7 and shown in Fig. 18. We can see that although the tractive energy consumptions for the interstations with a smaller train mass (e.g., CQN–CQ and CQ-YZ, see Fig. 18) of the stochastic model are a bit increased in comparison with the deterministic model, the total tractive energy consumption is reduced. It is because the stochastic model decreases the running time at inter-stations with a smaller train mass and allots the saved time to the inter-stations with a larger train mass. The results show that, in comparison with the deterministic model, the developed stochastic model: (1) can reduce the total tractive energy consumption during peak hours, off-peak hours and night hours; and (2) decreases the running time at inter-stations with a smaller train mass and reallocates the saved time to the inter-stations with a larger train mass, such that the total tractive energy consumption can be reduced.

X. Yang et al. / Transportation Research Part B 91 (2016) 424–445

443

Table 7 Comparisons on tractive energy consumption (unit: kW · h) between the deterministic and stochastic models. Inter-station

Deterministic model Peak hours

Off-peak hours

Stochastic model Night hours

Peak hours

Off-peak hours

Night hours

YZ–CQ CQ–CQN CQN–JH JH–TJ TJ–RC RC–RJ RJ–WY WY–WHY WHY–YZQ YZQ–JG JG–XHM XHM–XC XC–SJZ

6.6693 6.5244 7.5216 9.1275 8.9313 8.1041 7.7645 8.7056 5.9267 10.0860 12.0441 7.0129 8.4333

6.6693 6.4525 7.4479 8.8095 8.5625 7.6059 7.3350 8.2284 5.4590 9.3524 11.0350 6.2434 7.7468

6.6693 6.3775 7.3681 8.6082 8.1898 7.0994 6.8737 7.5815 5.1222 8.5850 10.0950 5.6131 7.0433

6.7870 6.5244 7.5216 9.2254 8.4466 7.9298 7.7645 8.7056 5.2775 9.2944 11.2526 6.1893 7.9406

6.7870 6.4525 7.4479 9.0057 8.0780 7.4299 7.3350 8.2284 5.0236 8.6048 10.2894 5.5486 7.2520

6.7870 6.3775 7.4234 8.8049 7.7610 6.9219 6.8737 7.5815 4.6136 7.9420 9.5392 4.9826 6.6742

SJZ–XC XC–XHM XHM–JG JG–YZQ YZQ–WHY WHY–WY WY–RJ RJ–RC RC–TJ TJ–JH JH–CQN CQN–CQ CQ–YZ Total

10.5099 7.3760 10.8283 11.5280 6.4479 9.1292 7.4857 8.2123 8.7448 8.9211 9.7614 6.0210 6.3260 218.1429

9.6850 6.6742 9.8069 10.5150 5.9549 8.4098 6.9683 7.9090 8.6575 9.0053 9.8459 6.0212 6.3260 206.7265

9.2464 6.4243 9.2731 9.9442 5.6383 8.0355 6.7325 7.6776 8.4595 8.8755 9.6944 6.0776 6.3260 197.6307

9.8643 6.7029 10.10 0 0 10.7514 5.8863 8.6220 6.9832 8.0975 8.7448 9.1642 9.9821 6.3967 6.6901 210.8447

9.1515 6.1335 9.1236 9.8345 5.3931 7.9568 6.5895 7.7942 8.6575 9.1020 9.9776 6.3969 6.6901 200.2836

8.7129 5.7462 8.6374 9.4072 5.1455 7.5244 6.3537 7.4472 8.4595 8.9724 9.9151 6.3274 6.6901 191.6212

(a) Tractive energy consumption for each inter-station during peak hours

15 The determistic model The stochastic model

Decreased

Tractive energy consumption (kW.h)

10

Increased

5

(b) Tractive energy consumption for each inter-station during off-peak hours

15 The determistic model The stochastic model 10

5

(c) Tractive energy consumption for each inter-station during night hours

15 The determistic model The stochastic model 10

5 J J J J Z Q Y C Y Y Z C Y Q C C G Q Q N H N H G M -C CQ N-J H-T J-R C-R J-W WH -YZ Q-J XHM -X -SJ Z-X X H M-J -YZ WH -W Y-R J-R C-T J-J CQ N-C Q-Y J T HM XC SJ - HY W T R R R Z Y YZ Q - CQ R C G H Q Q C H Y Y J J C X C JG X X W WH W YZ Fig. 18. Comparisons on tractive energy consumption between the deterministic and stochastic models.

444

X. Yang et al. / Transportation Research Part B 91 (2016) 424–445

7. Conclusion The main contribution of this paper is to develop a two-phase stochastic programming model for the integrated timetable and speed profile optimization problem that explicitly consider the uncertain train mass and the variable tractive force, braking force and basic running resistance. Moreover, we analyze the convergence property and formulation complexity of the developed model. The key points can be summarized as follows: (1) the dwell time at each station and the running time at each inter-station are the same during each of the three periods; (2) the frequency of trains is increased at peak hours and the headway is correspondingly reduced; (3) the train mass is a stochastic variable with a different probability distribution during each of the three different periods; and (4) the optimal journey for each given train mass at each interstation run is determined by a genetic algorithm and the expected energy consumption at each inter-station run in each period is determined by averaging over the entire probability mass distribution for that period. Based on the real-world operation data from the Beijing Metro Yizhuang Line, two numerical examples were conducted to illustrate the effectiveness of the proposed model. The results show that, during peak hours, off-peak hours and night hours, the total tractive energy consumptions can be reduced by: (1) 10.66%, 9.94% and 9.13% in comparison with the current timetable and speed profile; and (2) 3.35%, 3.12% and 3.04% in comparison with the deterministic model. One of our future research directions is to minimize the total net energy consumption instead of the total tractive energy consumption by considering the regenerative braking energy. The conversion, transmission and utilization of regenerative braking energy among multiple trains should be calculated accurately by modeling the whole electrical network of a metro system. Fast analytic algorithms are also worthy for future study to determine the optimal solution. In addition, we plan to test our numerical results on the real-world operations in the Beijing Metro Yizhuang Line. Based on the tested results, we will further adjust the model to obtain better achievements on energy conservation and efficiency. Acknowledgment The authors are grateful to three anonymous referees and Professor Hai Yang for their constructive comments and suggestions to improve the quality and clarity of the paper. This work was supported by the China National Funds for Distinguished Young Scientists (No. 71525002), the National Natural Science Foundation of China (Nos. 71210 0 01 and 61403020), the Research Foundation of State Key Laboratory of Rail Traffic Control and Safety (No. RCS2015ZZ003), the Beijing Laboratory of Urban Rail Transit, the Beijing Key Laboratory of Urban Rail Transit Automation and Control, the Hong Kong Polytechnic University, and the Chang Jiang Chair Professorship funded by the Ministry of Education in China. References Albrecht, A., Howlett, P., Pudney, P., Vu, X., 2013. Energy-efficient train control: from local convexity to global optimization and uniqueness. Automatica 49 (10), 3072–3078. Albrecht, A., Howlett, P., Pudney, P., Vu, X., Zhou, P., 2015a. The key principles of optimal train control-part 1: formulation of the model, strategies of optimal type, evolutionary lines, location of optimal switching points. Transp. Res. Part B Methodol. doi:10.1016/j.trb.2015.07.023. Albrecht, A., Howlett, P., Pudney, P., Vu, X., Zhou, P., 2015b. The key principles of optimal train control-part 2: existence of an optimal strategy, the local energy minimization principle, uniqueness, computational techniques. Transp. Res. Part B. Methodol. doi:10.1016/j.trb.2015.07.024. Bagloee, S.A., Ceder, A.A., 2011. Transit-network design methodology for actual-size road networks. Transp. Res. Part B Methodol. 45 (10), 1787–1804. ´ N., Goverde, R.M., Quaglietta, E., Roberti, R., 2016. An integrated micro-macro approach to robust railway timetabling. Transp. Res. Part B Bešinovic, Methodol. 87, 14–32. Birge, J.R., Louveaux, F., 2011. Introduction to Stochastic Programming. Springer Science & Business Media, New York. Bocharnikov, Y.V., Tobias, A.M., Robe, C., 2010. Reduction of train and net energy consumption using genetic algorithms for trajectory optimisation. In: Proceedings of IET Conference on Railway Traction Systems. IEEE Publisher, Birmingham, UK, pp. 32–36. Cacchiani, V., Caprara, A., Toth, P., 2010. Scheduling extra freight trains on railway networks. Transp. Res. Part B Methodol. 44 (2), 215–231. Chen, A., Kim, J., Lee, S., Kim, Y., 2010. Stochastic multi-objective models for network design problem under demand uncertainty. Expert Syst. Appl. 37 (2), 1608–1619. Chen, A., Subprasom, K., 2007. Analysis of regulation and policy of private toll roads in a build-operate-transfer scheme under demand uncertainty. Transp. Res. Part A Policy Pract. 41 (6), 537–558. Dai, L., Chen, C.H., Birge, J.R., 20 0 0. Convergence properties of two-stage stochastic programming. J. Optim. Theory Appl. 106 (3), 489–509. Ding, Y., Liu, H.D., Bai, Y., Zhou, F.M., 2011. A two-level optimization model and algorithm for energy-efficient urban train operation. J. Transp. Syst. Eng. Inf. Technol. 11 (1), 96–101. Fernández-Rodríguez, A., Fernández-Cardador, A., Cucala, A.P., Domínguez, M., Gonsalves, T., 2015. Design of robust and energy-efficient ATO speed profiles of metropolitan lines considering train load variations and delays. IEEE Trans. Intell. Transp. Syst. 16 (4), 2061–2071. Ghoseiri, K., Szidarovszky, F., Asgharpour, M.J., 2004. A multi-objective train scheduling model and solution. Transp. Res. Part B Methodol. 38 (10), 927–952. Gong, C., Zhang, S., Zhang, F., Jiang, J., Wang, X., 2014. An integrated energy-efficient operation methodology for metro systems based on a real case of shanghai metro line one. Energies 7 (11), 7305–7329. Howlett, P., 1988. Existence of an optimal strategy for the control of a train. School of Mathematics Report 3, University of South Australia. Howlett, P., 1990. An optimal strategy for the control of a train. J. Aust. Math. Soc. Ser. B 31, 454–471. Howlett, P., 20 0 0. The optimal control of a train. Ann. Oper. Res. 98 (1–4), 65–87. Howlett, P., Pudney, P., 1995. Energy-Efficient Train Control. Springer, London. Howlett, P., Pudney, P., Vu, X., 2009. Local energy minimization in optimal train control. Automatica 45 (11), 2692–2698. Ishikawa, K., 1968. Application of optimization theory for bounded state variable problems to the operation of trains. Bull. JSME 11 (47), 857–865. Kaniovski, Y.M., King, A.J., Wets, R.J., 1995. Probabilistic bounds (via large deviations) for the solutions of stochastic programming problems. Ann. Oper. Res. 56 (1), 189–208. Ke, B.R., Chen, M.C., Lin, C.L., 2009. Block-layout design using MAX-MIN ant system for saving energy on mass rapid transit systems. IEEE Trans. Intell. Transp. Syst. 10 (2), 226–235. Ke, B.R., Lin, C.L., Yang, C.C., 2012. Optimisation of train energy-efficient operation for mass rapid transit systems. IET Intell. Transp. Syst. 6 (1), 58–66. Khmelnitsky, E., 20 0 0. On an optimal control problem of train operation. IEEE Trans. Automat. Control 45 (7), 1257–1266.

X. Yang et al. / Transportation Research Part B 91 (2016) 424–445

445

Kim, M., Schonfeld, P., 2014. Integration of conventional and flexible bus services with timed transfers. Transp. Res. Part B Methodol. 16, 76–97. King, A.J., Rockafellar, R.T., 1993. Asymptotic theory for solutions in statistical estimation and stochastic programming. Math. Oper. Re. 18 (1), 148–162. Lee, D.H., Milroy, I.P., Tyler, K., 1982. Application of pontryagin’s maximum principle to the semi-automatic control of rail vehicles. In: The Second Conference on Control Engineering, Newcastle, Australia, pp. 233–236. Li, X., Lo, H.K., 2014a. An energy-efficient scheduling and speed control approach for metro rail operations. Transp. Res. Part B: Methodol. 64, 73–89. Li, X., Lo, H.K., 2014b. Energy minimization in dynamic train scheduling and control for metro rail operations. Transp. Res. Part B: Methodol. 70, 269–284. Li, X., Yang, X., 2013. A stochastic timetable optimization model in subway systems. Int. J. Uncertainty Fuzziness Knowl. Based Syst. 21 (1), 1–15. Liu, R.R., Golovitcher, I.M., 2003. Energy-efficient operation of rail vehicles. Transp. Res. Part A Policy Pract. 37 (10), 917–932. Lu, S., Hillmansen, S., Ho, T.K., Roberts, C., 2013. Single-train trajectory optimization. IEEE Trans. Intell. Transp. Syst. 14 (2), 743–750. Meng, L., Zhou, X., 2014. Simultaneous train rerouting and rescheduling on an n-track network: a model reformulation with network-based cumulative flow variables. Transp. Res. Part B Methodol. 67, 208–234. Milroy, I.P., 1981. Minimum-energy control of rail vehicles. In: The Railway Engineering Conference, Sydney, Australia, pp. 103–114. Nasri, A., Fekri Moghadam, M., Mokhtari, H., 2010. Timetable optimization for maximum usage of regenerative energy of braking in electrical railway systems. In: International Symposium on Power Electronics, Electrical Drives, Automation and Motion, Pisa, Italy, pp. 1218–1221. Niu, H., Zhou, X., 2013. Optimizing urban rail timetable under time-dependent demand and oversaturated conditions. Transp. Res. Part C Emerg. Technol. 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. Transp. Res. Part B Methodol. 76, 117–135. Online News, 2015. Survey on the development of the urban rail transit in China in 2014. (in Chinese). [online] Available at: http://www.tig-energy.com/ newsview.asp?classid=10&id=202. Peña Alcaraz, M., Fernández-Cardador, A., Cucala, A.P., Ramos, A., Pecharromán, R.R., 2012. Optimal underground timetable design based on power flow for maximizing the use of regenerative-braking energy. Proc. Ins. Mech. Eng. Part F J. Rail Rapid Transit 226 (4), 397–408. Ramos, A., Pena, M., Fernández-Cardador, A., Cucala, A.P., 2007. Mathematical programming approach to underground timetabling problem for maximizing time synchronization. In: International Conference on Industrial Engineering and Industrial Management, Madrid, Spain, pp. 88–95. Shiryayev, A.N., 1984. Probability. Springer Verlag, New York. Su, S., Li, X., Tang, T., Gao, Z., 2013. A subway train timetable optimization approach based on energy-efficient operation strategy. IEEE Trans. Intell. Transp. Syst. 14 (2), 883–893. Tuyttens, D., Fei, H.Y., Mezmaz, M., Jalwan, J., 2013. Simulation-based genetic algorithm towards an energy-efficient railway traffic control. Math. Probl. Eng. Article ID 805410, 1–12. Xiao, X., 2014. The Record of Vehicle Condition for Beijing Metro Yizhuang Line. Beijing Mass Transit Railway Operation Corporation Limited Technical report. (in Chinese). Xu, X., Li, K., Yang, L., 2015. Scheduling heterogeneous train traffic on double tracks with efficient dispatching rules. Transp. Res. Part B Methodol. 78, 364–384. Xun, J., Cao, F., Yang, X., 2014. Statistics of the Numbers of Passengers on the Train at Each Inter-Station During Different Operational Periods for the Beijing Metro Yizhuang Line. Beijing Jiaotong University Technical report. (in Chinese). 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. Transp. Res. Part C Emerg. Technol. 57, 13–29. Yang, X., Li, X., Gao, Z., Wang, H., Tang, T., 2013. A cooperative scheduling model for timetable optimization in subway systems. IEEE Trans. Intell. Transp. Syst. 14 (1), 438–447. Yang, X., Li, X., Ning, B., Tang, T., 2015b. An optimization method for train scheduling with minimum energy consumption and travel time in metro rail systems. Transportmetrica B Transp. Dynamics 3 (2), 79–98. Yang, X., Li, X., Ning, B., Tang, T., 2016. A survey on energy-efficient train operation for urban rail transit. IEEE Trans. Intell. Transp. Syst. 17 (1), 2–13. Yang, X., Ning, B., Li, X., Tang, T., 2014. A two-objective timetable optimization model in subway systems. IEEE Trans. Intell. Transp. Syst. 15 (5), 1913–1921. Yu, C., Ma, J., Lo, H.K., Yang, X., 2015. Optimization of mid-block pedestrian crossing network with discrete demands. Transp. Res. Part B Methodol. 73, 103–121. Zhang, L., 2013. Energy Consumption for the Beijing Metro in 2012. Beijing Mass Transit Railway Operation Corporation Limited Technical report. (in Chinese). Zhang, L., 2014. The Operation Data for the Beijing Metro Yizhuang Line. Beijing Mass Transit Railway Operation Corporation Limited Technical report. (in Chinese). Zhou, X., Zhong, M., 2007. Single-track train timetabling with guaranteed optimality: branch-and-bound algorithms with enhanced lower bounds. Transp. Res. Part B Methodol. 41 (3), 320–341.