- Email: [email protected]
Generation of classes of robust periodic railway timetables
Computers & Operations Research 33 (2006) 2283 – 2299 www.elsevier.com/locate/cor
Generation of classes of robust periodic railway timetables Michiel A. Odijka,1 , H. Edwin Romeijnb,∗ , Hans van Maarenc a ORTEC Consultants B.V., P.O. Box 490, 2800 AL Gouda, The Netherlands b Department of Industrial and Systems Engineering, University of Florida, 302-B Weil Hall, P.O. Box 116595,
Gainesville, FL 32611-6595, USA c Faculty of Information Technology and Systems, Delft University of Technology, Delft, The Netherlands
Available online 17 March 2005
Abstract In this paper we discuss the problem of randomly sampling classes of fixed-interval railway timetables from a so-called timetable structure. Using a standard model for the timetable structure, we introduce a natural partitioning of the set of feasible timetables into classes. We then define a new probability distribution where the probability of each class depends on the robustness of the timetables in that class. Due to the difficulty of sampling directly from this distribution, we propose a heuristic sampling method and illustrate using practical data that our method indeed favors classes containing robust timetables over others. 䉷 2005 Elsevier Ltd. All rights reserved.
1. Introduction Under the name Rail 21, large parts of the Dutch passenger railway network are (or will be) under (re)construction. The goal of this endeavor is to increase the capacity of the network in anticipation of an increase in the number of passengers in the near future. One of the biggest challenges within this program is the design of expensive railway infrastructure in such a way that it indeed increases the network capacity where such an increase is needed in the future. The railway infrastructure has two natural components:
∗ Corresponding author. Tel.: +1 352 392 3088; fax: +1 352 392 2537.
E-mail addresses: [email protected] (M.A. Odijk), [email protected]fl.edu (H.E. Romeijn), [email protected] (Hans van Maaren). 1 The work of this author was supported by Railned, Utrecht, The Netherlands. 0305-0548/$ - see front matter 䉷 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.cor.2005.02.004
2284
M.A. Odijk et al. / Computers & Operations Research 33 (2006) 2283 – 2299
infrastructure within stations, i.e., tracks and platforms, and infrastructure outside stations, i.e., tracks that connect stations. Decisions about such infrastructure are often made at a higher planning level than decisions about individual station infrastructure and are therefore outside the scope of this paper. In the remainder, we assume that the infrastructure outside the stations is known and fixed and focus is on infrastructure planning within a station. To deal with the uncertainties associated with long-term network capacity planning, the Dutch nonprofit organization for railway capacity planning and allocation called Railned performs extensive scenario studies. The outcomes of these studies are then used to specify railway infrastructure (see, e.g., [1,2]). Scenarios are generally based on origin–destination matrices, i.e., numbers of travellers between railway stations, that are derived from forecast data of demographical, economical, political, and/or geographical nature. Thus, when designing infrastructure, several candidate origin–destination matrices are often taken as starting point. However, designing infrastructure directly from such matrices is generally not possible. Therefore, from each origin–destination matrix less abstract objects should be extracted that can be used for the purpose of specifying railway infrastructure. On the basis of lengthy experience within Railned it was found out that only timetables provide a proper and necessary understanding of the interactions in a dense and complex railway network. Timetables specify arrival and departure times of trains and can be used to specify the required infrastructure in terms of numbers of platforms, tracks, switches, etc. As a first step towards generating timetable scenarios, origin–destination matrices are used to deduce specific sets of formal conditions on the aspired timetables, which specify a timetable structure. More formally, in the context of this paper where we deal with a given railway station, a timetable structure is a system of mathematical constraints on the arrival and departure times of trains at that station. These constraints can model various timetable requirements, e.g., headway times, connections between trains, etc. Any timetable that satisfies these constraints then possesses the desired properties. In this paper, we do not deal with the problem of how to obtain a timetable structure from an origin–destination matrix, nor do we deal with the mathematically complex problem of assuring that the timetable structure is consistent and feasible. For these issues we refer to Hooghiemstra [3] and Odijk [4–6]. Finding a feasible timetable from a timetable structure induces (i) a combinatorial sequencing problem, and (ii) a timing problem. Each solution to the sequencing problem yields a timing problem and each solution to a timing problem yields a timetable. (For details on how to solve the timing problem see, e.g., [5].) Since each timetable belongs to precisely one sequence, the solutions to the sequencing problem partition the set of timetables into timetable classes in a natural way. If a timetable class contains many different timetables, we say that the timetables in this class are robust. For robust timetable classes it is expected that a timetable that is globally feasible (i.e., not only locally at the current station but also globally with respect to the entire network can be found). In addition, robust timetable classes have the property that slight perturbations to the input data can be dealt with by modifying the timetable within its class. (See also [7,8].) We will introduce a probability distribution that assigns higher probability to classes containing robust timetable classes. This distribution is motivated by the fact that timetable planners strongly prefer robust timetables. The problem of randomly generating classes from the proposed distribution is very hard. In particular, the problem of deciding whether a single class exists is already NP-complete (see [9,10]). Therefore, we will provide a heuristic method for sampling approximately from this distribution. The remainder of this paper is organized as follows. In Section 2, we start by discussing the mathematical modeling of timetable structures. In addition, we define the concept of timetable classes, and derive several relevant properties. Finally, we propose a volume and robustness distribution that guide the random
M.A. Odijk et al. / Computers & Operations Research 33 (2006) 2283 – 2299
2285
generation of timetable classes. In Section 3 we develop a heuristic method for sampling according to the robustness distribution, while the performance of our approach is illustrated in Section 4 on a real-life example. We briefly summarize the contributions of the paper in Section 5.
2. Characterization of robust classes of timetables 2.1. The timetable structure model In The Netherlands, as in many other countries, passenger railway timetables are periodic and regular, i.e., the timetable repeats itself regularly in a cyclic fashion, and if a line is serviced more than once per cycle trains call at (approximately) equidistant moments. In this paper, we will use a model first introduced by Serafini and Ukovich [10] as the Periodic Event Scheduling Problem in the context of job shop scheduling. This model was applied to the railway timetabling problem by Odijk [5]. Examples of other applications of this model that have been reported in the literature are airline scheduling (see [11]) and traffic light scheduling (Serafini and Ukovich [12]). This model is currently used by Railned and is very suitable to model timetable periodicity and regularity (see also [13,3]). Formally, a railway timetable specifies an arrival and departure time for all trains that call on the station. These times are often referred to as events, and are denoted by = (0 , . . . , N ). To model driving times, headway times, connections between trains, etc., we can formulate so-called interval constraints on the event times. In particular, such constraints bound the difference between two event times in a cyclic fashion. Denoting the cycle length (which is assumed to be integral) of a timetable by T , this means that an interval constraint bounds the difference between two event times modulo T . Note that a typical choice for T is 60 min, corresponding to an hourly timetable pattern. A typical interval constraint can be written as ij i − j + Tp ij uij ,
(1)
where ij and uij are given integers satisfying 0 < uij − ij < T . The event times i and j as well as the integral variable pij are the unknowns. Some examples of how such interval constraints can be used (within a cyclic schedule with cycle length T = 60) are: • Connection/transfer constraints: Train i should leave the railway station between 2 and 5 min after train j arrives; if i denotes the departure time of train i and j denotes the arrival time of train j , this constraint can be represented as 2 i − j + 60pij 5. • Infrastructure constraints: Trains i and j cannot leave within two minutes of one another; if i and j now denote the departure times of trains i and j , respectively, this constraint can be represented as 2 i − j + 60pij 58. • Frequency constraints: Trains i and j belong to the same family, i.e., they have the same origin and destination, and it is required that these two trains visit the station in a semi-regular fashion; if i and j now denote the arrival (or departure) times of trains i and j , respectively, this constraint can be
2286
M.A. Odijk et al. / Computers & Operations Research 33 (2006) 2283 – 2299
represented as 28 i − j + 60pij 32. We can represent a set of M interval constraints via a simple directed graph G whose vertices correspond to the events. An arc from node i to node j then represents the fact that an interval constraint (1) exists. Denoting the vertex-arc incidence matrix of G by the (N + 1) × M matrix B, the set of timetables can be written in vector notation as Y = (, p) ∈ RN +1 × ZM | B + Tp u , where and u are the vectors of lower and upper bounds, and p is the integral vector p expressing the periodicity of the timetables. This vector is often called the phase shift. Finding a timetable in Y for a given phase shift p (or determining that none exists) is relatively easy (see [14], where an algorithm with an O(N 3 ) running time is described). However, finding a full timetable (, p) in Y is NP-complete in the strong sense (see [9,10]). Serafini and Ukovich [10] developed a branch-and-bound method for solving this problem, whereas Odijk [5] developed a cutting plane method. A particularly powerful solver was developed by Voorhoeve [15] and is currently used by Railned to construct timetables. 2.2. Partitioning the set of timetables into classes It is easy to see that the set of timetables in Y can be written as Y= (, p)| ∈ RN +1 , B + Tp u . p∈ZM
We will next discuss several sources of redundancy in this set which will lead to a reduction in the size of the problem. First, note that adding a constant to all event times in a cyclic timetable does not change its characteristics but simply rotates it in time. Thus, without loss of generality we may arbitrarily fix one event time (say 0 = 0) and redefine the system characterizing the timetable structure accordingly. In addition, it can be shown that, without loss of generality, we can fix (to zero) the elements of the phase shift vector p corresponding to some arbitrary (but fixed) spanning tree T in the graph G (see [10]). Without loss of generality we will assume that we have numbered the interval constraints in such a way that T corresponds to the last N arcs in G (i.e., the last N columns in B, numbered M − N + 1, . . . , M). Reflecting the elimination of both redundancies in the original system, we can now redefine the set Y as follows: (, p)| ∈ RN , B˜ + T Ep u , Y= p∈ZM−N
where B˜ denote the matrix B from which the first row (corresponding to vertex 0) is deleted and E is a matrix that appends N zeroes to the end of a vector of length M − N . Note that, with a slight abuse of notation, we denote the reduced vectors of event times and phase shifts by and p, respectively.
M.A. Odijk et al. / Computers & Operations Research 33 (2006) 2283 – 2299
2287
Finally, observe that 0 < u − < T , which implies that a feasible value of uniquely specifies the corresponding phase shift p. In other words, we can define the set of timetables purely in terms of the event times . A timetable class is then defined to be the set of timetables that are consistent with p, i.e., K(p) = ∈ RN | B˜ + T Ep u . If we denote the set of feasible phase shifts, i.e., the set of vectors p with the property that there exists some such that (, p) ∈ Y , by , the set of all feasible timetables can be written as K(p). p∈
Therefore, the set partitions the timetable structure into a collection of disjoint timetable classes. In the remainder of this paper, we will also consider a relaxation of the phase shift vector to real values. To this end, we will define the set of relaxed timetables as Y¯ = (, x) ∈ RN × RM−N | B˜ + T Ex u = (, x)| ∈ RN , B˜ + T Ex u . x∈RM−N
Furthermore, we denote the set of timetables consistent with the relaxed phase shift x by N ˜ ¯ K(x) = ∈ R | B + T Ex u and the set of feasible relaxed phase shifts by ¯ . The set of all feasible relaxed timetables can then be written as ¯ K(x). ¯ x∈
Note that the set ¯ partitions the relaxed timetable structure into a collection of disjoint timetable classes. 2.3. Properties In this section, we will study properties of the set of (relaxed) timetables as well as the set of feasible (relaxed) phase shifts. Proposition 2.1. The set ¯ is a full-dimensional polytope in RM−N . Moreover, this polytope is centrally symmetric around the point 1 F ( + u), 2T where F = I − (B˜ 2−1 B˜ 1 ) and (B˜ 1 B˜ 2 ) partitions B˜ into its first M − N columns and the remaining N columns.
2288
M.A. Odijk et al. / Computers & Operations Research 33 (2006) 2283 – 2299
Proof. We will first show that the set of relaxed timetables Y¯ is full-dimensional and bounded. First, rewrite the inequalities defining the set Y¯ B˜ + T Ex u as
B˜ 1 B˜ 2
T I M−N 0
x
u,
(2)
by noting that the last N columns of B˜ correspond to the N arcs of the spanning tree T. Since B˜ 2 is invertible, the entire constraint matrix is invertible. The boundedness of Y then follows from the boundedness of and u, and the full-dimensionality follows from the fact that < u. We now immediately conclude that ¯ is a full-dimensional polytope in RM−N since it is the projection of Y¯ onto the M − N coordinates corresponding to x. Finally, the central symmetry of ¯ follows from the fact that 1 I − (B˜ 2−1 B˜ 1 ) ( + u) − x ∈ ¯ , T whenever x ∈ ¯ . We will therefore refer to the set ¯ as the phase shift polytope. Moreover, since ⊆ ¯ , we obtain that there are only a finite number of timetable classes: Corollary 2.2. The set is finite. If the timetable class K(p) corresponding to some phase shift p is of full dimension, then a timetable somewhere in the ‘center’ of K(p) lies away from the boundary of K(p) and is therefore robust compared to timetables close to the boundary. As timetable planners are most interested in robust timetables, it is thus natural to generate phase shifts in a way that is biased towards timetable classes that contain many timetables, i.e., timetable classes with larger volumes. However, it may be the case that no phase shift p exists for which K(p) is full-dimensional, so that the volume of all sets K(p) is equal to zero. In that case, it seems intuitive that we should bias the generation towards phase shifts having the highest dimensionality and, among those, towards ones having the highest volume in that dimension. The following proposition provides insight into the dimensionality of K(x). We will denote the interior of ¯ by int(¯ ) and the boundary of ¯ by j¯ . Proposition 2.3. Let x ∈ ¯ . Then the polytope K(x) is full-dimensional if and only if x ∈ int(¯ ). Moreover, if x ∈ j¯ , then all implicit equalities are of the form i − j = ij − T (Ex)ij i − j = uij − T (Ex)ij .
Proof. Let C be the matrix whose rows are the incidence vectors of all elementary circuits in G and their reverse (where elementary circuits are circuits that are not decomposable into circuits with fewer arcs).
M.A. Odijk et al. / Computers & Operations Research 33 (2006) 2283 – 2299
2289
In addition, define the operator U that, when applied to a circuit c, yields the value 1 + u c − l c− , U (c) = T where c+ = max(c, 0) and c− = max(−c, 0), and when applied to the matrix C (which can be viewed as a column vector of elementary circuits) yields a vector with each element corresponding to an elementary circuit. Following Odijk [5], the set ¯ can be written as ¯ = x ∈ RM−N |CEx U (C) . We will first show that if x ∈ j¯ , then the polytope K(x) contains implicit equalities of the form given in the proposition and is therefore not full-dimensional. By the above characterization of ¯ , there exists an elementary circuit c in G for which c Ex = U (c). Now consider the system defining K(x): − T Ex B˜ u − T Ex.
(3)
Consider the upper bounds in (3) corresponding to forward arcs in c and the lower bounds in (3) corresponding to backward arcs in c, i.e., eij B˜ eij (u − T Ex)
if cij = 1,
−eij B˜ − eij ( − T Ex)
if cij = −1.
(4) (5)
Adding all constraints (4) and (5) yields (c+ − c− ) B˜ (c+ ) (u − T Ex) − (c− ) ( − T Ex) = T U (c) − T c Ex = 0. However, note also that (c+ − c− ) B˜ = c B˜ = 0, since c is a circuit, so that all inequalities in (4) and (5) hold at equality. Conversely, let x ∈ ¯ and suppose that K(x) has an implicit equality, that is, system (3) implies an equality, say y1 (u − T Ex) − y2 (l − T Ex) = 0, y1 B˜ − y2 B˜ = 0, for some (y1 , y2 ) ∈ R2M + not equal to 0. Define ˜ | B(z − z ) = 0 . L = (z1 , z2 ) ∈ R2M 1 2 + Since the extremal rays of L encode the elementary circuits of G and z1 (u − T Ex) − z2 ( − T Ex) 0,
2290
M.A. Odijk et al. / Computers & Operations Research 33 (2006) 2283 – 2299
for all (z1 , z2 ) ∈ L (see [14]), we conclude that there exists an elementary circuit c of G for which (c+ ) (u − T Ex) − (c− ) ( − T Ex) = 0 or, equivalently, c Ex =
1 + (u c − c− ) = U (c). T
Hence, x ∈ j¯ . In the next subsection we will use the above results to propose a continuous probability distribution on the set of relaxed phase shifts which we will call the volume distribution. In addition, we will use this continuous distribution to develop a discrete probability distribution on the set of phase shifts which we will call the robustness distribution. By construction, this probability distribution favors timetable classes containing robust timetables. 2.4. The volume and robustness distributions As mentioned above, a timetable class is robust if the volume of the corresponding set K(p) is large. In order to construct a suitable probability distribution that reflects the preference for such timetable classes, we first define the volume function on ¯ as follows: v(x) = (K(x))
for x ∈ ¯ ,
where the function denotes N-dimensional geometric volume. The following proposition derives some properties of the volume function. Proposition 2.4. The volume function v is a continuous and bounded function on ¯ that is centrally symmetric around the point 1 F ( + u), 2T where F is defined as in Proposition 2.1. Moreover, the maximum of v is attained at this point. Finally, v(x) = 0 if and only if x ∈ j¯ . Proof. We will first show that the function v 1/N is concave. From the definition of the set K(x) it follows that, for all x, x˜ ∈ ¯ and 0 1, K(x) + (1 − )K(x) ˜
= + (1 − )˜ | B˜ + T Ex u, B˜ ˜ + T E x˜ u ⊆ + (1 − )˜ | B˜ ( + (1 − )˜) + T (Ex + (1 − )E x) ˜ u = | B˜ + T (Ex + (1 − )E x) ˜ u ˜ . = K (x + (1 − )x)
M.A. Odijk et al. / Computers & Operations Research 33 (2006) 2283 – 2299
2291
We thus have v 1/N (x + (1 − )x) ˜ = 1/N (K(x + (1 − )x)) ˜ 1/N (K(x) + (1 − )K(x)) ˜ 1/N 1/N (K(x)) + (1 − ) (K(x)) ˜ 1/N 1/N = v (x) + (1 − )v (x) ˜
(6)
where inequality (6) follows from the Brunn–Minkowski Theorem (see, e.g., [16]). Continuity of v now follows from the concavity of v 1/N . Given this result, v is also bounded since its domain ¯ is compact. Since K(x) is a translation of K((1/T )F (+u)−x) the volumes of these two sets are equal. Therefore, it follows that the function v is centrally symmetric around the center of ¯ , the point (1/T )F ( + u). By the fact that v 1/N is concave we conclude that (1/T )F ( + u) in fact maximizes v on ¯ . Finally, it is easy to see that v(x) = 0 if and only if K(x) is not full-dimensional. Therefore, Proposition 2.3 says that v(x) = 0 if and only if x ∈ j¯ . Motivated by these results, we now define the volume distribution V on ¯ as the distribution having density proportional to the volume of the underlying timetable class. More formally, we define the volume distribution as v(y) dy for all measurable sets A ⊆ ¯ . V(A) = A ¯ v(y) dy Note that the results of Propositions 2.1 and 2.4 imply that A v(y) dy is well-defined and that therefore V is a properly defined probability distribution. Using the idea behind the volume distribution defined above, we may now define a discrete probability distribution R on . In particular, we define R as the conditionalization of V onto . If there exists some phase shift p ∈ such that K(p) is of full dimension, i.e., (K(p)) > 0, the robustness distribution is defined as R(p) =
v(p) ˜ p∈ ˜ v(p)
for p ∈ .
(7)
It is clear that the flexibility of timetables in classes with smaller dimension is far smaller than of those with larger dimension. An important property of the robustness distribution that reflects this is the fact that R(p) = 0 if the timetable class K(p) is not full-dimensional. Unfortunately, in general a full-dimensional timetable class K(p) may not exist. To account for these situations, we generalize the definition of R as follows. Let denote the dimension of the timetable class with highest dimension among all timetable classes. Then, let the function v (p) denote the -dimensional geometric volume of the timetable-class K(p) for p ∈ (where, of course, vN = v). Then we extend the definition of R as follows: R(p) =
vn (p) ˜ p∈ ˜ v (p)
for p ∈ .
Clearly, as in the special case of = N discussed above, this distribution has the property that R only assigns positive probability to phase shifts having the largest dimensionality among all phase shifts.
2292
M.A. Odijk et al. / Computers & Operations Research 33 (2006) 2283 – 2299
3. Generating classes of robust timetables The obvious approach towards sampling phase shifts according to robustness distribution is to (i) enumerate all phase shifts, (ii) compute the probability distribution R, and (iii) sample straightforwardly. Unfortunately, in practical applications enumerating all phase shifts can be computationally intractable. Even in cases where enumeration of the phase shifts is possible, the computational burden of the second task of computing the probability distribution R will prevent us from applying this approach. Therefore, we turn towards a heuristic sampling approach that is at least biased towards the generation of classes of robust timetables. In particular, we propose to start by sampling a relaxed phase shift x from the volume distribution. In the next section we will show that it is possible to do this very efficiently. Next, we round x to a phase shift in using an appropriately chosen distance measure. In cases where all phase shifts can be enumerated this rounding step can be done by simply computing the distance of x to each phase shift and finding the nearest one. However, we also propose an optimization approach to this rounding step to deal with cases where enumeration of the set of phase shifts is not practical. 3.1. Sampling from the volume distribution Recall from Section 2.2 the definition of the set of timetables Y . This set can equivalently be written as Y = (, x) ∈ RN × RM−N | x ∈ ¯ , ∈ K(x) . The next lemma shows that sampling a point x from the volume distribution V is equivalent to sampling a point (, x) uniformly from Y and discarding the component . Lemma 3.1. Let the random variable (, x) be uniformly distributed in Y . Then the random variable x is distributed in ¯ according to the volume distribution V. Proof. Let A ⊆ ¯ be measurable and let D ⊆ Y be such that A is its projection on ¯ , i.e. D = (, x) ∈ RM | x ∈ A, ∈ K(x) . Then P r(x ∈ A) = P r((, x) ∈ D) d(, y) = D d(, y) Y A K(y) d dy = d dy ¯ K(y) v(y) dy = A ¯ v(y) dy = V(A).
M.A. Odijk et al. / Computers & Operations Research 33 (2006) 2283 – 2299
2293
We can sample uniformly from Y by means of the inverse transformation technique (see e.g. [17] or [18]) as stated in the following proposition. Proposition 3.2. Let z be distributed uniformly on the hypercube H = {z ∈ RM | z u} and define the random variable x=
1 F z, T
where F is defined as in Proposition 2.1. Then, x is distributed according to the volume distribution V on ¯ . Proof. Following Eq. (2), the set Y can be written as
B˜ T I M−N Y = (z, x) ˜ 1 ∈H , x B2 0 where (B˜ 1 B˜ 2 ) partitions B˜ into its first M − N columns and the remaining N columns (see Proposition 2.1). Since a linear transformation preserves uniformity, the vector B˜ 1 T I M−N −1 z 0 B˜ 2 is uniformly distributed in Y . Now note that, ⎞ ⎛ −1 −1 ˜ 0 B B˜ 1 T I M−N 2 =⎝ 1 ⎠ 1 ˜ B2 0 −1 ˜ ˜ IM−N − B2 B1 T T and recall that F = I − B˜ 2−1 B˜ 1 , which proves the desired result. 3.2. Rounding to a feasible phase shift Given a uniformly generated relaxed phase shift, in the rounding step we wish to find the closest integral phase shift to it. We might simply define ‘closest’ to be in terms of the Euclidean norm on R M−N . However, if we recall the original definition of the timetable classes before eliminating the redundancy in the system, we see that a more appropriate measure of the distance between two (relaxed) phase shifts x and x˜ would be the distance between the two sets − T x B˜ u − T x
(8)
˜ − T x˜ B˜ u − T x.
(9)
and
2294
M.A. Odijk et al. / Computers & Operations Research 33 (2006) 2283 – 2299
This notion of distance removes any undesirable dependence on the choice of spanning tree T used to reduce the dimensionality of the phase shift vector and the choice of event time to reduce the dimension of the timetable. In particular, this means that we can define the distance between x, x˜ ∈ ¯ as (x, x) ¯ = min Ex − E x˜ + B˜ z2 |z ∈ RN . This optimization problem in fact says that the distance between x and x˜ is the Euclidean norm of the projection of the vector E x˜ − Ex = E(x˜ − x) onto the linear space {B˜ z | z ∈ RN }. It follows that an explicit expression for the norm is given by (x, x) ˜ = QE(x − x) ˜ 2,
where ˜ Q = I − B˜ (B˜ B˜ )−1 B. It is easy to verify that is indeed a well defined distance measure on ¯ (see, e.g., [19]). The rounding step of our heuristic can now be formulated as an integer quadratic programming problem: min (x, p) = QE(x − p)2 s.t. p ∈ .
(QP)
This problem is at least as difficult as the so-called Closest Vector problem (see [20]). As mentioned before, if the timetable structure allows for complete enumeration of the phase shifts, (QP) can be solved by inspection. Unfortunately, such approach is not always tractable. We may then solve (QP) either exactly or heuristically. If exact optimization of (QP) is too time consuming, we could solve a slight modification of (QP) that uses the Manhattan norm · 1 rather than the Euclidean norm, which reduces (QP) to a mixed-integer linear programming problem. As a computationally more attractive alternative, we describe a local search heuristic for (QP) that is currently used by Railned. Basically, the algorithm iteratively builds a list p1 , p2 , . . . , of feasible phase shifts. We initialize this list by applying any of the available algorithms that finds a feasible timetable structure p 1 (see, e.g., the ones in [10,15,5]). Then, in iteration n, we try to find a feasible phase shift pn+1 ∈ to add to the list that is closer to x than any of the phase shifts found so far. Here, we may either use the distance or its variant using the Manhattan norm. If such a phase shift cannot be found we stop the algorithm with p = pn . In our approach, we try to find an improving phase shift as follows. We first compute 1 + − u fi − fi , i = 1, . . . , M − N uˆ i = T and ˆi =
1 + fi − u fi− , T
i = 1, . . . , M − N,
where fi is the ith fundamental circuit of G relative to the spanning tree T. Note that ˆ p uˆ for all p ∈ and that ˆ and uˆ have to be calculated only once in the course of the algorithm. Choose a parameter
M.A. Odijk et al. / Computers & Operations Research 33 (2006) 2283 – 2299
2295
0 1 which determines the magnitude of the neighborhood in the local search algorithm. Then, in each iteration, sample vectors p and p¯ as follows:
n pi − 1 with probability , pi = i = 1, . . . , M − N, pin otherwise, and
p¯ i =
pin + 1 pin
with probability , otherwise,
i = 1, . . . , M − N.
Finally, try to find a feasible solution p = p n+1 to the system B + T Ep u, p p p, ¯ ˆ p u, ˆ QE(x − p)
min QE(x − p m ) − ,
m=1,...,n
p ∈ ZM−N , where > 0 is some small threshold improvement. If is not too large, then it can be expected that pi = pin = p¯ i for many i = 1, . . . , M − N, which significantly reduces the size of the problem. If desired, the procedure can be repeated with several values of and . One major advantage of this heuristic is that, in subsequent executions of the sampling heuristic, we may reuse the list of phase shifts previously generated, thereby speeding up the algorithm. 4. An example In this section we illustrate our approach to sampling timetable classes by considering the situation at the Dutch railway station Arnhem CS as an example. Details on this problem instance can be found in [5] and are therefore omitted from this paper. The timetable structure has N = 23 free event times and M = 54 interval constraints. We enumerated all elements in , which yielded 1703 different phase shifts. This enumeration of phase shifts is very time-consuming (on the order of minutes per phase shift). However, this step will not be performed in practice, and is here simply included for testing purposes. Further study of the constraints revealed circuit in the constraint graph that corresponds to a set of interval constraints that have to be satisfied as equalities, which implied that the maximum time table class dimension is N − 3 = 20. In particular, it was found that out of the 1703 phase shifts precisely 1168 phase shifts correspond to a timetable class of dimension N − 3, while the remaining 535 have even smaller dimension. We denote these sets by + and − , respectively. Recall that the robustness distribution R assigns positive probability to only the phase shifts in + . To test our approach, we generated 17,030 samples from V (see Proposition 3.2) and rounded these to integral phase shifts using the heuristic described in Section 3.2. The computation time required for generating a sample from V is negligible (on the order of seconds). However, the rounding phase is more time consuming: to avoid the need for tuning the value of we chose to randomly sample this parameter.
2296
M.A. Odijk et al. / Computers & Operations Research 33 (2006) 2283 – 2299
Fig. 1. Results for the robustness distribution applied to Arnhem CS for all 1703 phase shifts together.
The rounding step took an average of 18 iterations and 138 s. We then performed a frequency count for each of the 1703 phase shifts. It turned out that approximately 13% of the phase shifts in + occurred with frequency zero, while this was the case for approximately 40% of the phase shifts in − . This suggests that the heuristic is indeed biased towards avoiding timetable classes with low dimension. Computing the volumes of the 1703 timetable classes is too time consuming. In order to assess the success of our approach in biasing the random generation towards timetable classes containing robust timetables, we determined the distance of each of the sampled phase shifts to the center of the polytope ¯ , which is also the relaxed phase shift at which the volume function v attains its maximum (see Propositions 2.1 and 2.4). These results are plotted in Fig. 1. In this figure, the size of the dot indicates the number of phase shifts that the dot corresponds to. To improve the readability of the figure, seven phase shifts were omitted since they had a frequency of over 85 (the highest frequency observed was 126), although these observations were of course taken into account in all calculations. The horizontal and vertical lines indicate the mean distance and frequency, respectively. This distance measure is based on the intuitive expectation that phase shifts at a large distance from the center of ¯ are assigned less probability mass by R than phase shifts for which this distance is small. This is motivated by the fact that the function v 1/N is concave and the center of ¯ maximizes v. To test this expectation, we have performed a linear regression of the distance on the observed frequency. This regression yielded a slope of −0.0034 with a t-value of 35, which means that the slope is very significantly negative. Figs. 2 and 3 show the data of Fig. 1 partitioned by the phase shifts in + and − , respectively. In this case, performing linear regression of the distance on observed frequency yielded a slope of −0.0030 with a t-value 31 for the data in Fig. 2 and a slope of −0.0032 with a t-value of 15 for the data in Fig. 3. In both cases the negative slope is quite significant.
M.A. Odijk et al. / Computers & Operations Research 33 (2006) 2283 – 2299
2297
Fig. 2. Results for the robustness distribution applied to Arnhem CS for the phase shifts with a timetable set of dimension N − 3.
Fig. 3. Results for the robustness distribution applied to Arnhem CS for the phase shifts with a timetable set of dimension less than N − 3.
Finally, we may observe from the latter two figures that the phase shifts that should be out of sight, i.e. the phase shifts in − , have a much smaller average frequency and a somewhat larger distance to the center of ¯ than the phase shifts in + .
2298
M.A. Odijk et al. / Computers & Operations Research 33 (2006) 2283 – 2299
5. Summary In this paper we studied the problem of randomly generating timetable classes with a bias towards classes that contain robust timetables. This problem is very relevant in performing capacity planning using a scenario analysis to take into account uncertainties in future demand for rail travel. We formalized the notion of a robust timetable and developed a heuristic approach to solve problem of randomly generating classes of robust timetables. Using a practical example, we have illustrated that our method shows promise as it indeed favors classes of robust timetables over other classes. References [1] Odijk MA, Maarschalkerweerd MSPAM, Weits EAG. Timetable—independent design and analysis of railway infrastructure. In: Allan J, Brebbia CA, Hill RJ, Sciutto G, Sone S, editors. Computers in railways V, vol. 1: railway systems and management. Proceedings of CompRail 96, Berlin, Germany; 1996. pp. 3–12. [2] Odijk MA, Zwaneveld PJ, Hooghiemstra JS, Kroon LG, Salomon M. Decision support systems help railned to search for ‘win-win’ solutions in railway network design. Interfaces 1999;29(2):15–32. [3] Hooghiemstra JS. Design of regular interval timetables for strategic and tactical railway planning. In: Allan J, Brebbia CA, Hill RJ, Sciutto G, Sone S, editors. Computers in railways V, vol. 1: railway systems and management. Proceedings of CompRail 96, Berlin, Germany; 1996. pp. 393–402. [4] Odijk MA. Uniform generation of regular interval timetables to facilitate and improve the design of rail infrastructure. In: Hertel G, editor. Schriftenreihe des Instituts für Verkehrssystemtheorie und Bahnverkehr 2. Proceedings of Das Zweite Eisenbahnbetriebswis- senschaftliches Kolloquium, TU Dresden, Dresden, Germany; 1995. pp. 115–26. [5] Odijk MA. A constraint generation algorithm for the construction of periodic railway timetables. Transportation Research B 1996;30(6):455–64. [6] Odijk MA. Sensitivity analysis of a railway stations track layout with respect to a given timetable. European Journal of Operational Research 1999;112(3):517–30. [7] Peeters L, Kroon L. A cycle based optimization model for the cyclic railway timetabling problem. In: Voss S, Daduna JR, editors. Computer-aided scheduling of public transport. Lecture notes in economics and mathematical systems 505. Berlin: Springer; 2001. pp. 275–96. [8] Liebchen C, Möhring RH. A case study in periodic timetabling. In: Proceedings of ATMOS 2002—algorithmic methods and models for optimization of railways, Malaga, Spain. Electronic Notes in Theoretical Computer Science 2002; 66(6). [9] Garey MR, Johnson DS. Computers and intractability—a guide to the theory of NP-completeness. New York: Freeman and Company; 1979. [10] Serafini P, Ukovich W. A mathematical model for periodic scheduling problems. SIAM Journal on Discrete Mathematics 1989;2(4):550–81. [11] Gertsbakh I, Serafini P. Periodic transportation schedules with flexible departure times. European Journal of Operational Research 1991;50:298–309. [12] Serafini P, Ukovich W. A mathematical model for the fixed-time traffic control problem. European Journal of Operational Research 1989;42:152–65. [13] van den Berg JHA, Odijk MA. DONS: Computer aided design of regular service time-tables. In: Murthy TKS, Brebbia CA, Mellitt B, Sciutto G, Sone S, editors. Computers in railways IV, vol. 2: railway operations. Proceedings of CompRail 94, Madrid, Spain; 1994. pp. 109–16. [14] Rockafellar RT. Network flows and monotropic optimization. New York: Wiley; 1984. [15] Voorhoeve M. Rail scheduling with discrete sets. Working paper. Eindhoven University of Technology, Eindhoven, The Netherlands. 1993. [16] Bonnesen T, Fenchel W. Theorie der Konvexen Körper. Berlin: Springer; 1935 (in German). [17] Rubinstein RY. Simulation and the Monte Carlo method. New York: Wiley; 1981.
M.A. Odijk et al. / Computers & Operations Research 33 (2006) 2283 – 2299
2299
[18] Smith RL. Efficient Monte Carlo procedures for generating points uniformly distributed over bounded regions. Operations Research 1984;32(6):1296–308. [19] Luenberger DG. Optimization by vector space methods. New York: Wiley; 1969. [20] Nemhauser GL, Wolsey LA. Integer and combinatorial optimization. New York: Wiley; 1988.
Generation of classes of robust periodic railway timetables Michiel A. Odijka,1 , H. Edwin Romeijnb,∗ , Hans van Maarenc a ORTEC Consultants B.V., P.O. Box 490, 2800 AL Gouda, The Netherlands b Department of Industrial and Systems Engineering, University of Florida, 302-B Weil Hall, P.O. Box 116595,
Gainesville, FL 32611-6595, USA c Faculty of Information Technology and Systems, Delft University of Technology, Delft, The Netherlands
Available online 17 March 2005
Abstract In this paper we discuss the problem of randomly sampling classes of fixed-interval railway timetables from a so-called timetable structure. Using a standard model for the timetable structure, we introduce a natural partitioning of the set of feasible timetables into classes. We then define a new probability distribution where the probability of each class depends on the robustness of the timetables in that class. Due to the difficulty of sampling directly from this distribution, we propose a heuristic sampling method and illustrate using practical data that our method indeed favors classes containing robust timetables over others. 䉷 2005 Elsevier Ltd. All rights reserved.
1. Introduction Under the name Rail 21, large parts of the Dutch passenger railway network are (or will be) under (re)construction. The goal of this endeavor is to increase the capacity of the network in anticipation of an increase in the number of passengers in the near future. One of the biggest challenges within this program is the design of expensive railway infrastructure in such a way that it indeed increases the network capacity where such an increase is needed in the future. The railway infrastructure has two natural components:
∗ Corresponding author. Tel.: +1 352 392 3088; fax: +1 352 392 2537.
E-mail addresses: [email protected] (M.A. Odijk), [email protected]fl.edu (H.E. Romeijn), [email protected] (Hans van Maaren). 1 The work of this author was supported by Railned, Utrecht, The Netherlands. 0305-0548/$ - see front matter 䉷 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.cor.2005.02.004
2284
M.A. Odijk et al. / Computers & Operations Research 33 (2006) 2283 – 2299
infrastructure within stations, i.e., tracks and platforms, and infrastructure outside stations, i.e., tracks that connect stations. Decisions about such infrastructure are often made at a higher planning level than decisions about individual station infrastructure and are therefore outside the scope of this paper. In the remainder, we assume that the infrastructure outside the stations is known and fixed and focus is on infrastructure planning within a station. To deal with the uncertainties associated with long-term network capacity planning, the Dutch nonprofit organization for railway capacity planning and allocation called Railned performs extensive scenario studies. The outcomes of these studies are then used to specify railway infrastructure (see, e.g., [1,2]). Scenarios are generally based on origin–destination matrices, i.e., numbers of travellers between railway stations, that are derived from forecast data of demographical, economical, political, and/or geographical nature. Thus, when designing infrastructure, several candidate origin–destination matrices are often taken as starting point. However, designing infrastructure directly from such matrices is generally not possible. Therefore, from each origin–destination matrix less abstract objects should be extracted that can be used for the purpose of specifying railway infrastructure. On the basis of lengthy experience within Railned it was found out that only timetables provide a proper and necessary understanding of the interactions in a dense and complex railway network. Timetables specify arrival and departure times of trains and can be used to specify the required infrastructure in terms of numbers of platforms, tracks, switches, etc. As a first step towards generating timetable scenarios, origin–destination matrices are used to deduce specific sets of formal conditions on the aspired timetables, which specify a timetable structure. More formally, in the context of this paper where we deal with a given railway station, a timetable structure is a system of mathematical constraints on the arrival and departure times of trains at that station. These constraints can model various timetable requirements, e.g., headway times, connections between trains, etc. Any timetable that satisfies these constraints then possesses the desired properties. In this paper, we do not deal with the problem of how to obtain a timetable structure from an origin–destination matrix, nor do we deal with the mathematically complex problem of assuring that the timetable structure is consistent and feasible. For these issues we refer to Hooghiemstra [3] and Odijk [4–6]. Finding a feasible timetable from a timetable structure induces (i) a combinatorial sequencing problem, and (ii) a timing problem. Each solution to the sequencing problem yields a timing problem and each solution to a timing problem yields a timetable. (For details on how to solve the timing problem see, e.g., [5].) Since each timetable belongs to precisely one sequence, the solutions to the sequencing problem partition the set of timetables into timetable classes in a natural way. If a timetable class contains many different timetables, we say that the timetables in this class are robust. For robust timetable classes it is expected that a timetable that is globally feasible (i.e., not only locally at the current station but also globally with respect to the entire network can be found). In addition, robust timetable classes have the property that slight perturbations to the input data can be dealt with by modifying the timetable within its class. (See also [7,8].) We will introduce a probability distribution that assigns higher probability to classes containing robust timetable classes. This distribution is motivated by the fact that timetable planners strongly prefer robust timetables. The problem of randomly generating classes from the proposed distribution is very hard. In particular, the problem of deciding whether a single class exists is already NP-complete (see [9,10]). Therefore, we will provide a heuristic method for sampling approximately from this distribution. The remainder of this paper is organized as follows. In Section 2, we start by discussing the mathematical modeling of timetable structures. In addition, we define the concept of timetable classes, and derive several relevant properties. Finally, we propose a volume and robustness distribution that guide the random
M.A. Odijk et al. / Computers & Operations Research 33 (2006) 2283 – 2299
2285
generation of timetable classes. In Section 3 we develop a heuristic method for sampling according to the robustness distribution, while the performance of our approach is illustrated in Section 4 on a real-life example. We briefly summarize the contributions of the paper in Section 5.
2. Characterization of robust classes of timetables 2.1. The timetable structure model In The Netherlands, as in many other countries, passenger railway timetables are periodic and regular, i.e., the timetable repeats itself regularly in a cyclic fashion, and if a line is serviced more than once per cycle trains call at (approximately) equidistant moments. In this paper, we will use a model first introduced by Serafini and Ukovich [10] as the Periodic Event Scheduling Problem in the context of job shop scheduling. This model was applied to the railway timetabling problem by Odijk [5]. Examples of other applications of this model that have been reported in the literature are airline scheduling (see [11]) and traffic light scheduling (Serafini and Ukovich [12]). This model is currently used by Railned and is very suitable to model timetable periodicity and regularity (see also [13,3]). Formally, a railway timetable specifies an arrival and departure time for all trains that call on the station. These times are often referred to as events, and are denoted by = (0 , . . . , N ). To model driving times, headway times, connections between trains, etc., we can formulate so-called interval constraints on the event times. In particular, such constraints bound the difference between two event times in a cyclic fashion. Denoting the cycle length (which is assumed to be integral) of a timetable by T , this means that an interval constraint bounds the difference between two event times modulo T . Note that a typical choice for T is 60 min, corresponding to an hourly timetable pattern. A typical interval constraint can be written as ij i − j + Tp ij uij ,
(1)
where ij and uij are given integers satisfying 0 < uij − ij < T . The event times i and j as well as the integral variable pij are the unknowns. Some examples of how such interval constraints can be used (within a cyclic schedule with cycle length T = 60) are: • Connection/transfer constraints: Train i should leave the railway station between 2 and 5 min after train j arrives; if i denotes the departure time of train i and j denotes the arrival time of train j , this constraint can be represented as 2 i − j + 60pij 5. • Infrastructure constraints: Trains i and j cannot leave within two minutes of one another; if i and j now denote the departure times of trains i and j , respectively, this constraint can be represented as 2 i − j + 60pij 58. • Frequency constraints: Trains i and j belong to the same family, i.e., they have the same origin and destination, and it is required that these two trains visit the station in a semi-regular fashion; if i and j now denote the arrival (or departure) times of trains i and j , respectively, this constraint can be
2286
M.A. Odijk et al. / Computers & Operations Research 33 (2006) 2283 – 2299
represented as 28 i − j + 60pij 32. We can represent a set of M interval constraints via a simple directed graph G whose vertices correspond to the events. An arc from node i to node j then represents the fact that an interval constraint (1) exists. Denoting the vertex-arc incidence matrix of G by the (N + 1) × M matrix B, the set of timetables can be written in vector notation as Y = (, p) ∈ RN +1 × ZM | B + Tp u , where and u are the vectors of lower and upper bounds, and p is the integral vector p expressing the periodicity of the timetables. This vector is often called the phase shift. Finding a timetable in Y for a given phase shift p (or determining that none exists) is relatively easy (see [14], where an algorithm with an O(N 3 ) running time is described). However, finding a full timetable (, p) in Y is NP-complete in the strong sense (see [9,10]). Serafini and Ukovich [10] developed a branch-and-bound method for solving this problem, whereas Odijk [5] developed a cutting plane method. A particularly powerful solver was developed by Voorhoeve [15] and is currently used by Railned to construct timetables. 2.2. Partitioning the set of timetables into classes It is easy to see that the set of timetables in Y can be written as Y= (, p)| ∈ RN +1 , B + Tp u . p∈ZM
We will next discuss several sources of redundancy in this set which will lead to a reduction in the size of the problem. First, note that adding a constant to all event times in a cyclic timetable does not change its characteristics but simply rotates it in time. Thus, without loss of generality we may arbitrarily fix one event time (say 0 = 0) and redefine the system characterizing the timetable structure accordingly. In addition, it can be shown that, without loss of generality, we can fix (to zero) the elements of the phase shift vector p corresponding to some arbitrary (but fixed) spanning tree T in the graph G (see [10]). Without loss of generality we will assume that we have numbered the interval constraints in such a way that T corresponds to the last N arcs in G (i.e., the last N columns in B, numbered M − N + 1, . . . , M). Reflecting the elimination of both redundancies in the original system, we can now redefine the set Y as follows: (, p)| ∈ RN , B˜ + T Ep u , Y= p∈ZM−N
where B˜ denote the matrix B from which the first row (corresponding to vertex 0) is deleted and E is a matrix that appends N zeroes to the end of a vector of length M − N . Note that, with a slight abuse of notation, we denote the reduced vectors of event times and phase shifts by and p, respectively.
M.A. Odijk et al. / Computers & Operations Research 33 (2006) 2283 – 2299
2287
Finally, observe that 0 < u − < T , which implies that a feasible value of uniquely specifies the corresponding phase shift p. In other words, we can define the set of timetables purely in terms of the event times . A timetable class is then defined to be the set of timetables that are consistent with p, i.e., K(p) = ∈ RN | B˜ + T Ep u . If we denote the set of feasible phase shifts, i.e., the set of vectors p with the property that there exists some such that (, p) ∈ Y , by , the set of all feasible timetables can be written as K(p). p∈
Therefore, the set partitions the timetable structure into a collection of disjoint timetable classes. In the remainder of this paper, we will also consider a relaxation of the phase shift vector to real values. To this end, we will define the set of relaxed timetables as Y¯ = (, x) ∈ RN × RM−N | B˜ + T Ex u = (, x)| ∈ RN , B˜ + T Ex u . x∈RM−N
Furthermore, we denote the set of timetables consistent with the relaxed phase shift x by N ˜ ¯ K(x) = ∈ R | B + T Ex u and the set of feasible relaxed phase shifts by ¯ . The set of all feasible relaxed timetables can then be written as ¯ K(x). ¯ x∈
Note that the set ¯ partitions the relaxed timetable structure into a collection of disjoint timetable classes. 2.3. Properties In this section, we will study properties of the set of (relaxed) timetables as well as the set of feasible (relaxed) phase shifts. Proposition 2.1. The set ¯ is a full-dimensional polytope in RM−N . Moreover, this polytope is centrally symmetric around the point 1 F ( + u), 2T where F = I − (B˜ 2−1 B˜ 1 ) and (B˜ 1 B˜ 2 ) partitions B˜ into its first M − N columns and the remaining N columns.
2288
M.A. Odijk et al. / Computers & Operations Research 33 (2006) 2283 – 2299
Proof. We will first show that the set of relaxed timetables Y¯ is full-dimensional and bounded. First, rewrite the inequalities defining the set Y¯ B˜ + T Ex u as
B˜ 1 B˜ 2
T I M−N 0
x
u,
(2)
by noting that the last N columns of B˜ correspond to the N arcs of the spanning tree T. Since B˜ 2 is invertible, the entire constraint matrix is invertible. The boundedness of Y then follows from the boundedness of and u, and the full-dimensionality follows from the fact that < u. We now immediately conclude that ¯ is a full-dimensional polytope in RM−N since it is the projection of Y¯ onto the M − N coordinates corresponding to x. Finally, the central symmetry of ¯ follows from the fact that 1 I − (B˜ 2−1 B˜ 1 ) ( + u) − x ∈ ¯ , T whenever x ∈ ¯ . We will therefore refer to the set ¯ as the phase shift polytope. Moreover, since ⊆ ¯ , we obtain that there are only a finite number of timetable classes: Corollary 2.2. The set is finite. If the timetable class K(p) corresponding to some phase shift p is of full dimension, then a timetable somewhere in the ‘center’ of K(p) lies away from the boundary of K(p) and is therefore robust compared to timetables close to the boundary. As timetable planners are most interested in robust timetables, it is thus natural to generate phase shifts in a way that is biased towards timetable classes that contain many timetables, i.e., timetable classes with larger volumes. However, it may be the case that no phase shift p exists for which K(p) is full-dimensional, so that the volume of all sets K(p) is equal to zero. In that case, it seems intuitive that we should bias the generation towards phase shifts having the highest dimensionality and, among those, towards ones having the highest volume in that dimension. The following proposition provides insight into the dimensionality of K(x). We will denote the interior of ¯ by int(¯ ) and the boundary of ¯ by j¯ . Proposition 2.3. Let x ∈ ¯ . Then the polytope K(x) is full-dimensional if and only if x ∈ int(¯ ). Moreover, if x ∈ j¯ , then all implicit equalities are of the form i − j = ij − T (Ex)ij i − j = uij − T (Ex)ij .
Proof. Let C be the matrix whose rows are the incidence vectors of all elementary circuits in G and their reverse (where elementary circuits are circuits that are not decomposable into circuits with fewer arcs).
M.A. Odijk et al. / Computers & Operations Research 33 (2006) 2283 – 2299
2289
In addition, define the operator U that, when applied to a circuit c, yields the value 1 + u c − l c− , U (c) = T where c+ = max(c, 0) and c− = max(−c, 0), and when applied to the matrix C (which can be viewed as a column vector of elementary circuits) yields a vector with each element corresponding to an elementary circuit. Following Odijk [5], the set ¯ can be written as ¯ = x ∈ RM−N |CEx U (C) . We will first show that if x ∈ j¯ , then the polytope K(x) contains implicit equalities of the form given in the proposition and is therefore not full-dimensional. By the above characterization of ¯ , there exists an elementary circuit c in G for which c Ex = U (c). Now consider the system defining K(x): − T Ex B˜ u − T Ex.
(3)
Consider the upper bounds in (3) corresponding to forward arcs in c and the lower bounds in (3) corresponding to backward arcs in c, i.e., eij B˜ eij (u − T Ex)
if cij = 1,
−eij B˜ − eij ( − T Ex)
if cij = −1.
(4) (5)
Adding all constraints (4) and (5) yields (c+ − c− ) B˜ (c+ ) (u − T Ex) − (c− ) ( − T Ex) = T U (c) − T c Ex = 0. However, note also that (c+ − c− ) B˜ = c B˜ = 0, since c is a circuit, so that all inequalities in (4) and (5) hold at equality. Conversely, let x ∈ ¯ and suppose that K(x) has an implicit equality, that is, system (3) implies an equality, say y1 (u − T Ex) − y2 (l − T Ex) = 0, y1 B˜ − y2 B˜ = 0, for some (y1 , y2 ) ∈ R2M + not equal to 0. Define ˜ | B(z − z ) = 0 . L = (z1 , z2 ) ∈ R2M 1 2 + Since the extremal rays of L encode the elementary circuits of G and z1 (u − T Ex) − z2 ( − T Ex) 0,
2290
M.A. Odijk et al. / Computers & Operations Research 33 (2006) 2283 – 2299
for all (z1 , z2 ) ∈ L (see [14]), we conclude that there exists an elementary circuit c of G for which (c+ ) (u − T Ex) − (c− ) ( − T Ex) = 0 or, equivalently, c Ex =
1 + (u c − c− ) = U (c). T
Hence, x ∈ j¯ . In the next subsection we will use the above results to propose a continuous probability distribution on the set of relaxed phase shifts which we will call the volume distribution. In addition, we will use this continuous distribution to develop a discrete probability distribution on the set of phase shifts which we will call the robustness distribution. By construction, this probability distribution favors timetable classes containing robust timetables. 2.4. The volume and robustness distributions As mentioned above, a timetable class is robust if the volume of the corresponding set K(p) is large. In order to construct a suitable probability distribution that reflects the preference for such timetable classes, we first define the volume function on ¯ as follows: v(x) = (K(x))
for x ∈ ¯ ,
where the function denotes N-dimensional geometric volume. The following proposition derives some properties of the volume function. Proposition 2.4. The volume function v is a continuous and bounded function on ¯ that is centrally symmetric around the point 1 F ( + u), 2T where F is defined as in Proposition 2.1. Moreover, the maximum of v is attained at this point. Finally, v(x) = 0 if and only if x ∈ j¯ . Proof. We will first show that the function v 1/N is concave. From the definition of the set K(x) it follows that, for all x, x˜ ∈ ¯ and 0 1, K(x) + (1 − )K(x) ˜
= + (1 − )˜ | B˜ + T Ex u, B˜ ˜ + T E x˜ u ⊆ + (1 − )˜ | B˜ ( + (1 − )˜) + T (Ex + (1 − )E x) ˜ u = | B˜ + T (Ex + (1 − )E x) ˜ u ˜ . = K (x + (1 − )x)
M.A. Odijk et al. / Computers & Operations Research 33 (2006) 2283 – 2299
2291
We thus have v 1/N (x + (1 − )x) ˜ = 1/N (K(x + (1 − )x)) ˜ 1/N (K(x) + (1 − )K(x)) ˜ 1/N 1/N (K(x)) + (1 − ) (K(x)) ˜ 1/N 1/N = v (x) + (1 − )v (x) ˜
(6)
where inequality (6) follows from the Brunn–Minkowski Theorem (see, e.g., [16]). Continuity of v now follows from the concavity of v 1/N . Given this result, v is also bounded since its domain ¯ is compact. Since K(x) is a translation of K((1/T )F (+u)−x) the volumes of these two sets are equal. Therefore, it follows that the function v is centrally symmetric around the center of ¯ , the point (1/T )F ( + u). By the fact that v 1/N is concave we conclude that (1/T )F ( + u) in fact maximizes v on ¯ . Finally, it is easy to see that v(x) = 0 if and only if K(x) is not full-dimensional. Therefore, Proposition 2.3 says that v(x) = 0 if and only if x ∈ j¯ . Motivated by these results, we now define the volume distribution V on ¯ as the distribution having density proportional to the volume of the underlying timetable class. More formally, we define the volume distribution as v(y) dy for all measurable sets A ⊆ ¯ . V(A) = A ¯ v(y) dy Note that the results of Propositions 2.1 and 2.4 imply that A v(y) dy is well-defined and that therefore V is a properly defined probability distribution. Using the idea behind the volume distribution defined above, we may now define a discrete probability distribution R on . In particular, we define R as the conditionalization of V onto . If there exists some phase shift p ∈ such that K(p) is of full dimension, i.e., (K(p)) > 0, the robustness distribution is defined as R(p) =
v(p) ˜ p∈ ˜ v(p)
for p ∈ .
(7)
It is clear that the flexibility of timetables in classes with smaller dimension is far smaller than of those with larger dimension. An important property of the robustness distribution that reflects this is the fact that R(p) = 0 if the timetable class K(p) is not full-dimensional. Unfortunately, in general a full-dimensional timetable class K(p) may not exist. To account for these situations, we generalize the definition of R as follows. Let denote the dimension of the timetable class with highest dimension among all timetable classes. Then, let the function v (p) denote the -dimensional geometric volume of the timetable-class K(p) for p ∈ (where, of course, vN = v). Then we extend the definition of R as follows: R(p) =
vn (p) ˜ p∈ ˜ v (p)
for p ∈ .
Clearly, as in the special case of = N discussed above, this distribution has the property that R only assigns positive probability to phase shifts having the largest dimensionality among all phase shifts.
2292
M.A. Odijk et al. / Computers & Operations Research 33 (2006) 2283 – 2299
3. Generating classes of robust timetables The obvious approach towards sampling phase shifts according to robustness distribution is to (i) enumerate all phase shifts, (ii) compute the probability distribution R, and (iii) sample straightforwardly. Unfortunately, in practical applications enumerating all phase shifts can be computationally intractable. Even in cases where enumeration of the phase shifts is possible, the computational burden of the second task of computing the probability distribution R will prevent us from applying this approach. Therefore, we turn towards a heuristic sampling approach that is at least biased towards the generation of classes of robust timetables. In particular, we propose to start by sampling a relaxed phase shift x from the volume distribution. In the next section we will show that it is possible to do this very efficiently. Next, we round x to a phase shift in using an appropriately chosen distance measure. In cases where all phase shifts can be enumerated this rounding step can be done by simply computing the distance of x to each phase shift and finding the nearest one. However, we also propose an optimization approach to this rounding step to deal with cases where enumeration of the set of phase shifts is not practical. 3.1. Sampling from the volume distribution Recall from Section 2.2 the definition of the set of timetables Y . This set can equivalently be written as Y = (, x) ∈ RN × RM−N | x ∈ ¯ , ∈ K(x) . The next lemma shows that sampling a point x from the volume distribution V is equivalent to sampling a point (, x) uniformly from Y and discarding the component . Lemma 3.1. Let the random variable (, x) be uniformly distributed in Y . Then the random variable x is distributed in ¯ according to the volume distribution V. Proof. Let A ⊆ ¯ be measurable and let D ⊆ Y be such that A is its projection on ¯ , i.e. D = (, x) ∈ RM | x ∈ A, ∈ K(x) . Then P r(x ∈ A) = P r((, x) ∈ D) d(, y) = D d(, y) Y A K(y) d dy = d dy ¯ K(y) v(y) dy = A ¯ v(y) dy = V(A).
M.A. Odijk et al. / Computers & Operations Research 33 (2006) 2283 – 2299
2293
We can sample uniformly from Y by means of the inverse transformation technique (see e.g. [17] or [18]) as stated in the following proposition. Proposition 3.2. Let z be distributed uniformly on the hypercube H = {z ∈ RM | z u} and define the random variable x=
1 F z, T
where F is defined as in Proposition 2.1. Then, x is distributed according to the volume distribution V on ¯ . Proof. Following Eq. (2), the set Y can be written as
B˜ T I M−N Y = (z, x) ˜ 1 ∈H , x B2 0 where (B˜ 1 B˜ 2 ) partitions B˜ into its first M − N columns and the remaining N columns (see Proposition 2.1). Since a linear transformation preserves uniformity, the vector B˜ 1 T I M−N −1 z 0 B˜ 2 is uniformly distributed in Y . Now note that, ⎞ ⎛ −1 −1 ˜ 0 B B˜ 1 T I M−N 2 =⎝ 1 ⎠ 1 ˜ B2 0 −1 ˜ ˜ IM−N − B2 B1 T T and recall that F = I − B˜ 2−1 B˜ 1 , which proves the desired result. 3.2. Rounding to a feasible phase shift Given a uniformly generated relaxed phase shift, in the rounding step we wish to find the closest integral phase shift to it. We might simply define ‘closest’ to be in terms of the Euclidean norm on R M−N . However, if we recall the original definition of the timetable classes before eliminating the redundancy in the system, we see that a more appropriate measure of the distance between two (relaxed) phase shifts x and x˜ would be the distance between the two sets − T x B˜ u − T x
(8)
˜ − T x˜ B˜ u − T x.
(9)
and
2294
M.A. Odijk et al. / Computers & Operations Research 33 (2006) 2283 – 2299
This notion of distance removes any undesirable dependence on the choice of spanning tree T used to reduce the dimensionality of the phase shift vector and the choice of event time to reduce the dimension of the timetable. In particular, this means that we can define the distance between x, x˜ ∈ ¯ as (x, x) ¯ = min Ex − E x˜ + B˜ z2 |z ∈ RN . This optimization problem in fact says that the distance between x and x˜ is the Euclidean norm of the projection of the vector E x˜ − Ex = E(x˜ − x) onto the linear space {B˜ z | z ∈ RN }. It follows that an explicit expression for the norm is given by (x, x) ˜ = QE(x − x) ˜ 2,
where ˜ Q = I − B˜ (B˜ B˜ )−1 B. It is easy to verify that is indeed a well defined distance measure on ¯ (see, e.g., [19]). The rounding step of our heuristic can now be formulated as an integer quadratic programming problem: min (x, p) = QE(x − p)2 s.t. p ∈ .
(QP)
This problem is at least as difficult as the so-called Closest Vector problem (see [20]). As mentioned before, if the timetable structure allows for complete enumeration of the phase shifts, (QP) can be solved by inspection. Unfortunately, such approach is not always tractable. We may then solve (QP) either exactly or heuristically. If exact optimization of (QP) is too time consuming, we could solve a slight modification of (QP) that uses the Manhattan norm · 1 rather than the Euclidean norm, which reduces (QP) to a mixed-integer linear programming problem. As a computationally more attractive alternative, we describe a local search heuristic for (QP) that is currently used by Railned. Basically, the algorithm iteratively builds a list p1 , p2 , . . . , of feasible phase shifts. We initialize this list by applying any of the available algorithms that finds a feasible timetable structure p 1 (see, e.g., the ones in [10,15,5]). Then, in iteration n, we try to find a feasible phase shift pn+1 ∈ to add to the list that is closer to x than any of the phase shifts found so far. Here, we may either use the distance or its variant using the Manhattan norm. If such a phase shift cannot be found we stop the algorithm with p = pn . In our approach, we try to find an improving phase shift as follows. We first compute 1 + − u fi − fi , i = 1, . . . , M − N uˆ i = T and ˆi =
1 + fi − u fi− , T
i = 1, . . . , M − N,
where fi is the ith fundamental circuit of G relative to the spanning tree T. Note that ˆ p uˆ for all p ∈ and that ˆ and uˆ have to be calculated only once in the course of the algorithm. Choose a parameter
M.A. Odijk et al. / Computers & Operations Research 33 (2006) 2283 – 2299
2295
0 1 which determines the magnitude of the neighborhood in the local search algorithm. Then, in each iteration, sample vectors p and p¯ as follows:
n pi − 1 with probability , pi = i = 1, . . . , M − N, pin otherwise, and
p¯ i =
pin + 1 pin
with probability , otherwise,
i = 1, . . . , M − N.
Finally, try to find a feasible solution p = p n+1 to the system B + T Ep u, p p p, ¯ ˆ p u, ˆ QE(x − p)
min QE(x − p m ) − ,
m=1,...,n
p ∈ ZM−N , where > 0 is some small threshold improvement. If is not too large, then it can be expected that pi = pin = p¯ i for many i = 1, . . . , M − N, which significantly reduces the size of the problem. If desired, the procedure can be repeated with several values of and . One major advantage of this heuristic is that, in subsequent executions of the sampling heuristic, we may reuse the list of phase shifts previously generated, thereby speeding up the algorithm. 4. An example In this section we illustrate our approach to sampling timetable classes by considering the situation at the Dutch railway station Arnhem CS as an example. Details on this problem instance can be found in [5] and are therefore omitted from this paper. The timetable structure has N = 23 free event times and M = 54 interval constraints. We enumerated all elements in , which yielded 1703 different phase shifts. This enumeration of phase shifts is very time-consuming (on the order of minutes per phase shift). However, this step will not be performed in practice, and is here simply included for testing purposes. Further study of the constraints revealed circuit in the constraint graph that corresponds to a set of interval constraints that have to be satisfied as equalities, which implied that the maximum time table class dimension is N − 3 = 20. In particular, it was found that out of the 1703 phase shifts precisely 1168 phase shifts correspond to a timetable class of dimension N − 3, while the remaining 535 have even smaller dimension. We denote these sets by + and − , respectively. Recall that the robustness distribution R assigns positive probability to only the phase shifts in + . To test our approach, we generated 17,030 samples from V (see Proposition 3.2) and rounded these to integral phase shifts using the heuristic described in Section 3.2. The computation time required for generating a sample from V is negligible (on the order of seconds). However, the rounding phase is more time consuming: to avoid the need for tuning the value of we chose to randomly sample this parameter.
2296
M.A. Odijk et al. / Computers & Operations Research 33 (2006) 2283 – 2299
Fig. 1. Results for the robustness distribution applied to Arnhem CS for all 1703 phase shifts together.
The rounding step took an average of 18 iterations and 138 s. We then performed a frequency count for each of the 1703 phase shifts. It turned out that approximately 13% of the phase shifts in + occurred with frequency zero, while this was the case for approximately 40% of the phase shifts in − . This suggests that the heuristic is indeed biased towards avoiding timetable classes with low dimension. Computing the volumes of the 1703 timetable classes is too time consuming. In order to assess the success of our approach in biasing the random generation towards timetable classes containing robust timetables, we determined the distance of each of the sampled phase shifts to the center of the polytope ¯ , which is also the relaxed phase shift at which the volume function v attains its maximum (see Propositions 2.1 and 2.4). These results are plotted in Fig. 1. In this figure, the size of the dot indicates the number of phase shifts that the dot corresponds to. To improve the readability of the figure, seven phase shifts were omitted since they had a frequency of over 85 (the highest frequency observed was 126), although these observations were of course taken into account in all calculations. The horizontal and vertical lines indicate the mean distance and frequency, respectively. This distance measure is based on the intuitive expectation that phase shifts at a large distance from the center of ¯ are assigned less probability mass by R than phase shifts for which this distance is small. This is motivated by the fact that the function v 1/N is concave and the center of ¯ maximizes v. To test this expectation, we have performed a linear regression of the distance on the observed frequency. This regression yielded a slope of −0.0034 with a t-value of 35, which means that the slope is very significantly negative. Figs. 2 and 3 show the data of Fig. 1 partitioned by the phase shifts in + and − , respectively. In this case, performing linear regression of the distance on observed frequency yielded a slope of −0.0030 with a t-value 31 for the data in Fig. 2 and a slope of −0.0032 with a t-value of 15 for the data in Fig. 3. In both cases the negative slope is quite significant.
M.A. Odijk et al. / Computers & Operations Research 33 (2006) 2283 – 2299
2297
Fig. 2. Results for the robustness distribution applied to Arnhem CS for the phase shifts with a timetable set of dimension N − 3.
Fig. 3. Results for the robustness distribution applied to Arnhem CS for the phase shifts with a timetable set of dimension less than N − 3.
Finally, we may observe from the latter two figures that the phase shifts that should be out of sight, i.e. the phase shifts in − , have a much smaller average frequency and a somewhat larger distance to the center of ¯ than the phase shifts in + .
2298
M.A. Odijk et al. / Computers & Operations Research 33 (2006) 2283 – 2299
5. Summary In this paper we studied the problem of randomly generating timetable classes with a bias towards classes that contain robust timetables. This problem is very relevant in performing capacity planning using a scenario analysis to take into account uncertainties in future demand for rail travel. We formalized the notion of a robust timetable and developed a heuristic approach to solve problem of randomly generating classes of robust timetables. Using a practical example, we have illustrated that our method shows promise as it indeed favors classes of robust timetables over other classes. References [1] Odijk MA, Maarschalkerweerd MSPAM, Weits EAG. Timetable—independent design and analysis of railway infrastructure. In: Allan J, Brebbia CA, Hill RJ, Sciutto G, Sone S, editors. Computers in railways V, vol. 1: railway systems and management. Proceedings of CompRail 96, Berlin, Germany; 1996. pp. 3–12. [2] Odijk MA, Zwaneveld PJ, Hooghiemstra JS, Kroon LG, Salomon M. Decision support systems help railned to search for ‘win-win’ solutions in railway network design. Interfaces 1999;29(2):15–32. [3] Hooghiemstra JS. Design of regular interval timetables for strategic and tactical railway planning. In: Allan J, Brebbia CA, Hill RJ, Sciutto G, Sone S, editors. Computers in railways V, vol. 1: railway systems and management. Proceedings of CompRail 96, Berlin, Germany; 1996. pp. 393–402. [4] Odijk MA. Uniform generation of regular interval timetables to facilitate and improve the design of rail infrastructure. In: Hertel G, editor. Schriftenreihe des Instituts für Verkehrssystemtheorie und Bahnverkehr 2. Proceedings of Das Zweite Eisenbahnbetriebswis- senschaftliches Kolloquium, TU Dresden, Dresden, Germany; 1995. pp. 115–26. [5] Odijk MA. A constraint generation algorithm for the construction of periodic railway timetables. Transportation Research B 1996;30(6):455–64. [6] Odijk MA. Sensitivity analysis of a railway stations track layout with respect to a given timetable. European Journal of Operational Research 1999;112(3):517–30. [7] Peeters L, Kroon L. A cycle based optimization model for the cyclic railway timetabling problem. In: Voss S, Daduna JR, editors. Computer-aided scheduling of public transport. Lecture notes in economics and mathematical systems 505. Berlin: Springer; 2001. pp. 275–96. [8] Liebchen C, Möhring RH. A case study in periodic timetabling. In: Proceedings of ATMOS 2002—algorithmic methods and models for optimization of railways, Malaga, Spain. Electronic Notes in Theoretical Computer Science 2002; 66(6). [9] Garey MR, Johnson DS. Computers and intractability—a guide to the theory of NP-completeness. New York: Freeman and Company; 1979. [10] Serafini P, Ukovich W. A mathematical model for periodic scheduling problems. SIAM Journal on Discrete Mathematics 1989;2(4):550–81. [11] Gertsbakh I, Serafini P. Periodic transportation schedules with flexible departure times. European Journal of Operational Research 1991;50:298–309. [12] Serafini P, Ukovich W. A mathematical model for the fixed-time traffic control problem. European Journal of Operational Research 1989;42:152–65. [13] van den Berg JHA, Odijk MA. DONS: Computer aided design of regular service time-tables. In: Murthy TKS, Brebbia CA, Mellitt B, Sciutto G, Sone S, editors. Computers in railways IV, vol. 2: railway operations. Proceedings of CompRail 94, Madrid, Spain; 1994. pp. 109–16. [14] Rockafellar RT. Network flows and monotropic optimization. New York: Wiley; 1984. [15] Voorhoeve M. Rail scheduling with discrete sets. Working paper. Eindhoven University of Technology, Eindhoven, The Netherlands. 1993. [16] Bonnesen T, Fenchel W. Theorie der Konvexen Körper. Berlin: Springer; 1935 (in German). [17] Rubinstein RY. Simulation and the Monte Carlo method. New York: Wiley; 1981.
M.A. Odijk et al. / Computers & Operations Research 33 (2006) 2283 – 2299
2299
[18] Smith RL. Efficient Monte Carlo procedures for generating points uniformly distributed over bounded regions. Operations Research 1984;32(6):1296–308. [19] Luenberger DG. Optimization by vector space methods. New York: Wiley; 1969. [20] Nemhauser GL, Wolsey LA. Integer and combinatorial optimization. New York: Wiley; 1988.