A bi-objective timetable optimization model incorporating energy allocation and passenger assignment in an energy-regenerative metro system

Transportation Research Part B 133 (2020) 85–113

Contents lists available at ScienceDirect

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

A bi-objective timetable optimization model incorporating energy allocation and passenger assignment in an energy-regenerative metro system Songpo Yang a,b,c, Feixiong Liao c,∗, Jianjun Wu a,∗, Harry J.P. Timmermans c,d, Huijun Sun e, Ziyou Gao e a

State Key Laboratory of Rail Traffic Control and Safety, Beijing Jiaotong University, Beijing, China Beijing Key Laboratory of Traffic Engineering, Beijing University of Technology, Beijing, China c Urban Planning Group, Eindhoven University of Technology, Eindhoven, the Netherlands d Department of Air Transportation Management, NUAA, Nanjing, China e School of Traffic and Transportation, Beijing Jiaotong University, Beijing, China b

a r t i c l e

i n f o

Article history: Received 9 October 2019 Revised 25 December 2019 Accepted 2 January 2020

Keywords: Passenger assignment Energy allocation Irregular timetable Block operation Local search

a b s t r a c t Complex passenger demand and electricity transmission processes in metro systems cause difficulties in formulating optimal timetables and train speed profiles, often leading to inefficiency in energy consumption and passenger service. Based on energy-regenerative technologies and smart-card data, this study formulates an optimization model incorporating energy allocation and passenger assignment to balance energy use and passenger travel time. The Non-Dominated Sorting Genetic Algorithm II (NSGA-II) is applied and the core components are redesigned to obtain an efficient Pareto frontier of irregular timetables for maximizing the use of regenerative energy and minimizing total travel time. Particularly, a parallelogram-based method is developed to generate random feasible timetables; crossover and local-search-driven mutation operators are proposed relying on the graphic representations of the domain knowledge. The suggested approach is illustrated using realworld data of a bi-directional metro line in Beijing. The results show that the approach significantly improves regenerative energy use and reduces total travel time compared to the fixed regular timetable. © 2020 Elsevier Ltd. All rights reserved.

1. Introduction Metro systems play an important role in urban economies and the social development of big cities due to their low fares, high loading capacity, punctuality, and low energy consumption. The metro is also seen as an environmental-friendly transport mode that has a large potential for alleviating traffic congestion and emissions. With these advantages, metro systems have expanded rapidly in big cities all over the world (Yang et al., 2015). In Beijing (China), for example, to date, 22 lines have been built to connect the primary points of interest. The metro system serves more than 12 million commuters on an average weekday. Besides the high travel time reliability, it is estimated that taking the metro rather than private car produces around 90% less emission for an average commuting trip during peak hours. ∗

Corresponding authors. E-mail addresses: [email protected] (S. Yang), [email protected] (F. Liao), [email protected] (J. Wu).

https://doi.org/10.1016/j.trb.2020.01.001 0191-2615/© 2020 Elsevier Ltd. All rights reserved.

86

S. Yang, F. Liao and J. Wu et al. / Transportation Research Part B 133 (2020) 85–113

In general, a well-designed metro system has the infrastructure and timetables developed in a passenger-oriented and eco-friendly manner. In particular, train scheduling aims at generating timetables for all service trains, specifying their arrival and departure times at each station of the metro lines. Regular timetables are commonly used by the current metro systems, in which the headway times, running times and dwelling times are all constant. In contrast, irregular timetables contain certain varied time components. The designed timetable affects train energy consumption and passenger travel time. Passenger travel time is closely related to service quality, while energy consumption, which is mainly for traction, has negative economic and environmental effects. In early studies (e.g., Howlett., 1996), traction energy is determined by the running phases and speed profiles of independent trains. In an energy-regenerative system (Yang et al., 2015, 2019a,b), energy generated by one decelerating train can be used by one accelerating train. Considering the above two objectives simultaneously is challenging for the following reasons. First, the relationship between passenger flow characteristics and travel times is complex (Shang et al., 2019). The majority of prior studies (e.g., Yang et al., 2013; Niu and Zhou, 2013; Li and Lo 2014b; Luan et al., 2018) are based on passenger characteristics obtained from extrapolation or surveys. The challenge is to develop a model realizing that passenger travel time is subject to passenger flow characteristics and that passenger assignment should deal with time-variant features. Second, the optimization of energy consumption for train operations is an inherently complex nonlinear problem (Yin et al., 2017). Due to the integrity of electricity transmission in a metro system, electricity supplies should be considered as a whole. Third, these two objectives are often in conflict with each other. For example, higher frequency and shorter running times reduce passenger travel time but increase energy consumption. Therefore, the trade-off between the two objectives should be considered when designing timetables. Energy-saving timetable optimization of metro systems has attracted much attention recently in the operations and transportation research communities (e.g., Yang et al., 2013; Li and Lo, 2014a; Albrecht et al., 2016a,b; Gupta et al., 2016; Howlett, 2016; Ye and Liu, 2016; Yang et al., 2019a). Several studies have focused on minimizing traction energy during a train operating period, which is usually composed of three phases: acceleration, coasting, and deceleration. For example, Howlett (1996) first addressed the running time on a single section, which has an influence on traction energy consumption. Albrecht et al. (2016a) developed an optimal control theory for a single train on a section with a long-fixed distance using continuous control approaches, which suggest driving strategies to minimize the total energy consumption and ensure on-time arrivals. Ye and Liu (2016) considered an optimization problem of train control and scheduling. The optimal control problem involved complex factors such as undulating track, variable speed limitations, running resistances, and maximum tractive/braking forces. The numerical results showed that the designed control strategy could reduce energy consumption by 0.27% compared with the current timetable. Furthermore, Albrecht et al. (2018) found the necessary conditions for minimizing the total traction energy of two trains on a non-level track, subject to the safety-compatible constraints, and optimized the constraints using Pontryagin principle (Hartl et al., 1995) to obtain the sequences of driving speeds of each train. As a result, they calculated the speed profiles separately since the two trains are logically independent in the problem. With the equipment of energy-regenerative devices, traction energy can be further reduced by coordinating the arrival and departure times of multiple trains. According to the three-working-phase mode, a few studies formulated train optimal control strategies to improve the use of regenerative energy involving variable speed limits, variable track gradients, and uncertain dwelling times. For example, Li and Lo (2014b) first formulated an integrated energy-efficient operational model to design the timetables and speed profiles to reduce net energy consumption. In their study, the authors divided a running section into several segments with constant speed limits. Then, they assumed that speed profiles consist of accelerationcruising-braking phases on a single segment according to the analysis of the optimal phase sequences. Yang et al. (2015) attempted to coordinate the arrival and departure times of trains, which are located in the same independent electricity supply interval for reducing energy use. Gupta et al. (2016) proposed a two-step linear optimization model to obtain a timetable aiming at minimizing the total energy consumption of a metro network. Yang et al. (2018a) formulated an integer programming model to minimize energy consumption. However, the above optimization models may not be able to achieve the desired effects as only one objective has been considered. In view of the importance of the multiple objectives, there are some studies considering passenger waiting times and system energy consumption simultaneously in bi-objective optimization models. For example, Ning et al. (2015) proposed an integrated model to optimize headways by adjusting the arrival times at stations to decrease the delay times and average waiting times. A dynamic programming algorithm was developed to calculate the switching points of the speed profile in an operating section. The results show that the average passenger waiting time is reduced by 23% and the energy consumption is reduced by 16%. Niu et al. (2015) developed a unified integer quadratic programming model with linear constraints. A number of mathematically tractable forms were used to characterize the total waiting times using general optimization models, which were implemented in a general algebraic modeling system in a corridor. Yin et al. (2017) proposed an integrated model to minimize energy consumption and passenger waiting time in a bi-directional metro line. An algorithm based on Lagrangian relaxation was suggested to find a good solution after employing the space-time network and combining the two objectives into one. Yang et al. (2018b) developed a bi-objective nonlinear programming model to determine the optimal timetables and speed profiles considering minimum energy consumption and passenger waiting time. A quadratic programming method was developed to calculate the switching points of train speed profiles. Sun et al. (2019) considered time-variant passenger demands at stations for minimizing passenger waiting time and energy consumption, for which a genetic algorithm was developed to find good solutions. In addition, Yang et al. (2016) and Scheepmaker et al. (2017)

S. Yang, F. Liao and J. Wu et al. / Transportation Research Part B 133 (2020) 85–113

87

Table 1 The main differences from the relevant studies.

n

Related objective

Yang et al. (2013)

Energy consumption Passenger waiting time Energy consumption Passenger travel time Energy consumption Energy consumption Energy consumption Passenger waiting time Energy consumption Passenger waiting time Energy consumption Passenger waiting time Energy consumption Passenger travel time

Wang et al. (2015) Yang et al. (2015) Gupta et al. (2016) Yin et al. (2017) Yang et al. (2018a) Sun et al. (2019) This study

Energy allocation feature

Solution outcome

Solution algorithm

One-to-one

Single solution

GA

One-to-one

Single solution

One-to-one One-to-one One-to-one

Single solution Single solution Single solution

One-to-one

Single solution

Iterative convex programming GA Optimization solver Lagrangian-based algorithm Lingo software

One-to-one

Single solution

GA

One-to-many

Pareto frontier

NSGA-II

conducted literature reviews that focused on energy consumption optimization for timetabling, speed profile generation, and the integration of both. Nevertheless, since passenger waiting time is only a portion of travel time, it may not be suitable to have one main objective for assessing the effectiveness of the timetables. Wang et al. (2015) proposed a real-time scheduling algorithm with the aim of minimizing total passenger travel time and energy consumption. The optimization problem was transformed into a real-valued nonlinear nonconvex problem. Subsequently, they developed an iterative convex programming approach and compared the performance with two heuristic algorithms. The results showed that the performance of sequential quadratic programming is better than the genetic algorithm in irregular timetable optimization, but its computational time is higher even when processing regular timetables. To summarize, the above studies have three limitations. First, most studies treated the two objectives simply using the weighted average method, which nonetheless requires perfect knowledge of the coefficients. Second, passenger assignment is not considered at a detailed level that allows calculating the proportions of passengers at different travel stages (i.e., waiting, on-train, and walking out). Third, a comprehensive energy allocation mechanism is lacking to capture the energy transmission patterns among multiple accelerating and decelerating trains. As an extension of the above bi-objective models, we develop a bi-objective optimization model incorporating passenger assignment and energy allocation. The Non-Dominated Sorting Genetic Algorithm II (NSGA-II) is applied to obtain an efficient Pareto frontier of irregular timetables for maximizing the use of regenerative energy and minimizing total travel time. Using the principle of parallelograms that the lengths of parallel edges are equal, we randomly generate a set of irregular timetables in strict accordance with the constraints. The crossover and mutation operators of the NSGA-II are designed based on the structure of the timetable solutions. The model is applied to a bi-directional metro line with the input of smart-card data. This study contributes to the state-of-the-art energy optimization of timetables and speed profiles with the consideration of energy regenerative technologies. First, using historical smart-card data, we formulate the train- and timetable-dependent OD tables. With the limited recorded passenger information, we develop an effective passenger assignment rule to calculate the travel times at different travel stages according to a timetable. Second, for calculating the use of regenerative energy, we formulate the dynamic energy allocation mechanism incorporating time, speed, and position dependencies. Third, to evaluate the bi-objective functions simultaneously, we decompose the decelerating phases using the multi-state space-time diagram (Zhou and Zhong, 2005; Zhang et al., 2011; Liao et al., 2010, 2011, 2013, 2017; Liao, 2016; Lu et al., 2018; Tong et al., 2019). Fourth, we redesign the core components of the NSGA-II algorithm to obtain an efficient Pareto frontier of irregular timetables for maximizing regenerative energy use and minimizing total travel time. Particularly, a parallelogram-based method is developed to generate random feasible timetables; crossover and local-search-driven mutation operators are proposed relying on the graphic representations of the domain knowledge. The main differences between the relevant studies and this paper are listed in Table 1. The remainder of this paper is organized as follows. Section 2 presents the problem description and basic considerations. In Section 3, the bi-objective optimization model is proposed for improving the regenerative energy use and reducing passenger travel time. In Section 4, we discuss an elaborated NSGA-II algorithm to solve the bi-objective optimization model. Section 5 presents the numerical experiments and discusses the effectiveness of the proposed algorithm for obtaining a quality Pareto frontier. Finally, conclusions and possibilities for further research are discussed in Section 6. 2. Problem formulation A metro line is considered as a closed-loop (Fig. 1(a)), which consists of stations as the nodes and sections between two adjacent stations as the arcs. There are N stations and N − 1 undirected arcs in Fig. 1(a). Without loss of generality, the trains

88

S. Yang, F. Liao and J. Wu et al. / Transportation Research Part B 133 (2020) 85–113

Fig. 1. Illustration of energy allocation in a bi-directional metro line.

are assumed to start from two depots located at the two ends of a metro line, indexed by terminus 1 and N respectively. Each train fleet has K trains starting from termini 1 and N in the upstream and downstream directions respectively. Within each train fleet, the train index k is ranged from 1 to K. After turning around, the trains run in the opposite direction and the turnaround index tr is updated as tr + 1. Particularly, to distinguish the train fleets, we introduce a matrix as the solution timetable. The odd rows represent the information of the train fleet starting from termini 1 at the beginning time (up-train fleet); the even rows represent the information of the train fleet starting from termini N at the beginning time (down-train fleet). In a train fleet, train k + 1 (k + 1 ≤ K) runs immediately after train k according to a timetable. Accounting for the bidirectional and turnaround properties, the design of an irregular timetable should consider the requirements of unbalanced passenger demands and circular operations. The substations built along the stations supply electrical power to ensure the normal operations of the trains. In a metro line with energy-regenerative devices, the moving trains interact through energy transmission. For example, in Fig. 1(b), the arrows of different colors represent the different moving phases of trains. The train-like icons of different colors run in opposite directions. The energy generated by the two decelerating trains in the downstream direction can be transferred to the two accelerating trains in the upstream direction, according to the status (time, speed, and position) of the trains. By coordinating the accelerating and decelerating trains within the maximum electricity transmission distance, the generated energy can be efficiently used. 2.1. Notation The following notation is used for formulating the proposed approach. General variables and parameters i, j Station index, i = 1, 2, …, N; j = 1, …, N. Passengers’ origin stations of the metro trips. io , jo k Train index, k = 1, 2, ..., K. Accelerating train index. k t A time instance that uses 1 s as the unit. τ A time instance that uses 0.1 s as the unit. tr Turnaround index, tr = 0, 1, ..., TR.  Decelerating phase and accelerating phase index, z, z = 1, 2, ..., z, z (TR + 1) • (N − 1). n Element index, n = 1, 2, ..., (TR+1) • (2N − 1). a Acceleration rate during the acceleration phase. Deceleration rate and substitute deceleration rate in mutation operator. b , b r Resistance rate. M Train mass. θ 1, θ 2 Conversion factors from electricity to kinetic energy and from kinetic energy to electricity. Length (distance) between i and i + 1. li,i + 1 (continued on next page)

S. Yang, F. Liao and J. Wu et al. / Transportation Research Part B 133 (2020) 85–113

T¯ DW , T DW T¯ , T H¯ , H T¯ TR , T TR TTR C TW λi ( t )

ωi,j

The upper and lower bounds of dwelling time at a station. The upper and lower bounds of running time on a section. The upper and lower bounds of headway between two consecutive trains. The upper and lower bounds of turnaround time at termini. Turnaround time. Train loading capacity. The walking-out time per traveler at a station. Arrival rate of passengers at i at time t. Proportion of passengers from i to j.

Intermediate variables EkA ,z

Traction energy of k in z

B,U B,D , Ek,z Ek,z

Regenerative energy by k of the up- and down-train fleet in z.

U,U U,D , Ek,z Ek,z

Use of energy regenerated by k of the up- and down-train fleet in z.

W,U B,U , Ptr,k,i, Ptr,k,i, j j W,U B,U Ptr,k,i , Ptr,k,i

The numbers of waiting and boarded passengers for k from i to j at times. The total numbers of waiting and boarded passengers for k at i.

L,U R,U , Ptr,k,i Ptr,k,i

The total numbers of alighted and stranded passengers for k at i.

B E , tk,z tk,z

The beginning and ending times of k in z.

Decision variables A D , ttr,k,i,i ttr,k,i,i +1 +1 A D , ttr,k,i,i ttr,k,i,i −1 −1 A D , ttr,k,N,N ttr,k,N,N A D , ttr,k, ttr,k, 1,1 1,1 WA WD , ttr,k,i,i ttr,k,i,i +1 +1

WA WD , ttr,k,i,i ttr,k,i,i −1 −1

89

Arrival and departure times of k at i after turning around tr times in the upstream direction, respectively. Arrival and departure times of k at i after turning around tr times in the downstream direction, respectively. Arrival and departure times of k at terminus N in the upstream direction, respectively. Arrival and departure times of k at terminus 1 in the downstream direction, respectively. Switching time points from accelerating to coasting and from coasting to decelerating of k after turning around tr times in the upstream direction, respectively. Switching time points from accelerating to coasting and from coasting to decelerating of k after turning around tr times in the downstream direction, respectively.

2.2. Basic considerations The following assumptions are made to facilitate the presentation of the modeling framework. (A1) The study focuses on the timetabling problem of a metro line. According to the metro smart-card data, the arrival rate of passengers is constant in a time interval and the time-dependent OD matrix is also constant (Gao et al., 2016). (A2) Passengers can board a train until it has reached its maximum loading capacity. Non-served passengers incur additional waiting times. (A3) There is no transfer to other metro lines. This assumption, albeit unlikely in big cities, is widely adopted in the existing studies of energy-efficient timetable optimization (e.g., Yang et al., 2013; Li and Lo 2014b; Niu et al., 2015; Yin et al., 2017; Yang et al., 2018a; Sun et al., 2019). (A4) Train running times and speed profiles may vary in the different sections and involve three phases in each section: acceleration-coasting-deceleration (Yang et al., 2015). (A5) Trains consume energy in the form of electricity. To simplify the process of calculating energy consumption, we assume that the conversion factors are constants (Yang et al., 2018a). 2.3. The bi-objective function In this subsection, we discuss the formulations of the two conflicting objectives: minimizing total passenger travel time and maximizing the use of regenerative energy. Less travel time contributes to improving the service quality. A higher utilization rate usually leads to less energy consumption. The calculation of the utilization rate of regenerative energy considers the coordination of multiple trains in a timetable, while the calculation of the passenger travel time considers passenger assignment. The first objective is to minimize total travel time. The travel time of a metro trip includes multiple time components. In this paper, a metro trip is recorded from the time arriving at the platform of an origin station to the time checking out of the destination station. We divide the total travel time into four parts: waiting time (WT), on-train time (OT), extra waiting time (EWT), and walking-out time (WOT). The walking-in time is not considered since it is fixed given the same passenger demands. The walking-out time is considered since the passenger demands at their destination stations always

90

S. Yang, F. Liao and J. Wu et al. / Transportation Research Part B 133 (2020) 85–113

change according to the irregular timetables and train capacity. They are closely related to the calculation of the stranded passengers for each train. Overall, the first objective is expressed as

min [W T +OT + EW T + W OT ]

(1)

The second objective is to maximize the use of regenerative energy. The reasons for this consideration are listed below. In this study, we have mainly focused on the regenerative energy flows in the space-time diagram. It is straightforward to calculate the utilization rate from the modeling perspective, given that utilization rate and energy consumption are strongly positively correlated. However, it is not to say that formulating the total energy consumption is not possible. From managerial and technological points of view, the utilization rate is important. Having the minimum total energy consumption may not necessarily be the sole operational objective. For example, here are two timetables resulting in traction energy consumption of 100 J and 90 J and use of regenerative energy of 50 J and 40 J, respectively. The total energy consumptions are both 50 J, but the utilization rates are 50% and 44.44% respectively. Thus, the second objective is to maximize the average utilization rate u of the regenerative energy:

max [u]

(2)

3. Optimization model In this section, we formulate the two objectives of the irregular timetabling problem. First, the basic train timetable constraints are listed. Second, the numbers and the travel times of passengers at each stage (i.e., waiting, on-train, stranded, and walking-out) are calculated. Then, we define the situations of energy exchange and specify an energy allocation mechanism when multiple accelerating trains share the regenerative energy. Finally, the complete bi-objective optimization model is presented. 3.1. Train timetable constraints With the cyclic and irregular properties of a timetable, it is necessary to distinguish the different running directions. Specifically, the upstream direction is defined as the trains running from station i to i + 1(1 ≤ i < i + 1 ≤ N); likewise, the downstream direction is defined as the trains running from i to i − 1(1 ≤ i − 1 < i ≤ N). The basic constraints related to a timetable are given as follows.



A A H ≤ ttr,k − ttr,k,i,i ≤ H¯ +1,i,i+1 +1 , i = N , 1 ≤ k ≤ K − 1, 0 ≤ tr ≤ T R D D H≤t −t ≤ H¯

(3.1)

A A H ≤ ttr,k − ttr,k,i,i ≤ H¯ +1,i,i−1 −1 , i = 1, 1 ≤ k ≤ K − 1, 0 ≤ tr ≤ T R D D H≤t −t ≤ H¯

(3.2)

A A ttr+1 ,1,1,2 − ttr,K,1,2 ≥ H , 0 ≤ tr < tr + 1 ≤ T R D D ttr+1,1,1,2 − ttr,K, 1,2 ≥ H

(4.1)

A A ttr+1 ,1,N,N−1 − ttr,K,N,N−1 ≥ H , 0 ≤ tr < tr + 1 ≤ T R D D ttr+1,1,N,N−1 − ttr,K,N,N−1 ≥H

(4.2)

A D ttr,k − ttr,k,i,i > 0, i = N +1,i,i+1 +1 , 1 ≤ k ≤ K − 1, 0 ≤ tr ≤ T R A D ttr,k − t > 0 +1,N,N tr,k,N,N

(5.1)

A D ttr,k − ttr,k,i,i > 0, i = 1 +1,i,i−1 −1 , 1 ≤ k ≤ K − 1, 0 ≤ tr ≤ T R A D ttr,k − t > 0 +1,1,1 tr,k,1,1

(5.2)



tr,k+1,i,i+1



tr,k+1,i,i−1

tr,k,i,i+1

tr,k,i,i−1

  

⎧ DW D A DW ⎨T ≤ ttr,k,i,i+1 − ttr,k,i,i+1 ≤ T¯ , i = N

D A 1 ≤ k ≤ K , 0 ≤ tr ≤ T R T DW ≤ ttr,k,i,i − ttr,k,i,i ≤ T¯ DW , i = 1, −1 −1 DW T ≤ tD − tA ≤ T¯ DW , i = N or i = 1

⎩ 

T≤ T≤



tr,k,i,i

A ttr,k,i +1,i+2 A ttr,k,i −1,i−2

tr,k,i,i

D − ttr,k,i,i ≤ T¯ , i ≤ N − 2 +1 , 1 ≤ k ≤ K , 0 ≤ tr ≤ T R D − ttr,k,i,i−1 ≤ T¯ , i ≥ 3

A D T ≤ ttr,k,N,N − ttr,k,N−1 ≤ T¯ ,N , 1 ≤ k ≤ K , 0 ≤ tr ≤ T R A D ¯ T ≤t −t ≤T



tr,k,1,1

(6)

(7.1)

(7.2)

tr,k,2,1

A D T TR ≤ ttr+1 − ttr,k,N,N ≤ T¯ TR ,k,N,N−1 , 1 ≤ k ≤ K , 0 ≤ tr < tr + 1 ≤ T R TR A D T ≤t −t ≤ T¯ TR tr+1,k,1,2

tr,k,1,1

(8)

S. Yang, F. Liao and J. Wu et al. / Transportation Research Part B 133 (2020) 85–113

⎧ ⎨t WD



tr,k,i,i+1

WA − ttr,k,i,i = +1



2r·li,i+1 C

−

A C



2

A D · ttr,k,i − ttr,k,i,i , i≤N−2 +1,i+2 +1

, 1 ≤ k ≤ K , 0 ≤ tr ≤ T R A 2 2r·li,i−1 ⎩t WD A WA D − t = − · t − t , i ≥ 3 C C tr,k,i,i−1 tr,k,i,i−1 tr,k,i−1,i−2 tr,k,i,i−1 ⎧   2 2r·lN−1,N ⎨t WD WA A D − ttr,k,N−1 = − CA · ttr,k,N,N − ttr,k,N−1 C tr,k,N−1,N ,N ,N  , 1 ≤ k ≤ K , 0 ≤ tr ≤ T R 2 ⎩t WD − t WA = 2r·l2,1 − A · t A − tD tr,k,2,1

tr,k,2,1

C

tr,k,1,1

C

91

(9.1)

(9.2)

tr,k,2,1

Constraint (3) describes the headway constraints of arrival-arrival and departure-departure times to guarantee the safety A D of consecutive trains in a metro line, where ttr,k,i,i and ttr,k,i,i are the arrival and departure times of train k at i after k +1 +1 ¯ turns around tr(0 ≤ tr ≤ TR) times. H and H are the lower and upper bounds of headway times between every two conA A secutive trains. At terminus N in the upstream direction, the arrival-arrival headway time (ttr,k − ttr,k,N,N ) of k should +1,N,N A A be in the range [H , H¯ ]. At terminus 1 in the downstream direction, the arrival-arrival headway (ttr,k − t ) should +1,1,1 tr,k,1,1 ¯ also be in the range [H , H ]. Similarly, the departure-departure headway times at termini should also satisfy the headway

A constraints. Particularly, if tr = 0, ttr,k,i,i is the arrival time when train k is for the first time running in the upstream +1 A direction; ttr,k,i,i represents running in the downstream direction. To guarantee safety between the last train of one train −1 fleet and the first train in an ensuing fleet at the termini, the headways should be greater than the minimum headway time, A D shown as constraint (4), where ttr,K,N,N−1 and ttr,K,N,N−1 are the arrival and departure times of the last train K in one train

A D fleet at terminus N in the downstream direction. ttr+1 and ttr+1 respectively denote the arrival and departure ,1,N,N−1 ,1,N,N−1 times at N of the first train in the ensuing train fleet in the downstream direction. To avoid two trains being at the same station at the same time, the arrival time of a train at a station should be always later than the previous train’s departure time, shown as constraint (5). Constraint (6) describes the dwelling time constraints to guarantee the trains dwelling at stations within a time window [T DW , T¯ DW ]. Constraint (7) describes the running time constraints with a time window [T , T¯ ]. Constraint (8) describes the turnaround times of trains at termini 1 and N in different directions, which should fall within time window[T TR , T¯ TR ]. Eq. (9.1) calculates the switching points of a train that runs on a section in both directions, where WA WD ttr,k,i,i is the switching time point of k between i and i + 1 from accelerating to coasting, ttr,k,i,i is the switching time +1 +1 point of k from coasting to decelerating, li,i + 1 is the length of the section, r is the resistance rate, and A and C are the intermediate parameters (the detailed formulations are described in Appendix A). Eq. (9.2) calculates the switching points of a train that runs on the last sections in both directions.

3.2. Passenger travel time In this subsection, we formulate passenger-flow constraints related to a train timetable. The four components of the total travel time function in both running directions are calculated one by one. By definition above, the passengers travel from station i to j (j > i) in the upstream direction; the passengers travel from station j (j > i) to i in the downstream direction. (1) Passenger demand constraints In a metro system, the passenger tap-in and tap-out information is recorded by the Automatic Fare Collection (AFC) System. After discretizing the tap-in data at each station in time taking integer values, the time-dependent arrival rates are obtained. The dynamic passenger demands can be calculated based on the arrival rates and OD demands. The arrival rate λi (t), defined as the number of arrival passengers per second, records the demand at time t station i. The OD matrix describes the proportion of passengers arriving at the platforms of the origin stations to reach the destination stations. We use ωi,j (j > i) to denote the proportion of passengers traveling in the upstream direction, and correspondingly ωj,i (j > i) to denote the proportion of passengers traveling in the downstream direction. Thus, the numbers of passengers in each stage are shown as follows. W ,U Ptr,k,i, j

D ttr,k,i,i +1



=

λi (t ) · ωi, j ;

W ,D Ptr,k, j,i

D ttr,k,i,i −1

=

D t =ttr,k +1 −1,i,i+1

W ,U R,U Ptr,k,i = Ptr,k + −1,i

W ,U R,U Ptr, = Ptr−1 + 1,i ,K,i

W ,U Ptr, = 1,i

N j=i+1

N j=i+1 N j=i+1



λ j (t ) · ω j,i

(10)

D t =ttr,k +1 −1,i,i−1

W ,U W ,D R,D Ptr,k,i, , k = 1, tr = 0; Ptr,k, = Ptr,k + j j −1, j

W ,U W ,D R,D Ptr, , tr = 0; Ptr, = Ptr−1 + 1,i, j 1, j ,K, j

W ,U W ,U Ptr, , tr = 0; Ptr, = 1,i, j 1, j

j−1 i=1

W ,U Ptr, , tr = 0 1, j,i

j−1 i=1

j−1

W ,D Ptr,k, , k = 1, tr = 0 j,i

(11.1)

i=1

W ,D Ptr, , tr = 0 1, j,i

(11.2)

(11.3)

92

S. Yang, F. Liao and J. Wu et al. / Transportation Research Part B 133 (2020) 85–113 R,U W ,U B,U Ptr,k,i = Ptr,k,i − Ptr,k,i

(12)

B,U Ptr,k,i =

W ,U Ptr,k,i ,C

min



−

i−1

B,U Ptr,k,q

+

q=1

i

 L,U Ptr,k,q

, i = 1

(13.1)

q=2

B,U W Ptr,k,i = min Ptr,k,i ,C , i = 1 i−1

L,U Ptr,k,i =

io =1 B,U Ptr,k,i, = j

B,U L,D Ptr,k, ; Ptr,k, = i ,i j o

W ,U Ptr,k,i, j W ,U Ptr,k,i

(13.2) N jo = j+1

B,U B,D · Ptr,k,i ; Ptr,k, = j,i

B,D Ptr,k, j ,j

(14)

o

W ,D Ptr,k, j,i W ,D Ptr,k, j

B,D · Ptr,k, j

(15)

Eq. (10) calculates the numbers of waiting passengers at stations in both directions. In the upstream direction, the time D D D difference between two consecutive trains k − 1 and k is ttr,k,i,i − ttr,k . Similarly, the time difference is ttr,k,i,i − +1 −1,i,i+1 −1

D ttr,k in the downstream direction. One second is used as the time unit when calculating the passenger numbers at each −1,i,i−1

W,U D stage. The passengers are accumulated from time ttr,k + 1. Ptr,k,i, is the number of waiting passengers from station i to −1,i,i+1 j

j for train k after turning around tr times, and λi (t)

•

W,D ωi,j is the passenger arrival rate at i for j at time t;Ptr,k, denotes the j,i

W,U number of waiting passengers traveling in the downstream direction. In addition, the number of waiting passengers Ptr,k,i at

W,D R,U i in the upstream direction and Ptr,k, at j in the downstream direction can be calculated by Eq. (11), where Ptr,k is the j −1,i

number of stranded passengers from k − 1 at i. Due to the train loading capacity, some passengers may not be able to board a train immediately after they arrive at the platform. Based on Eq. (11), we can calculate the number of stranded passengers R,D who are left on the platform in the upstream direction by Eq. (12). Similarly, the number of stranded passengers Ptr,k, in j B,U the downstream direction can be calculated. Therefore, the number of boarded passengers Ptr,k,i at i for k in the upstream

L,U direction can be calculated by Eq. (13), where C is the train capacity and Ptr,k,i is the number of alighted passengers from B,D k at i. Similarly, the number of boarded passengers Ptr,k, in the downstream direction can be calculated. Eq. (14) calculates j B,U the number of alighted passengers at i in the upstream direction and j in the downstream direction, where Ptr,k, i B,D Ptr,k, j

o, j

o ,i

and

are the numbers of boarded passengers for k from origin station io (io < i) to i in the upstream direction and from

origin station jo (jo > j) to j in the downstream direction, respectively. Eq. (15) calculates the number of boarded passengers for each train and each OD pair according to a proportional function, which, however, cannot guarantee the first-in-first-out property in the metro system (Wang et al., 2015; Gao et al., 2016). (2) Travel time functions According to Eq. (1), we calculate the time expenses of the four components of the total travel time respectively, including waiting time (WT), on-train time (OT), walking-out time (WOT), and extra waiting time (EWT), where WT = WTU + WTD , OT = OTU + OTD ,WOT = WOTU + WOTD , and EWT = EWTU + EWTD with the consideration of both directions. The detailed functions in the upstream direction are given as follows.

⎛

WTU =

TR K tr=0 k=1

U

OT =

TR K tr=0 k=1

W OT U =

N−1

⎝



D ttr,k,i,i +1

N





i=1 j=i+1 t =t D



R,U D Ptr,k + λi (t ) · ωi, j · ttr,k,i,i +1 − t −1,i



⎞

⎠, P0R,,0U,1 = 0

+1 tr,k−1,i,i+1

N−2 N−1

B,U Ptr,k,i, j

·



A ttr,k, j, j+1

i=1 j=i+1

TR K N−1 N

D − ttr,k,i,i +1



+

N−1

B,U Ptr,k,i,N

·



A ttr,k,N,N

D − ttr,k,i,i +1

 

(16)

(17)

i=1



B,U Ptr,k,i, · TW j



(18)

tr=0 k=1 i=1 j=i+1

EW T U =

N−1

,U PTRR,K,i ·L

(19)

i=1

Eq. (16) calculates the total passenger waiting times on the platforms to board trains in the upstream direction. A Eq. (17) calculates the total passenger on-train times, where ttr,k, is the arrival time of train k at passenger destij, j+1 nation station j. When passengers leave their destination stations at the check-out gates, the metro trips are completed. Eq. (18) calculates the total passenger walking-out times at destination stations, where TW is the walking-out time per traveler. If passengers are left at the platforms after the last train has finished its service, extra waiting times will be added as

S. Yang, F. Liao and J. Wu et al. / Transportation Research Part B 133 (2020) 85–113

93

Fig. 2. The different train links in a space-time diagram.

,U a penalty, which is formulated as Eq. (19), where L is a positive number and PTRR,K,i is the number of passengers who cannot

catch the last train K at i. Similar to the above calculation process, WTD , OTD ,WOTU and EWTD can also be calculated in the downstream direction. 3.3. Utilization of regenerative energy In this subsection, we formulate the utilization rate of regenerative energy when energy transmits between multiple trains and introduce the energy exchanging constraints. (1) Utilization function of regenerative energy To estimate the magnitude and direction of electricity transmission, we use a space-time diagram to visualize the complicated energy flow process, based on which the allocation mechanism is modeled. We introduce several energy exchange situations between accelerating and decelerating trains and develop suitable strategies for calculating the use of regenerative energy. In Fig. 2(a), the ticks show the preferred overlap between the accelerating and decelerating trains; the cross shows that there is no energy transmission; the circle shows that a portion of regenerative energy will be transmitted. In Fig. 2(b), the solid and dashed lines show the running trains in different directions; the red and green curves show the flows of the traction and regenerative energy; the gray shadows show the different decelerating indices; the different curve directions correspond to different train running directions. In Fig. 2(c), the different shadows show the different relative position changes in each decelerating phase. In Fig. 2(d), the shadows show the different allocation modes at each time step in a decelerating phase. A decelerating index refers to a decelerating phase of a train, labeled as z = 1, …, (TR + 1) • (N − 1). In a single running direction, each train decelerates N − 1 times. During the deceleration phases, if there are accelerating trains within the maximum electricity transmission distance, the regenerative energy can be used by the accelerating trains. Thus, the time and positional information of accelerating trains in a decelerating index is important. The larger overlapping time and the closer relative distances between the accelerating and decelerating trains, the more regenerative energy will be utilized. For example, Fig. 2(a) illustrates the energy exchange between a decelerating train (the green solid line) and multiple accelerating trains (in red lines) in the space-time diagram. The different dashed lines represent the different beginning times of the accelerating trains in a decelerating phase. Fig. 2(b) shows that the traction (red curves) and regenerative (green curves) energy exchange during the time instances in the shadow areas. Fig. 2(c) displays four different situations of interactions in

94

S. Yang, F. Liao and J. Wu et al. / Transportation Research Part B 133 (2020) 85–113

the space-time diagram, depending on the relative moving positions and speeds of the accelerating and decelerating trains, which are the determinants of exchange magnitude. Therefore, once there are multiple accelerating trains in a single decelerating phase, an energy allocation mechanism is required. For example, in Fig. 2(d), the solid and dashed lines show the different running directions; the different color lines represent the three different running phases; the different shadow shapes represent the different allocation rates. If regenerative energy is used by two different trains (red and black shapes) simultaneously, the energy allocation mechanism will affect the utilization rates (shadow areas with crosses). The average utilization rate function can be obtained as



u=

K Z

U,U Ek,z

+

k=1 z=1

K Z



U,D Ek,z

k=1 z=1

K Z k=1 z=1

B,U Ek,z

+

K Z



B,D Ek,z

(20)

k=1 z=1

B,U U,U where Ek,z is the total energy regenerated by train k of the up-train fleet in decelerating phase z, and Ek,z is the amount in B,U U,D B,D Ek,z that is actually consumed. Ek,z and Ek,z are the use of regenerative energy and regenerative energy by train k of the down-train fleet in decelerating phase z.

(2) Energy exchanging constraints In a decelerating phase, the regenerative energy that is not used by the accelerating trains is wasted; in case of insufficient regenerative energy, the accelerating trains use traction energy. In addition, the accelerating trains out of the decelerating phases cannot use regenerative energy. The magnitude of regenerative energy of the up-train fleet is determined by the utilization coefficients of the accelerating trains, overlapping times, and the magnitudes of accelerating and decelerating U,U B,U energy of the up-train fleet in a decelerating phase. Ek,z and Ek,z are formulated as follows. t E,U

U,U Ek,z =

k,z



B,U τ =tk,z +1 k ∈ψk,z



min EkA ,z (τ ),

 B,U μk,k ,z (τ ) · Ek,z (τ )

ψ¯ k,z = ψ¯ k,z / k , ∀k ∈ ψk,z ⎧     ⎨E B,U (τ ) = M · θ2 · vk,z (τ + 1 ) 2 − vk,z (τ ) 2 /2, t B,U ≤ t ≤ t E,U k,z k,z k,z     ⎩E A (τ ) = M · vk ,z (τ + 1 ) 2 − vk ,z (τ ) 2 /2θ1 , t B ≤ t ≤ t E k ,z k ,z k ,z μk,k ,z = αk,k ,z · σk,k ,z ⎧   dk,k ,z (τ ) ⎨ αk,k ,z = β1 · exp −β2 · , dk,k ,z (τ ) ≤ dmax dmax ⎩ αk,k ,z = 0 , dk,k ,z (τ ) > dmax dk,k ,z (τ ) = | pk,z (τ ) − pk ,z (τ )|

 σk,k ,z (τ ) = αk,k ,z



(21)

(22)

(23)

(24)

(25)

(26)

 αk,k ,z

(27)

k ∈ψk,z

U,U Eq. (21) calculates the use of regenerative energy Ek,z in z from k of the up-train fleet based on mechanical power equation (e.g., Wang et al., 2015; Yang et al., 2015; Yin et al., 2017), where μk,k , z (τ ) is the utilization coefficient of accelerating B,U E,U train k in z at time τ , k is an accelerating train that can use the regenerative energy from k, tk,z and tk,z respectively

represent the beginning and ending times of k in z. tkB ,z and tkE ,z respectively represent the beginning and ending times of k in z of k. τ takes 0.1 s as the time unit and thus τ + 1 corresponds to the next time step after τ . During this time, if several accelerating trains exist simultaneously, we introduce a set ψ k,z that includes those accelerating train indices. Eq. (22) is a condition to ensure that one accelerating train receives energy from at most one decelerating train at one time, while one decelerating train may still transmit energy to multiple accelerating trains, where ψ¯ k,z is the complementary set of ψ k,z regarding k and z. The manipulation means that for any accelerating train in set ψ k,z , the accelerating train index will be B,U removed for further consideration. Ek,z and EkA ,z are the regenerative energy by the up-train fleet k and the traction energy

by k in z respectively, which can be calculated by Eq. (23), where M is the train mass, vk,z (τ ) is the speed of k at τ in z, θ 1 and θ 2 are the conversion factors from electricity to kinetic energy and from kinetic energy to electricity respectively. Particularly, for calculating the utilization energy, we need first to obtain the utilization coefficients according to Eq. (24). Then, by comparing the traction energy and available utilization energy for each accelerating train at each time step, the smaller value is selected. Thus, we obtain the utilization energy of one accelerating train at each time step. Finally, to obtain

S. Yang, F. Liao and J. Wu et al. / Transportation Research Part B 133 (2020) 85–113

95

the total utilization energy in a decelerating phase, we sum up the utilization energy of all accelerating trains appearing in z. As a result, the utilization coefficient μk,k , z (τ ) of k in z can be obtained by Eq. (25), where α k,k , z is the electricity decay factor, σ k,k , z is the proportion rate of energy allocation. To calculate the electricity decay factor, the resistance is a key factor, which depends on the length (from its source to the end) and the material itself. In a basic electric circuit, the resistance can be represented by = (ρ / ) · d, where ρ is the resistance rate and is the cross-sectional area of the material. In the metro system, the resistance is used to estimate the energy loss of the overhead connecting line materials (Liu et al., 2018). Therefore, based on the relative distance between accelerating and decelerating trains, we propose a nonlinear function to calculate α k,k , z by Eq. (26), where β 1 and β 2 are two resistance coefficients, dmax is the maximum distance that the instant power can reach by an electrical line. Since the energy exchanges occur at the same time, the relative distance dk,k , z (τ ) between k and k at τ can be calculated by Eq. (26), where pk (τ ) and pk (τ ) are positional information of k and k at τ . When several accelerating trains appear at the same time in z, it resembles the formation of a simple parallel circuit. To decide how much regenerative energy is distributed to the accelerating trains, an allocation rate σ k,k , z determined by the decay factors and the number of accelerating trains at τ in z is calculated by Eq. (27). Similarly, the use of regenerative U,D B,D energy Ek,z and the regenerative energy Ek,z for the down-train fleet can be obtained. 3.4. Bi-objective optimization model We conclude a bi-objective optimization model incorporating the energy allocation and the passenger assignment as



min [W T +OT + EW T + W OT ] max [u] s.t. Constraints (3 ) − (15 ) and (21 ) − (27 )

(28)

where constraints (9), (15), (21), (22), (23), (26), and (27) are all non-linear.

4. Solution procedure Either of the objectives Eqs. (1) and (2) can be categorized as a large-scale mixed-integer programming problem and thus is NP-hard (Kang et al., 2015; Corman et al., 2017; Guo et al., 2017; Wang et al., 2017; Yin et al., 2017; Van Thielen et al., 2018). For example, for the instance referred in Section 3.3 with discrete parameters (time, speed, and utilization rates, etc.), there are 20 × 7200 = 144,0 0 0 integer variables for representing speed, regenerative energy or traction energy in the discrete-time domain if we set the system time (t) and electricity transmission time (τ ) as 1 s and 0.1 s respectively for 20 trains running in a 2-h period. When trains run through 12 sections for 2 cycles, there are 12 × 4 = 48 decelerating phases. If the minimum running time is 100 s, there are 48 × 10 0 0 = 48,0 0 0 binary variables that represent whether the energy is exchanged or not. As we discretize time, these numbers refer to the number of decision variables in the temporal dimension. Thus, the optimization problem is difficult and needs an efficient approach to solve. In the literature, there are two common categories of methods to optimize energy-efficient timetables when considering the two objectives: genetic algorithms (e.g., Yang et al., 2013) and Lagrangian-based algorithms (e.g., Yin et al., 2017). The objectives are highly nonlinear and the constraints are not easy to be transformed into linear conditions. In this study, we modify a fast probability-based evolutionary method, i.e. NSGA-II (Deb et al., 2002) to solve the bi-objective optimization problem. It has been widely used in multi-objective optimization problems (e.g., Pacheco et al., 2013; Yu et al., 2015; Liu et al., 2016). The detailed flowchart of the NSGA-II is shown in Fig. 3. In the algorithm, the first step involves generating a number of random feasible irregular timetables (chromosomes in genetic algorithm terms), in which the running, headway, and dwelling times tend to be different. After evaluation, the solutions will be sorted by non-dominated sorting. The solutions are then selected in a tournament for the crossover and mutation operators. Given the special structure of the solutions, we redesign the crossover and mutation operators to generate offspring solutions, after which the old and new candidate solutions are combined. The ranks and crowding distances are calculated by non-dominated sorting. Finally, we choose the same population size of solutions for the next iteration if the termination condition is not satisfied.

4.1. Generation of an initial solution A solution includes detailed time and positional information of the trains. The blocks of a solution are shown as

X = [X1 , ..., XB ] where X is a space-time matrix and B is the block number (B = TR + 1). X1 is a component of X, which includes the trains running from terminus 1 to N (up-train fleet) and thereafter from N to 1 (down-train fleet). For example, below are two

96

S. Yang, F. Liao and J. Wu et al. / Transportation Research Part B 133 (2020) 85–113

Fig. 3. The main loop of the NSGA-II algorithm.

blocks in a solution to show the relationship between the upstream and downstream directions,

⎡x U

1,1 xD 1,1

⎢ ⎢ ⎢.. ⎢. ⎢ ⎢. . . X1 = ⎢ ⎢. . . ⎢ ⎢. ⎢.. ⎢ ⎣x U

K,1

xD K,1

... ... .. . xU k,n xD k,n .. . ... ...

... ... .. . xU k,n+1 xD k,n+1 .. . ... ...

⎤

⎡

xD xU 1,2N 1,2N−1 D x1,2N−1 ⎥ ⎢xU1,2N ⎥ ⎢ .. ⎥ ⎢.. ⎥ ⎢. . ⎥ ⎢ ⎥ ⎢ ... ⎥; X2 = ⎢. . . ⎥ ⎢. . . ... ⎥ ⎢ ⎥ ⎢. .. ⎥ ⎢.. . ⎥ ⎢ ⎦ ⎣xD xU K,2N K,2N−1 xU xD K,2N−1 K,2N

... ... .. . xD k,n xU k,n .. . ... ...

... ... .. . xD k,n+1 xU k,n+1 .. . ... ...

⎤

xD 1,4N−2 xU 1,4N−2 ⎥ ⎥ .. ⎥ ⎥ . ⎥ ⎥ ... ⎥ ⎥ ... ⎥ ⎥ .. ⎥ . ⎥ ⎦ xD K,4N−2 xU K,4N−2

where xU refers to train k at element index n running in the upstream direction; xD refers to train k at n running in the k,n k,n downstream direction. n is an element index that alternately represents a station or section. More specifically, the different element indices of x consist of different time and positional information. In addition, the element sequence guarantees the regular structure of each block. In X1 , the element index n is ranged from 1 to 2N − 1, in which the odd element indices represent the stations and even indices represent the sections; in X2 , the element index is ranged from 2N to 4N − 2, in which the odd element indices represent the sections and even indices represent the stations. For example, if tr = 0, two adjacent element indices correspond to the arrival-departure information at the stations and switching-point information on a section respectively, stated as

xU = k,n



xU = k,n+1

 

U,A U,D tk,n , tk,n ;



,A ,D pUk,n , pUk,n

  U,WD

U,WA tk,n , tk,n+1 ; +1

  U,WD

,WA pUk,n , pk,n+1 +1

, n = 1, 3, ..., 2N − 1 , n = 1, 3, ..., 2N − 3

(29)

U,A U,D where tk,n and tk,n are the arrival and departure times at the stations in the upstream direction. Accordingly, the variables

,A ,D U,WA ,WA ,WD pU and pU are equal to the station position. tk,n , pU , t U,WD and pU are the time and positional information of k,n k,n +1 k,n+1 k,n+1 k,n+1 the switching points on the sections. The rest rows represent information in the same manner. To generate an initial feasible population subject to the constraints (Section 3.1), we extend the principle of a parallelogram for guaranteeing the diversity and feasibility of a solution. Note that when a train turns around at the terminus in the upstream direction, the element indices will change accordingly. For example, in Fig. 4, the green circle and ring are the arrival and departure times of train k at station N − 1; the black circles and rings are the switching points between stations N − 1 and N. To obtain the complete feasible solutions satisfying the space-time constraints, we develop a random generation algorithm, labeled as Algorithm 1 in Appendix B.

S. Yang, F. Liao and J. Wu et al. / Transportation Research Part B 133 (2020) 85–113

97

Fig. 4. The element index changes at terminus N.

Fig. 5. Randomly generate a feasible solution.

Note that if element index n corresponds to a running section, n − 1 and n + 1 correspond to two adjacent stations. In Algorithm 1, the generation of the headway-dwelling and headway-running time pairs should satisfy the following constraints: D A DW Hk,n − TkDW = Hk,n − Tk,n −1 +1,n−1 −1 −1 A D Hk,n+1 − Tk+1,n = Hk,n−1 − Tk,n A D DW H ≤ Hk,n , Hk,n ≤ H¯ ; T ≤ Tk,n ≤ T¯ ; T DW ≤ Tk,n ≤ T¯ DW

(30)

A D where Hk,n and Hk,n are the arrival and departure time headways respectively. TkDW is the dwelling time of k + 1 +1 −1 +1,n−1 at n − 1. Tk + 1, n is the running time of k + 1 in n. A diagram (Fig. 5) is used to show the generation process of parallel A A DW solutions. To be specific, the arrival-arrival headway Hk,n of k at n − 1 is equal to tkA+1,n−1 − tk,n ; the dwelling time Tk,n −1 −1 −1 D A D D of k at n − 1 is equal to tk,n − tk,n . The departure-departure headway Hk,n of k at n − 1 is equal to tkD+1,n−1 − tk,n ; −1 −1 −1 −1

A D DW , running time T the running time Tk,n of k in n is equal to tk,n − tk,n . In Fig. 5, the dwelling time Tk,n k,n and headway +1 −1 −1 DW ¯ DW A time H are randomly generated in [T , T ], [T , T¯ ] and [H , H¯ ] respectively. Thus, under constraint Eq. (3), the first k,n−1

equation of the right-hand side of Eq. (28) is calculated as a constant value. For example, if n = 2 and k = 1, we have t1A,1 = 0. Then, after generating random T1DW ∈ [T DW , T¯ DW ],T1,2 ∈ [T , T¯ ], and H1A,1 ∈ [H , H¯ ], we have t1D,1 = t1A,1 + T1DW , t1A,3 = t1D,1 +T1,2 ,1 ,1 A and t2A,1 = t1A,1 +H1A,1 . Similarly, another pair (Hk,n ,T ) can also be obtained with the method. +1 k+1,n In line with the time series of the trains, we also obtain the positional series to make the timetable complete. When calculating the running distances of the train in different directions, the distances of trains are increasing in the upstream direction, and decreasing in the downstream direction. The method to obtain the positional information of the switching points follows the same rules above.

4.2. Evaluation of objectives After generating an initial timetable solution, we introduce the evaluation of two objectives below. Due to the complexity of each nonlinear objective, a fast evaluation method is needed to improve the performance of the whole algorithm. (1) Passenger travel time

98

S. Yang, F. Liao and J. Wu et al. / Transportation Research Part B 133 (2020) 85–113

To calculate the objective function of passenger travel time (Eq. (1)), we have formulated the time components in Section 3.2. Combining the initial timetable, we assign the passengers in the metro line to the trains at each stage. For example, the matrix of waiting passengers for k on each OD pair can be calculated, according to the timetable-oriented OD matrix

⎡

,U PW k

where

W,U Pk,i, j

Pk,W1,,U2

0 ⎢. ⎢ .. =⎢ ⎢. ⎣ .. 0

···

0 ..

···

. ···

U Pk,W1,,N .. . W ,U Pk,N−1 ,N 0

⎤ ⎥ ⎥ ⎥ ⎥ ⎦

is the number of waiting passengers in the upstream direction, which can be calculated based on Eq. (16).

The matrices of the boarded (PBk ,U ), alighted (PLk,U ), and stranded (PRk ,U ) passengers of each train can be calculated similarly according to Eqs. (10-15). Passengers who travel in the downstream direction with k are represented by a lower triangular matrix. After knowing the numbers of passengers and the entry times in a metro line, we evaluate the total time expenses at each stage according to the timetable. (2) Use of regenerative energy Similar to the modeling of energy allocation in Section 3.3, the energy transmission process should be formulated at the microscopic level, in which the time resolution should be high (usually less than 1 s per time step) Yang et al., 2015). According to the Eqs. (21-(27), all the relative states (speed, distance, traction, and regenerative energy) of each train at every time step should be calculated. To decompose this problem, we first develop a fast heuristic method to find the valid decelerating phases and separate the initial solution into two independent time series based on the beginning and ending points of the accelerating and decelerating trains. This operation allows distinguishing how many accelerating trains reside in a decelerating phase, which is shown as follows. Step 1: to simplify the matrixes constructed in Section 4.1, we establish the train accelerating and decelerating point matrixes XA and XD respectively as

⎡xBA

1,1

⎢.. ⎢. ⎢ XA = ⎢. . . ⎢. ⎣..

xBA K,1

⎡xBD

1,1

⎢.. ⎢. ⎢ XD = ⎢. . . ⎢. ⎣..

xBD K,1

xEA 1,1 .. . ... .. . xEA K,1

... .. . xBA k ,z .. . ...

xED 1,1 .. . ... .. . xED K,1

... .. . xBD k,z .. . ...

... .. . xEA k ,z .. . ... ... .. . xED k,z .. . ...

xBA 1,(T R+1)·(N−1 ) .. . ... .. . xBA K,(T R+1)·(N−1 ) xBD 1,(T R+1)·(N−1 ) .. . ... .. . xBD K,(T R+1)·(N−1 )

⎤

xEA 1,(T R+1)·(N−1 ) .. ⎥ ⎥ . ⎥ ... ⎥; ⎥ .. ⎦ . EA xK,(T R+1)·(N−1)

⎤

xED 1,(T R+1)·(N−1 ) .. ⎥ ⎥ . ⎥ ... ⎥ ⎥ .. ⎦ . ED xK,(T R+1)·(N−1)

where xBA includes the beginning time and positional information of accelerating train k in accelerating phase z (1 ≤ k ,z z ≤ (TR + 1)

(tkB ,z ,

pBk ,z );

•

(N − 1)) and xEA includes the ending time and positional information, which are shown as xBA = k ,z k ,z

xEA k ,z

= (tkE ,z , pEk ,z ). tkB ,z and tkE ,z respectively represent the beginning and ending times of accelerating train

B , pB ); xED = (t E , pE ). xBD represents k in z . Similarly, the decelerating matrix XD can also be obtained, where xBD = (tk,z k,z k,z k,z k,z k,z k,z

the beginning time and positional information of train k in decelerating phase z(1 ≤ z ≤ (TR + 1) B information.tk,z

E (ignoring tk,z

•

(N − 1)) and xEA repk,z

resents the ending time and positional and the running direction for the sake of convenience) respectively represent the beginning and ending times of decelerating train k in decelerating phase z. Step 2: to generate an overlapping matrix, we find the beginning or ending time points of accelerating trains under the B , t E ] within maximum electricity transmission distance. Then, the index information of decelerating phase constraints [tk,z k,z the accelerating trains is recorded. If two or more decelerating trains are simultaneously qualified to transmit the regenerative energy to the same accelerating train, only one decelerating train is chosen. Particularly, the train with the earlier beginning time of the decelerating phase is chosen; in rare cases, the one with the smaller train index is chosen if the beginning times are the same. Next, we discretize time and obtain the discrete positions, speed, and energy consumption of the trains appearing in the decelerating phases. Finally, the relevant parameters of the utilization rates at each time step are calculated according to

S. Yang, F. Liao and J. Wu et al. / Transportation Research Part B 133 (2020) 85–113

99

Fig. 6. Block operations.

Eq. (24). To evaluate the average utilization rate, the pseudo-code of the calculation procedure is listed as Algorithm 2 in Appendix B. 4.3. Non-dominated sorting After evaluating the two objectives of the solutions, we use the classic non-dominated ranking method of NSGA-II to rank the solutions. After that, the crowding distances are calculated in each rank. The pseudo-code of the calculation procedure is described as Algorithm 3 in Appendix B (Deb et al., 2002). 4.4. Selection After the solutions are sorted by the non-dominated sorting procedure, a tournament-based selection operator is applied to create the mating pool. In the selection, the best solutions in low (better) ranks have large probabilities to be selected (Deb et al., 2002). The pseudo-code of the calculation procedure is listed as Algorithm 4 in Appendix B. 4.5. Crossover The crossover operator is used in genetic algorithms to produce new timetable solutions from the parent population. The aim is to maintain some chromosome blocks that are potential parts of quality solutions while seeking a certain degree of diversity. After the selection operator, each pair of parent solutions will be operated to generate two offspring solutions based on a certain probability. Based on the structure of a timetable solution (Section 4.1), a timetable solution is partitioned into a number of blocks, which is determined by the turnaround numbers. The information on multiple trains running in the two directions respectively before turning around is considered as a block. The detailed rules can be found in Section 4.1. For example, if two trains run in both directions and turn around at both termini, the timetable solution will be divided into two blocks, which is shown in Fig. 6(a). Thus, if a train turns around TR times, the block number B will be TR+1. The crossover is operated on these B blocks. In the mating pool, two parent solutions are selected for each crossover operator. A probability is set to determine the crossover blocks, which may include BC , 1 ≤ BC < B, blocks. We introduce a probability parameter pC to determine how many blocks are exchanged (Fig. 6(b)). In addition, pC1 , pC2 , and pC3 ( pC1 < pC2 < pC3 ) are three probability values. If pC ≤ pC1 , two blocks at the same position from two parent solutions will be selected to be exchanged; if pC1 < pC ≤ pC2 , four blocks will be exchanged; and if pC2 < pC ≤ pC3 , six blocks will be exchanged. If the difference between the arrival time of the last train and the departure time of the first train is less than the minimum headway time H , the algorithm regenerates the exchanged blocks until the time difference is greater than H . After exchanging the blocks, if the turnaround time of a train is greater than the upper bound time TTR,1 , we reduce the turnaround time; if the turnaround time is less than the lower bound time TTR,2 , we increase the turnaround time. Following this principle, we first calculate the turnaround times of each train between every two consecutive blocks and the headway times between every two consecutive trains. If the turnaround time is greater than TTR,1 , the modified time (Tim,1 ) is equal to the turnaround time (TiTR ) minus the difference between the current headway time (the smallest headway time between two trains) and the minimum headway time (Eq. (31.1)). Then, all the time points of the second train in the latter block are minus the modified time. If the turnaround time is less than TTR,2 , the modified time (Tim,2 ) is equal to the turnaround

100

S. Yang, F. Liao and J. Wu et al. / Transportation Research Part B 133 (2020) 85–113

Fig. 7. Mutation operations on speed profiles.

time (TiTR ) minus the difference between the current headway time (the largest headway time between two trains) and the maximum headway time (Eq. (31.2)). Also, all the time points of the train in the latter block are plus the modified time. The equations are given as follows.

Tim,1 = TiTR − (min (hi,i−1 ) − H ), i = 1, TiTR > T TR,1

(31.1)

Tim,2 = TiTR − max (hi,i−1 ) − H¯ , i = 1, TiTR < T TR,2

(31.2)





where hi,i − 1 is the headway time set at each station between train i and i − 1. These modification mechanisms are applied to all the trains. 4.6. Mutation Recent studies have shown that mutation operators based on local search contribute to speeding up the convergence of the iterative process (Lara et al., 2010). In this subsection, we introduce a local search method based on graphic representations for the mutation operator in the NSGA-II algorithm. This method consists of three heuristic rules: (1) if the ending time and position of the accelerating phase of an accelerating train is located in a decelerating phase, we move backward the speed profiles of the accelerating train to reduce the running time (from the solid lines to the dashed lines in Fig. 7(a)); (2) if the beginning time and position of the accelerating train is located in a decelerating phase, we move forward the speed profile of the accelerating train in parallel (from the solid lines to the dashed lines in Fig. 7(b)); and (3) if multiple accelerating trains use the regenerative energy within a decelerating phase, we change the decelerating rates of accelerating trains (from the solid lines to the dashed lines in Fig. 7(c)). If an accelerating train satisfies scenarios 1 and 2 simultaneously, scenario 2 is applied. If scenario 3 is satisfied, rule 3 is applied even though there are multiple accelerating trains satisfying either scenario 1 or scenario 2. Rule 3 is implemented on multiple accelerating trains simultaneously regardless of the adjustments that have already been made on other accelerating trains. These heuristic rules deal with the relevant variables to create an adjustment set. This set guarantees that the offspring solutions comply with the basic constraints (Eqs. (3–7)). The detailed adjustments of mutation operations are formulated as follows. (1) Moving backward The moving backward operation will reduce the running times of accelerating trains between two adjacent stations. This operation prolongs the overlapping time to reduce regenerative energy loss and meanwhile guarantees no extra time expenses for passengers. Specifically, we increase the dwelling times at the departure stations and keep the arrival times the same (Fig. 8(a)). Note that when changing the running times, the operation should also comply with the constraints of bounds related to headway, dwelling time, and running time, so the departure time can be added with the minimal value in the set tempB = {t empB1 , t empB2 , t empB3 , t empB4}. Each item in the set is calculated as follows. A WA tempB1 = ttr,k,i,i +1 − ttr,k ,i ,i +1

(32)

D D tempB2 = ttr,k  +1,i ,i +1 − H − ttr,k ,i ,i +1

(33)

tempB3

= H¯

D − ttr,k  ,i ,i +1

D + ttr,k  −1,i ,i +1

(34)

A ¯ DW − t D    tempB4 = ttr,k  ,i ,i +1 + T tr,k ,i ,i +1

(35)

tempB5

(36)

=

A ttr,k  ,i +1,i +2

−T

D − ttr,k  ,i ,i +1

S. Yang, F. Liao and J. Wu et al. / Transportation Research Part B 133 (2020) 85–113

Distance

Station i+1

Maximum dwelling time Decelerating ending point Minimum running time Switching time point Minimum headway Maximum headway time

Train k

Station i+1

Minimum dwelling time Switching time point Maximum dwelling time Accelerating beginning point Train k Minimum headway Maximum headway time

Train k+1

The departure time

Distance

Station i

The departure time

Train k+1

Station i

Time

(a) Moving backward

101

Time

(b) Moving forward

Fig. 8. Mutation operations within the flexible domain.

Fig. 9. The stability of the NSGA-II algorithm.

A D where ttr,k,i,i is the arriving time of train k at station i (ending time of the decelerating phase). ttr,k  ,i ,i +1 is the departure +1

WA   time of accelerating train k at i . ttr,k  ,i ,i +1 is the switching time points from accelerating to coasting of k between i and

i + 1, which means that the accelerating trains can use the regenerative energy within the maximum electricity transmission distance. Eq. (32) aims to make the above two time points as close as possible, which are shown by the red and blue points. Eq. (33) guarantees that the headway between k and k +1 is greater than the minimum headway time when putting off the departure time of k . Eq. (34) guarantees that the headway between k and k − 1 is less than the maximum headway time, which is shown by the brown point. Eq. (35) guarantees that the dwelling time at an adjusted station i is less than the maximum dwelling time. The maximum dwelling time point is shown by the yellow point. Eq. (36) guarantees that the running time of train k is greater than the minimum running time, which is shown by the black point in Fig. 8(a). D As a result, the modified departure time ttr,k  ,i ,i +1 is equal to the current departure time plus the minimal value in the above adjustment set 



D D B ttr,k  ,i ,i +1 = ttr,k ,i ,i +1 + min temp



(37)

The modified and real departure times are shown as green and purple points (Fig. 9(a)). (2) Moving forward The moving forward operation means that a train departs and arrives earlier. By reducing the dwelling time at a station, the passenger travel time is reduced when the number of boarded passengers is fixed. Specifically, decreasing and increasing the dwelling times at stations i and i + 1 respectively, this operation should comply with the constraints of bounds related to headway and dwelling times, so the departure and arrival times can be added with the minimal value in the set

102

S. Yang, F. Liao and J. Wu et al. / Transportation Research Part B 133 (2020) 85–113

tempF = {t empF1 , t empF2 , t empF3 , t empF4}. The items are calculated by the following equations WD D tempF1 = ttr,k,i,i +1 − ttr,k ,i ,i +1

(38)

DW D A tempF2 = ttr,k − ttr,k  ,i ,i +1 − T  ,i ,i +1

(39)

D D tempF3 = ttr,k  ,i ,i +1 − H − ttr,k −1,i ,i +1

(40)

D D tempF4 = H¯ − ttr,k  +1,i ,i +1 + ttr,k ,i ,i +1

(41)

A ¯ DW − t D    tempF5 = ttr,k  ,i +1,i +2 + T tr,k ,i +1,i +2

(42)

A A tempF6 = ttr,k  ,i +1,i +2 − H − ttr,k −1,i +1,i +2

(43)

A A tempF7 = H¯ − ttr,k  +1,i +1,i +2 + ttr,k ,i +1,i +2

(44)

WD where ttr,k,i,i is the switching time point from coasting to decelerating of decelerating train k between i and i + 1 (begin+1 ning time of decelerating phase). Eq. (38) aims to make the above two time points as close as possible, which are shown by the red and blue points in Fig. 8(b). Eq. (39) guarantees that the dwelling time at i is greater than the minimum dwelling time. Eq. (40) ensures that the departure-departure headway between two adjacent trains is greater than the minimum headway time when decreasing the departure time of k , which is shown by the black point. Eq. (41) ensures that the headway between two adjacent trains is less than the maximum headway time, which is shown by the brown point. The D  departure time ttr,k  −1,i ,i +1 is shown by the purple point. Eq. (42) specifies that the dwelling time at i +1 is less than the maximum dwelling time, shown by the yellow point (Fig. 9(b)). Eq. (43) ensures that the arrival-arrival headway between two adjacent trains is greater than the minimum headway time when decreasing the arrival time of k , which is shown by the green point. Eq. (44) ensures that the departure-departure headway between two adjacent trains is less than the maximum headway time. D A  Finally, the modified departure time ttr,k  ,i ,i +1 and arrival time ttr,k ,i +1,i +2 of k are reduced by the minimum values in the adjustment set: 





D D F ttr,k  ,i ,i +1 = ttr,k ,i ,i +1 − min temp   A A ttr,k ,i +1,i +2 = ttr,k ,i +1,i +2 − min tempF

(45)

(3) Change the decelerating rates The operation of changing the decelerating rates of accelerating trains does not change the accelerating rate, compared with the previous operations. Herein, the decelerating rate will decrease to b . Since the running times stay the same, there are no new constraints introduced so that no loss of travel time for the passengers. This operation has two advantages: (1) it enlarges the overlapping times between accelerating and decelerating trains in one decelerating phase; and (2) it prolongs the train decelerating times for generating non-zero decelerating phases. As illustrated by the graphic representations, these three operations guarantee no increase in passenger travel time and improve the use of regenerative energy in each decelerating phase. Particularly, the operations of moving backward and forward may reduce passenger travel time. In each iteration of the NSGA-II algorithm, the local search method based on these three heuristic rules applies to all the decelerating phases. 5. Numerical example In this section, a numerical example is presented to test the performance of the developed model and the NSGA-II algorithm. The example concerns the timetabling optimization of a metro line (Yizhuang Line, Beijing) during the morning peak hours (7:0 0–9:0 0 am). The proposed algorithm is coded in MATLAB R2018b and run on a personal computer with Intel(R) Core(TM) i5-8700 @2.3 GHz CPU and 8 GB RAM using Microsoft Windows 10 Enterprise. 5.1. Experiment description The Yizhuang line is a bi-directional metro line of 23.23 km in length and includes 14 stations (13 running sections). Note that the last station (on the right-hand side), Yizhuang railway station (YZ), has been out-of-service since 2010. Table 2 in Appendix C shows the initial inputs for obtaining the random feasible timetables, including the distances of the sections, the maximum and minimum running times on each section, and the maximum and minimum dwelling times at each station. In the metro line, the train fleet size in each direction is K = 10, and the initial turnaround time at the termini is T TR = T¯ TR = 300 s. Table 3 shows the values and units of the relevant train operating parameters. In this numerical example, we

S. Yang, F. Liao and J. Wu et al. / Transportation Research Part B 133 (2020) 85–113

103

Table 3 Values and units of the relevant train operating parameters. Parameter

Value

Unit

Parameter

Value

Unit

M

287,080 0.7 0.8 3 80 1 90 300

kg – – – s – s s

a b b r dmax

0.8 0.5 0.8 0.05 3500 2 180 60

m/s2 m/s2 m/s2 m/s2 m – s s

θ1 θ2 TR TW

β1

H TTR,1

β2 H¯

use the smart-card data of time-variant arrival rates and OD demands, which are listed in Tables 4 and 5 in Appendix D. The time horizon is extended by two hours after the peak hours for ensuring the completeness of the optimization model. We consider two train loading capacities as 1200 and 600 respectively. According to the input of passenger demands, they correspond to two situations respectively, i.e., “no congestion” and “congestion”. The large valueLis set as 10 0 0. The settings for the train operating parameters are given in Table 3. In the following, we consider the first trips in both directions starting at 7:00 am, which is shown in the timetable as 0 for coding convenience. Given the time frame, it is found there are four blocks (B = 4) in a timetable solution. The ending time of the last train is not a fixed value due to the headway fluctuations. For the NSGA-II algorithm, the basic parameters are set as: population size (pop) is 30; the maximum number of generations (gen) is 20 or the algorithm terminates when the Pareto frontier is unchanged after 5 iterations; the crossover and mutation probabilities (pr and pm ) are 50% and 100% respectively; the block numbers 1, 2 and 3 are associated with probabilities pC1 = 25%, pC2 = 75%, and pC3 = 1 respectively; the tournament size is 30. 5.2. Computational results Based on the above setup, we run the NSGA-II procedure on the real-world case study. First, for the fixed passenger demands, we test the efficiency of our model. Then, the convergence and stability of the algorithm are analyzed. Thirdly, we randomly generate a population of feasible timetables and compare the irregular and regular timetables intuitively. Next, the bi-objective optimization model is solved and a Pareto frontier is derived. With the Pareto frontier, an optimal solution is selected for illustration, for which the effects on passenger travel time and energy consumption are further discussed. For considering a fixed passenger demand, we select a fixed time period from 7:00 am to 9:00 am by setting the arrival rates of passengers as 0 after 9:00 am. Meanwhile, we discard the sub-solutions after 9:00 am (that is, only the train timetable before 9:00 am is considered). Considering the congestion situation (train capacity is 600 p/t), the average passenger travel time is reduced by 8.5% (from 5.7728 × 107 s to 5.2789 × 107 s. Meanwhile, the utilization rate of regenerative energy is increased by 20.06%. We run the NSGA-II procedure ten times and each takes on average 0.64 h, of which around 75% of the time is used to evaluate the bi-objective functions. The coefficients of variance of utilization rates and travel times in the first ranks (including the best solutions that are non-dominated) are 0.0437 and 0.0127 respectively, which indicate the stability of the developed algorithm (Fig. 9). (1) Irregular and regular timetables For illustrating the timetabling outcomes, we display two timetables that are generated randomly by the method discussed in Section 4.1. In Fig. 10, the different colors show the different trains, which are labeled in the legend on the righthand side of the figure. The regular timetable has constant time components (e.g., running time, headway, and dwelling time) between two fixed trains (Fig. 10(a)); differently, the time components vary in the irregular space-time diagram (Fig. 10(b)). (2) Passenger travel time To show the improvement in passenger travel time without considering congestion (train loading capacity is 1200), we arbitrarily choose an optimized timetable from the Pareto solutions. According to the formulations in Section 3.2, the waiting passengers are assigned to adjust the arrival and departure times. For example, there are two trains running in the upstream direction after turning around three times, and these two trains are from the initial and optimized solutions respectively (Fig. 11). In Fig. 11(a and b), the horizontal axis denotes the station indexed from 1 to 12. Station 13 is the origin when running in the downstream direction and all the arriving passengers board the trains in the opposite direction. Therefore, the waiting time is 0. The vertical axis represents the accumulative numbers of boarded passengers. In Fig. 11(c and d), we show the changes in waiting times at each station, recorded from the time passenger arriving at the platform to the time boarding a train, according to Eq. (10). Note that different colors in Fig. 12 are used to represent different destination stations. When adjusting the train schedule, the numbers of boarded passengers who take stations 2 and 3 as destinations

104

S. Yang, F. Liao and J. Wu et al. / Transportation Research Part B 133 (2020) 85–113

Fig. 10. Regular and irregular timetables.

Fig. 11. The boarded numbers and waiting times of passengers on each O-D pair.

are reduced, so are the waiting times. Taken together, the passenger waiting times for this train is reduced by 2.70% (from 9.476 × 104 s to 9.218 × 104 s. To further test the effects on passenger travel time, we consider the “congestion” situation. When decreasing the train loading capacity by half to 600, congestions exist at most of the stations. Fig. 12 shows the waiting, on-train, extra waiting, and walking-out times of a train at each station respectively. Note that Fig. 12(a–c) only shows 12 stations regarding the calculation subjects. For example, the passenger waiting time at station 4 has decreased significantly to the blue circles (Fig. 12(a)). The proposed strategies (Section 4.6) offer more available spaces on the trains for most congested stations via adjusting the train departure times (Fig. 12(b)). As shown in Fig. 12(c), congestion exists in the initial solution from station 4 to 10 shown by the orange stars. In the optimized solution, the extra waiting times at almost all stations decrease. In Fig. 12(d), as expected, the total walking-out times are little changed because the model only adjusts the proportions at each station. In total, passenger travel time can be saved by 226.25 h considering all the trains in both directions. (3) Energy utilization According to the energy allocation mechanism discussed in Section 3.3, the relevant variables (speed, time, distance, allocation rates, regenerative and traction energies) are discrete in each decelerating phase. When the decelerating and accelerating trains are found in a decelerating phase, we discretize the above variables every 0.1 s. We obtain two categories of decelerating phases to show the changes in the variables. One is to include only one decelerating train and one

S. Yang, F. Liao and J. Wu et al. / Transportation Research Part B 133 (2020) 85–113

105

Fig. 12. Travel time consists of four items.

Fig. 13. Train relevant variables in a decelerating phase.

accelerating train in a decelerating phase; another is to include one decelerating train and multiple accelerating trains. In the first category, in Fig. 13(a, b), the blue curves show the distance and speed changes of the decelerating train, and the red curves show the changes of the accelerating train. The near-continuous curves show that the time resolution is sufficient. Regarding the second category, since the differences are caused by the relative speeds and distances between the decelerating train and multiple accelerating trains, the energy allocation mechanism takes effect to determine the proportion of energy-assigning rates. Due to the similarity of the speed and distance changes, we only illustrate the curves of utilization energy volumes and rates. For example, seven accelerating trains use the regenerative energy generated by one decelerating train. In Fig. 13(c, d), the different colors show the allocating utilization results of the different accelerating trains. The utilization rate of regenerative energy (u) for each accelerating train is calculated by Eq. (20) and the utilized regenerative U ) is calculated by Eq. (21). When only one accelerating train shares the regenerative energy, the utilization rates energy (Ek,z will increase significantly as shown in the top right corner in Fig. 13(d). Fig. 14 shows the comparison of decelerating phases of an initial and the above selected optimized timetable respectively. Each includes 960 possible decelerating phases considering that 20 trains run on 12 sections four times bi-directionally. The odd numbers in the vertical axis represent the trains in the upstream direction, and the even numbers represent the trains in the downstream direction. The horizontal axis represents the whole decelerating phases generated by all trains

106

S. Yang, F. Liao and J. Wu et al. / Transportation Research Part B 133 (2020) 85–113

Fig. 14. Utilization rates in each decelerating phase.

Fig. 15. Comparison of the energy consumption between initial and optimized solutions.

after finishing two cycles, and the vertical axis represents train numbers in both directions. The different colors show the different utilization rates of regenerative energy, which are remarked by the legend on the right-hand side of the figure. If the regenerative energy dissipates completely, we call it a zero-decelerating phase, in which the utilization rate is 0 shown by the white background. According to the different colors of the grids in Fig. 14, there are 300 and 319 non-zero decelerating phases in Fig. 14(a) and (b) respectively. In sum, the average utilization rate is increased from 13.68% to 21.65% after applying the optimization model. To illustrate that our model reduces the total energy consumption in a consistent way of enhancing utilization rate, we choose a set of solutions from the initial population and compare the results with the optimized solutions in one algorithm run. The traction, regenerative, and total energy are shown in Fig. 15. In Fig. 15, the red, blue, and green colors correspond to the traction, total, and regenerative energy, respectively. The stars represent the results calculated with the initial solutions; the circles represent the results calculated with the optimized solutions. As shown by the optimized solutions, the average traction energy is increased from 3.6455 × 1010 J to 3.7262 × 1010 J and the total energy is decreased from 3.5252 × 1010 J to 3.3755 × 1010 J. Thus, the utilization rate is increased from 13.68% to 21.65% (the utilized regenerative energy is increased from 1.2036 × 109 J to 3.5073 × 109 ). It shows that the maximum utilization rate and minimum total energy consumption are consistent objectives in the example.

S. Yang, F. Liao and J. Wu et al. / Transportation Research Part B 133 (2020) 85–113

107

Fig. 16. Pareto frontier.

Fig. 17. An irregular timetable.

(4) Pareto frontier To explore the effectiveness of the optimization model, the Pareto frontier of one arbitrary run is shown in Fig. 16. The horizontal axis represents the total passenger travel time, and the vertical represents the utilization rate of regenerative energy. For a clear illustration, we convert the objective of maximizing the utilization rate to minimizing the rate gap to 1 (100%). These two subfigures illustrate the tradeoff between the two objectives (total travel time and transformed utilization rate) and the convergence process of the NSGA-II algorithm. Fig. 16(a) shows the values of every five generations. In Fig. 16(b), the black star in the top right corner denotes the result of the currently used regular timetable; the blue stars represent the bi-objective values calculated using the optimized irregular timetables in the first rank. Comparing the results between regular and irregular timetables, the average utilization rate of regenerative energy is increased from 0.06% to 21.65%, and the average passenger travel time is reduced by 15.21%. Fig. 17 shows a timetable solution in the first rank of an arbitrary run. It takes around 5 s to generate a solution by Algorithm 1. This formulation of the timetable cannot only greatly enhance the flexibility of metro systems, but also increases the convenience for evaluating the objectives. In addition, the irregular timetables include many intersecting links in the space-time network at station 1 and 13 in both directions, during which there may be some use of regenerative

108

S. Yang, F. Liao and J. Wu et al. / Transportation Research Part B 133 (2020) 85–113

energy. However, when calculating based on a timetable in a single direction, there is no use of regenerative energy since the regenerative energy is neglected by the turnaround trains. All in all, according to the above numerical results, the suggested NSGA-II algorithm can produce a set of quality Pareto frontier solutions for the bi-objective optimization problem. These results confirm that our proposed model and algorithm are capable of capturing energy exchanges among the trains, increasing energy re-use, and reducing the total travel time sufficiently. 6. Conclusions In this study, we consider the joint optimization of energy consumption and passenger travel time considering passenger assignment and energy allocation in an energy regenerative metro system. To achieve an optimal tradeoff between the passenger travel time and the use of regenerative energy, we suggest an integrated bi-objective optimization model. With this model, we can optimize the arrival, departure, and dwelling times simultaneously to generate an irregular timetable of multiple trains in a long time period. An NSGA-II based algorithm is developed in combination with the space-time characteristics of train timetables to derive the Pareto frontier efficiently. The efficiency and effectiveness of the proposed model are demonstrated by a numerical example and the sensitivity analysis of train loading capacity. With an arbitrary optimized solution from the Pareto frontier, the use of regenerative energy is improved up to 20.09% and the total travel time is saved up to 3717 h compared with a fixed regular timetable. Our future research will focus on the following aspects. First, although we have considered time-dependent passenger demand, passenger demand in real conditions is more complicated. For example, the arrival rates may be counted at the check-in gates and the walking times should be considered as variable values. Extending our research to consider passenger choice behavior is also meaningful when calculating travel times. Addressing these issues requires dedicated data. Second, we explored an energy allocation problem only for a single line. However, electricity transmission is a complicated issue when considering the instantaneous distribution of electric energy among multiple lines. To make the electricity distribution more efficient and reliable, we should consider more general functions of electricity attenuation and allocation mechanisms in a full metro system. Third, with the consideration of the above aspects, efficient crossover and mutation operators should be developed accordingly to obtain quality Pareto frontier solutions. Fourth, we will consider an efficient way of constructing the solution structures in case some trains finish operations at a terminus. Fifth, we will test the differences between adjusting the deceleration and acceleration phases and the many-to-many energy transmission mechanism. CRediT authorship contribution statement Songpo Yang: Conceptualization, Methodology, Writing - original draft. Feixiong Liao: Conceptualization, Methodology, Writing - review & editing, Supervision. Jianjun Wu: Conceptualization, Methodology, Writing - review & editing, Supervision. Harry J.P. Timmermans: Conceptualization, Writing - review & editing. Huijun Sun: Writing - review & editing. Ziyou Gao: Writing - review & editing. Acknowledgments We are grateful to the Editor-in-Chief and four anonymous referees whose comments have led to a significant improvement of the paper. This work was supported jointly by the China National Funds for Distinguished Young Scientists (Grant no. 71525002), the National Natural Science Foundation of China (Grant no. 71890972, 71890970, 91846202, 71621001), the State Key Laboratory of Rail Traffic Control and Safety (Grant no. RCS2020ZZ001), the Fundamental Research Funds for the Central Universities (No. 2018RC010), the Project Funded by China Postdoctoral Science Foundation (Grant no. 2019M660373) and the Netherlands Organization for Scientific Research (NWO). Appendix A. Calculation of switching points Yang et al. (2018b) obtained the coasting time function w.r.t. the distance and running time on section i by using a strict mathematical derivation, shown as

!

WD ttr,k,i,i +1

WA − ttr,k,i,i +1

=

2 2r · ltr,k,i,i+1 A  A D − · ttr,k,i − ttr,k,i,i +1,i+2 +1 C C

A D where ttr,k,i − ttr,k,i,i is the running time of train k between station i and i + 1 after turning around tr times, li,i + 1 is +1,i+2 +1 the distance between i and i + 1. r is the resistance rate. A and C are two intermediate parameter functions. The relationships w.r.t. accelerating rate a, resistance rate r and decelerating rate b are stated as





 2   · b · r − b2     2 b 2

A = (Aa )2 · a · r + a2 + Ab



C = a·r+a

 2

· ( B a )2 + b · r − b

· B

S. Yang, F. Liao and J. Wu et al. / Transportation Research Part B 133 (2020) 85–113

where the intermediate variables Aa , Ba , Ab , Bb are stated as

Aa = 1 −

a a+r a a+r , Ba = − 1, Ab = , Bb = a+b a+b a+b a+b

The detail calculation process can be found in (Li et al., 2014b; Yang et al., 2018a,b). Appendix B. Pseudo-codes of different modules

Algorithm 1 Randomly generate a feasible timetable. (1) (2) (3) (4) (5) (6) (7) (8) (9) (10) (11) (12) (13)

For each k ≤ K do For each n ≤ (TR+1) • (2N − 1) do If k = 1 and n = 1 A DW A = 0, generate random Tk,n ∈ [T DW , T¯ DW ]; Tk,n+1 ∈ [T , T¯ ];Hk,n ∈ [H , H¯ ] tk, 1 D A DW A D A A = tk,n + Tk,n , tk,n = tk,n +Tk,n+1 , tkA+1,n = tk,n +Hk,n Set tk,n +2 Else if n=(tr+1) • (2N − 1) and k = 1 A D D A DW = tk,n +T TR , tk,n = tk,n +Tk,n tk,n +1 +1 +1 +1 Else A A A ¯ ], t A ∈ [ H = tk,n +Hk,n ; , H Generate random Hk,n −1 k+1,n−1 −1 −1 D A ) and (Hk,n , TkDW , Tk+1,n ) with constraints (29) Generate random pairs (Hk,n −1 +1,n−1 +1 End End End

Algorithm 2 Valuation of utilization energy. (1) (2) (3) (4) (5) (6) (7) (8)

For each k ≤ K do For each z ≤ Z do For each k ∈ ψ k,z Calculating τ k,z , τ k ,z , pk,z , pk ,z , vk,z , vk ,z , Ek,z , Ek ,z in 0.1 s U , μk,k ,z Calculating dk,k ,z , αk,k ,z , σk,k ,z , Ek,z End End End

109

110

S. Yang, F. Liao and J. Wu et al. / Transportation Research Part B 133 (2020) 85–113

Algorithm 3 Non-dominated sorting. i. Non-dominated ranks (1) For each s= 1:pop do (2) Ss = ∅ and ns = 0

s and g are solutions; Ss contains solutions dominated by s; ns is the number of solutions dominating s.

(3) For each g= 1:pop do (4) If ps ≺pg then (5) S s = S s ∪{ s } (6) Else if pg ≺ ps then (7) ns = ns + 1 (8) End (9) If ns = =0 then (10) ranks = 1; F1 =F1 ∪{s} (11) End (12) End (13) Initialize i = 1 (14) While Fi is not empty do (15) .. (16) For each si ∈ Fi (17) For each g ∈ Ss (18) ng = ng − 1 (19) If ng = =0 then (20) ranksi+1 = i+1; Q = Q∪{g} (21) End (22) End (23) End (24) i = i + 1; Fi = Q (25) End ii. Crowding distance calculation (1) While Fi is not empty do (2) Set the Fi (dj )= 0; omax = 2 (3) For each o = 1: omax do (4) I = sort(Fi ,o)

ps dominates pg .

s belongs to the first rank.

i is a rank counter. Q is used to store the solutions of the next rank. si is a solution that belongs to rank i. g belongs to the first rank. si + 1 belongs to the next rank.

Initialize the crowding distance dj . j corresponds to the j-th solutions in rank Fi . omax is the maximum number of objectives. Sort using objective values from small to large, where o is the objective index. To ensure the boundary points are always selected, wheren is the number of solutions in the rank.

(5)

I(di ) = ∞ and I(dn ) = ∞

(6) (7)

For each k= 2:n − 1 do I (dk ) = I (dk ) + (I (k + 1 ).o − I (k − 1 ).o)/( f omax − f omin )

(8) (9) (10) (11)

I(k).o is the value of the o-th objective function of the k-th solution in I, f omax and f omin are the maximum and minimum values of theo-th objective .

End End i=i+1 End

Algorithm 4 Selection. (1) (2) (3) (4) (5) (6) (7) (8) (9) (10) (11) (12) (13) (14) (15)

Initialize i = 1 While Fi is not empty do ranks1 = 1 when .. belongs in Fi i=i+1 End For each s= 1:pop do For each g= 1:pop do If ranks Fi (dg ) then s ≺n g End End End

≺n is a dominance comparison operator (Deb et al., 2002)

S. Yang, F. Liao and J. Wu et al. / Transportation Research Part B 133 (2020) 85–113

111

Appendix C. The bounds of the train timetable

Table 2 Section index, length, upper and lower bounds of dwelling and running times in the Yizhuang Line.

Index and sections

Length (m)

Minimum dwelling time (s)

Maximum dwelling time (s)

Minimum running time (s)

Maximum running time (s)

1 SJ-XC 2 XC-XH 3 XH-JG 4 JG-YZQ 5 YZQ-WH 6 WH-WY 7 WY-RJ 8 RJ-RC 9 RC-TJ 10 TJ-JH 11 JH-CQN 12 CQN-CQ 13 Turnaround

2631 1247 2366 1983 992 1539 1280 1354 2337 2265 2086 1287 –

20 20 20 20 20 20 20 20 20 20 20 20 –

40 40 40 40 40 40 40 40 40 40 40 40 –

120 120 120 120 120 120 120 120 120 120 120 120 60

150 150 150 150 150 150 150 150 150 150 150 150 540

13 CQ- CQN 12 CQN-JH 11 JH-TJ 10 TJ-RC 9 RC-RJ 8 RJ-WY 7 WY-WH 6 WH-YZQ 5 YZQ-JG 4 JG-XH 3 XH-XC 2 XC-SJ 1 Turnaround

1287 2086 2265 2337 1354 1280 1539 992 1983 2366 1247 2631 –

20 20 20 20 20 20 20 20 20 20 20 20 –

40 40 40 40 40 40 40 40 40 40 40 40 –

120 120 120 120 120 120 120 120 120 120 120 120 60

150 150 150 150 150 150 150 150 150 150 150 150 540

.

Appendix D. Passenger arrival rates and OD matrix Table 4 represents the passenger arrival rates at each station during a given time span, which uses passenger per second as the unit; Table 5 is an OD matrix that describes the proportion of passengers traveling from the origin stations to the destination stations. Table 4 Passenger arrival rates fluctuations with time. Time Station

7:00– 7:15

7:15– 7:30

7:30– 7:45

7:45– 8:00

8:00– 8:15

8:15– 8:30

8:30– 8:45

8:45– 9:00

9:00– 9:15

9:15– 9:30

9:30– 9:45

9:45– 10:00

10:00– 10:15

10:15– 10:30

10:30– 10:45

1 2 3 4 5 6 7 8 9 10 11 12 13

0.35 0.55 0.54 1.3 0.84 0.35 0.303 0.247 1.063 2.617 0.65 0.2 0.907

0.467 0.667 0.6 1.6 1 0.397 0.367 0.337 1.133 2.84 0.747 0.247 0.997

0.603 0.833 0.633 2.02 1.127 0.393 0.437 0.4 1.3 3.073 0.827 0.277 1.197

0.7 0.813 0.7 2.433 1.2 0.433 0.533 0.467 1.493 3.18 0.967 0.19 1.3

0.8 0.9 0.667 2.533 1.25 0.483 0.61 0.573 1.573 3.107 0.94 0.163 1.26

0.883 1.033 0.687 2.43 1.3 0.613 0.737 0.617 1.687 3.15 0.827 0.24 1.343

0.687 1.043 0.633 2.133 1.147 0.55 0.567 0.5 1.6 3.167 0.733 0.293 1.383

0.6 0.867 0.467 2.013 1.133 0.49 0.487 0.433 1.49 3.017 0.657 0.177 1.243

0.35 0.55 0.54 1.3 0.84 0.35 0.303 0.247 1.063 2.617 0.65 0.2 0.907

0.467 0.667 0.6 1.6 1 0.397 0.367 0.337 1.133 2.84 0.747 0.247 0.997

0.603 0.833 0.633 2.02 1.127 0.393 0.437 0.4 1.3 3.073 0.827 0.277 1.197

0.7 0.813 0.7 2.433 1.2 0.433 0.533 0.467 1.493 3.18 0.967 0.19 1.3

0.8 0.9 0.667 2.533 1.25 0.483 0.61 0.573 1.573 3.107 0.94 0.163 1.26

0.883 1.033 0.687 2.43 1.3 0.613 0.737 0.617 1.687 3.15 0.827 0.24 1.343

0.687 1.043 0.633 2.133 1.147 0.55 0.567 0.5 1.6 3.167 0.733 0.293 1.383

112

S. Yang, F. Liao and J. Wu et al. / Transportation Research Part B 133 (2020) 85–113 Table 5 OD matrix.

1 2 3 4 5 6 7 8 9 10 11 12 13

1

2

3

4

5

6

7

8

9

10

11

12

13

0 0.048 0.091 0.13 0.13 0.136 0.136 0.182 0.19 0.2 0.235 0.25 0.235

0 0 0.045 0.043 0.043 0.045 0.091 0.091 0.143 0.1 0.176 0.125 0.118

0.048 0.048 0 0.043 0.043 0.045 0.045 0.045 0.048 0.1 0.059 0.063 0.118

0.048 0.048 0.045 0 0.043 0.045 0.045 0.045 0.048 0.05 0.059 0.063 0.059

0.048 0.048 0.045 0.043 0 0.045 0.045 0.045 0.048 0.05 0.059 0.063 0.059

0.048 0.048 0.045 0.043 0.043 0 0.045 0.045 0.048 0.05 0.059 0.063 0.059

0.048 0.048 0.045 0.043 0.043 0.045 0 0.045 0.048 0.05 0.059 0.063 0.059

0.048 0.048 0.045 0.043 0.043 0.045 0.045 0 0.048 0.05 0.059 0.063 0.059

0.048 0.048 0.045 0.043 0.043 0.045 0.045 0.045 0 0.05 0.059 0.063 0.059

0.19 0.19 0.182 0.174 0.174 0.136 0.136 0.091 0.048 0 0.059 0.063 0.059

0.19 0.19 0.182 0.174 0.174 0.182 0.136 0.136 0.095 0.05 0 0.063 0.059

0.19 0.19 0.182 0.174 0.174 0.182 0.182 0.182 0.19 0.15 0.059 0 0.059

0.048 0.048 0.045 0.043 0.043 0.045 0.045 0.045 0.048 0.1 0.059 0.063 0

References Albrecht, A., Howlett, P., Pudney, P., Vu, X., Zhou, P., 2016a. 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 94, 482–508. Albrecht, A., Howlett, P., Pudney, P., Vu, X., Zhou, P., 2016b. 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 94, 509–538. Albrecht, A., Howlett, P., Pudney, P., Vu, X., Zhou, P., 2018. The two-train separation problem on non-level track-driving strategies that minimize total required tractive energy subject to prescribed section clearance times. Transp. Res. Part B 94, 482–508. 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. Transp. Res. Part E 105, 213–239. Deb, K., Pratap, A., Agarwal, S., Meyarivan, T., 2002. A fast and elitist multi-objective genetic algorithm: nsga-ii. IEEE Trans. Evol. Comput. 6 (2), 182–197. Gao, Y., Kroon, L., Schmidt, M., Yang, L., 2016. Rescheduling a metro line in an over-crowded situation after disruptions. Transp. Res. Part B 93, 425–449. Guo, X., Sun, H., Wu, J., Jin, J., Zhou, J., Gao, Z., 2017. Multiperiod-based timetable optimization for metro transit networks. Transp. Res. Part B 96, 46–67. Gupta, S.D., Tobin, J.K., Pavel, L., 2016. A two-step linear programming model for energy-efficient timetables in metro railway networks. Transp. Res. Part B 93, 57–74. Hartl, R.F., Sethi, S.P., Vickson, R.G., 1995. A survey of the maximum principles for optimal control problems with state constraints. SIAM Rev. 37 (2), 181–218. Howlett, P., 1996. Optimal strategies for the control of a train. Automatica 32 (4), 519–532. Howlett, P., 2016. A new look at the rate of change of energy consumption with respect to journey time on an optimal train journey. Transp. Res. Part B 94, 387–402. Kang, L., Wu, J., Sun, H., Zhu, X., Gao, Z., 2015. A case study on the coordination of last trains for the beijing subway network. Transp. Res. Part B 72, 112–127. Lara, A., Sanchez, G., Coello, C.A.C., Schütze, O., 2010. HCS: a new local search strategy for memetic multiobjective evolutionary algorithms. IEEE Trans. Evol. Comput. 14 (1), 112–132. Li, X., Lo, H.K., 2014a. An energy-efficient scheduling and speed control approach for metro rail operations. Transp. Res. Part B 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 70, 269–284. Liao, F., 2016. Modeling duration choice in space–time multi-state supernetworks for individual activity-travel scheduling. Transp. Res. Part C 69, 16–35. Liao, F., Arentze, T.A., Timmermans, H.J.P., 2010. Supernetwork approach for multi-modal and multi-activity travel planning. Transp. Res. Rec. 2175 (1), 38–46. Liao, F., Arentze, T.A., Timmermans, H.J.P., 2011. Constructing personalized transportation network in multi-state supernetworks: a heuristic approach. Int. J. Geogr. Inf. Sci. 25 (11), 1885–1903. Liao, F., Arentze, T.A., Timmermans, H.J.P., 2013. Incorporating space–time constraints and activity-travel time profiles in a multi-state supernetwork approach to individual activity-travel scheduling. Transp. Res. Part B 55, 41–58. Liao, F., Arentze, T.A., Molin, E.J.E., Bothe, W., Timmermans, H.J.P., 2017. Effects of integrated land-use transport scenarios on travel patterns: a multi-state supernetwork application. Transportation 44, 1–25. Liu, M., Lee, C.Y., Zhang, Z., Chu, C., 2016. Bi-objective optimization for the container terminal integrated planning. Transp. Res. Part B 93, 720–749. Liu, P., Yang, L., Gao, Z., Huang, Y., Li, S., Gao, Y., 2018. Energy-efficient train timetable optimization in the subway system with energy storage devices. IEEE Trans. Intell. Transp. Syst. 16 (3), 1469–1478. Lu, G., Zhou, X., Mahmoudi, M., Shi, T., Peng, Q., 2018. Optimizing resource recharging location-routing plans: a resource-space-time network modeling framework for railway locomotive refueling applications. Comput. Ind. Eng. doi:10.1016/j.cie.2018.03.015. Luan, X., Wang, Y., Schutter, B.D., Meng, L., Lodewijks, G., Corman, F., 2018. Integration of real-time traffic management and train control for rail networks– Part 2: extensions towards energy-efficient train operations. Transp. Res. Part B 115, 72–74. Ning, B., Xu, J., Gao, S., Zhang, L., 2015. An integrated control model for headway regulation and energy saving in urban rail transit. IEEE Trans. Intell. Transp. Syst. 16 (3), 1469–1478. Niu, H., Zhou, X., 2013. Optimizing urban rail timetable under time-dependent demand and oversaturated conditions. Transp. Res. Part C 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 76, 117–135. Pacheco, J., Caballero, R., Laguna, M., Molina, J., 2013. Bi-objective bus routing: an application to school buses in rural areas. Transp. Sci. 47 (3), 397–411. Scheepmaker, G.M., Goverde, M.P., Kroon, L.G., 2017. Review of energy-efficient train control and timetabling. Eur. J. Oper. Res. 257 (2), 355–376. Shang, P., Li, R., Guo, J., Xian, K., Zhou, X., 2019. Integrating lagrangian and eulerian observations for passenger flow state estimation in an urban rail transit network: a space-time-state hyper network-based assignment approach. Transp. Res. Part B 122, 1–19. Sun, H., Wu, J., Ma, H., Yang, X., Gao, Z., 2019. A bi-objective timetable optimization model for urban rail transit based on the time-dependent passenger volume. IEEE Trans. Intell. Transp. Syst. 20 (2), 604–615. Tong, L., Pan, Y., Shang, P., Guo, J., Xian, K., Zhou, X., 2019. Open-source public transportation mobility simulation engine DTALite-S: a discretized space–time network-based modeling framework for bridging multi-agent simulation and optimization. Urban Rail. Transit. 5 (1), 1–16. Van Thielen, S., Corman, F., Vansteenwegen, P., 2018. Considering a dynamic impact zone for real-time railway traffic management. Transp. Res. Part B 111, 39–59.

S. Yang, F. Liao and J. Wu et al. / Transportation Research Part B 133 (2020) 85–113

113

Wang, L., Yang, L., Gao, Z., Huang, Y., 2017. Robust train speed trajectory optimization: A stochastic constrained shortest path approach. Front. Eng. Manag. 4 (4), 408–417. Wang, Y., Tang, T., Ning, B., Van den Boom, T., Schutter, B.D., 2015. Passenger-demands-oriented train scheduling for an urban rail transit network. Transp. Res. Part C 60, 1–23. Yang, S., Wu, J., Sun, H., Yang, X., Gao, Z., Chen, A., 2018a. Bi-objective nonlinear programming with minimum energy consumption and passenger waiting time for metro systems, based on the real-world smart-card data. Transportmetrica B 6 (4), 302–319. Yang, S., Wu, J., Yang, X., Liao, F., Li, D., Wei, Y., 2019a. Analysis of energy consumption reduction in metro systems using rolling stop-skipping patterns. Comput. Ind. Eng. 127, 129–142. Yang, S., Wu, J., Yang, X., Sun, H., Gao, Z., 2018b. Energy-efficient timetable and speed profile optimization with multi-phase speed limits: theoretical analysis and application. Appl. Math. Model. 56, 32–50. Yang, X., Chen, A., Li, X., Ning, B., Tang, T., 2015. An energy-efficient scheduling approach to improve the utilization of regenerative energy for metro systems. Transp. Res. Part C 57, 13–29. Yang, X., Chen, A., Wu, J., Gao, Z., Tang, T., 2019b. An energy-efficient rescheduling approach under delay perturbations for metro systems. Transportmetrica B 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 Trans. Intell. Transp. Syst. 14 (1), 438–447. 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. Ye, H., Liu, R., 2016. A multiphase optimal control method for multi-train control and scheduling on railway lines. Transp. Res. Part B 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: mixed-integer linear programming approaches. Transp. Res. Part B 97, 182–213. Yu, C., Ma, W., Lo, H.K., Yang, X., 2015. Optimization of mid-block pedestrian crossing network with discrete demands. Transp. Res. Part B 73, 103–121. Zhang, J., Liao, F., Arentze, T., Timmermans, H., 2011. A multimodal transport network model for advanced traveler information systems. Procedia Comput. Sci. 5, 912–919. Zhou, X., Zhong, M., 2005. Bicriteria train scheduling for high-speed passenger railroad planning applications. Eur. J. Oper. Res. 167 (3), 752–771.