Integrating line planning, timetabling, and vehicle
schedulingInteger programming formulation and analysis∗
Marco Lübbecke

Christian Puchert
Anita Schöbel

Philine Schiewe

February 27, 2019

Abstract
Line planning, timetabling, and vehicle scheduling are three importantstages of public transport planning which highly depend on one another.
It is hence beneficial to solve them in an integrated way instead of sequentially. In this paper we present an integer linear programming formulation
for the integrated line planning, timetabling, and vehicle scheduling problem. Due to the inherent complexity of line planning and timetabling, it is
not possible to solve the integrated model directly. We hence analyze thestructure of the resulting integer program with respect to decomposition
approaches and show that there exist decompositions that are superior
to the canonical decomposition into the planning stages line planning,timetabling, and vehicle scheduling.
Keywords Line planning - Timetabling - Vehicle scheduling - Integrated
public transport planning - Integer programming - Decomposition∗ This work was partially funded by DFG research unit FOR 2083.

M. Lübbecke, C. Puchert
RWTH Aachen University
Kackertstraße 752072 Aachen, Germany
e-mail: luebbecke@or.rwth-aachen.de, puchert@or.rwth-aachen.de
P. Schiewe
University of Goettingen
Lotzestr. 16-18
37083 Göttingen, Germany
e-mail: p.schiewe@math.uni-goettingen.de
A. Schöbel
Technical University of Kaiserslautern
Gottlieb-Daimler-Straße, 67663 Kaiserslautern, Germany
e-mail: schoebel@mathematik.uni-kl.deandFraunhofer Institute for Industrial Mathematics ITWM
Fraunhofer Platz 1
67663 Kaiserslautern, Germany

1

1

Introduction

When planning public transport, three important and well researched stages
are line planning, timetabling, and vehicle scheduling. They are usually solved
sequentially, see, e.g. (Bussieck et al, 1997; Huisman et al, 2005; Desaulniers andHickman, 2007; Guihaire and Hao, 2008). Thus, at first a line plan is constructedthat covers the infrastructure network. Afterwards, arrival and departure times
for all lines and all stops are determined in the timetabling stage minimizing
the travel time of the passengers. For the lines and the corresponding timetablea vehicle schedule is constructed such that the operational costs are minimized.
Due to the nature of the sequential planning process, the solution quality andeven the feasibility of the later problems highly depend on the solution of theearlier problems and the quality of the resulting public transport supply is very
likely to benefit from an integrated approach. Already in (Bussieck et al, 1997)
it is therefore suggested to consider integrated problems spanning more than
one planning stage.
Recently, integration has been the focus of many publications in public
transport planning. The integration of timetabling and passenger routing is
considered in (Siebert and Goerigk, 2013; Schmidt, 2014; Schmidt and Schöbel,
2015a,b; Borndörfer et al, 2016; Gattermann et al, 2016; Robenek et al, 2017)
while line planning and timetabling are considered integratedly in (Schmidt,
2005; Liebchen, 2008; Rittner and Nachtigall, 2009; Kaspi and Raviv, 2013;
Burggraeve et al, 2017; Schiewe, 2018). Integrating timetabling and vehicle
scheduling is considered in (van den Heuvel et al, 2008; Guihaire and Hao,
2010; Ibarra-Rojas and Rios-Solis, 2011; Cadarso and Marı́n, 2012; Petersen
et al, 2013; Schmid and Ehmke, 2015; Yue et al, 2017; Fonseca et al, 2018).
While the integration of several planning stages leads to better solutions, italso increases the computational challenge. Therefore, many approaches to integrating problems in public transport planning are heuristic, such as iterative
approaches in (Kinder, 2008; Siebert and Goerigk, 2013; Burggraeve et al, 2017)or meta-heuristic ones in (van den Heuvel et al, 2008; Guihaire and Hao, 2010;
Petersen et al, 2013; Schmid and Ehmke, 2015; Robenek et al, 2017; Yue et al,
2017; Fonseca et al, 2018). A general scheme for designing iterative algorithms
for integrating line planning, timetabling, and vehicle scheduling is proposedin (Schöbel, 2017), of which parts are implemented in (Michaelis and Schöbel,
2009; Pätzold et al, 2017; Schiewe and Schiewe, 2018).
While the integration of two planning stages is subject of ongoing research,
the integration of all three problems line planning, timetabling, and vehicle
scheduling has not been handled before. We propose an integrated formulation
for all three problems which also includes passenger routing in Section 2 ofthis paper and consider different decompositions in regard of their solvability in
Section 4.

2

Integrating line planning, timetabling, and vehicle scheduling with passenger routing

In this section we shortly describe the usual sequential approach to public transport planning and introduce an IP model for the integrated problem.

2

2.1

Sequential solution

While the sequential planning stages are thoroughly researched, we only give a
brief introduction to each of them. Note that there are many different models
for each of the planning stages which are described in the overview articles
referenced below.
Line planning. The goal of line planning is to cover an infrastructure networkby lines such that a certain travel quality for the passengers is guaranteed andthe costs for the infrastructure provider is not too high. This is usually doneby routing all passengers along shortest paths and assigning lower and upper
frequency bounds on the edges of the infrastructure network that make sure thatthese passengers paths can be realized. For an overview, see (Schöbel, 2012).
Timetabling. We consider periodic timetabling, i.e., the timetable is constructed for a fixed period and then repeated. This is usually modeled as periodic event scheduling problem (PESP), as introduced in (Serafini and Ukovich,1989), which is intrinsically hard to solve, see e.g. (Liebchen, 2007). To avoid integrating passenger routing, it is often assumed that passengers travel on fixed
paths independently from the timetable. For an overview, see (Lusby et al,2011).
Vehicle scheduling. In comparison to line planning and timetabling, vehicle scheduling is an easier problem as its basic form can be solved by a flow
formulation. The goal is to minimize the operational costs by scheduling theoperation of lines by vehicles such that additional costs arising from empty relocation trips between line trips are minimized. It is an aperiodic problem asvehicle schedules do not need to repeat with the same period as the timetable
does. For an overview, see (Bunte and Kliewer, 2009).
Passenger routing. Passenger routing is used in different stages of the sequential approach. It is used to find lower frequency bounds in the line planningstage as well as for defining weights in the timetabling stage. In the objectivefunction it is used to evaluate timetables with respect to the passengers’ travel
times. We have to consider passenger routing in the integrated model, because
assigning passengers to transfers beforehand is impossible when no line plan isfixed.

2.2

Integer programming (IP) model

Our goal is to find a public transport supply consisting of a line plan, a timetable,
a routing of the passengers, and a vehicle schedule that is feasible and minimizes
a weighted sum of the operational costs and the passengers’ travel times. Here,the operational costs contain duration-based and distance-based costs for trips
covering lines and relocations between them as well as costs depending on thenumber of vehicles needed. The structure of the integrated IP model is thefollowing:3

objective: minimize β·travel time + γ·operational costsline planningLfeasibility problem
coupling constraints

LPpass. routing
flow problem

P

coupling constraints

PT

timetabling
PESP

T

coupling constraints

TV

vehicle sched.
flow problem
coupling constraints

VLT, LV

Figure 1: Structure of the integrated line planning, timetabling, and vehicle
scheduling problem with passenger routing.
We abbreviate the four planning stages by L, P, T and V, and hence classify the
coupling constraints as LP, LT, LV, PT and TV. These classes are again usedin the IP model to simplify the presentation of the constraints.
We now develop the integrated integer program in more detail. In order to keep
the complexity of the model manageable, we introduce the following assumptions:• A public transportation network (PTN) (V, E) consisting of stops V anddirect connections (e.g., tracks) E between the nodes is given.
• Lines are chosen from a fixed line pool L0 and are operated with a frequency of one if used. This is encoded in the Boolean variable fl , l ∈ L0 .
• Lower and upper frequency bounds femin , femax on the edges e ∈ E areprovided. The lower frequency bounds femin guarantee a prespecified service quality while the upper frequency bounds femax deal with operational
issues such as headway constraints.
• The timetable is periodic with period length T .
• origin-destination (OD) pairs are routed as a unit on shortest paths according to the travel time.
• The number of periods considered for vehicle scheduling is pmax , we write
the set of considered periods as P = {1, . . . , pmax }.
• The number of vehicles is not limited.
• There is one fixed depot where all vehicles start and end their journey.
• The minimal turnover time between the last event of line l1 and the firstevent of line l2 is given as Ll1 ,l2 .

4The model is based on an event-activity network (EAN) N 0 = (E 0 , A0 ), see(Serafini and Ukovich, 1989; Odijk, 1996), which is derived from the PTN andcontains events and activities for all lines in the pool. The events represent
departures and arrivals of vehicles at stops while the activities represent driving
and waiting of vehicles as well as transferring of passengers between different
lines.
00
E 0 = Earr
∪ Edepwhere0
Earr= {(v, l, arr) : v ∈ l ∩ V, l ∈ L0 }0Edep
= {(v, l, dep) : v ∈ l ∩ V, l ∈ L0 }A0 = A0drive ∪ A0wait ∪ A0trans where
A0drive = {((v1 , l, dep), (v2 , l, arr)) : {v1 , v2 } ∈ l ∩ E, l ∈ L0 }A0wait = {((v, l, arr), (v, l, dep)) : v ∈ l ∩ V, l ∈ L0 }A0trans = {((v, l1 , arr), (v, l2 , dep)) : v ∈ l1 ∩ l2 ∩ V, l1 , l2 ∈ L0 }These events have to be scheduled according to the lower and upper bounds
La , Ua , respectively, on the duration of the activities a ∈ A0 . Therefore, thevariables πi for the periodic time of the events i ∈ E 0 and za for the modulo
parameters on the activities a ∈ A0 are introduced. The auxiliary variables
ya are used to decide whether all lines corresponding to activity a ∈ A0 are
operated. Note that A0 (l1 , l2 ) is the set of activities a = (i, j) such that event i
belongs to line l1 and event j belongs to line l2 while A0 (l) is the set of activities
a = (i, j) such that i or j belongs to l .
To correctly model the passenger routing, the network has to be extended
further to include source and target nodes for all OD pairs which correspond to
nodes in the underlying infrastructure network as well as activities connecting
these special events to the rest of the EAN. These new events need not be scheduled in the timetable. Along the lines of (Schmidt, 2014), we get an extended
event-activity network N̄ = (Ē, Ā) with
0Ē = E 0 ∪ EOD
0EOD= {(u, v, source), (u, v, target) : u, v ∈ V }Ā = A0 ∪ A0to ∪ A0from
A0to = {((u, v, source), (u, l, dep)) : u ∈ l ∩ V, u, v ∈ V }A0from = {((v, l, arr), (u, v, target)) : v ∈ l ∩ V, u, v ∈ V }.
Analogously to the definition of A0 (l), we define Ā(l) as the set of activities
a = (i, j) ∈ Ā such that event i or event j belongs to line l.
As passengers are routed in the extended EAN N̄ in a flow model, we introduce a variable pu,v
for each combination of OD pair (u, v) with u, v ∈ V andaactivity a ∈ Ā indicating whether the OD pair uses the activity.
For vehicle scheduling we introduce Boolean variables x(p1 ,l1 ),(p2 ,l2 ) to indicate whether the p2 -th driving of line l2 is done by the vehicle that directly before
that did the p1 -th driving of line l1 . Similarly, Boolean variables xdepot,(p,l) indicate whether the p-th driving of line l is done by a new vehicle from the depot
and x(p,l),depot indicates whether the vehicle that did the p-th driving of line l is5going to the depot. To correctly describe the p-th driving of line l the following
variables are used: dl is the time it takes in the timetable to get from the firstevent in the line (first(l)) to the last event in the line (last(l)), sp,l , ep,l is thestart or end time of the p-th driving of line l, respectively.
The IP model can now be formulated in the following way, specifying thegeneral structure given in Figure 1.
min β ·

XXpu,va · (πj − πi + za · T )(1)

a=(i,j)∈A0u,v∈V+ γ1 ·

X

Cu,v ·
fl · dl(2)

fl · lengthl

(3)

l∈L0

+ γ2 ·Xl∈L0+ γ3 ·X X  X X
p1 ∈P l1 ∈L0s(p2 ,l2 ) − e(p1 ,l1 ) · x(p1 ,l2 ),(p2 ,l2 )p2 ∈P l2 ∈L0+ xdepot,(p1 ,l1 ) · Ldepot,l1 + x(p1 ,l1 ),depot · Ll1 ,depot
X X  X X

+ γ4 ·
x(p1 ,l1 ),(p2 ,l2 ) · Dl1 ,l2
p1 ∈P l1 ∈L0



(4)

p2 ∈P l2 ∈L0+ xdep,(p1 ,l1 ) · Ddep,l1 + x(p1 ,l1 ),dep · Dl1 ,dep
XX
+ γ5 ·
xdepot,(p,l)
p∈P l∈L0

6



(5)
(6)

X

fl ≥ femin

e∈E

(L1)

fl ≤ femax

e∈E

(L2)

l∈L0 :
e∈l

X
0

l∈L :
e∈l

πj − πi + za · T ≥ ya · L a
πj − πi + za · T ≤ Ua + M · (1 − ya )Au,vya = fl1 · fl2u,v· (pa )a∈Ā = bu,vfl ≥ pu,v
a
X

dl =

a = (i, j) ∈ A0

(T1)

0

a = (i, j) ∈ A

(T2)

0

a ∈ A (l1 , l2 )(LT1)u, v ∈ V(P1)u, v ∈ V, a ∈ Ā(l)

(LP1)

0

(πj − πi + za T )l ∈ L(TV1)

a=(i,j)∈A0 (l,l)

sp,l = p · T + πfirst(l)ep,l = p · T + πfirst(l) + dl

p ∈ P, l ∈ L0

(TV2)

0

(TV3)

p ∈ P, l ∈ L

sp2 ,l2 − ep1 ,l1 ≥ x(p1 ,l1 ),(p2 ,l2 ) · Ll1 ,l2− M 0 · (1 − x(p1 ,l1 ),(p2 ,l2 ) )
X X
fl2 =
x(p1 ,l1 ),(p2 ,l2 )

p1 , p2 ∈ P, l1 , l2 ∈ L0

(V1)

p2 ∈ P, l2 ∈ L0

(LV1)

p1 ∈ P, l1 ∈ L0

(LV2)

p1 ∈P l1 ∈L0+ xdepot,(p2 ,l2 )
X X
fl1 =
x(p1 ,l1 ),(p2 ,l2 )p2 ∈P l2 ∈L0+ x(p1 ,l1 ),depot
x(p,l),• ≤ fl0p ∈ P, l ∈ L(LV3)x•,(p,l) ≤ fl

0

(LV4)

p ∈ P, l ∈ L

πi ∈ {0, . . . , T − 1}, za ∈ Z, ya , fl , pu,va ∈ {0, 1}, dl , sp,l , ep,l ∈ N,
x(p1 ,l1 ),(p2 ,l2 ) , xdepot,(p,l) , x(p,l),depot ∈ {0, 1}

The objective function, see (1) to (6), minimizes a weighted sum of the travel
time of the passengers, see (1), and the operational costs, see (2) to (6). While
the travel time in (1) is weighted with β, the duration of trips in (2) and theirlength in (3) are weighted by γ1 and γ2 , respectively. The duration of relocations
between trips in (4) and their length in (5) are weighted by γ3 and γ4 , respectively while the number of vehicles in (6) is weighted by γ5 . Linearizing thequadratic terms in (1) and (4) leads to adding further variables and constraints.
Equations (L1) and (L2) are the standard feasibility constraints for line
planning, guaranteeing a prespecified service quality as well as a rough upperbound on the costs. Equations (T1), (T2), and (LT1) guarantee the feasibilityof the timetable for the lines which are operated. The passenger flow is modeledby (P1) where a node-arc-incidence matrix is used and equation (LP1) ensures7

that only arcs belonging to operated lines are used by passengers. Equations(TV1), (TV2), and (TV3) model the correct (aperiodic) time for the start andend of the trips which is used in equation (V1) to ensure that the time between
two trips which are operated directly after one another is sufficiently large. The
flow of the vehicles is modeled in equations (LV1) and (LV2) while equations(LV3) and (LV4) ensure that only trips belonging to operated lines are used.
Note that constraint (LT1) can easily be linearized as product of two Boolean
variables.
Note that the formulation given above can furthermore be extended to include the time slice model introduced in (Gattermann et al, 2016) in order to
distribute the favored departure times of the passengers. This is omitted here
to not further complicate the model.

3

Dantzig-Wolfe decompositions

In this section we apply a generic column generation approach to the integrated
line planning, timetabling, and vehicle scheduling problem.
The structure presented in Section 2.2 can be exploited using the so-called
Dantzig-Wolfe decomposition (DWD) (Dantzig and Wolfe, 1960): The problem
is reformulated according to the given structure where each block is represented
as a subproblem. Furthermore, a master problem has the task to select feasible
solutions from each subproblem such that the coupling constraints are satisfied.
Due to the exponentially large number of variables, this master problem issolved by column generation: variables are generated dynamically when solving
the linear relaxation. Embedding this in a branch-and-bound algorithm yields
branch-and-price. For an overview on column generation and branch-and-price,see e.g., (Desaulniers et al, 2005; Vanderbeck and Wolsey, 2010).
The above problem structure, consisting of the subproblems line planning,
passenger routing, timetabling, and vehicle scheduling, seems to be the “canonical” one for applying a DWD. However, any structure that subdivides thecoefficient matrix into blocks and coupling constraints is theoretically suitablefor DWD. Here, two different blocks are independent from one another as theyneither share variables nor constraints. (If linking variables are present, i.e.,variables that are shared by two or more blocks, one can reformulate the problem by adding for each such variable a copy for each block it appears in, andthen introducing coupling constraints that state that the variable copies must
attain the same values, a so-called Langrangian decomposition.) Thus, a broad
variety of structures exist that might be used to decompose the problem. The
questions that arise are:• What other decomposition structures do exist, apart from the canonicalone?
• Is the canonical structure suited best for applying a DWD and performing
branch-and-price? Or, if there exist other structures, does this decomposition-based solution approach perform better on them?
• Are there any properties that can serve as indicators of a good performance?

8

To find other decompositions than the canonical one, we use several structure
detection algorithms, some of them described in (Bergner et al, 2015). Formally,a structure detection algorithm tries to find a mapping C → N0 , where C is theset of constraints. A constraint that is mapped to 0 is a coupling constraint, i.e.,a constraint that belongs to no block, and thus is part of the master problem.
Depending on the detection algorithm, the mapping either already guarantees
that constraints mapped to the same integer form a block or blocks have to beformed by moving variables to the set of linking variables.
A key feature of the detection that we use is that the algorithms are allowed
to determine partial structures C → N0 ∪ {open}; i.e., a constraint can be left
undecided (mapping it to open), and a partial structure that contains undecided constraints can then be completed by another algorithm. In the end, eachcomplete structure (i.e. a structure with no undecided constraints left) usually
has been detected by a combination of several of these detection algorithms,see (Gleixner et al, 2018) for more details. This increases the number of foundstructures and the chance to find suitable decompositions, but leaves the challenge to choose a “meaningful” one with which branch-and-price is expected toperform best.
Very roughly, structure detection proceeds in the following steps:
1. Constraint classifiers determine partitions of C, e.g., according to thenumber of variables and their coefficients. With these partitions, potential
candidates for the number of blocks are determined.
2. Then, partial decompositions are built that only assign certain constraints
to be coupling constraints, but leave the remainder open. This is done in
the following ways:• by the above mentioned constraint classifiers;• by analyzing the densities of the constraints: Constraints with a highnumber of variables are assigned as coupling constraints.
• by graph partitioning: The coefficient matrix A ∈ Rm×n is modeled
as a hypergraph in two different ways:– hyper row graph: Each node represents a column j, and a hyperedge {j : aij 6= 0} for each row i is introduced;– hyper row-column graph: Each node represents a matrix entry(i, j) with aij 6= 0, and each row i is represented by a hyperedge{(i, j) : aij 6= 0} containing its nonzero entries; analogously, there
is a hyperedge for each column j containing its nonzero entries.
Then, graph partitioning algorithms are applied on these graphs.
These graph partitioners yield complete decompositions as well as
partial decompositions which again only assign coupling constraints.
3. The partial decompositions are completed by looking for connected components on the remaining constraints.
4. Last, a postprocessing routine checks if coupling constraints can be assigned to blocks: If a coupling constraint only contains variables of one
block, it will be moved to this block.
A overview on the detection algorithms we use is given in Table 1.
9

Step

LetterAlgorithm

1/2
2
2

c
a
r

constraint classificationgraph partitioning on the hyper row-column graph
graph partitioning on the hyper row graph

3
3

C
d

searching connected components
detection by constraint densities

4p

postprocessing

Table 1: Overview on detection algorithms. The algorithms c, a, and r derivea partial decomposition according to Step 2. They can be followed by the algorithms C and d, see Step 3. The postprocessing algorithm p can be performed
after each of the other algorithms.

4

Computational experiments

The structure detection is implemented in the generic branch-and-price solver
GCG (Gamrath and Lübbecke, 2010) which we use in a development version
based on version 2.1.4. GCG is an extension to SCIP, used in version 4.0, see(Gleixner et al, 2018), a solver for mixed integer programs that also serves as aframework for branch-price-and-cut.
We applied the above structure detection scheme on a small example instancedepicted in Figure 2. Since SCIP comes with various presolving routines which
may change the problem formulation and in particular add new constraints, thedetection scheme was applied twice: first on the original IP formulation, then,after presolving, again on potentially newly added constraints. In total, thisyielded 75 decompositions.
We evaluated each decomposition w.r.t. its computational performance withinbranch-and-price: Therefore, we tried to solve the root LP relaxation within atime limit of one hour. Note that this is the LP relaxation of the master problem,i.e., the subproblems are solved with integrality constraints but the combinationof columns found by the subproblems to a solution of the master problem canbe rational instead of integral. The computations were performed on an Intel(R) Core(TM) i7-2600 CPU at 3.6 GHz, with 16 GB RAM and 8 MB cache,running on openSUSE Leap 42.2 with Linux kernel 4.4. The results are shownin Table 2; for each decomposition, it shows the involved algorithms; moreover,
the relative block and border area, the time and number of LP iterations needed
to solve the root LP relaxation, and the gap between the dual bound and the
optimal solution value of the integrated IP. Here, the optimal value of the IP
can be used to compute the gap as the integrated problem is small enough suchthat it can be solved to optimality by commercial solvers.

4.1

Canonical decomposition

At first, we consider the “canonical” decomposition structure which uses the
subproblems line planning, timetabling, vehicle scheduling, and passenger routing from the sequential process as blocks. It is depicted in Figure 3.
Figure 3a is a reordering of the schematic representation of the matrix struc-

10

v1

v2

v3

v4

Figure 2: Data set used for testing. The solid lines represent the infrastructure
network while the dashed lines represent the lines of the line pool.
coupling constraints
LT, TV, LV, LP, TP
line planning L
feasibility prob.
timetabling
PESP

Tvehic. sched. V
flow problem
pass. routing P
flow problem(a) Schematic representation of the matrix structure, reordering of Figure 1.
0

2000

4000

6000

8000

10000
0

1000

2000

3000

4000

5000

(b) Actual matrix structure. The variables are numbered on the x-axis and the constraints on the y-axis. The small blocks (dark green area) represent the subproblems
timetabling, vehicle scheduling, and passenger routing, with an additional block forline planning which consists of 12 constraints and cannot be seen. The large dark blue
area represents the large number of coupling constraints.

Figure 3: Canonical decomposition of the integrated line planning, timetabling,
and vehicle scheduling problem with passenger routing.
ture given in Figure 1 while Figure 3b represents the actual matrix structure

11

for the instance given in Figure 2. The large area at the top represents thecoupling constraints which clearly make up most of the coefficient matrix thus
making it hard to find subproblems which can be solved independently. Also
note the large number of variables which do not occur in any of the blocks.
These are auxiliary variables used for linearizations of constraints (LT1) and of
the objective and are not explicitly mentioned in the model.
When solving the problem in this canonical form by SCIP, the LP at theroot node of the branch-and-price tree cannot even be solved within the time
limit. Due to the then poor lower bound, the gap is still at 5182.32 % which isfar from optimal.

4.2

Influence of detection algorithms

Therefore, we now consider other decompositions found by GCG. In the following, we denote each combination of detection algorithms by the letters in
Table 1. E.g., apCp means applying graph partitioning and postprocessing
on the un-presolved problem, then searching connected components and again
performing postprocessing on newly added constraints after presolving.
Figure 4 shows the solvability of the matrix structures found by the (combinations of) different algorithms indicated by the gap after solving the root nodeof the branch-and-price tree.

5000

gap in %

4000

3000

2000

1000

r

cCpdCp

cCpCp

apdCp

apdC

apCp

apC

ap

adCp

adC

aCp

aC

a

0

Figure 4: Box and whiskers plot of the performance of the different decomposition algorithms. The algorithms listed on the x-axis are combinations of thedetection algorithms given in Table 1. The boxes mark the 25th to 75th percentile while the whiskers mark the minimal and maximal values. The medianis depicted by a red line. Here, the performance is measured as the gap after
solving the root node of the branch-and-price tree which is given on the y-axis.
Figure 4 suggests that graph partitioning algorithms on hyper row-column
graphs combined with connected components are better suited for the integrated

12

line planning, timetabling, and vehicle scheduling problem then algorithms using
constraint classification or graph partitioning on hyper row graphs.
Especially the algorithms apC, apCp, apdC and apdCp lead to good structures. A typical example of a decomposition found by these algorithms is depicted in Figure 5. Such decompositions are called arrowhead matrices due totheir shape. Intuitively, these decompositions seem to be easier to solve due tothe low number of coupling constraints and variables combined with independent blocks of reasonable sizes.
0
500
1000
1500
2000
2500
3000
3500
4000
4500
0

500

1000

1500

2000

2500

Figure 5: Example for an arrowhead matrix found by the algorithm apC with26 blocks. The variables are numbered on the x-axis and the constraints on the
y-axis. The dark blue area represents the coupling constraints, the purple area
represents the linking variables, and the dark green area represents the blocks.

4.3

Influence of the number of blocks

Figure 6 shows the influence of the number of blocks on the solvability. While
good decompositions could be found for a large span of the number of blocks,
high and low numbers of blocks can also lead to bad decompositions while a
medium number of around 20 to 30 blocks is more promising.

13

5000

gap in %

4000

3000

2000

1000

0
0

20

40

60

80

100

120

number of blocks
a
aC

aCp
adC

adCp
ap

(a) Algorithms a, aC, aCp, adC, adCp and ap, see Table 1.

200
180
160

gap in %

140
120
100
80
60
40
5

10

15

20

25

30

35

40

number of blocks
apC

apCp

apdC

apdCp

(b) Algorithms apC, apCp, apdC, apdCp, see Table 1.

Figure 6: Influence of the number of blocks on the performance of column
generation. Here, the performance is measured as the gap after solving the root
node of the branch-and-price tree which is given on the y-axis.
Note that the scale of Figure 6a and 6b varies as Figure 6b only contains“good” algorithms which lead to a gap of less than 200 percent. Figure 6a also

14

shows an effect which occurred for all decompositions considered here: The gap
is either acceptably small (up to 200%) or the root node LP could not be solved,
leading to a gap of several thousand percent.
Examples for decompositions with a large gap are given in Figure 7. Theyeither feature many very small blocks (Figure 7a and 7b) or few large blocks(Figure 7d) or a combination of both (Figure7c). Medium-sized blocks seem to
be more promising.
0

0

500

500

1000

1000

1500

1500

2000

2000

2500

2500

3000

3000

3500

3500

4000

4000

4500

4500
0

500

1000

1500

2000

2500

0

(a) Algorithm aCp, 106 blocks.

500

1000

1500

2000

2500

(b) Algorithm cCpdCp, 77 blocks.

0

0

500

500

1000

1000

1500

1500

2000

2000

2500

2500

3000

3000

3500

3500

4000

4000

4500

4500
0

500

1000

1500

2000

2500

0

(c) Algorithm cCpCp, 24 blocks.

500

1000

1500

2000

2500

(d) Algorithm r, 3 blocks.

Figure 7: Decompositions with a large gap. The variables are numbered on
the x-axis and the constraints on the y-axis. The dark blue areas represent the
coupling constraints, the purple areas represent the linking variables and thedark green areas represent the blocks.

4.4

Block and border scores

To further characterize the decompositions we consider the block score
PK
block =

k=1 mk · nk

m·n

(7)

and the border score
PK
m0 · n + ( k=1 mk ) · n0
border =
,
m·n

(8)

where K is the number of blocks, m and n are the total number of constraints
and variables, respectively, mk and nk are the number of constraints and variables in block k, respectively, m0 is the number of coupling constraints and n0
the number of linking variables.
In an IP model structure, these two scores indicate the relative block andborder area, respectively. Our expectation is that decompositions with smaller
15

scores lead to a better computational performance, see (Bergner et al, 2015;
Gleixner et al, 2018).
Figure 8 shows the block and border scores for each decomposition. The
canonical decomposition which is depicted by a star differentiates itself by avery high border score compared to all other decompositions considered here
and one of the lowest block scores. This could already be seen in Figure 3 inthe large number of coupling constraints. Decompositions for which the root
LP can be solved such that the gap is less than 200 % all feature both low
block and border scores. Nevertheless, low block and border scores are only anindicator and no guarantee for good solvability, see Bergner et al (2015) for ageneral discussion of such measures.

0.9

0.8

0.7

border score

0.6

0.5

0.4

0.3

0.2

0.1

0.0
0.0

0.1

0.2

0.3

0.4

0.5

block scoregap ≥ 900 %canonical decomposition

gap < 100 %gap < 200 %Figure 8: Comparing different decompositions regarding the block score and theborder score, see equations (7) and (8), respectively.

5

OutlookIn this paper, we developed an integrated formulation for line planning, timetabling, and vehicle scheduling with passenger routing. We analyzed the structure of the resulting integer program with respect to possible decompositions.
The analysis of different matrix structure and their influence on the solvability of the integrated line planning, timetabling, and vehicle scheduling problem
presented here can only be a first step to understanding this problem. Although
we show that the canonical decomposition is far from being an optimal one to16

solve the integrated problem it is hard to determine which decompositions are
suited better to this end. The best indicator so far is a small block score and asmall border score at the same time. Therefore, extensive computational studies are needed to confirm this correlation and find further indicators for good
decompositions.
To improve the efficiency of solving different decompositions, the combination with existing heuristics for problems in public transport planning would beinteresting. Whenever a subproblem is detected that can be interpreted as oneof the problems stemming from the original planning process (i.e., line planning, timetabling, or vehicle scheduling), well researched specialized heuristics
can be used to solve this subproblem. From this we expect a speed-up of thecomputations, leading to better results for solving the integrated problem.
In addition to better understanding the decompositions, it would be interesting to further extend the scope of the integrated problem, e.g., by adding
robustness aspects. It is of particular interest which robustness concept makes
sense in this context (for timetabling, see (Goerigk and Schöbel, 2010) for astudy of different robustness concepts). In our setting, robustness can be considered from a passengers’ and an operators’ point of view. While robust passenger paths depend mainly on transfers which should not break due to small
delays, robust vehicle schedules need should allow for a certain number of delayed vehicles before new ones have to be scheduled.
Furthermore, the integrated model could be extended by additional planningstages. On the one hand, adding network design aspects such as stop location to
the formulation extends the scope of the model in regard to planning completely
new public transport supplies. On the other hand, adding crew schedulingleads to an even better approximation of the operational costs. This, however,increases the computational challenge of the model even more and has to be
combined with solution approaches for reducing the problem size.

Acknowledgements
We would like to thank Michael Bastubbe from RWTH Aachen University fordesigning and implementing the structure detection scheme in GCG, as well as
contributing a number of detection algorithms.

Appendix
Det.

K

block

border# LP iters

Time (s)Gap (%)

a
a
a
a
a
a
a

2
4
8
16
32
3
10

0.50
0.24
0.12
0.06
0.03
0.33
0.10

0.01
0.02
0.03
0.04
0.05
0.01
0.04

274
1913
16564
29172
48081
268
40097

3600.00
3600.00
3600.00
223.2097.90
3600.00
3600.00

4784.504784.504784.5054.30
54.30
4784.504784.50

Continue next page

17

Det.

K

block

border# LP iters

Time (s)Gap (%)

a
a
a
aC
aC
aC
aC
aC
aC
aC
aC
aC
aCp
aCp
aCp
aCp
aCp
aCp
aCp
aCp
adC
adC
adC
adC
adC
adC
adC
adC
adC
adCp
adCp
adCp
adCp
adCp
adCp
adCp
adCp
ap
ap
ap
ap
ap
apC
apC
apC
apC

6
20
100
10
14
22
38
9
12
16
26
106
10
22
38
9
12
16
26
106
5
9
17
33
4
7
11
21
101
5
17
33
4
7
11
21
101
16
32
10
6
100
38
12
16
26

0.16
0.05
0.01
0.27
0.13
0.07
0.03
0.32
0.17
0.10
0.05
0.01
0.27
0.07
0.03
0.32
0.17
0.10
0.05
0.01
0.27
0.13
0.07
0.03
0.32
0.17
0.10
0.05
0.01
0.27
0.07
0.03
0.32
0.17
0.10
0.05
0.01
0.06
0.03
0.10
0.16
0.01
0.03
0.17
0.11
0.05

0.03
0.05
0.10
0.06
0.07
0.09
0.11
0.04
0.12
0.12
0.14
0.18
0.05
0.09
0.11
0.04
0.10
0.11
0.12
0.16
0.05
0.07
0.09
0.11
0.04
0.12
0.12
0.13
0.18
0.05
0.09
0.11
0.04
0.10
0.10
0.11
0.16
0.04
0.05
0.04
0.03
0.10
0.11
0.11
0.08
0.12

18243
46721
1420136276
147856
187441251646273
20195554483625537191589003276
224273265512273
2909775776239686333211795
271
108036138008206704268
2854435743881876435
1232117271
154553259472268
2834976990987213622822062
33182
44427
38957
23019
32615821285930333522811811797093600.00
160.80481.103600.00
2853.20636.70739.303600.00
1415.90904.402829.101463.903600.00
635.90727.403600.00
1639.60
1019.201266.703423.60
3600.00
3600.00
585.00609.603600.00
1414.801030.902067.101703.803600.00
543.40805.803600.00
1813.00
1081.20
1072.60
3518.00
234.7099.00
2621.403600.00
385.80620.801840.00
843.701773.904784.5054.30
4784.504784.5054.40
113.20190.804784.50
57.20
56.80
59.80
4784.504784.50
113.20190.804784.50
57.20
56.60
59.70
4784.504784.504784.50
113.20190.804784.50
57.20
56.80
59.80
1064.30
4784.50
113.20190.804784.50
57.20
56.50
59.70
1050.20
54.30
54.30
54.30
4784.50
96.70
190.80
56.90
55.10
58.80

Continue next page

18

Det.
apC
apCp
apCp
apCp
apCp
apdC
apdC
apdC
apdC
apdC
apdCp
apdCp
apdCp
apdCp
cCpCp
cCpCp
cCpCp
cCpCp
cCpCp
cCpdCpcCpdCpr

K

blockborder# LP iters

Time (s)Gap (%)

106
0.01
0.12
1642579
2619.60
38
0.03
0.11
247731760.4012
0.17
0.103325211776.50
26
0.05
0.10
5777821052.80
106
0.01
0.12
18004972891.40
33
0.03
0.11
235631670.507
0.17
0.11330598
1882.80
11
0.11
0.08
348440836.6021
0.05
0.111327333
1828.00
101
0.01
0.12
740547
1826.1070.17
0.102792362162.70
11
0.11
0.08
315679908.8021
0.05
0.10571431
978.20
101
0.01
0.12
994508
1950.80
24
0.50
0.05
292
3600.00
24
0.49
0.05
902
3600.00
24
0.11
0.15
308
3600.00
24
0.49
0.05
472
3600.00
24
0.11
0.15
322
3600.00
19
0.10
0.15
392
3600.00
19
0.10
0.15
1460
3600.00
3
0.22
0.33
10928
3600.00
Table 2: Detailed results for each decomposition

910.40190.8056.90
58.80
910.40190.80
56.90
55.10
58.80
910.40
56.90
55.10
58.80
910.404784.504784.504784.504784.504784.504784.504784.504784.50

References
Bergner M, Caprara A, Ceselli A, Furini F, Lübbecke ME, Malaguti E, Traversi
E (2015) Automatic Dantzig–Wolfe reformulation of mixed integer programs.
Mathematical Programming 149(1-2):391–424
Borndörfer R, Hoppmann H, Karbstein M (2016) Passenger routing for periodic
timetable optimization. Public Transport DOI 10.1007/s12469-016-0132-0
Bunte S, Kliewer N (2009) An overview on vehicle scheduling models. Public
Transport 1(4):299–317
Burggraeve S, Bull S, Vansteenwegen P, Lusby R (2017) Integrating robust
timetabling in line plan optimization for railway systems. Transportation Research Part C: Emerging Technologies 77:134–160
Bussieck M, Winter T, Zimmermann U (1997) Discrete optimization in public
rail transport. Mathematical Programming 79(1-3):415–444
Cadarso L, Marı́n Á (2012) Integration of timetable planning and rolling stock
in rapid transit networks. Annals of Operations Research 199(1):113–135
Dantzig GB, Wolfe P (1960) Decomposition principle for linear programs. Operations Research 8(1):101–111

19

Desaulniers G, Hickman M (2007) Public transit. Handbooks in Operations
Research and Management Science 14:69–127
Desaulniers G, Desrosiers J, Solomon MM (eds) (2005) Column Generation.
Springer
Fonseca J, van der Hurk E, Roberti R, Larsen A (2018) A matheuristic fortransfer synchronization through integrated timetabling and vehicle scheduling. Transportation Research Part B: Methodological 109:128–149
Gamrath G, Lübbecke ME (2010) Experiments with a generic Dantzig-Wolfe
decomposition for integer programs. In: Experimental Algorithms, Springer,
LNCS, vol 6049, pp 239–252
Gattermann P, Großmann P, Nachtigall K, Schöbel A (2016) Integrating passengers’ routes in periodic timetabling: A SAT approach. In: ATMOS
2016, OpenAccess Series in Informatics (OASIcs), vol 54, pp 1–15, URL
http://drops.dagstuhl.de/opus/volltexte/2016/6527
Gleixner A, Bastubbe M, Eifler L, Gally T, Gamrath G, Gottwald RL, Hendel
G, Hojny C, Koch T, Lübbecke ME, Maher SJ, Miltenberger M, Müller B,
Pfetsch ME, Puchert C, Rehfeldt D, Schlösser F, Schubert C, Serrano F,
Shinano Y, Viernickel JM, Walter M, Wegscheider F, Witt JT, Witzig J(2018) The SCIP Optimization Suite 6.0. ZIB-Report 18-26, Zuse-Institute
BerlinGoerigk M, Schöbel A (2010) An empirical analysis of robustness concepts for timetabling. In: Erlebach T, Lübbecke M (eds) Proceedings of ATMOS10, Schloss Dagstuhl–Leibniz-Zentrum fuer Informatik,
Dagstuhl, Germany, OpenAccess Series in Informatics (OASIcs), vol 14,
pp 100–113, DOI http://dx.doi.org/10.4230/OASIcs.ATMOS.2010.100, URL
http://drops.dagstuhl.de/opus/volltexte/2010/2753
Guihaire V, Hao J (2008) Transit network design and scheduling: A global
review. Transportation Research Part A: Policy and Practice 42(10):1251–
1273
Guihaire V, Hao JK (2010) Transit network timetabling and vehicle assignment
for regulating authorities. Computers & Industrial Engineering 59(1):16–23
van den Heuvel A, van den Akker J, van Kooten M (2008) Integrating
timetabling and vehicle scheduling in public bus transportation. Tech. rep.,
Utrecht University
Huisman D, Kroon L, Lentink R, Vromans M (2005) Operations research inpassenger railway transportation. Statistica Neerlandica 59(4):467–497
Ibarra-Rojas O, Rios-Solis Y (2011) Integrating synchronization bus timetablingand single-depot single-type vehicle scheduling. In: ORP3 Meeting, Cadiz
Kaspi M, Raviv T (2013) Service-oriented line planning and timetabling forpassenger trains. Transportation Science 47(3):295–311
Kinder M (2008) Models for periodic timetabling. Master’s thesis, Technische
Universität Berlin20

Liebchen C (2007) Periodic timetable optimization in public transport. Springer
Liebchen C (2008) Linien-, Fahrplan-, Umlauf- und Dienstplanoptimierung: Wieweit können diese bereits integriert werden? In: Heureka’08
Lusby R, Larsen J, Ehrgott M, Ryan D (2011) Railway track allocation: models
and methods. OR Spectrum 33(4):843–883
Michaelis M, Schöbel A (2009) Integrating line planning, timetabling, and vehicle scheduling: A customer-oriented approach. Public Transport 1(3):211–232
Odijk M (1996) A constraint generation algorithm for the construction of periodic railway timetables. Transportation Research Part B: Methodological30(6):455–464
Pätzold J, Schiewe A, Schiewe P, Schöbel A (2017) Look-Ahead Approaches
for Integrated Planning in Public Transportation. In: D’Angelo G, Dollevoet
T (eds) 17th Workshop on Algorithmic Approaches for Transportation Modelling, Optimization, and Systems (ATMOS 2017), Schloss Dagstuhl–LeibnizZentrum fuer Informatik, Dagstuhl, Germany, OpenAccess Series in Informatics (OASIcs), vol 59, pp 17:1–17:16
Petersen H, Larsen A, Madsen O, Petersen B, Ropke S (2013) The Simultaneous
Vehicle Scheduling and Passenger Service Problem. Transportation Science
47(4):603–616
Rittner M, Nachtigall K (2009) Simultane Liniennetz- und Fahrlagenoptimierung. Der Eisenbahningenieur
Robenek T, Azadeh S, Maknoon Y, Bierlaire M (2017) Hybrid cyclicity: Combining the benefits of cyclic and non-cyclic timetables. Transportation Research Part C: Emerging Technologies 75:228–253
Schiewe A, Schiewe P (2018) An Iterative Approach for Integrated Planningin Public Transportation. Tech. rep., Georg-August-Universität Göttingen,working Paper
Schiewe P (2018) Integrated optimization in public transport planning. PhD
thesis, Georg-August-Universität Göttingen
Schmid V, Ehmke JF (2015) Integrated timetabling and vehicle scheduling withbalanced departure times. OR Spectrum 37(4):903–928
Schmidt D (2005) Linien- und Taktfahrplanung - Ein integrierter Optimierungsansatz. Master’s thesis, Technische Universität Berlin, in German
Schmidt M (2014) Integrating Routing Decisions in Public Transportation Problems, Optimization and Its Applications, vol 89. Springer
Schmidt M, Schöbel A (2015a) The complexity of integrating routing decisions
in public transportation models. Networks 65(3):228–243
Schmidt M, Schöbel A (2015b) Timetabling with passenger routing. OR Spectrum 37:75–97

21

Schöbel A (2012) Line planning in public transportation: models and methods.
OR Spectrum 34(3):491–510
Schöbel A (2017) An eigenmodel for iterative line planning, timetabling and vehicle scheduling in public transportation. Transportation Research C 74:348–
365
Serafini P, Ukovich W (1989) A mathematical model for periodic schedulingproblems. SIAM Journal on Discrete Mathematics 2(4):550–581
Siebert M, Goerigk M (2013) An experimental comparison of periodic
timetabling models. Computers & Operations Research 40(10):2251–2259
Vanderbeck F, Wolsey LA (2010) Reformulation and decomposition of integer
programs. In: 50 Years of Integer Programming 1958-2008, Springer, pp 431–
502
Yue Y, Han J, Wang S, Liu X (2017) Integrated Train Timetabling and Rolling
Stock Scheduling Model Based on Time-Dependent Demand for Urban RailTransit. Computer-Aided Civil and Infrastructure Engineering 32(10):856–873

22