# Metro Timetabling for Time-Varying Passenger Demand and Congestion at Stations

Authors: Keping Li; Hangfei Huang; Paul Schonfeld
Journal: Journal of Advanced Transportation
 (2018)

---

## Abstract

For the train timetabling problem (TTP) in a metro system, the operator-oriented and passenger-oriented objectives are both important but partly conflicting. This paper aims to minimize both objectives by considering frequency (in the line planning stage) and train cost (in the vehicle scheduling stage). Time-varying passenger demand and train capacity are considered in a nonsmooth, nonconvex programming model, which is transformed into a mixed integer programming model with a discrete time-space graph (DTSG). A novel dwell time determining process considering congestion at stations is proposed, which turns the dwell times into dependent variables. In the solution approach, we decompose the TTP into a subproblem for optimizing segment travel times (OST) and a subproblem for optimizing departure headways from the shunting yard (OH). Branch-and-bound and frequency determining algorithms are designed to solve OST. A novel rolling optimization algorithm is designed to solve OH. The numerical experiments include case studies on a short metro line and Beijing Metro Line 4, as well as sensitivity analyses. The results demonstrate the predictive ability of the model, verify the effectiveness and efficiency of the proposed approach, and provide practical insights for different scenarios, which can be used for decision-making support in daily operations.


---

## Body

## 1. Introduction

With increasing concerns about urban congestion and climate change, urban metro rail transportation receives increasing attention due to its high capacity, punctuality and sustainability. Metro timetabling is the problem of assigning precise utilization times for infrastructure resources to every train in the metro system [1]. In large cities, a metro system is often the busiest public transportation system. To meet travel demand and passenger satisfaction, a timetable should be as compact and flexible as possible. Meanwhile, since many metro systems need government subsidies to cover operating expenses, the planners focus more on operating cost than on passenger factors when designing a timetable [2]. This motivates some recent studies to combine these two conflicting objectives in the train timetabling problem [3, 4].Two types of timetables are used in metro systems: cyclic and noncyclic. In a cyclic timetable, the departure times of trains are scheduled in equally spaced cycles (e.g., a half hour cycle). In a railway system, the cyclic timetables are assumed to be preferred by both operators and passengers because such timetables are easy to operate and remember. However, in a metro system, passengers usually arrive randomly at stations and wait for the next available train, without checking precise timetables beforehand. Besides, as pointed out by Robenek et al. [5], a cyclic timetable provides an inefficient operation as there is a mismatch between the supply (determined by the timetables) and the demand (characterized by the time-varying demand). On the other hand, a noncyclic timetable imposes no special rules on the departure times of trains [6]. The noncyclic timetables are more flexible regarding time-varying passenger demand, especially when the demand is large [2]. Thus, the noncyclic timetabling problem is worth exploring [2, 3, 7].The planning process in public transportation can be generally split into three stages: line planning, timetabling, and vehicle scheduling [8]. The line planning stage determines the route/station layouts and frequency. Timetabling is based on the output of line planning, after which the vehicle scheduling can be designed. It is pointed out by Schöbel [8] that going through all these stages sequentially (i.e., independently) leads to unsatisfactory solutions. Since they are interrelated, important factors and objectives at the line planning stage and the vehicle scheduling stage should be accounted for in the timetabling problem.This paper proposes a rolling optimization algorithm to obtain noncyclic timetables considering time-varying passenger demand and the effects of congestion at stations. It integrates the objectives in the line planning (frequency), timetabling (conflicting objectives including passenger wait/in-vehicle time and energy), and vehicle scheduling (train cost). The main contributions of this paper are as follows.(1) This paper specifically models the dynamic evolution of passenger loads on trains at each station, by considering passenger arrival rates, limited train capacity and actual passenger alighting/boarding rates associated with congestion. This is an extension of existing studies [9]. Particularly, Niu and Zhou [10] and Niu et al. [11] propose the concept of “effective loading time” that represents the actual time interval within which the arriving passengers can board a train. The detailed boarding and alighting processes are not discussed in their model. Wang et al. [3] set a lower bound (which accounts for details including the number of boarding/alighting passengers and their boarding/alighting rates) for dwell times but consider the dwell times as variables to be optimized, rather than parameters that depend on arrival/departure times and passenger demand. Impacts of gradually increasing passenger demand are analyzed in Robenek et al. [2]; however, in modelling they do not consider the effects of passenger congestion.(2) This paper integrates objectives from different planning stages in the train timetabling problem, for which solution approaches are scarce [8]. Particularly, Schöbel [8] considers the line planning, timetabling and vehicle scheduling in an integrated way, with passenger- and cost-oriented objectives. However, since his eigenmodel is general and resulting approaches are generic, important objectives described in this paper are not considered in his work, that is, passenger wait time, dwell time, and train capacity. These objectives are included in [2, 5], whereas the frequency and dwell times are fixed in their studies. Besides, energy is not considered in [2, 5].(3) This paper proposes a novel and effective solution approach, which decomposes the master timetabling model into two subproblems. The first subproblem is solved by branch-and-bound and frequency determining algorithms, and the second one is solved by a rolling optimization algorithm which optimizes the departure headways from the shunting yard (which are equivalent to arrival times of trains at their first station) and considers interdependent variables in each rolling step. The design of decision variables in this paper is similar to [2, 5], but the models and solution approaches are different. Compared to Wang et al. [3], which addresses the similar timetabling problem with traditional approaches such as sequential quadratic programming and a genetic algorithm, the approach proposed in this paper performs better computationally in terms of run time and solution quality.The remainder of this paper is as follows. Section2 reviews relevant studies on the train timetabling problem and summarizes the research gap in the literature. Section 3 first describes the problem, introduces formulation of different cost functions, and constructs the train timetabling model. Then the discrete time-space graph is proposed to transform the model into a mixed integer program (MIP). Section 4 discusses the solution approach for the optimization model. Section 5 presents different numerical experiments and sensitivity analyses on a short metro line and on Beijing Metro Line 4, to demonstrate the predictive ability of the model as well as the effectiveness and efficiency of the proposed approach. Section 6 summarizes the works done in this paper and presents the limitations and potential topics in future studies.


## 2. Literature Review

The train timetabling problem (TTP) aims to schedule trains to transport passengers (or goods) without conflicts, by specifying train arrival and departure times at stations. It is interrelated with other planning stages, that is, the line planning stage and the vehicle scheduling stage. We organize the related literature from the perspectives of operator-oriented objectives, passenger-oriented objectives, and their integration.


### 2.1. Operator-Oriented Objectives

In the literature, operator-oriented objectives can be found only for the noncyclic version of the TTP [2]. Most studies focus on the operation indicators, for example, energy [17, 18], train delay [19–21], and train travel time [22–24].The total train delay and train travel time are usually considered in the TTP for railway networks, which aims to find a feasible timetable by minimizing the profit loss resulting from changes to the ideal/planned timetable. Energy is the focus and main cost for metro system operators. Most energy is consumed in train operations [25, 26], and thus obtaining an energy-efficient train timetable receives considerable attention [17, 18, 27]. Note that energy is affected by the travel time in the segment, and the segment travel time is determined by the arrival and departure times at adjacent stations, that is, the timetable. Hence, the energy-efficient train timetabling problem may be extended to include other objectives. For example, Yang et al. [23] consider both energy and train travel time as optimization objectives. Some recent studies also account for regenerative braking. Yang et al. [28] develop a scheduling approach to coordinate arrivals and departures of trains within the same electricity supply zones, thereby effectively recovering the regenerative energy. Generally, these studies focus solely on the timetabling stage, where frequency (in line planning) is fixed and train cost (in-vehicle scheduling) is not considered.The efficient use of railway rolling stock (vehicle scheduling) is an important objective pursued by a railway agency or company because of intensive capital investment in rolling stock [29]. To this end, Lai et al. [29] develop an optimization model to improve the efficiency of rolling stock usage considering necessary regulations and practical constraints, where a hybrid heuristic process is designed to improve solution quality and efficiency. Haahr et al. [27] use CPLEX and a column and row generation approach to assign rolling stock units to timetable services in passenger railways, prepare daily schedules, and check their real-time applicability by testing different disruption scenarios. These vehicle scheduling studies require the trips as input data. Based on Schöbel [8], in a metro system, every trip describes the operation of a train between the start and end time of the line at its first and its last stations (given from the timetable). The objective function then aims at minimizing the number of trains and the costs for their movements.Particularly, Schmid and Ehmke [30] and Schöbel [8] demonstrate that the integration of timetabling and vehicle scheduling is more beneficial than the sequential planning process. However, the studies on the integration of timetabling and vehicle scheduling usually adopt relatively general and generic models, and the TTP model is often simplified [31]. For example, Cadarso and Marin [32] focus on shunting operations and their timetable is calculated by readjusting frequencies, whereas specific details such as the dwell times and segment travel times are not considered. This motivates us to consider train cost in the timetabling stage, while focusing on the specified TTP formulation.In addition, the conflict-free timetable is also pursued by the operators [33, 34]. It is worth mentioning that the discrete time-space graph proposed in Caprara et al. [6] is a directed multigraph in which nodes correspond to arrivals/departures at a certain station and at a given time stamp. This graph is widely used to formulate the TTP and derive different integer programming models that correspond to specific objectives. For example, Caprara et al. [6] use the graph to derive an integer linear programming (ILP) model with Lagrangian relaxation, which is embedded within a heuristic algorithm. Cacchiani et al. [12] use the graph to formulate an ILP with linear programming (LP) relaxation. Cacchiani et al. [35] use the time-space graph-based LP relaxation of an ILP to derive a dual bound in the TTP for a set of stations in an urban area interconnected by tracks, thus aiming to resolve the conflicts and evaluate the capacity saturation.


### 2.2. Passenger-Oriented Objectives

For a metro system, passenger-oriented timetables that consider reliability and reduction of passenger time are most desirable [11]. In some studies, the passenger demand is considered stable [36]. Assuming that passengers prefer easily memorable timetables, such timetables are usually cyclic. They are designed on the basis of a period event scheduling problem (PESP), aiming to minimize passenger travel and wait time [15, 37] and to maximize network stability [38]. However, as mentioned above, metro passengers arrive randomly and do not remember the timetable. Besides, cyclic timetables are less flexible than noncyclic ones in accounting for the passenger demand [2].Recent studies consider time-varying demands due to the growing concerns for service level and congestion at stations. Such timetables are noncyclic, with objectives of improving passenger satisfaction [3, 7, 10]. Particularly, Barrena et al. [7] propose two nonlinear programming (NP) formulations to generalize noncyclic train timetables on a single line, which are solved by a fast adaptive large neighborhood search metaheuristic. The objective is to minimize passenger wait times at stations. Niu et al. [11] construct mathematically rigorous and algorithmically tractable nonlinear mixed integer programming (MIP) models for both real-time scheduling and medium-term planning applications, to jointly synchronize effective passenger loading time windows and train arrival and departure times at each station. Their work aims to minimize the total waiting times of passengers at stations. Robenek et al. [5] define a timetable as a set of departure times of every train from its origin station on every line and consider four attributes in passenger satisfaction: the in-vehicle-time, the waiting time at transfers, the number of transfers, and the schedule passenger delay. They describe the hybrid timetables through additional constraints that are imposed on the original passenger centric train timetabling problem [2], so that the passengers would obtain the same level of service as a cyclic timetable, with more flexibility. A specifically defined simulated annealing heuristic is proposed to solve the problem. Luan et al. [16] integrate the TTP and preventive maintenance time slots (PMTSs) planning problem on a general railway network. They propose an MILP model to jointly optimize train routes, orders and passing times at each station, and PMTSs. The objective is to minimize the sum of the absolute arrival time deviations of real trains at destinations between the ideal and actual timetables. A Lagrangian-based label correcting algorithm is designed for solving the time-dependent least cost path problem.In addition, the delay management (DM) problem determines whether trains should wait for a delayed feeder train or should depart on time, while considering passenger-oriented objectives [13]. Some recent studies integrate the macroscopic DM with the microscopic TTP [39]. Particularly, Dollevoet et al. [14] formulate an integer programming (IP) model and propose an iterative optimization approach to solve such a bilevel problem that the macroscopic level is the delay management and the microscopic level is the train timetabling. The optimization approach repeats a process that uses DM to find a solution and uses TTP to validate it, until a feasible solution is found. It should be noted that the graph-based models are adopted in their work, which simplify the formulations: the TTP is formulated with the alternative graph [40], and the DM is formulated with an event-activity network [15].


### 2.3. Integrated Objectives

Since the operator- and passenger-oriented objectives are partly conflicting, focusing on only one side would yield undesirable solutions; for example, the best possible service for passengers may also be the most costly alternative for the operator [2]. Thus, many TTP studies integrate these objectives jointly [22, 41].Wang et al. [3] aim to minimize passenger travel time and energy in their nonsmooth nonconvex programming (NSNCP) model. They propose an event-driven model that involves arrival and departure events and passenger arrival rates change events and use nonlinear programming approaches and evolutionary algorithms to solve the problem.Qi et al. [42] focus on a line planning problem with budget constraints and evaluate the station layout with train timetable indicators. The objectives include constructing cost for additional tracks, total travel time, and network capacity. A mixed integer linear programming (MILP) model is formulated to address the TTP, and then a bilevel programming (BP) model is designed to address the integrated problem. Both CPLEX and a local search-based heuristic are developed to solve the models.Robenek et al. [2] model the PCTTP as an MILP model with objectives of maximizing the operator’s profit while maintaining passenger satisfaction. The model uses the output of the line planning problem as the input of the train timetabling problem, to address the gap found between them. CPLEX is used to construct a Pareto Frontier of the contradictory objectives.Li et al. [43] combine dynamic train regulation and passenger flow control to minimize the timetable and the headway deviations for metro lines, thereby reducing the operator profit loss and the passenger delay. The model predictive control and a quadratic programming algorithm are proposed to solve the problem, and stability analysis is conducted to verify the system convergence.However, the above studies do not consider the impacts of frequency. Schöbel [44] reviews some approaches to model and solve the line planning problem, which demonstrate the impacts of frequency to both operator- and passenger-oriented objectives. Particularly, Samanta and Jha [45] consider different objectives in the line planning stage which include minimizing user cost, operator cost, and location cost and maximizing the ridership or the service coverage of the line. The objectives are separately formulated with NP models, and a GA is developed to optimize the objective functions subject to constraints on the number and spacing of stations. Fu et al. [46] propose a BP model to determine line frequencies and the stopping patterns for a mix of fast and slow train lines, with the objectives of minimizing passenger travel time and the number of transfers and maximizing train capacity occupancy. The heuristics solve the problem by combining additional conflicting demand-supply factors. Cascetta and Coppola [47] specify both the frequency-based assignment models and the timetable-based models in the demand forecasting, to analyze their differences in modal split estimates and flows on individual trains. The results show that when passenger demand is time-varying or the timetable is noncyclic, the difference is significant. Thus, the integration of frequency and timetabling is worth investigating.


### 2.4. Research Gap

Table1 highlights some publications related to the TTP. All the graph-based models use integer programs, including the model in this paper. Although different objectives are considered, very few studies integrate the objectives in all three stages. Particularly, Wang et al. [3] only consider the objectives in one stage, that is, timetabling stage; Robenek et al. [2] integrate passenger satisfaction and operating profit (including train cost) but do not consider frequency and energy; Schöbel [8] takes into account objectives in all stages with an iterating framework; however, crucial factors such as passenger wait time, dwell time, and train capacity are missing. Kroon et al. [15] integrate the objectives in the TTP and vehicle scheduling but omit energy.Table 1
Summary of publications related to the TTP.
Work Model Objectives Decision variable Solution approach Caprara et al. [6] ILP Maximize profit Arrival/departure times Lagrangian-based heuristic  Cacchiani et al. [12] ILP Maximize profit A full timetable of a train Column generation, decomposition, and heuristic  Dollevoet et al. [13] IP Minimize total passenger delay If a connection is maintained/used Modified Dijkstra  Niu and Zhou [10] IP Minimize passenger wait time Arrival/departure times GA  Barrena et al. [7] LP Minimize passenger wait time Arrival/departure times Branch and bound  Dollevoet et al. [14] IP Minimize(1) passenger arrival times, (2) maximize consecutive delay Connection and arrival/departure times Iterative optimization approach  Kroon et al. [15] MIP Minimize the number of trains, passenger transfer time, and total travel time Arrival/departure times CPLEX  Niu et al. [11] MIP Minimize passenger wait time Arrival/departure times GAMS  Wang et al. [3] NSNCP Minimize passenger travel time and energy Arrival/departure times Rolling horizon with SQP and GA  Robenek et al. [2] MILP Maximize operating profit while maintaining passenger satisfaction Departure times at first station CPLEX  Robenek et al. [5] MILP Maximize passenger satisfaction Departure times at first station Simulated Annealing (SA) heuristic  Schöbel [8] Eigenmodel Passenger- and operator-oriented objectives in all three planning stages Frequency, A/D times, trip assignment Iterating algorithm  Luan et al. [16] MILP Minimize absolute arrival time deviations A/D times, route selection Lagrangian relaxation  Yin et al. [4] MILP Minimize passenger wait time and energy A/D times, train control Lagrangian-based heuristic  This paper MIP Minimize passenger wait and travel time, energy, and train cost Segment travel times and arrival times at the first station Branch and bound, Rolling optimizationApart from the studies focusing on the delay management, most TTP-related studies optimize arrival and departure times of trains at each station; that is, the dwell times are optimized in the solution approach. However, the dwell times are heavily affected by the number of boarding and alighting passengers and their boarding and alighting rates, associated with congestion at stations. Thus, they should be considered as parameters that are dependent on other decision variables, that is, the departure times from the first station [2, 5] and the segment travel times. In this paper, we account for these characteristics and introduce specific equations to determine the dwell times in different scenarios. Such a process for determining dwell times has not been found in the literature.Finally, heuristics and exact approaches are designed to find solutions based on the properties of the models. Particularly, rolling horizon [3] and iterative optimization [8, 14] are used to solve problems that have similar properties to our proposed model. However, the boundary constraints of the rolling horizon limit the solution flexibility of our model; the iterative optimization considers different stages with same importance, while this paper focuses only on the timetabling stage. Besides, the algorithms used in the rolling horizon and iterative optimization are relatively generic [3, 8]. Thus, a more specific algorithm is desirable for better solutions. In this regard, we develop a novel rolling optimization approach to solve the model, which combines the benefits of rolling horizon and iterative optimization.


## 2.1. Operator-Oriented Objectives

In the literature, operator-oriented objectives can be found only for the noncyclic version of the TTP [2]. Most studies focus on the operation indicators, for example, energy [17, 18], train delay [19–21], and train travel time [22–24].The total train delay and train travel time are usually considered in the TTP for railway networks, which aims to find a feasible timetable by minimizing the profit loss resulting from changes to the ideal/planned timetable. Energy is the focus and main cost for metro system operators. Most energy is consumed in train operations [25, 26], and thus obtaining an energy-efficient train timetable receives considerable attention [17, 18, 27]. Note that energy is affected by the travel time in the segment, and the segment travel time is determined by the arrival and departure times at adjacent stations, that is, the timetable. Hence, the energy-efficient train timetabling problem may be extended to include other objectives. For example, Yang et al. [23] consider both energy and train travel time as optimization objectives. Some recent studies also account for regenerative braking. Yang et al. [28] develop a scheduling approach to coordinate arrivals and departures of trains within the same electricity supply zones, thereby effectively recovering the regenerative energy. Generally, these studies focus solely on the timetabling stage, where frequency (in line planning) is fixed and train cost (in-vehicle scheduling) is not considered.The efficient use of railway rolling stock (vehicle scheduling) is an important objective pursued by a railway agency or company because of intensive capital investment in rolling stock [29]. To this end, Lai et al. [29] develop an optimization model to improve the efficiency of rolling stock usage considering necessary regulations and practical constraints, where a hybrid heuristic process is designed to improve solution quality and efficiency. Haahr et al. [27] use CPLEX and a column and row generation approach to assign rolling stock units to timetable services in passenger railways, prepare daily schedules, and check their real-time applicability by testing different disruption scenarios. These vehicle scheduling studies require the trips as input data. Based on Schöbel [8], in a metro system, every trip describes the operation of a train between the start and end time of the line at its first and its last stations (given from the timetable). The objective function then aims at minimizing the number of trains and the costs for their movements.Particularly, Schmid and Ehmke [30] and Schöbel [8] demonstrate that the integration of timetabling and vehicle scheduling is more beneficial than the sequential planning process. However, the studies on the integration of timetabling and vehicle scheduling usually adopt relatively general and generic models, and the TTP model is often simplified [31]. For example, Cadarso and Marin [32] focus on shunting operations and their timetable is calculated by readjusting frequencies, whereas specific details such as the dwell times and segment travel times are not considered. This motivates us to consider train cost in the timetabling stage, while focusing on the specified TTP formulation.In addition, the conflict-free timetable is also pursued by the operators [33, 34]. It is worth mentioning that the discrete time-space graph proposed in Caprara et al. [6] is a directed multigraph in which nodes correspond to arrivals/departures at a certain station and at a given time stamp. This graph is widely used to formulate the TTP and derive different integer programming models that correspond to specific objectives. For example, Caprara et al. [6] use the graph to derive an integer linear programming (ILP) model with Lagrangian relaxation, which is embedded within a heuristic algorithm. Cacchiani et al. [12] use the graph to formulate an ILP with linear programming (LP) relaxation. Cacchiani et al. [35] use the time-space graph-based LP relaxation of an ILP to derive a dual bound in the TTP for a set of stations in an urban area interconnected by tracks, thus aiming to resolve the conflicts and evaluate the capacity saturation.


## 2.2. Passenger-Oriented Objectives

For a metro system, passenger-oriented timetables that consider reliability and reduction of passenger time are most desirable [11]. In some studies, the passenger demand is considered stable [36]. Assuming that passengers prefer easily memorable timetables, such timetables are usually cyclic. They are designed on the basis of a period event scheduling problem (PESP), aiming to minimize passenger travel and wait time [15, 37] and to maximize network stability [38]. However, as mentioned above, metro passengers arrive randomly and do not remember the timetable. Besides, cyclic timetables are less flexible than noncyclic ones in accounting for the passenger demand [2].Recent studies consider time-varying demands due to the growing concerns for service level and congestion at stations. Such timetables are noncyclic, with objectives of improving passenger satisfaction [3, 7, 10]. Particularly, Barrena et al. [7] propose two nonlinear programming (NP) formulations to generalize noncyclic train timetables on a single line, which are solved by a fast adaptive large neighborhood search metaheuristic. The objective is to minimize passenger wait times at stations. Niu et al. [11] construct mathematically rigorous and algorithmically tractable nonlinear mixed integer programming (MIP) models for both real-time scheduling and medium-term planning applications, to jointly synchronize effective passenger loading time windows and train arrival and departure times at each station. Their work aims to minimize the total waiting times of passengers at stations. Robenek et al. [5] define a timetable as a set of departure times of every train from its origin station on every line and consider four attributes in passenger satisfaction: the in-vehicle-time, the waiting time at transfers, the number of transfers, and the schedule passenger delay. They describe the hybrid timetables through additional constraints that are imposed on the original passenger centric train timetabling problem [2], so that the passengers would obtain the same level of service as a cyclic timetable, with more flexibility. A specifically defined simulated annealing heuristic is proposed to solve the problem. Luan et al. [16] integrate the TTP and preventive maintenance time slots (PMTSs) planning problem on a general railway network. They propose an MILP model to jointly optimize train routes, orders and passing times at each station, and PMTSs. The objective is to minimize the sum of the absolute arrival time deviations of real trains at destinations between the ideal and actual timetables. A Lagrangian-based label correcting algorithm is designed for solving the time-dependent least cost path problem.In addition, the delay management (DM) problem determines whether trains should wait for a delayed feeder train or should depart on time, while considering passenger-oriented objectives [13]. Some recent studies integrate the macroscopic DM with the microscopic TTP [39]. Particularly, Dollevoet et al. [14] formulate an integer programming (IP) model and propose an iterative optimization approach to solve such a bilevel problem that the macroscopic level is the delay management and the microscopic level is the train timetabling. The optimization approach repeats a process that uses DM to find a solution and uses TTP to validate it, until a feasible solution is found. It should be noted that the graph-based models are adopted in their work, which simplify the formulations: the TTP is formulated with the alternative graph [40], and the DM is formulated with an event-activity network [15].


## 2.3. Integrated Objectives

Since the operator- and passenger-oriented objectives are partly conflicting, focusing on only one side would yield undesirable solutions; for example, the best possible service for passengers may also be the most costly alternative for the operator [2]. Thus, many TTP studies integrate these objectives jointly [22, 41].Wang et al. [3] aim to minimize passenger travel time and energy in their nonsmooth nonconvex programming (NSNCP) model. They propose an event-driven model that involves arrival and departure events and passenger arrival rates change events and use nonlinear programming approaches and evolutionary algorithms to solve the problem.Qi et al. [42] focus on a line planning problem with budget constraints and evaluate the station layout with train timetable indicators. The objectives include constructing cost for additional tracks, total travel time, and network capacity. A mixed integer linear programming (MILP) model is formulated to address the TTP, and then a bilevel programming (BP) model is designed to address the integrated problem. Both CPLEX and a local search-based heuristic are developed to solve the models.Robenek et al. [2] model the PCTTP as an MILP model with objectives of maximizing the operator’s profit while maintaining passenger satisfaction. The model uses the output of the line planning problem as the input of the train timetabling problem, to address the gap found between them. CPLEX is used to construct a Pareto Frontier of the contradictory objectives.Li et al. [43] combine dynamic train regulation and passenger flow control to minimize the timetable and the headway deviations for metro lines, thereby reducing the operator profit loss and the passenger delay. The model predictive control and a quadratic programming algorithm are proposed to solve the problem, and stability analysis is conducted to verify the system convergence.However, the above studies do not consider the impacts of frequency. Schöbel [44] reviews some approaches to model and solve the line planning problem, which demonstrate the impacts of frequency to both operator- and passenger-oriented objectives. Particularly, Samanta and Jha [45] consider different objectives in the line planning stage which include minimizing user cost, operator cost, and location cost and maximizing the ridership or the service coverage of the line. The objectives are separately formulated with NP models, and a GA is developed to optimize the objective functions subject to constraints on the number and spacing of stations. Fu et al. [46] propose a BP model to determine line frequencies and the stopping patterns for a mix of fast and slow train lines, with the objectives of minimizing passenger travel time and the number of transfers and maximizing train capacity occupancy. The heuristics solve the problem by combining additional conflicting demand-supply factors. Cascetta and Coppola [47] specify both the frequency-based assignment models and the timetable-based models in the demand forecasting, to analyze their differences in modal split estimates and flows on individual trains. The results show that when passenger demand is time-varying or the timetable is noncyclic, the difference is significant. Thus, the integration of frequency and timetabling is worth investigating.


## 2.4. Research Gap

Table1 highlights some publications related to the TTP. All the graph-based models use integer programs, including the model in this paper. Although different objectives are considered, very few studies integrate the objectives in all three stages. Particularly, Wang et al. [3] only consider the objectives in one stage, that is, timetabling stage; Robenek et al. [2] integrate passenger satisfaction and operating profit (including train cost) but do not consider frequency and energy; Schöbel [8] takes into account objectives in all stages with an iterating framework; however, crucial factors such as passenger wait time, dwell time, and train capacity are missing. Kroon et al. [15] integrate the objectives in the TTP and vehicle scheduling but omit energy.Table 1
Summary of publications related to the TTP.
Work Model Objectives Decision variable Solution approach Caprara et al. [6] ILP Maximize profit Arrival/departure times Lagrangian-based heuristic  Cacchiani et al. [12] ILP Maximize profit A full timetable of a train Column generation, decomposition, and heuristic  Dollevoet et al. [13] IP Minimize total passenger delay If a connection is maintained/used Modified Dijkstra  Niu and Zhou [10] IP Minimize passenger wait time Arrival/departure times GA  Barrena et al. [7] LP Minimize passenger wait time Arrival/departure times Branch and bound  Dollevoet et al. [14] IP Minimize(1) passenger arrival times, (2) maximize consecutive delay Connection and arrival/departure times Iterative optimization approach  Kroon et al. [15] MIP Minimize the number of trains, passenger transfer time, and total travel time Arrival/departure times CPLEX  Niu et al. [11] MIP Minimize passenger wait time Arrival/departure times GAMS  Wang et al. [3] NSNCP Minimize passenger travel time and energy Arrival/departure times Rolling horizon with SQP and GA  Robenek et al. [2] MILP Maximize operating profit while maintaining passenger satisfaction Departure times at first station CPLEX  Robenek et al. [5] MILP Maximize passenger satisfaction Departure times at first station Simulated Annealing (SA) heuristic  Schöbel [8] Eigenmodel Passenger- and operator-oriented objectives in all three planning stages Frequency, A/D times, trip assignment Iterating algorithm  Luan et al. [16] MILP Minimize absolute arrival time deviations A/D times, route selection Lagrangian relaxation  Yin et al. [4] MILP Minimize passenger wait time and energy A/D times, train control Lagrangian-based heuristic  This paper MIP Minimize passenger wait and travel time, energy, and train cost Segment travel times and arrival times at the first station Branch and bound, Rolling optimizationApart from the studies focusing on the delay management, most TTP-related studies optimize arrival and departure times of trains at each station; that is, the dwell times are optimized in the solution approach. However, the dwell times are heavily affected by the number of boarding and alighting passengers and their boarding and alighting rates, associated with congestion at stations. Thus, they should be considered as parameters that are dependent on other decision variables, that is, the departure times from the first station [2, 5] and the segment travel times. In this paper, we account for these characteristics and introduce specific equations to determine the dwell times in different scenarios. Such a process for determining dwell times has not been found in the literature.Finally, heuristics and exact approaches are designed to find solutions based on the properties of the models. Particularly, rolling horizon [3] and iterative optimization [8, 14] are used to solve problems that have similar properties to our proposed model. However, the boundary constraints of the rolling horizon limit the solution flexibility of our model; the iterative optimization considers different stages with same importance, while this paper focuses only on the timetabling stage. Besides, the algorithms used in the rolling horizon and iterative optimization are relatively generic [3, 8]. Thus, a more specific algorithm is desirable for better solutions. In this regard, we develop a novel rolling optimization approach to solve the model, which combines the benefits of rolling horizon and iterative optimization.


## 3. Model Formulation

In this section, we present a mixed integer programming formulation for the TTP. It is a multiobjective optimization problem. The weighted-sum scalarization of different objective functions in terms of generalized monetary cost is used to resolve the tradeoffs among them [8].


### 3.1. Problem Formulation

Table2 summarizes the parameters and variables used in this paper. Some are also explained further in the text. The model’s simplifying assumptions are listed below and also explained in the problem description.Table 2
Parameters and variables.
Name Description Units Purpose I Set of stations - Set K Set of trains - Set U i Set of time stamps in which trains can arrive at stationi - Set W i Set of time stamps in which trains departure from stationi - Set  θ Shunting yard of the metro line - Symbol i ,  j Station index - Symbol I Total number of stations - Symbol k Train index - Symbol K Total number of trains within a given period/Frequency - Symbol t Time index - Symbol u k Arrival node of traink - Symbol w k Departure node of traink - Symbol π ( u k ) Time stamp associated with nodeuk (or wk) - Symbol Δ ( u k ,  w k ′ ) Time duration between two nodesuk and wk′ (k=k′ is allowed) - Symbol R k Path of traink - Symbol T Time duration - Symbol a The largest integer that is smaller than real numbera - Symbol a The smallest integer that is larger than real numbera - Symbol b ⃑ A vector (represented by any bold alphabet) - Symbol M A matrix (represented by any double-struck alphabet) - Symbol M a : b ,  x : y A sub-matrix whose elements correspond to those in matrixM, from rows a to b, and from columns from x to y - Symbol M Sum of all (numerical) elements in matrixM - Symbol  c e Average energy cost ¥/Kwh Parameter c o Average train operating cost ¥ /(train⋅km) Parameter c p t Average passenger travel time cost ¥/(passenger⋅hour) Parameter c p w Average passenger wait time cost ¥/(passenger⋅hour) Parameter c t r a Average capital cost of one train ¥/(train⋅hour) Parameter C t o t a l Total cost of the whole system ¥ Parameter C t o t a l e Total cost of energy consumption ¥ Parameter C t o t a l p Total cost of passenger time ¥ Parameter C t o t a l t r a Total cost of train (rolling stock) ¥ Parameter d i ¯ Maximum allowable dwell time at stationi Second Parameter d i _ Minimum allowable dwell time at stationi Second Parameter e i Energy consumption per unit mass for one train in segmenti Kwh/Kg Parameter f t Tractive force per unit mass m / s 2 Parameter f r Running resistance per unit mass m / s 2 Parameter F s Fleet size (the number of trains) Trains Parameter h m i n Minimum headway Second Parameter M k Net weight of a train Kg/train Parameter M p Average weight of a passenger Kg/passenger Parameter n d Number of doors per train - Parameter N k Capacity of traink Passengers Parameter P k ,  i a Number of passengers alighting traink at station i Passengers Parameter P k ,  i b Number of passengers boarding traink at station i Passengers Parameter P k ,  ( i ,  j ) b Number of passengers boarding traink at station i with destination j Passengers Parameter P i ( t 1 ,  t 2 ) Number of passengers arriving at stationi during [t1,t2] Passengers Parameter P ( i ,  j ) t 1 ,  t 2 Number of passengers arriving at stationi with destination j during [t1,t2] Passengers Parameter P i ,  t w Number of passengers waiting at stationi at time t Passengers Parameter P k ,  t v Number of passengers in traink at time t Passengers Parameter s l ( i ) Length of segmenti (between station i and i+1) Meter Parameter s t m a x ( i ) Maximum segment travel time for segmenti Second Parameter s t m i n ( i ) Minimum segment travel time for segmenti Second Parameter t k ,  i 0 Arrival time of traink at station i Second Parameter t k ,  i a Alighting end time for traink at station i Second Parameter t k ,  i b Boarding end time for traink at station i/Departure time Second Parameter t e n d p End time of passenger demand for a period Second Parameter t e n d t End time of a complete timetable for a period Second Parameter T r o u n d k Round trip time (total travel time from the departure from the shunting yard to the arrival at the shunting yard) of traink Second Parameter T p t Total passenger travel time Second Parameter T p w Total passenger wait time Second Parameter T t u r n Turnaround time Second Parameter v k ( t ) Speed of traink at time t m/s Parameter x k ( t ) Position of traink at time t Meter Parameter τ i ,  j ( t ) Passenger arrival rate at stationi with destination j at time t Passengers/second Parameter λ i ( t ) Passenger arrival rate at stationi at time t Passengers/second Parameter ζ a Passenger alighting rate Passengers/second Parameter ζ b Passenger boarding rate Passengers/second Parameter χ x Parameter for energy approximation - Parameter χ y Parameter for energy approximation - Parameter  s t ( i ) Segment travel time of segmenti Second Decision t k ,  1 0 Arrival time of traink at the first station Second DecisionAssumption 1.
Within one period, trains travel through each segment with the same speed profile; that is, the travel times for different trains in a segment are equal.Assumption 2.
We assume that the last train arrives at the first station at the end time of the period, to ensure that no late-arriving passengers will be left behind.Assumption 3.
Passenger arrival rateλi(t) should not exceed total passenger boarding rate ζb·nd. If λit>ζb·nd, the queue of waiting passengers increases even if a train is dwelling at the station. In practice, passenger entrance flow control∗ is taken to avoid such a case.R e m a r k ∗. At most busy metro stations in China, passenger entrance flow control is adopted during peak hours to ensure operation safety, and thus the passenger arrival rate is a restricted number, which is expected not to exceed passenger boarding rate. Assumption 3 is introduced for consistency with practical experience.Assumption 4.
All passengers in traink who have reached their destination are assumed to exit the train at its arrival time of station i, that is, tk,i0. All waiting passengers at station i able to board train k are assumed to board the train at its departure time, that is, tk,ib. In other words, the number of in-vehicle passengers Pk,tv changes only at the arrival/departure times of train k at stations.This paper focuses on a bidirectional metro line with a single track, as shown in Figure1. Properties of such a metro line are described in [48]. The timetabling problem is studied within a given period. The time-varying passenger demand is given as input, and a series of trains (whose frequency is to be optimized) are scheduled to transport the passengers.Figure 1
A bidirectional metro line with 2I stations.The arrival times at the first station are decision variables that must be optimized. The departure times at the first station are determined by the dwell times. As mentioned earlier, this paper obtains the dwell times based on trains’ conditions when they arrive at stations. In this regard, optimizing arrival times at the first station is equivalent to optimizing departure times at the first station, as in [2, 5]. For the second station, the arrival times are determined by the segment travel time in the first segment. Practically, for the metro operators, uniform segment travel times are preferred within one period, because they not only make the system more stable and easier to manage but also simplify synchronized passenger transfers with connecting transit lines. Assumption 1 is introduced to consider this preference. Thus, within one period, the segment travel time in each segment needs to be optimized only once (for all trains within this period). The speed profile optimization is outside the scope of this paper. The departure times at the second station are also determined by the dwell time determining process.Analogously, for subsequent train arrivals and departures, the arrival times are determined by the previous segment travel time, and the departure times are determined by the dwell times. Thus, within a period, the arrival times at the first station and the segment travel times determine the train timetable.To represent the TTP, a directed multigraph called discrete time-space graph (DTSG) is proposed [6]. As shown in Figure 2, θ represents the shunting yard of the metro line and Ui  (Wi) corresponds to the set of time stamps in which trains can arrive at (depart from) station i. The nodes in Ui  (Wi) are called arrival (departure) nodes. Unlike the original graph that is designed for a railway network [6, 12], this paper considers a single-track metro line, and thus the constraints are different and implicitly imposed in the definition of the structure of the DTSG.Figure 2
Representation of discrete time-space graph.The continuous time is discretized into evenly spaced time intervals (1 second in this paper), and the period extends from 0 totendt. Note that the values of tendt and tendp are different. tendp is the last time stamp of the passenger demand. The period mentioned above refers to the one extending from 0 to tendp. Assumption 2 imposes that the last train departs from the first station at time tendp. However, if the DTSG only considers the period from 0 to tendp, the arrival/departures of the last train at subsequent stations are not considered, and the passenger-oriented objectives are hard to formulate. Thus, the DTSG considers the period from 0 to tendt. The relation between tendt and tendp is set as:(1)tendt=tendp+∑i=12I-1sti+∑i=12Idi¯,where st(i) is the segment travel time of segment i and di¯ is the maximum dwell time of station i. This setting ensures that even the slowest train can be completely included in the DTSG.The construction of the DTSG is noncyclic. The objective functions of passenger time, energy, and train cost are associated with the chosen arcs and are described below.


### 3.2. Passenger Time

#### 3.2.1. Passenger Arriving, Alighting, and Boarding

(1) Time-Varying Demand Information. The passenger demand at each station can be expressed by τi,j(t), which is defined as(2)τi,jt≥0,If1≤i<j≤IorI≤i<j≤2I;τi,jt=0,Otherwise;∀i,j∈I.Passenger arrival rates at stations areλit=∑j=i+12Iτi,j(t), ∀i,j∈I. In this paper, demands are regarded as deterministic values, which are given by τi,j(t). Thus, the number of arriving passengers at station i with destination j during time interval [t1,t2] is(3)Pi,jt1,t2=∑t=t1t2τi,jt,∀i,j∈I.The number of arriving passengers at stationi within time interval [t1,t2] is the sum of arriving passengers with different destinations:(4)Pit1,t2=∑j=i+12IPi,jt1,t2=∑j=i+12I∑t=t1t2τi,jt=∑t=t1t2λit,∀i,j∈I.(2) Determination of Dwell Arcs. Letting nd be the number of doors per train, the average number of alighting (boarding) passengers per door when train k arrives at (departs from) station i is Pk,ia/nd  (Pk,ib/nd). In addition, we denote the average number of boarding passengers per second (boarding rate) as ζb and alighting rate as ζa, which vary with passenger density in trains. For simplicity, we consider that boarding passengers wait until all alighting passengers get off the train; that is, the boarding process starts at time tk,ia.The dwell time is the sum of alighting and boarding times. The calculation of alighting time is solely based on the number of alighting passengers:(5)tk,ia=tk,i0+Pk,iand·ζa,∀k∈K,i∈I.The calculation of boarding time is slightly more complex. Figure3 shows the relations of train remaining capacity, waiting passengers, boarding rate, and max/min dwell times. The horizontal axis represents time, and the vertical axis represents the number of passengers. These four factors determine the boarding time for train k at station i, and in Figure 3 they are represented by (1) remaining capacity of train k at station i: Nk-Pk,tk,iav, a horizontal line; (2) number of waiting passengers at station i at time t, a curve Pi,tw; (3) maximum allowable number of boarding passengers, a diagonal function fkt with the slope of nd·ζb; (4) minimum and maximum dwell times at station i, that is, tk,i0+di_ and tk,i0+di¯. In Figure 3, the diagonal P=fkt intersects the horizontal line P=Nk-Pk,tk,iav at time t1k,i and intersects the curve P=Pi,tw at time t2k,i. The curve P=Pi,tw intersects the horizontal line P=Nk-Pk,tk,iav at time t3k,i. Note that the time is discrete, and Assumption 3 ensures that the slope of P=fkt always exceeds the slope of P=Pi,tw. Thus, we have(6)t1k,i=max⁡t∣t-tk,ia·nd·ζb≤Nk-Pk,tk,iav,∀k∈K,i∈I,(7)t2k,i=maxt∣t-tk,ia·nd·ζb≤Pi,tw,∀k∈K,i∈I,(8)t3k,i=max⁡t∣Pi,tw≤Nk-Pk,tk,iav,∀k∈K,i∈I.Figure 3
Relations of remaining capacity, waiting passengers, boarding rate, and max/min dwell times.Practically, in an oversaturated scenario, a train departs from a station as soon as the number of in-train passengers reaches its capacity. Based on this strategy, (6) shows that the boarding passengers who arrive after t1k,i exceed train capacity, and they will be left behind; (7) shows that before t2k,i there are passengers waiting at station, and after t2k,i all passengers board the train; (8) points that the train is full at t3k,i, and passengers arriving after t3k,i should wait for the next train.The departure time of traink at station i  (tk,ib) is determined based on the relations among t1k,i, t2k,i, t3k,i, tk,i0+di_, and tk,i0+di¯. Based on Assumption 3, we consider two cases: (1)  t1k,i>t2k,i and (2)  t1k,i≤t2k,i. Note that di_ and di¯ should be considered in both cases. In case (1), we can further derive that t3k,i>t1k,i>t2k,i (which is the solid curve in Figure 3), where the number of boarding passengers per second is determined by passenger arrival rate. Thus, the departure time is determined as tk,ib=min{tk,i0+di¯,max⁡t2k,i,tk,i0+di_,t3k,i}. In case (2), we can further derive that t3k,i≤t1k,i≤t2k,i (which is the dashed curve in Figure 3), where the number of boarding passengers per second is determined by the actual passenger boarding rate. Thus, the departure time is determined as tk,ib=min{tk,i0+di¯,t1k,i}. The departure time of train k at station i is expressed as(9)tk,ib=min⁡tk,i0+di¯,max⁡t2k,i,tk,i0+di_,t3k,i,t1k,i>t2k,imin⁡tk,i0+di¯,t1k,i,t1k,i≤t2k,i.Finally, we consider the relations between alighting/boarding rates and congestion at stations. Here the congestion is described by the number of in-vehicle passengers. Practically, the greater the passenger density in a train is, the slower the passengers can get off/on the train. Here, we design a monotonically decreasing functionρ(Pk,tv)≤1 with respect to the number of in-vehicle passengers, to approximate actual alighting/boarding rates according to the passenger density in trains. For simplicity, in this paper, the value of ρ(Pk,tv) during the dwell process is fixed and is determined by the number of in-vehicle passengers at its arrival time. That is, ζa and ζb remain constant during the dwell process:(10)ζa=ρPk,tk,i0v·ζaidealζb=ρPk,tk,i0v·ζbideal;∀k∈K,i∈I,where ζaideal and ζbideal are the ideal alighting and boarding rates in the uncongested condition.(3) Passenger Alighting and Boarding. Once the dwell arcs are determined, the alighting and boarding times for train k at station i are determined. Based on the results of dwell arcs, three cases should be considered: (1)  t1k,i<t2k,i; (2)  t1k,i≥t2k,i and t2k,i≤tk,i0+di_; (3)  t1k,i≥t2k,i and t2k,i>tk,i0+di_. In case (1), waiting passengers board the train at rate ζb. In case (2), actual passenger boarding rate is equal to the passenger arrival rate. In case (3), the actual boarding rate is ζb. Thus, the number of boarding passengers for train k at station i is(11)Pk,ib=Pi,tk,ibw;ift1k,i≥t2k,i,t2k,i≤tk,i0+di_;tk,ib-tk,ia·ζb·nd;otherwise;∀k∈K,i∈I.Wang et al. [3] point out that, in a busy metro line, passengers with different destinations are well mixed at each station. Thus, the ratio between the numbers of boarding passengers with destination j and of total boarding passengers is equivalent to the ratio between τi,j(t) and λi(t):(12)Pk,i,jb=∑t=tk-1,ibtk,ibτi,jt∑t=tk-1,ibtk,ibλit·Pk,ib;∀k∈K,j∈I\1,i=1,2,…,j-1.The number of alighting passengers from traink at station j consists of all previously boarding passengers in train k with destination j:(13)Pk,ja=∑i=1j-1Pk,i,jb;∀k∈K,j=2,…,I,I+2,…,2I.There are no alighting passengers at the first station in each direction (station 1 andI+1).


#### 3.2.2. Passenger Wait and Travel Time

(1) Passenger Wait Time. Total passenger wait time at a station can be calculated by summing up the products of the cumulative number of passengers and their time. Let klast be the index of last train that has departed from station i before time t, that is, tklast,ib≤t<tklast+1,ib. The cumulative number of passengers at station i at time t is the difference between cumulative numbers of arrival passengers and boarding passengers within the period [0,t]:(14)Pi,tw=Pi0,t-∑k=1klastPk,ib;∀i∈I.Since the continuous time is discretized into 1 second intervals, the passenger wait time between timest and t+1 is Pi,tw·t+1-t=Pi,tw. Thus, the total passenger wait time is(15)Tpw=∑i∈I∑t=0tendtPi,tw.(2) Passenger Travel Time. The number of in-train passengers in train k at time t is the difference between cumulative numbers of boarding passengers and alighting passengers. Let ilast be the index of the last station that train k has visited before time t. We consider two cases: (1) train k is dwelling at station ilast, that is, tk,ilast0≤t<tk,ilastb; (2) train k is in the segment between station ilast and ilast+1, that is, tk,ilastb≤t<tk,ilast+10. Based on Assumption 4, the boarding passengers at station ilast are not considered in case (1) but are included in case (2). Thus, the number of in-vehicle passengers in train k at time t is(16)Pk,tv=∑i=1ilast-1Pk,ib-∑i=1ilastPk,ia;tk,ilast0≤t<tk,ilastb;∑i=1ilastPk,ib-Pk,ia;tk,ilastb≤t<tk,ilast+10;∀k∈K.Analogously, the total passenger travel time is(17)Tpt=∑k∈K∑t=0tendtPk,tv.Finally, the total cost of passenger time is formulated as(18)Ctotalp=cpw·Tpw3600+cpt·Tpt3600.The cost unit is yuan(¥).


### 3.3. Energy

In general, for each segmenti, with given segment travel time st(i) and segment length sl(i), the speed profile determines energy use [49]. Thus, given a speed profile generating method, the energy per unit mass of different trains in segment i is the same. We denote the energy per unit mass for one train in segment i as ei. Its calculation usually considers the following train dynamic equations [48]:(19)dvktdt=ft-fr,dxktdt=vkt;∀t∈tk,ib,tk,i+10,i,i+1∈I,fr=α·vk2t+β·vkt+γ,vktk,i0=vktk,ib=0.The first two equations are train motion equations, wherevkt and xk(t), respectively, represent traveling speed and distance. The third one is the Davis resistance function where Davis parameters are α, β, and γ. The last one requires that every train stops at each station. An optimal train control strategy can be obtained based on (19), in which energy per unit mass is(20)ei=∫tk,ibtk,i+10ft·vktdt;∀i,i+1∈I,k∈K.However, in (20), the relation between energy and decision variable st(i) is unclear. It is also difficult to integrate the optimal train control with other objective functions [18]. Since the speed profile optimization is not considered, for simplicity, we use the method in [50] to generate practical speed profiles. Based on these speed profiles, we then adopt a linear piecewise curve to approximate the energy per unit mass of train k in segment i [48]:(21)ei=χx·sti+χy·sli;∀i∈I.We introduce here an estimation procedure for the values ofχx and χy, which is based on mean squared error (MSE). First, we set several feasible values (e.g., z=1,2,…,100) of slz and st(z). Then, we obtain optimal speed profiles based on (19) and (20), thereby obtaining the optimal energy with respect to slz and st(z). Finally, we define the optimal energy per unit mass as ez, and the estimated energy per unit mass obtained with (21) as ez~. The MSE is formulated as(22)MSEχx,χy=1100∑z=1100ez~-ez2.The optimal values ofχx and χy can be obtained when MSEχx,χy reaches its minimum and are found as follows:(23)∂MSEχx,χy∂χx=0,∂MSEχx,χy∂χy=0.Withχx∗ and χy∗ obtained, the total energy can be calculated as(24)Etotal=∑i=12I-1∑k=1Kχx∗·sti+χy∗·sli·Mk+Mp·Pk,tk,ibv.Thus, the total cost of energy consumption is(25)Ctotale=ce·Etotal=ce·∑i=12I-1∑k=1Kχx∗·sti+χy∗·sli·Mk+Mp·Pk,tk,ibv.Note that when all tracks are straight and level,χx and χy suffice for estimating energy. However, if grades and curves are considered, we may have ei=χxi·st(i)+χyi·sli,∀i∈I, where χxi and χyi are specified parameters to describe the characteristics in segment i. The parameter estimation procedure is the same.


### 3.4. Train Cost

In this paper, the train cost is estimated with empirical equations with respect to the interest rate, train price, and average round trip time, where long-term requirements (e.g., the consideration of the highest peak during a year and demand evolution over the planning horizon) are not considered.Train cost consists of train capital cost and operating cost (excluding energy). Train capital cost is related to the price and economic life of a train, as well as interest rate. LetCannual be the annual cost per train, Ctrain be the current price per train, r be the interest rate, and Tls be the lifespan of a train. The relations among these four factors can be expressed as(26)Cannual=Ctrain·r·1+rTls1+rTls-1.The capital cost per train hourctra can be then calculated by dividing the annual cost Cannual by the product of average working hours per day (Dwork) and working days per year (Ywork)(27)ctra=CannualDwork·Ywork.In addition, a train can be reassigned after it finishes a round trip and returns to the shunting yard, and the round trip timeTroundk can be obtained as(28)Troundk=tk,2Ib-tk,10+Tturn,∀k∈K.Thus, the actual period for train running is∑k=1KTroundk  (∑k=1KTroundk≤tendt). The estimated fleet size Fs within period [0,∑k=1KTroundk] depends on the round trip time Troundk and the departure headway from the shunting yard (which is equivalent to the arrival time of train k at the first station tk,10):(29)Fs=∑k=1KTroundk/K∑k=2Ktk,10-tk-1,10/K-1,where the term (∑k=1KTroundk)/K represents the average round trip time and the term ∑k=2K(tk,10-tk-1,10)/(K-1) represents the average headway. The total cost of trains is formulated as(30)Ctotaltra=ctra·Fs·∑k=1KTroundk3600+co·K·∑i=12I-1sli1000.The first termctra·Fs·(∑k=1KTroundk/3600) represents the capital cost of Fs trains for the actual running period [0,∑k=1KTroundk]. The second term co·K·∑i=12I-1sli represents the operating cost with frequency K (trains) on a line with the length of the sum of 2I-1 segments (Km). Both decision variables (segment travel times and arrival times at the first station) affect the train cost.


### 3.5. Mixed Integer Programming Formulation

The optimization model aims to minimize the weighted-sum of different cost functions, so that the conflicting objectives are generalized in terms of monetary cost. In this paper, the weights of different functions are reflected by their average costs, that is,ce, co, cpt, cpw, and ctra. For example, an operator-oriented timetable (neglecting passenger satisfaction) can be obtained by setting cpt=0 and cpw=0; a passenger-oriented timetable can be obtained by setting ce=co=ctra=0. The weights for integrated objectives can be set similarly, with larger cost coefficients for larger weights. The model is formulated as follows:Minimize(31)Ctotal=Ctotalp+Ctotale+Ctotaltra.Subject to(32)tK,10=tendp,(33)stI=Tturn,(34)stmini≤sti≤stmaxi;∀i∈I,(35)tk,ib-tk-1,ib≥hmin;∀k=2,…,K,i∈I,(36)Pk,tv≤Nk;∀k∈K,t∈0,tendt.This model is a nonsmooth, nonconvex programming model, where the nonsmoothness is caused by the limited train capacity, and the nonconvexity is caused by the determination of dwell times [3]. Constraint (32) is equivalent to Assumption 2. Constraint (33) considers the turnaround from station I to I+1 as traveling in segment I, with fixed travel time Tturn. Constraint (34) ensures that the segment travel time should be optimized within a reasonable range [stmini,stmaxi]. Constraint (35) imposes the minimum headway constraint for both departures and arrivals (based on Assumption 1). Constraint (36) requires that the number of in-train passengers must not exceed train capacity.In Caprara et al. [6], the time-space graph is used to derive an integer linear programming model, thereby allowing a considerably faster solution. Analogously, this paper transforms the nonsmooth, nonconvex programming model into a mixed integer programming (MIP) model based on the DTSG. First, we construct two (2I×tendt) matrixes U and W, respectively, representing the selection of arrival and departure nodes in the DTSG:(37)Ui,t=1;ift=πuk,uk∈Ui,k∈K;0;otherwise;Wi,t=1;ift=πwk,wk∈Wi,k∈K;0;otherwise;∀i∈I,t∈0,tendt,where Ui,t=1  (Wi,t=1) if a train arrives at (departs from) station i at time t. As mentioned earlier, the value of U1,t is determined by πuk=tk,10, for ∀uk∈U1,  k∈K. The value of Wi,t is determined by Ui,t and the dwell arcs, that is, Δuk,wk=tk,ib-tk,i0, for ∀uk∈Ui,  wk∈Wi,  i∈I,  k∈K. The value of Ui+1,t is determined by Wi,t and the travel time in segment i, that is, Δwk,uk=st(i), for ∀wk∈Wi,  uk∈Ui+1,  and i∈I\{2I},  k∈K. Thus, the decision variables determine U and W.Based on the structures ofU and W, as well as Assumption 4, we introduce two (2I×tendt) matrixes PA and PB, respectively, representing passenger alighting and boarding information (number and time):(38)PAi,t=Pk,ia;t=πuk,uk∈Ui,k∈K;0;otherwise;PBi,t=Pk,ib;t=πwk,wk∈Wi,k∈K;0;otherwise;∀i∈I,t∈0,tendt.The cumulative numbers of waiting passengersQ  (2I×tendt) and in-vehicle passengers V  (K×tendt) can be constructed based on PA and PB:(39)Qi,t=Qi,t-PBi,1:t;∀i∈I,t∈0,tendt,(40)Vk,t=Vk,t-1+Pk,ib;ift=πwk,wk∈Wi,i∈I;Vk,t-1-Pk,ia;ift=πuk,uk∈Ui,i∈I;Vk,t-1;otherwise;∀k∈K,t∈1,tendt.In (39), Qi,t=Qi,0+∑q=0tλi(t) represents the total waiting passengers at station i, where Qi,0 is the number of waiting passengers at the start of this period. According to the definition of cumulative passengers, we have Qi,t=Pi,tw. Thus, the total passenger wait time can be expressed by Q:Tpw=∑i∈I∑t=0tendtPi,tw=∑i∈I∑t=0tendtQi,t=Q  (since Qi,t≥0). In (40), the initial number of in-vehicle passengers is zero, that is, Vk,0=0. The number of in-vehicle passengers changes only at arrival/departure times, which is consistent with Assumption 4. Analogously, the total passenger travel time can be expressed as Tpt=V. Thus, the total cost for passenger time can be rewritten as(41)Ctotalp=13600cpw·Q+cpt·V.For energy cost, the energy per unit mass in segmenti can be rewritten as(42)ei=χx·Δwk,uk+χy·sli,wk∈Wi,uk∈Ui+1;∀i∈I\2I.The total energy can be rewritten as(43)Etotal=∑i=12I-1∑k=1K∑t=0tendtϕt·χx·Δwk,uk+χy·sli·Mk+Mp·Vk,t,wk∈Wi,uk∈Ui+1,where ϕt=1 if t=πwk; otherwise ϕt=0. The total energy cost can be rewritten correspondingly.For train cost, the round trip timeTroundk can be rewritten as(44)Troundk=Δuk,wk+Tturn,uk∈U1,wk∈W2I;∀k∈K.The departure headway from the shunting yard can be rewritten as(45)tk,10-tk-1,10=Δuk-1,uk,uk-1,uk∈U1;∀k∈K\1.The total cost of trains can be rewritten correspondingly.Finally, the last departing train constraint (32) is equivalent to (46)πuK=tendt,uK∈U1.The turnaround constraint (33) is equivalent to(47)Δwk,uk=Tturn,wk∈WI,uk∈UI+1;∀k∈K.The segment travel time window constraint (34) is equivalent to(48)stmini≤Δwk,uk≤stmaxi,wk∈Wi-1,uk∈Ui;∀k∈K,i∈I\1.The minimum headway constraint (35) is equivalent to(49)Δwk-1,wk≥hmin,wk-1,wk∈Wi;k∈K\1,i∈I.The train capacity constraint (36) is equivalent to(50)Vk,t≤Nk,∀k∈K,t∈0,tendt.The selection of arcs determines the arrival nodesuk and departure nodes wk, thereby determining every element in the optimization model. Thus, the model is transformed into an MIP. The next section discusses how its solution can be found.


## 3.1. Problem Formulation

Table2 summarizes the parameters and variables used in this paper. Some are also explained further in the text. The model’s simplifying assumptions are listed below and also explained in the problem description.Table 2
Parameters and variables.
Name Description Units Purpose I Set of stations - Set K Set of trains - Set U i Set of time stamps in which trains can arrive at stationi - Set W i Set of time stamps in which trains departure from stationi - Set  θ Shunting yard of the metro line - Symbol i ,  j Station index - Symbol I Total number of stations - Symbol k Train index - Symbol K Total number of trains within a given period/Frequency - Symbol t Time index - Symbol u k Arrival node of traink - Symbol w k Departure node of traink - Symbol π ( u k ) Time stamp associated with nodeuk (or wk) - Symbol Δ ( u k ,  w k ′ ) Time duration between two nodesuk and wk′ (k=k′ is allowed) - Symbol R k Path of traink - Symbol T Time duration - Symbol a The largest integer that is smaller than real numbera - Symbol a The smallest integer that is larger than real numbera - Symbol b ⃑ A vector (represented by any bold alphabet) - Symbol M A matrix (represented by any double-struck alphabet) - Symbol M a : b ,  x : y A sub-matrix whose elements correspond to those in matrixM, from rows a to b, and from columns from x to y - Symbol M Sum of all (numerical) elements in matrixM - Symbol  c e Average energy cost ¥/Kwh Parameter c o Average train operating cost ¥ /(train⋅km) Parameter c p t Average passenger travel time cost ¥/(passenger⋅hour) Parameter c p w Average passenger wait time cost ¥/(passenger⋅hour) Parameter c t r a Average capital cost of one train ¥/(train⋅hour) Parameter C t o t a l Total cost of the whole system ¥ Parameter C t o t a l e Total cost of energy consumption ¥ Parameter C t o t a l p Total cost of passenger time ¥ Parameter C t o t a l t r a Total cost of train (rolling stock) ¥ Parameter d i ¯ Maximum allowable dwell time at stationi Second Parameter d i _ Minimum allowable dwell time at stationi Second Parameter e i Energy consumption per unit mass for one train in segmenti Kwh/Kg Parameter f t Tractive force per unit mass m / s 2 Parameter f r Running resistance per unit mass m / s 2 Parameter F s Fleet size (the number of trains) Trains Parameter h m i n Minimum headway Second Parameter M k Net weight of a train Kg/train Parameter M p Average weight of a passenger Kg/passenger Parameter n d Number of doors per train - Parameter N k Capacity of traink Passengers Parameter P k ,  i a Number of passengers alighting traink at station i Passengers Parameter P k ,  i b Number of passengers boarding traink at station i Passengers Parameter P k ,  ( i ,  j ) b Number of passengers boarding traink at station i with destination j Passengers Parameter P i ( t 1 ,  t 2 ) Number of passengers arriving at stationi during [t1,t2] Passengers Parameter P ( i ,  j ) t 1 ,  t 2 Number of passengers arriving at stationi with destination j during [t1,t2] Passengers Parameter P i ,  t w Number of passengers waiting at stationi at time t Passengers Parameter P k ,  t v Number of passengers in traink at time t Passengers Parameter s l ( i ) Length of segmenti (between station i and i+1) Meter Parameter s t m a x ( i ) Maximum segment travel time for segmenti Second Parameter s t m i n ( i ) Minimum segment travel time for segmenti Second Parameter t k ,  i 0 Arrival time of traink at station i Second Parameter t k ,  i a Alighting end time for traink at station i Second Parameter t k ,  i b Boarding end time for traink at station i/Departure time Second Parameter t e n d p End time of passenger demand for a period Second Parameter t e n d t End time of a complete timetable for a period Second Parameter T r o u n d k Round trip time (total travel time from the departure from the shunting yard to the arrival at the shunting yard) of traink Second Parameter T p t Total passenger travel time Second Parameter T p w Total passenger wait time Second Parameter T t u r n Turnaround time Second Parameter v k ( t ) Speed of traink at time t m/s Parameter x k ( t ) Position of traink at time t Meter Parameter τ i ,  j ( t ) Passenger arrival rate at stationi with destination j at time t Passengers/second Parameter λ i ( t ) Passenger arrival rate at stationi at time t Passengers/second Parameter ζ a Passenger alighting rate Passengers/second Parameter ζ b Passenger boarding rate Passengers/second Parameter χ x Parameter for energy approximation - Parameter χ y Parameter for energy approximation - Parameter  s t ( i ) Segment travel time of segmenti Second Decision t k ,  1 0 Arrival time of traink at the first station Second DecisionAssumption 1.
Within one period, trains travel through each segment with the same speed profile; that is, the travel times for different trains in a segment are equal.Assumption 2.
We assume that the last train arrives at the first station at the end time of the period, to ensure that no late-arriving passengers will be left behind.Assumption 3.
Passenger arrival rateλi(t) should not exceed total passenger boarding rate ζb·nd. If λit>ζb·nd, the queue of waiting passengers increases even if a train is dwelling at the station. In practice, passenger entrance flow control∗ is taken to avoid such a case.R e m a r k ∗. At most busy metro stations in China, passenger entrance flow control is adopted during peak hours to ensure operation safety, and thus the passenger arrival rate is a restricted number, which is expected not to exceed passenger boarding rate. Assumption 3 is introduced for consistency with practical experience.Assumption 4.
All passengers in traink who have reached their destination are assumed to exit the train at its arrival time of station i, that is, tk,i0. All waiting passengers at station i able to board train k are assumed to board the train at its departure time, that is, tk,ib. In other words, the number of in-vehicle passengers Pk,tv changes only at the arrival/departure times of train k at stations.This paper focuses on a bidirectional metro line with a single track, as shown in Figure1. Properties of such a metro line are described in [48]. The timetabling problem is studied within a given period. The time-varying passenger demand is given as input, and a series of trains (whose frequency is to be optimized) are scheduled to transport the passengers.Figure 1
A bidirectional metro line with 2I stations.The arrival times at the first station are decision variables that must be optimized. The departure times at the first station are determined by the dwell times. As mentioned earlier, this paper obtains the dwell times based on trains’ conditions when they arrive at stations. In this regard, optimizing arrival times at the first station is equivalent to optimizing departure times at the first station, as in [2, 5]. For the second station, the arrival times are determined by the segment travel time in the first segment. Practically, for the metro operators, uniform segment travel times are preferred within one period, because they not only make the system more stable and easier to manage but also simplify synchronized passenger transfers with connecting transit lines. Assumption 1 is introduced to consider this preference. Thus, within one period, the segment travel time in each segment needs to be optimized only once (for all trains within this period). The speed profile optimization is outside the scope of this paper. The departure times at the second station are also determined by the dwell time determining process.Analogously, for subsequent train arrivals and departures, the arrival times are determined by the previous segment travel time, and the departure times are determined by the dwell times. Thus, within a period, the arrival times at the first station and the segment travel times determine the train timetable.To represent the TTP, a directed multigraph called discrete time-space graph (DTSG) is proposed [6]. As shown in Figure 2, θ represents the shunting yard of the metro line and Ui  (Wi) corresponds to the set of time stamps in which trains can arrive at (depart from) station i. The nodes in Ui  (Wi) are called arrival (departure) nodes. Unlike the original graph that is designed for a railway network [6, 12], this paper considers a single-track metro line, and thus the constraints are different and implicitly imposed in the definition of the structure of the DTSG.Figure 2
Representation of discrete time-space graph.The continuous time is discretized into evenly spaced time intervals (1 second in this paper), and the period extends from 0 totendt. Note that the values of tendt and tendp are different. tendp is the last time stamp of the passenger demand. The period mentioned above refers to the one extending from 0 to tendp. Assumption 2 imposes that the last train departs from the first station at time tendp. However, if the DTSG only considers the period from 0 to tendp, the arrival/departures of the last train at subsequent stations are not considered, and the passenger-oriented objectives are hard to formulate. Thus, the DTSG considers the period from 0 to tendt. The relation between tendt and tendp is set as:(1)tendt=tendp+∑i=12I-1sti+∑i=12Idi¯,where st(i) is the segment travel time of segment i and di¯ is the maximum dwell time of station i. This setting ensures that even the slowest train can be completely included in the DTSG.The construction of the DTSG is noncyclic. The objective functions of passenger time, energy, and train cost are associated with the chosen arcs and are described below.


## 3.2. Passenger Time

### 3.2.1. Passenger Arriving, Alighting, and Boarding

(1) Time-Varying Demand Information. The passenger demand at each station can be expressed by τi,j(t), which is defined as(2)τi,jt≥0,If1≤i<j≤IorI≤i<j≤2I;τi,jt=0,Otherwise;∀i,j∈I.Passenger arrival rates at stations areλit=∑j=i+12Iτi,j(t), ∀i,j∈I. In this paper, demands are regarded as deterministic values, which are given by τi,j(t). Thus, the number of arriving passengers at station i with destination j during time interval [t1,t2] is(3)Pi,jt1,t2=∑t=t1t2τi,jt,∀i,j∈I.The number of arriving passengers at stationi within time interval [t1,t2] is the sum of arriving passengers with different destinations:(4)Pit1,t2=∑j=i+12IPi,jt1,t2=∑j=i+12I∑t=t1t2τi,jt=∑t=t1t2λit,∀i,j∈I.(2) Determination of Dwell Arcs. Letting nd be the number of doors per train, the average number of alighting (boarding) passengers per door when train k arrives at (departs from) station i is Pk,ia/nd  (Pk,ib/nd). In addition, we denote the average number of boarding passengers per second (boarding rate) as ζb and alighting rate as ζa, which vary with passenger density in trains. For simplicity, we consider that boarding passengers wait until all alighting passengers get off the train; that is, the boarding process starts at time tk,ia.The dwell time is the sum of alighting and boarding times. The calculation of alighting time is solely based on the number of alighting passengers:(5)tk,ia=tk,i0+Pk,iand·ζa,∀k∈K,i∈I.The calculation of boarding time is slightly more complex. Figure3 shows the relations of train remaining capacity, waiting passengers, boarding rate, and max/min dwell times. The horizontal axis represents time, and the vertical axis represents the number of passengers. These four factors determine the boarding time for train k at station i, and in Figure 3 they are represented by (1) remaining capacity of train k at station i: Nk-Pk,tk,iav, a horizontal line; (2) number of waiting passengers at station i at time t, a curve Pi,tw; (3) maximum allowable number of boarding passengers, a diagonal function fkt with the slope of nd·ζb; (4) minimum and maximum dwell times at station i, that is, tk,i0+di_ and tk,i0+di¯. In Figure 3, the diagonal P=fkt intersects the horizontal line P=Nk-Pk,tk,iav at time t1k,i and intersects the curve P=Pi,tw at time t2k,i. The curve P=Pi,tw intersects the horizontal line P=Nk-Pk,tk,iav at time t3k,i. Note that the time is discrete, and Assumption 3 ensures that the slope of P=fkt always exceeds the slope of P=Pi,tw. Thus, we have(6)t1k,i=max⁡t∣t-tk,ia·nd·ζb≤Nk-Pk,tk,iav,∀k∈K,i∈I,(7)t2k,i=maxt∣t-tk,ia·nd·ζb≤Pi,tw,∀k∈K,i∈I,(8)t3k,i=max⁡t∣Pi,tw≤Nk-Pk,tk,iav,∀k∈K,i∈I.Figure 3
Relations of remaining capacity, waiting passengers, boarding rate, and max/min dwell times.Practically, in an oversaturated scenario, a train departs from a station as soon as the number of in-train passengers reaches its capacity. Based on this strategy, (6) shows that the boarding passengers who arrive after t1k,i exceed train capacity, and they will be left behind; (7) shows that before t2k,i there are passengers waiting at station, and after t2k,i all passengers board the train; (8) points that the train is full at t3k,i, and passengers arriving after t3k,i should wait for the next train.The departure time of traink at station i  (tk,ib) is determined based on the relations among t1k,i, t2k,i, t3k,i, tk,i0+di_, and tk,i0+di¯. Based on Assumption 3, we consider two cases: (1)  t1k,i>t2k,i and (2)  t1k,i≤t2k,i. Note that di_ and di¯ should be considered in both cases. In case (1), we can further derive that t3k,i>t1k,i>t2k,i (which is the solid curve in Figure 3), where the number of boarding passengers per second is determined by passenger arrival rate. Thus, the departure time is determined as tk,ib=min{tk,i0+di¯,max⁡t2k,i,tk,i0+di_,t3k,i}. In case (2), we can further derive that t3k,i≤t1k,i≤t2k,i (which is the dashed curve in Figure 3), where the number of boarding passengers per second is determined by the actual passenger boarding rate. Thus, the departure time is determined as tk,ib=min{tk,i0+di¯,t1k,i}. The departure time of train k at station i is expressed as(9)tk,ib=min⁡tk,i0+di¯,max⁡t2k,i,tk,i0+di_,t3k,i,t1k,i>t2k,imin⁡tk,i0+di¯,t1k,i,t1k,i≤t2k,i.Finally, we consider the relations between alighting/boarding rates and congestion at stations. Here the congestion is described by the number of in-vehicle passengers. Practically, the greater the passenger density in a train is, the slower the passengers can get off/on the train. Here, we design a monotonically decreasing functionρ(Pk,tv)≤1 with respect to the number of in-vehicle passengers, to approximate actual alighting/boarding rates according to the passenger density in trains. For simplicity, in this paper, the value of ρ(Pk,tv) during the dwell process is fixed and is determined by the number of in-vehicle passengers at its arrival time. That is, ζa and ζb remain constant during the dwell process:(10)ζa=ρPk,tk,i0v·ζaidealζb=ρPk,tk,i0v·ζbideal;∀k∈K,i∈I,where ζaideal and ζbideal are the ideal alighting and boarding rates in the uncongested condition.(3) Passenger Alighting and Boarding. Once the dwell arcs are determined, the alighting and boarding times for train k at station i are determined. Based on the results of dwell arcs, three cases should be considered: (1)  t1k,i<t2k,i; (2)  t1k,i≥t2k,i and t2k,i≤tk,i0+di_; (3)  t1k,i≥t2k,i and t2k,i>tk,i0+di_. In case (1), waiting passengers board the train at rate ζb. In case (2), actual passenger boarding rate is equal to the passenger arrival rate. In case (3), the actual boarding rate is ζb. Thus, the number of boarding passengers for train k at station i is(11)Pk,ib=Pi,tk,ibw;ift1k,i≥t2k,i,t2k,i≤tk,i0+di_;tk,ib-tk,ia·ζb·nd;otherwise;∀k∈K,i∈I.Wang et al. [3] point out that, in a busy metro line, passengers with different destinations are well mixed at each station. Thus, the ratio between the numbers of boarding passengers with destination j and of total boarding passengers is equivalent to the ratio between τi,j(t) and λi(t):(12)Pk,i,jb=∑t=tk-1,ibtk,ibτi,jt∑t=tk-1,ibtk,ibλit·Pk,ib;∀k∈K,j∈I\1,i=1,2,…,j-1.The number of alighting passengers from traink at station j consists of all previously boarding passengers in train k with destination j:(13)Pk,ja=∑i=1j-1Pk,i,jb;∀k∈K,j=2,…,I,I+2,…,2I.There are no alighting passengers at the first station in each direction (station 1 andI+1).


### 3.2.2. Passenger Wait and Travel Time

(1) Passenger Wait Time. Total passenger wait time at a station can be calculated by summing up the products of the cumulative number of passengers and their time. Let klast be the index of last train that has departed from station i before time t, that is, tklast,ib≤t<tklast+1,ib. The cumulative number of passengers at station i at time t is the difference between cumulative numbers of arrival passengers and boarding passengers within the period [0,t]:(14)Pi,tw=Pi0,t-∑k=1klastPk,ib;∀i∈I.Since the continuous time is discretized into 1 second intervals, the passenger wait time between timest and t+1 is Pi,tw·t+1-t=Pi,tw. Thus, the total passenger wait time is(15)Tpw=∑i∈I∑t=0tendtPi,tw.(2) Passenger Travel Time. The number of in-train passengers in train k at time t is the difference between cumulative numbers of boarding passengers and alighting passengers. Let ilast be the index of the last station that train k has visited before time t. We consider two cases: (1) train k is dwelling at station ilast, that is, tk,ilast0≤t<tk,ilastb; (2) train k is in the segment between station ilast and ilast+1, that is, tk,ilastb≤t<tk,ilast+10. Based on Assumption 4, the boarding passengers at station ilast are not considered in case (1) but are included in case (2). Thus, the number of in-vehicle passengers in train k at time t is(16)Pk,tv=∑i=1ilast-1Pk,ib-∑i=1ilastPk,ia;tk,ilast0≤t<tk,ilastb;∑i=1ilastPk,ib-Pk,ia;tk,ilastb≤t<tk,ilast+10;∀k∈K.Analogously, the total passenger travel time is(17)Tpt=∑k∈K∑t=0tendtPk,tv.Finally, the total cost of passenger time is formulated as(18)Ctotalp=cpw·Tpw3600+cpt·Tpt3600.The cost unit is yuan(¥).


## 3.2.1. Passenger Arriving, Alighting, and Boarding

(1) Time-Varying Demand Information. The passenger demand at each station can be expressed by τi,j(t), which is defined as(2)τi,jt≥0,If1≤i<j≤IorI≤i<j≤2I;τi,jt=0,Otherwise;∀i,j∈I.Passenger arrival rates at stations areλit=∑j=i+12Iτi,j(t), ∀i,j∈I. In this paper, demands are regarded as deterministic values, which are given by τi,j(t). Thus, the number of arriving passengers at station i with destination j during time interval [t1,t2] is(3)Pi,jt1,t2=∑t=t1t2τi,jt,∀i,j∈I.The number of arriving passengers at stationi within time interval [t1,t2] is the sum of arriving passengers with different destinations:(4)Pit1,t2=∑j=i+12IPi,jt1,t2=∑j=i+12I∑t=t1t2τi,jt=∑t=t1t2λit,∀i,j∈I.(2) Determination of Dwell Arcs. Letting nd be the number of doors per train, the average number of alighting (boarding) passengers per door when train k arrives at (departs from) station i is Pk,ia/nd  (Pk,ib/nd). In addition, we denote the average number of boarding passengers per second (boarding rate) as ζb and alighting rate as ζa, which vary with passenger density in trains. For simplicity, we consider that boarding passengers wait until all alighting passengers get off the train; that is, the boarding process starts at time tk,ia.The dwell time is the sum of alighting and boarding times. The calculation of alighting time is solely based on the number of alighting passengers:(5)tk,ia=tk,i0+Pk,iand·ζa,∀k∈K,i∈I.The calculation of boarding time is slightly more complex. Figure3 shows the relations of train remaining capacity, waiting passengers, boarding rate, and max/min dwell times. The horizontal axis represents time, and the vertical axis represents the number of passengers. These four factors determine the boarding time for train k at station i, and in Figure 3 they are represented by (1) remaining capacity of train k at station i: Nk-Pk,tk,iav, a horizontal line; (2) number of waiting passengers at station i at time t, a curve Pi,tw; (3) maximum allowable number of boarding passengers, a diagonal function fkt with the slope of nd·ζb; (4) minimum and maximum dwell times at station i, that is, tk,i0+di_ and tk,i0+di¯. In Figure 3, the diagonal P=fkt intersects the horizontal line P=Nk-Pk,tk,iav at time t1k,i and intersects the curve P=Pi,tw at time t2k,i. The curve P=Pi,tw intersects the horizontal line P=Nk-Pk,tk,iav at time t3k,i. Note that the time is discrete, and Assumption 3 ensures that the slope of P=fkt always exceeds the slope of P=Pi,tw. Thus, we have(6)t1k,i=max⁡t∣t-tk,ia·nd·ζb≤Nk-Pk,tk,iav,∀k∈K,i∈I,(7)t2k,i=maxt∣t-tk,ia·nd·ζb≤Pi,tw,∀k∈K,i∈I,(8)t3k,i=max⁡t∣Pi,tw≤Nk-Pk,tk,iav,∀k∈K,i∈I.Figure 3
Relations of remaining capacity, waiting passengers, boarding rate, and max/min dwell times.Practically, in an oversaturated scenario, a train departs from a station as soon as the number of in-train passengers reaches its capacity. Based on this strategy, (6) shows that the boarding passengers who arrive after t1k,i exceed train capacity, and they will be left behind; (7) shows that before t2k,i there are passengers waiting at station, and after t2k,i all passengers board the train; (8) points that the train is full at t3k,i, and passengers arriving after t3k,i should wait for the next train.The departure time of traink at station i  (tk,ib) is determined based on the relations among t1k,i, t2k,i, t3k,i, tk,i0+di_, and tk,i0+di¯. Based on Assumption 3, we consider two cases: (1)  t1k,i>t2k,i and (2)  t1k,i≤t2k,i. Note that di_ and di¯ should be considered in both cases. In case (1), we can further derive that t3k,i>t1k,i>t2k,i (which is the solid curve in Figure 3), where the number of boarding passengers per second is determined by passenger arrival rate. Thus, the departure time is determined as tk,ib=min{tk,i0+di¯,max⁡t2k,i,tk,i0+di_,t3k,i}. In case (2), we can further derive that t3k,i≤t1k,i≤t2k,i (which is the dashed curve in Figure 3), where the number of boarding passengers per second is determined by the actual passenger boarding rate. Thus, the departure time is determined as tk,ib=min{tk,i0+di¯,t1k,i}. The departure time of train k at station i is expressed as(9)tk,ib=min⁡tk,i0+di¯,max⁡t2k,i,tk,i0+di_,t3k,i,t1k,i>t2k,imin⁡tk,i0+di¯,t1k,i,t1k,i≤t2k,i.Finally, we consider the relations between alighting/boarding rates and congestion at stations. Here the congestion is described by the number of in-vehicle passengers. Practically, the greater the passenger density in a train is, the slower the passengers can get off/on the train. Here, we design a monotonically decreasing functionρ(Pk,tv)≤1 with respect to the number of in-vehicle passengers, to approximate actual alighting/boarding rates according to the passenger density in trains. For simplicity, in this paper, the value of ρ(Pk,tv) during the dwell process is fixed and is determined by the number of in-vehicle passengers at its arrival time. That is, ζa and ζb remain constant during the dwell process:(10)ζa=ρPk,tk,i0v·ζaidealζb=ρPk,tk,i0v·ζbideal;∀k∈K,i∈I,where ζaideal and ζbideal are the ideal alighting and boarding rates in the uncongested condition.(3) Passenger Alighting and Boarding. Once the dwell arcs are determined, the alighting and boarding times for train k at station i are determined. Based on the results of dwell arcs, three cases should be considered: (1)  t1k,i<t2k,i; (2)  t1k,i≥t2k,i and t2k,i≤tk,i0+di_; (3)  t1k,i≥t2k,i and t2k,i>tk,i0+di_. In case (1), waiting passengers board the train at rate ζb. In case (2), actual passenger boarding rate is equal to the passenger arrival rate. In case (3), the actual boarding rate is ζb. Thus, the number of boarding passengers for train k at station i is(11)Pk,ib=Pi,tk,ibw;ift1k,i≥t2k,i,t2k,i≤tk,i0+di_;tk,ib-tk,ia·ζb·nd;otherwise;∀k∈K,i∈I.Wang et al. [3] point out that, in a busy metro line, passengers with different destinations are well mixed at each station. Thus, the ratio between the numbers of boarding passengers with destination j and of total boarding passengers is equivalent to the ratio between τi,j(t) and λi(t):(12)Pk,i,jb=∑t=tk-1,ibtk,ibτi,jt∑t=tk-1,ibtk,ibλit·Pk,ib;∀k∈K,j∈I\1,i=1,2,…,j-1.The number of alighting passengers from traink at station j consists of all previously boarding passengers in train k with destination j:(13)Pk,ja=∑i=1j-1Pk,i,jb;∀k∈K,j=2,…,I,I+2,…,2I.There are no alighting passengers at the first station in each direction (station 1 andI+1).


## 3.2.2. Passenger Wait and Travel Time

(1) Passenger Wait Time. Total passenger wait time at a station can be calculated by summing up the products of the cumulative number of passengers and their time. Let klast be the index of last train that has departed from station i before time t, that is, tklast,ib≤t<tklast+1,ib. The cumulative number of passengers at station i at time t is the difference between cumulative numbers of arrival passengers and boarding passengers within the period [0,t]:(14)Pi,tw=Pi0,t-∑k=1klastPk,ib;∀i∈I.Since the continuous time is discretized into 1 second intervals, the passenger wait time between timest and t+1 is Pi,tw·t+1-t=Pi,tw. Thus, the total passenger wait time is(15)Tpw=∑i∈I∑t=0tendtPi,tw.(2) Passenger Travel Time. The number of in-train passengers in train k at time t is the difference between cumulative numbers of boarding passengers and alighting passengers. Let ilast be the index of the last station that train k has visited before time t. We consider two cases: (1) train k is dwelling at station ilast, that is, tk,ilast0≤t<tk,ilastb; (2) train k is in the segment between station ilast and ilast+1, that is, tk,ilastb≤t<tk,ilast+10. Based on Assumption 4, the boarding passengers at station ilast are not considered in case (1) but are included in case (2). Thus, the number of in-vehicle passengers in train k at time t is(16)Pk,tv=∑i=1ilast-1Pk,ib-∑i=1ilastPk,ia;tk,ilast0≤t<tk,ilastb;∑i=1ilastPk,ib-Pk,ia;tk,ilastb≤t<tk,ilast+10;∀k∈K.Analogously, the total passenger travel time is(17)Tpt=∑k∈K∑t=0tendtPk,tv.Finally, the total cost of passenger time is formulated as(18)Ctotalp=cpw·Tpw3600+cpt·Tpt3600.The cost unit is yuan(¥).


## 3.3. Energy

In general, for each segmenti, with given segment travel time st(i) and segment length sl(i), the speed profile determines energy use [49]. Thus, given a speed profile generating method, the energy per unit mass of different trains in segment i is the same. We denote the energy per unit mass for one train in segment i as ei. Its calculation usually considers the following train dynamic equations [48]:(19)dvktdt=ft-fr,dxktdt=vkt;∀t∈tk,ib,tk,i+10,i,i+1∈I,fr=α·vk2t+β·vkt+γ,vktk,i0=vktk,ib=0.The first two equations are train motion equations, wherevkt and xk(t), respectively, represent traveling speed and distance. The third one is the Davis resistance function where Davis parameters are α, β, and γ. The last one requires that every train stops at each station. An optimal train control strategy can be obtained based on (19), in which energy per unit mass is(20)ei=∫tk,ibtk,i+10ft·vktdt;∀i,i+1∈I,k∈K.However, in (20), the relation between energy and decision variable st(i) is unclear. It is also difficult to integrate the optimal train control with other objective functions [18]. Since the speed profile optimization is not considered, for simplicity, we use the method in [50] to generate practical speed profiles. Based on these speed profiles, we then adopt a linear piecewise curve to approximate the energy per unit mass of train k in segment i [48]:(21)ei=χx·sti+χy·sli;∀i∈I.We introduce here an estimation procedure for the values ofχx and χy, which is based on mean squared error (MSE). First, we set several feasible values (e.g., z=1,2,…,100) of slz and st(z). Then, we obtain optimal speed profiles based on (19) and (20), thereby obtaining the optimal energy with respect to slz and st(z). Finally, we define the optimal energy per unit mass as ez, and the estimated energy per unit mass obtained with (21) as ez~. The MSE is formulated as(22)MSEχx,χy=1100∑z=1100ez~-ez2.The optimal values ofχx and χy can be obtained when MSEχx,χy reaches its minimum and are found as follows:(23)∂MSEχx,χy∂χx=0,∂MSEχx,χy∂χy=0.Withχx∗ and χy∗ obtained, the total energy can be calculated as(24)Etotal=∑i=12I-1∑k=1Kχx∗·sti+χy∗·sli·Mk+Mp·Pk,tk,ibv.Thus, the total cost of energy consumption is(25)Ctotale=ce·Etotal=ce·∑i=12I-1∑k=1Kχx∗·sti+χy∗·sli·Mk+Mp·Pk,tk,ibv.Note that when all tracks are straight and level,χx and χy suffice for estimating energy. However, if grades and curves are considered, we may have ei=χxi·st(i)+χyi·sli,∀i∈I, where χxi and χyi are specified parameters to describe the characteristics in segment i. The parameter estimation procedure is the same.


## 3.4. Train Cost

In this paper, the train cost is estimated with empirical equations with respect to the interest rate, train price, and average round trip time, where long-term requirements (e.g., the consideration of the highest peak during a year and demand evolution over the planning horizon) are not considered.Train cost consists of train capital cost and operating cost (excluding energy). Train capital cost is related to the price and economic life of a train, as well as interest rate. LetCannual be the annual cost per train, Ctrain be the current price per train, r be the interest rate, and Tls be the lifespan of a train. The relations among these four factors can be expressed as(26)Cannual=Ctrain·r·1+rTls1+rTls-1.The capital cost per train hourctra can be then calculated by dividing the annual cost Cannual by the product of average working hours per day (Dwork) and working days per year (Ywork)(27)ctra=CannualDwork·Ywork.In addition, a train can be reassigned after it finishes a round trip and returns to the shunting yard, and the round trip timeTroundk can be obtained as(28)Troundk=tk,2Ib-tk,10+Tturn,∀k∈K.Thus, the actual period for train running is∑k=1KTroundk  (∑k=1KTroundk≤tendt). The estimated fleet size Fs within period [0,∑k=1KTroundk] depends on the round trip time Troundk and the departure headway from the shunting yard (which is equivalent to the arrival time of train k at the first station tk,10):(29)Fs=∑k=1KTroundk/K∑k=2Ktk,10-tk-1,10/K-1,where the term (∑k=1KTroundk)/K represents the average round trip time and the term ∑k=2K(tk,10-tk-1,10)/(K-1) represents the average headway. The total cost of trains is formulated as(30)Ctotaltra=ctra·Fs·∑k=1KTroundk3600+co·K·∑i=12I-1sli1000.The first termctra·Fs·(∑k=1KTroundk/3600) represents the capital cost of Fs trains for the actual running period [0,∑k=1KTroundk]. The second term co·K·∑i=12I-1sli represents the operating cost with frequency K (trains) on a line with the length of the sum of 2I-1 segments (Km). Both decision variables (segment travel times and arrival times at the first station) affect the train cost.


## 3.5. Mixed Integer Programming Formulation

The optimization model aims to minimize the weighted-sum of different cost functions, so that the conflicting objectives are generalized in terms of monetary cost. In this paper, the weights of different functions are reflected by their average costs, that is,ce, co, cpt, cpw, and ctra. For example, an operator-oriented timetable (neglecting passenger satisfaction) can be obtained by setting cpt=0 and cpw=0; a passenger-oriented timetable can be obtained by setting ce=co=ctra=0. The weights for integrated objectives can be set similarly, with larger cost coefficients for larger weights. The model is formulated as follows:Minimize(31)Ctotal=Ctotalp+Ctotale+Ctotaltra.Subject to(32)tK,10=tendp,(33)stI=Tturn,(34)stmini≤sti≤stmaxi;∀i∈I,(35)tk,ib-tk-1,ib≥hmin;∀k=2,…,K,i∈I,(36)Pk,tv≤Nk;∀k∈K,t∈0,tendt.This model is a nonsmooth, nonconvex programming model, where the nonsmoothness is caused by the limited train capacity, and the nonconvexity is caused by the determination of dwell times [3]. Constraint (32) is equivalent to Assumption 2. Constraint (33) considers the turnaround from station I to I+1 as traveling in segment I, with fixed travel time Tturn. Constraint (34) ensures that the segment travel time should be optimized within a reasonable range [stmini,stmaxi]. Constraint (35) imposes the minimum headway constraint for both departures and arrivals (based on Assumption 1). Constraint (36) requires that the number of in-train passengers must not exceed train capacity.In Caprara et al. [6], the time-space graph is used to derive an integer linear programming model, thereby allowing a considerably faster solution. Analogously, this paper transforms the nonsmooth, nonconvex programming model into a mixed integer programming (MIP) model based on the DTSG. First, we construct two (2I×tendt) matrixes U and W, respectively, representing the selection of arrival and departure nodes in the DTSG:(37)Ui,t=1;ift=πuk,uk∈Ui,k∈K;0;otherwise;Wi,t=1;ift=πwk,wk∈Wi,k∈K;0;otherwise;∀i∈I,t∈0,tendt,where Ui,t=1  (Wi,t=1) if a train arrives at (departs from) station i at time t. As mentioned earlier, the value of U1,t is determined by πuk=tk,10, for ∀uk∈U1,  k∈K. The value of Wi,t is determined by Ui,t and the dwell arcs, that is, Δuk,wk=tk,ib-tk,i0, for ∀uk∈Ui,  wk∈Wi,  i∈I,  k∈K. The value of Ui+1,t is determined by Wi,t and the travel time in segment i, that is, Δwk,uk=st(i), for ∀wk∈Wi,  uk∈Ui+1,  and i∈I\{2I},  k∈K. Thus, the decision variables determine U and W.Based on the structures ofU and W, as well as Assumption 4, we introduce two (2I×tendt) matrixes PA and PB, respectively, representing passenger alighting and boarding information (number and time):(38)PAi,t=Pk,ia;t=πuk,uk∈Ui,k∈K;0;otherwise;PBi,t=Pk,ib;t=πwk,wk∈Wi,k∈K;0;otherwise;∀i∈I,t∈0,tendt.The cumulative numbers of waiting passengersQ  (2I×tendt) and in-vehicle passengers V  (K×tendt) can be constructed based on PA and PB:(39)Qi,t=Qi,t-PBi,1:t;∀i∈I,t∈0,tendt,(40)Vk,t=Vk,t-1+Pk,ib;ift=πwk,wk∈Wi,i∈I;Vk,t-1-Pk,ia;ift=πuk,uk∈Ui,i∈I;Vk,t-1;otherwise;∀k∈K,t∈1,tendt.In (39), Qi,t=Qi,0+∑q=0tλi(t) represents the total waiting passengers at station i, where Qi,0 is the number of waiting passengers at the start of this period. According to the definition of cumulative passengers, we have Qi,t=Pi,tw. Thus, the total passenger wait time can be expressed by Q:Tpw=∑i∈I∑t=0tendtPi,tw=∑i∈I∑t=0tendtQi,t=Q  (since Qi,t≥0). In (40), the initial number of in-vehicle passengers is zero, that is, Vk,0=0. The number of in-vehicle passengers changes only at arrival/departure times, which is consistent with Assumption 4. Analogously, the total passenger travel time can be expressed as Tpt=V. Thus, the total cost for passenger time can be rewritten as(41)Ctotalp=13600cpw·Q+cpt·V.For energy cost, the energy per unit mass in segmenti can be rewritten as(42)ei=χx·Δwk,uk+χy·sli,wk∈Wi,uk∈Ui+1;∀i∈I\2I.The total energy can be rewritten as(43)Etotal=∑i=12I-1∑k=1K∑t=0tendtϕt·χx·Δwk,uk+χy·sli·Mk+Mp·Vk,t,wk∈Wi,uk∈Ui+1,where ϕt=1 if t=πwk; otherwise ϕt=0. The total energy cost can be rewritten correspondingly.For train cost, the round trip timeTroundk can be rewritten as(44)Troundk=Δuk,wk+Tturn,uk∈U1,wk∈W2I;∀k∈K.The departure headway from the shunting yard can be rewritten as(45)tk,10-tk-1,10=Δuk-1,uk,uk-1,uk∈U1;∀k∈K\1.The total cost of trains can be rewritten correspondingly.Finally, the last departing train constraint (32) is equivalent to (46)πuK=tendt,uK∈U1.The turnaround constraint (33) is equivalent to(47)Δwk,uk=Tturn,wk∈WI,uk∈UI+1;∀k∈K.The segment travel time window constraint (34) is equivalent to(48)stmini≤Δwk,uk≤stmaxi,wk∈Wi-1,uk∈Ui;∀k∈K,i∈I\1.The minimum headway constraint (35) is equivalent to(49)Δwk-1,wk≥hmin,wk-1,wk∈Wi;k∈K\1,i∈I.The train capacity constraint (36) is equivalent to(50)Vk,t≤Nk,∀k∈K,t∈0,tendt.The selection of arcs determines the arrival nodesuk and departure nodes wk, thereby determining every element in the optimization model. Thus, the model is transformed into an MIP. The next section discusses how its solution can be found.


## 4. Solution Approach

In this section, we decompose the TTP into two subproblems. The first subproblem optimizes the segment travel times (OST), aiming at minimizing the costs for energy and passenger segment travel times. A branch-and-bound algorithm and a frequency determining algorithm are introduced to solve the OST. The second subproblem optimizes departure headways from the shunting yard (OH), aiming to minimize passenger wait time, dwell time and train cost. A novel rolling optimization algorithm is designed to solve the OH.


### 4.1. Decomposition

The rationale of the decomposition is explained as follows. The objectives in the optimization model include passenger wait time (at their original stations), passenger travel time (segment travel times and station dwell times), energy, and train cost. Based on (43), the total energy is related to the energy per unit mass and the number of passengers in each segment. The energy per unit mass is solely determined by segment travel times st(i). Since passengers board at their origins and alight at their destinations regardless of which train they take, the sum of in-vehicle passengers in one segment i is determined by passenger demand. In this regard, the total energy can be expressed as(51)Etotal=Mk·K·∑i=12I-1ei+Mp·∑i=12I-1∑j=i+12I∑r=1i∑t=0tendpτr,jt·ei.Therefore, the total energy cost depends only on segment travel timesst(i). The passenger segment travel times also depend only on segment travel times.On the other hand, passenger wait time, dwell time, and train cost depend on the train arrival and departure times at stations. Based on Assumption1, train arrival/departure times at stations are affected by the departure headways from the shunting yard. Thus, it is reasonable to use segment travel times as input and consider the costs for passenger wait time, dwell time, and trains as functions of the departure headways.In general, the OST subproblem is formulated as follows:Minimize(52)CtotalOST=∑i=12I-1ce·Mk·K·ei+ce·Mp·ei+cpt·sti·∑j=i+12I∑r=1i∑t=0tendpτr,jt.Subject to constraints (47) and (48).It should be noted that since segment travel times are given inputs, minimizing the cost of passenger dwell time is equivalent to minimizing the cost of total travel time. Since we already obtain the in-train passenger matrixV, there is no need to further derive dwell times at stations. Thus, the objectives in the OH subproblem include the total passenger time cost (both Q and V) instead of costs of passenger wait and dwell times. The OH subproblem is formulated as follows (with segment travel times as input):Minimize(53)CtotalOH=Ctotalp+Ctotaltra.Subject to constraints (46), (49) and (50).


### 4.2. Solution Approach for OST

The solution approach for OST consists of branch-and-bound (B&B) and frequency determining algorithms. On one hand, the B&B algorithm needs frequency as input. On the other hand, the frequency determining algorithm needs segment travel times obtained by the B&B to evaluate the performance. Thus, we iteratively evaluate the performance with different frequencies (each frequency is used as input for B&B to find a solution of segment travel times) and choose the one with the best performance. The segment travel times are then determined correspondingly.


#### 4.2.1. Branch-and-Bound (B&B) Algorithm

The OST subproblem is an integer linear programming problem with2I-1 decision variables (st(i), for i=1,2,…,2I-1). Here we use a generic B&B, which is widely used to solve linear programming problems [7]. Generally, a B&B needs an input model formulated as follows:(54)minc⃑·x⃑T,(55)s.t.A·x⃑T≤b⃑TAeq·x⃑T=beq⃑T,where c⃑ is the coefficient vector and x⃑T is the decision variable vector. Both c⃑ and x⃑T have 2I-1 elements. A and Aeq are coefficient matrixes, and b⃑T and beq⃑T are constant vectors representing the bounds. Here A is a (4I-2)×(2I-1) matrix, and Aeq is a 1×(2I-1) matrix. Correspondingly, b⃑T has 2I-1 elements, and beq⃑T has one element. We define a vector with 2I-1 elements as(56)ϵ⃑j=0,0,…,1,…,0,0,j=1,…,2I-1.The decision variable vector can be represented asx⃑i=st(i) for i=1,2,…,2I-1, where x⃑i is the i-th element in the vector x⃑. We then define ℸi=∑j=i+12I∑r=1i∑t=0tendpτr,jt for i=1,2,…,2I-1 as the total number of passengers that travel through segment i. Based on (52), c⃑i can be represented as(57)c⃑i=ce·Mk·K·χx·sli+ℸi·ce·Mp·χx·sli+cpt;i=1,2,…,2I-1.Constraints (55) are represented as follows. Based on constraint (47), we have Aeq=ϵ⃑I and beq⃑=Tturn. Based on constraint (48), we have Ai,:=-ϵ⃑i, for i=1,2,…,2I-1; Ai,:=ϵ⃑i, for i=2I,2I+1,…,4I-2, where Ai,: represents the i-th row of matrix A. Analogously, b⃑i=-stmini, for i=1,2,…,2I-1; b⃑i=stmaxi, for i=2I,2I+1,…,4I-2. Since we use the standard B&B algorithm [40], the detailed algorithm description is not shown here.


#### 4.2.2. Frequency Determining Algorithm

Here we consider cyclic timetables (train departure headways from the shunting yard are constant, e.g., 200 s) to discuss the frequency, because they have following properties:(1) The cyclic timetable with minimum cost (i.e., the best cyclic timetable) can be considered as an upper bound of the noncyclic timetable. (2) If passenger demand is stable, the best cyclic timetable is a good approximation to an optimal noncyclic one. (3) The cyclic timetables are easy to obtain.First, the range of frequency should be estimated. We define the peak passenger volume in one segment asqpeak, which is the maximum value among the numbers of passengers traveling through the segment, that is, qpeak=max⁡(ℸi∣i=1,2,…,2I-1). The minimum frequency Kmin is determined by the train capacity Nk and qpeak:Kmin=qpeak/Nk. The maximum frequency Kmax is determined by the minimum headway: Kmax=tendp/hmin+1.Second, for each frequencyK∈[Kmin,Kmax], we use K as input of the B&B, to obtain the segment travel times. Since the frequency for a cyclic timetable determines the departure headways from the shunting yard (h=tendp/(K-1), where h is the departure headway), the cyclic timetable that corresponds to frequency K can be obtained. Then, the total cost of the timetable can be calculated with (52).Finally, we obtain different cyclic timetables with different costs, each corresponding to a frequencyK∈[Kmin,Kmax]. We select the timetable with minimum cost as the best cyclic timetable, and its corresponding segment travel times as the OST solution within this period. The detailed algorithm is as shown in Algorithm 1.Algorithm 1:Frequency determining algorithm.
Step  1: Estimate the frequency range, i.e., [Kmin,Kmax]. Set the initial frequency of the cyclic timetable as K=Kmin, the best frequency found during the process asKbest=0, and the minimum total cost found during the process as TCmin=+∞. Go to Step 2. Step  2: For frequency K, obtain the segment travel times through the B&B algorithm, and obtain the uniform departure headways from the shunting yardh=tendp/(K-1). Go to Step 3. Step  3: Calculate total cost TCK corresponding to the frequency K. If TCK<TCmin, then set TCmin=TCK, and Kbest=K. Go to Step4. Step  4: If K=Kmax, then go to Step 5. Otherwise, set K=K+1, and return to Step 2. Step  5: Output Kbest and TCmin. The segment travel times within this period are the solution of the B&B algorithm that corresponds to frequencyKbest.


### 4.3. Solution Approach for OH

Here we design a rolling optimization (RO) algorithm that optimizes the arrival time at the first station for one train at a time. The main difference between the RO and the rolling horizon (RH) approach is in the RH, the bound conditions (i.e., the start and end times of each horizon) require some of the original variables to have fixed, known values [3], while in the RO, each decision variable may still be optimized. Thus, the rolling optimization is more flexible.


#### 4.3.1. The Rolling Optimization (RO) Algorithm

Based on (53), the objectives include the total costs of passenger time and trains. To obtain these costs, we need the segment travel times as input. Here, instead of considering the whole period, we consider a much shorter period tp≤hmin+∑i=12I-1st(i)+∑i=12Idi¯, e.g., 600 s, for solving the OST. For example, for each train k≥2, its previous train arrives at the first station at time tk-1,10 (the first train arrives at time 0), and we consider the period [tk-1,10,tk-1,10+tp] to solve the OST (the constraint of tp ensures that tk-1,10+tp≤tendt). In this way, the key ideas of “iterative optimization” (the frequencies, segment travel times, and headways are interdependent) and “rolling horizon” (the smaller period studied in each step) are adopted, and future passenger demand is considered.Then, we discuss the costs for one traink≥2 if its arrival time at the first station is given. Since each train is scheduled in chronological order, the waiting passenger matrix Q and the in-train passenger matrix V can be constructed chronologically. For each train k≥2, the departure time of its preceding train at the first station is tk-1,1b, and train k departs from the first station at time tk,1b. We define the passenger wait time with respect to train k as Tpw(k), which is formulated as(58)Tpwk=∑i=12I∑t=tk-1,1btk,1bQi,t;∀k∈K.Analogously, we define the passenger travel time for traink as Tpt(k), and we obtain(59)Tptk=∑t=tk,10tk,2IbVk,t;∀k∈K.The estimated fleet size of traink (Fsk) depends on its round trip time (Troundk) and the departure headway from the shunting yard (tk,10-tk-1,10):(60)Fsk=Troundktk,10-tk-1,10;∀k∈K.Based on (18), (30), (53), (58), (59), and (60), we define the total costs in the OH for train k as CtotalOH(k), which can be expressed as(61)CtotalOHk=cpw·Tpwk+cpt·Tptk+ctra·Fsk·Troundk3600+co·∑i=12I-1sli1000;∀k∈K.The value ofCtotalOHk is a function of the arrival time for train k at the first station. Since the arrival time for the previous train is given, we consider the value of CtotalOHk as a function of the headway h between train k-1 and train k (h=tk,10-tk-1,10), for consistency with other TTP studies.Next we discuss how the optimal value forh can be obtained. Based on (58), (59), and (60), both the passenger time and train cost depend on passenger demand. Since the passenger demand is time-varying, the relation between CtotalOHk and h cannot be described with a fixed function.The expected objective function of OH is shown in Figure4. Its value first decreases (the headway starts from hmin) to a minimum value and then increases monotonically. We specify that this function reaches its minimum value at headway h∗. Since the DTSG divided the continuous time into seconds, we iteratively calculate the objective function of OH with respect to h (starting from h=hmin, and h=h+1 in each following step), until the optimal headway h∗ is found, or the total cost increases as h increases (in case that h∗ cannot be reached). It should be noted that sometimes a headway is not feasible due to the model constraints (e.g., a short headway might violate the minimum headway constraint). In such cases, we define the value of CtotalOHk to be +∞ if it corresponds to an infeasible headway.Figure 4
Total cost for train k with respect to headway.Finally, for each traink≥2, we can obtain its optimal headway h∗ (between train k-1 and train k). By rolling optimization of departure headways for trains, the timetable can be obtained. Let H⃑ be the vector of optimal headways for all trains k≥2. The detailed rolling optimization (RO) algorithm is described in Algorithm 2. The arrival and departure times of the first train should be obtained beforehand, and t1,10=0.Algorithm 2:Rolling optimization algorithm.
Step  1: Set H⃑=(0,0,…,0,0). Set the period tp for OST solution approach. Set k=2, h=hmin, and go to Step 2. Step  2: Obtain the segment travel times for train k by OST solution approach. Obtain the arrival and departure times at stations for traink. Go to Step 3. Step  3: Calculate the value of CtotalOHk with respect to h, and calculate the value of CtotalOH′k with respect to h′=h+1. Go to Step4. Step  4: If CtotalOH′k>CtotalOHk, then go to Step 5; Otherwise, set h=h′, and return to Step 2. Step  5: Set the optimal headway for train k as h∗=h, and set H⃑k-1=h∗. Go to Step 6. Step  6: If H⃑≥tendp, then the rolling optimization is finished, and thus the timetable for the period [0,tendt] can be obtained. Otherwise, setk=k+1, h=hmin, and return to Step 2.


#### 4.3.2. Other Approaches

As pointed out by Wang et al. [3], other approaches such as a pattern search [51] or a genetic algorithm [36] can be applied to solve the TTP. Here we briefly describe the generic structures of these two approaches.A genetic algorithm (GA) is an iterative heuristic that seeks an optimal individual among a population of solutions in each generation, and uses selection, crossover, mutation and possibly other operations to obtain new generations. Such process is repeated until the termination condition is satisfied. Since the size of the chromosome should be preset, the frequency is considered here, a given and unchangeable input. Then, the segment travel times can be obtained with the given frequency. An individual consists of headways for all trains (K-1 headways), and its objective function is calculated based on (31). The initial individuals are generated as the cyclic headways. Random values are then added or subtracted to the headways, considering constraints (46) and (49). The details of standard GA [3, 36] are not shown here.Compared with the proposed solution approach (RO) in this paper, generic GA (GGA) has following properties:(1) In terms of input, both GGA and RO require passenger demand as input. However, GGA requires fixed frequency as input, while RO determines optimal frequency during the process. (2) In terms of optimization process, GGA optimizes the solution (of all trains) in a random direction, while RO optimizes one train at a time. (3) In terms of computation complexity, GGA calculates the objective function of the whole system for each chromosome, while RO calculates the objective of one train at each step. These differences result in different computation performances, which will be further demonstrated in the numerical experiments.For the numerical experiments, we set the parameter values in the generic GA as follows. The population size is 20, the maximum iteration number is 1000, the crossover rate is 0.7, and the mutation rate is 0.1. The algorithm terminates if the maximum number of iterations is reached or the best fitness value in this generation is the same as in the previous 100 generations ahead.The pattern search (PS) algorithm is proposed to solve unconstrained optimization problems [52], which varies one variable (headway) at a time by steps of the same magnitude. When no such increase or decrease in any one variable further decreases the objective function, it decreases the step size and repeats the process until the steps are deemed sufficiently small. Here the unconstrained objective is constructed with the Augmented Lagrangian (AL) method, which adds the constraints as penalty terms to the objective function. The penalty term is the product of a Lagrangian multiplier and a constraint equation (e.g., from g(x)≥0 to gx-y2=0, where y is a new variable). Let the number of constraints be m. Set the multiplier that corresponds to the j-th constraint as ωj, and set a large enough parameter σ. By minimizing the objective function with respect to y (so that the intermediate variable y does not influence the objective function), we can obtain the typical form of the penalty term as(62)12σ∑j=1mmax⁡0,ωj-σ·gjx2-ωj2.The parameter valueωj is optimized in each iterative step, until the termination condition is satisfied. In each step of AL, we use the PS to find a minimum value of the unconstrained objective function. In this way, we integrate AL and PS and define such approach as AL + PS. The same drawback exists that the frequency is preset and considered as input. The initial solution is the cyclic headways. For the numerical experiments, we set the values of parameters in the AL + PS as follows. The initial step size for PS is 32 s. The speeding factor for PS is 1. The shrinking factor for PS is 0.5. The allowance error for both PS and AL is 1 s (since the continuous time is discretized into seconds). The augmented factor for AL is 1.2. The large enough parameter σ in AL is 200.


## 4.1. Decomposition

The rationale of the decomposition is explained as follows. The objectives in the optimization model include passenger wait time (at their original stations), passenger travel time (segment travel times and station dwell times), energy, and train cost. Based on (43), the total energy is related to the energy per unit mass and the number of passengers in each segment. The energy per unit mass is solely determined by segment travel times st(i). Since passengers board at their origins and alight at their destinations regardless of which train they take, the sum of in-vehicle passengers in one segment i is determined by passenger demand. In this regard, the total energy can be expressed as(51)Etotal=Mk·K·∑i=12I-1ei+Mp·∑i=12I-1∑j=i+12I∑r=1i∑t=0tendpτr,jt·ei.Therefore, the total energy cost depends only on segment travel timesst(i). The passenger segment travel times also depend only on segment travel times.On the other hand, passenger wait time, dwell time, and train cost depend on the train arrival and departure times at stations. Based on Assumption1, train arrival/departure times at stations are affected by the departure headways from the shunting yard. Thus, it is reasonable to use segment travel times as input and consider the costs for passenger wait time, dwell time, and trains as functions of the departure headways.In general, the OST subproblem is formulated as follows:Minimize(52)CtotalOST=∑i=12I-1ce·Mk·K·ei+ce·Mp·ei+cpt·sti·∑j=i+12I∑r=1i∑t=0tendpτr,jt.Subject to constraints (47) and (48).It should be noted that since segment travel times are given inputs, minimizing the cost of passenger dwell time is equivalent to minimizing the cost of total travel time. Since we already obtain the in-train passenger matrixV, there is no need to further derive dwell times at stations. Thus, the objectives in the OH subproblem include the total passenger time cost (both Q and V) instead of costs of passenger wait and dwell times. The OH subproblem is formulated as follows (with segment travel times as input):Minimize(53)CtotalOH=Ctotalp+Ctotaltra.Subject to constraints (46), (49) and (50).


## 4.2. Solution Approach for OST

The solution approach for OST consists of branch-and-bound (B&B) and frequency determining algorithms. On one hand, the B&B algorithm needs frequency as input. On the other hand, the frequency determining algorithm needs segment travel times obtained by the B&B to evaluate the performance. Thus, we iteratively evaluate the performance with different frequencies (each frequency is used as input for B&B to find a solution of segment travel times) and choose the one with the best performance. The segment travel times are then determined correspondingly.


### 4.2.1. Branch-and-Bound (B&B) Algorithm

The OST subproblem is an integer linear programming problem with2I-1 decision variables (st(i), for i=1,2,…,2I-1). Here we use a generic B&B, which is widely used to solve linear programming problems [7]. Generally, a B&B needs an input model formulated as follows:(54)minc⃑·x⃑T,(55)s.t.A·x⃑T≤b⃑TAeq·x⃑T=beq⃑T,where c⃑ is the coefficient vector and x⃑T is the decision variable vector. Both c⃑ and x⃑T have 2I-1 elements. A and Aeq are coefficient matrixes, and b⃑T and beq⃑T are constant vectors representing the bounds. Here A is a (4I-2)×(2I-1) matrix, and Aeq is a 1×(2I-1) matrix. Correspondingly, b⃑T has 2I-1 elements, and beq⃑T has one element. We define a vector with 2I-1 elements as(56)ϵ⃑j=0,0,…,1,…,0,0,j=1,…,2I-1.The decision variable vector can be represented asx⃑i=st(i) for i=1,2,…,2I-1, where x⃑i is the i-th element in the vector x⃑. We then define ℸi=∑j=i+12I∑r=1i∑t=0tendpτr,jt for i=1,2,…,2I-1 as the total number of passengers that travel through segment i. Based on (52), c⃑i can be represented as(57)c⃑i=ce·Mk·K·χx·sli+ℸi·ce·Mp·χx·sli+cpt;i=1,2,…,2I-1.Constraints (55) are represented as follows. Based on constraint (47), we have Aeq=ϵ⃑I and beq⃑=Tturn. Based on constraint (48), we have Ai,:=-ϵ⃑i, for i=1,2,…,2I-1; Ai,:=ϵ⃑i, for i=2I,2I+1,…,4I-2, where Ai,: represents the i-th row of matrix A. Analogously, b⃑i=-stmini, for i=1,2,…,2I-1; b⃑i=stmaxi, for i=2I,2I+1,…,4I-2. Since we use the standard B&B algorithm [40], the detailed algorithm description is not shown here.


### 4.2.2. Frequency Determining Algorithm

Here we consider cyclic timetables (train departure headways from the shunting yard are constant, e.g., 200 s) to discuss the frequency, because they have following properties:(1) The cyclic timetable with minimum cost (i.e., the best cyclic timetable) can be considered as an upper bound of the noncyclic timetable. (2) If passenger demand is stable, the best cyclic timetable is a good approximation to an optimal noncyclic one. (3) The cyclic timetables are easy to obtain.First, the range of frequency should be estimated. We define the peak passenger volume in one segment asqpeak, which is the maximum value among the numbers of passengers traveling through the segment, that is, qpeak=max⁡(ℸi∣i=1,2,…,2I-1). The minimum frequency Kmin is determined by the train capacity Nk and qpeak:Kmin=qpeak/Nk. The maximum frequency Kmax is determined by the minimum headway: Kmax=tendp/hmin+1.Second, for each frequencyK∈[Kmin,Kmax], we use K as input of the B&B, to obtain the segment travel times. Since the frequency for a cyclic timetable determines the departure headways from the shunting yard (h=tendp/(K-1), where h is the departure headway), the cyclic timetable that corresponds to frequency K can be obtained. Then, the total cost of the timetable can be calculated with (52).Finally, we obtain different cyclic timetables with different costs, each corresponding to a frequencyK∈[Kmin,Kmax]. We select the timetable with minimum cost as the best cyclic timetable, and its corresponding segment travel times as the OST solution within this period. The detailed algorithm is as shown in Algorithm 1.Algorithm 1:Frequency determining algorithm.
Step  1: Estimate the frequency range, i.e., [Kmin,Kmax]. Set the initial frequency of the cyclic timetable as K=Kmin, the best frequency found during the process asKbest=0, and the minimum total cost found during the process as TCmin=+∞. Go to Step 2. Step  2: For frequency K, obtain the segment travel times through the B&B algorithm, and obtain the uniform departure headways from the shunting yardh=tendp/(K-1). Go to Step 3. Step  3: Calculate total cost TCK corresponding to the frequency K. If TCK<TCmin, then set TCmin=TCK, and Kbest=K. Go to Step4. Step  4: If K=Kmax, then go to Step 5. Otherwise, set K=K+1, and return to Step 2. Step  5: Output Kbest and TCmin. The segment travel times within this period are the solution of the B&B algorithm that corresponds to frequencyKbest.


## 4.2.1. Branch-and-Bound (B&B) Algorithm

The OST subproblem is an integer linear programming problem with2I-1 decision variables (st(i), for i=1,2,…,2I-1). Here we use a generic B&B, which is widely used to solve linear programming problems [7]. Generally, a B&B needs an input model formulated as follows:(54)minc⃑·x⃑T,(55)s.t.A·x⃑T≤b⃑TAeq·x⃑T=beq⃑T,where c⃑ is the coefficient vector and x⃑T is the decision variable vector. Both c⃑ and x⃑T have 2I-1 elements. A and Aeq are coefficient matrixes, and b⃑T and beq⃑T are constant vectors representing the bounds. Here A is a (4I-2)×(2I-1) matrix, and Aeq is a 1×(2I-1) matrix. Correspondingly, b⃑T has 2I-1 elements, and beq⃑T has one element. We define a vector with 2I-1 elements as(56)ϵ⃑j=0,0,…,1,…,0,0,j=1,…,2I-1.The decision variable vector can be represented asx⃑i=st(i) for i=1,2,…,2I-1, where x⃑i is the i-th element in the vector x⃑. We then define ℸi=∑j=i+12I∑r=1i∑t=0tendpτr,jt for i=1,2,…,2I-1 as the total number of passengers that travel through segment i. Based on (52), c⃑i can be represented as(57)c⃑i=ce·Mk·K·χx·sli+ℸi·ce·Mp·χx·sli+cpt;i=1,2,…,2I-1.Constraints (55) are represented as follows. Based on constraint (47), we have Aeq=ϵ⃑I and beq⃑=Tturn. Based on constraint (48), we have Ai,:=-ϵ⃑i, for i=1,2,…,2I-1; Ai,:=ϵ⃑i, for i=2I,2I+1,…,4I-2, where Ai,: represents the i-th row of matrix A. Analogously, b⃑i=-stmini, for i=1,2,…,2I-1; b⃑i=stmaxi, for i=2I,2I+1,…,4I-2. Since we use the standard B&B algorithm [40], the detailed algorithm description is not shown here.


## 4.2.2. Frequency Determining Algorithm

Here we consider cyclic timetables (train departure headways from the shunting yard are constant, e.g., 200 s) to discuss the frequency, because they have following properties:(1) The cyclic timetable with minimum cost (i.e., the best cyclic timetable) can be considered as an upper bound of the noncyclic timetable. (2) If passenger demand is stable, the best cyclic timetable is a good approximation to an optimal noncyclic one. (3) The cyclic timetables are easy to obtain.First, the range of frequency should be estimated. We define the peak passenger volume in one segment asqpeak, which is the maximum value among the numbers of passengers traveling through the segment, that is, qpeak=max⁡(ℸi∣i=1,2,…,2I-1). The minimum frequency Kmin is determined by the train capacity Nk and qpeak:Kmin=qpeak/Nk. The maximum frequency Kmax is determined by the minimum headway: Kmax=tendp/hmin+1.Second, for each frequencyK∈[Kmin,Kmax], we use K as input of the B&B, to obtain the segment travel times. Since the frequency for a cyclic timetable determines the departure headways from the shunting yard (h=tendp/(K-1), where h is the departure headway), the cyclic timetable that corresponds to frequency K can be obtained. Then, the total cost of the timetable can be calculated with (52).Finally, we obtain different cyclic timetables with different costs, each corresponding to a frequencyK∈[Kmin,Kmax]. We select the timetable with minimum cost as the best cyclic timetable, and its corresponding segment travel times as the OST solution within this period. The detailed algorithm is as shown in Algorithm 1.Algorithm 1:Frequency determining algorithm.
Step  1: Estimate the frequency range, i.e., [Kmin,Kmax]. Set the initial frequency of the cyclic timetable as K=Kmin, the best frequency found during the process asKbest=0, and the minimum total cost found during the process as TCmin=+∞. Go to Step 2. Step  2: For frequency K, obtain the segment travel times through the B&B algorithm, and obtain the uniform departure headways from the shunting yardh=tendp/(K-1). Go to Step 3. Step  3: Calculate total cost TCK corresponding to the frequency K. If TCK<TCmin, then set TCmin=TCK, and Kbest=K. Go to Step4. Step  4: If K=Kmax, then go to Step 5. Otherwise, set K=K+1, and return to Step 2. Step  5: Output Kbest and TCmin. The segment travel times within this period are the solution of the B&B algorithm that corresponds to frequencyKbest.


## 4.3. Solution Approach for OH

Here we design a rolling optimization (RO) algorithm that optimizes the arrival time at the first station for one train at a time. The main difference between the RO and the rolling horizon (RH) approach is in the RH, the bound conditions (i.e., the start and end times of each horizon) require some of the original variables to have fixed, known values [3], while in the RO, each decision variable may still be optimized. Thus, the rolling optimization is more flexible.


### 4.3.1. The Rolling Optimization (RO) Algorithm

Based on (53), the objectives include the total costs of passenger time and trains. To obtain these costs, we need the segment travel times as input. Here, instead of considering the whole period, we consider a much shorter period tp≤hmin+∑i=12I-1st(i)+∑i=12Idi¯, e.g., 600 s, for solving the OST. For example, for each train k≥2, its previous train arrives at the first station at time tk-1,10 (the first train arrives at time 0), and we consider the period [tk-1,10,tk-1,10+tp] to solve the OST (the constraint of tp ensures that tk-1,10+tp≤tendt). In this way, the key ideas of “iterative optimization” (the frequencies, segment travel times, and headways are interdependent) and “rolling horizon” (the smaller period studied in each step) are adopted, and future passenger demand is considered.Then, we discuss the costs for one traink≥2 if its arrival time at the first station is given. Since each train is scheduled in chronological order, the waiting passenger matrix Q and the in-train passenger matrix V can be constructed chronologically. For each train k≥2, the departure time of its preceding train at the first station is tk-1,1b, and train k departs from the first station at time tk,1b. We define the passenger wait time with respect to train k as Tpw(k), which is formulated as(58)Tpwk=∑i=12I∑t=tk-1,1btk,1bQi,t;∀k∈K.Analogously, we define the passenger travel time for traink as Tpt(k), and we obtain(59)Tptk=∑t=tk,10tk,2IbVk,t;∀k∈K.The estimated fleet size of traink (Fsk) depends on its round trip time (Troundk) and the departure headway from the shunting yard (tk,10-tk-1,10):(60)Fsk=Troundktk,10-tk-1,10;∀k∈K.Based on (18), (30), (53), (58), (59), and (60), we define the total costs in the OH for train k as CtotalOH(k), which can be expressed as(61)CtotalOHk=cpw·Tpwk+cpt·Tptk+ctra·Fsk·Troundk3600+co·∑i=12I-1sli1000;∀k∈K.The value ofCtotalOHk is a function of the arrival time for train k at the first station. Since the arrival time for the previous train is given, we consider the value of CtotalOHk as a function of the headway h between train k-1 and train k (h=tk,10-tk-1,10), for consistency with other TTP studies.Next we discuss how the optimal value forh can be obtained. Based on (58), (59), and (60), both the passenger time and train cost depend on passenger demand. Since the passenger demand is time-varying, the relation between CtotalOHk and h cannot be described with a fixed function.The expected objective function of OH is shown in Figure4. Its value first decreases (the headway starts from hmin) to a minimum value and then increases monotonically. We specify that this function reaches its minimum value at headway h∗. Since the DTSG divided the continuous time into seconds, we iteratively calculate the objective function of OH with respect to h (starting from h=hmin, and h=h+1 in each following step), until the optimal headway h∗ is found, or the total cost increases as h increases (in case that h∗ cannot be reached). It should be noted that sometimes a headway is not feasible due to the model constraints (e.g., a short headway might violate the minimum headway constraint). In such cases, we define the value of CtotalOHk to be +∞ if it corresponds to an infeasible headway.Figure 4
Total cost for train k with respect to headway.Finally, for each traink≥2, we can obtain its optimal headway h∗ (between train k-1 and train k). By rolling optimization of departure headways for trains, the timetable can be obtained. Let H⃑ be the vector of optimal headways for all trains k≥2. The detailed rolling optimization (RO) algorithm is described in Algorithm 2. The arrival and departure times of the first train should be obtained beforehand, and t1,10=0.Algorithm 2:Rolling optimization algorithm.
Step  1: Set H⃑=(0,0,…,0,0). Set the period tp for OST solution approach. Set k=2, h=hmin, and go to Step 2. Step  2: Obtain the segment travel times for train k by OST solution approach. Obtain the arrival and departure times at stations for traink. Go to Step 3. Step  3: Calculate the value of CtotalOHk with respect to h, and calculate the value of CtotalOH′k with respect to h′=h+1. Go to Step4. Step  4: If CtotalOH′k>CtotalOHk, then go to Step 5; Otherwise, set h=h′, and return to Step 2. Step  5: Set the optimal headway for train k as h∗=h, and set H⃑k-1=h∗. Go to Step 6. Step  6: If H⃑≥tendp, then the rolling optimization is finished, and thus the timetable for the period [0,tendt] can be obtained. Otherwise, setk=k+1, h=hmin, and return to Step 2.


### 4.3.2. Other Approaches

As pointed out by Wang et al. [3], other approaches such as a pattern search [51] or a genetic algorithm [36] can be applied to solve the TTP. Here we briefly describe the generic structures of these two approaches.A genetic algorithm (GA) is an iterative heuristic that seeks an optimal individual among a population of solutions in each generation, and uses selection, crossover, mutation and possibly other operations to obtain new generations. Such process is repeated until the termination condition is satisfied. Since the size of the chromosome should be preset, the frequency is considered here, a given and unchangeable input. Then, the segment travel times can be obtained with the given frequency. An individual consists of headways for all trains (K-1 headways), and its objective function is calculated based on (31). The initial individuals are generated as the cyclic headways. Random values are then added or subtracted to the headways, considering constraints (46) and (49). The details of standard GA [3, 36] are not shown here.Compared with the proposed solution approach (RO) in this paper, generic GA (GGA) has following properties:(1) In terms of input, both GGA and RO require passenger demand as input. However, GGA requires fixed frequency as input, while RO determines optimal frequency during the process. (2) In terms of optimization process, GGA optimizes the solution (of all trains) in a random direction, while RO optimizes one train at a time. (3) In terms of computation complexity, GGA calculates the objective function of the whole system for each chromosome, while RO calculates the objective of one train at each step. These differences result in different computation performances, which will be further demonstrated in the numerical experiments.For the numerical experiments, we set the parameter values in the generic GA as follows. The population size is 20, the maximum iteration number is 1000, the crossover rate is 0.7, and the mutation rate is 0.1. The algorithm terminates if the maximum number of iterations is reached or the best fitness value in this generation is the same as in the previous 100 generations ahead.The pattern search (PS) algorithm is proposed to solve unconstrained optimization problems [52], which varies one variable (headway) at a time by steps of the same magnitude. When no such increase or decrease in any one variable further decreases the objective function, it decreases the step size and repeats the process until the steps are deemed sufficiently small. Here the unconstrained objective is constructed with the Augmented Lagrangian (AL) method, which adds the constraints as penalty terms to the objective function. The penalty term is the product of a Lagrangian multiplier and a constraint equation (e.g., from g(x)≥0 to gx-y2=0, where y is a new variable). Let the number of constraints be m. Set the multiplier that corresponds to the j-th constraint as ωj, and set a large enough parameter σ. By minimizing the objective function with respect to y (so that the intermediate variable y does not influence the objective function), we can obtain the typical form of the penalty term as(62)12σ∑j=1mmax⁡0,ωj-σ·gjx2-ωj2.The parameter valueωj is optimized in each iterative step, until the termination condition is satisfied. In each step of AL, we use the PS to find a minimum value of the unconstrained objective function. In this way, we integrate AL and PS and define such approach as AL + PS. The same drawback exists that the frequency is preset and considered as input. The initial solution is the cyclic headways. For the numerical experiments, we set the values of parameters in the AL + PS as follows. The initial step size for PS is 32 s. The speeding factor for PS is 1. The shrinking factor for PS is 0.5. The allowance error for both PS and AL is 1 s (since the continuous time is discretized into seconds). The augmented factor for AL is 1.2. The large enough parameter σ in AL is 200.


## 4.3.1. The Rolling Optimization (RO) Algorithm

Based on (53), the objectives include the total costs of passenger time and trains. To obtain these costs, we need the segment travel times as input. Here, instead of considering the whole period, we consider a much shorter period tp≤hmin+∑i=12I-1st(i)+∑i=12Idi¯, e.g., 600 s, for solving the OST. For example, for each train k≥2, its previous train arrives at the first station at time tk-1,10 (the first train arrives at time 0), and we consider the period [tk-1,10,tk-1,10+tp] to solve the OST (the constraint of tp ensures that tk-1,10+tp≤tendt). In this way, the key ideas of “iterative optimization” (the frequencies, segment travel times, and headways are interdependent) and “rolling horizon” (the smaller period studied in each step) are adopted, and future passenger demand is considered.Then, we discuss the costs for one traink≥2 if its arrival time at the first station is given. Since each train is scheduled in chronological order, the waiting passenger matrix Q and the in-train passenger matrix V can be constructed chronologically. For each train k≥2, the departure time of its preceding train at the first station is tk-1,1b, and train k departs from the first station at time tk,1b. We define the passenger wait time with respect to train k as Tpw(k), which is formulated as(58)Tpwk=∑i=12I∑t=tk-1,1btk,1bQi,t;∀k∈K.Analogously, we define the passenger travel time for traink as Tpt(k), and we obtain(59)Tptk=∑t=tk,10tk,2IbVk,t;∀k∈K.The estimated fleet size of traink (Fsk) depends on its round trip time (Troundk) and the departure headway from the shunting yard (tk,10-tk-1,10):(60)Fsk=Troundktk,10-tk-1,10;∀k∈K.Based on (18), (30), (53), (58), (59), and (60), we define the total costs in the OH for train k as CtotalOH(k), which can be expressed as(61)CtotalOHk=cpw·Tpwk+cpt·Tptk+ctra·Fsk·Troundk3600+co·∑i=12I-1sli1000;∀k∈K.The value ofCtotalOHk is a function of the arrival time for train k at the first station. Since the arrival time for the previous train is given, we consider the value of CtotalOHk as a function of the headway h between train k-1 and train k (h=tk,10-tk-1,10), for consistency with other TTP studies.Next we discuss how the optimal value forh can be obtained. Based on (58), (59), and (60), both the passenger time and train cost depend on passenger demand. Since the passenger demand is time-varying, the relation between CtotalOHk and h cannot be described with a fixed function.The expected objective function of OH is shown in Figure4. Its value first decreases (the headway starts from hmin) to a minimum value and then increases monotonically. We specify that this function reaches its minimum value at headway h∗. Since the DTSG divided the continuous time into seconds, we iteratively calculate the objective function of OH with respect to h (starting from h=hmin, and h=h+1 in each following step), until the optimal headway h∗ is found, or the total cost increases as h increases (in case that h∗ cannot be reached). It should be noted that sometimes a headway is not feasible due to the model constraints (e.g., a short headway might violate the minimum headway constraint). In such cases, we define the value of CtotalOHk to be +∞ if it corresponds to an infeasible headway.Figure 4
Total cost for train k with respect to headway.Finally, for each traink≥2, we can obtain its optimal headway h∗ (between train k-1 and train k). By rolling optimization of departure headways for trains, the timetable can be obtained. Let H⃑ be the vector of optimal headways for all trains k≥2. The detailed rolling optimization (RO) algorithm is described in Algorithm 2. The arrival and departure times of the first train should be obtained beforehand, and t1,10=0.Algorithm 2:Rolling optimization algorithm.
Step  1: Set H⃑=(0,0,…,0,0). Set the period tp for OST solution approach. Set k=2, h=hmin, and go to Step 2. Step  2: Obtain the segment travel times for train k by OST solution approach. Obtain the arrival and departure times at stations for traink. Go to Step 3. Step  3: Calculate the value of CtotalOHk with respect to h, and calculate the value of CtotalOH′k with respect to h′=h+1. Go to Step4. Step  4: If CtotalOH′k>CtotalOHk, then go to Step 5; Otherwise, set h=h′, and return to Step 2. Step  5: Set the optimal headway for train k as h∗=h, and set H⃑k-1=h∗. Go to Step 6. Step  6: If H⃑≥tendp, then the rolling optimization is finished, and thus the timetable for the period [0,tendt] can be obtained. Otherwise, setk=k+1, h=hmin, and return to Step 2.


## 4.3.2. Other Approaches

As pointed out by Wang et al. [3], other approaches such as a pattern search [51] or a genetic algorithm [36] can be applied to solve the TTP. Here we briefly describe the generic structures of these two approaches.A genetic algorithm (GA) is an iterative heuristic that seeks an optimal individual among a population of solutions in each generation, and uses selection, crossover, mutation and possibly other operations to obtain new generations. Such process is repeated until the termination condition is satisfied. Since the size of the chromosome should be preset, the frequency is considered here, a given and unchangeable input. Then, the segment travel times can be obtained with the given frequency. An individual consists of headways for all trains (K-1 headways), and its objective function is calculated based on (31). The initial individuals are generated as the cyclic headways. Random values are then added or subtracted to the headways, considering constraints (46) and (49). The details of standard GA [3, 36] are not shown here.Compared with the proposed solution approach (RO) in this paper, generic GA (GGA) has following properties:(1) In terms of input, both GGA and RO require passenger demand as input. However, GGA requires fixed frequency as input, while RO determines optimal frequency during the process. (2) In terms of optimization process, GGA optimizes the solution (of all trains) in a random direction, while RO optimizes one train at a time. (3) In terms of computation complexity, GGA calculates the objective function of the whole system for each chromosome, while RO calculates the objective of one train at each step. These differences result in different computation performances, which will be further demonstrated in the numerical experiments.For the numerical experiments, we set the parameter values in the generic GA as follows. The population size is 20, the maximum iteration number is 1000, the crossover rate is 0.7, and the mutation rate is 0.1. The algorithm terminates if the maximum number of iterations is reached or the best fitness value in this generation is the same as in the previous 100 generations ahead.The pattern search (PS) algorithm is proposed to solve unconstrained optimization problems [52], which varies one variable (headway) at a time by steps of the same magnitude. When no such increase or decrease in any one variable further decreases the objective function, it decreases the step size and repeats the process until the steps are deemed sufficiently small. Here the unconstrained objective is constructed with the Augmented Lagrangian (AL) method, which adds the constraints as penalty terms to the objective function. The penalty term is the product of a Lagrangian multiplier and a constraint equation (e.g., from g(x)≥0 to gx-y2=0, where y is a new variable). Let the number of constraints be m. Set the multiplier that corresponds to the j-th constraint as ωj, and set a large enough parameter σ. By minimizing the objective function with respect to y (so that the intermediate variable y does not influence the objective function), we can obtain the typical form of the penalty term as(62)12σ∑j=1mmax⁡0,ωj-σ·gjx2-ωj2.The parameter valueωj is optimized in each iterative step, until the termination condition is satisfied. In each step of AL, we use the PS to find a minimum value of the unconstrained objective function. In this way, we integrate AL and PS and define such approach as AL + PS. The same drawback exists that the frequency is preset and considered as input. The initial solution is the cyclic headways. For the numerical experiments, we set the values of parameters in the AL + PS as follows. The initial step size for PS is 32 s. The speeding factor for PS is 1. The shrinking factor for PS is 0.5. The allowance error for both PS and AL is 1 s (since the continuous time is discretized into seconds). The augmented factor for AL is 1.2. The large enough parameter σ in AL is 200.


## 5. Numerical Experiments

In this section, we present several numerical experiments to test the effectiveness and efficiency of the proposed solution approach. All experiments are performed with Matlab on a PC with 2 GHz Intel Core7 and 8 GB memory. Parameters in the examples are shown in Table3. We adopt here the same energy calculation method as in [48]. To save space, we do not show the detailed calculations and directly provide the estimated values of parameters χx and χy in Table 3.Table 3
Parameter values in the experiments.
Parameter Value M k (Kg) 1.99 × 10 5 M p  (Kg) 75 N k (Passenger) 1290 c e (¥/Kwh) 0.79 c o (¥ /(train⋅Km)) 20 c p t (¥/(passenger⋅hour)) 10 c p w (¥/(passenger⋅hour)) 20 c t r a (¥/(train⋅hour)) 800 h m i n (s) 100 f r (m/s) 1.36 × 10 - 4 ⋅ v k 2 t + 1.45 × 10 - 2 ⋅ v k t + 0.08 χ x - 1.944 × 10 - 4 χ y 0.4176


### 5.1. A Small Case Study

Here we consider a short metro line, as shown in Figure5. The segment lengths are 1800 m, 1600 m, and 2000 m. For simplicity, the minimum and maximum dwell times for each station are all set as 30 s and 90 s, respectively.Figure 5
Representation of a short metro line.The baseline passenger demand is shown in Table4. The “Boarding (Alighting) Sum” represents the number of total boarding (alighting) passengers per second at a station. Four cases are considered here. As shown in Figure 6, the curves represent demand variance over time, and “multiplier” represents the ratio between actual and baseline passenger demand. For example, in case 1, during period [1200,2400] (s), the multiplier is 1.6, and passenger demand from station 2 to station 3 at t∈[1200,2400] is 0.96 (0.6×1.6) passengers/second. The four cases have different passenger demand distributions.Table 4
Baseline passenger demands for a short metro line (passengers per second).
Station index 1 2 3 4 5 6 7 8 Boarding sum 1 0 0.3 0.3 0.45 0 0 0 0 1.05 2 0 0 0.6 0.3 0 0 0 0 0.9 3 0 0 0 0.3 0 0 0 0 0.3 4 0 0 0 0 0 0 0 0 0 5 0 0 0 0 0 0.45 0.75 0.3 1.5 6 0 0 0 0 0 0 0.6 0.3 0.9 7 0 0 0 0 0 0 0 0.45 0.45 8 0 0 0 0 0 0 0 0 0 Alighting sum 0 0.3 0.9 1.05 0 0.45 1.35 1.05 -Figure 6
Passenger demand distributions in different cases.The computational results of different cases are shown in Table5 and Figure 7. Specifically, in Table 5, the noncyclic result obtained by RO is denoted as NK=KN, which schedules KN trains within a given time duration. The cyclic result obtained by Algorithm 1 is denoted as UK=KU, which schedules KU trains within a given time duration. Particularly, the best cyclic result is denoted as UK=KUbest. The GGA result is denoted as GAK=KN, which uses frequency KN from RO result as input. The objectives of different results are calculated, including passenger wait time (Tpw), passenger travel time (Tpt), total energy (Etotal), fleet size (Fs), and total cost (Ctotal). Figure 7 specifically shows the numerical results in passenger satisfaction.Table 5
Comparison between cyclic and non-cyclic timetables.
Time duration (s) Timetable type  T p w  (106) T p t  (106) E t o t a l   (103) F s C t o t a l   (104) Gap (%) 0-10800 Non-cyclic N K = 60 7.2330 14.384 18.069 7 12.320 - Case 1 Non-cyclic G A K = 60 7.2786 14.402 18.069 7 12.329 0.07 Cyclic U K = 53 b e s t 8.3360 14.523 16.148 6 12.496 1.41 Cyclic U K = 60 7.3561 14.406 18.069 7 + ∞ -  0-4800 Non-cyclic N K = 24 2.3171 4.9470 7.1411 6 4.3586 - Case 2 Non-cyclic G A K = 24 2.5489 4.9758 7.1411 6 4.4954 3.04 Cyclic U K = 24 b e s t 2.6102 4.9988 7.1411 6 4.5407 4.01 Cyclic U K = 23 2.7619 5.0347 6.8667 6 4.5681 4.59 Cyclic U K = 25 2.4989 4.9859 7.4156 7 4.5441 4.08  0-7200 Non-cyclic N K = 45 4.3948 11.292 13.615 8 8.8022 - Case 3 Non-cyclic G A K = 45 4.5726 11.328 13.615 8 8.9245 1.37 Cyclic U K = 48 b e s t 4.3975 11.319 14.438 8 9.0243 2.46 Cyclic U K = 45 4.7268 11.417 13.615 8 9.0345 2.57  0-10800 Non-cyclic N K = 60 6.0280 14.301 18.069 7 11.621 - Case 4 Non-cyclic G A K = 60 6.3113 14.560 18.069 7 11.854 1.97 Cyclic U K = 60 b e s t 6.4638 14.330 18.069 7 11.875 2.14 Cyclic U K = 59 6.5981 14.381 17.795 7 11.899 2.34 Cyclic U K = 61 6.3572 14.327 18.344 7 11.884 2.21Figure 7
Passenger costs for different results in different cases.Table5 and Figure 7 demonstrate that, under time-varying passenger demand, noncyclic timetables obtained by RO outperform the other timetables, while GGA results outperform cyclic ones obtained by Algorithm 1. This result verifies that the noncyclic timetable improves passenger satisfaction. Particularly, the passenger wait time (Tpw), passenger travel time (Tpt), and total cost (Ctotal) have their lowest values in the RO results, which demonstrates that the RO method is effective in solving the proposed train timetabling problem. In case 1, the result of UK=60 is infeasible because it violates the minimum headway constraint at intermediate stations (the dwell times are different from UK=60 to NK=60), and thus we have Ctotal=+∞. It is interesting that NK=60 and UK=60 require the same energy and fleet size. This result verifies the reasonableness of model decomposition in this paper from two perspectives: (1) When all passengers reach their destinations, energy is solely determined by segment travel times (they are the same in NK=60 and UK=60). (2) With given segment travel times, the estimated fleet size is solely determined by headways. Similar results can be observed in other cases.The gaps in Table5 shows that case 2 has the largest gap, case 3 and 4 have similar gaps, and case 1 has the smallest gap. Note that, in case 2, passenger demand changes drastically; in cases 3 and 4, demand patterns are quite similar; in case 1, within certain subperiods, passenger demand is the same, which is the steadiest among all cases. Therefore, we can conclude that as the passenger demand varies more significantly over time, the gap between the cyclic and noncyclic timetables increases.Figure8 illustrates the obtained train timetables in different cases, where the departure headways vary in accordance with the fluctuating passenger demand. For example, in cases 2, 3, and 4, the headways are small in the peak period and are larger before and after the peak period. In case 1, the small headways occur during subperiods [1200,2400] (s) and [7200,8400] (s), where the multipliers of passenger demand are largest (1.6 and 1.8). In Figure 8, the circles represent the dwell times that exceed the minimum dwell time (i.e., 30 s). We observe that the excess dwell times (circles) mostly occur at station 7. This is consistent with the OD pairs we present in Table 4. Specifically, the difference between the cumulative boarding sum and the cumulative alighting sum reaches its maximum value at station 7 (1.5+0.9-0.45=1.95). In this regard, the alighting and boarding rates are smallest at station 7. In addition, the total number of alighting and boarding passengers is largest (i.e., 1.8 passengers per second) at station 7, and thus the excess dwell times here are reasonable.Figure 8
Obtained timetables for different cases.
(a)
Obtained timetable for case 1 (NK=60) (b)
Obtained timetable for case 2 (NK=24) (c)
Obtained timetable for case 3 (NK=45) (d)
Obtained timetable for case 4 (NK=60)Moreover, despite the symmetry of passenger demand distributions, the headways at the ends of periods are larger than at their starts. This can be explained by that passenger demand increases at the start and decreases to zero at the end, which leads to different cost changes. Specifically, the passenger time cost increases more significantly when demand increases, while the train cost does not change with demand. Thus, the headway tends to be smaller at the start of periods than at the ends. In fact, it is often seen in actual operations that the headways between the last trains at the end of day exceed those between the first trains in the morning.Finally, Table6 shows the computational performances for different methods, including RO method, GGA, and AL + PS method. As described above, GGA and AL + PS require a preset fixed frequency, and thus we use the same frequency from RO result KN as input. We terminate the algorithm when its computation time reaches 3 hours and output the best solution obtained by then.Table 6
Computational performances for different methods.
RO method Generic GA AL + PS (Case 1) Time (s) 225.5 8299.4 10800 Total cost (¥) 12.32 × 10 4 12.33 × 10 4 12.42 × 10 4  (Case 2) Time (s) 86.1 4658.2 10800 Total cost (¥) 4.3586 × 10 4 4.4954 × 10 4 4.4872 × 10 4  (Case 3) Time (s) 120.6 5009.3 10800 Total cost (¥) 8.8022 × 10 4 8.9245 × 10 4 8.9567 × 10 4  (Case 4) Time (s) 240.717 8229.1 10800 Total cost (¥) 11.621 × 10 4 11.854 × 10 4 11.774 × 10 4Table6 demonstrates that RO is very efficient, obtaining a timetable very fast. In comparison, GGA takes more than 1.5 hours to obtain a timetable, and AL + PS takes more than 3 hours in all cases. In each case, the timetable obtained with the RO method has the lowest total cost, which further verifies its effectiveness. The GGA timetables perform better than AL + PS in cases 1 and 3 but are outperformed in cases 2 and 4.Combining all the computational results in this case study, we can conclude some disadvantages of GGA and AL + PS compared to RO.(1) Although the timetables obtained by them are noncyclic, train departure headways in their results differ greatly from one another without clear correspondence to the passenger demand. In comparison, RO timetables are more regular (headways change gradually) and capture well the passenger demand patterns. (2) Due to the computation complexity, the computation time for running GGA and AL + PS once is far greater than RO. (3) They both require a preset frequency as input. If we do not know the best frequency beforehand, which is true in most real-world scenarios, we must run these methods with different frequencies and select the best solution, thereby further increasing substantially computation time.In general, the RO method is superior to GGA and AL + PS in obtaining a noncyclic timetable with time-varying passenger demand.


### 5.2. Case Study of Beijing Metro Line 4

This section aims to test the applicability and effectiveness of RO in a real-world case study based on Beijing’s Metro Line 4 (BML4). The infrastructure information and some operational parameters of BML4 are shown in Table7. The turnaround time is 120s. To satisfy the large passenger demand, we set train capacity as Nk=2160 passengers, which is a maximum overloaded capacity. If not mentioned, the other parameters are the same as in Table 3. The baseline passenger demand data are based on the research conducted by Duan et al. [53], which are presented in Figure 9, where zx,y, ∀x<y represents morning peak demand and zx,y, ∀x>y represents evening peak demand.Table 7
Real-world infrastructure data of Beijing Metro Line 4.
Station Minimum dwell time (s) Maximum dwell time (s) Segment Minimum travel time (s) Maximum travel time (s) Segment length (m) 1/48 AHQB 30 90 1 124 170 1363 2/47 BGM 30 90 2 114 156 1251 3/46 XY 30 90 3 152 209 1672 4/45 YMY 30 90 4 118 162 1295 5/44 BDDM 30 90 5 81 111 887 6/43 ZGC 30 90 6 82 113 900 7/42 HDHZ 30 90 7 97 133 1063 8/41 RMDX 30 90 8 96 131 1051 9/40 WGC 30 90 9 151 207 1658 10/39 NL 30 90 10 138 190 1517 11/38 BJZ 30 90 11 131 180 1441 12/37 XZM 30 90 12 93 128 1025 13/36 XJK 30 90 13 100 138 1100 14/35 PAL 30 90 14 100 138 1100 15/34 XS 30 90 15 79 109 869 16/33 LJHT 30 90 16 92 126 1011 17/32 XD 30 90 17 74 102 815 18/31 XWM 30 90 18 105 144 1152 19/30 CSK 30 90 19 109 150 1200 20/29 TRT 30 90 20 149 205 1643 21/28 BJSRS 30 90 21 135 185 1480 22/27 MJP 30 90 22 75 103 827 23/26 JMX 30 90 23 90 124 989 24/25 GYXQ 30 90 - - - -Figure 9
Baseline passenger demand for BML4.Here we consider three types of scenarios whose passenger demands are large and unsteady, which are difficult to manage: (a) period of 3 hours (i.e., 10800 s) including morning peak; (b) period of 3 hours including evening peak; and (c) a shopping festival hosted in the shopping centers next to the XD station (whose index is 8/41). Representations of passenger demands are the same as in the small case study (as shown in Figure10): the actual passenger demand is the product of baseline passenger demand (Figure 9) and the multipliers (Figure 10). To test RO with more instances, we design three subcases of scenario (a) with different peak multipliers, and two subcases of scenario (c): without the special event (c-1) or with it (c-2). In case (c-2), the passenger demand at station XD is tripled.Figure 10
Passenger demand patterns in different cases.It should be noted that the real operation data of BML4 is dynamic and confidential. Besides, the operator today prefers to schedule trains uniformly. Therefore, in this experiment, we use cyclic timetables with best performance to approximate the real-world timetables. The notations are the same as in the small case study.Table8 and Figure 11 show the computational results in different cases. Note that, for a major metro line such as BML4, the computation times of GGA and AL + PS are extremely large, and in many cases they do not yield a satisfactory solution within a reasonable time. Thus, we do not present their results in this experiment. The results of UK=46 in case (a-1), UK=54 in case (a-2), and UK=54 in case (a-3) are infeasible, and thus their total costs are +∞. Figure 11 shows that considerable costs are reduced by the noncyclic timetables obtained by RO. Interestingly, in the cases where cyclic timetables are feasible, that is, cases (a-1), (b), (c-1), and (c-2), the frequency of the noncyclic timetable is always larger than the best cyclic ones. This is because the passenger demand in the real-world fluctuates so drastically that cyclic timetables cannot capture its dynamic patterns, and the minimum headway constraint makes it hard to obtain feasible cyclic timetables. In this regard, adding more trains in a cyclic timetable may be either uneconomic (i.e., the increased train operating cost exceeds the reduced passenger time value) or infeasible.Table 8
Computational results for BML4.
Scenario   T p w  (108) T p t  (108) E t o t a l  (104) F s C t o t a l  (106) Gap (%) Case (a-1) Non-cyclic N K = 46 3.8665 3.0555 7.4898 31 3.1788 - 3869 s Cyclic U K = 44 b e s t 5.3237 3.0943 7.2385 30 3.9925 20.38 Cyclic U K = 46 4.9805 3.0985 7.4898 31 + ∞ -  Case (a-2) Non-cyclic N K = 54 11.426 4.3898 9.2418 37 7.7869 - 1544 s Cyclic U K = 54 10.653 4.4060 9.2418 37 + ∞ -  Case (a-3) Non-cyclic N K = 54 19.568 4.7953 9.4566 38 12.426 - 1133 s Cyclic U K = 54 18.891 4.816 9.4958 38 + ∞ -  Case (b) Non-cyclic N K = 38 2.2416 2.1501 5.9624 25 1.9902 - 3932 s Cyclic U K = 30 b e s t 4.4170 2.1632 4.9379 21 3.1748 37.31  Case (c-1) Non-cyclic N K = 30 0.5482 1.0256 4.3421 19 0.7004 - 6848 s Cyclic U K = 22 b e s t 1.2744 1.0946 3.3695 14 1.0957 36.08  Case (c-2) Non-cyclic N K = 41 1.1415 1.5061 6.0564 26 1.2342 - 4068 s Cyclic U K = 25 b e s t 3.8100 1.5533 3.9911 17 2.6459 53.53Figure 11
Total costs for different results in different cases.In comparison, the noncyclic timetables obtained by RO respond well to demand changes. For example, in case (a-1), the result ofNK=46 substantially reduces passenger wait time compared to UK=44best and UK=46. The adaptation for time-varying demand is further demonstrated in cases (c-1) and (c-2). In case (c-1), the demand distribution pattern is similar to cases (a-1) and (b), and the cost reduction of the noncyclic timetable is closer to the other cases (36.08%). In case (c-2), the special event leads to a more heterogeneous demand distribution, and the cost reduction reaches 53.35%. Besides, the computation times (shown below each case) are acceptable for such a major and busy metro line.However, in cases (a-2) and (a-3), the passenger wait times for noncyclic timetables are larger than for cyclic ones. This occurs because the RO method considers the objective function that corresponds to an infeasible headway as+∞. With a large demand, the first departure headway that is feasible may greatly exceed h∗ (as can be seen in Figure 4), whereas the cyclic timetable does not consider the feasibility and divides headways uniformly. In other words, the increase of departure headways (leading to larger Tpw) is due to safety issues, which is consistent with the actual operational requirements that put safety ahead of costs. In addition, with parameters given in this experiment, the largest frequency for BML4 is KN=54, as can be concluded by comparing cases (a-2) and (a-3). Therefore, real-time dispatching and passenger flow control methods should be applied with passenger demand at or above the levels in case (a-2).Figures12, 13, and 14 illustrate the noncyclic timetables obtained in case (a-1), case (b), and case (c-2), respectively. The circles here represent the occurrences of maximum dwell times. These figures illustrate that the headways vary in accordance with demand patterns, which is consistent with the results in the small case study. Particularly in Figure 14, the dense headways in the middle are set to transport many passengers due to the special event. This shows the flexibility of noncyclic timetables in special scenarios.Figure 12
Timetable for morning peak (case a-1).Figure 13
Timetable for evening peak (case b).Figure 14
Timetable for a special event (case c-2).In addition, from the distribution of the circles, we can see that case (a-1) is the busiest scenario among these figures. The many circles in case (a-1) help explain the results in cases (a-2) and (a-3) (which have larger demands than case (a-1)). With so many trains taking maximum dwell times to load passengers, real-time dispatching and passenger flow control methods are necessary. In case (b), the maximum dwell times mostly occur at stations 7, 11, and 12. Specifically, the neighborhood of station 7 (HDHZ) includes famous schools, fashion shopping streets, and dense office buildings. The neighborhood of station 11 (BJZ) includes a major bus hub, the Beijing Zoo, and one of the largest wholesale markets. The neighborhood of station 12 (XZM), which has transfers with two other busy metro lines, includes a railway station, large shopping malls, and office buildings. They are well known for the large numbers of passengers during peak hours. Thus, we find that the model is reasonable and consistent with real-world scenarios. Analogously, in case (c-2), the circles are mostly seen at stations near the station XD, where the special event is held.It should be noted that the last two trains have an unusually small headway in case (b). The reason is that we assume that the last train must arrive at the first station at the end of the period. It is expected that if there is still demand after 10800 s, the headways will be optimized by our approaches.


### 5.3. Additional Experiments

In this section, we implement three sets of additional experiments to investigate the relations between the values of costs and the objectives. The first experiment tests the sensitivity of passenger time and peak frequency to unit passenger time cost, based on the demand in case (b) of BML4 study. In practice, the unit time value of passenger travel time is often estimated at half of passenger wait time. Thus, we setcpt=0.5·cpw and obtain the results that correspond to different values of cpw, as shown in Figure 15. Results show that as time value increases, total passenger wait time tends to decrease, while the frequency within the peak hour tends to increase. They demonstrate the tradeoffs between passenger-oriented objectives (service level) and operator-oriented objectives (train cost). It may be expected that if the average value of time continues growing, the system manager should add more trains to improve passenger satisfaction.Figure 15
Sensitivity to unit passenger time cost.
(a)
Passenger time with different unit time costs (b)
Peak hour frequency with different unit time costsIn addition, neither passenger wait time nor peak hour frequency changes are strictly monotonic. This occurs because the minimum headway constraints would somehow limit the flexibility of the RO method, especially based on a major metro line (as described above), where very different distributions of passenger demand lead to varying dwell times. Thus, in the real-world operations, safety issues may limit improvements in service level. The timetables respond well to different values of time.Besides, we find that total passenger travel time barely changes. In fact, all the segment travel times reach minimum values. The reason will be specifically discussed in the second experiment.The second experiment tests how unit energy cost affects passenger time and train round trip time, based on the passenger demand in case (c-1) of BML4. Figure16(a) demonstrates that when unit energy cost changes within a practical range, neither passenger wait time nor travel time is affected. The minimum segment travel times are obtained in all cases. This can be seen in daily operations that, with large enough passenger demand, the metro operators would run trains at the fastest speed, that is, with minimum segment travel times.Figure 16
Sensitivity to unit energy cost.
(a)
Passenger time with different unit energy costs (b)
Round trip time with different unit time costsTo extend the discussion, we run more experiments to show the relations between larger unit energy costs and round trip time and obtain the results shown in Figure16(b). Although such large unit energy costs are unrealistic, we can see that the round trip time increases; that is, the average speed decreases, as energy cost increases. In other words, as the ratio between unit energy cost and passenger demand increases, round trip time increases. This result is consistent with the practical operation strategy that trains should run faster with more passengers and slower when they carry fewer passengers (where the unit energy cost is fixed, and thus the ratio depends on the passenger demand).The third experiment aims to test the sensitivity of the passenger wait time and the total cost to train capital cost, based on the passenger demand in case (c-1). Figure17 illustrates that, with the increased train capital cost, the passenger wait time and total cost tend to increase. The tradeoffs between passenger wait time and train cost, that is, the passenger-oriented and operator-oriented objectives, is again demonstrated. The fluctuation of passenger wait time is caused by the minimum headway constraints, as explained for the first experiment. In addition, the monotonically increasing total cost demonstrates that technologies which increase the train’s economic life help reduce total cost in a metro system.Figure 17
Sensitivity to train capital cost.Generally, we obtain reasonable and practical results with the proposed solution approach. These results are consistent with practical experiences and can be used as decision-making support in daily operations. Since similar results are obtained with different passenger demand patterns and on different metro lines, the predictive ability of the model is demonstrated. The effectiveness of the proposed approach is further verified.


## 5.1. A Small Case Study

Here we consider a short metro line, as shown in Figure5. The segment lengths are 1800 m, 1600 m, and 2000 m. For simplicity, the minimum and maximum dwell times for each station are all set as 30 s and 90 s, respectively.Figure 5
Representation of a short metro line.The baseline passenger demand is shown in Table4. The “Boarding (Alighting) Sum” represents the number of total boarding (alighting) passengers per second at a station. Four cases are considered here. As shown in Figure 6, the curves represent demand variance over time, and “multiplier” represents the ratio between actual and baseline passenger demand. For example, in case 1, during period [1200,2400] (s), the multiplier is 1.6, and passenger demand from station 2 to station 3 at t∈[1200,2400] is 0.96 (0.6×1.6) passengers/second. The four cases have different passenger demand distributions.Table 4
Baseline passenger demands for a short metro line (passengers per second).
Station index 1 2 3 4 5 6 7 8 Boarding sum 1 0 0.3 0.3 0.45 0 0 0 0 1.05 2 0 0 0.6 0.3 0 0 0 0 0.9 3 0 0 0 0.3 0 0 0 0 0.3 4 0 0 0 0 0 0 0 0 0 5 0 0 0 0 0 0.45 0.75 0.3 1.5 6 0 0 0 0 0 0 0.6 0.3 0.9 7 0 0 0 0 0 0 0 0.45 0.45 8 0 0 0 0 0 0 0 0 0 Alighting sum 0 0.3 0.9 1.05 0 0.45 1.35 1.05 -Figure 6
Passenger demand distributions in different cases.The computational results of different cases are shown in Table5 and Figure 7. Specifically, in Table 5, the noncyclic result obtained by RO is denoted as NK=KN, which schedules KN trains within a given time duration. The cyclic result obtained by Algorithm 1 is denoted as UK=KU, which schedules KU trains within a given time duration. Particularly, the best cyclic result is denoted as UK=KUbest. The GGA result is denoted as GAK=KN, which uses frequency KN from RO result as input. The objectives of different results are calculated, including passenger wait time (Tpw), passenger travel time (Tpt), total energy (Etotal), fleet size (Fs), and total cost (Ctotal). Figure 7 specifically shows the numerical results in passenger satisfaction.Table 5
Comparison between cyclic and non-cyclic timetables.
Time duration (s) Timetable type  T p w  (106) T p t  (106) E t o t a l   (103) F s C t o t a l   (104) Gap (%) 0-10800 Non-cyclic N K = 60 7.2330 14.384 18.069 7 12.320 - Case 1 Non-cyclic G A K = 60 7.2786 14.402 18.069 7 12.329 0.07 Cyclic U K = 53 b e s t 8.3360 14.523 16.148 6 12.496 1.41 Cyclic U K = 60 7.3561 14.406 18.069 7 + ∞ -  0-4800 Non-cyclic N K = 24 2.3171 4.9470 7.1411 6 4.3586 - Case 2 Non-cyclic G A K = 24 2.5489 4.9758 7.1411 6 4.4954 3.04 Cyclic U K = 24 b e s t 2.6102 4.9988 7.1411 6 4.5407 4.01 Cyclic U K = 23 2.7619 5.0347 6.8667 6 4.5681 4.59 Cyclic U K = 25 2.4989 4.9859 7.4156 7 4.5441 4.08  0-7200 Non-cyclic N K = 45 4.3948 11.292 13.615 8 8.8022 - Case 3 Non-cyclic G A K = 45 4.5726 11.328 13.615 8 8.9245 1.37 Cyclic U K = 48 b e s t 4.3975 11.319 14.438 8 9.0243 2.46 Cyclic U K = 45 4.7268 11.417 13.615 8 9.0345 2.57  0-10800 Non-cyclic N K = 60 6.0280 14.301 18.069 7 11.621 - Case 4 Non-cyclic G A K = 60 6.3113 14.560 18.069 7 11.854 1.97 Cyclic U K = 60 b e s t 6.4638 14.330 18.069 7 11.875 2.14 Cyclic U K = 59 6.5981 14.381 17.795 7 11.899 2.34 Cyclic U K = 61 6.3572 14.327 18.344 7 11.884 2.21Figure 7
Passenger costs for different results in different cases.Table5 and Figure 7 demonstrate that, under time-varying passenger demand, noncyclic timetables obtained by RO outperform the other timetables, while GGA results outperform cyclic ones obtained by Algorithm 1. This result verifies that the noncyclic timetable improves passenger satisfaction. Particularly, the passenger wait time (Tpw), passenger travel time (Tpt), and total cost (Ctotal) have their lowest values in the RO results, which demonstrates that the RO method is effective in solving the proposed train timetabling problem. In case 1, the result of UK=60 is infeasible because it violates the minimum headway constraint at intermediate stations (the dwell times are different from UK=60 to NK=60), and thus we have Ctotal=+∞. It is interesting that NK=60 and UK=60 require the same energy and fleet size. This result verifies the reasonableness of model decomposition in this paper from two perspectives: (1) When all passengers reach their destinations, energy is solely determined by segment travel times (they are the same in NK=60 and UK=60). (2) With given segment travel times, the estimated fleet size is solely determined by headways. Similar results can be observed in other cases.The gaps in Table5 shows that case 2 has the largest gap, case 3 and 4 have similar gaps, and case 1 has the smallest gap. Note that, in case 2, passenger demand changes drastically; in cases 3 and 4, demand patterns are quite similar; in case 1, within certain subperiods, passenger demand is the same, which is the steadiest among all cases. Therefore, we can conclude that as the passenger demand varies more significantly over time, the gap between the cyclic and noncyclic timetables increases.Figure8 illustrates the obtained train timetables in different cases, where the departure headways vary in accordance with the fluctuating passenger demand. For example, in cases 2, 3, and 4, the headways are small in the peak period and are larger before and after the peak period. In case 1, the small headways occur during subperiods [1200,2400] (s) and [7200,8400] (s), where the multipliers of passenger demand are largest (1.6 and 1.8). In Figure 8, the circles represent the dwell times that exceed the minimum dwell time (i.e., 30 s). We observe that the excess dwell times (circles) mostly occur at station 7. This is consistent with the OD pairs we present in Table 4. Specifically, the difference between the cumulative boarding sum and the cumulative alighting sum reaches its maximum value at station 7 (1.5+0.9-0.45=1.95). In this regard, the alighting and boarding rates are smallest at station 7. In addition, the total number of alighting and boarding passengers is largest (i.e., 1.8 passengers per second) at station 7, and thus the excess dwell times here are reasonable.Figure 8
Obtained timetables for different cases.
(a)
Obtained timetable for case 1 (NK=60) (b)
Obtained timetable for case 2 (NK=24) (c)
Obtained timetable for case 3 (NK=45) (d)
Obtained timetable for case 4 (NK=60)Moreover, despite the symmetry of passenger demand distributions, the headways at the ends of periods are larger than at their starts. This can be explained by that passenger demand increases at the start and decreases to zero at the end, which leads to different cost changes. Specifically, the passenger time cost increases more significantly when demand increases, while the train cost does not change with demand. Thus, the headway tends to be smaller at the start of periods than at the ends. In fact, it is often seen in actual operations that the headways between the last trains at the end of day exceed those between the first trains in the morning.Finally, Table6 shows the computational performances for different methods, including RO method, GGA, and AL + PS method. As described above, GGA and AL + PS require a preset fixed frequency, and thus we use the same frequency from RO result KN as input. We terminate the algorithm when its computation time reaches 3 hours and output the best solution obtained by then.Table 6
Computational performances for different methods.
RO method Generic GA AL + PS (Case 1) Time (s) 225.5 8299.4 10800 Total cost (¥) 12.32 × 10 4 12.33 × 10 4 12.42 × 10 4  (Case 2) Time (s) 86.1 4658.2 10800 Total cost (¥) 4.3586 × 10 4 4.4954 × 10 4 4.4872 × 10 4  (Case 3) Time (s) 120.6 5009.3 10800 Total cost (¥) 8.8022 × 10 4 8.9245 × 10 4 8.9567 × 10 4  (Case 4) Time (s) 240.717 8229.1 10800 Total cost (¥) 11.621 × 10 4 11.854 × 10 4 11.774 × 10 4Table6 demonstrates that RO is very efficient, obtaining a timetable very fast. In comparison, GGA takes more than 1.5 hours to obtain a timetable, and AL + PS takes more than 3 hours in all cases. In each case, the timetable obtained with the RO method has the lowest total cost, which further verifies its effectiveness. The GGA timetables perform better than AL + PS in cases 1 and 3 but are outperformed in cases 2 and 4.Combining all the computational results in this case study, we can conclude some disadvantages of GGA and AL + PS compared to RO.(1) Although the timetables obtained by them are noncyclic, train departure headways in their results differ greatly from one another without clear correspondence to the passenger demand. In comparison, RO timetables are more regular (headways change gradually) and capture well the passenger demand patterns. (2) Due to the computation complexity, the computation time for running GGA and AL + PS once is far greater than RO. (3) They both require a preset frequency as input. If we do not know the best frequency beforehand, which is true in most real-world scenarios, we must run these methods with different frequencies and select the best solution, thereby further increasing substantially computation time.In general, the RO method is superior to GGA and AL + PS in obtaining a noncyclic timetable with time-varying passenger demand.


## 5.2. Case Study of Beijing Metro Line 4

This section aims to test the applicability and effectiveness of RO in a real-world case study based on Beijing’s Metro Line 4 (BML4). The infrastructure information and some operational parameters of BML4 are shown in Table7. The turnaround time is 120s. To satisfy the large passenger demand, we set train capacity as Nk=2160 passengers, which is a maximum overloaded capacity. If not mentioned, the other parameters are the same as in Table 3. The baseline passenger demand data are based on the research conducted by Duan et al. [53], which are presented in Figure 9, where zx,y, ∀x<y represents morning peak demand and zx,y, ∀x>y represents evening peak demand.Table 7
Real-world infrastructure data of Beijing Metro Line 4.
Station Minimum dwell time (s) Maximum dwell time (s) Segment Minimum travel time (s) Maximum travel time (s) Segment length (m) 1/48 AHQB 30 90 1 124 170 1363 2/47 BGM 30 90 2 114 156 1251 3/46 XY 30 90 3 152 209 1672 4/45 YMY 30 90 4 118 162 1295 5/44 BDDM 30 90 5 81 111 887 6/43 ZGC 30 90 6 82 113 900 7/42 HDHZ 30 90 7 97 133 1063 8/41 RMDX 30 90 8 96 131 1051 9/40 WGC 30 90 9 151 207 1658 10/39 NL 30 90 10 138 190 1517 11/38 BJZ 30 90 11 131 180 1441 12/37 XZM 30 90 12 93 128 1025 13/36 XJK 30 90 13 100 138 1100 14/35 PAL 30 90 14 100 138 1100 15/34 XS 30 90 15 79 109 869 16/33 LJHT 30 90 16 92 126 1011 17/32 XD 30 90 17 74 102 815 18/31 XWM 30 90 18 105 144 1152 19/30 CSK 30 90 19 109 150 1200 20/29 TRT 30 90 20 149 205 1643 21/28 BJSRS 30 90 21 135 185 1480 22/27 MJP 30 90 22 75 103 827 23/26 JMX 30 90 23 90 124 989 24/25 GYXQ 30 90 - - - -Figure 9
Baseline passenger demand for BML4.Here we consider three types of scenarios whose passenger demands are large and unsteady, which are difficult to manage: (a) period of 3 hours (i.e., 10800 s) including morning peak; (b) period of 3 hours including evening peak; and (c) a shopping festival hosted in the shopping centers next to the XD station (whose index is 8/41). Representations of passenger demands are the same as in the small case study (as shown in Figure10): the actual passenger demand is the product of baseline passenger demand (Figure 9) and the multipliers (Figure 10). To test RO with more instances, we design three subcases of scenario (a) with different peak multipliers, and two subcases of scenario (c): without the special event (c-1) or with it (c-2). In case (c-2), the passenger demand at station XD is tripled.Figure 10
Passenger demand patterns in different cases.It should be noted that the real operation data of BML4 is dynamic and confidential. Besides, the operator today prefers to schedule trains uniformly. Therefore, in this experiment, we use cyclic timetables with best performance to approximate the real-world timetables. The notations are the same as in the small case study.Table8 and Figure 11 show the computational results in different cases. Note that, for a major metro line such as BML4, the computation times of GGA and AL + PS are extremely large, and in many cases they do not yield a satisfactory solution within a reasonable time. Thus, we do not present their results in this experiment. The results of UK=46 in case (a-1), UK=54 in case (a-2), and UK=54 in case (a-3) are infeasible, and thus their total costs are +∞. Figure 11 shows that considerable costs are reduced by the noncyclic timetables obtained by RO. Interestingly, in the cases where cyclic timetables are feasible, that is, cases (a-1), (b), (c-1), and (c-2), the frequency of the noncyclic timetable is always larger than the best cyclic ones. This is because the passenger demand in the real-world fluctuates so drastically that cyclic timetables cannot capture its dynamic patterns, and the minimum headway constraint makes it hard to obtain feasible cyclic timetables. In this regard, adding more trains in a cyclic timetable may be either uneconomic (i.e., the increased train operating cost exceeds the reduced passenger time value) or infeasible.Table 8
Computational results for BML4.
Scenario   T p w  (108) T p t  (108) E t o t a l  (104) F s C t o t a l  (106) Gap (%) Case (a-1) Non-cyclic N K = 46 3.8665 3.0555 7.4898 31 3.1788 - 3869 s Cyclic U K = 44 b e s t 5.3237 3.0943 7.2385 30 3.9925 20.38 Cyclic U K = 46 4.9805 3.0985 7.4898 31 + ∞ -  Case (a-2) Non-cyclic N K = 54 11.426 4.3898 9.2418 37 7.7869 - 1544 s Cyclic U K = 54 10.653 4.4060 9.2418 37 + ∞ -  Case (a-3) Non-cyclic N K = 54 19.568 4.7953 9.4566 38 12.426 - 1133 s Cyclic U K = 54 18.891 4.816 9.4958 38 + ∞ -  Case (b) Non-cyclic N K = 38 2.2416 2.1501 5.9624 25 1.9902 - 3932 s Cyclic U K = 30 b e s t 4.4170 2.1632 4.9379 21 3.1748 37.31  Case (c-1) Non-cyclic N K = 30 0.5482 1.0256 4.3421 19 0.7004 - 6848 s Cyclic U K = 22 b e s t 1.2744 1.0946 3.3695 14 1.0957 36.08  Case (c-2) Non-cyclic N K = 41 1.1415 1.5061 6.0564 26 1.2342 - 4068 s Cyclic U K = 25 b e s t 3.8100 1.5533 3.9911 17 2.6459 53.53Figure 11
Total costs for different results in different cases.In comparison, the noncyclic timetables obtained by RO respond well to demand changes. For example, in case (a-1), the result ofNK=46 substantially reduces passenger wait time compared to UK=44best and UK=46. The adaptation for time-varying demand is further demonstrated in cases (c-1) and (c-2). In case (c-1), the demand distribution pattern is similar to cases (a-1) and (b), and the cost reduction of the noncyclic timetable is closer to the other cases (36.08%). In case (c-2), the special event leads to a more heterogeneous demand distribution, and the cost reduction reaches 53.35%. Besides, the computation times (shown below each case) are acceptable for such a major and busy metro line.However, in cases (a-2) and (a-3), the passenger wait times for noncyclic timetables are larger than for cyclic ones. This occurs because the RO method considers the objective function that corresponds to an infeasible headway as+∞. With a large demand, the first departure headway that is feasible may greatly exceed h∗ (as can be seen in Figure 4), whereas the cyclic timetable does not consider the feasibility and divides headways uniformly. In other words, the increase of departure headways (leading to larger Tpw) is due to safety issues, which is consistent with the actual operational requirements that put safety ahead of costs. In addition, with parameters given in this experiment, the largest frequency for BML4 is KN=54, as can be concluded by comparing cases (a-2) and (a-3). Therefore, real-time dispatching and passenger flow control methods should be applied with passenger demand at or above the levels in case (a-2).Figures12, 13, and 14 illustrate the noncyclic timetables obtained in case (a-1), case (b), and case (c-2), respectively. The circles here represent the occurrences of maximum dwell times. These figures illustrate that the headways vary in accordance with demand patterns, which is consistent with the results in the small case study. Particularly in Figure 14, the dense headways in the middle are set to transport many passengers due to the special event. This shows the flexibility of noncyclic timetables in special scenarios.Figure 12
Timetable for morning peak (case a-1).Figure 13
Timetable for evening peak (case b).Figure 14
Timetable for a special event (case c-2).In addition, from the distribution of the circles, we can see that case (a-1) is the busiest scenario among these figures. The many circles in case (a-1) help explain the results in cases (a-2) and (a-3) (which have larger demands than case (a-1)). With so many trains taking maximum dwell times to load passengers, real-time dispatching and passenger flow control methods are necessary. In case (b), the maximum dwell times mostly occur at stations 7, 11, and 12. Specifically, the neighborhood of station 7 (HDHZ) includes famous schools, fashion shopping streets, and dense office buildings. The neighborhood of station 11 (BJZ) includes a major bus hub, the Beijing Zoo, and one of the largest wholesale markets. The neighborhood of station 12 (XZM), which has transfers with two other busy metro lines, includes a railway station, large shopping malls, and office buildings. They are well known for the large numbers of passengers during peak hours. Thus, we find that the model is reasonable and consistent with real-world scenarios. Analogously, in case (c-2), the circles are mostly seen at stations near the station XD, where the special event is held.It should be noted that the last two trains have an unusually small headway in case (b). The reason is that we assume that the last train must arrive at the first station at the end of the period. It is expected that if there is still demand after 10800 s, the headways will be optimized by our approaches.


## 5.3. Additional Experiments

In this section, we implement three sets of additional experiments to investigate the relations between the values of costs and the objectives. The first experiment tests the sensitivity of passenger time and peak frequency to unit passenger time cost, based on the demand in case (b) of BML4 study. In practice, the unit time value of passenger travel time is often estimated at half of passenger wait time. Thus, we setcpt=0.5·cpw and obtain the results that correspond to different values of cpw, as shown in Figure 15. Results show that as time value increases, total passenger wait time tends to decrease, while the frequency within the peak hour tends to increase. They demonstrate the tradeoffs between passenger-oriented objectives (service level) and operator-oriented objectives (train cost). It may be expected that if the average value of time continues growing, the system manager should add more trains to improve passenger satisfaction.Figure 15
Sensitivity to unit passenger time cost.
(a)
Passenger time with different unit time costs (b)
Peak hour frequency with different unit time costsIn addition, neither passenger wait time nor peak hour frequency changes are strictly monotonic. This occurs because the minimum headway constraints would somehow limit the flexibility of the RO method, especially based on a major metro line (as described above), where very different distributions of passenger demand lead to varying dwell times. Thus, in the real-world operations, safety issues may limit improvements in service level. The timetables respond well to different values of time.Besides, we find that total passenger travel time barely changes. In fact, all the segment travel times reach minimum values. The reason will be specifically discussed in the second experiment.The second experiment tests how unit energy cost affects passenger time and train round trip time, based on the passenger demand in case (c-1) of BML4. Figure16(a) demonstrates that when unit energy cost changes within a practical range, neither passenger wait time nor travel time is affected. The minimum segment travel times are obtained in all cases. This can be seen in daily operations that, with large enough passenger demand, the metro operators would run trains at the fastest speed, that is, with minimum segment travel times.Figure 16
Sensitivity to unit energy cost.
(a)
Passenger time with different unit energy costs (b)
Round trip time with different unit time costsTo extend the discussion, we run more experiments to show the relations between larger unit energy costs and round trip time and obtain the results shown in Figure16(b). Although such large unit energy costs are unrealistic, we can see that the round trip time increases; that is, the average speed decreases, as energy cost increases. In other words, as the ratio between unit energy cost and passenger demand increases, round trip time increases. This result is consistent with the practical operation strategy that trains should run faster with more passengers and slower when they carry fewer passengers (where the unit energy cost is fixed, and thus the ratio depends on the passenger demand).The third experiment aims to test the sensitivity of the passenger wait time and the total cost to train capital cost, based on the passenger demand in case (c-1). Figure17 illustrates that, with the increased train capital cost, the passenger wait time and total cost tend to increase. The tradeoffs between passenger wait time and train cost, that is, the passenger-oriented and operator-oriented objectives, is again demonstrated. The fluctuation of passenger wait time is caused by the minimum headway constraints, as explained for the first experiment. In addition, the monotonically increasing total cost demonstrates that technologies which increase the train’s economic life help reduce total cost in a metro system.Figure 17
Sensitivity to train capital cost.Generally, we obtain reasonable and practical results with the proposed solution approach. These results are consistent with practical experiences and can be used as decision-making support in daily operations. Since similar results are obtained with different passenger demand patterns and on different metro lines, the predictive ability of the model is demonstrated. The effectiveness of the proposed approach is further verified.


## 6. Conclusions

This paper investigates the train timetabling problem (TTP) with time-varying passenger demand and accomplishes the following:(1) It accounts for both the passenger-oriented and operator-oriented objectives. Factors in the line planning stage, the timetabling stage, and the vehicle scheduling stage are integrated in our model and solution approach. The noncyclic timetable obtained in this paper outperforms the cyclic timetable, the noncyclic timetables obtained by the generic GA, and a classic direct optimization approach.(2) The dwell time determining process is specifically designed to consider both congestion at stations and passenger demand. Besides, limited train capacity and varying passenger alighting/boarding rates are accounted for. Thus, the noncyclic timetable obtained in this paper is more practical in depicting passenger loading evolution. For this reason, some timetables with otherwise good properties would violate the minimum headway constraint at intermediate stations and thus be infeasible. Passenger flow control methods and real-time dispatching can be employed to make these timetables applicable.(3) Based on a general scheduling scheme, the TTP is decomposed into the OST and OH subproblems. The solution approach includes a branch-and-bound (B&B) algorithm, a frequency determining algorithm, and a novel rolling optimization (RO) method. Then, computational results based on a short metro line and the Beijing Metro Line 4 are obtained. They verify the effectiveness and efficiency of the proposed solution approach and demonstrate that the noncyclic timetables obtained in this paper respond well to passenger demand characteristics. Since the results are based on different cases with different passenger distributions, the predictive ability of the proposed solution approach can also be demonstrated.(4) Additional experiments are presented whose results demonstrate the following points: (a) As passenger time values increase, passenger wait time may be further reduced by scheduling more trains or by shifting departure times for riders. (b) Average train speed decreases as passenger demand increases, due to longer dwell times at stations. This dwell time effect outweighs the speed increase justified by more riders per train. (c) Longer train economic life not only reduces total system cost, but also indirectly improves the service level.This paper has some limitations. Specifically, the linear piecewise approximation simplifies the energy calculation. Assumption1 is consistent with a common strategy in practice but does not account for the relevant effects of driver behaviors. Assumption 4 may be extended to a more accurate model which distributes the passenger boarding and alighting over the station dwell times. A more systematic comparison with real data of BML4 in the case study should be conducted to fully validate RO method.These drawbacks can be addressed by a more complex model which considers driver behaviors (by estimating the magnitude of their relevance), line layouts, and physics of real train run process, as well as simultaneously optimizing departure headways and segment travel times under continuous time. Future studies can investigate the improvement of our solution approach, to reduce the computation time needed for such a complex model. Besides, it would also be interesting to adjust the model’s parameters with respect to different real-world operation data, to obtain more realistic timetables.



---

*Source: 3690603-2018-07-03.xml*