A delay propagation algorithm for large-scale railway traffic networks

Transportation Research Part C 18 (2010) 269–287

Contents lists available at ScienceDirect

Transportation Research Part C journal homepage: www.elsevier.com/locate/trc

A delay propagation algorithm for large-scale railway traffic networks Rob M.P. Goverde * Department of Transport and Planning, Delft University of Technology, P.O. Box 5048, 2600 GA Delft, The Netherlands

a r t i c l e

i n f o

Article history: Received 26 September 2007 Received in revised form 7 January 2010 Accepted 21 January 2010

Keywords: Railway timetable Railway traffic Train delay Max-plus algebra Timed event graph

a b s t r a c t In scheduled railway traffic networks a single delayed train may cause a domino effect of secondary delays over the entire network, which is a main concern to planners and dispatchers. This paper presents a model and an algorithm to compute the propagation of initial delays over a periodic railway timetable. The railway system is modelled as a linear system in max-plus algebra including zero-order dynamics corresponding to delay propagation within a timetable period. A timed event graph representation is exploited in an effective graph algorithm that computes the propagation of train delays using a bucket implementation to store the propagated delays. The behaviour of the delay propagation and the convergence of the algorithm is analysed depending on timetable properties such as realisability and stability. Different types of delays and delay behaviour are discussed, including primary and secondary delays, structural delays, periodic delay regimes, and delay explosion. A decomposition method based on linearity is introduced to deal with structural and initial delays separately. The algorithm can be applied to large-scale scheduled railway traffic networks in real-time applications such as interactive timetable stability analysis and decision support systems to assist train dispatchers. Ó 2010 Elsevier Ltd. All rights reserved.

1. Introduction Train punctuality is a main performance indicator for railways. European railways are typically operated according to a master timetable. This timetable represents a conflict-free coordination of train paths and includes slack time to manage train delays. Moreover, buffer times between train paths prevent that slightly delayed trains immediately hinder following trains. Nevertheless, depending on the distribution of running time supplements and buffer times, a single delayed train may cause a domino effect of secondary delays over the entire network due to train connections and route conflicts. This delay propagation behaviour can be studied by a mathematical model that computes for a given set of initial delays all arrival and departure delays of any train at any station, and thus also the settling time after which all delays have been absorbed by available recovery time. A delay propagation model can be used by timetable planners to study the robustness of timetable designs to delays, and by dispatchers to recognise the impact of delays and to evaluate dispatching proposals. A scheduled railway system can be modelled efficiently as a discrete-event dynamic system using max-plus algebra (Braker, 1993; Heidergott and De Vries, 2001; Goverde, 2007). In the case of periodic timetables the dynamics of train sequences and synchronisations are simple recursive linear equations in max-plus algebra, which can be viewed as the linear state-space representation of a timed event graph (Baccelli et al., 1992). The deterministic max-plus model assumes that each event waits until all required resources are available and maintains the fixed order of all events according to the given timetable. Furthermore, late trains use available recovery times and propagation of delays is reduced by available buffer times. The max-plus recursions can be used to effectively compute the propagation of train delays over the railway network

* Tel.: +31 15 2783178; fax: +31 15 2783179. E-mail address: [email protected] 0968-090X/$ - see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.trc.2010.01.002

270

R.M.P. Goverde / Transportation Research Part C 18 (2010) 269–287

in both time and space, which is the topic of this paper. The bucket-based graph algorithm presented in this paper has been implemented in the tool PETER (Performance Evaluation of Timed Events in Railways) (Goverde and Odijk, 2002; Goverde, 2005), which has been developed as a decision support to railway planners for evaluating and improving timetable stability. The algorithm can also be applied in a decision support system to assist train dispatchers in evaluating various dispatching measures to manage train delays (Goverde, 1998; Heidergott and De Vries, 2001; De Schutter et al., 2002). Goverde (2007) gives a general overview of the max-plus stability theory of railway timetables using a polynomial matrix representation to deal with zero-order terms explicitly, and includes a straightforward recursive procedure to compute delay propagation over the network. The zero-order terms correspond to constraints between events falling in the same basic timetable period and are necessary to compute the delay propagation from and over all events in the system, rather than in a reduced system with aggregated processes and removed internal events. The present paper covers an in-depth discussion of the delay propagation behaviour of max-plus linear systems, again including zero-order terms explicitly. The model used in this paper also deals with arrival, departure, and passing events explicitly rather than departures only. The main new contribution is an efficient delay propagation algorithm with a performance proportional to the number of delayed events. The underlying theory generalises results from Braker (1993) to higher-order systems. Also a marked precedence graph is introduced which generalises the precedence graph associated to a max-plus matrix of a pure first-order max-plus system to higher-order systems (including zero-order terms), as a static counterpart of timed event graphs. The interplay between timetable properties and delay propagation behaviour is investigated thoroughly and used to prove convergence properties of the delay propagation algorithm. Structural delays due to infeasible schedules are covered in depth and a decomposition method based on system linearity is proposed to effectively study the impact of different delay scenarios and separate structural delays from delay propagation of initial delays. Furthermore, the paper discusses efficient data structures for large-scale sparse max-plus linear systems, the event list for selecting the next event in the algorithm, and delay lists. The paper is structured as follows. Section 2 explains the modelling of scheduled railway systems in max-plus algebra and timed event graphs. The delay propagation model and underlying theory of the algorithm are considered in Section 3. The bucket-based delay propagation algorithm is presented in Section 4 and its convergence with respect to timetable properties is analysed in Section 5, along with a discussion of timetable and model properties on delay propagation behaviour. Section 6 ends with conclusions and open extensions. 2. The scheduled railway traffic model The propagation of train delays over the railway network depends on timetable and infrastructure constraints, and in particular on the interconnection structure and the distribution of time supplements and buffer times over the network. The variables of interest are the train event times at stations, which are connected by activities or processes (train runs, stops, transfers, etc.). This paper assumes that trains are operated according to a periodic railway timetable as is common in European railways. The length of a timetable period, the cycle time T, is a common multiple of the line cycle times which allows the daily timetable to be a repetition of a basic timetable period. Often T ¼ 60 min according to a basic hour timetable. A periodic railway timetable gives the scheduled arrival and departure times within a basic period of length T for each train line at all served stations. In addition, also the scheduled passing times at stations where a train does not stop are given, as well as the scheduled passing times at e.g., junctions of merging railway lines, movable bridges, and any other ‘timetable point’ in the network where a minimum headway time has to be satisfied corresponding to a safe separation distance on conflicting train routes. The timetable periodicity can be exploited in a discrete-event dynamic system by modelling event times as a function of the period. An event i corresponds to a triple i ¼ ðEi ; Li ; Si Þ, where Ei is the event type (arrival, departure, passage), Li is the associated train line, and Si is a location (station, junction, signal). The event time xi ðkÞ denotes the kth occurrence time of a periodic event i. The counter k represents the period (hour) in which the event is scheduled to occur. If the initial timetable period is ½0; TÞ then the kth period is ½k  T; ðk þ 1Þ  TÞ. Note that k  T ¼ kT is k times T, and ½0; TÞ is the semi-open interval 0 6 x < T. Hence, time T lies in the next period ½T; 2  TÞ. If all event times and process times are integers (and so is T) then the initial period is ½0; T  1, the next is ½T; 2T  1, etc. 2.1. Timetable constraints A timetable gives constraints to the train event times. In particular, a passenger train is usually not allowed to depart early from a station. Let di ðkÞ denote the scheduled event time of event i in period k. If event i may not occur before its scheduled departure time then the event times must satisfy the timetable constraint

xi ðkÞ P di ðkÞ:

ð1Þ

If an early event time is not prohibited, such as early arrival and passing times at stations, then a redundant timetable constraint can be defined by setting di ðkÞ equal to the previous scheduled (departure) event time of line Li plus the subsequent trip time to event i. In principle, we could also just set di ðkÞ ¼ 1 when early events are allowed but to compute a possible delay we still require a finite reference event time di ðkÞ. Hence, without loss of generality we assume throughout that each event i has finite scheduled event times di ðkÞ for which (1) must be satisfied. In conjunction with the deterministic

R.M.P. Goverde / Transportation Research Part C 18 (2010) 269–287

271

precedence constraints defined below this implies that punctual trains use their scheduled running times while delayed trains speed up and exploit running time supplements to recover from delays (assuming that the scheduled running times have a positive running time supplement). In a periodic timetable with cycle time T, the scheduled event time of event i in successive periods k 2 N is given by 0 di ðkÞ ¼ di ð0Þ þ k  T, where di ð0Þ ¼ di 2 ½0; TÞ is an initial scheduled event time. So for example, if T ¼ 60 min and di ð0Þ ¼ 11 min, then di ð1Þ ¼ 71 min; di ð2Þ ¼ 131 min; di ð3Þ ¼ 191 min, etc. 2.2. Precedence constraints The events are connected to each other by activities or processes. An event can only occur if each preceding process has been completed. For instance, a train departure can occur only if the physical train has arrived and the dwell process is completed, possible feeder trains have arrived and transferring passengers have boarded the train, and possible conflicting routes of preceding trains have been released. So before an event i is ready to occur it must satisfy a number of precedence constraints. In general, these precedence constraints take the form

xi ðkÞ P aij þ xj ðk  lij Þ;

ð2Þ

for all predecessor events j of i. Here, aij 2 RP0 is the minimum process time from event j to i, and lij 2 N0 :¼ N [ f0g is the period shift denoting whether both events i and its predecessor j are scheduled in the same period ðlij ¼ 0Þ, or whether the preceding event j is scheduled in a previous period ðlij P 1Þ. Thus, if lij ¼ 0 then the scheduled process from j to i falls entirely within the same period, while for lij ¼ 1 the process crosses a period boundary, and in general a process from j to i crosses lij periods. 0 0 The period shift lij 2 N0 from an event j to i depends on the scheduled event times dj and di , and the scheduled process 0 0 time a0ij P 0. Note that the scheduled process time may be larger than the cycle time T, while dj ; di 2 ½0; TÞ. For a given timetable the period shifts can be determined as (Goverde, 2005, 2007) 0

lij ¼

0

a0ij þ dj  di T

ð3Þ

: 0

0

Note that this period shift formula generates integers lij 2 N0 , since dj þ a0ij ¼ di mod T. The scheduled process time a0ij is the sum of the minimum process time aij and a time supplement or buffer time sij . The latter can be computed as 0 0 sij ¼ di  dj  aij mod T provided that 0 6 sij < T. A supplement or buffer time can however also be negative or larger than T, see Section 3.2. For instance, when an arriving train at the end station has a minimum turn-around time that exceeds the scheduled departure time for the return trip by 1 min. Then a planner can either choose a generous layover buffer time of 59 min or spare rolling stock by allowing 1 min departure delay. If on the other hand, the layover process is ready just 1 min before the scheduled departure time for the return trip, then a planner may decide that 1 min layover buffer time is too small for a stable train circulation and calculates 61 min buffer time, which requires an additional train to take over the return trip. Hence, deriving the correct scheduled process times takes some care. In particular for rolling stock connections where rolling stock is (cyclically) exchanged between different train lines. The precedence constraint (2) can be used for all line schedule constraints (successive train stops and runs), synchronization constraints (passenger transfers, crew transfers, rolling stock connections), and infrastructure constraints (minimum headway at open tracks and conflicting routes). Note that an infrastructure constraint also defines a train order at a conflicting route: event j initiates a process that must be finished before event i is allowed to pass the conflict point. Infrastructure constraints may occur between approaching or departing trains at a station, between approaching and departing trains with conflicting inbound and outbound routes, a fast train overtaking a slow train, opposite trains meeting at sidings on single-track lines, etc. The minimum headway time depends on the safety and signalling systems and may also contain running times to the conflict point (Goverde, 2005). Without loss of generality, we may assume that all these minimum headway times are nonnegative. 2.3. Max-plus linear systems Assuming that an event occurs immediately after all conditions are satisfied we obtain the equations

  xi ðkÞ ¼ max maxðaij þ xj ðk  lij ÞÞ; di ðkÞ ; j

i ¼ 1; . . . ; n;

ð4Þ

where j ranges over the predecessors of i, and n is the total number of events. Hence, each event waits until the last preceding process is completed and the scheduled event time has been reached. If the last preceding process has finished before the scheduled event time then the event occurs at the scheduled event time, while otherwise the event occurs as soon as possible with some delay. The max operations make (4) nonlinear equations in the discrete event times xi ðkÞ. It can however be written and analysed as a dynamic linear system in the so-called max-plus algebra. Max-plus algebra is the set Rmax ¼ R [ feg of real numbers extended with e :¼ 1 on which addition and multiplication operations are defined as a  b :¼ maxða; bÞ and a  b :¼ a þ b for all a; b 2 Rmax . Note that in max-plus algebra e is the neutral element of addition: a  e ¼ e  a ¼ a and e :¼ 0 is the neutral element of multiplication: a  e ¼ e  a ¼ a for all

272

R.M.P. Goverde / Transportation Research Part C 18 (2010) 269–287

a 2 Rmax . Let aij ¼ e if there is no precedence relation from j to i. Then the nonlinear Eq. (4) are formulated in max-plus algebra as a system of recursive linear equations n

xi ðkÞ ¼ aðaij  xj ðk  lij ÞÞ  di ðkÞ;

i ¼ 1; . . . ; n;

ð5Þ

j¼1 n

where aj¼1 aj ¼ a1  . . .  an denotes repeated ‘addition’. Note the resemblance with (4). The max-plus operations are extended to matrix addition and multiplication in a standard way:

½A  Bij :¼ aij  bij ¼ maxðaij ; bij Þ n

½A  Bij :¼ a aik  bkj ¼ max ðaik þ bkj Þ k¼1

k¼1;...;n

nn for any two square matrices A; B 2 Rnn max . The set Rmax of n  n matrices with entries in Rmax and the above addition and multiplication is an algebraic structure called an idempotent semiring (or dioid). This means that the operations satisfy a number of algebraic properties (associativity, distributivity, neutral elements), but in comparison to a field a semiring may not have nn , which additive or multiplicative inverses. Instead, addition in max-plus algebra is idempotent, i.e., A  A ¼ A for all A 2 Rmax is a powerful property that compensates the lack of inverses. For more details on max-plus algebra, see e.g., Baccelli et al. (1992), Cuninghame-Green (1979), Heidergott et al. (2006). > We will now formulate the system in matrix notation. Define the event time vector xðkÞ ¼ ðx1 ðkÞ; . . . ; xn ðkÞÞ and timetable > vector dðkÞ ¼ ðd1 ðkÞ; . . . ; dn ðkÞÞ , with superscript > denoting the transpose, and denote the maximum period shift over all processes by p :¼ maxðj;iÞ lij . Now collect the process times associated to all processes with period shift l in a matrix Al and fill the remaining entries with e, i.e., define the matrices A0 ; . . . ; Ap by ½Alij ij ¼ aij for all processes ðj; iÞ and ½Al ij ¼ e otherwise. If there are parallel processes between the same events with the same period shift then only the maximum process time has to be added to the matrices. The system of recursive linear Eq. (5) can now be written in matrix form as

p

xðkÞ ¼ a Al  xðk  lÞ  dðkÞ;

k 2 N;

ð6Þ

l¼0

where the maximum period shift p 2 N is called the order of the system. Eq. (6) defines a discrete-event dynamic system (DEDS). A periodic timetable is defined as a linear system in max-plus algebra as well by

dðkÞ ¼ T  dðk  1Þ;

dð0Þ ¼ d0 ;

which has solution dðkÞ ¼ d0  T k , where a power in max-plus algebra denotes repeated ‘multiplication’, i.e., 0 0 0 di ðkÞ ¼ di  T k ¼ di þ k  T. Note that di denotes the ith entry of vector d0 . The general pth-order state-space equations contain (many) zero-order terms. These implicit terms can be eliminated from the dynamic model to obtain a purely recursive state-space realisation (Baccelli et al., 1992). However, we then also lose the possibility of taking into account initial delays in the first period. Moreover, this procedure aggregates the propagation of delays within a timetable period. We will therefore keep the zero-order terms in the max-plus model to correctly represent and compute the propagation of delays over all events, see also Goverde (2007). Note that the event time vector xðkÞ is on both sides of the Eq. (6), which seems to define an implicit system. However, if A0 is strictly lower triangular (or there is a perturbation of event numbering such that A0 is strictly lower triangular) then the system is recursive. So under this assumption, the successive event time vectors fxðkÞjk 2 Ng can be computed recursively by (6) for a given initial timetable d0 and initial conditions xðlÞ ¼ xl for 1  p 6 l 6 1. Example. As an example consider the small railway transportation network of Fig. 1. This network consists of two main stations and three train lines: line 1 serves the regional area around station S1 , line 2 circulates between the two stations, and line 3 serves the regional area around station S2 . Both stations have two arrival and departure events. The departures are numbered 1–4 and the arrivals 5–8, corresponding to the nodes 1–8 in Fig. 1. The shown arc weights are the minimum process times in minutes: the minimum running and minimum dwell times (black arcs), minimum transfer times (dark gray arcs), and minimum departure headway times (light gray arcs). The (italic) node weights are the scheduled event times of the initial timetable d0 ¼ ð31; 30; 0; 1; 22; 57; 27; 57Þ> . The timetable cycle time is T ¼ 60. Note that all scheduled running times have 1 min running time supplement. Finally, the period shifts are depicted as tokens on the arcs of processes that cross the hour with respect to the given timetable. For example, line 1 departs at 31 after the hour from station S1 (node 1)

Fig. 1. Example railway network.

R.M.P. Goverde / Transportation Research Part C 18 (2010) 269–287

273

and returns at 22 in the next hour (node 5). An arc without token implies a period shift 0. The period shifts are computed by (3) for the given timetable and scheduled process times. For instance, the scheduled running time of train line 1, arc (1, 5), is given by the minimum running time plus running time supplement a0ij ¼ 50 þ 1 ¼ 51, and therefore l51 ¼ ð51 þ 31  22Þ=60 ¼ 60=60 ¼ 1. The tokens can be viewed as the processes that are active just before the whole hour. So if the trains run according to schedule then at each whole hour a train runs on line 1, the trains of lines 2 and 3 dwell in station S2 , and passengers may transfer between the trains of line 1 and 2 in station S2 . The precedence constraints for the example network are given as follows:

x5 ðkÞ P 50 þ x1 ðk  1Þ; x6 ðkÞ P 26 þ x2 ðkÞ;

x1 ðkÞ P 2 þ x7 ðkÞ;

x1 ðkÞ P 1 þ x2 ðkÞ;

x1 ðkÞ P 2 þ x5 ðkÞ;

x3 ðkÞ P 2 þ x6 ðk  1Þ; x2 ðkÞ P 2 þ x5 ðkÞ;

x4 ðkÞ P 1 þ x3 ðkÞ:

x8 ðkÞ P 55 þ x4 ðkÞ;

x7 ðkÞ P 26 þ x3 ðkÞ;

x3 ðkÞ P 2 þ x8 ðk  1Þ;

x4 ðkÞ P 2 þ x8 ðk  1Þ;

x2 ðkÞ P 2 þ x7 ðkÞ;

x4 ðkÞ P 2 þ x6 ðk  1Þ;

The 1st array contains the successive schedule constraints for line 1 and line 3, the 2nd array the schedule constraints for line 2, the 3rd array the transfer constraints, and the last array the two infrastructure constraints. The state vector xðkÞ ¼ ðx1 ðkÞ; . . . ; x8 ðkÞÞ> contains both departure and arrival events. The order of the system is p ¼ 1, since the maximum period shift is 1. The max-plus linear system is now given by

xðkÞ ¼ A0  xðkÞ  A1  xðk  1Þ  dðkÞ;

ð7Þ

with

The matrices A0 and A1 correspond to the processes scheduled within an hour and crossing the hour, respectively, and are partitioned into four parts corresponding to the departure and arrival events. The lower-left part contains the running times, the upper-right part contains all dwell and transfer times from arrivals to departures, the upper-left part contains minimum headway times between departures, and the lower-right part contains minimum headway times between arrivals (none in this example). Minimum headway times between arrivals and departures can also be defined in the lower-left and upper-right parts. 2.4. Marked precedence graphs and sparse matrix representations For large-scale railway networks the matrices A0 ; . . . ; Ap of the max-plus linear system (6) are typically very large and sparse. An efficient sparse representation can however be used corresponding to precedence graphs associated to max-plus nn is a weighted digraph with nodes V ¼ f1; . . . ; ng matrices. A precedence graph GðAÞ ¼ ðV; E; wÞ associated to a matrix A 2 Rmax and arcs E ¼ fðj; iÞjaij – eg with arc weight wðj; iÞ ¼ aij (Baccelli et al., 1992). We can generalise a precedence graph of a matrix p to a marked precedence graph of matrices A0 ; . . . ; Ap associated to a max-plus linear system xðkÞ ¼ al¼0 Al  xðk  lÞ as follows. Definition (Marked precedence graph). A marked precedence graph GðA0 ; . . . ; Ap Þ ¼ ðV; E; l; wÞ associated to a max-plus linear p system xðkÞ ¼ al¼0 Al  xðk  lÞ with Al 2 Rnn max for l ¼ 0; . . . ; p is a biweighted digraph with nodes V ¼ f1; . . . ; ng and for each ½Al ij – e an arc ðj; iÞ with marking lij ¼ l 2 N0 and holding time wði; jÞ ¼ ½Al ij . An efficient representation of the max-plus linear system is obtained by an adjacency list. An adjacency list adjðGÞ of a digraph G is a list of all nodes, where each node contains a list of (links to) all outgoing arcs. The arcs are defined in an arc list of the marked precedence graph, i.e., a list of tuples

ðj; i; l; ½Al ij Þ;

ð8Þ

for each ½Al ij – e, or equivalently a list of all ðj; i; lij ; aij Þ. The successors of a node j can thus be found in the list adjðjÞ. A marked precedence graph has an intuitive graphical representation: the arcs are marked by tokens denoting the associated period shifts. The minimum process times (holding times) are added to the marked graph as arc weights. Note that an adjacency list allows multiple (or parallel) arcs from a node j to node i corresponding to multiple processes from event j to i with possible different period shifts and process times. If there are parallel arcs with the same period shift then only one arc with the maximum process time is of interest. So in the matrix notation only the maximum process time is specified. For a scheduled max-plus linear system (6) the marked precedence graph can be extended to GðA0 ; . . . ; Ap ; d0 Þ with additional node

274

R.M.P. Goverde / Transportation Research Part C 18 (2010) 269–287

weights corresponding to the initial timetable d0 . The node weights can be included in the adjacency list. We already saw an example of a marked precedence graph in Fig. 1. The marked precedence graph clearly visualises the precedence constraints and the interconnection structure they generate. Example. An arc list of the matrices A1 and A0 of the example max-plus linear system (7) corresponding to the marked precedence graph shown in Fig. 1 is

ð9Þ

Here, the arcs ðj; i; lij ; aij Þ are shown in columns instead of rows for presentational reasons only. The precedence graph GðA0 Þ can also be seen as the subgraph of the marked precedence graph GðA0 ; . . . ; Ap Þ with zero marking. If this graph is acyclic then A0 is strictly lower triangular (or there is a perturbation of event numbering such that A0 is strictly lower triangular) so that the system (6) is recursive and explicitly solvable. On the other hand, if GðA0 Þ contains a circuit then the system cannot be solved explicitly. Hence, in the following we assume that GðA0 Þ is acyclic. Note that if GðA0 Þ is not acyclic then there is a circuit of events and processes within one timetable period with no well-defined precedence order. 2.5. Timed event graphs An arc list with entries (8) also defines a timed event graph (TEG) of which the max-plus linear system (6) is a state-space realisation (Baccelli et al., 1992; Goverde, 2005). A TEG is a decision-free Petri net (Murata, 1989), see Fig. 2. It consists of transitions (depicted as bars) corresponding to the events, places (depicted as circles) corresponding to the processes, and arcs connecting each place with exactly one input and one output transition. Moreover, an initial marking and holding times are assigned to the places. The initial marking (depicted as tokens in the places) corresponds to the period shifts of the processes and can be interpreted as one or more active processes (at the beginning of each period), while the holding time of a place is the minimum process time that a token is occupied in the place. A place in a TEG can also be viewed as a (timed) marked arc. In fact, in the literature a TEG is also called a timed marked graph (Murata, 1989). A timed marked graph is drawn equivalently to the marked precedence graph as introduced in the previous section, although also the Petri net graphical notation is used. The big difference between the marked precedence graph and a timed event (or marked) graph is that the former is just a static graph, while in the latter also the dynamics of event occurrences and movements of tokens can be animated and analysed, see below. Hence, the initial marking in the definition of a timed event graph. The dynamic behaviour of a TEG is defined by the movements of tokens over the places governed by the occurrences of events obeying the following two-step firing rule: (1) a transition i is enabled (ready to occur) if each upstream place contains a token and the associated holding times have expired and (2) a firing (occurrence) of an enabled transition i removes one token from each upstream place and adds one token to each downstream place. A sequence of firings satisfying the firing rule is called a legal firing sequence (feasible event sequence). An animation of the movement of tokens (active processes) over the network is obtained by successively applying the firing rule to the TEG. As an example, Fig. 3 shows the example TEG just after the hour, when event three has fired. Note that the tokens of the input places of event three are removed and tokens are added to its output

Fig. 2. Timed event graph of the example network (initial marking).

Fig. 3. Timed event graph of the example network after event three occurred.

275

R.M.P. Goverde / Transportation Research Part C 18 (2010) 269–287

places. The execution of a TEG corresponds to the max-plus linear system where an event occurs (‘fires’) after the process time of the last predecessor event has finished. Hence, there is a one-to-one correspondence between TEGs and max-plus linear systems (Cohen et al., 1985). In fact, one iteration of the max-plus recursion (6) gives the next event time vector xðkÞ corresponding to a legal firing sequence in the TEG where all transitions fire once and the TEG has returned to the initial marking. In addition to the arcs (8) of the precedence relations also the timetable has to be included in the timed event graph. This can be done in (at least) two ways. The first is by transforming the timetable constraint xðkÞ P dðkÞ ¼ T k  d0 with d0 2 ½0; TÞn into precedence constraints of the form (2). Define a new ‘clock’ event 0 with as initial event time the beginning of the initial period x0 ð0Þ ¼ 0, and add the precedence constraints 0

x0 ðkÞ P T þ x0 ðk  1Þ and xi ðkÞ P di þ x0 ðkÞ; k

i ¼ 1; . . . ; n: 0

0

Then, x0 ðkÞ ¼ T ¼ k  T gives the beginning of the kth period and xi ðkÞ P di þ x0 ðkÞ ¼ di þ k  T equals the original constraint (1) for a periodic timetable. Hence, the timetable can be included into   the timed event graph by adding an extra ‘clock’ 0 transition 0, a ‘clock’ place ð0; 0; 1; TÞ, and the timetable places 0; i; 0; di for all i ¼ 1; . . . ; n (Goverde, 2005), see Fig. 4. A timetable can also be incorporated into the timed event graph by extending the firing rule. So instead of adding an extra event and n þ 1 places as described above, we extend firing rule (1) as follows: (1’) A transition i is enabled if each upstream place contains a token, the associated holding times have expired, and the scheduled event time has been reached. To execute such a scheduled timed event graph either the scheduled event times must be assigned to the transitions or just the initial scheduled event time with an updating scheme for the next scheduled event time after each firing. The latter implies also an extension of the 2nd firing rule. Assuming a periodic timetable with cycle time T this becomes: (2’) A firing of an enabled transition i removes one token from each upstream place, adds one token to each downstream place, and updates the 0 di þ T. The timed event graph is initialised with di ¼ di , i ¼ 1; . . . ; n. The algorithm of Secnext scheduled event time as di tion 4 is an implementation of this extended firing rule with a legal firing sequence (event list) derived in Section 3.3. A scheduled timed event graph thus has two additional parameters: the cycle time T and the initial timetable vector d0 2 ½0; TÞn which assigns weights to the transitions. In the sequel, a scheduled timed event graph (or marked precedence graph) is represented as GðA0 ; . . . ; Ap ; d0 ; TÞ ¼ ðV; E; l; w; d0 ; TÞ with ðV; E; l; wÞ defined by the adjacency list adjðGÞ linking to the arc list (8). A timed event graph is called live if each transition can fire from the initial marking no matter what firing sequence is chosen (Murata, 1989). A timed event graph that is not live is called deadlocked. A necessary and sufficient condition for a live TEG is that the subgraph GðA0 Þ of the places with zero initial marking is acyclic (Murata, 1989). Suppose there is a zero-token circuit in a timed event graph. By definition, a transition i1 on this circuit can only fire if all preceding places contain a token. Since i1 lies on the zero-token circuit there is a preceding place ði2 ; i1 Þ without token. Hence, the predecessor i2 must fire first before i1 is enabled. But i2 has the same situation, and ultimately there is a transition iq on the circuit with preceding zero-token place ðiq ; i1 Þ requiring that i1 has to fire before iq . Hence, none of the transitions on a zero-token circuit will ever be enabled and fire, and eventually any transition accessible from this circuit will never be enabled again: they are deadlocked. Hence, a live TEG corresponds to a recursive max-plus linear system (6). 3. Delay propagation 3.1. The max-plus delay propagation model Assume that at some reference time t0 (for instance the current time) we know or have a prediction of the delays of all events scheduled up to t0 . Without loss of generality we assume that t 0 occurs in the first period of the max-plus model, t 0 2 ½T; 2TÞ. Then the initial conditions to the max-plus linear system (6) are the event time vectors of the last p periods

Fig. 4. Scheduled timed event graph of the example network, the gray part shows the clock transition (event 0) and clock place (T ¼ 60), and the timetable places to each transition.

276

R.M.P. Goverde / Transportation Research Part C 18 (2010) 269–287

x1p ; . . . ; x0 2 Rnmax such that xðlÞ ¼ xl for l ¼ 1  p; . . . ; 0, and an initial vector x1 2 Rnmax that partially determines the state xð1Þ  > for the events scheduled up to t 0 . Writing x1 ¼ x11 ; . . . ; x1n , we have xi ð1Þ ¼ x1i for those events 1 6 i 6 n that have a scheduled event time di ð1Þ 2 ½T; t0 , since these event times are known at time t 0 . The event times of the remaining events i are still unknown at time t0 and depend on both the delayed events of preceding processes that were initiated before time t0 as well as on the delay propagation within the first period over ðt0 ; 2TÞ. Since we have no information on the event times scheduled after t0 they are given an initial estimate x1i ¼ di ð1Þ, corresponding to the scheduled event times, and thus xi ð1Þ P x1i . Thus, we denote by xð1 Þ ¼ x1 the partially known state vector of the first period at time t0 . Let d0 be a finite initial timetable vector with all entries falling in a timetable period of length T and let zðkÞ ¼ xðkÞ  dðkÞ be an output vector containing the delays in period k. Then the delay propagation model can be written as the scheduled max-plus linear system

8 p > > > xðkÞ ¼ a Al  xðk  lÞ  dðkÞ > > > l¼0 < dðkÞ ¼ T  dðk  1Þ > > > zðkÞ ¼ xðkÞ  dðkÞ > > > : xð1 Þ ¼ x1 ; xð0Þ ¼ x0 ; . . . ; xð1  pÞ ¼ x1p ; dð0Þ ¼ d0 :

ð10Þ

Note that the first dynamic equation implies xðkÞ P dðkÞ and therefore zi ðkÞ P 0 for all 1 6 i 6 n. The initial conditions can be given in terms of initial delay vectors z1 ; z0 ; . . . ; z1p 2 RnP0 of delays that reach events in the first (current) timetable period or beyond. The initial state vectors are then computed as xl ¼ dðlÞ þ zl with dðlÞ ¼ d0 þ l  T for 1  p 6 l 6 1. The delay propagation of any given initial delay scenario can be computed from (10) using matrix–vector multiplications in max-plus algebra (Goverde, 2007). For large-scale sparse networks, however, this direct procedure requires high memory usage and a huge amount of unnecessary computations. In Section 4 we therefore give an efficient graph algorithm with a running time that is bounded by a constant times the number of delayed events, which is the best one can wish for. The next section first deals with structural delays due to unrealisable (or unfeasible) process times. 3.2. Realisability and structural delays A scheduled process time a0ij consists of a minimum (or nominal) process time aij plus a process time supplement or buffer time sij , i.e., a0ij ¼ aij þ sij . Time supplements are used to recover from existing train delays while buffer times prevent or reduce secondary delays propagated from one train to another due to conflicting routes or waiting at stations to secure connections. However, a minimum process time can also be larger than the scheduled process time, i.e., sij < 0. This happens for instance for short stops where the minimum dwell time plus supplement is rounded down to minutes so that trains do not have to wait unnecessary on their scheduled departure time. The same applies to trains passing stations (and junctions, etc.) without dwelling: the running time supplement is spread even over the train run or at the end and intermediate passing times are rounded down to avoid wasting recovery time. Also timetable flaws cause unrealisable processes when process times are estimated too low or standard norms are applied that are not valid everywhere. Temporary speed restrictions and local route changes may also cause larger running times than originally scheduled, and a delay propagation model could be applied to analyse the impact of these changed operational conditions on the delay propagation. Note that the unrealisable processes explain that a0ij is used in the period shift formula (3) instead of aij : otherwise extra trains would be added incorrectly to the model to guarantee feasibility. The used minimum process times may depend on the application of the model. For delay forecasting, the minimum process times aij must be representative process times. Then a time supplement or buffer time is available for delay recovery if aij < a0ij , and otherwise an existing delay is increased or a primary delay is generated when the train is not too early. For timetable evaluation during the timetable design different objectives may ask for different choices, such as the scheduled process times, calculated process times, standard norms, statistics from realization data, or a mixture of these. Let A0l be the matrix of the scheduled process times associated to Al for 0 6 l 6 p. A timetable fdðkÞjk 2 N0 g or the corresponding scheduled max-plus linear system (6) is called realisable if all (minimum) process times are smaller or equal to their scheduled process time, i.e., Al 6 A0l for 0 6 l 6 p. A necessary and sufficient condition for a realisable timetable is p

dðkÞ P a Al  dðk  lÞ for all k 2 N;

ð11Þ

l¼0

i.e., di ðkÞ P ½Al ij  dj ðk  lÞ for all i; j; l; k. Hence, if there are no initial delays and the scheduled max-plus linear system is realisable then no (primary) delays will be generated, i.e., xðkÞ ¼ dðkÞ for all k 2 N. In the case of a periodic timetable dðkÞ ¼ d0  T k Eq. (11) becomes (Goverde, 2007) p

d0 P a Al  T l  d0 :

ð12Þ

l¼0

which is independent from k. In a realisable periodic timetable each process thus has a nonnegative slack time 0

0

sij ¼ di  dj  ½Al ij þ l  T:

ð13Þ

R.M.P. Goverde / Transportation Research Part C 18 (2010) 269–287

277

If event j has a delay z > 0 in some period k then this delay is propagated to each successor event i for which the slack time sij < z. Event i then receives a delay z  sij from event j. Note however that event i may be delayed more if there is another path from j to i with less total slack time, or if it receives a larger delay from another predecessor event. In a realisable timetable a delay can grow by propagation of another dominating delay. It can however be proven that delays do not increase over intervals of p periods, with p the order of the system. Denote the maximum entry in a vector z 2 Rnmax as n

kzk :¼ a zi : i¼1

Lemma 1. Let (10) be a realisable scheduled max-plus linear system. Then the maximum delay is nonincreasing over intervals of p periods, i.e., p

kzðkÞk 6 a kzðk  lÞk

for all k 2 N:

l¼0

Moreover, if there exists an integer ks 2 N such that zðks þ lÞ ¼ 0 for all 1 6 l 6 p then zðkÞ ¼ 0 for all k P ks þ 1. Proof. Recall zðkÞ ¼ xðkÞ  dðkÞ P 0 and by realisability ½Al ij 6 di ðkÞ  dj ðk  lÞ. Then n

kzðkÞk ¼ aðxi ðkÞ  di ðkÞÞ; i¼1

!  p n  ¼ a a a ð½Al ij þ xj ðk  lÞÞ  di ðkÞ  di ðkÞ ; n

i¼1

l¼0 j¼1 p

n

!

n

6 a a aððdi ðkÞ  dj ðk  lÞ þ xj ðk  lÞÞ  di ðkÞÞ  di ðkÞ ; i¼1

l¼0 j¼1 p

n

n

¼ a a aðxj ðk  lÞ  dj ðk  lÞÞ  0; i¼1 l¼0 j¼1 p

¼ a kzðk  lÞk : l¼0

It follows immediately that after an interval of p periods without delays there will be no more delays ever after. Note that for the particular case of p ¼ 1 the sequence fkzðkÞk jk 2 Ng is nonincreasing. h If the scheduled system is unrealisable then there are processes with slack time sij < 0, i.e., aij > a0ij , which always will generate a primary delay to the downstream event. These structural primary delays are computed as p

^z0 ¼

a Al  T l  d0  d0

!  0;

ð14Þ

l¼0

    0 0 where 0 ¼ ð0; . . . ; 0Þ> 2 Rnmax , i.e., ^z0i ¼ max maxðj;i;lÞ dj þ ½Al ij  l  T  di ; 0 P 0. The structural primary delays ^z0 may also propagate further over the network and cause structural secondary delays to other trains as well. The (primary and secondary) structural delays will occur repetitively in each period independent of other delays. All (primary and secondary) structural delays can be computed by applying the delay propagation model (10) with the structural primary delays ^z0 as initial conditions, i.e., zl ¼ ^z0 for all 1  p 6 l 6 1. Note that for unrealisable timetables ks ¼ 1 in Lemma 1, unless the structural delays are neglected. We come back to the structural delays in Section 5.2. 3.3. Feasible event list The adjacency list representation of the marked precedence graph associated to the max-plus linear system (10) can be used to recursively compute the delay propagation of the initial delays, similar to the execution of a timed event graph. The efficiency of this procedure depends on how the next event is selected. Moreover, the output is the sequence of delay vectors fzðkÞjk 2 Ng, or just the sequence of delayed events fzi ðkÞ > 0j1 6 i 6 n; k 2 Ng (note that any event not in this list has zero delay). To determine this output we can discard events that are not delayed and have no adjacent unrealisable process time, since these events will not propagate or generate a delay. In this section, we derive the structure of a (dynamic) event list that defines the successive events to be scanned for delay generation or propagation. The event list is feasible if all delays can be computed by passing the event list once. The event list is updated during the algorithm and therefore dynamic. The delay propagation model (10) is recursively solvable iff GðA0 Þ is acyclic as explained from different perspectives in Sections 2.3, 2.4 and 2.5. In this paper the initial marking (period shift) is computed using (3). So this formula defines the 0 0 matrix A0 . In particular, any process ðj; iÞ with scheduled process time a0ij ¼ 0 has scheduled event times dj ¼ di , and there0 0 fore a period shift lij ¼ 0 by (3). Conversely, any process ðj; iÞ with scheduled event times dj ¼ di and a period shift lij ¼ 0

278

R.M.P. Goverde / Transportation Research Part C 18 (2010) 269–287

has scheduled process time a0ij ¼ 0 by (3). Denote by A00 the matrix of scheduled process times for processes with zero period     shift. Note that GðA0 Þ and G A00 only differ in the arc weights. Now, define the subgraph G0 ðA00 Þ ¼ ðV 0 ; E0 Þ # G A00 as

o n o n 0 0 E0 ¼ ðj; iÞja0ij ¼ 0 ¼ ðj; iÞjlij ¼ 0; di ¼ dj ; V 0 ¼ f1 6 i 6 nj9j : ði; jÞ 2 E0 or ðj; iÞ 2 E0 g:

  The graph G0 A00 has only arcs with arc weight zero, has no isolated nodes (i.e., nodes without predecessor or successor),   and each path consists of nodes sharing a common scheduled event time. The graph G0 ðA0 Þ has the same structure of G0 A00 except maybe for the arc weights. In particular, they are equal if A0 contains realisable process times only. This graph plays a central role in the event list description. First, we prove the following lemma. Lemma 2. Consider the max-plus linear system (10). Then GðA0 Þ is acyclic if and only if G0 ðA0 Þ is acyclic. Proof. We give a constructive proof that if there is a circuit in GðA0 Þ then it satisfies the conditions of G0 ðA00 Þ and is therefore contained in both, since G0 ðA00 Þ # GðA00 Þ. First, consider a circuit n in GðA0 ; . . . ; Ap Þ with at least one arc ðj; iÞ 2 n with a0ij > 0. P Then also the circuit weight wðnÞ ¼ ðj;iÞ2n a0ij > 0 since all process times are assumed nonnegative. Then

X ðj;iÞ2n

lij ¼

X a0ij þ d0j  d0i ðj;iÞ2n

T

¼

1 T

X ðj;iÞ2n

a0ij þ

!  X 0 wðnÞ 0 ¼ dj  di > 0: T ðj;iÞ2n

Here, the fact that the scheduled event times on a circuit cancel each other out: if n ¼ ði1 ; . . . ; iq Þ then  we used  0 0 0 0 0 0 0 0 dj  di ¼ di1  di2 þ di2  di3 þ . . . þ diq  di1 ¼ 0. Hence, a necessary condition for a zero-marking circuit is that all 0 0 arcs have zero process time. Now assume wðnÞ ¼ 0. Then a0ij ¼ 0 for each ðj; iÞ 2 n, and therefore di ¼ dj for each ðj; iÞ 2 n, i.e.,  all scheduled event times on the circuit are equal. The period shift for each arc on n becomes      lij ¼ a0ij þ d0j  d0i =T ¼ 0=T ¼ 0. Hence, G0 A00 is exactly the subgraph of G A00 satisfying the conditions for a (zero-marking) circuit. It follows that GðA0 Þ has a (zero-marking) circuit iff G0 ðA0 Þ has a circuit, and therefore also GðA0 Þ is acyclic iff G0 ðA0 Þ is acyclic. h P

ðj;iÞ2n

The following theorem is the main theoretical result of this paper and the basis of the developed algorithm. Visiting the nodes of an acyclic digraph in topological order means that for each arc ðj; iÞ node j is visited before node i (Aho et al., 1983). For example, the arc list (9) of the marked precedence graph in is topologically ordered with respect to GðA0 Þ. Theorem 1. Let (10) be a realisable max-plus linear system with acyclic G0 ðA0 Þ ¼ ðV 0 ; E0 Þ. Then a feasible event list for computing the delay propagation fzðkÞjk 2 Ng is a list of delayed events only, ordered by scheduled event time, and breaking ties in topological order for delayed events i 2 V 0 and randomly otherwise. Proof. By definition of a realisable timetable (11), a punctual event cannot generate a delay to successors. Conversely, only a delayed event can propagate a delay to successors. Hence, we can discard all events that are not reached by a delay and only have to scan delayed events. We now prove feasibility of the event ordering. If V 0 – £ then G0 ðA0 Þ is an acyclic digraph by assumption, and thus a topological order of the nodes in V 0 exists. Let i 2 f1; . . . ; ng be an arbitrary event and assume first that i R V 0 . The event time xi ðkÞ can be computed in any period k 2 N only if the event times of all predecessors are available. If di ðkÞ is the associated scheduled event time then each predecessor j associated to some marked arc ðj; i; lÞ has scheduled event time dj ðk  lÞ 6 di ðkÞ. Moreover, from (3) it follows that if di ðkÞ ¼ dj ðkÞ then either a0ij > 0 and therefore l ¼ lij > 0, or a0ij ¼ 0 and so i; j 2 V 0 . Hence, for i R V 0 we must have dj ðk  lÞ < di ðkÞ and the scheduled order equals the precedence order. If for some j 2 f1; . . . ; ng we have dj ðkÞ ¼ di ðkÞ then j cannot be a predecessor of i and vice versa since this would violate the condition dj ðk  lÞ < di ðkÞ. Hence, all events with equal scheduled event time but not in V 0 can be computed independently and thus in random order. Now suppose i 2 V 0 . Then there may exist a predecessor j 2 V 0 with a0ij ¼ lij ¼ 0 and dj ðkÞ ¼ di ðkÞ. If so, then the event time xj ðkÞ of this predecessor must be computed before xi ðkÞ can be computed. Analogously, since also j 2 V 0 the event time of any predecessor s 2 V 0 , if one exists, must be computed prior to xj ðkÞ, etc. Hence, by visiting the events in V 0 upstream of i, if any, in topological order the correct precedence order is obtained for computing the event times one at a time. h   G0 A00 can be constructed by going through the adjacency list of GðA00 Þ once and adding arcs with a0ij ¼  simply   0 to the (initially empty) arc list of G A00 , which takes OðmÞ time for m arcs. This assumes however that an arc list for G A00 is available or that the arc listof GðA0 Þ contains the scheduled process times as additional arc weights. Another approach is to con0 0 struct G0 ðA0 Þ (not G A00 ) by passing through the adjacency list of GðA0 Þ once, checking for each arc with lij ¼ 0 if di ¼ dj , and if so adding the arc to the adjacency list of G0 ðA0 Þ. If d0 is included in the adjacency list or implemented as an array indexed by event number then this procedure also takes OðmÞ time. This graph extraction procedure either returns an adjacency list of G0 ðA0 Þ or returns an empty list denoting E0 ¼ £. A necessary and sufficient condition for E0 ¼ £, and thus V 0 ¼ £, is that A00 has only positive entries a0ij > 0. n o n o 0 0 Corollary 1. Let (10) be a realisable scheduled max-plus linear system with E0 ¼ ðj; iÞja0ij ¼ 0 ¼ ðj; iÞjlij ¼ 0; di ¼ dj ¼ £. Then a feasible event list for computing fzðkÞjk 2 Ng is a list of delayed events only, ordered by scheduled event time and breaking ties in random order.

R.M.P. Goverde / Transportation Research Part C 18 (2010) 269–287

279

We may renumber the nodes of the acyclic digraph G0 ðA0 Þ in topological order such that j < i for all arcs ðj; iÞ 2 E0 . Topological sorting of an acyclic digraph can be done in linear Oðn þ mÞ time using depth-first search (Aho et al., 1983), where n and m are the number of nodes and arcs in the acyclic graph, respectively. In the following we therefore assume that the nodes of the subgraph G0 ðA0 Þ are topologically ordered. Note that if GðA0 Þ is topologically ordered then so is G0 ðA0 Þ and so the nodes may also be topologically ordered with respect to GðA0 Þ, although GðA0 Þ is typically considerably larger than G0 ðA0 Þ. Corollary 2. Let (10) be a realisable scheduled max-plus linear system with acyclic G0 ðA0 Þ and events topologically ordered with respect to G0 ðA0 Þ. Then a feasible event list for computing fzðkÞjk 2 Ng is a list of delayed events only, ordered by scheduled event time and breaking ties in ascending event index order. In the following we assume without loss of generality that either E0 ¼ £ or that the events are topologically sorted with respect to G0 ðA0 Þ. 4. A bucket-based delay propagation algorithm Theorem 1 and its corollaries suggest a graph algorithm where the event list is defined by the delayed events maintained in buckets of equal scheduled event time. Here, a bucket Bt is a linked list (of variable size) of all delayed events ði; k; zi ðkÞÞ with scheduled event time di ðkÞ ¼ t. The buckets themselves are elements of an array B ¼ ðBt Þ, or more precisely, B is an array of (headers for) lists indexed by scheduled event time. A new delayed event ði; k; zi ðkÞÞ with scheduled event time di ðkÞ is inserted in bucket Bdi ðkÞ and an existing delayed event can be found by searching for it in the bucket indexed by the associated scheduled event time. This approach of sorting events is known as bucket sort (Aho et al., 1983). We here tacitly assume that the scheduled event times are specified in whole minutes, which is generally true for railway timetables. If the timetable precision is higher, say a tenth of a minute or a second, then the method is easily adapted. If the indices of array B are restricted to nonnegative integers then it is convenient to choose the initial timetable vector d0 such that 0 d0 2 ½ðp  1Þ  T; p  TÞn , where p is the order of the max-plus linear system. Then di ð1  pÞ ¼ di þ ð1  pÞ  T 2 ½0; TÞ for all 1 6 i 6 n and thus any delayed event ði; k; zi ðkÞÞ, k P 1  p, has at least scheduled event time 0. Algorithm 1 gives the pseudocode for computing the propagation of any initial delay scenario in a scheduled max-plus linear system. The required input is an adjacency list of the associated marked precedence graph (or TEG), the initial timetable vector d0 2 ½ðp  1Þ  T; p  TÞn , the cycle time T, and a list of initial delays Z 0 . We assume that the initial delayed events are specified by tuples of the form ði; k; zi ðkÞÞ 2 Z 0 , where

Z 0 ¼ fði; k; zi ðkÞÞjzi ðkÞ > 0; 1  p 6 k 6 1g:

ð15Þ

Here i 2 f1; . . . ; ng is the delayed event number, k 2 Z the period in which the delay occurs, and zi ðkÞ > 0 the delay. The output is a list of all delays

Z ¼ fði; k; zi ðkÞÞjzi ðkÞ > 0; 1 6 i 6 n; k P 1  pg: Note that the delay list can be implemented or returned as a sequence of vectors (or arrays) fzðkÞjk P 1  pg, with zi ðkÞ ¼ 0 for punctual events. Algorithm 1. Delay propagation Input: Adjacency list of GðA0 ; . . . ; Ap Þ, timetable vector d0 , cycle time T, initial delay list Z 0 . Output: Delay list Z. 1 t 1;dmax 1; /*Initialisation */ 2 foreach ði; k; zi ðkÞÞ 2 Z 0 do 0

3 d di þ k  T; 4 Bd Bd [ fði; k; zi ðkÞÞg; 5 t minðt; dÞ; dmax maxðdmax ; dÞ; 6 while t 6 dmax 7 while Bt – £ do 8 ðj; k; zj ðkÞÞ arg minfjjðj; k; zj ðkÞÞ 2 Bt g; 9 Bt Bt n fðj; k; zj ðkÞÞg; Z Z [ fðj; k; zj ðkÞÞg; 10 foreach ðj; i; l; ½Al ij Þ 2 adjðjÞ do 11 12

d z

/*Main loop*/ /*Inner loop*/ /*Delay selection*/ /*Scan*/

0

di þ ðk þ lÞ  T; t þ ½Al ij þ zj ðkÞ  d;

13 if z > 0 then d; Bd fi; k þ l; zg; 14 if d > dmax then dmax 15 else if ði; k þ l; Þ 2 Bd then zi ðk þ lÞ maxðzi ðk þ lÞ; zÞ; 16 else Bd Bd [ fði; k þ l; zÞg; 17 t t+1; 18 return Z;

/*propagate*/

/*Terminate*/

280

R.M.P. Goverde / Transportation Research Part C 18 (2010) 269–287

Dutch railway timetable 2007 250

Count

200

150

100

50

0

0

5

10

15

20

25

30

35

40

45

50

55

60

Scheduled event time (min) Fig. 5. Histogram of the scheduled event times in the Dutch basic hour railway timetable 2007.

Algorithm 1 uses two counters t and dmax denoting the current scheduled event time and the maximum scheduled event time of delayed events encountered so far, respectively. Initially, t (resp. dmax ) is set to a very large (resp. negative) value (line 1). Then each initial delay is sorted in buckets and the counters t and dmax are updated according to the scheduled event times of the sorted delayed events (lines 2–5). At the end of the initialisation t is the smallest scheduled event time amongst the initial delayed events, Bt is the nonempty bucket of smallest index, and Bdmax is the nonempty bucket with largest index. The main program (lines 6–17) contains an outer and an inner loop. The outer loop iterates over the counter t and stops as soon as t exceeds dmax which implies that all buckets are empty. The algorithm then terminates by returning the computed delay list Z, which contains all initial and secondary delayed events in scheduled order (line 18). For each current value of the counter t the inner loop (lines 7–16) is called, which explores the delayed events within bucket Bt . If bucket Bt is empty the program returns to the outer loop, increases the current scheduled event time t by one and starts a next iteration. On the other hand, if Bt is nonempty then the event with minimal index j is selected (line 8), deleted from the bucket and added to the final list of delays (line 9). Lines 10–16 scan each outgoing arc of j. If the current delay zj ðkÞ exceeds the slack on an outgoing arc then the delay is propagated to the successor event i and one of three options occurs: if the scheduled event time di ðkÞ is larger than any scheduled event time encountered so far then dmax is updated, the bucket array B is expanded with empty buckets up to Bdi ðkÞ , and the new (reduced) delay is inserted into Bdi ðkÞ (line 14). Otherwise, if the successor i is already in bucket Bdi ðkÞ then the delay is updated (line 15), and if i is not yet in the bucket then it is inserted (line 16). If the adjacency list of event j is exhausted, a new iteration of the inner-loop starts by finding the next event in bucket Bt or establishing that the bucket has become empty. 4.1. Data structures The bucket implementation is particular efficient if the scheduled event times are evenly distributed over a finite discrete set (in each period). In general, railway timetables of large-scale networks show a rather discrete uniform distribution of scheduled event times over the integers from 0 to 59 because of the diverse running times over the various stretches between station stops. Fig. 5 shows a histogram of the hourly scheduled event times in the national Dutch railway timetable of the year 2007. This periodic timetable contains 11,401 integer scheduled event times ranging from 0 to 59 minutes, with 2779 departure times, 2779 arrival times, and 5843 passing times. The histogram indeed suggests that the events are spread evenly over the integers ½0; 59. The Chi-Square goodness-of-fit test with hypothesised discrete uniform distribution U½0; 59 confirms the good fit (v2 ¼ 6:5099, 11 d.f., p ¼ 0:8373). Fig. 5 also shows that in this case the size of a bucket will never exceed 245 events (for scheduled event time 38), where this upper bound corresponds to the (rare) case that in some period all events over the network with scheduled event time 38 are delayed. The bucket implementation is therefore an effective way to decompose the list of delayed events during the course of the algorithm for quick access. The efficiency of Algorithm 1 also depends on an appropriate data structure of the buckets. In particular, each bucket should be implemented as an ordered list with respect to increasing event number or more generally as a priority queue (Aho et al., 1983). An ordered list implementation of the buckets improves finding, updating, and inserting items in a bucket as opposed to an unordered list regardless of the event selection method (line 8). Hence, even if events in a bucket may be selected randomly the ordered list implementation of the buckets is the preferred data structure. The buckets have a dynamic size and the maximum size differs from bucket to bucket. Most buckets contain only a few events and even ‘full’ buckets have limited size as discussed above. Furthermore, the total number of events and therefore the range of the event numbering can be very large. Given these characteristics, a linked ordered (or linear) list is the preferred data structure for each bucket (Knuth, 1998). In a linked ordered list finding the minimal element in a bucket Bt

R.M.P. Goverde / Transportation Research Part C 18 (2010) 269–287

281

(cf. line 8) just reduces to returning the first element in Bt , which takes Oð1Þ time, and so does deleting this minimal element from the ordered list (line 9). Finding a given event by its number (line 15) or establishing that the event is not yet in the list takes Oðnt Þ time, where nt is the maximal size of bucket t. When an event has been located in a bucket, updating the delay takes Oð1Þ time. On the other hand, as soon as the list search meets a higher numbered event it is known that the current event is not yet in the list and so it is inserted between the located higher numbered event and its predecessor (line 16) in Oð1Þ time. In practice, nt is a reasonable small number that can be bounded by a constant. For a uniform discrete distribution of scheduled event times U½0; T  1 the average maximum size of the buckets is dn=Te. For the Dutch railway timetable 2007 the average maximum bucket size is d11401=60e ¼ 191, while the maximum bucket size is Bt ¼ 245 for t ¼ 38 ðmod 60Þ. However, in practice the bucket sizes are typically much smaller, since (1) only areas reached by the initial delays are considered and (2) punctual events are not added to the buckets. For the Dutch railway timetables the linked ordered list implementation performs well. The amount of work per delayed event j is therefore bounded by a constant times the number of successors of node j. Moreover, the number of successors of each event (the outdegree) is also typically limited to a small number (<10). Hence, the computation time of the algorithm is bounded by a constant times the number of scanned delayed events.

5. Timetable properties and their impact on delay propagation and algorithm convergence 5.1. Realisable and stable timetables The convergence of Algorithm 1 depends on the stability of the scheduled railway system. A scheduled max-plus linear system is stable if any initial delay in fz1p ; . . . ; z1 g settles after a finite number of periods. Formally, a scheduled max-plus linear system is stable if for any fz1p ; . . . ; z1 g there exists a ks ðz1p ; . . . ; z1 Þ 2 N such that zðkÞ ¼ 0 for all k P ks þ 1. The minimal such integer ks is called the settling time associated to the initial delays fz1p ; . . . ; z1 g, measured in number of periods. Assume p P 2 and consider any process ðj; iÞ with period shift lij ¼ p P 2. If event j is delayed in period k by more than the slack time, zj ðkÞ > sij , then the successor i will be delayed in period k þ p, zi ðk þ pÞ > 0. However, if the slack times of all other processes are sufficient to settle the delays up to period k, then zðk þ lÞ ¼ 0 for l ¼ 1; . . . ; p  1. Hence, reaching a delay vector zðkÞ ¼ 0 does not necessarily mean that all delays have settled unless p ¼ 1. From Lemma 1 it follows that the settling time (or settling period) ks is determined from (10) as

ks ¼ minfk 2 Njxðk þ lÞ ¼ dðk þ lÞ; 1 6 l 6 pg:

ð16Þ

Note that ks is defined as the last period with zðks Þ > 0, so that all delays have settled after ks periods. Applying the recursive delay propagation model (10) directly thus requires p  1 more iterations after reaching a zero delay vector to be sure that all delays have settled. This is of course not needed in delay propagation Algorithm 1, since this algorithm looks ahead by maintaining a delay list with delayed events that still have to be scanned. As soon as this list is empty there are and will be no more (unstructural) delays and the period of the last delayed event is the settling time. A necessary and sufficient condition for a scheduled max-plus linear system to be stable is that each circuit in the marked precedence graph (or TEG) GðA0 ; . . . ; Ap Þ has some positive slack time available for delay recovery. Or equivalently, the system is stable if the maximum cycle mean k is smaller than T (Braker, 1993; Goverde, 2007). Here, a cycle mean is defined for any circuit in a marked precedence graph (or TEG) as the ratio of the sum of holding times and the marking of the circuit. The maximum cycle mean is defined as

k ¼ max n2C

wðnÞ ; lðnÞ

ð17Þ

P P where C is the set of all elementary circuits in GðA0 ; . . . ; Ap Þ; wðnÞ ¼ ðj;iÞ2n aij is the circuit weight, and lðnÞ ¼ ðj;iÞ2n lij is the circuit marking. This k equals the largest (generalised) eigenvalue of the state matrices in max-plus algebra, which is the p solution k to the (generalised) max-plus eigenvalue problem al¼0 Al  kl  v ¼ v , or maxðj;i;lÞ ð½Al ij  l  k þ v j Þ ¼ v i for all 1 6 i 6 n (Baccelli et al., 1992; Goverde, 2005; Goverde, 2007). For the first-order representation xðkÞ ¼ A  xðk  1Þ  dðkÞ the eigenvalue problem looks more familiar as A  v ¼ k  v . The eigenvector v 2 Rnmax represents a realisable initial timetable allowing periods with minimal cycle time k (Braker, 1993; Goverde, 2007). Efficient algorithms are available for solving the (generalised) eigenvalue problem, and in particular graph algorithms based on the maximum cycle mean (17) (CochetTerrasson et al., 1998; Dasdan, 2004). A scheduled max-plus linear system is called stable if k < T, critical if k ¼ T, and unstable if k > T. Stability is concerned with the internal delay propagation behaviour over circuits in the network. In a stable timetable all circuits contain some positive slack time and therefore each delay is decreasing over any circuit, although it could be dominated by a propagated delay from another circuit. The next lemma shows that the maximum delay in a stable timetable is strictly decreasing over circuit traversals until it settles.

282

R.M.P. Goverde / Transportation Research Part C 18 (2010) 269–287

Lemma 3. Let (10) be a stable scheduled max-plus linear system with acyclic GðA0 Þ, C the set of all elementary circuits in GðA0 ; . . . ; Ap Þ, and zðkÞ > 0. If the timetable dðkÞ ¼ d0  T k is also realisable then the maximum delay is decreasing over circuit traversals,

kzðkÞk > a kzðk þ lðnÞÞk :

ð18Þ

n2C

~ ~0  T k such that (18) holds If the timetable dðkÞ ¼ d0  T k is unrealisable then there is a realisable (reference) timetable dðkÞ ¼d ~ true with zðkÞ ¼ xðkÞ  dðkÞ. Proof. Assume first that dðkÞ ¼ d0  T k is realisable. Let zðkÞ > 0 and C be the set of all elementary circuits in GðA0 ; . . . ; Ap Þ. Then

a kzðk þ lðnÞÞk ¼ a a xi ðk þ lðnÞÞ  di ðk þ lðnÞÞ n2C i2n

n2C

¼ a a ððxi ðkÞ þ wðnÞÞ  di ðk þ lðnÞÞÞ  di ðk þ lðnÞÞ n2C i2n

¼ a aðxi ðkÞ þ wðnÞ  di ðk þ lðnÞÞÞ  0 n2C i2n

¼ a a ðxi ðkÞ þ wðnÞ  di ðkÞ  lðnÞ  T Þ  0 n2C i2n

6 a a ðzi ðkÞ þ lðnÞ  k  lðnÞ  T Þ  0 n2C i2n n

< a a zi ðkÞ  0 ¼ a zi ðkÞ ¼ kzðkÞk : n2C i2n

i¼1

The inequality follows from (17) and the strict inequality follows from the stability condition k < T. So, the maximum delay is strictly decreasing after circuit traversals until it settles. Also delayed source events (events without predecessor) are strictly decreasing after reaching a circuit. The assumption of realisability is required to avoid that an unrealisable process time just before event i causes a structural delay to zi ðk þ lÞ for all l 2 N, which hides the reduction or settling of the earlier delay zi ðkÞ. If dðkÞ ¼ d0  T k is unrealisable then we can take an eigenvector to define a realisable reference timetable, ~0 ¼ v mod T. Note that the eigenvalue equation is stronger than the realisability condition (12), since d v ¼ apl¼0 Al  kl  v P apl¼0 Al  T l  v , or in conventional notation v i ¼ maxðj;i;lÞ ð½Al ij  l  k þ v j Þ > maxðj;i;lÞ ð½Al ij  l  Tþ v j Þ P ½Al ij þ v j  l  T for all ði; j; l; ½Al ij Þ. h By Lemma 3 any stable scheduled max-plus linear system admits a realisable timetable. Hence, structural delays due to unrealisable process times are unnecessary by defining a suitable timetable (for instance an eigenvector). However, timetables are usually given in minutes and rounding feasible event times to minutes may cause unrealisable process times. Nevertheless, a working timetable with a higher precision of seconds or deci-minutes for internal usage by e.g., railway employee, automatic route setting, and other control systems, can eliminate structural delays. Theorem 2. Let (10) be a stable and realisable scheduled max-plus linear system with acyclic G0 ¼ ðV 0 ; E0 Þ and events topologically ordered with respect to G0 ¼ ðV 0 ; E0 Þ. Then for any initial delay scenario fz1p ; . . . ; z1 g the settling time ks is finite, and Algorithm 1 computes the finite delay list Z, or the delay propagation fzðkÞj1 6 k 6 ks g, in finite time proportional to the number of delays, with the settling time ks depending on fz1p ; . . . ; z1 g. Proof. The theorem follows directly from Lemmas 1, 3 and Theorem 1: by realisability delays are nonincreasing over intervals of p periods (1). Moreover, by stability all delays are decreasing with positive time steps each time they have propagated over a circuit using the positive circuit slack time Lemma 3. Since delays are nonnegative they necessarily settle in finite time, and once settled no new delays are generated by definition of realisability. It follows that the initial delays settle in finite time ks and therefore the total number of delays jZj is finite. The algorithm proceeds by scanning each delayed event exactly once as an immediate consequence of Theorem 1 and its Corollaries 1 and 2, and therefore terminates after a finite number of iterations. h The amount of work done by Algorithm 1 can be compared to that of a straightforward application of the max-plus recursion as given in Goverde (2007, Th. 7). Consider the delay propagation model (10) and assume a finite settling period ks . The direct recursive method using matrix multiplication and addition requires Oðpn3 þ pks n2 Þ operations but also Oðn2 Þ space to store the matrices, which is thus only applicable to small networks. However, a sparse matrix–vector multiplication can be done in OðmÞ time using the arc list representation of sparse matrices. This procedure simply corresponds to scanning all events per period in topological order until all delays are zero for p successive periods. A sparse version of the direct recursive method then requires Oððks þ pÞ  ðm þ nÞÞ operations for ks þ p iterations of a sparse matrix multiplication and a vector addition. Moreover, an additional Oðm0 þ n0 Þ operations are required to sort the arcs in topological order with respect to GðA0 Þ, with n0 and m0 the number of (connected) nodes and arcs of GðA0 Þ, respectively. With the order p a small integral constant (usually, p ¼ 1) and ks depending on the initial conditions, the computational complexity of the sparse direct recursive method is Oðks m þ ks nÞ. Algorithm 1 requires OðjZjÞ operations, with jZj the number of delays, and also Oðm0 þ n0 Þ operations to

283

R.M.P. Goverde / Transportation Research Part C 18 (2010) 269–287 Table 1 Delay propagation behaviour depending on timetable properties.

Stable Critical Unstable a

Realisable

Unrealisablea

Settling Settling or periodic –

Settling Settling or periodic Settling, periodic, or exploding

Structural delays of unreached events discarded.

sort the arc list in topological order with respect to GðA0 Þ. Typically, jZj  ks ðm þ nÞ. Note that jZj ¼ ks n if all events are delayed for all ks periods and suddenly have settled in period ks þ 1. Even in this improbable case, the direct method requires pðm þ nÞ more operations to be sure that all delays have settled. A well designed railway timetable is both realisable and stable, by which initial delays are guaranteed to settle in finite time. Nevertheless, in a preliminary design or in special circumstances these conditions may not be satisfied. Therefore, it is also interesting to study how the algorithm –and delay propagation in general– behaves if one or both of the conditions are not satisfied. Table 1 gives an overview of the five different possible combinations. Note that a realisable timetable cannot be unstable. Theorem 2 above dealt with the case that both stability and realisability are satisfied. The other possibilities are discussed in the next sections. 5.2. Structural delays In this section we assume that the timetable is stable but not realisable. Recall from Section 3.2 that if a timetable is unrealisable then there are process times ½Al ij exceeding the scheduled process time, 0

0

½Al ij > di  dj þ l  T;

ð19Þ

and thus causing structural delays. Nevertheless, if the timetable is stable then downstream of the unreliable processes there is enough slack available for compensation so that on any circuit there is some recovery time guaranteeing stability. This situation typically occurs due to a timetabling trick where the scheduled time–distance trajectory between two stations is ‘bent’ so that it fits the timetable: on paper the timetable then looks feasible but in practice the train will always have some arrival or passing delay which is however compensated for ‘further down the line’. In fact, these timetabling practices make the timetable stability problem nontrivial. The structural primary delays ^z0 are computed using (14). The secondary delays caused by these structural primary delays can be computed using the delay propagation Algorithm 1 with the structural primary delays ^z0 as initial conditions, i.e., zl ¼ ^z0 for all 1  p 6 l 6 1. By stability the delays settle and ^z ¼ zð1Þ is the structural delay vector containing all structural primary delays in period 1 and the secondary delays propagated by the structural primary delays over all periods 1  p 6 l 6 1. Alternatively, the structural delays can be computed using Algorithm 1 with initial conditions zl ¼ e for all 1  p 6 l 6 1. Clearly the algorithm terminates (for stable timetables) and the output ^z ¼ zð1Þ gives the structural delays occurring repetitively each period due to unrealisable process times. Yet another alternative is running Algorithm 1 with iniks zðkÞ. Recall that the settling time ks is finite due to stability. Note that line 12 in Algotial condition z1 ¼ e and taking ^z ¼ al¼0 rithm 1 generates a primary delay z > 0 for any arc ðj; i; l; ½Al ij Þ satisfying the unrealisability condition (19) even though zj ð0Þ ¼ 0. Algorithm 1 only computes delays generated or propagated by events reached from the initial delays. So structural delays become visible only if the delay propagation reaches the preceding unrealisable process times. Furthermore, the structural delays are no longer repetitive when the delays have reduced to a level so that they settle before returning to the unrealisable process. This algorithm behaviour is desired to study properly the propagation of initial delay scenarios without disguising the settling behaviour by structural delays. Clearly the structural delays can be analysed separately as above. Of course, for delay forecasting the impact of both initial and structural delays must be computed. For this, however, the structural delays only have to be computed once because of the linearity in max-plus algebra. n o  Lemma 4. Let x1 ðkÞ ; . . . ; fxq ðkÞg be q 2 N event time sequences satisfying (6) with (different) initial conditions xi1p ; . . . ; xi1 q

for 1 6 i 6 q, respectively, and let fxðkÞg be the event time sequence satisfying (6) with initial condition xl ¼ ai¼1 xil for 1  p 6 l 6 1. Then q

xðkÞ ¼ a xi ðkÞ for all k 2 N: i¼1

Proof. The theorem is a direct consequence of the linearity of the system equations in max-plus algebra. Let x1 ðÞ; x2 ðÞ 2 Rnmax be two event time sequences satisfying (6) with (possibly) different initial conditions x1l and x2l for 1  p 6 l 6 1, respectively. Define xðkÞ ¼ x1 ðkÞ  x2 ðkÞ. Then

284

R.M.P. Goverde / Transportation Research Part C 18 (2010) 269–287

xðkÞ ¼ x1 ðkÞ  x2 ðkÞ; p

p

¼ aðAl  x1 ðk  lÞ  dðkÞÞ  aðAl  x2 ðk  lÞ  dðkÞÞ; l¼0

l¼0

  ¼ a Al  x1 ðk  lÞ  aðAl  x2 ðk  lÞÞ  dðkÞ; p

p

l¼0

l¼0

p

¼ aðAl  x1 ðk  lÞ  Al  x2 ðk  lÞÞ  dðkÞ; l¼0 p

¼ aðAl  ðx1 ðk  lÞ  x2 ðk  lÞÞÞ  dðkÞ; l¼0 p

¼ a Al  xðk  lÞ  dðkÞ: l¼0

The 3rd equality is due to commutativity and idempotency of addition, the 4th due to commutativity of addition, and the 5th due to distributivity. The general case of any finite sum is obtained by taking the two sequences xðkÞ ¼ ðx1 ðkÞ      xq1 Þ  xq ðkÞ. h From Lemma 4 it follows that the combined delay propagation of different sets of initial delays can be obtained by the finite (max-plus) ‘sum’ of the independent delay vectors. In particular, many scenarios consisting of combinations of initial delays can be decomposed in computing the delay propagation of single initial delays and combining the results. Also initial delays and structural delays can be decomposed this way as a direct consequence of Lemma 4, which gives the following theorem. Theorem 3. Let (10) be a stable but unrealisable scheduled max-plus linear system with acyclic GðA0 Þ. Then the delay propagation of both initial delays and structural delays is obtained by the maximum of the structural delays ^z and the delay propagation fzðkÞjk 2 Ng computed using Algorithm 1 with the initial delays z1p ; . . . ; z1 , i.e.,

zall ðkÞ ¼ zðkÞ  ^z;

k 2 N:

Moreover, the delays caused by the initial delays in excess of the structural delays are given by ðzðkÞ  ^zÞ  0 for k 2 N. Suppose that the initial delays do not reach events that start a process with an unrealisable process time. Then the same reasoning as in the proof of Theorem 2 applies and the delays settle in finite time. On the other hand, if a process with an unrealisable process time is reached by some delay sequence then this activity will always generate a primary delay. However, due to the stability condition any delay generated by such a process is reduced again, and even if a delay sequence reaches the process a second time then this delay is smaller than the first time, until at some point the delay is so small that it settles before it can return again. Hence, also for stable but unrealisable systems delay propagation Algorithm 1 terminates although the unrealisable processes would generate new structural delays periodically, which are thus discarded but can be added using Theorem 3. 5.3. Critical systems and periodic delay regimes Now assume that the timetable is realisable but not stable. This implies that the system is critical with k ¼ T, i.e., there are circuits where the processes require exactly a multiple of the cycle time by which no slack is available for delay recovery. If an initial delay starts on or reaches such a critical circuit then delays will circulate over the events in the circuit forever and also keep propagating delays from the critical events to adjacent processes out of the circuit. Hence, the system will converge to a periodic delay regime where all initial delays in the stable parts of the network have settled but the maximum delays that have reached the critical circuits will never settle, and periodically keep spreading delays into the network. This periodic behaviour can be analysed separately by starting the delay propagation algorithm with a delay at some critical event. However, the delay impact area of this critical circuit now depends on the size of the initial delay of the critical event, since the larger this delay the further it will propagate into the network. Clearly, the delay algorithm will not terminate automatically in this case. However, by computing the maximum eigenvalue k before starting the algorithm an additional stop condition of a maximum settling period (say 24 periods) can be defined if k P T. Note that if the initial delays do not reach a critical circuit then the delays settle and the delay propagation algorithm does terminate within a finite number of periods. Indeed, in this case it is interesting whether initial delays reach the critical circuits or not. If the timetable is critical and unrealisable then the same behaviour is observed as in the realisable case, except that also structural delays may reach the critical circuits. In that case, a structural periodic delay regime will always be present. Moreover, this periodic delay regime may get worse if large enough initial delays reach the critical circuits. Nevertheless, since structural delays can be computed separately by Theorem 3, Algorithm 1 may terminate in finite time if the initial delays do not reach the critical circuit either directly or via unrealisable processes.

R.M.P. Goverde / Transportation Research Part C 18 (2010) 269–287

285

Fig. 6. Delay propagation of IC 1900 The Hague-Venlo in railway timetable 2007 (PETER). (a) Propagation over 1st hour (after 0:21:00). (b) Propagation over 2nd hour. (c) Propagation over 3rd hour. (d) Propagation over 4th hour (last delay 3:49:36).

5.4. Unstable systems and delay explosion The final case is an unstable timetable. Note that a realisable timetable can never be unstable since by definition of realisability all process times fit their scheduled process times. Hence, an unstable scheduled railway system implies an unrealisable timetable, with unrealisable process times contained in critical (or unstable) circuits with cycle mean k > T. A delay on an unstable circuit increases after each circulation over the unstable circuit. Hence, initial delays reaching an unstable circuit will explode and the ever increasing delayed events on these circuits will seed a delay propagation that finally will reach all events reachable from the unstable circuits (unless dispatching measures intervene). Obviously, also in this case the delay propagation algorithm must be guarded against non-convergence by setting a maximum number of periods after which the algorithm stops. Initial delays may still settle if they do not reach unstable circuits. More disastrously, unstable circuits necessarily contain unrealisable process times. Hence, the unstable circuits generate structural delays which grow after each circulation of the unstable circuit and thus reach further and further in the network with larger and larger structural secondary delays. In the limit all events in the subgraphs reachable from an unstable circuit have a structural growing delay. Clearly, the structural delays in unstable scheduled systems really hide the delay propagation behaviour of initial delays. Since it is now more important than ever to investigate whether initial delays reach unstable areas or not, decomposition of the structural delays from the initial delays is definitely necessary. The unstable components can also be found using max-plus spectral analysis, see Goverde (2005) and Goverde (2007). Table 1 summarises the possible delay behaviour depending on the timetable properties as discussed in this and the previous sections. If a timetable is unrealisable delays may increase, if it is critical periodic delay regimes may manifest, and for unstable timetables delays may explode. The classification in Table 1 includes the structural delays in the last column (unrealisable timetable), but the classification for initial delays only is based on whether or not the initial delays reach critical or

286

R.M.P. Goverde / Transportation Research Part C 18 (2010) 269–287

unstable circuit(s). Unrealisable or unstable areas may also appear by for example temporary speed reductions or maintenance work. In this case, the delay propagation algorithm can help finding optimal dispatching decisions such as cancelling connections or even trains (or stops) to prevent that delays reach the unstable areas (Braker, 1993; Goverde, 1998; Heidergott and De Vries, 2001). 5.5. Implementation The implemented delay propagation algorithm in the software PETER (Goverde and Odijk, 2002; Goverde, 2005) typically takes less than a second to compute the propagation of initial delays over the dense Dutch national railway timetable. If not yet determined, the algorithm first calls the policy iteration algorithm (Cochet-Terrasson et al., 1998) to compute the maximum eigenvalue k and determine whether the timetable is stable, and if not ðk P TÞ the delay propagation is terminated after 24 periods (default value for T ¼ 60 min). The policy algorithm also terminates within seconds, so that the overall algorithm typically takes a few seconds on a modern PC (Intel CPU 2.80 GHz, 512 MB RAM). The maximal computation time is of course obtained for unstable timetables and many large initial delays reaching the entire network. Note that it is not necessary to compute the eigenvalue in advance: the stop criterion of a default maximum number of periods can be implemented anyhow. A stability analysis can follow in case the algorithm terminates with a warning message that the maximum number of periods had been reached. In PETER, the eigenvalue usually has been computed already before the delay propagation algorithm is applied. This value is then input to the delay propagation algorithm and when it is not yet known it is computed just for additional information. When many initial delay scenarios are evaluated the eigenvalue is of course only determined once after which it is known to the system. Fig. 6a–d shows the PETER visualisation of the propagation of 10 min departure delay of the intercity 1900 from The Hague Central to Venlo in the basic Dutch railway timetable 2007. The visualisations show a projection of the largest delay in a timetable point per period. The colours show the size of the delays according to the legend on the right. The initial delay of 10 min is orange (darkest in B/W), white is no delay. Note that in successive periods the delays reach further over the network but with decreasing delay intensity. There are 257 delayed events (11 arrivals, 142 departures, and 103 passings), the maximum delay is 10:28, and the settling time is ks ¼ 4 periods (hours), although in the fourth period all delays are below 3 min which is registered as punctual in the Netherlands. There are 146 secondary delayed trains (4 transfer connections, 1 turn-around, and 141 headway and route conflicts), although a train can get multiple secondary delays each time it is reached by a larger delay. The cumulative and average secondary delay are 519:40 and 3:33 min, respectively. The other 110 delayed events are follow-up delays at timetable points of already delayed trains. The timetable is unstable due to several single-track lines with unrealisable running times between meeting stations (in this model). There is an unstable circuit in the northeast with eigenvalue 61:02 which reaches only a local area, and an unstable circuit in the east with eigenvalue 60:23 which has access to the main part of the network. The delay from The Hague does however not reach the unstable circuit and thus settles. The number of nodes and arcs in the model are n ¼ 11 401 and m ¼ 89 643, respectively, and the order is p ¼ 1. So a direct computation of the delay propagation using the max-plus recursion would require a constant times ðks þ pÞðm þ nÞ ¼ 4  ð89 643 þ 11 401Þ ¼ 404 176 operations, while Algorithm 1 requires only a constant times jZj ¼ 257 operations. The topological sorting of the arc list is required by both and is not counted here. The direct matrix implementation requires a constant times 1:5  1012 operations. 6. Conclusion This paper introduced an effective delay propagation algorithm based on a timed event graph representation of a scheduled railway system. The model takes into account running time supplements and dwell buffer times to recover from delays, and buffer times to reduce delay propagation to connecting or conflicting trains. The behaviour of the delay propagation has been analysed depending on realisability and stability properties of a timetable, which also led to convergence properties of the algorithm. Using an effective bucket data structure to store the delays during the algorithm, the amount of work per delayed event is minimised to a constant times the number of its successors. The computation time of the algorithm is therefore proportional to the number of delayed events, which depends on the initial delays. The delay propagation algorithm presented in this paper assumes deterministic process times (state matrices) and a perip odic timetable. The algorithm is easily adapted to period-varying systems xðkÞ ¼ al¼0 Al ðkÞxðk  lÞ  dðkÞ, where the state matrices Al depend on the period k. In particular, some train services may not run in off-peak periods or some connections may be cancelled. Also aperiodic timetables pose no problem as long as the scheduled max-plus linear system is realisable p for each k, i.e., dðkÞ P al¼0 Al ðkÞ  dðk  lÞ. Of course, the size of the input data increases if it depends on the periods. The efficiency of the algorithm enables its utilisation in real-time applications such as interactive timetable stability analysis and decision support systems to assist train dispatchers. The fast computational performance is crucial for repeated application of the algorithm in multiple timetable comparisons or dispatching variants to determine the best alternative given a set of actual delays. Moreover, the fast performance allows an embedding the algorithm in simulation optimization methods or model-predictive control frameworks to automatically compute optimal dispatching decisions. This is a topic of current research.

R.M.P. Goverde / Transportation Research Part C 18 (2010) 269–287

287

Acknowledgments This research was partially funded by the Dutch Technology Foundation STW (DWV.6072). The author thanks the three anonymous referees for their helpful comments. References Aho, A.V., Hopcroft, J.E., Ullman, J.D., 1983. Data Structures and Algorithms. Addison-Wesley, Reading. Baccelli, F.L., Cohen, G., Olsder, G.J., Quadrat, J.P., 1992. Synchronization and Linearity: An Algebra for Discrete Event Systems. Wiley, Chichester. Braker, J.G., 1993. Algorithms and Applications in Timed Discrete Event Systems. PhD Thesis, Delft University of Technology, Delft. Cochet-Terrasson, J., Cohen, G., Gaubert, S., Mc Gettrick, M., Quadrat, J.-P., 1998. Numerical computation of spectral elements in max-plus algebra. In: Lafay, J.-F. (Ed.), Proceedings of the IFAC Conference on System Structure and Control 1998, Nantes, France, pp. 699–706 . Cohen, G., Dubois, D., Quadrat, J.-P., Viot, M., 1985. A linear-system-theoretic view of discrete-event processes and its use for performance evaluation in manufacturing. IEEE Transactions on Automatic Control AC 30 (3), 210–220. Cuninghame-Green, R.A., 1979. Minimax Algebra. Lecture Notes in Economics and Mathematical Systems, 166. Springer, Berlin. Dasdan, A., 2004. Experimental analysis of the fastest optimum cycle ratio and mean algorithms. ACM Transactions on Design Automation of Electronic Systems 9 (4), 385–418. De Schutter, B., Van den Boom, T., Hegyi, A., 2002. A model predictive control approach for recovery from delays in railway systems. Transportation Research Record 1793, 15–20. Goverde, R.M.P., 1998. Synchronization control of scheduled train services to minimise passenger waiting times. In: Proceedings of the 4th TRAIL Year Congress, Part 2. TRAIL Research School, Delft. Goverde, R.M.P., 2005. Punctuality of Railway Operations and Timetable Stability Analysis. PhD thesis, Delft University of Technology, Delft. %3chttp:// repository.tudelft.nl/view/ir/uuid:a40ae4f1-1732-4bf3-bbf5-fdb8dfd635e7%3e. Goverde, R.M.P., 2007. Railway timetable stability analysis using max-plus system theory. Transportation Research Part B 41 (2), 179–201. Goverde, R.M.P., Odijk, M.A., 2002. Performance evaluation of network timetables using PETER. In: Allan, J., Hill, R.J., Brebbia, C.A., Sciutto, G., Sone, S. (Eds.), Computers in Railways, vol. VIII. WIT Press, Southampton, pp. 731–740. Heidergott, B., De Vries, R., 2001. Towards a control theory for transportation networks. Discrete Event Dynamic Systems 11 (4), 371–398. Heidergott, B., Olsder, G.J., Van der Woude, J., 2006. Max Plus at Work: Modelling and Analysis of Synchronized Systems. Princeton University Press, Princeton, New Jersey. Knuth, D.E., 1998. The Art of Computer Programming Volume 3: Sorting and Searching, 2nd ed. Addison-Wesley, Boston. Murata, T., 1989. Petri nets: properties, analysis and applications. Proceedings of the IEEE 77 (4), 541–580.