A combined simulation-optimization approach for minimizing travel time and delays in railway timetables

Transportation Research Part B 126 (2019) 192–212

Contents lists available at ScienceDirect

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

A combined simulation-optimization approach for minimizing travel time and delays in railway timetables Johan Högdahl∗, Markus Bohlin, Oskar Fröidh Department of Civil and Architectural Engineering KTH Royal Institute of Technology SE-100 44, Stockholm, Sweden

a r t i c l e

i n f o

Article history: Received 18 September 2018 Revised 27 February 2019 Accepted 8 April 2019 Available online 17 June 2019 Keywords: Railroad Robustness Optimization Simulation Punctuality

a b s t r a c t Minimal travel time and maximal reliability are two of the most important properties of a railway transportation service. This paper considers the problem of finding a timetable for a given set of departures that minimizes the weighted sum of scheduled travel time and expected delay, thereby capturing these two important socio-economic properties of a timetable. To accurately represent the complex secondary delays in operational railway traffic, an approach combining microscopic simulation and macroscopic timetable optimization is proposed. To predict the expected delay in the macroscopic timetable, a surrogate function is formulated, as well as a subproblem to calibrate the parameters in the model. In a set of computational experiments, the approach increased the socio-economic benefit by 2–5% and improved the punctuality by 8–25%. © 2019 Elsevier Ltd. All rights reserved.

1. Introduction To make the best use of the limited resources in society, existing and planned railway transportation services need to be evaluated and compared by policy makers and infrastructure managers, typically using socio-economic evaluation. This evaluation, by means of a cost benefit analysis (CBA), has the aim to optimize welfare in society in a manageable size, defined by the project or problem studied. For railway services in general, a demand of both passenger and freight services is met by a preplanned supply described as train paths in the timetable. In many aspects, the timetable can be seen as the foundation of railway operations, and in comparison to infrastructure improvements, improved planning of the services give benefits in the short term to a low cost since capacity utilization is managed towards improved efficiency, or increased welfare. Furthermore, socio-economic value is important as it serves as the basis for strategic investment decisions and prioritizations. A central term in most socio-economic timetable calculations is reliability, including the passenger’s and freight transport users value of reliability, and the reliability of the service itself (i.e. the consistency of the planned timetable when operated). The topic of passengers’ value of reliability has been investigated in several works. For example, Bates et al. (2001) and Carrion and Levinson (2012) gather findings about travellers’ value of reliability for both public and private transport and present an overview of valuations, and Concas and Kolpakov (2009) conclude that especially the effect of changes requires further research. De Jong et al. (2007), Tseng (2008), and Chang (2010) all show the importance of including reliability in socio-economic analysis. Further, Parbo et al. (2016) state that reliability is more important than travel time for services with low frequency. Altogether, there is a clear need to explicitly include reliability in any evaluation of timetable quality, and hence in any approach for timetable planning. ∗

Corresponding author. E-mail address: [email protected] (J. Högdahl).

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

J. Högdahl, M. Bohlin and O. Fröidh / Transportation Research Part B 126 (2019) 192–212

193

The method in this paper is centered around a simple idea, which is: when adjusting a timetable to improve punctuality, events should be assigned time supplements proportionally to their expected delays. This is formalized into a method that, for a given initial timetable, combine micro-simulation (to estimate the expected delays) with macroscopic mixed integer linear programming to reschedule it, such that the weighted sum of scheduled travel time and expected delay is minimized given that the sequence of the trains is considered to be fixed. The benefit of the combined approach is that the effects of random disturbances on a microscopic level is included in the planning problem on the macroscopic level, and by keeping the timetable planning on a macroscopic level and considering the sequence of the trains as fixed it is possible to apply the approach on larger networks. 1.1. Related work The current state of practice in railway timetable analysis is micro-simulation, which main advantage is the ability to handle complex and detailed models, while methods in the scientific community for timetable planning (often referred to as the train timetabling problem) are most often based on macro-level linear and integer programming. In the following section an overview of the scientific literature in the area is given. Railway timetabling and robust railway timetabling have been an active field in railway operations research for the last decade and recent surveys on railway timetabling models have been conducted by Harrod (2012), Cacchiani and Toth (2012), Caimi et al. (2017), where the two latter also includes robust timetabling. Various robustness definitions have also been surveyed by Andersson (2014, p. 11) and Dewilde (2014, p. 35) and both authors emphasize the importance of margins in the timetable for handling small everyday-delays that occur in rail operations. Stochastic programming has been utilized by Kroon et al. (2008), who describe a two-stage stochastic optimization model with recourse to improve the robustness of a periodic timetable. They assume that an operating timetable is known, including fixed routing and platforming at stations which remain unchanged, and combine a timetabling model with a recourse model to minimize the sum of weighted delays. A related approach is light robustness, which is based on the principles of robust programming and aims to generate solutions that are robust against uncertain parameters, while at the same time assuring that the deterioration of the objective function value is bounded. Light robustness was introduced as a general method and applied to timetabling by Fischetti and Monaci (2009). In their paper, they found that the method achieves similar robustness as stochastic programming but in shorter time, and concluded that light robustness is more suitable than stochastic programming to solve large-scale scenarios. Cacchiani et al. (2012) later proposed a Lagrangian heuristic method for generation of robust solutions, based on iteratively solving a relaxed problem of the train timetabling problem (TTP), with additional parameters that controls robustness in combination with a scheme for dynamically increasing their weight in the objective. In a set of experiments they showed results of comparable quality with the solutions obtained by light robustness, but with much shorter computing times, which prove effectiveness of the approach. Schöbel and Kratz (2009) and Schlechte and Borndörfer (2010) both considered the robustness problem as a multiobjective optimization problem, in which efficiency and robustness are two conflicting objectives to be optimized. Schöbel and Kratz present three different delay management strategies, propose corresponding robustness measures for each strategy, and derive necessary conditions for the pareto-optimal solutions. Schlechte and Borndörfer propose a robustness function based on buffer time allocation and to scalarize the multi-objective problem, which is solved by an approach based on column generation, to yield pareto-optimal solutions. Chow and Pavlides (2018) consider timetable optimization from the customers perspective and propose a stochastic multi-objective optimization problem which take journey times, waiting times, punctuality, and crowdedness of trains as objective functions. The optimization problem is scalarized and solved with an iterative heuristic algorithm which in each iteration generates a set of candidate timetables (which are evaluated with Monte Carlo simulation) using a genetic algorithm together with Dijkstra’s algorithm. Sels et al. (2016) also consider timetable optimization from a customer (passenger) perspective and formulate the problem as to minimize the total expected passenger time (including the effects of primary and secondary delays as expected values of stochastic variables). By including both the travel time and the expected delay in the objective function, the need for artificial upper limits on the time supplements is solved, which solved infeasability issues of earlier attempts of robust PESP scheduling. There have also been some attempts to increase robustness by distributing timetable margins without directly considering simulated or actual performance. Khoshniyat and Peterson (2017) propose such an approach based on the observation that on average delays are proportional to the traveled distance. In this work, a strategy is proposed where the minimum headway increases along with each train’s journey. A related approach is the concept of Robustness in Critical Points (RCP), proposed by Andersson et al. (2013). Critical points refer to time-sensitive dependencies between pairs of trains, for instance meetings or overtakings. The RCP measure for each time-space point is then defined for two trains in order as the sum of runtime margin for the first train, runtime margin for the second train and the headway margin between the trains. A mixed integer linear programming model was proposed to redistribute the margin time in an existing timetable. Applied to a realistic example the model achieved a 20% reduction in total delay at the final destination compared to a non-optimized timetable (Solinen et al., 2017). However, it was also observed that the performance of some trains did not improve (or

194

J. Högdahl, M. Bohlin and O. Fröidh / Transportation Research Part B 126 (2019) 192–212

even became worse) by increasing the RCP. Jovanovic´ et al. (2017) extended the operational criteria of RCP with a commercial criterion and modeled the problem of redistributing buffer times as a multi-dimensional knapsack problem. A common characteristic with the above-mentioned methods is that none include the effects of stochasticity explicitly in the model. Approaches that combine optimization with detailed micro-level simulation have received less attention than other approaches. Lee et al. (2017) described the robustness problem as a trade-off between efficient travel times and small delays. They proposed an iterative heuristic method that is based on essentially three steps: (1) solve a LP-problem to allocate time supplements and buffer times such that the deviation from the desired target value of each activity is minimized; (2) simulate the obtained timetable; (3) update the target values of the LP-problem based on the simulation output. The output after iterating the three-step approach is a timetable for one value of a global parameter which determines how large time supplements that are allowed without generating a penalty in the objective. By varying the value of the global parameter a set of alternative timetables is generated, and among those a timetable with a desired balance between time supplements and average delays can be chosen. Burggraeve and Vansteenwegen (2017) focus on robust timetabling of complex station areas which is addressed using an iterative four-step approach consisting of the following steps: (1) find an optimal route plan for each train such that the traffic is spread evenly at the station area; (2) timetable each trains route by maximizing the weighted minimum buffer time between each train; (3) simulate the obtained timetable; (4) calculate new dwell and running times such that a predetermined percentile of the train movements, with respect to the simulated realizations of the timetable, is expected to be on-time. Repeat the steps until the approach converge or until either the time or the number of iterations run out. Contrary to the approach in this paper, the expected delay is not directly minimized in the optimization model. Instead, as in Lee et al. (2017), a heuristic search of the time supplement and buffer time-space is conducted, which give a set of alternative timetables from which the best one is chosen. 1.2. Relation to previous work Table 1 summarizes the relation between this paper and previous work. In contrast to existing approaches, the model in this paper minimize the weighted sum of scheduled travel time and the predicted value of the expected delay. This has the benefit that the use of time supplements is distributed more efficiently, so that trains that operate with large average delays get larger supplements compared with trains that operate with small average delays. From a model perspective it also has the benefit (as in Sels et al. (2016)) that artificial constraints on running and dwell times are not necessary (which might lead to issues with infeasibility, Sels et al. (2016)) and, simultaneously, that inefficient use of time supplements is avoided. Compared with Burggraeve and Vansteenwegen (2017) and Lee et al. (2017), who propose approaches that combine micro-simulation with optimization, a prediction of the expected delay is included in the objective function of the optimization model, and the expected delay time can therefore be minimized directly. This has the benefit that a separate search for time supplements and buffer times (either as target or actual values) does not have to be conducted. In our approach, the parameters in the objective function are instead calibrated by solving a separate problem off-line. The objective function in this paper is close to the one used by Sels et al. (2016), who minimize the total expected passenger journey time (which is the sum of scheduled journey time and expected journey delay time). However, it differs in

Table 1 Relation between this paper and previous work. Objective = Objective value; Interaction = Knock-on interaction model between trains in the planning model; Order = Train order; Sim. scope = Scope of the simulation; Traffic = Train traffic type considered; Topology = Topology of the infrastructure (line, network or a station area); RCP = Robustness in Critical Points; WD = Weighted delays; VTT1 = Value of travel time including journey time and delay time; VTT2 = Value of travel time including journey time, waiting time, punctuality, and crowdedness on trains; ETT = Expected travel time; HW = Headway; ESD = Empirical secondary delays; XSD = Exponential delays; Dev. = Deviation from ideal timetable; TTS = Total time supplements; WBT = Max-min weighted buffer time; WER = Weighted sum of efficiency and estimated robustness; MO = multi-objective optimization of efficiency and robustness measure. Optimization model

Evaluation

Approach

Objective

Interaction

Order

Sim. scope

Traffic

Topology

Solinen et al. (2017) Jovanovic´ et al. (2017) Khoshniyat and Peterson (2017) Fischetti and Monaci (2009) Kroon et al. (2008) Cacchiani et al. (2012) Schlechte and Borndörfer (2010) Schöbel and Kratz (2009) Sels et al. (2016) Chow and Pavlides (2018) Lee et al. (2017) Burggraeve and Vansteenwegen (2017) This paper

RCP RCP Dev. Dev. WD WER WER MO ETT VTT2 TTS + Dev. WBT VTT1

HW HW HW HW HW HW HW HW HW + XSD HW HW HW HW + ESD

Fixed Fixed Flex. Fixed Fixed Fixed Flex. Flex. Flex. Flex Fixed Flex. Fixed

Macro Macro Micro Micro Micro

Mix Mix Mix Mix P Mix Mix P Mix Mix P Mix

Line Line Line Line Network Line Line Network Line Line Station Line

J. Högdahl, M. Bohlin and O. Fröidh / Transportation Research Part B 126 (2019) 192–212

195

an important way since the objective in this paper is to minimize the weighted sum of scheduled travel time and expected delay time (where the weighting parameter correspond to the valuation of decreasing the expected delay time compared to decreasing the scheduled in-vehicle travel time). This implies that the timetable generated by the approach in this paper may have non-minimal expected travel time (depending on the weighting parameter) but optimal scheduled travel time with respect to the valuation of reducing delays compared to reducing the scheduled travel time (i.e. the balance between scheduled in-vehicle travel time and average delay is socioeconomically optimal). Moreover, the approach in our paper is different from the approach of Sels et al. (2016), since they rely solely on optimization. Finally, using a micro-simulation approach opens up for evaluation of the effect of different infrastructure layouts, signal systems and dispatching strategies. An evaluation of this type is however outside the scope of this paper. 1.3. Contributions The contributions of the paper are the following. 1. The train timetabling problem with scheduled travel time and passenger and freight delay minimization using observed arrival times is formally defined (Section 2.1). In particular, the model includes both the number of passengers and knock-on effects (Section 2.1.2). 2. The stochastic train timetabling calibration problem, which minimizes the root mean square error of the predicted average delay for each event in an optimal timetable solution, is formally defined (Section 2.2). 3. The two optimization problems are combined with micro-simulation to improve the robustness of a non-periodic timetable, for which the sequence of the trains is considered as fixed. In the resulting method, actual arrival times of a timetable, obtained from micro-simulation, are used as input to the timetabling problem, which results in a timetable that minimizes a weighted sum of travel and expected delay time (Section 2.3). 4. A thorough calibration and combined optimization-simulation study on a 383 km section of the Western Main Line in Sweden is used to evaluate the method. The results of the study showed that socio-economic benefits and punctuality could be increased, in particular for the case where the initial punctuality was comparable to the real-world punctuality on the line (Section 3). 1.4. Paper outline This paper consists of four parts. In the next section we propose the train timetabling problem of travel time and delay minimization and also the calibration problem; in the third section we describe the computational experiments and the obtained results; and finally, in the fourth section the conclusions of this article and future work are stated. Distributions of the entry delays and dwell time extensions that have been used in the experiments, and also the settings that have been used with the micro-simulation tool Railsys (Bendfeldt et al., 20 0 0), are given in the appendix. 2. A simulation-based approach for timetable optimization The nominal train timetabling problem (TTP) can be described as follows: Consider a railway network that consists of a number of stations connected with tracks, and a set of railway lines that we wish to timetable. The problem is then to find a timetable such that it satisfies the operational constraints and utilizes the infrastructure as efficient as possible. A practical issue with the solution to the nominal TTP is that it may not yield a timetable which is robust against delays, and for this reason additional time supplements and buffer times have to be included in the timetable. 2.1. The timetabling problem of travel time and delay minimization To address the issue of how supplements should be included and distributed in the timetable we propose to extend the nominal TTP into a new problem such that the overall aim is to find the timetable which minimizes the following two conflicting objectives: scheduled travel time (i.e. planned capacity utilization) and expected delay time (i.e. expected additional capacity utilization). This is formally defined as the timetabling problem of travel time and delay minimization (TTD), which we, in its scalar form, formulate as follows:

minimize subject to

f (t ) = α F (t ) + (1 − α )G(t ) t ∈ X (x ),

(TTD)

where t is the timetable, x is the train sequence, F(t) is the scheduled travel time, G(t) is the total predicted average delay (which ideally is equal to the total expected delay), α is the delay valuation coefficient (i.e. the valuation of decreasing the average delay time compared to decreasing the scheduled in-vehicle travel time), and X(x) is the set of all feasible timetables for the train sequence x. In this article we will study a special case of the TTD, namely the problem of improving an initial timetable t¯ given that the sequence of the trains is considered to be fixed. This problem is called the timetabling problem of travel time and delay

196

J. Högdahl, M. Bohlin and O. Fröidh / Transportation Research Part B 126 (2019) 192–212

minimization with fixed sequence of the trains (TTD-fixed), and is formulated as follows:

minimize subject to

f (t ) = α F (t ) + (1 − α )G(t t ∈ X (x¯ ),

| x¯ )

(TTD-fix)

where x¯ is the fixed sequence of the trains as given in the initial timetable t¯, and G(t | x¯ ) is the total predicted average delay. 2.1.1. Operational constraints of a one-direction single-track line We will focus on a one-direction single-track line (i.e. one of the tracks on a double-track line), connecting two endpoint stations, with a number of intermediate stations. Several models exist to formulate this as an optimization problem and we base our model on the TTP-formulation used by Fischetti et al. (2009). Let S be the set of stations, As be the set of arrival events for station s ∈ S, Ds be the set of departure events for station s ∈ S, T be the set of trains, and tih be the time of event i ∈ [1, Nh − 1] for train h, where Nh is the last event for train h and the events are sorted such that tih ≤ tih+1 . The operational constraints that must be considered are the following: 1. The time difference between two consecutive events of the same train must be larger than either the minimum running time or the minimum dwell time depending on the events, i.e. h tih+1 − tih ≥ di,i +1 ,

i = 1, . . . , Nh − 1,

∀h ∈ T ;

2. Two trains arriving at the same station must be separated by a minimum arrival headway

|

tih

− t kj

|≥

∀

a s,

tih , t kj

∈ As ,

(1)

as ,

i.e.

∀s ∈ S;

(2)

3. Two trains departing from the same station must be separated by a minimum departure headway ds , i.e.

|tih − t kj | ≥ ds , ∀tih , t kj ∈ Ds , ∀s ∈ S;

(3)

4. If one train departs earlier than another train from one station, then the first train must also arrive earlier than the latter train at the following station, i.e.

∀tih , t kj ∈ Ds , ∀s ∈ S,

xh,k = xh,k , i, j i+1, j+1

(4)

where xh,k is a binary variable which is equal to 1 if and only if event i for train h occurs before event j for train k. i, j Due to the absolute values in Eqs. (2) and (3), the constraints above will not give a MIP-formulation when included in an optimization model. To overcome this, Fischetti et al. (2009) simulate the effect of the absolute values by utilizing binary variables and Big-M constraints, which is a common technique for linearizing absolute value constraints. By this approach we obtain

⎧h k h,k ⎨ti − t j ≥  − Mxi, j h k k h |ti − t j | ≥  ⇔ t j − ti ≥  − Mxk,h j,i ⎩ h,k k,h

(5)

xi, j + x j,i = 1,

where M is a sufficiently large constant. 2.1.2. The objective function of TTD-fix The objective function consists of two parts, the scheduled travel time F(t) and total the predicted average delay G(t | x¯ ), and they will be treated separately in this section. The scheduled travel time is calculated as the sum of the travel time multiplied by a weighting parameter (e.g. the number of passengers) between each OD-pair for all trains, i.e:

F (t ) =





  λhi, j t hj − tih ,

(6)

h∈T (i, j )∈ODh

where (i, j) is a trip from departure event i to arrival event j, ODh is the set of all trips with train h (i.e. all its OD-pairs), and λhi, j is the weight (e.g. number of passengers) of trip i to j with train h. When calculating the predicted average delay of the timetable t it is assumed that K observations yr (with r = 1, . . . , K being the observation number) of the initial timetable t¯ is available and that the departure time from the first station is unchanged. Then, for the optimization problem with a fixed sequence of trains, TTD-fix, we propose to model the total predicted average delay G(t | x¯ ) as follows below. First, let the total predicted average delay be defined as the sum of the predicted average delay for each train, i.e:

G(t

| x¯ ) =



gh (t

| x¯ )

h∈T

where gh (t | x¯ ) is the predicted average delay for each train.

(7)

J. Högdahl, M. Bohlin and O. Fröidh / Transportation Research Part B 126 (2019) 192–212

197

Similarly, let the predicted average delay of train h be the sum of the predicted average delay multiplied with a weighting parameter (e.g. number of disembarking passengers) at each arrival event j, i.e:

gh (t

| x¯ ) =



λhj ghj (t | x¯ ),

(8)

j∈Eah

where Eah is the set of all arrival events for train h, and λhj and ghj (t | x¯ ) is a weighting parameter (e.g. number of passengers that disembark at that event) and the predicted average delay, respectively, at event j for train h. Next, assume that a train driver always try to operate the train at the highest possible speed. Then a small change in the timetable (e.g. adjusting the time supplement up and to event j for train h or adjusting the headway to the previous train) does not affect the train driver’s behavior and thus, does not change the actual arrival time, but only the average deviation from the timetable. Therefore, the predicted average delay at the arrival event j for train h is calculated as the sum of the predicted new delay (i.e. the predicted delay after adjusting the arrival time) multiplied with the probability (estimated by simulating the initial timetable) that this new delay occur, i.e:

ghj (t

| x¯ ) =

 ω ∈

phj,ω mhj,ω (t

| x¯ ),

(9)

h j

where hj is the set of observed delays (w.r.t. the initial timetable t¯) at event j for train h, phj,ω is the empirical probability that event j of train h is delayed by ω, and mhj,ω (t | x¯ ) is the predicted new delay of event j for train h (i.e. given that the new timetable t would have been used instead of the initial timetable t¯ in the situation where the delay ω initially was observed). In Eq. (9), the empirical probability phj,ω is computed as the number of times the delay ω was observed divided with the total number of observations, i.e:

phj,ω =

K  1  h 1ω y j,r − t¯hj , K

(10)

r=1

where yhj,r denotes the observed time of event j for train h in the r:th replication of the simulation of the initial timetable

t¯, K is the number of replications, and 1ω (· ) is the indicator function, which is 1 if and only if yhj,r − t¯hj = ω, else 0. Regarding the new delay function mhj,ω (t | x¯ ) in Eq. (9), assume that: (1) if the new arrival time t hj is increased, then the average delay is expected to decrease and vice versa; (2) if the headway to the proceeding train is decreased, then the average delay is assumed to increase (due to greater risk of knock-on delays) and vice versa; and (3) a driver does not change driving behavior due to the advertised arrival time. The new delay is then computed as follows:

mhj,ω (t

     | x¯ ) = max 0, t¯hj + ω − t hj + βφ t hj , tlk ,

(11)

where β is the knock-on gain and φ (t hj , tlk ) is called the train dependency penalty function which penalizes decreased buffer times between adjacent trains, and is defined as follows:

φ



t hj , tlk





=

tlk − t¯lk , 0,

if t hj − tlk < b else,

(12)

where event l for train k directly precedes event j of train h, and b is the train independence threshold (i.e. when the headway between two trains is so large that it is unlikely that the following train receives a knock-on delay from the leading train). An objection to the above calculation of the new delay is that if one train travels between stations A and B and if the scheduled departure time from A is increased (i.e. pushed later in time compared with the initial timetable) and the scheduled arrival time at station B is unchanged, then intuitively the delay of arrival at B should increase - which it does not in the above model. However, this case would also indicate that the initial time supplement between A and B is used to reduce a delay that occurred before the departure from A and that if the departure delay would have been avoided in the first place (e.g. by pushing the departure from A later in time), then the time supplement between A and B would have remained unused, and hence it should be possible to reduce it. However, as the proposed model clearly does not capture the effects of reducing time supplements it is of interest to investigate alternative models in future research. Finally, by combining Eqs. (7)–(9) and (11) the total predicted average delay is expressed as

G(t

| x) =

  h∈T j∈Eah ω∈h j



  h λhj phj,ω max 0, t j + ω − t hj − βφ t hj , tlk .

This function and its resulting value will be referred to as the predicted average delay later in this article.

(13)

198

J. Högdahl, M. Bohlin and O. Fröidh / Transportation Research Part B 126 (2019) 192–212

2.1.2.1. Linearization. The predicted average delay G(t | x¯ ), as expressed in Eq. (13), is not a linear function and will not yield a mixed integer linear programming formulation when used in TTD-fix. To obtain a MIP we must replace the max -function and φih (tih , t kj ) in the objective with linear expressions.

First, we linearize the function φ (t hj , tlk ). This is done by introducing the binary variable zhj , which is equal to 1 if t hj − tlk <

b and otherwise 0, and a state variable h,k . By utilizing zhj and h,k we can now write φ (t hj , tlk ) as j,l j,l

    φ t hj , tlk = tlk − t¯lk zhj = h,k . j,l

Let

h,k j,l

(14)

replace φ (t hj , tlk ) in Eq. (13), and to simulate the if-else statement in Eq. (12), include the following set of constraints

in the problem TTD-fix:

t hj − tlk − b ≥ −Mzhj   t hj − tlk − (b − ) ≤ M 1 − zhj

h,k ≥ −Mzhj j,l h,k  j,l ≤ Mzhj     h,k ≥ tlk − t¯lk − M 1 − zhj j,l     h,k ≤ tlk − t¯lk + M 1 − zhj , j,l

(15)

where ≥ 0 is a small constant simulating the strict inequality in Eq. (12) and M is a sufficiently large constant. We now have the following expression for the predicted average delay:

G(t

| x¯ ) =

 

h∈T j∈Eah ω∈h j

   λhj phj,ω max 0, t¯hj + ω − t hj + β h,k . j,l

(16)

Next, to obtain a linear formulation of Eq. (16) the max -function must be replaced with linear expressions. To begin, we note that



minimize

f (x ) = max {0, x} ⇔

minimize s. t.

y x y

≤ ≥

y 0.

By introducing yhj,ω ≥ 0 and together with Eq. (15) add





≤ yhj,ω t¯hj + ω − t hj + β h,k j,l

(17)

to the constraints of problem TTD-fix, we can thus write the predicted average delay G(t | x¯ ) as

G(t

| x¯ ) =

  h∈T

j∈Eah

ω ∈

λhj phj,ω yhj,ω .

(18)

h j

2.1.3. Full timetabling model By including Eq. (1) together with the linearized expressions of Eqs. (2) and (3) by using Eq. (5) and having xh,k = x¯h,k j,l j,l (i.e. to fix the sequence of the trains to the sequence in the initial timetable, which make Eq. (4) redundant) and t1h = t¯1h (i.e. enforcing that the initial departure time of each train is unchanged), the operational constraints is included. By having F(t) and G(t | x¯ ) as formulated in Eqs. (6) and (18), respectively, and including Eqs. (15) and (17), we obtain a complete MILP formulation of TTD-fix such that for a given timetable on a one-direction single-track line the weighted sum of scheduled travel time and predicted average delay (which ideally is equal or close to the expected value) is minimized. This complete formulation of the linearized TTD-fix is presented in Appendix A. 2.2. Calibration of model parameters The predicted average delay depend on the parameters b and β , which have to be calibrated. For this reason a separate calibration problem is formulated as the following optimization problem:

minimize subject to

L(b, β | α ) b, β ≥ 0.

(PC)

where b and β are the train independence threshold and the knock-on gain, respectively, and α is the delay valuation coefficient (pre-specified by the user). The objective function is defined as the root mean square error (RMSE) of the predicted average delay per event with respect to the observed average delay. This is calculated as follows:



L(b,

β | α) =

 1 |Ea |

(h, j )∈Ea



Hˆ hj (α

| b, β ) − H hj t(∗α

2 | b, β )

,

(19)

J. Högdahl, M. Bohlin and O. Fröidh / Transportation Research Part B 126 (2019) 192–212

199

Table 2 A random search method to approximately solve the calibration problem. Algorithm 1: A random search method to approximately solve the calibration problem PC Select the initial timetable t¯, the sample space , the number of iterations N, the delay valuation coefficient α , and set the minimal objective function value Lmin = ∞. Simulate the initial timetable t¯ to generate observed arrival and departure times, y¯ , of the initial timetable. Randomly choose b and β from , generate t(α | b, β ) by solving TTD-fix with parameters α , b, and β , and simulate the obtained timetable t(α | b, β ) . Evaluate the objective function L(b, β | α ) according to equation (19). If L(b, β | α ) < Lmin , then update Lmin := L(b, β | α ) and the currently best parameter values b∗ := b and β ∗ := β . Terminate if the procedure has been repeated N times, else return to Step 2.

Step 0: Step 1: Step 2: Step 3: Step 4:

Start Select an initial timetable and a delay valuation coefficient α

Simulate the initial timetable N times

Calibrate the predicted average delay model by solving the calibration problem (PC)

Solve the timetabling problem of travel time and delay minimization with fixed sequencing (TTD-fix) Optimized timetable

Fig. 1. An overview of the combined simulation-optimization timetable improvement approach.

where Hˆ hj (α | b, β ) is the predicted average delay of event (h, j) obtained by solving TTD-fix with the parameters α , β , and b, and where H hj (t(∗α | b, β ) ) is the observed average delay of event (h, j) obtained by simulating optimized timetable t(∗α | (which is the optimal solution of TTD-fix solved with the parameters α , b and β ). The observed average delay of event (h, j) is computed as



H hj t(∗α |

 b, β )

=

 R  1 max 0, yhj,r − t(∗α | R



h b, β ) j

,

b, β )

(20)

r=1

where R is the number of replications of the simulation of timetable t(∗α |



replication r, and

t(∗α |

h

b, β )

b, β )

, yhr, j is the observed time of event (h, j) in

is the scheduled time of event (h, j). j

The calibration problem objective function is stochastic and requires simulation in order to be evaluated. Therefore, a random search method (see for instance Spall (2005)) is applied to approximately solve the calibration problem, described in Table 2. 2.3. The combined simulation-optimization approach for minimizing travel time and delays in timetables with fixed train order By combining the calibration problem with the TTD-fix an optimized timetable is obtained by using the following approach (also illustrated in Fig. 1): 1. 2. 3. 4.

Select an initial timetable and a delay valuation coefficient α . Simulate the initial timetable N times. Decide the model parameters of the predicted average delay model by solving the calibration problem PC, see Table 2. Generate an optimized timetable by solving the optimization problem TTD-fix using the arrival times obtained in step 2 as observations of the initial timetable.

3. Computational experiments To demonstrate the performance of the proposed approach, a series of computational experiments have been conducted, which together with, the experiment setup, are described in this section.

200

J. Högdahl, M. Bohlin and O. Fröidh / Transportation Research Part B 126 (2019) 192–212

3.1. Experiment setup The experiments have been conducted by optimization and simulation of rail operations on the section between Gnesta (63 km southwest of Stockholm) and Partille (near Gothenburg) on the Swedish Western Main Line (which is one of Sweden’s most heavily utilized lines and connects the capital Stockholm with the second-largest city Gothenburg). The section between Gnesta and Partille is an electrified double-track line of 383 km and has 53 stations, including passing tracks. By this choice, the multiple track sections and the Stockholm commuter services on the Western Main Line are excluded. The delay valuation coefficient α in the formulation of TTD and TTD-fix has an important meaning in the model, as it corresponds to the delay time factor (i.e. the value of decreasing the average delay time compared with decreasing the scheduled in-vehicle travel time) used when computing the value of travel time in socio-economic analyses. Table 3 show a compilation of delay time factors used in different countries (Warg and Bohlin, 2016, p. 257). To show the general applicability of TTD-fix, the experiment will be repeated for each delay time factor in Table 3. The remainder of this section will describe the setup of the initial timetable, the simulation model and the optimization model. 3.1.1. Initial timetable The initial timetable is based on the national Swedish timetable of 2015, from which we have selected all southbound trains that operated on any section of the Swedish Western Main Line between Gnesta and Partille on Wednesday, September 23, from 05:00 to 12:00. This selection included 87 trains, which can be divided into the following categories: • • • •

9 freight trains (varying train weight, 70–100 km/h top speed). 17 commuter trains (frequent stops, Alingsås-Gothenburg only, 140 km/h top speed). 40 Intercity/regional trains (some stops, 160–200 km/h top speed). 21 fast trains (interregional, none or few stops, 200 km/h top speed).

The timetable of the selected trains was not conflict-free, since it contained some trains with too small headways (compared to the recommended values specified by the Swedish Transport Administration, see the end of Section 3.1.3) and some trains with infeasible running times, according to the simulation software Railsys, which has been used in the experiments. Therefore, in order to make a fair comparison with the optimized timetables (i.e. the solution of TTD-fix), it was necessary to construct a feasible initial timetable. This was carried out by first generating a macroscopically feasible timetable by solving the TTP (using the 87 trains in the national timetable as a temporary initial timetable) with the following principles: 1. The running time between two stations shall be the largest value of the minimum technical running time, which is calculated by Railsys, and the originally scheduled running time in the national timetable. 2. The dwell times shall be equal to the scheduled dwell times in the national timetable. 3. The headway shall be larger than the minimum headways according to the Swedish Transport Administration’s recommendations (see the end of Section 3.1.3). 4. The objective is to minimize the translation in time of each train. In order to arrive at the initial timetable used in the experiments, the next step was to solve the remaining platforming conflicts and then check that the modified timetable still complied with principles 1–3 above. Fig. 2 shows the graphical timetable of the generated initial timetable, in which we see that we have mixture a of trains operating on different parts of the line.

Table 3 Delay time factors (i.e. the value of decreasing the average delay time compared with decreasing the scheduled in-vehicle travel time) that are used in some countries, together with the corresponding value of α (i.e. the weighting parameter in (TTD)). Data source: (Warg and Bohlin, 2016, p. 257). Delay time factor

Corresponding α -value

Countries

3.7 3.5 3.0 2.8 2.7 2.1 2.0

0.213 0.222 0.250 0.263 0.270 0.303 0.333

New South Wales, Australia Sweden Denmark, UKa Norway ( < 50km) Germany Norway, ( > 50km) Denmarkb , UK

a b

National standard value of the United Kingdom National standard value of Denmark

J. Högdahl, M. Bohlin and O. Fröidh / Transportation Research Part B 126 (2019) 192–212

201

Fig. 2. Graphical timetable for the initial timetable. The trains are color coded as follows: solid green is for commuter trains; dash-dot blue is for Intercity/regional trains; dash red is for fast trains; and solid black is for freight trains. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

3.1.2. Simulation model setup The simulations have been carried out in Railsys version 9.8.25 (Bendfeldt et al., 20 0 0), which is a commercial software developed by Rail Management Consultants GmbH1 used to model, simulate, and evaluate railway traffic. It is the standard simulation tool used by the capacity planners at the Swedish Transport Administration, and they have provided us with the complete simulation model for the Swedish railway, including the national timetable of 2015. The simulations have been carried out such that traffic on branch lines, i.e. trains that enters and leaves the main line at the same station, is excluded, meaning that only the trains in the initial timetable (see Section 3.1.1) are simulated. The assumption simplifies the computational experiments and can be considered justifiable as the branch line traffic has lower priority, and its inclusion would therefore only impact the analysis marginally. In order to guarantee that all trains arrive at their end station, the simulation period has been set from 05.00 to 18:00 (although the last train is scheduled to arrive at its terminal station at 14:47). In each simulation cycle, randomly generated primary delays are added to the entry, dwell and running times of each train. The setup of the primary delays and the dispatching functionality will be described in the remaining two parts of this section, and as an addition, the non-default settings that have been used (and have influence on the simulation results) are listed in Table C.12 in Appendix C. 3.1.2.1. Primary delays. In the experiments, the following types of primary delays have been used: entry delays, which are added to the scheduled starting time; dwell time extensions, which are added to the minimum dwell time; and running time extensions (which are added in-between stations). Seven delay scenarios, which are based on the same distributions but differ in the average level of running time extensions, have been defined for the experiments. They are listed in Table 4 and have been constructed as follows: •

1

Entry delays are generated for 80% of the trains and they are randomly chosen in the range of 0–6 minutes using the “entry” distribution given in Fig. B.8 (see Appendix B). URL: http://www.rmcon.de/railsys-en/.

202

J. Högdahl, M. Bohlin and O. Fröidh / Transportation Research Part B 126 (2019) 192–212 Table 4 Brief summary of the delay scenarios used in the experiments. The percentage of the average running time extension is computed with respect to the minimum running time for each train class (divided into fast trains, intercity/regional trains, commuter trains, and freight trains).

•

Scenario

Entry delays and dwell time extensions

Average running time extensions (%)

D0 D1 D5 D8 D10 D15 D20

See Appendix B

0 1 5 8 10 15 20

Dwell time extensions are randomly chosen at each stop for all trains using one of the distributions shown in Fig. B.8. Which distribution that is used for a specific train at a certain stop depends on the train’s train type and the station of the stop. A detailed description is provided in Table B.11 (see Appendix B).

The running time extensions have been randomly generated from negative exponential distributions, where one distribution per line section (that connects two stations) and train category have been defined. Each distribution has been generated with the following parameters: • • •

The average running time extension is μ percent of the minimum running time for the specific train category. The proportion of trains that get a running time extension is 100%. The maximum running time extension is 10 minutes (which is approximately in the order of 10–100 times larger than the average running time extension).

The scenarios listed in Table 4 have been defined such that each scenario have a common value of μ (specified in the rightmost column of the table), which determines the average percentage of the minimum running times that are added as running time extensions.

3.1.2.2. Dispatching. The dispatching functionality in Railsys has been disabled. In practice, this implements a no-wait dispatching policy where a movement authority to enter the next block section is granted if it is not currently occupied. If several trains are queuing for a movement authority to enter the same block section, access is granted using a first-in first-out serving policy of the queue.

3.1.3. Optimization model setup The optimization model contains some parameters which must be either chosen by the user or imported. In this section the choice of OD-pair weights is described as well as how the minimum running time and minimum dwell times have been chosen. Also included here is the definition of minimum headways and the chosen value of the Big-M constant, which is used in Eqs. (5) and (15), and , which is used in Eq. (15). In practice, it is desired to have detailed knowledge of the complete OD-matrix. However, this is often considered a commercial trade secret and therefore not available. OD-matrix estimation could be applied but lies outside the scope of this article. For this reason, OD-pair weights have been generated using the following simplified principles: • •

For freight trains we have assumed that all freight is loaded at the first station and is heading to the final station. For passenger trains we have assumed that all passengers board the train at the first station, that 50% of them travel to the end station, and that the remaining passengers disembark the train evenly distributed at the intermediate stations.

The minimum technical running times, as calculated by Railsys, have been used as the minimum running times in the h optimization model (i.e. the parameter di,i if i and i + 1 are events at different stations), and as minimum dwell times (i.e. +1

h the parameter di,i if i and i + 1 are events at the same station) the scheduled dwell times of the initial timetable have +1 been used. According to the recommendations by the Swedish Transport Administration for the 2015 national timetable (Trafikverket, 2014) minimum headways, a and d , shall be as follows: • • •

Gnesta - Hallsberg: 5 minutes Hallsberg - Alingsås: 4 minutes Alingsås - Partille: 5 minutes

Finally, the model parameters and M have been set to = 1 second and M = 105 .

J. Högdahl, M. Bohlin and O. Fröidh / Transportation Research Part B 126 (2019) 192–212

203

3.1.4. Calibration problem setup As mentioned in Section 2.2, it is necessary to calibrate the parameters b (i.e. the train independence threshold) and β (i.e. the knock-on gain) for the predicted average delay G (t | x¯ ) in the problem TTD-fix. Since the calibration objective function, L(·), is stochastic and not explicitly defined (due to it depending on simulation output), a random search for the optimal parameter values for b and β has been carried out (as described in Table 2). To restrict the size of the search space the following assumptions have been made: (1) it is not expected that trains which are separated by a large headway should affect each other; and (2) it is not expected that reduced headway should have larger influence on the expected delay than increased running time supplements. Therefore, the sampling has been limited to the region where the train independence threshold b ∈ [0, 60] minutes and the knock-on gain β ∈ [0, 1.5]. The sampling has been repeated such that 100 samples have been evaluated for each delay valuation coefficient in Table 3, with on average 10% and 20% of the minimum running time as running time extensions (i.e. delay scenarios D10 and D20, which are described in Section 3.1.2.1). 3.2. Results This section presents the results from the computational experiments and is divided into the following parts: Sections 3.2.1 and 3.2.2 is devoted to calibration and evaluation of the calibration; Section 3.2.3 is focused to evaluate optimality of the timetables that this approach yield; and Section 3.2.4 evaluates the operational performance of the optimized timetables. 3.2.1. Parameter calibration The purpose of the calibration is to find values of the train independence threshold b and the knock-on gain β such that the prediction error of the predicted average delay G (t | x¯ ), which is used in the objective function of the timetabling problem TTD-fix, is minimized. As a reminder, the meaning of the parameters b and β is as follows: •

•

The train independence threshold, b, defines the threshold for when a change in the scheduled headway between two consecutive trains affects the average delay of the following train. The knock-on gain, β , defines how much changes in scheduled headway between two trains affect the average delay, compared to changes in the running time supplements.

As was described in Section 3.1.4, the calibration problem has been approximately solved by random sampling of the search space. Fig. 3 shows the outcome of the sampling, where the figures show the result using the delay scenarios with on average 10% and 20% of the minimum running time as running time extensions, respectively. From this figure we see that the calibration objective function tend to have larger values for the scenario with on average 20% running time extensions compared to the scenario with on average only 10% running time extensions. This means that if the primary delays are large, the prediction error of the predicted average delay will be large (in comparison with a scenario with smaller primary delays), which is consistent with expectations. Another observation we make from Fig. 3 is that for the scenario with on average 10% running time extensions, the best parameter values of b and β are located in a region where they affect the predicted average delay little or not at all. However, for the scenario with on average 20% running time extensions, the best parameter values of b and β are sufficiently large, which means that they will have impact on the predicted average delay. This indicates that the magnitude of the primary delays in the scenario with on average 10% running time extensions may be too small to result in headway dependencies between the trains that this model can capture. For this reason, the remaining analysis in this section will focus on the results from the calibration using the scenario with on average 20% running time extensions. Table 5 shows the correlation between the samples in Fig. 3, divided into groups based on which value of the delay valuation coefficient is used when the calibration objective function is evaluated. From this table we see that the correlation

Table 5 The correlation between the samples in Fig. 3, divided into groups based on which value of the delay valuation coefficient is used when the calibration objective function is evaluated. The correlation coefficients in this table are statistically significant ( p < .001) for all pairs of delay valuation coefficients, except for the pair with the delay valuation coefficients α = 0.222 and α = 0.270, for which a p-value of p = .492 were obtained.

0.213 0.222 0.250 0.263 0.270 0.323 0.333

0.213

0.222

0.250

0.263

0.270

0.323

0.333

1.00 0.63 0.64 0.64 0.63 0.66 0.66

1.00 0.61 0.62 0.47 0.68 0.59

1.00 0.59 0.60 0.67 0.68

1.00 0.62 0.60 0.58

1.00 0.63 0.62

1.00 0.61

1.00

204

J. Högdahl, M. Bohlin and O. Fröidh / Transportation Research Part B 126 (2019) 192–212

Delay scenario: D10

1.5

5

4.5

1 4

0.5

3.5

3

0 0

10

20

30

40

50

60

Threshold, b [min]

2.5

Delay scenario: D20

1.5

RMSE [min]

Knock-on gain,

Minimum for each Avg. minimium for all

2

Knock-on gain,

Minimum for each Avg. minimium for all 1.5

1

1

0.5 0.5

0

0

0

10

20

30

40

50

60

Threshold, b [min] Fig. 3. The search space of the calibration problem for the delay scenarios D10 (the upper subplot) and D20 (the lower subplot). Each point shows the value of the loss function of the calibration problem. The best parameter value for each delay valuation coefficient is marked with a black diamond, and the average of those points is marked with the red star. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Table 6 The best values of b and β for each delay valuation coefficient. Delay valuation coefficient

b (min)

β

0.213 0.222 0.250 0.263 0.270 0.323 0.333 Average:

34.8 28.2 47.0 42.2 25.1 27.3 21.8 32.4

0.19 0.16 0.38 0.23 0.11 0.10 0.38 0.22

coefficients are in the range of 0.47 − 0.68, which means that there is a moderate positive correlation between the groups. The correlation coefficients are statistically significant ( p < .001) except between the delay valuation coefficients α = 0.222 and α = 0.270, where the correlation coefficient is 0.47 and the p-value is only 0.492. Nonetheless, this indicates that if a specific value of b and β is good for one choice of the delay valuation coefficient, then it is likely that the same b and β is also a good choice for any other value of the delay valuation coefficient. This means that the same values of b and β may be used to study the effects on different values of the delay time factor (i.e. the value of decreasing the average delay time compared with decreasing the scheduled in-vehicle travel time). Table 6 summarizes the best values of b and β for each delay valuation coefficient, together with their respective average best values. Since the samples are noisy, the best values of b and β for each delay valuation coefficient will also be noisy. Therefore, we will use the average best values of b and β in the remainder of this article (unless something else is stated). This gives b = 32.4 minutes and β = 0.22, which is interpreted as follows. The value of the train independence threshold b implies that when the headway between two consecutive trains is less than 32.4 minutes, then a decreased headway to the train ahead will lead to the average delay for the following train to

J. Högdahl, M. Bohlin and O. Fröidh / Transportation Research Part B 126 (2019) 192–212

205

Average error [%]

40

30

20

10

0 0

2

4

6

8

10

12

14

16

18

Average running time extension compared with the minimum running times [%] Fig. 4. The relation between delay scenario used in the simulation and the average percentage error of the predicted average delay and its alternative model. Each mark show the average percentage error for the optimal timetables obtained by solving TTD-fix using the delay valuation coefficients α listed in Table 3. The lines show the regression line of each data set.

increase (and vice versa). This, rather large, value of the train independence threshold is likely a consequence of that large primary delays have been used in the simulation (on average the running time extension for each train has been about 20% of its minimal running time). Since large delays might cause delay propagation even between pair of trains that are separated with a large headway, it appears reasonable that the average delay of a following train might be affected even for small changes between trains that in the initial timetable was separated with a large headway. The effect of a decreased headway to the train ahead, with respect to the average delay for the following train, is determined by the knock-on gain β . The value β = 0.22 means, approximately, that for each 4-minute decrease of the headway (e.g. by adding a running time supplement to the train ahead), 1 additional minute of running time supplement for the following train is required in order to not increase its average delay. 3.2.2. Prediction error of the predicted average delay In this section the precision (here evaluated as the percentage error) of the predicted average delay is evaluated and compared with a related (and simplified) model to calculate the predicted average delay, in which the contribution of the headway between trains is excluded, which is calculated as follows:

Galt (t

| x¯ ) =

 

h∈T j∈Eah ω∈h j

   λhj phj,ω max 0, t¯hj + ω − t hj ,

(21)

where Galt (t | x¯ ) will be referred to as the predicted average delay of the alternative model. Fig. 4 show the relation between the delay scenario that has been used in the simulation and the average percentage error for the optimal timetables obtained by solving TTD-fix using the delay valuation coefficients α listed in Table 3. The blue and red marks show the percentage error for the predicted average delay and its alternative model, respectively, and the lines show the corresponding regression lines. An interpretation of the differences in the inclination for the regression lines in Fig. 4 is that the predicted average delay is roughly 50% less sensitive to variations in the magnitude of primary delays. This is due to the higher precision that the predicted average delay achieves for the delay scenario with on average 20% running time extensions, which is a consequence of the additional information in the model obtained by including the relation between headway and delay time. The increased robustness against variations in the magnitude of primary delays also means that this model seems to be less sensitive to imperfections in the modeling of delay distributions. 3.2.3. Evaluation of the combined simulation-optimization approach By solving the timetable problem (TTD-fix) for a certain choice of parameters b and β we obtain, for the resulting timetable, a prediction of the minimal weighted sum of scheduled travel time and average delay. However a different choice of parameters b and β , will yield a different solution with a different predicted optimal value. To evaluate if the timetable generated by solving the optimization problem TTD-fix with the parameters b and β that was obtained by the calibration actually minimizes the weighted sum of scheduled travel time and observed average delay (i.e. instead of the predicted average delay) a random search according to Table 7 has been carried out, where the values of b and β has been sampled in the region b ∈ [0, 1.5] and β ∈ [0, 60] min (i.e. the same region that was used during the

206

J. Högdahl, M. Bohlin and O. Fröidh / Transportation Research Part B 126 (2019) 192–212

100

min sequence max sequence random sequence b = 28.2, =0.16

TPI

95

90

85 0

50

100

150

200

250

300

350

400

Iteration Fig. 5. The results from the random search described in Table 7. TPI denote the timetable performance index, here computed as the weighted sum of scheduled travel time and observed average delay. As a comparison, the TPI for the timetable generated with the parameter values found by the calibration is included (the dashed purple line). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Table 7 A random search method to minimize the scheduled travel time and observed average delay. Algorithm 2: A random search method to minimize the weighted sum of scheduled travel time and observed average delay. Step 0: Step 1: Step 2: Step 3:

Step 4:

Select the initial timetable t¯, the sample space , the number of iterations N, the delay valuation coefficient α , and set the minimal objective function value Lmin = ∞. Simulate the initial timetable t¯ to generate observed arrival and departure times of its operation. Randomly choose b and β from , generate a new current timetable, t, by solving TTD-fix with parameters α , b, and β , and simulate its operation. Compute the weighted sum of scheduled travel time and observed average delay L(t ) = F (t ) + cD (t ), such that F(t) and D(t) is the scheduled travel time and the observed average delay, respectively, of the current timetable, and c is the delay time factor. If L(t ) < Lmin , then update Lmin := L. Terminate if the procedure has been repeated N times, else return to Step 2.

calibration). The result of the random search is a sequence of timetables which all minimize the optimization problem TTDfix for a random value of the parameters b and β . Each timetable is evaluated by computing the timetable performance index (TPI), i.e. the weighted sum of scheduled travel time and observed average delay, which is a measure to evaluate timetable quality proposed by Warg (2016). Fig. 5 shows the result of the random search described in Table 7, where the random sequence of the TPI is given by the thin yellow line, the maximum sequence and the minimum sequence of the TPI is given by the tick red and blue lines, respectively, and the TPI for the timetable generated by solving the optimization problem TTD-fix with the parameter values found by the calibration is given by the dashed purple line. From this figure it is concluded that, for the sampled values of b and β (N = 448), the proposed method generates a near optimal timetable, where the gap is 1.2% from the timetable with minimal TPI (which, in this case, was obtained for b=15.6 minutes and β = 0.27). The main reason for this gap, which is expected, is that the calibration is based on minimizing a different objective, which is: the root mean square error of the prediction error of the predicted average delay for each event. As the error for each event is squared, the calibration objective might favor solutions with moderate errors in all events over a solution with many small and some large errors, which then is likely to yield different optimal parameter values compared with the approach described in Table 7. 3.2.4. Operational performance of the optimized timetables The choice of delay valuation coefficients (see Table 3) corresponds to valuating decreased delay time between 2.0 and 3.7 times higher than reducing the scheduled in-vehicle travel time. It is therefore expected that the scheduled travel times of the optimized timetables should be larger compared with the initial timetable, that the observed average delay should decrease, and that the punctuality should improve. The performance of an optimized timetable compared with the initial timetable can be evaluated by utilizing the timetable performance index (TPI), as proposed by Warg (2016) as a measure for timetable quality evaluation. Let TPIx denote the TPI of timetable x, and calculate it as the weighted sum of scheduled travel time and observed average delay, i.e:

T P Ix = T Tx + cDx ,

(22)

where TTx is the scheduled travel time of timetable x, c is the delay time factor (i.e. the value of decreasing the average delay time compared with decreasing the scheduled in-vehicle travel time), and Dx is the observed average delay. The performance difference (in percent) for an optimized timetable o, compared with the initial timetable i is then

J. Högdahl, M. Bohlin and O. Fröidh / Transportation Research Part B 126 (2019) 192–212

8

Delay cost factor 3.7 3.5 3.0 2.8 2.7 2.5 2.0

6

TPI [%]

207

4 2 0 D0

D1

D5

D8

D10

D15

D20

Delay scenario Fig. 6. Difference (percentage) in TPI between the optimized timetables and the initial timetable for the set of delay scenarios and different delay time factors (i.e. the value of decreasing the average delay time compared with decreasing the scheduled in-vehicle travel time).

Punctuality [%]

100 95 90 85 80 75 70 D0

D1

D5

D8

D10

D15

D20

Delay scenario Average for the optimized timetables

Initial timetable

Fig. 7. End station punctuality for the different delay scenarios. The bars for the optimized timetables show average punctuality with respect to the different values of the delay valuation coefficient in Table 3.

T PIo =

100 (T PIi − T PIo ). T P Ii

(23)

Fig. 6 shows the performance difference of the optimized timetables compared with the initial timetable for each delay scenario (D0-D20, as defined in Table 4). Each color represents a single delay time factor (between 2.0 and 3.7), and each bar shows the performance difference obtained by optimizing the initial timetable with the corresponding delay valuation coefficient α and a specific delay scenario. This figure shows that, on average for each delay scenario, the TPI is reduced by between 0.04% and 4.64%. This means that the improvement in terms of reduction of delays is larger than the simultaneous deterioration of the scheduled travel times. We also observe that the TPI reduction is proportional to both the cost of delay and the magnitude of the primary delays, which is logical as the potential delay reduction is proportional to the total average delay and the benefit of reducing delays is proportional to the value of the delay time factor. The concluding interpretation is that this approach yields socio-economic improvements of the timetable in terms of the balance between scheduled travel time and average delay. In particular, the proposed approach is suitable in cases where the value of the delay time factor is high, for instance in Sweden (delay time factor 3.5) or New South Wales in Australia (delay time factor 3.7), in combination with running time extensions of either 15% or 20% of the minimum running times. In this case, improvements between 3.25% and 7.43% are observed. Fig. 7 shows end station punctuality for each delay scenario (D0-D20, see Table 4) for the optimized timetables and for the initial timetable. This figure shows that punctuality of the optimized timetables are slightly higher for scenario D1-D10, and substantially improved for scenario D15 and D20. This is consistent with expectations as it should be hard to further improve the punctuality when it is very high (which is the case for the delay scenarios with up to on average 10% running time extensions with respect to the minimum running times). Another observation we make is that when the punctuality is moderate or poor, as in the scenarios D15 and D20, there is a large improvement (from 90.3% and 75.4% to 97.1% and 93.4%, respectively). Put together, this indicates that the potential improvement of this approach is inverse proportional to the punctuality of the initial timetable, meaning that the benefit of this approach is larger for an initial timetable with poor punctuality than for an initial timetable with good punctuality. This means that poor initial punctuality might be a good indicator for when this approach is suitable to apply in railway timetable planning.

208

J. Högdahl, M. Bohlin and O. Fröidh / Transportation Research Part B 126 (2019) 192–212

To relate the results from the analysis of punctuality, it may be mentioned that the average punctuality of all passenger trains in Sweden is about 90%, and has been so since 2013 (Transportanalys, 2018). However, for long distance passengers trains and freight trains the punctuality is worse and for that reason the scenarios that appears to be closest to real world is scenario D15 and D20. As the socio-economic improvement seems to be proportional to the delay time factor and the magnitude of the primary delays, and since the initial punctuality is inversely proportional to the magnitude of primary delays, this approach is best suited for railway lines with poor punctuality in countries which use a high delay time factor. 4. Conclusions and future work The main contribution of this article is a model that is able to use observed times of events, such as arrival times obtained either from simulation or real operations, to predict the average delay as the timetable is optimized. It differs from existing simulation-based optimization methods for railway timetabling since the effects of secondary delays have been included explicitly. The formulation of the problem in this article also differs from the usual train timetabling problem as the objective here is to minimize the weighted sum of scheduled travel time and predicted average delay, to improve the socio-economic benefit of a given timetable. To demonstrate the approach and to validate the model a thorough computational study, on a large part of the Swedish Western Main Line, has been carried out, and based on these results we draw the following conclusions: •

•

•

•

•

By including the effects of secondary delays in the predicted average delay G (t | x¯ ), the percentage error as well as the sensitivity to modeling imperfections in the primary delays may be reduced by between 0 and 50%. For the set of timetables that are optimal solutions to the optimization problem TTD-fix, the approach in this paper yield a timetable that seems to be near optimal with respect to the weighted sum of scheduled travel time and observed average delay. This approach yields socio-economic benefits and increased punctuality compared to the initial timetable. On average we observed socio-economic improvements in the range of 0–5%, depending on the magnitude of the primary delays, and increased punctuality in the range of 0–25%, depending on the initial punctuality. For the delay scenarios where the punctuality of the initial timetable is comparable to the real-world punctuality we observed socio-economic improvement on average in the range of 2–5% and improved punctuality on average in the range of 8–25%. The most substantial improvements by this approach are obtained on railway lines with poor punctuality in countries that use a high delay time factor in socio-economic assessments of the transport system.

The impact of parameters and assumptions in the simulation model has not been studied in this article. However, the findings of this article are based on comparisons between similar alternatives, and even though some of the obtained results might differ we anticipate that the overall conclusions in this paper will still hold. Some issues that will be addressed in future work include studying the effects of different dispatching strategies and including train connections in the timetable (e.g. for passenger exchanges or circulation of rolling stock). Further extension of the simulation methodology could include a systematic approach to non-independent disturbances, such as: modeling malfunction of trains and infrastructure; using randomized train sets with different weights and lengths (i.e. randomize the trains’ characteristics); and allowing freight train departures before schedule (common praxis in Sweden). The following topics could be addressed in future research in order to increase the generality and applicability of the approach: extend the optimization model so that it can be applied for simultaneous optimization in both directions, so that it can be used on sections with single as well as multiple tracks, and so that it can take into account traffic from branch lines; increase the degree of freedom in the optimization model (in this article the first departure time of each train, as well as the train order, is predetermined); and include constraints from the operators’ perspective, such as circulation times, etc. As discussed in Section 2.1.2 it is also of interest to investigate alternative formulations of the predicted average delay model. Acknowledgments This research was funded by Swedish Transport Administration (Trafikverket) grant TRV 2016/5090. The authors would like to express their gratitude to Magnus Wahlborg (Trafikverket) for providing the necessary Railsys infrastructure models and timetable data. This article is an extension of Högdahl et al. (2017). Appendix A. The complete linearized optimization model This section presents the complete linearized optimization model used in the experiments. For a description of each set, variable and parameter in the model, see Tables A.8, A.9, and A.10, respectively. When (h, j) is used in this section, this is shorthand for the j:th event of train h.

J. Högdahl, M. Bohlin and O. Fröidh / Transportation Research Part B 126 (2019) 192–212

Table A.8 List of all sets in the model. Set

Description

T S As Ds Ea Eah ODh

The The The The The The The The The

hj Pjh

set set set set set set set set set

of of of of of of of of of

all trains. all stations. all arrival events at station s ∈ S. all departure events from station s ∈ S. all arrival events. all arrival events for train h all trips for train h. all observed delays of event (h, j). the event that directly precedes event (h, j).

Table A.9 List of all variables in the model. Variable

Type

Description

t hj yhj,ω h,k j,l

Decision variable State variable



State variable

zhj

State variable

The scheduled time of event (h, j). The outcome delay time of event (h, j) that has probability phj,ω . If t hj − tlk < b, then this variable measures the change in headway between those events caused by rescheduling event (k, j). Binary variable that is 1 if and only if the headway to the train ahead is less than threshold value b.

Table A.10 List of all parameters in the model. Parameter

Description

x¯h,k j,l

A binary parameter that is 1 if and only if t hj < tlk in the initial timetable. The delay valuation coefficient. The weight (e.g. number of passengers) of the trip with train h that starts with event (h, i) and ends with event (h, j).

α λhi, j λhj

phj,ω

A weighting parameter for the arrival event j for train h (e.g the number of disembarking passengers). The probability that event (h, j) was delayed by ω.

dhj, j+1

The minimum time separation between consecutive events of the same train.

as as

The minimum headway between two arrival events at station s. The minimum headway between two departure events from station s.

b

The train independence threshold.

β

The knock-on gain.

M t¯1h

A parameter for simulating the strict less than relation in Eq. (12). A Big-M constant. The departure time from the first station for train h in the initial timetable.

209

210

J. Högdahl, M. Bohlin and O. Fröidh / Transportation Research Part B 126 (2019) 192–212

min f (t ) = α





h∈T (i, j )∈ODh

+ (1 − α )

  λhi, j t hj − tih

  

h∈T j∈Eah ω∈h j

λhj phj,ω yhj,ω

s.t. t hj+1 − t hj ≥ dhj, j+1 , t hj

− tlk ≥

∀(h, j ), (k, l ) ∈ As,

t hj − tlk ≥ ds − Mx¯h,k , j,l

(k, l ) = (h, j ), ∀s ∈ S, ∀(h, j ), (k, l ) ∈ Ds,

tlk − t hj ≥ ds − Mx¯k,h , l, j

(k, l ) = (h, j ), ∀s ∈ S, ∀(h, j ), (k, l ) ∈ Ds,

− tlk

− b≥

−Mzhj ,





(k, l ) = (h, j ), ∀s ∈ S, (k, l ) ∈ Pjh , ∀(h, j ) ∈ Ea,

+ (b − )≥ −M 1 − zhj , (k, l ) ∈ Pjh ,



h,k ≥ j,l

−



 h

−Mzhj ,

h,k ≥ j,l



−Mzhj ,

 k

M 1 − z j − tlk − t¯l ≥ −

 k

tlk − t¯l



∀h ∈ T ,

(k, l ) = (h, j ), ∀s ∈ S, ∀(h, j ), (k, l ) ∈ As,

+ tlk



j = 1, . . . , Nh − 1,

Mx¯h,k , j,l

tlk − t hj ≥ as − Mx¯k,h , l, j

t hj −t hj

 − a s



 h

h,k , j,l

+ M 1 − z j ≥ h,k , j,l



− t¯hj + ω + t hj β h,k ≥ −yhj,ω , j,l t1h = t¯1h , zhj ∈ {0, 1} yhj,ω ≥ 0

∀(h, j ) ∈ Ea,

(k, l ) ∈

Pjh ,

∀(h, j ) ∈ Ea,

(k, l ) ∈

Pjh ,

∀(h, j ) ∈ Ea,

(k, l ) ∈ Pjh , ∀(h, j ) ∈ Ea, (k, l ) ∈ Pjh , ∀(h, j ) ∈ Ea,

∀ω ∈ hj , (k, l ) ∈ Pjh, ∀(h, j ) ∈ Ea , ∀h ∈ T , ∀(h, j ) ∈ Ea , ∀ω ∈ hj , ∀(h, j ) ∈ Ea .

Appendix B. Entry delays and dwell time extensions Entry delays and dwell time extensions have been applied using the distributions given in Fig. B.8. The distribution entitled “entry” is an entry delay distribution and the rest are distributions for the dwell time extensions. Entry delays have been applied at the first station for 80% of the trains and dwell time extensions have been applied for each train at each stop according to Table B.11.

J. Högdahl, M. Bohlin and O. Fröidh / Transportation Research Part B 126 (2019) 192–212

small

small commuter

0.1

prob. density

prob. density

0.04

empirical pdf avg. delay

0.02

empirical pdf avg. delay

0.05

0

0 0

100

200

300

400

0

20

40

delay [s]

80

100

large commuter

0.04

prob. density

prob. density

60

delay [s]

medium

0.02

empirical pdf avg. delay

0.01

empirical pdf avg. delay

0.02

0

0 0

100

200

300

400

0

20

40

delay [s]

60

80

100

delay [s]

large

10-3

6

prob. density

0.02

prob. density

211

empirical pdf avg. delay

0.01

0

entry empirical pdf avg. delay

4 2 0

0

100

200

300

400

0

100

delay [s]

200

300

delay [s]

Fig. B.8. The probability density functions of the dwell time extensions that have been used in the simulations. The red dotted line shows the mean value. The distributions entitled “small commuter” and “large commuter” are only used for commuter trains. The distribution “entry” is only used for entry delays and is applied for 80% of the trains. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Table B.11 Dwell time distributions used for each station and train type. Abbreviations are as follows: S = small, M = medium, L = large, SC = small commuter, LC = large commuter. Station

Freight

Commuter

Regional

Fast

Alingsås Aspedalen Falköping Floby Floda Flen Gårdsjö Gnesta Hallsberg pbg Herrljunga Jonsered V Katrineholm Lerum Laxå Laxå Östra Skövde Stenkullen Stenstorp Töreboda Västra bodarna Vårgårda Vingåker Väring

L

LC SC

M

M

M M S S S S M

L

M

M

M

L

L L S L

SC L SC L L

S M LC S M SC M S

L

L

212

J. Högdahl, M. Bohlin and O. Fröidh / Transportation Research Part B 126 (2019) 192–212

Appendix C. Railsys settings Table C.12 lists some important settings that have been made which have influence on the simulation result. Table C.12 Settings that have been used in the simulation. Setting

Value

Comment

Portion of allowances that can be used to reduce lateness Vmax from lateness [s] Vnormal below lateness [s] Routing methods

100% 0 sec 0 sec Disabled

Trains are able to run at maximum speed. Maximum speed is always used when trains are late. Maximum speed is maintained until the lateness is completely removed. All routing methods have been disabled.

References Andersson, E.V., 2014. Assessment of robustness in railway traffic timetables Lic. thesis. Linköping University, Sweden. Andersson, E.V., Peterson, A., Törnquist Krasemann, J., 2013. Quantifying railway timetable robustness in critical points. J. Rail Transp. Plann. Manage. 3 (3), 95–110. Bates, J., Polak, J., Jones, P., Cook, A., 2001. The valuation of reliability for personal travel. Transp. Res. Part E, 37 (2–3), 191–229. Bendfeldt, J., Mohr, U., Muller, L., 20 0 0. Railsys, a system to plan future railway needs. In: Computers in Railways VII. WIT Press, pp. 249–255. Burggraeve, S., Vansteenwegen, P., 2017. Optimization of supplements and buffer times in passenger robust timetabling. J. Rail Transp. Plann. Manage., 7 (3), 171–186. Cacchiani, V., Caprara, A., Fischetti, M., 2012. A lagrangian heuristic for robustness, with an application to train timetabling. Transp. Sci., 46 (1), 124–133. Cacchiani, V., Toth, P., 2012. Nominal and robust train timetabling problems. Eur. J. Oper. Res., 219 (3), 727–737. Caimi, G., Kroon, L., Liebchen, C., 2017. Models for railway timetable optimization: applicability and applications in practice. J. Rail Transp. Plann. Manage., 6 (4), 285–312. Carrion, C., Levinson, D., 2012. Value of travel time reliability: a review of current evidence. Transp. Res. part A, 46 (4), 720–741. Chang, J.S., 2010. Assessing travel time reliability in transport appraisal. J. Transp. Geo., 18 (3), 419–425. Chow, A.H., Pavlides, A., 2018. Cost functions and multi-objective timetabling of mixed train services. Transp. Res. Part A, 113, 335–356. Concas, S., Kolpakov, A., 2009. Synthesis of Research on Value of Time and Value of Reliability. Technical Report. Florida Department of Transportation, USA. Dewilde, T., 2014. Improving the Robustness of a Railway System in Large and Complex Station Areas. KU Leuven, Belgium. Fischetti, M., Monaci, M., 2009. Light robustness. In: Robust and Online Large-Scale Optimization. Lecture Notes in Computer Science. Springer, pp. 61–84. Fischetti, M., Salvagnin, D., Zanette, A., 2009. Fast approaches to improve the robustness of a railway timetable. Transp. Sci., 43 (3), 321–335. Harrod, S.S., 2012. A tutorial on fundamental model structures for railway timetable optimization. Surveys Oper. Res. Manage. Sci., 17 (2), 85–96. Högdahl, J., Bohlin, M., Fröidh, O., 2017. Combining optimization and simulation to improve railway timetable robustness. In: Paper presented at 7th International Conference on Railway Operations Modelling and Analysis - RailLille2017, Lille. ´ P., Kecman, P., Bojovic, ´ N., Mandic, ´ D., 2017. Optimal allocation of buffer times to increase train schedule robustness. Eur. J. Oper. Res., 256 (1), Jovanovic, 44–54. De Jong, G., Tseng, Y., Kouwenhoven, M., Verhoef, E., Bates, J., 2007. The Value of Travel Time and Travel Time Reliability. Technical Report. The Netherlands Ministry of Transport, Public Works and Water Management. Khoshniyat, F., Peterson, A., 2017. Improving train service reliability by applying an effective timetable robustness strategy. J. Intell. Transp. Syst., 21 (6), 525–543. Kroon, L., Maróti, G., Helmrich, M.R., Vromans, M., Dekker, R., 2008. Stochastic improvement of cyclic railway timetables. Transpo. Res. Part B, 42 (6), 553–570. Lee, Y., Lu, L.-S., Wu, M.-L., Lin, D.-Y., 2017. Balance of efficiency and robustness in passenger railway timetables. Transp. Res. Part B, 97, 142–156. Parbo, J., Nielsen, O.A., Prato, C.G., 2016. Passenger perspectives in railway timetabling: a literature review. Transp. Rev., 36 (4), 500–526. Schlechte, T., Borndörfer, R., 2010. Balancing efficiency and robustness–a bi-criteria optimization approach to railway track allocation. In: Multiple Criteria Decision Making for Sustainable Energy and Transportation Systems. Springer, pp. 105–116. Schöbel, A., Kratz, A., 2009. A bicriteria approach for robust timetabling. In: Robust and Online Large-Scale Optimization. Lecture Notes in Computer Science. Springer, pp. 61–84. Sels, P., Dewilde, T., Cattrysse, D., Vansteenwegen, P., 2016. Reducing the passenger travel time in practice by the automated construction of a robust railway timetable. Transp. Res. Part B, 84, 124–156. Solinen, E., Nicholson, G., Peterson, A., 2017. A microscopic evaluation of railway timetable robustness and critical points. J. Rail Transp. Plan. Manage., 7 (4), 207–223. Spall, J.C., 2005. Introduction to Stochastic Search and Optimization: Estimation, Simulation, and Control, 65. John Wiley & Sons. Trafikverket, 2014. Riktlinjer tåthet mellan tåg: Planeringsförutsättningar: Tågplan T15 [Guidelines for train density: planning prerequisites: Train plan T15]. Transportanalys, 2018. Punktlighet på järnväg 2017, Statistik 2018:7 [Train performance 2017, Statistics 2018:7]. Tseng, Y.-Y., 2008. Valuation of travel time reliability in passenger transport. Vrije University, Amsterdam, The Netherlands. Warg, J., 2016. Timetable evaluation with focus on quality for travellers Lic. thesis. KTH Royal Institute of Technology, Stockholm, Sweden. Warg, J., Bohlin, M., 2016. The use of railway simulation as an input to economic assessment of timetables. J. Rail Transp. Plann. Manage, 6 (3), 255–270.