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.
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. Section 2 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
Table 1 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.
Apart 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
Table 2 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.
Assumption 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 should not exceed total passenger boarding rate . If , the queue of waiting passengers increases even if a train is dwelling at the station. In practice, passenger entrance flow is taken to avoid such a case.
. 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 train who have reached their destination are assumed to exit the train at its arrival time of station , that is, . All waiting passengers at station able to board train are assumed to board the train at its departure time, that is, . In other words, the number of in-vehicle passengers changes only at the arrival/departure times of train at stations.
This paper focuses on a bidirectional metro line with a single track, as shown in Figure 1. 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.
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 corresponds to the set of time stamps in which trains can arrive at (depart from) station . The nodes in 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.
The continuous time is discretized into evenly spaced time intervals (1 second in this paper), and the period extends from 0 to . Note that the values of and are different. is the last time stamp of the passenger demand. The period mentioned above refers to the one extending from 0 to . Assumption 2 imposes that the last train departs from the first station at time . However, if the DTSG only considers the period from 0 to , 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 . The relation between and is set as:where is the segment travel time of segment and is the maximum dwell time of station . 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 , which is defined as
Passenger arrival rates at stations are , . In this paper, demands are regarded as deterministic values, which are given by . Thus, the number of arriving passengers at station with destination during time interval is
The number of arriving passengers at station within time interval is the sum of arriving passengers with different destinations:
(2) Determination of Dwell Arcs. Letting be the number of doors per train, the average number of alighting (boarding) passengers per door when train arrives at (departs from) station is . In addition, we denote the average number of boarding passengers per second (boarding rate) as and alighting rate as , 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 .
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:
The calculation of boarding time is slightly more complex. Figure 3 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 at station , and in Figure 3 they are represented by remaining capacity of train at station : , a horizontal line; number of waiting passengers at station at time , a curve ; maximum allowable number of boarding passengers, a diagonal function with the slope of ; minimum and maximum dwell times at station , that is, and . In Figure 3, the diagonal intersects the horizontal line at time and intersects the curve at time . The curve intersects the horizontal line at time . Note that the time is discrete, and Assumption 3 ensures that the slope of always exceeds the slope of . Thus, we have
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 exceed train capacity, and they will be left behind; (7) shows that before there are passengers waiting at station, and after all passengers board the train; (8) points that the train is full at , and passengers arriving after should wait for the next train.
The departure time of train at station is determined based on the relations among , , , , and . Based on Assumption 3, we consider two cases: and . Note that and should be considered in both cases. In case , we can further derive that (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 . In case , we can further derive that (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 . The departure time of train at station is expressed as
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 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 during the dwell process is fixed and is determined by the number of in-vehicle passengers at its arrival time. That is, and remain constant during the dwell process:where and 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 at station are determined. Based on the results of dwell arcs, three cases should be considered: ; and ; and . In case , waiting passengers board the train at rate . In case , actual passenger boarding rate is equal to the passenger arrival rate. In case , the actual boarding rate is . Thus, the number of boarding passengers for train at station is
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 and of total boarding passengers is equivalent to the ratio between and :
The number of alighting passengers from train at station consists of all previously boarding passengers in train with destination :
There are no alighting passengers at the first station in each direction (station 1 and ).
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 be the index of last train that has departed from station before time , that is, . The cumulative number of passengers at station at time is the difference between cumulative numbers of arrival passengers and boarding passengers within the period :
Since the continuous time is discretized into 1 second intervals, the passenger wait time between times and is . Thus, the total passenger wait time is
(2) Passenger Travel Time. The number of in-train passengers in train at time is the difference between cumulative numbers of boarding passengers and alighting passengers. Let be the index of the last station that train has visited before time . We consider two cases: train is dwelling at station , that is, ; train is in the segment between station and , that is, . Based on Assumption 4, the boarding passengers at station are not considered in case but are included in case . Thus, the number of in-vehicle passengers in train at time is
Analogously, the total passenger travel time is
Finally, the total cost of passenger time is formulated as
The cost unit is yuan .
3.3. Energy
In general, for each segment , with given segment travel time and segment length , the speed profile determines energy use [49]. Thus, given a speed profile generating method, the energy per unit mass of different trains in segment is the same. We denote the energy per unit mass for one train in segment as . Its calculation usually considers the following train dynamic equations [48]:
The first two equations are train motion equations, where and , 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
However, in (20), the relation between energy and decision variable 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 in segment [48]:
We introduce here an estimation procedure for the values of and , which is based on mean squared error (MSE). First, we set several feasible values (e.g., ) of and . Then, we obtain optimal speed profiles based on (19) and (20), thereby obtaining the optimal energy with respect to and . Finally, we define the optimal energy per unit mass as , and the estimated energy per unit mass obtained with (21) as . The MSE is formulated as
The optimal values of and can be obtained when reaches its minimum and are found as follows:
With and obtained, the total energy can be calculated as
Thus, the total cost of energy consumption is
Note that when all tracks are straight and level, and suffice for estimating energy. However, if grades and curves are considered, we may have , where and are specified parameters to describe the characteristics in segment . 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. Let be the annual cost per train, be the current price per train, be the interest rate, and be the lifespan of a train. The relations among these four factors can be expressed as
The capital cost per train hour can be then calculated by dividing the annual cost by the product of average working hours per day and working days per year
In addition, a train can be reassigned after it finishes a round trip and returns to the shunting yard, and the round trip time can be obtained as
Thus, the actual period for train running is . The estimated fleet size within period [] depends on the round trip time and the departure headway from the shunting yard (which is equivalent to the arrival time of train at the first station ):where the term represents the average round trip time and the term represents the average headway. The total cost of trains is formulated as
The first term represents the capital cost of trains for the actual running period . The second term represents the operating cost with frequency (trains) on a line with the length of the sum of 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, , , , , and . For example, an operator-oriented timetable (neglecting passenger satisfaction) can be obtained by setting and ; a passenger-oriented timetable can be obtained by setting . The weights for integrated objectives can be set similarly, with larger cost coefficients for larger weights. The model is formulated as follows:
Minimize
Subject to
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 to as traveling in segment , with fixed travel time . Constraint (34) ensures that the segment travel time should be optimized within a reasonable range . 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 matrixes and , respectively, representing the selection of arrival and departure nodes in the DTSG:where if a train arrives at (departs from) station at time . As mentioned earlier, the value of is determined by , for , . The value of is determined by and the dwell arcs, that is, , for , , , . The value of is determined by and the travel time in segment , that is, , for , , and , . Thus, the decision variables determine and .
Based on the structures of and , as well as Assumption 4, we introduce two matrixes and , respectively, representing passenger alighting and boarding information (number and time):
The cumulative numbers of waiting passengers and in-vehicle passengers can be constructed based on and :
In (39), represents the total waiting passengers at station , where is the number of waiting passengers at the start of this period. According to the definition of cumulative passengers, we have . Thus, the total passenger wait time can be expressed by since . In (40), the initial number of in-vehicle passengers is zero, that is, . 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 . Thus, the total cost for passenger time can be rewritten as
For energy cost, the energy per unit mass in segment can be rewritten as
The total energy can be rewritten aswhere if ; otherwise . The total energy cost can be rewritten correspondingly.
For train cost, the round trip time can be rewritten as
The departure headway from the shunting yard can be rewritten as
The total cost of trains can be rewritten correspondingly.
Finally, the last departing train constraint (32) is equivalent to
The turnaround constraint (33) is equivalent to
The segment travel time window constraint (34) is equivalent to
The minimum headway constraint (35) is equivalent to
The train capacity constraint (36) is equivalent to
The selection of arcs determines the arrival nodes and departure nodes , 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 . 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 is determined by passenger demand. In this regard, the total energy can be expressed as
Therefore, the total energy cost depends only on segment travel times . 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 Assumption 1, 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
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 matrix , 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 and ) instead of costs of passenger wait and dwell times. The OH subproblem is formulated as follows (with segment travel times as input):
Minimize
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 with decision variables (, for ). 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:where is the coefficient vector and is the decision variable vector. Both and have elements. and are coefficient matrixes, and and are constant vectors representing the bounds. Here is a matrix, and is a matrix. Correspondingly, has elements, and has one element. We define a vector with elements as
The decision variable vector can be represented as for , where is the -th element in the vector . We then define for as the total number of passengers that travel through segment . Based on (52), can be represented as
Constraints (55) are represented as follows. Based on constraint (47), we have and . Based on constraint (48), we have , for ; , for , where represents the -th row of matrix . Analogously, , for ; , for . 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: The cyclic timetable with minimum cost (i.e., the best cyclic timetable) can be considered as an upper bound of the noncyclic timetable. If passenger demand is stable, the best cyclic timetable is a good approximation to an optimal noncyclic one. The cyclic timetables are easy to obtain.
First, the range of frequency should be estimated. We define the peak passenger volume in one segment as , which is the maximum value among the numbers of passengers traveling through the segment, that is, . The minimum frequency is determined by the train capacity and . The maximum frequency is determined by the minimum headway: .
Second, for each frequency , we use 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 (, where is the departure headway), the cyclic timetable that corresponds to frequency 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 frequency . 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.
| ||||||||||||||||||||
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 , e.g., 600 s, for solving the OST. For example, for each train , its previous train arrives at the first station at time (the first train arrives at time ), and we consider the period to solve the OST (the constraint of ensures that ). 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 train if its arrival time at the first station is given. Since each train is scheduled in chronological order, the waiting passenger matrix and the in-train passenger matrix can be constructed chronologically. For each train , the departure time of its preceding train at the first station is , and train departs from the first station at time . We define the passenger wait time with respect to train as , which is formulated as
Analogously, we define the passenger travel time for train as , and we obtain
The estimated fleet size of train () depends on its round trip time () and the departure headway from the shunting yard ():
Based on (18), (30), (53), (58), (59), and (60), we define the total costs in the OH for train as , which can be expressed as
The value of is a function of the arrival time for train at the first station. Since the arrival time for the previous train is given, we consider the value of as a function of the headway between train and train (), for consistency with other TTP studies.
Next we discuss how the optimal value for 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 and cannot be described with a fixed function.
The expected objective function of OH is shown in Figure 4. Its value first decreases (the headway starts from ) to a minimum value and then increases monotonically. We specify that this function reaches its minimum value at headway . Since the DTSG divided the continuous time into seconds, we iteratively calculate the objective function of OH with respect to (starting from , and in each following step), until the optimal headway is found, or the total cost increases as increases (in case that 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 to be if it corresponds to an infeasible headway.
Finally, for each train , we can obtain its optimal headway (between train and train ). By rolling optimization of departure headways for trains, the timetable can be obtained. Let be the vector of optimal headways for all trains . 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 .
| ||||||||||||||||||
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 ( 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: 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. In terms of optimization process, GGA optimizes the solution (of all trains) in a random direction, while RO optimizes one train at a time. 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 to , where is a new variable). Let the number of constraints be . Set the multiplier that corresponds to the -th constraint as , and set a large enough parameter . By minimizing the objective function with respect to (so that the intermediate variable does not influence the objective function), we can obtain the typical form of the penalty term as
The parameter value 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 Table 3. 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 and in Table 3.
5.1. A Small Case Study
Here we consider a short metro line, as shown in Figure 5. 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.
The baseline passenger demand is shown in Table 4. 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 (s), the multiplier is 1.6, and passenger demand from station 2 to station 3 at is 0.96 () passengers/second. The four cases have different passenger demand distributions.
The computational results of different cases are shown in Table 5 and Figure 7. Specifically, in Table 5, the noncyclic result obtained by RO is denoted as , which schedules trains within a given time duration. The cyclic result obtained by Algorithm 1 is denoted as , which schedules trains within a given time duration. Particularly, the best cyclic result is denoted as . The GGA result is denoted as , which uses frequency from RO result as input. The objectives of different results are calculated, including passenger wait time (), passenger travel time (), total energy (), fleet size (), and total cost (). Figure 7 specifically shows the numerical results in passenger satisfaction.
Table 5 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 (), passenger travel time (), and total cost () 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 is infeasible because it violates the minimum headway constraint at intermediate stations (the dwell times are different from to ), and thus we have . It is interesting that and require the same energy and fleet size. This result verifies the reasonableness of model decomposition in this paper from two perspectives: When all passengers reach their destinations, energy is solely determined by segment travel times (they are the same in and ). 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 Table 5 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.
Figure 8 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, , and , 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 (s) and (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 (). 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.
(a) Obtained timetable for case 1 ()
(b) Obtained timetable for case 2 ()
(c) Obtained timetable for case 3 ()
(d) Obtained timetable for case 4 ()