Single-line rail rapid transit timetabling under dynamic passenger demand

Transportation Research Part B 70 (2014) 134–150

Contents lists available at ScienceDirect

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

Single-line rail rapid transit timetabling under dynamic passenger demand Eva Barrena a,b,⇑, David Canca c, Leandro C. Coelho a,d, Gilbert Laporte a,b a

Interuniversity Research Center on Network Enterprise, Logistics and Transportation (CIRRELT), Canada HEC Montréal, 3000 chemin de la Côte-Sainte-Catherine, Montréal H3T 2A7, Canada c School of Engineering, University of Seville, Avenida de los Descubrimientos s/n, 41092 Seville, Spain d Faculté des sciences de l’administration, Université Laval, 2325 rue de la Terrasse, Québec G1K 7P4, Canada b

a r t i c l e

i n f o

Article history: Received 9 December 2013 Received in revised form 25 August 2014 Accepted 27 August 2014

Keywords: Rail rapid transit Demand-based timetable Passenger welfare Adaptive large neighborhood search metaheuristic

a b s t r a c t Railway planning is a complex activity which is usually decomposed into several stages, traditionally network design, line design, timetabling, rolling stock, and staffing. In this paper, we study the design and optimization of train timetables for a rail rapid transit (RRT) line adapted to a dynamic demand environment, which focuses on creating convenient timetables for passengers. The objective is to minimize the average passenger waiting time at the stations, thus focusing on passenger welfare. We first propose two mathematical programming formulations which generalize the non-periodic train timetabling problem on a single line under a dynamic demand pattern. We then analyze the properties of the problem before introducing a fast adaptive large neighborhood search (ALNS) metaheuristic in order to solve large instances of the problem within short computation times. The algorithm yields timetables that may not be regular or periodic, but are adjusted to a dynamic demand behavior. Through extensive computational experiments on artificial and real-world based instances, we demonstrate the computational superiority of our ALNS compared with a truncated branch-and-cut algorithm. The average reduction in passenger waiting times is 26%, while the computational time of our metaheuristic is less than 1% of that required by the alternative CPLEX-based algorithm. Out of 120 open instances, we obtain 84 new best known solutions and we reach the optimum on 10 out of 14 instances with known optimal solutions. Ó 2014 Elsevier Ltd. All rights reserved.

1. Introduction Railway planning is a complex activity which is usually decomposed into several stages, including network design, line design, timetabling, rolling stock, and staffing (Cordeau et al., 1998; Guihaire and Hao, 2008). Traditionally, these problems have been optimized based on the operator’s perspective. In this paper, we provide a new approach, which takes the operator’s restrictions into account, but maximizes the passenger welfare, measured as the average waiting time (AWT) at the stations. The train timetabling problem solved in this paper consists of determining departure and arrival times of each train service to and from each station along a single line rail rapid transit (RRT) line over a long planning horizon characterized by a dynamic passenger arrival pattern. ⇑ Corresponding author at: HEC Montréal, 3000 chemin de la Côte-Sainte-Catherine, Montréal H3T 2A7, Canada. Tel.: +1 514 343 6111x8707; fax: +1 514 343 7121. E-mail addresses: [email protected] (E. Barrena), [email protected] (D. Canca), [email protected] (L.C. Coelho), [email protected] (G. Laporte). http://dx.doi.org/10.1016/j.trb.2014.08.013 0191-2615/Ó 2014 Elsevier Ltd. All rights reserved.

E. Barrena et al. / Transportation Research Part B 70 (2014) 134–150

135

A service is defined as a trip from an origin to a final destination. Here, a train refers to the service it operates. We consider the case of a two-directional two-track rail rapid transit system, where trains in opposite directions stop at different platforms, which means that departure and arrival times can be designed independently for each direction without interaction at stations (Vuchic, 2007, pp. 304–305). This is a typical design observed in many RRT systems around the world. We analyze the properties of a mathematical model and, on this basis, we develop an efficient heuristic for the timetabling problem, yielding solutions that are adapted to the demand pattern. The proposed heuristic works with small discrete time intervals, as fine as needed. In what follows, we review the relevant literature before positioning our contributions. When designing railway timetables for passengers it is essential to respect the track capacities, and it is common to optimize an objective function relevant to the operator or to the infrastructure manager, such as minimizing the deviations from an ideal timetable, or minimizing the services’ running times (Cacchiani and Toth, 2012). This can be achieved in two ways. The first is to construct a periodic timetable, i.e., a timetable that is constantly repeated, for example with departures at 03, 21 and 46 min every hour. Periodic timetabling models typically follow the Periodic Event Scheduling Problem formulation of Serafini and Ukovich (1989). Periodic timetables have the advantage of being easily memorized by passengers and have proved their ability to deal with large-scale railway networks (Kroon et al., 2009). The alternative, which consists of constructing a non-periodic timetable, is appropriate when demand cannot be assumed to be constant over time or when the service operates in long corridors with a high track density. Several different integer linear programming models have been proposed for the non-periodic timetabling problem. The model of Carey (1994) uses binary variables to describe the precedences between services, and continuous variables to represent departure and arrival times. It minimizes the total service running time. Caprara et al. (2002) have used a time discretization and a representation of the problem on a time–space graph. These authors have applied their model to a single one-way line. Their objective consists of minimizing the deviation with respect to an ideal timetable. Vansteenwegen and Van Oudheusden (2006) have studied the problem of improving passenger service on a small part of the Belgian railway network, taking into account waiting times and delays. They define a set of ideal train buffer times to ensure connections and minimize deviations with respect to them. None of the above contributions considers passenger demand in the design of the timetables. However, this problem has been studied by Ceder (2001, 2009) in a bus timetabling context. This author does not consider a demand curve as a function of time, but estimates the demand from the bus loads and aims to obtain the desired passenger loads instead of minimizing the AWT. Heuristics have been extensively employed to solve variants of the train timetabling problem, including greedy heuristics with different objectives (Cai et al., 1998; Yuan et al., 2008), mathematical programming-based methods to minimize interchange waiting times (Wong et al., 2008), graph theory-based algorithms (Liebchen and Möhring, 2002), cut-based algorithms (Liebchen, 2005), and sequential decomposition (Caimi et al., 2009). Local search algorithms include a hybrid of tabu search and branch-and-bound (Corman et al., 2010), a combination of insertion, backtracking and dynamic route selection (Burdett and Kozan, 2010b), and a hybridization of simulated annealing with particle swarm optimization (Jamili et al., 2012). Population-based algorithms include ant colony optimization (Chen et al., 2011; Reimann and Leal, 2013) and genetic algorithms (Chung et al., 2009; Nachtigall and Voget, 1996; Tormos et al., 2008). Constraint satisfaction techniques have been used by Odijk (1996), Abril et al. (2006), and Ingolotti et al. (2006), whereas methods based on linear programming relaxations and duality-based methods have been frequently applied (Brännlund et al., 1998; Cacchiani et al., 2008; Cacchiani et al., 2010a; Cacchiani et al., 2010b; Caprara et al., 2006; Lee and Chen, 2009; Zhou and Zhong, 2007). Other recent heuristics are those of Burdett and Kozan (2010a) and of Furini and Kidd (2013), which are based on a disjunctive graph representation of the timetable and on relaxed dynamic programming, respectively. Most of these papers optimize an objective function relevant to the service provider. From the infrastructure manager point of view, one commonly used objective is to minimize the deviation from a solution proposed by the operator. This objective is frequently used in the periodic case. From the users’ perspective, the objective of minimizing total travel time has been frequently considered. However, these contributions do not explicitly take variable passenger demand into consideration. Passenger demand is commonly assumed to be constant and is represented by a set of origin–destination (OD) matrices corresponding to different hours of a typical day. With respect to demand behavior, the following contributions have focused on its dynamic structure. Hänseler et al. (2012) investigate ways of computing dynamic OD matrices for train stations. Specifically, these authors have designed an algorithm to predict the OD demand as a function of the train timetabling. Yano and Newman (2001) propose a dynamic programming algorithm for the timetabling of trains used to transport containers arriving dynamically at the origin. Cordone and Redaelli (2011) determine the optimal periodic timetable when the train competes with an alternative mode that offers services between the same origin and destination. They consider a logit modal split model to calculate the demand captured by each transportation mode. Finally, Niu and Zhou (2013) study the train timetabling problem under a dynamic demand scenario, and are interested in how the system behaves under congestion, namely when the demand exceeds the capacity. These authors propose a local search heuristic and a dynamic programming algorithm in order to optimize the timetable for a single station. A genetic algorithm is developed and computation experiments are performed on a real case with 13 stations. In this paper, we consider a generalization of the non-periodic train timetabling problem under a dynamic demand pattern (Canca et al., 2014), which arises in an RRT context. We first provide two new non-linear formulations to precisely describe this problem in mathematical terms. The objective is to minimize passenger AWT at the stations, and the solutions are train timetables adjusted to a dynamic demand pattern which may neither be regular nor periodic. Barrena et al. (2014) have already proposed equivalent formulations and an exact branch-and-cut algorithm applicable to this model, but it can

136

E. Barrena et al. / Transportation Research Part B 70 (2014) 134–150

only solve relatively small instances to optimality. Here, we take advantage of Riemann’s sums theory to analyze the formulations of the problem and thus extract useful information in order to design a solution algorithm that minimizes passenger AWT. More specifically, we propose an adaptive large neighborhood search (ALNS) metaheuristic capable of solving larger and more realistic instances of the problem. We compare the solutions obtained by our ALNS algorithm with those of a truncated branch-and-cut algorithm (Barrena et al., 2014), both in terms of passenger AWT and of computing time. The remainder of this paper is organized as follows. In Section 2 we formally describe the problem and present two mathematical formulations for it. In Section 3 we exploit the mathematical properties of the second model to obtain useful information for the design of the solution algorithm. We then describe our ALNS metaheuristic in Section 4. We present the results of extensive computational experiments in Section 5, and conclusions in Section 6. 2. Problem description and mathematical models We focus on constructing train timetables adapted to a dynamic demand pattern using the same criteria as in Barrena et al. (2014) who represented them in the form of time–space diagrams, as shown in Fig. 1. The horizontal axis represents the planning horizon, and the vertical axis the distance from each station to the first one. Fig. 1(a) illustrates a regular timetable, i.e., the headway between consecutive trains is constant, whereas Fig. 1(b) depicts a non-regular timetable. Let S ¼ f1; . . . ; ng be the ordered set of stations defining a two-track railway line. The planning horizon is discretized into time intervals of length d. Thus, the time instant t 2 T ¼ f0; 1; . . . ; pg corresponds to dt time units elapsed since the beginning of the planning horizon. The discretization constant d represents the length of the smallest time interval considered in the t problem and therefore, from now on we will consider it as the time unit which can be as small as desired. Let dij be the passenger demand between stations i; j 2 S; i < j during the interval ðt  1; t. Let lij be the length of the segment between stations i and j; hmin the minimum headway, i.e., the minimum amount of time required between the departure of two consecutive trains at each station, wmin and wmax the minimum and maximum allowed dwell time at stations, and smin and smax the inverse of the maximum and minimum traveling speed of a train. The inverse of the speed is used to avoid non-linearities in the constraints of the problem. The operator has a set M ¼ f1; . . . ; mg of trains available. The aim of the problem is to determine train departure times at stations and train speeds on rail segments in order to minimize the passenger AWT at the stations. In what follows, we present the demand pattern and two non-linear formulations for this problem. 2.1. Demand function and passenger average waiting time In order to determine an optimal timetable under a dynamic demand behavior, we consider the cumulative demand function for each pair of stations. Several authors deal with variable demand, describing it through functions defined from different perspectives. Hurdle (1973) starts from a stepwise cumulative demand function and assumes that the passengers’ arrival pattern can be reasonably approximated by a differentiable function. The author does not require the algebraic expression of this function, but only that it can be drawn. Unlike in Hurdle (1973), our original data are not given in a regular form, the steps height does not correspond to a passenger, but data given at arbitrary points in time define the steps heights, each corresponding to a different number of passengers. We also approximate the cumulative demand by a smooth function which is a sum of sigmoidal functions (Canca et al., 2014). Daganzo (1997) presents a cumulative plot approach as an efficient way to represent input–output data, but recognizes that in practical applications, data are not always available and one usually has to construct the plot from limited information. This is what happens in our case, where data are irregularly given and we use a sum of sigmoidal functions to approximate the cumulative demand for each origin–destination pair (Canca et al., 2014). As Daganzo (1997) states, Newell (1971) was the author who demonstrated the full potential of cumulative plots as an analysis and thinking tool in connection with his queueing theory work. Starting from a given cumulative demand function, Newell (1971) proposes a formula for the waiting time of passengers, which takes into consideration the time elapsed between the passenger arrival time and the

(a) Regular timetable

(b) Non-regular timetable

Fig. 1. Time–space diagrams of train timetables for a one line corridor. Source: Barrena et al. (2014).

E. Barrena et al. / Transportation Research Part B 70 (2014) 134–150

137

departure of the next train after the arrival time of each passenger. This total waiting time is given by the area between the demand curve and a step function. In order to obtain an easy-to-manage objective function in the optimization problem, we consider an approximation of the area under the cumulative demand curve by calculating the area of the triangles formed by the cumulative demand values between the departure times of two consecutive trains. This is a reasonable approximation, taking into account the fact that the cumulative demand is a monotone increasing function and it is not an exact description of the demand since it is the result of an approximation from demand data given at random points in time. This simplification counts all passengers arriving between the departure of two trains and considers that each passenger arriving in this interval waits on average half of the interval time.

2.2. Formulation 1 The first formulation uses binary variables yk ; k 2 M, equal to one if and only if train k is used in the solution. Let the launch time of train k at station i be represented by an integer variable xik . Let Di be the cumulative function of the outgoing P P t0 demand at station i, i.e., Di ðtÞ ¼ tt0 ¼0 j2S;j>i dij . Note that we combine the demand over different destinations since we measure the passenger welfare in terms of the AWT. We consider a dynamic demand pattern where the original demand funct tions dij have certain peaks at rush hours, and the cumulative demand functions Di can therefore be considered as a sum of sigmoidal functions (Canca et al., 2014), each sigmoidal function representing the outgoing demand from station i to each of the others stations j 2 S; j > i. These variables give rise to a non-linear objective function expressed in terms of the launch times of each pair of consecutive trains. In order to write the following model, we consider two dummy trains, one at the beginning and another at the end of the planning horizon. The problem can then be formulated as follows:

ðF1Þ minimize

 1X X  Di ðxik Þ  Di ðxik1 Þ ðxik  xik1 Þ 2 i2S k2M[fmþ1g

ð1Þ

subject to

xi0 ¼ 0 i 2 S

ð2Þ

ximþ1 ¼ p i 2 S

ð3Þ

pð1  yk Þ 6 xik 6 xikþ1

xik

6 p i 2 S;

i 2 S;

xkiþ1 ¼ xik þ ski li;iþ1 þ wkiþ1 yk wmin 6 wki 6 yk wmax yk smin 6 ski 6 yk smax xikþ1 P xik þ hmin ykþ1 xik 2 Nþ

i 2 S;

k2M

ð4Þ

k 2 M n fmg i 2 S n fng; i 2 S; i 2 S;

k2M

k2M

i 2 S n fng;

k2M

yk 2 f0; 1g k 2 M:

ð5Þ

k2M

k 2 M n fmg

ð6Þ ð7Þ ð8Þ ð9Þ ð10Þ ð11Þ

In this formulation, the objective function (1) minimizes the total passenger AWT, which must be divided by the total P demand DT ¼ i2S;t2T Di ðtÞ in order to obtain the AWT per passenger. Since the total demand is a constant, we do not consider it in the objective function and from now on we will refer to the objective function (1) as the passenger AWT. Constraints (2) and (3) define the dummy train launch times whose role is to ensure that the objective function takes all waiting times and demands into account. Constraints (4) link the variables xik and yk . Specifically, if train k is not used, i.e., yk ¼ 0, then it is launched at the end of the planning horizon, i.e., xik ¼ p. On the other hand, if train k is used, i.e., yk ¼ 1, then xik can take any value smaller than p, i.e., the train must be launched before the end of the planning horizon. Constraints (5) order the trains in terms of launch times, i.e., together with constraints (4), they ensure that trains with lower index values are launched first. This is a necessary condition for this formulation, but it also breaks equivalent solutions (solutions differing only by train indices and having the same objective function). Constraints (6) ensure that if train k is launched from station i at time xik , then this train has to be launched from the next station i þ 1 at time xik , plus the time required to travel from i to i þ 1, plus the dwell time at station i þ 1. Constraints (7) set the bounds for dwell times of launched trains and consider zero dwell times for those that are not launched. Constraints (8) are similar to (7) and bound the inverse of the speed ski of train k leaving station i, and set it to zero if train k is unlaunched. Note that constraints (6)–(8) ensure that if train k is not launched, then xik ¼ p for all stations i. Constraints (9) play a stronger role than (5) and (6) by ensuring that no two trains are launched from the same station at the same time, and by also considering the minimum headway hmin . Integrality and binary conditions on the variables are enforced through constraints (10) and (11). Note that we declare xik as integer variables since they are used in the objective function as the argument of the demand function defined over a discrete set.

138

E. Barrena et al. / Transportation Research Part B 70 (2014) 134–150

2.3. Formulation 2 Our second formulation is similar to the first but disregards speed and waiting times variables. The previous formulation contains a number of equivalent solutions since the same launch time at a station i þ 1 can be obtained by the launch time at the previous station i, and different combinations of arrival and dwell times at station i þ 1. This situation can be avoided by not considering the speeds ski and dwell times wki of each train k, but constraining the launch times at stations to lie within the interval in which a service must be launched in order to yield a feasible solution. In this way, one would obtain the same launch times at stations as in F1, avoiding equivalent solutions. This is achieved in our second formulation:

ðF2Þ minimize

 1X X  Di ðxik Þ  Di ðxik1 Þ ðxik  xik1 Þ 2 i2S k2M[fmþ1g

ð12Þ

subject to 2, 3, 4 and 5, 9, 10 and 11, and to

xiþ1 6 xik þ yk ðsmax li;iþ1 þ wmax Þ i 2 S n fng; k

k2M

ð13Þ

xik þ yk ðsmin li;iþ1 þ wmin Þ 6 xiþ1 k

k 2 M:

ð14Þ

i 2 S n fng;

In this formulation, the objective function (12) is the same as (1). Constraints (13) and (14) enforce the bounds on the launch time from station i þ 1 according to the launch time from station i, and the minimum and maximum speeds and waiting times. 3. Analysis of the mathematical properties of the problem Formulation 2 is easier to solve than Formulation 1 under a continuous demand pattern, as in Canca et al. (2014). However, the cumulative demand functions Di are the sum of sigmoidal functions, which means that the objective function (12) is non-linear and non-convex, and therefore remains beyond the reach of exact algorithms for non-trivial instances. The train timetabling problem studied by Caprara et al. (2002) is NP-hard, which explains why it cannot be solved exactly for large instances. The proof of NP-hardness is based on a polynomial reduction of the NP-hard Max-Independent Set Problem (MISP) (Garey and Johnson, 1979). In the case of Caprara et al. (2002), a MISP instance is reduced to an initial (ideal) timetable obtaining the Max-Independent set. Our case is different because we do not have an ideal timetable to be fitted as much as possible. Suppose an instance where demand peak is concentrated in a short time period. Then, relaxing the safety headway constraints and allowing identical departure and arrival times for different trains will yield a probably infeasible ‘‘ideal’’ timetable minimizing waiting time where an important number of trains will be incompatible. Then, our problem consists of designing a new timetable, as close as possible to the ‘‘ideal’’ one, but taking into consideration all the constraints previously relaxed. Moreover, a new set of constraints will be imposed, forbidding overtaking at stations. The resulting problem is very similar to the problem described in Caprara et al. (2002). The version of the problem we consider can only be solved exactly for instances with three stations and at most 300 time units in the planning horizon (Barrena et al., 2014). Therefore the use of a heuristic becomes necessary for larger and more realistic instances. In this section we analyze the mathematical properties of the problem, which then provides some directions and motivations for the heuristic to be presented in Section 4. Specifically, in this section we apply Riemann’s theory (Riemann, 1867) to gain insights into the problem. To this end, we first introduce some definitions and properties. Definition 3.1. An n-order partition of an interval ½a; b is a finite set P ¼ fx0 ; x1 ; . . . ; xn g such that a ¼ x0 < x1 < x2 < . . . < xn1 < xn ¼ b. Intervals ½xi ; xiþ1 ; i ¼ 0; 1; . . . ; n  1, are called intervals of the partition P. If all these intervals have the same length, the partition is called regular. The diameter or mesh size of the partition P is defined as dðPÞ ¼ maxfjxkþ1  xk j; i ¼ 0; 1; . . . ; n  1g. Let P and P 0 be two partitions of ½a; b. Then P 0 is a refinement of P if P # P 0 .

Lemma 3.1. Let P # P 0 be two partitions of a closed interval ½a; b. It follows directly that dðP 0 Þ 6 dðPÞ. Definition 3.2. Let f ðxÞ be a real function defined on a closed support ½a; b 2 R, and let P ¼ fx0 ; x1 ; . . . ; xn g be an n-order partition of ½a; b. The Riemann lower sum of f for the partition P is defined as

sðP; f ; ½a; bÞ ¼

n X

mk ðxk  xk1 Þ;

k¼1

where mk ¼ inf ff ðxÞ : x 2 ½xk1 ; xk g; 0 6 k 6 n  1, and the Riemann upper sum of f is defined as

SðP; f ; ½a; bÞ ¼

n X Mk ðxk  xk1 Þ; k¼1

where M k ¼ supff ðxÞ : x 2 ½xk1 ; xk g; 0 6 k 6 n  1.

139

E. Barrena et al. / Transportation Research Part B 70 (2014) 134–150

Geometrically, sðP; f ; ½a; bÞ is the sum of the areas of the rectangles based on the intervals ½xk ; xkþ1  lying below the graph of f, and SðP; f ; ½a; bÞ is the sum of the areas of the rectangles lying above the graph of f. Theorem 3.1. Let f : ½a; b ! R be a real function bounded in the variable domain, let P and P 0 be partitions of ½a; b such that P # P 0 . Then

sðP; f ; ½a; bÞ 6 sðP 0 ; f ; ½a; bÞ and

SðP 0 ; f ; ½a; bÞ 6 SðP; f ; ½a; bÞ: Proof. See Apostol (1957) and Riemann (1867).

h

3.1. Analysis of the problem and of the objective function We analyze the problem based on our second formulation (Section 2.3), and we observe that at each station i 2 S, the departure time variables xik of the trains k 2 M induce an ðm þ 1Þ-order partition P i of the time horizon ½0; p, that is, P i ¼ f0 ¼ xi0 ; xi1 ; . . . ; ximþ1 ¼ pg. The problem consists of finding the partition of the planning horizon at each station such that the objective function (12) is minimized, while ensuring feasibility. We now analyze the objective function by applying Riemann’s theory in order to obtain insightful information that will drive the design of a heuristic algorithm. We first rewrite the objective function (12) as

ðOFÞ

 1X X  1X Di ðxik Þ  Di ðxik1 Þ ðxik  xik1 Þ ¼ 2 i2S k2M[fmþ1g 2 i2S

X k2M[fmþ1g

Di ðxik Þðxik



xik1 Þ



X

! Di ðxik1 Þðxik



xik1 Þ

:

k2M[fmþ1g

The function Di ðtÞ represents the cumulative demand at each station i 2 S and is therefore monotonically increasing. It follows that for each interval ½xik1 ; xik , the Riemann lower and upper sums are equal to the Riemann left and right sums, respectively, i.e., mk ¼ Di ðxik1 Þ and M k ¼ Di ðxik Þ. The objective function can then be written as

ðOFÞ ¼

1X ðSðP i ; Di ; ½0; pÞ  sðP i ; Di ; ½0; pÞÞ: 2 i2S

ð15Þ

This means that for all i 2 S, the objective function (15) is minimized when the differences SðP i ; Di ; ½0; pÞ  sðP i ; Di ; ½0; pÞ are minimized. The problem then consists of finding the best partitions P i ¼ f0 ¼ xi0 ; xi1 ;    ; ximþ1 ¼ pg of the planning horizon at each station i 2 S such that the difference between the upper and lower sums of cumulative demand functions Di ðtÞ is minimized, while ensuring feasibility. The order (m þ 1) of these partitions is given by the number of available train services. Fig. 2 provides a graphical representation of the problem encountered for one station. The cumulative demand functions Di ðtÞ are the sum of sigmoidal functions (see Canca et al. (2014)), and the areas with maximum slope are an indicator of the

Fig. 2. Cumulative demand function at one station i 2 S. In light grey, the Riemann lower sums; in dark, the double of the objective function at station i, represented by SðP i ; Di ; ½0; pÞ  sðP i ; Di ; ½0; pÞ.

140

E. Barrena et al. / Transportation Research Part B 70 (2014) 134–150

peak hour demand. The light rectangles represent the Riemann lower sums. Above these, we show in dark the difference between the Riemann upper and lower sums. These dark areas are a visual representation of twice our objective function at each station. 3.2. Effect of train insertions and removals In this section we prove the effect of train insertions and removals in the passenger AWT. Theorem 3.2. Let P i ; i 2 S be a partition of the time horizon ½0; p at each station i, representing a timetable with a set M0  M of train services. When a train k 2 M is inserted in the timetable, the passenger AWT decreases. In contrast, the passenger AWT 0 increases when a train k 2 M0 is removed from the timetable.

Proof. When a train k 2 M is inserted in the timetable, it departs from every station i 2 S and therefore, a point xik is added to every partition P i . So, a refined partition P 0i of P i is obtained at each station and, according to Theorem 3.1, the following inequalities hold:

sðP i ; Di ; ½0; pÞ 6 sðP 0i ; Di ; ½0; pÞand SðP 0i ; Di ; ½0; pÞ 6 SðP i ; Di ; ½0; pÞ: Since Di ðtÞ P 0 for all t 2 T , the Riemann upper and lower sums are also positive, and we can conclude that

SðP 0i ; Di ; ½0; pÞ  sðP 0i ; Di ; ½0; pÞ 6 SðP i ; Di ; ½0; pÞ  sðP i ; Di ; ½0; pÞ;

for all i 2 S:

This implies that the AWT (Eq. (15)) also satisfies

 1X 1 X SðP 0i ; Di ; ½0; pÞ  sðP 0i ; Di ; ½0; pÞ 6 ðSðP i ; Di ; ½0; pÞ  sðP i ; Di ; ½0; pÞÞ: 2 i2S 2 i2S Analogously, when a train is removed, it yields coarser partitions and therefore the objective function value increases. h In fact, when the number of intervals of the partition tends to infinity, the Riemann upper and lower sums become equal and this value is called the Riemann integral (Apostol, 1957). In the train timetabling problem, this means that when the number of trains tends to infinity, the passenger AWT tends to zero, which is obvious. 3.3. Selection of the partition or timetable According to Theorem 3.2, the objective function value decreases with every train insertion so that an infinite number of trains would make the passenger AWT equal to zero. Since there are limitations on the number of trains, which are given by the minimum headway, and by the available number of trains, a good partition should be selected in order to obtain a small AWT. In order to derive good partitions, we first show that a regular partition, i.e., a regular timetable in the terms of the problem being studied, is not necessarily optimal, especially when the demand is non-regular. Counter-example 3.1. Consider the following cumulative demand function for a problem with only one departure station

DðtÞ ¼

8 > < > :

0

if t 2 ½0; 100;

t  100 if t 2 ð100; 200; 100

if t 2 ð200; 900;

which comes from a non-regular demand function. According to this demand, no passenger arrives during the first interval ½0; 100, one passenger per time unit arrives during the interval ð100; 200, and no passenger arrives during the interval ð200; 900. Let P 1 ¼ f0; 150; 300; 450; 600; 750; 900g be a regular six-order partition of the time horizon ½0; p, and P 2 ¼ f0; 100; 125; 150; 175; 200; 900g be a non-regular six-order partition of ½0; p. The points in the partition indicate the trains departure time and they are represented by vertical dashed lines in Fig. 3. In this figure, the dark areas represent the differences between Riemann upper and lower bounds for partitions P 1 and P 2 , whose values are

SðP 1 ; D; ½0; pÞ  sðP 1 ; D; ½0; pÞ ¼ 15000; and

SðP 2 ; D; ½0; pÞ  sðP 2 ; D; ½0; pÞ ¼ 2500: This means that the objective function value, i.e., the AWT, obtained with the regular partition is six times larger than the one obtained with the selected non-regular one.

E. Barrena et al. / Transportation Research Part B 70 (2014) 134–150

(a) Regular timetable P 1

141

(b) Non-regular timetable P 2

Fig. 3. Graphical representation of SðP l ; D; ½0; 900Þ  sðP l ; D; ½0; 900Þ; l ¼ 1; 2.

So, a regular timetable does not necessarily minimize the difference between the Riemann upper and lower sums, and consequently the passenger AWT. In  order to reduce the value of the objective function (12), we attempt to reduce the largest sums. Since ðxik  xik1 Þ P 0 and Di ðxik Þ  Di ðxik1 Þ P 0, their product can be large due to large individual values of each of them. In what follows, we analyze each case separately. 3.3.1. Large values of the partition interval lengths ðxik  xik1 Þ According to Definition 3.1, the maximum of the interval lengths at each station is the diameter of the partition P i . Following Lemma 3.1, refining the partition reduces the diameter, which in turn decreases the value of the objective function. In order to refine the partition, one needs to add a point to it at some time instant in the interval corresponding to the largest diameter over the stations and obtain a partition P 0i so that dðP 0i Þ 6 dðP i Þ. We will later exploit this observation in Section 4.2.2 when designing the operators of our local search metaheuristic.   3.3.2. Large values of the cumulative demand Di ðxik Þ  Di ðxik1 Þ According to Theorem 3.1, the finer the partition, the smaller the difference between the Riemann upper and lower sums. This means that refining the partition within the partition interval having the highest increment in demand will reduce the difference in the sums and thus decrease the value of the objective function. This observation is exploited in Section 4.2.2. 3.4. Insights and directions After having analyzed the problem in terms of its mathematical properties, we have gained some insights into how the optimization algorithm should work if passenger welfare is to be maximized. According to Theorem 3.2, the algorithm should insert trains in order to reduce the value of the objective function, i.e., the passenger AWT. As demonstrated in Sections 3.3.1 and 3.3.2, insertions in certain areas of the planning horizon have a much stronger effect in reducing the AWT than in some other areas. 4. Adaptive large neighborhood search metaheuristic What is needed to solve our problem is an algorithm that removes trains from the time–space diagram and inserts new ones. One such algorithm is the adaptive large neighborhood search (ALNS) metaheuristic (Ropke and Pisinger, 2006) which is based on destroy and repair operators, and appears to be an ideal framework for our problem. In our context, destroy operators remove trains, while repair operators insert them back. As per Theorem 3.2, every time an insertion operator is selected, the objective function value is reduced, while it increases whenever a train is removed. This means that if only the improving moves were to be accepted by the heuristic, the removals operators would never be used, and the algorithm would be dysfunctional. This can be remedied by using a probabilistic acceptance criterion such as the one used in simulated annealing since it accepts deteriorating solutions with some probability and allows a better exploration of the search space. The ALNS metaheuristic enables us to simultaneously select good partitions P i at each station i 2 S for the demand functions Di , while ensuring feasibility. In our problem, feasibility is defined by the minimum and maximum train speed and dwell time at stations, which restricts the set of the

142

E. Barrena et al. / Transportation Research Part B 70 (2014) 134–150

possible departure times from consecutive stations, as well as by the limitation in the headway between consecutive trains. We have designed a powerful ALNS metaheuristic in order to maximize the passenger welfare, capable of handling a large number of origin–destination pairs and of simultaneously determining departure times, speeds, stopping times, while minimizing the passenger AWT. The destroy and repair operators are randomly selected at each iteration. Basically, train services are removed and inserted in the initial time space diagram by means of destroy and repair operators. A roulette wheel mechanism controls the choice of the operators, with a probability that depends on their past performance. More concretely, to each operator i are associated a score pi and a weight xi whose values depend on the past performance of the operator. Then, P given h operators, operator j will be selected with probability xj = hi¼1 xi . Initially, all weights are set to one and all scores are set to zero. The search is divided into segments of u iterations each, and the weights are computed by taking into account the performance of the operators during the last segment. At each iteration, the score of the selected operator is increased by r1 if the operator identifies a new best solution, by r2 if it identifies a solution better than the incumbent, and by r3 if the solution is not better but is still accepted. After u iterations, the weights are updated by considering the scores obtained in the last segment as follows: let oij be the number of times operator i has been used in the last segment j. The updated weights are then

xi :¼



xi

if oij ¼ 0

ð1  gÞxi þ gpi =oij

if oij – 0;

ð16Þ

where g 2 ½0; 1 is called the reaction factor and controls how quickly the weight adjustment reacts to changes in the operator performance. The scores are reset to zero at the end of each segment. As in other ALNS implementations (Coelho et al., 2012; Ropke and Pisinger, 2006), we use an acceptance criterion based on simulated annealing. Given a solution s, a neighbor solution s0 is always accepted if zðs0 Þ < zðsÞ, and is accepted with prob0 ability eðzðs ÞzðsÞÞ=s otherwise, where zðsÞ is the solution cost and s > 0 is the current temperature. The temperature starts at sstart and is decreased by a cooling rate factor / at each iteration, where 0 < / < 1. We now describe the main features of our algorithm. 4.1. Initial solution The algorithm can be initialized from an empty solution or from an arbitrary solution, for example, a regular solution as in Fig. 1(a). In our implementation, we start with an empty solution. Later, we assess the effect of providing an initial solution on the performance of the algorithm. 4.2. List of operators When designing the operators, we have paid close attention to the mathematical properties of a solution, as explained in Section 3, and to the need to diversify the search. Given a solution, destroy and repair operators will delete and insert train services in the time–space diagram. Restrictions on the headway, dwell times and speeds are considered by each operator in order to maintain the feasibility of the solution. We now list the destroy and repair operators we have developed. In what follows, q1 is the number of train services removed from the solution and q2 is the number of train services inserted in the solution at each iteration. The parameters q1 and q2 are integers randomly drawn from the interval ½1; mr  ðr ¼ 1; 2Þ, where m1 ¼ ms ; m2 ¼ m  ms , and ms is the number of train services of the incumbent solution s. More precisely, qr follows a semiqffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi triangular distribution with a negative slope, i.e., qr ¼ bmr  ð1  uÞðmr  1Þ2 þ 0:5c, where u is a random value in the interval ½0; 1. 4.2.1. Destroy operators 1. Randomly remove q1 train services This operator randomly selects q1 trains and removes them. It is useful for refining the solution since it does not change the solution much when q1 is small, which happens frequently due to the shape of the probability distribution of q1 . However, it still yields a major transformation of the solution when q1 is large. 2. Remove q1 train services from the smallest interval This operator identifies the two consecutive trains with the smallest interval and removes the earlier one, that is, it removes the train k  1 corresponding to the smallest values of the partition interval lengths ðxik  xik1 Þ (see Section 3.3.1). This procedure is repeated q1 times. 3. Remove q1 train services with the smallest demand This operator removes the train with the smallest passenger demand in one of its tracks, that is, it removes the train k   yielding the smallest Di ðxik Þ  Di ðxik1 Þ (see Section 3.3.2). It is applied q1 times.

E. Barrena et al. / Transportation Research Part B 70 (2014) 134–150

143

4.2.2. Repair operators 1. Randomly insert q2 train services This operator randomly inserts q2 trains. Each insertion is achieved by randomly selecting a time instant from the planning horizon at the first station and inserting a train starting at this time instant with random speed and dwell times at the next stations, while ensuring feasibility. Feasibility is ensured by respecting the speed and dwell bounds: starting from station i 2 S n fng and according to constraints (13) and (14), we calculate the feasible interval, i.e., the interval whose lower bound is given by the maximum speed and minimum dwell time and whose upper bound is given by the minimum speed and maximum dwell time. Any time instant within this interval is a feasible departure time from the next station i þ 1. 2. Insert q2 train services in the largest interval This operator inserts a train in the largest inter-departure interval, that is, in the interval corresponding to the largest values of the partition interval lengths ðxik  xik1 Þ (see Section 3.3.1). Inserting trains in this interval helps reduce the AWT and at the same time avoid large intervals of the planning horizon without train departures. This is important to guarantee service quality since, even if there are not many passengers waiting at a large interval, it is important to provide minimum services for them. Insertions are performed at a random time instant at the first station of this interval, and speed and dwell times at the following stations are randomly selected within their bounds while ensuring feasibility. This process is repeated q2 times. 3. Insert q2 train services near those having the largest demand This operator inserts a train just before the train with the most loaded track, that is, before the train k yielding the largest   Di ðxik Þ  Di ðxik1 Þ (see Section 3.3.2). Inserting trains in these high-demand areas helps reduce the AWT and at the same time implies that train frequency is increased at demand peak hours, which are normally the rush hours, and therefore it is convenient that passengers do not have to wait long during these periods. The insertion is carried out by assigning the departure of the new train from the first station just hmin time units before the departure of the most loaded train from first station. The new train runs parallel to the most loaded one, i.e., with the same speed and the same dwell time at the same stations. This process is repeated q2 times. 4.3. Parameter settings and pseudocode We now describe the parameter settings we have used in our ALNS implementation. These were set after an early tuning phase. The maximum number of iterations imax depends on the starting temperature sstart and on the cooling rate /. We have set these parameters as follows:

sstart ¼ 60; 000 / ¼ ð0:01=sstart Þ1=imax :

ð17Þ ð18Þ

This makes the cooling rate a function of the desired number of iterations, adjusting accordingly the probability that the ALNS mechanism will accept worsening solutions. The stopping criterion is satisfied when the temperature reaches 0.01. In our implementation, the maximum number of iterations imax was set to 70,000, the segment length u was set to 200 iterations, and the reaction factor g was set to 0.7, thus defining the new weights by 70% of the performance on the last segment and 30% of the last weight value. The scores are updated with r1 = 10, r2 = 5 and r3 = 2. Algorithm 1 shows the pseudocode of our ALNS implementation. In our experiments, we consider 20 reheatings, that is, after finishing the algorithm, the temperature is set to sstart and the process is repeated 20 times. 4.4. Possible extensions The aim of the paper is to efficiently solve a specific RRT problem. However, the proposed ALNS can be adapted to handle other features such as overtaking at stations (a common problem in case of metropolitan railways, at least in central stations) by maintaining a list with the train departure order at each station. Several train types can be also managed by considering different speed intervals, as well as a different stopping pattern for each train type. Note that in our formulation, a train can travel on any segment at any speed and can maintain a different dwell time at each station. An ALNS performs a feasibility test before applying an insertion operator. In this procedure, overtaking at segments is avoided. For specific segments with parallel tracks, this test can be relaxed, allowing the algorithm to manage multi-track segments. In the case of two-directional one-track segments, the feasibility test can be modified by comparing trains in both directions. 5. Computational experiments We now provide some implementation specifications, we describe the instances used and we present the results of extensive computational experiments. All computations were performed on a grid of Intel Xeon™ processors running at 2.66 GHz

144

E. Barrena et al. / Transportation Research Part B 70 (2014) 134–150

with up to 24 GB of RAM installed per node, with the Scientific Linux 6.1 operating system and using a single thread. Our algorithm was coded in C++ and uses a single processor and a maximum of only 1 GB of memory. Algorithm 1. ALNS metaheuristic 1: Initialize: set all weights equal to 1 and all scores equal to 0. 2: sbest s initial solution; s sstart ; reheatings 0. 3: while (s > 0.01 and reheatings < 20) do 4: s0 s. 5: Select a destroy and a repair operator using the roulette-wheel mechanism based on the current weights. Apply the operators to s0 and update the number of times they are used. 6: if zðs0 Þ < zðsÞ then 7: s s0 ; 8: if zðsÞ < zðsbest Þ then 9: sbest s; 10: update the score for the operators used with r1 ; 11: else 12: update the score for the operators used with r2 ; 13: end if 14: else 15: if s0 is accepted by the simulated annealing criterion then 16: s s0 ; 17: update the scores for the operators used with r3 . 18: end if 19: end if 20: if the iteration count is a multiple of u then 21: update the weights of all operators and reset their scores. 22: end if 23: s /s; 24: if s 6 0:01 then 25: reheatings reheatings þ 1. 26: end if 27: end while 28: return sbest ; 5.1. Set of instances We have used the set of instances initially proposed by Barrena et al. (2014) in order to compare the results of our metaheuristic with those obtained by their branch-and-cut algorithm. A real-world based instance was obtained for the line C5 of the Madrid Metropolitan Railway (Canca et al., 2014; Barrena et al., 2014), which qualifies as an RRT (Vuchic, 2005, page 595). The Spanish railway operator RENFE provided us surveys on the demand data between six stations of line C5: Móstoles-Soto, Móstoles, Las Retamas, Alarcón, San José, and Cuatro Vientos. We have used these survey data to adjust the 15 cumulative demand functions among these stations, and we have used these functions as inputs for our model. The set of artificial instances was generated according to the following parameters:         

number of stations n: 3, 6, 10; horizon p: 200, 400, 600, 800, 1200 min; discretization constant d: 1, 2, 4 min; maximum number of trains m: 5, 10; maximum inverse speed of the trains smin : 0.0015 min/m (speed = 40 km/h); minimum inverse speed of the trains smax : 0.00075 min/m (speed = 80 km/h); minimum headway hmin : 12 min; minimum stopping time at the stations wmin : 4 min; maximum stopping time at the stations wmax : 12 min.

These instances will be referred to as TT-n-p-d-m, e.g., TT-3-800-2-10, corresponding to a train timetabling instance with three stations, a planning horizon of 800 min, a discretization constant of two minutes, and a maximum of 10 trains. The C5 instances will refer to the six stations previously mentioned and their corresponding demand functions, varying the rest of parameters. These will be referred to as C5-p-d-m. The set of instances as well as their solutions are available on http:// www.leandro-coelho.com.

145

E. Barrena et al. / Transportation Research Part B 70 (2014) 134–150

Fig. 4. Cost of the solution accepted at the end of each iteration for instance TT-3-200-1-5.

Table 1 Summary of the average waiting times and running time for the two algorithms. Regular

Barrena et al. (2014)

AWT

AWT

Time (s)

ALNS AWT

Time (s)

3 stations 6 stations 10 stations Madrid C5

45.19 60.06 60.56 52.35

10.54 32.09 49.18 49.78

6656 10,803 10,803 10,802

10.04 21.09 27.15 45.71

33 85 148 73

Average

54.54

35.40

9766

25.99

84.75

Table 2 Summary of computational results of a regular timetable, of Barrena et al. (2014), and of our ALNS on instances with three stations. Instance

Regular

Barrena et al. (2014)

AWT

UB1

LB

Time (s)

AWT

ALNS UB2

Time (s)

AWT

ALNS vs. Barrena et al. (2014) Improvement (%)

TT-3-200-1-5⁄ TT-3-200-2-5⁄ TT-3-200-4-5⁄ TT-3-400-1-5 TT-3-400-2-5⁄ TT-3-400-4-5⁄ TT-3-600-1-5 TT-3-600-2-5⁄ TT-3-600-4-5⁄ TT-3-800-1-5 TT-3-800-2-5 TT-3-800-4-5⁄ TT-3-1000-1-5 TT-3-1000-2-5 TT-3-1000-4-5⁄

19.39 18.51 15.81 44.05 43.50 39.75 44.99 43.33 39.44 87.36 89.26 97.84 112.70 111.55 108.02

37,368 35,210 29,428 102,418 98,412 87,260 167,941 162,008 150,748 290,061 282,996 270,468 393,820 294,992 283,516

37,368 35,210 29,428 52,149 98,412 87,260 31,144 162,008 150,748 29,381 65,078 270,468 25,130 50,677 283,516

366 28 1 10,800 476 35 10,801 7979 906 10,802 10,801 3480 10,802 10,801 4627

6.27 6.90 6.93 9.08 9.72 9.73 12.34 12.90 13.07 17.01 17.60 17.85 23.10 18.30 18.62

37,368 35,210 29,428 102,456 98,412 87,260 167,175 162,008 150,748 290,068 282,996 270,468 301,576 294,992 283,516

20 10 5 32 17 8 49 25 13 60 31 15 74 39 20

6.27 6.90 6.93 9.08 9.72 9.73 12.29 12.90 13.07 17.01 17.60 17.85 17.69 18.30 18.62

0.00 0.00 0.00 -0.04 0.00 0.00 0.46 0.00 0.00 0.00 0.00 0.00 23.42 0.00 0.00

TT-3-200-1-10⁄ TT-3-200-2-10⁄ TT-3-200-4-10⁄ TT-3-400-1-10 TT-3-400-2-10 TT-3-400-4-10⁄ TT-3-600-1-10 TT-3-600-2-10 TT-3-600-4-10⁄ TT-3-800-1-10 TT-3-800-2-10 TT-3-800-4-10 TT-3-1000-1-10 TT-3-1000-2-10 TT-3-1000-4-10

6.48 8.17 8.17 22.48 19.24 19.24 24.47 21.75 21.75 40.43 43.65 43.66 55.08 52.76 52.76

28,961 26,366 20,744 59,715 55,714 45,592 99,034 83,512 70,460 167,218 168,932 127,528 205,169 158,180 130,460

28,961 26,366 20,744 29,887 54,734 45,592 17,418 43,803 70,460 19,581 29,814 54,541 16,550 26,658 53,162

153 18 2 10,801 10,800 140 10,801 10,801 8649 10,803 10,801 10,801 10,805 10,801 10,801

4.86 5.42 5.48 5.29 5.94 6.04 7.28 7.13 7.17 9.81 10.91 9.47 12.03 10.28 9.65

29,032 26,502 21,332 61,717 56,338 45,592 89,538 84,062 70,584 141,642 135,556 119,832 143,477 137,028 119,736

31 16 8 49 26 13 62 32 16 77 42 22 92 46 24

4.87 5.44 5.58 5.47 5.99 6.04 6.58 7.17 7.18 8.31 8.95 9.02 8.42 9.04 9.02

-0.25 -0.52 -2.83 -3.35 -1.12 0.00 9.59 -0.66 -0.18 15.30 19.76 6.03 30.07 13.37 8.22

Average

45.19

137,808

64,875

6656

10.54

129,188

32.47

10.04

3.91

146

E. Barrena et al. / Transportation Research Part B 70 (2014) 134–150

Fig. 5. Percentage AWT improvement on instances with six stations with respect to a regular timetable.

Table 3 Summary of computational results of a regular timetable, of Barrena et al. (2014), and of our ALNS on instances with six stations. Instance

Regular

Barrena et al. (2014)

AWT

UB1

LB

Time (s)

AWT

ALNS UB2

Time (s)

AWT

ALNS vs. Barrena et al. (2014) Improvement (%)

TT-6-400-1-5 TT-6-400-2-5 TT-6-400-4-5 TT-6-600-1-5 TT-6-600-2-5 TT-6-600-4-5 TT-6-800-1-5 TT-6-800-2-5 TT-6-800-4-5 TT-6-1000-1-5 TT-6-1000-2-5 TT-6-1000-4-5 TT-6-1200-1-5 TT-6-1200-2-5 TT-6-1200-4-5

38.94 39.42 39.42 38.76 39.27 39.29 91.73 92.21 92.22 93.19 93.73 93.72 145.75 146.20 146.14

697,886 618,322 575,780 1,583,040 1,446,760 1,091,000 4,501,100 2,331,550 1,815,690 2,655,010 2,281,690 1,904,200 8,653,510 3,079,810 1,973,860

31,357 89,743 235,328 64,565 105,349 260,442 44,831 52,831 216,420 48,904 63,335 244,483 39,394 67,542 133,449

10,801 10,800 10,801 10,802 10,800 10,800 10,803 10,801 10,801 10,804 10,801 10,801 10,806 10,802 10,801

17.70 16.66 16.58 29.12 27.57 22.03 71.63 38.05 30.83 42.25 37.26 32.24 137.72 49.94 33.35

657,026 642,086 610,176 1,272,434 1,252,694 1,205,644 1,930,834 1,908,152 1,853,772 2,009,144 1,987,872 1,936,144 2,055,965 2,035,826 1,984,188

66 34 18 109 54 29 114 65 31 151 78 38 181 84 47

16.67 17.26 17.45 23.41 24.01 24.13 30.73 31.32 31.44 31.98 32.59 32.75 32.72 33.35 33.51

5.85 3.84 5.97 19.62 13.41 10.51 57.10 18.16 2.10 24.33 12.88 1.68 76.24 33.90 0.52

TT-6-400-1-10 TT-6-400-2-10 TT-6-400-4-10 TT-6-600-1-10 TT-6-600-2-10 TT-6-600-4-10 TT-6-800-1-10 TT-6-800-2-10 TT-6-800-4-10 TT-6-1000-1-10 TT-6-1000-2-10 TT-6-1000-4-10 TT-6-1200-1-10 TT-6-1200-2-10 TT-6-1200-4-10

18.32 18.81 18.81 25.18 25.67 25.68 44.29 44.77 44.78 51.55 52.07 52.07 49.60 50.11 50.11

541,427 299,506 264,824 1,289,140 738,628 585,832 1,964,960 1,292,330 1,123,040 2,301,980 2,803,080 1,360,370 2,420,880 2,290,650 1,339,860

31,841 54,620 112,558 43,125 85,181 162,597 18,149 59,877 159,611 21,797 54,501 118,285 27,191 50,684 90,199

10,802 10,800 10,801 10,804 10,801 10,800 10,805 10,802 10,801 10,807 10,803 10,802 10,812 10,803 10,801

13.73 8.59 8.71 23.72 14.57 12.76 31.27 21.54 19.83 36.64 45.55 23.60 38.53 37.40 23.28

367,925 349,606 299,864 725,086 680,844 652,088 998,887 972,616 924,604 1,022,806 1,000,818 954,772 1,023,010 1,011,344 945,752

95 53 29 126 60 33 161 84 46 196 113 54 223 119 61

9.33 9.86 9.59 13.34 13.51 13.97 15.90 16.46 16.68 16.28 16.90 17.16 16.28 17.07 17.02

32.05 16.73 13.23 43.75 7.82 11.31 49.17 24.74 17.67 55.57 64.30 29.82 57.74 55.85 29.41

Average

60.06

1,860,857

92,940

10,802

32.09

1,175,733

21.09

22.12

85.07

5.2. Computational results We first analyze the performance of Algorithm 1 with respect to its parameters, after having performed various preliminary tests to adjust their settings. The following results are typical of those obtained across all instances. For presentation reasons we have chosen to illustrate them on instance TT-3-200-1-5 because this is a small instance for which the optimal solution is known (Barrena et al., 2014). In Fig. 4 we depict the cost of the incumbent solution at the end of each iteration. This figure illustrates how the algorithm converges and also how it can escape from local optima by accepting worsening

147

E. Barrena et al. / Transportation Research Part B 70 (2014) 134–150

Fig. 6. Percentage AWT improvement on instances with 10 stations with respect to a regular timetable.

Table 4 Summary of computational results of a regular timetable, of Barrena et al. (2014), and of our ALNS on instances with 10 stations. Instance

Regular

Barrena et al. (2014)

AWT

UB1

LB

Time (s)

AWT

UB2

Time (s)

AWT

43.59 44.16 45.69 49.64 50.37 50.72 78.24 78.75 76.97 100.55 101.10 100.51 136.81 137.30 133.25

3,033,290 2,833,700 2,531,620 7,519,990 6,895,040 4,721,520 16,622,200 12,715,800 8,954,050 19,846,900 20,319,300 10,387,000 29,476,600 29,013,600 22,723,500

237,065 615,530 1,116,720 767,088 994,635 1,299,070 178,477 194,557 404,503 96,539 120,569 381,954 155,406 115,835 214,828

10,801 10,800 10,800 10,803 10,801 10,801 10,804 10,801 10,800 10,806 10,802 10,800 10,808 10,802 10,801

21.61 21.17 20.01 42.18 39.62 28.43 76.40 59.37 43.08 91.14 94.18 49.60 135.36 134.05 106.13

3,808,505 3,778,528 3,583,692 4,510,632 4,480,222 4,301,624 7,828,656 7,931,410 7,584,952 8,463,409 8,387,908 8,142,144 8,897,358 8,891,232 8,720,024

120 62 36 157 86 46 214 116 57 258 138 74 350 170 86

27.13 27.89 27.49 25.30 26.10 26.08 35.98 37.41 36.80 38.87 39.46 39.31 40.86 41.77 41.96

25.56 33.34 41.56 40.02 35.02 8.89 52.90 37.63 15.29 57.36 58.72 21.61 69.82 69.35 61.63

TT-10-400-1-10 TT-10-400-2-10 TT-10-400-4-10 TT-10-600-1-10 TT-10-600-2-10 TT-10-600-4-10 TT-10-800-1-10 TT-10-800-2-10 TT-10-800-4-10 TT-10-1000-1-10 TT-10-1000-2-10 TT-10-1000-4-10 TT-10-1200-1-10 TT-10-1200-2-10 TT-10-1200-4-10

30.11 33.84 33.75 28.38 28.86 28.74 37.08 36.79 36.53 44.33 45.18 45.21 53.80 53.53 53.13

3,983,120 3,292,210 1,285,320 4,816,760 4,314,360 2,498,310 8,006,310 7,689,280 6,324,230 9,627,130 9,578,170 8,321,170 11,534,900 10,365,300 9,601,930

195,417 364,183 534,844 614,681 750,225 998,377 28,024 110,549 178,630 7267.72 117,209 289,559 14,218 97,221 194,502

10,802 10,800 10,800 10,805 10,801 10,800 10,808 10,802 10,801 10,813 10,803 10,800 10,818 10,805 10,801

28.37 24.43 11.14 27.02 25.17 15.99 36.80 36.30 31.01 44.21 44.92 40.14 52.97 48.53 46.00

3,022,055 3,012,880 2,798,428 2,959,429 2,919,640 2,750,996 4,136,948 4,040,168 3,859,848 4,399,797 4,368,718 4,173,512 4,550,840 4,486,944 4,300,716

169 105 50 203 114 60 259 155 81 325 175 93 374 203 94

21.53 22.44 21.91 16.60 17.35 17.40 19.02 19.55 19.71 20.20 21.03 21.13 20.90 21.58 21.71

24.13 8.48 117.72 38.56 32.33 10.11 48.33 47.46 38.97 54.30 54.39 49.84 60.55 56.71 55.21

Average

60.56

9,961,087

379,589

10,803

49.18

5,169,707

147.67

27.15

28.97

TT-10-400-1-5 TT-10-400-2-5 TT-10-400-4-5 TT-10-600-1-5 TT-10-600-2-5 TT-10-600-4-5 TT-10-800-1-5 TT-10-800-2-5 TT-10-800-4-5 TT-10-1000-1-5 TT-10-1000-2-5 TT-10-1000-4-5 TT-10-1200-1-5 TT-10-1200-2-5 TT-10-1200-4-5

ALNS

ALNS vs. Barrena et al. (2014) Improvement (%)

solutions. In our tests, we have executed up to 70,000 iterations, but in Fig. 4 we only show the results corresponding to the first 30,000 iterations since the algorithm has already converged at this point. It can be observed that, as expected, at the beginning of the running process, the algorithm accepts a relatively large number of worsening solutions while it becomes more conservative towards the end, and finally converges to the optimal solution. After executing 20 reheatings, the running time for the presented instance TT-3-200-1-5 was equal to 20 s. In Fig. 4, we show the results of the first run since the algorithm already yields the optimal solution. The running time is then equal to 1 s.

148

E. Barrena et al. / Transportation Research Part B 70 (2014) 134–150

Table 5 Summary of computational results of a regular timetable, of Barrena et al. (2014), and of our ALNS on C5 instances Instance

Regular

Barrena et al. (2014)

AWT

UB1

LB

Time (s)

AWT

ALNS UB2

Time (s)

AWT

ALNS vs. Barrena et al. (2014) Improvement (%)

C5-400-1-5 C5-400-2-5 C5-400-4-5 C5-600-1-5 C5-600-2-5 C5-600-4-5 C5-800-1-5 C5-800-2-5 C5-800-4-5 C5-1000-1-5 C5-1000-2-5 C5-1000-4-5 C5-1200-1-5 C5-1200-2-5 C5-1200-4-5

32.04 32.27 32.17 49.99 50.26 49.87 69.22 69.08 68.03 84.98 85.44 84.63 109.00 108.70 107.73

351,578 366,436 368,016 713,138 712,808 697,884 1,420,300 1,295,520 1,297,920 2,232,260 2,117,930 2,069,540 3,273,840 3,350,500 2,961,000

44,711 62,644 144,576 30,773 62,346 100,063 44,429 69,211 130,709 37,783 73,783 132,404 34,870 67,301 127,146

10,801 10,801 10,801 10,801 10,801 10,800 10,802 10,801 10,801 10,804 10,801 10,801 10,806 10,802 10,800

28.07 28.87 29.33 45.87 44.74 43.94 69.22 61.77 61.57 84.94 78.77 76.83 109.00 108.70 96.02

347,774 365,958 367,672 657,123 685,316 691,816 1,167,838 1,204,986 1,222,216 1,920,051 1,985,946 2,009,336 2,716,630 2,796,192 2,815,152

60 32 17 84 43 23 105 53 28 129 66 34 154 80 44

27.77 28.84 29.31 42.27 43.05 43.58 56.91 57.53 58.09 73.06 73.92 74.66 90.45 90.88 91.39

1.08 0.13 0.09 7.85 3.86 0.87 17.78 6.99 5.83 13.99 6.23 2.91 17.02 16.54 4.93

C5-400-1-10 C5-400-2-10 C5-400-4-10 C5-600-1-10 C5-600-2-10 C5-600-4-10 C5-800-1-10 C5-800-2-10 C5-800-4-10 C5-1000-1-10 C5-1000-2-10 C5-1000-4-10 C5-1200-1-10 C5-1200-2-10 C5-1200-4-10

17.56 18.96 18.98 25.73 26.26 26.30 35.13 35.34 35.44 44.44 44.68 44.79 54.50 54.43 54.63

205,070 221,142 194,616 400,063 396,106 389,040 720,960 732,108 697,328 1,167,870 1,189,540 1,177,340 1,637,000 1,662,340 1,562,580

21,974 38,235 56,767 15,834 38,197 71,623 20,380 32,511 92,559 18,472 30,750 108,988 11,019 22,438 101,524

10,802 10,801 10,800 10,803 10,801 10,801 10,807 10,802 10,801 10,808 10,802 10,801 10,811 10,803 10,801

16.37 17.82 16.45 25.73 25.31 25.38 35.13 35.34 34.00 44.44 44.68 44.57 54.50 54.43 51.62

197,638 200,196 197,020 358,608 372,610 368,204 649,808 667,874 667,392 1,063,456 1,092,802 1,089,628 1,489,795 1,531,472 1,527,444

88 47 26 113 59 32 136 71 38 163 85 45 189 101 54

15.78 16.23 16.63 23.07 23.86 24.13 31.67 32.33 32.63 40.47 41.13 41.40 49.60 50.23 50.50

3.62 9.47 1.24 10.36 5.93 5.36 9.87 8.77 4.29 8.94 8.13 7.45 8.99 7.87 2.25

Average

52.35

1,186,059

61,467

10,802

49.78

1,080,932

45.71

6.87

73.3

Regarding the use of an initial solution, we have run all the tests with and without an initial solution and no significant differences were observed in the results. These findings are typical of other ALNS implementations for different problems and indicate that the algorithm is quite robust since its performance is not affected by the quality of the initial solution. We present in Table 1 a summary of the AWT per passenger and of the running time over all instances, along with a comparison with other algorithms. The rows represent the different instance sizes containing three, six and 10 stations, as well as the real-world based case of the line C5 of Madrid Metropolitan Railway. The columns present the AWT when using a regular timetable, the best upper bound on the AWT obtained by the branch-and-cut algorithm of Barrena et al. (2014), which was truncated after 3 h of computation, and the solutions of our ALNS algorithm. We note that most instances could not be solved to optimality within the time limit used in Barrena et al. (2014), and optimal solutions were obtained only for some instances with three stations. It can be observed that the AWT is considerably improved with respect to a regular timetable, and that our ALNS algorithm can significantly improve the previous best known AWT. Our ALNS implementation can solve the problem within much shorter computation times than the truncated branch-and-cut algorithm. On average, the AWT is reduced by 26.58%, and the ALNS metaheuristic finishes within less than 1% of the running time required by the truncated branchand-cut algorithm, most of the time converging much earlier. In order to compare the results of our ALNS algorithm with those obtained by branch-and-cut (Barrena et al., 2014), we provide in Tables 2–5 detailed results for instances with three, six, and 10 stations, and the line C5 of Madrid Metropolitan Railway, respectively. These tables provide the AWT per passenger for a regular timetable, as well as for the timetables obtained in Barrena et al. (2014) and with our ALNS metaheuristic. We also show the upper bound (UB) on the objective function, and the running time needed by both algorithms, as well as the lower bound obtained in Barrena et al. (2014). In order to compare the two methods, we provide in the last column the Improvement (%) statistic, which represents the improvement of the upper bound UB2 obtained with the ALNS over the upper bound UB1 obtained by Barrena et al. (2014). It is computed as Improvement (%) ¼ 100ðUB1  UB2 Þ=UB1 . Out of 106 open instances without a known optimal solution, we obtain 84 new best known solutions. The results of Tables 2–5 were obtained with 20 reheatings, that is, after finishing the algorithm, the temperature s is reset to its initial value and the process is repeated 20 times. We have also stored the solutions after each reheating and we have observed that the same solutions are often obtained after the first reheating. So, we could have avoided the 20 reheatings and reduced the computational time by a factor of 20 without much impact on the quality of the results.

E. Barrena et al. / Transportation Research Part B 70 (2014) 134–150

149

In Table 2, if we disregard the cases where the optimum is reached (marked with ⁄), the resulting average improvement then becomes 7.33%. If we only consider the instances with known optimal solutions, our ALNS metaheuristic is able to find the optimum in more than 70% of the cases within only a fraction of the computational time used by the exact branch-andcut algorithm. When ALNS does not reach the optimum, the worst optimality gap is only 2.83%. The average results of Tables 2–4 show that, as expected, ALNS is able to obtain better improvements on the larger instances than on the smaller ones. Moreover, we see from Tables 2–5 that ALNS tends to provide better improvements on those instances having a larger planning horizon, a larger number of trains, and smaller time units, i.e., on the most difficult instances. In order to derive a better insight into our results, we depict in Figs. 5 and 6 the AWT improvement obtained with ALNS and with the algorithm of Barrena et al. (2014) over a regular timetable for the cases of six and 10 stations. Once again, it can be observed that the improvement of ALNS over Barrena et al. (2014) is remarkable, especially on the larger instances. These figures also clearly show that there is nearly no difference between the ALNS results obtained with and without an initial solution, except for the smaller instances, where an initial solution provided by a regular timetable enables the ALNS metaheuristic to reach marginally better final solutions. 6. Conclusions We have proposed two non-linear programming formulations for a single-line timetabling problem arising in an RRT system under dynamic demand, with the objective of minimizing the AWT per passenger. Through a suitable analysis of the second formulation, we have derived principles that should guide a heuristic in removing trains from a timetable and in introducing new ones. We have concluded that the choice of an ALNS metaheuristic which effectively combines these two operators together with a simulated annealing acceptance criterion should be highly appropriate for our problem. Our ALNS algorithm maximizes passenger welfare, while respecting all constraints and limitations of the problem. To the best of our knowledge, this is the first time that this type of timetabling problem has been efficiently solved under a dynamic demand behavior. In comparison, the exact algorithm developed in Barrena et al. (2014) can only be used to solve very small instances, and yields very large optimality gaps for medium and large instances. By performing extensive computational experiments on real-world based and randomly generated instances, we have compared our ALNS implementation against a state-of-art truncated branch-and-cut algorithm. We have succeeded in reducing the passenger AWT by 26% by using less than 1% of the computation time required by the latter algorithm. Out of the 120 open instances, we were able to improve the best known solution in 84 cases and to reach the optimum for 10 out of 14 instances with known optimal solutions. These clearly confirm the efficiency and superiority of our ALNS metaheuristic. Acknowledgements This work was partly supported by the Natural Sciences and Engineering Research Council of Canada (NSERC) under grant 39682-10, and by the Excellence Program of the Andalusian Government under grant P09-TEP-5022. This support is gratefully acknowledged. We also thank the Calcul Québec for providing parallel computing facilities. Thanks are due to the Associate Editor and to the referees for their valuable comments. References Abril, M., Salido, M.A., Barber, F., Ingolotti, L., Tormos, P., Lova, A., 2006. Distributed constraint satisfaction problems to model railway scheduling problems. In: Allan, J., Brebbia, C.A., Rumsey, A.F., Scuitto, G., Sone, S., Goodman, C.J. (Eds.), Proceedings of the X International Conference Computers in Railways, Computer System Design and Operation in the Railway and Other Transit Systems. WIT Press, Southampton, UK, pp. 289–297, Prague, Czech Republic. Apostol, T.M., 1957. Mathematical Analysis. Addison-Wesley, Reading, MA. Barrena, E., Canca, D., Coelho, L.C., Laporte, G., 2014. Exact formulations and algorithm for the train timetabling problem with dynamic demand. Computers & Operations Research 44, 66–74. Brännlund, U., Lindberg, P.O., Nou, A., Nilsson, J.E., 1998. Railway timetabling using Lagrangian relaxation. Transportation Science 32 (4), 358–369. Burdett, R.L., Kozan, E., 2010a. A disjunctive graph model and framework for constructing new train schedules. European Journal of Operational Research 200 (1), 85–98. Burdett, R.L., Kozan, E., 2010b. A sequencing approach for creating new train timetables. OR Spectrum 32 (1), 163–193. Cacchiani, V., Toth, P., 2012. Nominal and robust train timetabling problems. European Journal of Operational Research 219 (3), 727–737. Cacchiani, V., Caprara, A., Toth, P., 2008. A column generation approach to train timetabling on a corridor. 4OR 6 (2), 125–142. Cacchiani, V., Caprara, A., Toth, P., 2010a. Non-cyclic train timetabling and comparability graphs. Operations Research Letters 38 (3), 179–184. Cacchiani, V., Caprara, A., Toth, P., 2010b. Solving a real-world train-unit assignment problem. Mathematical Programming 124 (1-2), 207–231. Cai, X., Goh, C.J., Mees, A.I., 1998. Greedy heuristics for rapid scheduling of trains on a single track. IIE Transactions 30 (5), 481–493. Caimi, G., Burkolter, D., Herrmann, T., Chudak, F., Laumanns, M., 2009. Design of a railway scheduling model for dense services. Networks and Spatial Economics 9 (1), 25–46. Canca, D., Barrena, E., Algaba, E., Zarzo, A., 2014. Design and analysis of demand-adapted railway timetables. Journal of Advanced Transportation 48 (2), 119–137. Caprara, A., Fischetti, M., Toth, P., 2002. Modeling and solving the train timetabling problem. Operations Research 50 (5), 851–861. Caprara, A., Monaci, M., Toth, P., Guida, P.L., 2006. A Lagrangian heuristic algorithm for a real-world train timetabling problem. Discrete Applied Mathematics 154 (5), 738–753. Carey, M., 1994. A model and strategy for train pathing with choice of lines, platforms, and routes. Transportation Research Part B 28B (5), 333–353. Ceder, A., 2001. Bus timetables with even passenger loads as opposed to even headways. Transportation Research Record 1760 (1), 3–9. A. Ceder. Public-transport automated timetables using even headway and even passenger load concepts. In: The 32nd Australian Transport Research Forum, pp. 1–17, Auckland, New Zealand, 2009.

150

E. Barrena et al. / Transportation Research Part B 70 (2014) 134–150

Chen, D., Lv, M., Ni, S., 2011. Study on initial schedule optimization model of intercity passenger trains based on ACO algorithm. International Journal of Advancements in Computing Technology 3 (4), 222–228. Chung, J.-W., Oh, S.-M., Choi, I.-C., 2009. A hybrid genetic algorithm for train sequencing in the korean railway. Omega 37 (3), 555–565. Coelho, L.C., Cordeau, J.-F., Laporte, G., 2012. Consistency in multi-vehicle inventory-routing. Transportation Research Part C 24, 270–287. Cordeau, J.-F., Toth, P., Vigo, D., 1998. A survey of optimization models for train routing and scheduling. Transportation Science 32 (4), 380–404. Cordone, R., Redaelli, F., 2011. Optimizing the demand captured by a railway system with a regular timetable. Transportation Research Part B 45 (2), 430– 446. Corman, F., D’Ariano, A., Pacciarelli, D., Pranzo, M., 2010. A tabu search algorithm for rerouting trains during rail operations. Transportation Research Part B 44 (1), 175–192. Daganzo, C.F., 1997. Fundamentals of Transportation and Traffic Operations. Pergamon-Elsevier, Oxford, UK. Furini, F., Kidd, M.P., 2013. A fast heuristic approach for train timetabling in a railway node. Electronic Notes in Discrete Mathematics 41 (5), 205–212. Garey, M.R., Johnson, D.S., 1979. Computers and Intractability: A Guide to the Theory of NP-Completeness. W.H. Freeman & Co., San Francisco. Guihaire, V., Hao, J.K., 2008. Transit network design and scheduling: A global review. Transportation Research Part A 42 (10), 1251–1273. F. Hänseler, B. Farooq, and M. Bierlaire. Preliminary ideas for dynamic estimation of pedestrian origin-destination demand within train stations. In: Proceedings of the Swiss Transport Research Conference, Ascona, Switzerland, 2012. Hurdle, V.F., 1973. Minimum cost schedules for a public transportation route. Transportation Science 7 (2), 109–157. Ingolotti, L., Lova, A., Barber, F., Tormos, P., Salido, M.A., Abril, M., 2006. New heuristics to solve the CSOP railway timetabling problem. In: Ali, M., Dapoigny, R. (Eds.), Advances in Applied Artificial Intelligence, Lecture Notes in Artificial Intelligence, vol. 4031. Springer, Berlin Heidelberg, pp. 400–409. Jamili, A., Shafia, M.A., Sadjadi, S.J., Tavakkoli-Moghaddam, R., 2012. Solving a periodic single-track train timetabling problem by an efficient hybrid algorithm. Engineering Applications of Artificial Intelligence 25 (4), 793–800. Kroon, L.G., Huisman, D., Abbink, E., Fioole, P.-J., Fischetti, M., Maróti, G., Schrijver, A., Steenbeek, A., Ybema, R., 2009. The new Dutch timetable: the OR revolution. Interfaces 39 (1), 6–17. Lee, Y., Chen, C.-Y., 2009. A heuristic for the train pathing and timetabling problem. Transportation Research Part B 43 (8), 837–851. Liebchen, C., 2005. A cut-based heuristic to produce almost feasible periodic railway timetables. In: Nikoletseas, S.E. (Ed.), Experimental and Efficient Algorithms, Lecture Notes in Computer Science, vol. 3503. Springer, Berlin Heidelberg, pp. 354–366. Liebchen, C., Möhring, R., 2002. A case study in periodic timetabling. Electronic Notes in Theoretical Computer Science 66 (6), 1–14. Nachtigall, K., Voget, S., 1996. A genetic algorithm approach to periodic railway synchronization. Computers & Operations Research 23 (5), 453–463. Newell, G.F., 1971. Dispatching policies for a transportation route. Transportation Science 5 (1), 91–105. Niu, H., Zhou, X., 2013. Optimizing urban rail timetable under time-dependent demand and oversaturated conditions. Transportation Research Part C 36, 212–230. Odijk, M.A., 1996. A constraint generation algorithm for the construction of periodic railway timetables. Transportation Research Part B 30 (6), 455–464. Reimann, M., Leal, J.E., 2013. Single line train scheduling with ACO. In: Middendorf, M., Blum, C. (Eds.), Evolutionary Computation in Combinatorial Optimization, Lecture Notes in Computer Science, vol. 7832. Springer, Berlin Heidelberg, pp. 226–237. B. Riemann. Über die Darstellbarkeit einer Function durch eine trigonometrische Reihe. Dieterich, 1867. URL http://eudml.org/doc/203787. Ropke, S., Pisinger, D., 2006. An adaptive large neighborhood search heuristic for the pickup and delivery problem with time windows. Transportation Science 40 (4), 455–472. Serafini, P., Ukovich, W., 1989. A mathematical model for periodic scheduling problems. SIAM Journal on Discrete Mathematics 2 (4), 550–581. Tormos, P., Lova, A., Barber, F., Ingolotti, L., Abril, M., Salido, M.A., 2008. A genetic algorithm for railway scheduling problems. In: Xhafa, F., Abraham, A. (Eds.), Metaheuristics for Scheduling in Industrial and Manufacturing Applications, Studies in Computational Intelligence, vol. 128. Springer, Berlin Heidelberg, pp. 255–276. Vansteenwegen, P., Van Oudheusden, D., 2006. Developing railway timetables which guarantee a better service. European Journal of Operational Research 173 (1), 337–350. Vuchic, V.R., 2005. Urban Transit, Operations, Planning and Economics. Wiley, Hoboken, NJ. Vuchic, V.R., 2007. Urban Transit Systems and Technology. Wiley, Hoboken, NJ. Wong, R.C.W., Yuen, T.W.Y., Fung, K.W., Leung, J.M.Y., 2008. Optimizing timetable synchronization for rail mass transit. Transportation Science 42 (1), 57– 69. Yano, C.A., Newman, A.M., 2001. Scheduling trains and containers with due dates and dynamic arrivals. Transportation Science 35 (2), 181–191. Yuan, Z., Fügenschuh, A., Homfeld, H., Balaprakash, P., Stützle, T., Schoch, M., 2008. Iterated greedy algorithms for a real-world cyclic train scheduling problem. In: Blesa, M.J., Blum, C., Cotta, C., Fernández, A.J., Gallardo, J.E., Roli, A., Sampels, M. (Eds.), Hybrid Metaheuristics, Lecture Notes in Computer Science, vol. 5296. Springer, Berlin Heidelberg, pp. 102–116. Zhou, X., Zhong, M., 2007. Single-track train timetabling with guaranteed optimality: branch-and-bound algorithms with enhanced lower bounds. Transportation Research Part B 41 (3), 320–341.