Math Meth Oper Res (2004) 59:205–220
DOI 10.1007/s001860300332

A fast algorithm for near cost optimal line plans
Michael R. Bussieck1 , Thomas Lindner2 , Marco E. Lübbecke3
1
GAMS Development Corporation, 1217 Potomac St. NW, Washington DC 20007(e-mail: MBussieck@gams.com)
2
Siemens Transportation Systems, Ackerstraße 22, D-38126 Braunschweig, Germany(e-mail: Thomas.Lindner@ts.siemens.de)3Technische Universität Berlin, Institut für Mathematik, Sekr. MA 6-1, Straße des 17. Juni 136,
D-10623 Berlin, Germany (e-mail: m.luebbecke@math.tu-berlin.de)

Manuscript received: May 2003/Final version received: November 2003

Abstract. We consider the design of line plans in public transport at a minimal total cost. Both, linear and nonlinear integer programming are adequate
and intuitive modeling approaches for this problem. We present a heuristic
variable ﬁxing procedure which builds on problem knowledge from bothtechniques. We derive and compare lower bounds from diﬀerent linearizations in order to assess the quality of our solutions. The involved integer
linear programs are strengthened by means of problem speciﬁc valid
inequalities. Computational results with practical data from the Dutch
Railways indicate that our algorithm gives excellent solutions within minutesof computation time.
Key words: Integer programming, Nonlinear integer programming, Cutting
planes, Public transport, Line planning
MSC (2000): 90C10, 90C30, 90C57, 90B06

1 Introduction
In a public transportation network, we refer to a line as a route, or itinerary, together with the speciﬁcation of how often the route is operated. Theline planning problem is the design of lines which meet several operational
constraints, most notably the passengers’ demand for transportation. The
resulting line plan, or route map, is visually known to anyone who evertraveled by bus, tram, subway, or train. Two objectives have been considered in the literature, and these reﬂect both sides of the story, namely
service versus cost aspects. While travelers demand for convenient, ideallydirect connections [2, 3], transportation companies are forced to make mosteﬃcient use of their resources [2, 6, 8, 11], not least for the reason of market

206

Michael R. Bussieck et al.

deregulations. Minimization of the total cost is the goal of our investigation
as well.
One main issue of our paper is to further the integration of two areas inmathematical programming all too often regarded separately, viz. nonlinear
and mixed integer optimization. While the latter is a well established technique in the planning of public rail transportation issues [4], the track
towards the topic via the former is practically unbeaten. In fact, mixed integer
nonlinear models and algorithms appear not to belong to the ‘‘traditional’’operations research active areas. This contrasts the tenor of the OR/MS
Today Special Issue on innovative education: Nonlinear optimization
deserves more emphasis [9].
In the cost optimal line planning context, Goossens et al. [8] demonstrate how to use lower bounds from a binary linear program [6] in abranch-and-cut algorithm. While the quality of lower bounds is excellent,primal solutions of matching quality cannot be found in reasonable time.
In fact, the emphasis of previous exact approaches is on lower bounds.
Most recently, Fischetti and Lodi [7] developed a general strategy for
obtaining feasible solutions to mixed integer programs, and also makesome computational progress on our particular problem. The main purpose of this paper is to justify that an appropriate combination of twomodels is suited to yield very high quality solutions within very shortcomputation time limits. After all, this and only this – is relevant topractitioners.
Our paper delivers updated material from [2] which has already been built
on by other authors [8]. We contribute linear and nonlinear integer programs,
on which we base our variable ﬁxing algorithm. Our successful use of anonlinear model is a new and surprising development in the design of line
plans.
2 Cost optimal line planning
In this section we formally describe the cost optimal line planning problem
and introduce our notation. We do not repeat the practical background, andwe refrain from comparing cost minimization with the direct traveler
approach; see e.g., [3, 6] for details on both issues.
The underlying rail network is represented by a graph G ¼ ðV ; EÞ. Edges
e 2 E model physical links which connect stations given by the set V ofvertices, cf. Figure 1 for an example. Only a predeﬁned subset R of paths,
and sometimes cycles, in G qualiﬁes as possible itineraries along which
trains can be operated. The design of a line l ¼ ðr; uÞ requires the choice ofa route r 2 R and a frequency u 2 F  Zþ , i.e., the number of times theroute is operated per hour. Unless otherwise stated, we assume that F is aninterval of (small) integers, including zero. The length of route r 2 R in
kilometers is denoted by kmr . A feasible set of lines, or line plan, must alsoobey lower and upper bounds Femin and Femax on the line frequency
requirement, i.e., on the total number of trains which traverse a link e 2 E
per hour. In addition, for each edge e 2 E we are given its traﬃc load Le ,that is, the total number of seats per hour which are to be provided on thelink. Our study is driven by the assumption that the cost of a line planessentially consist of personnel cost and the cost of rolling stock. ThisA fast algorithm for near cost optimal line plans

207

Fig. 1. The Netherlands’ InterRegio network with 86 vertices and 113 edges

motivates making the number of coaches of each train a decision variable
for every particular line. This number must range within given bounds c
and c. It is important to note that we consider a single type of coach only,with a single seat capacity cap. Therefore, it is more natural to express thetraﬃc load on an edge in terms of coaches, rather than in terms of of seats.
That is, we deﬁne Ke :¼ dLe =cape, the number of coaches to be provided onlink e.
We denote by Ctfix and Ccfix , respectively, the hourly ﬁxed capital andpersonnel cost incurred per train/per car. The respective operational or variable cost, that is, the per kilometer cost are denoted by Ctvar and Ccvar . The
cost optimal line planning problem is to design a feasible line plan at minimum total cost. This problem is strongly N P-hard [2].
3 Primal solutions from a nonlinear model
Perhaps the most intuitive model formulation for our problem is a nonlinear
mixed integer program (MINLP). Such a formulation was already given in
Claessens’ master’s thesis [5] but was later dropped for the publication [6]because all presented results came out of a linearized model. We develop new

208

Michael R. Bussieck et al.

ideas to overcome some of the technical obstacles which motivates us torevisit the MINLP approach.
Nonlinear mixed integer programming is a comparatively recent areawith the ﬁrst general purpose algorithms appearing in the early 1990s(GAMS/DICOPT [10]). Recently, several implementations emerged (branchand-bound, outer approximation, extended cutting plane) due to morestable nonlinear programming codes (FILTERBB, SBB, AlphaECP,
mittlp, BARON (global), other MINLP solvers embedded in modelinglanguages like LINGO and MINOPT1 ). An increase of interest in MINLP
models and their solution is reﬂected in the growing number (absolute andrelative) of entries at the NEOS2 site. Within six months (09/01–02/02),
their MINLP section got about 2,500 user submissions (private communication).
Practically relevant nonlinear models have some drawbacks comparedto their linear counterparts: There is practically no proof of global optimality due to the non-convexity of practical models; the starting point is
critical; in order to prevent algorithms from failing, problem speciﬁcknowledge is essential for a good model formulation, and thus, for algorithmic guidance. We have chosen our nonlinear approach in view of these
considerations.
For each route r 2 R we introduce a variable xr 2 Zþ which reﬂects its
frequency and an integer variable c  yr  c for the number of coaches per
train. The requirement that the line plan provide suﬃcient capacity for eachlink in the network then naturally writes down asXxr  yr  Ke
8e 2 E:
ð1Þ
r2R;r3e

Similar quadratic terms appear in the objective function. However, for theactual number of trains per hour needed to operate a line at a given frequency we have to count more carefully. It depends on the total duration ofa route, including its minimum turn-around time needed e.g., for cleaning
the train, maintenance, and changing the crew. A division of this durationin minutes by 60 yields a fractional estimate Tr of the number of trainsneeded to operate r 2 R once per hour. The overall requirement on thatroute is then dTr  xr e trains. Although undeﬁned derivatives of such discontinuous terms are handled by solvers, modelers are advised against using
them. Instead, deﬁning variables zr 2 Zþ , r 2 R via
Tr  x r  zr

8r 2 Rð2Þ

we theoretically can get over this problem. Because of our minimization
objective there is no need for an explicit upper bound on these variables. Bymeans of a translation we eliminate the non-trivial lower bounds on yr . Thatis, new variables yr count the number of additional coaches in excess of theminimal number c on route r 2 R. The complete MINLP model reads asfollows:

1

Solver information available on the web under http://www.gamsworld.org/minlp/
links:htm, http://at8.abo.fi/~hasku/mittlp/, and http://www.lindo.com/
table=lgosolvet:html
2
URL http://www-fp.mcs.anl.gov/otc/Guide/ SoftwareGuide=A fast algorithm for near cost optimal line plans

minimizeX

209

zr  ðCtfix þ ðc þ yr Þ  Ccfix Þ þ kmr  xr  ðCtvar þ ðc þ yr Þ  Ccvar Þ

r2R

subject toX

xr  Femin

8e 2 E

r2R; r3eX

xr  Femax 8e 2 E

r2R; r3eX

xr  ðc þ yr Þ  Ke

8e 2 ETr  x r  zr

8r 2 R

r2R; r3eðNLPÞyr  c  c 8r 2 R
xr 2 F

8r 2 Ry r 2 Zþ8r 2 Rz r 2 Zþ8r 2 R

Solving this model leaves us with the dilemma of local versus global optimality. Implementations of global optimization algorithms like BARON do
not ﬁnd useful bounds and suﬀer from the same diﬃculties as linear branchand-bound codes for ﬁnding integer solutions for real-world sized instances.
Even large scale local optimization solvers perform poorly on these instances. Our ﬁrst remedy to make this model useful for the cost optimal line
planing problem simpliﬁes the calculation of the number of trains required.
Since
dTr  xr e  dTr e  xr

ð3Þ

holds for xr 2 Zþ we may substitute zr for the right hand side of (3), thereby
getting rid of the additonal variables as well as their deﬁning inequalities. This
modiﬁcation aﬀects the ﬁxed cost term in the objective function only. For ourdata sets, the eﬀect is mild since we are not given a Ctfix . Still, the level of error
induced by the overestimation clearly depends on the instance. With hourly
frequencies, i.e., F ¼ f0; 1g, which are common in German and Dutch InterCity networks this approach is exact. For our instances, quite a lot of
inequalities (3) are strict, but comparing the exact cost with the modiﬁed cost
for all connections in our solutions, we obtain an average relative error of at
most 3%. This is explained by the fact that, in our experiments, the best
found solutions have almost all lines with a frequency of one anyway. Notethat an obvious advantage of this formulation is its small number of variables.
4 Lower bounds from linearizations
Since (NLP) is a non convex model, the ‘‘lower bound’’ obtained from solving
its continuous relaxation with available local optimization solvers is mathematically useless. As a remedy to this shortcoming we present three linearizations. The ﬁrst essentially is a proposal by Claessens et al. [6] with an210

Michael R. Bussieck et al.

enormous number of variables. We introduce two new models which make up
for this defect.
4.1 A binary linear program
The fact that all nonlinear terms in (NLP) are of the form xr yr and all variables are integer suggests a linearization by considering only the meaningfulvalues of the respective variable products. More precisely, introduce a binary
variable zr;u;c which assumes a value of one if and only if route r 2 R is
operated with frequency u 2 F using c 2 C :¼ fc; . . . ; cg coaches [6]. With anadditional constraint
PwhichPensures that exactly one combination is selected
we replace xr yr by u2F c2C u  c  zr;u;c . Note that this substitution eveneliminates the discontinuous term from the objective function. The modelreads
X X X
minimizedTr  ue  ðCtfix þ c  Ccfix Þ þ kmr  u  ðCtvar þ c  Ccvar Þ  zr;u;cr2R u2F c2C

subject to

X XX

u  zr;u;c  Femin 8e 2 E

r2R; r3e u2F c2C

X XX

u  zr;u;c  Femax 8e 2 E

ðBLPÞr2R; r3e u2F c2C

X XX

u  c  zr;u;c  Ke

8e 2 E

r2R; r3e u2F c2C

XX

zr;u;c  1

8r 2 Ru2F c2C

zr;u;c 2 f0;1g 8r 2 R;u 2 F ;c 2 C
Model (BLP) is more ﬂexible than (NLP) with respect to important
operational constraints which result in holes in the domains of variables.
For instance, F ¼ f1; 2; 4g is a reasonable choice where explicitly 3 2= F.
Moreover, a technical particularity of the rolling stock used by the Dutch
Railways does not allow for trains made up of exactly ﬁve coaches, while
C ¼ f3; 4; 6; 7; . . .g is allowed. Such restrictions entail additional constraintsin (NLP) but are easily incorporated in (BLP) simply by dropping therespective variables.
4.2 An integer linear program based on frequencies
A drawback of (BLP) clearly is the comparatively large number of variables.
Using integer variables, we are able to introduce a more compact new formulation. Similar to (NLP) we state the model in the formulation withadditional coaches. Binary variables xr;u indicate whether a route/frequency
combination is chosen, and if so, variables yr;u reﬂect the number of coaches
more than c.

A fast algorithm for near cost optimal line plans

minimizeXX

211

dTr  ue  ðxr;u  Ctfix þ ðc  xr;u þ yr;u Þ  Ccfix Þ

r2R u2F

þ kmr  u  ðxr;u  Ctvar þ ðc  xr;u þ yr;u Þ  Ccvar Þ
X X
subject tou  xr;u  Femin 8e 2 E
r2R; r3e u2F

X

X

X

u  xr;u  Femax

8e 2 E

r2R; r3e u2FXu  ðc  xr;u þ yr;u Þ  Ke

8e 2 E

r2R; r3e u2FðILPÞ

yr;u  ðc  cÞ  xr;u  0 8r 2 R; u 2 FXxr;u  1 8r 2 R
u2F

xr;u 2 f0; 1g 8r 2 R; u 2 F
yr;u 2 Zþ

8r 2 R; u 2 F

4.3 A third linearization based on capacities
Similarly, one might be tempted to come up with a third linearization. Here, abinary variable xr;c selects the combination of a route and the corresponding
capacity while an integer variable yr;c indicates the frequency of the respective
combination. Even though the constraints of this model are similar to (ILP)
their structure is slightly more simple. Most notably, we can state the bounds
on the line frequency requirement with constraints having 0/1 coeﬃcients
only. This happens at the cost of the discontinuous term reappearing in theobjective function, and the necessary workaround in spirit of Equation (2).
We may drop the ceiling when using this model for the purpose of lower
bounding only.
XX
dTr  yr;c e  ðCtfix þ c  Ccfix Þ þ kmr  yr;c  ðCtvar þ c  Ccvar Þ
minimizer2R c2CX X
subject toyr;c  Femin 8e 2 E
r2R; r3e c2CXX

yr;c  Femax 8e 2 E

r2R; r3e c2CXXc  yr;c  Ke 8e 2 E

ðILP-cÞ

r2R; r3e c2C

yr;c  umax  xr;c  0 8r 2 R; c 2 CXxr;c  1 8r 2 R
c2C

xr;c 2 f0; 1g 8r 2 R; c 2 C
yr;c 2 F 8r 2 R; c 2 C
This model has another interesting property in the case that only one frequency (besides zero) is admissible. Again, we may substitute the discontinuous term in the objective function for dTr e  yr;c but we can also completely

212

Michael R. Bussieck et al.

drop the yr;c variables including their deﬁning constraints from the formulation
altogether. This is because for F ¼ f0; 1g, i.e., umax ¼ 1, together withPx
c2C r;c  1, r 2 R we have yr;c  1  xr;c anyway and we substitute xr;c for
yr;c . The remaining model is (BLP) with ﬁxed u  1! The emphasis withrespect to important decisions in this model clearly lies on capacities rather
than on frequencies. In a sense, among the three proposals this is the linearization closest to (NLP).
5 Valid inequalities
In strengthening our linearizations by valid inequalities we aim at improving
the quality of lower bounds. We present valid inequalities for (ILP) only. Our
argumentation is based on the problem rather than on models, and all resultsof this section immediately carry over to the other linearizations. We prefergiving the intuition behind our results only. The reason for doing so is todemonstrate that rather fancy looking cutting planes are nothing more (andless) but good problem knowledge. Detailed proofs and some advanced
preprocessing techniques canPbe found in [2].
Let us abbreviate KE0 :¼ e2E0 Ke . We will now generalize the following
idea. Let e; f 2 E be the only edges incident to v 2 V , and let Ke < Kf . Whenno line via f ends in v, the number of coaches running via e is at least Kf .
Proposition 1. Let E0  E, f 2 EnE0 , arE0 ¼ jE0 \ rj, and Kf > KE0 . Then,X
X
XX
xr;u þ
arE0  u  ðc  xr;u þ yr;u Þ  Kf ð4ÞðKf  KE0 Þ
r2R;arE0 ¼0;r3f u2F

r2R u2F

is a valid inequality for (ILP).
The following two classes are particularly useful in our experiments. A
tentative explanation is that they exploit observations on capacities and frequencies. Consider an edge e 2 E with c  Femin < Ke < cðFemin þ 1Þ. Then, thenumber of trains via e is at least Femin þ 1 or the number of additional coachesof lines via e is at least n :¼ Ke  c  Femin .
Proposition 2. For every edge e 2 E with c  Femin < Ke < cðFemin þ 1ÞX X
n  u  xr;u þ minfn; ug  yr;u  n  ðFemin þ 1Þð5Þ

r2R;r3e u2F

is a valid inequality for (ILP).
The line frequency requirement impacts the total number of coaches on anedge due to the lower bound c on the number of coaches per train. Suppose
the line plan contains a line ðr; uÞ with e 2 r. Independently of the particulardemand Ke there must be at least cðFemin  uÞ coaches of other lines via e.
This is generalized in the following result which is the most eﬀective in our
experience.
00
Proposition
P
P 3. Let e 2 E and R  R with e 2 r for all r 2 R and
r2R0
u2F xr;u  1 in any feasible solution to (ILP). Then,A fast algorithm for near cost optimal line plans

X

X

213

u  ðc  xr;u þ yr;u Þ  Ke 1 

r2RnR0 ;r3e u2F

XX

!
xr;u

r2R0 u2F

þc

XX
r2R

ðFemin  uÞ  xr;u

ð6Þ

0 u2F

is a valid inequality for (ILP).
Our ﬁnal proposal involves edges e 2 E with a large demand Ke but with asmall minimum line frequency requirement Femin . When the number of trains
via e is close to Femin the number of coaches in the operated lines must be close
to c. The following is a generalization of this statement. Observe that in ourmodels we count coaches in addition to c.
Proposition 4. Let e 2 E and ðr ; u Þ be a line with e 2 r . Deﬁne f :¼ Ke cðFemin  maxf1; Femin  Femax þ u gÞ  u c. If f  0, then0
1Xu  yr ;u  @Femin u  xr;u Adf=u e  0ð7Þðr;uÞ6¼ðr ;u Þ;r3e

is a valid inequality for (ILP).
We close this section with a few words on separation, i.e., on identifying
valid inequalities which are violated by a given fractional solution to (ILP).
We cannot give a general scheme for separating (4) for all E0  E which
actually is N P-hard [2]. However, our experiments show that inequalities
with jE0 j ¼ 1 dominate those with larger jE0 j. The number of inequalities (5)and (7) is polynomial by deﬁnition. Also, in (6) we consider the special case ofjR0 j ¼ 1 only. Hence, for our instances, we check only a polynomial numberof inequalities. In particular, we separate valid inequalities of Propositions 1–
4 in a cut-and-branch fashion, i.e., in the root node of the branch-and-boundtree only, see Algorithm 1. In the remainder of this paper, we refer to thislatter model (ILP) augmented by valid inequalities as model (CUT). For abranch-and-cut approach based on (BLP) see [8].
6 Turning models into an algorithm
We already remarked that optimally solving our nonlinear model is not an
option. Even our simpliﬁed version of (NLP) using (3) represents a challenge
to existing MINLP codes. The problem is still far from being trivial. This canbe deduced from the performance of the until recently only method to solve
MINLPs: A linearization at the (local) optimum point of the relaxed MINLP
results in a mixed integer program and is solved as a second step in the outerapproximation algorithm of DICOPT. Solving only the MIP by CPLEX 7.5
takes over four hours.
Lagrangian relaxation yields unsatisfactory results on the original MINLP
model by Claessens [5]. She proposes an iterative rounding heuristic along thefollowing lines. Promising routes with xr yr larger than a given threshold areincorporated in the line plan by bounding xr  1 from below. On the other
hand, routes with small xr yr value are eliminated from further inspection.
When all routes are either deleted or incorporated, all xr variables are214

Michael R. Bussieck et al.

rounded up to the next integer. The resulting program is an integer linear
program in the yr variables. An alternative to this approach is to ﬁx yrvariables instead. In the resulting integer linear program one has ﬁxed
capacities on each route and suﬃcient seat capacity needs to be provided foreach link at minimal cost. This problem remains N P-hard even for
Femin  Femax  1, e 2 E [2]. All these earlier proposals failed in ﬁnding goodprimal solutions.
For large scale mixed integer programs in general, and for our (integer)
linearizations in particular, we cannot hope for an optimal solution within
practically relevant time bounds. One would terminate a branch-and-boundalgorithm when a predeﬁned computation time limit is reached. A drawback
is the risk of not having found any primal feasible solution at that time.
Modern commercial branch-and-bound codes optionally apply sophisticatedstrategies and invest quite some time in the root node in order to come up
with heuristic solutions. Commonly used variable rounding fails for our
instances. The latest CPLEX 8.0 fares better in that it implements various MIPemphasis preferences to accomodate, among others, a bias towards ﬁnding
feasible solutions. Fischetti and Lodi [7] develop a general two-level
branching framework which controls feasibility at a higher ‘‘strategic’’ level.
Computationally, this method outperforms CPLEX’ defaults. A considerable
improvement as compared to relying on general strategies is to develop a
problem speciﬁc algorithm which exploits the models we introduced earlier.
This is what we do now.
Fixing all lines according to Claessens’ procedure is a too restrictive
measurement, and makes permanent decisions too early. Our idea is to use
model (NLP) to get a hint to promising lines and pass this information tomodel (ILP), still leaving enough freedom to work on these. At ﬁrst, we invest
in a meaningful starting point for (NLP). To this end, we solve the linear
programming relaxation of (ILP). Variables in (NLP) are initialized according to their transformed counterparts in the linearized model’s optimal
solution. A warm or hot start with a former good solution is possible as well.
Then, the continuous relaxation of (NLP) is solved. In the solution, a variablexr close to zero is interpreted as r being an unpromising route, and all
variables xr? ;u or lines involving such a route are eliminated from model (ILP)
by ﬁxing these variables to zero. All other variables are left untouched. Wethen succeed in ﬁnding an integer solution to the partially ﬁxed (ILP) model
within three minutes.
From this point on, we have a feasible solution and we may invest in either
improving it or assessing its quality. Thus, we add valid inequalities to (ILP)
as discussed at the end of the preceding section. We remove the variable ﬁxing
again, and continue with a standard branch-and-bound on model (ILP). We
summarize the whole procedure, we term (FIX), in Algorithm 1. Note thatthe variable ﬁxing is a heuristic, even though for all our instances we quickly
ﬁnd an integer solution. In one sentence, it is a very careful initialization of
branch-and-bound.
Our richness of models allows for several alternatives. Instead of using(ILP) we could run (BLP), or combinations of the models. Also, the ﬁnal
model to be solved need not be a linearization. We could feed our integral
solution obtained from the ﬁxed linearization into model (NLP) as a better
starting point. The subsequent nonlinear programming based branch-andbound phase performs as is well-known for linear mixed integer programs,

A fast algorithm for near cost optimal line plans

215

except that ‘‘bounds’’ in each node are obtained from the relaxed nonlinear
integer program (NLP).
7 Computational results
We implemented all our models in the GAMS modeling language [1], version
20.6, using the combination of CONOPT2 and SBB for solving the nonlinear
models, and CPLEX 7.5 for the solution of the linear (integer) models. Our
code is publicly available.3 We deem it important that the reader be in aposition to get ﬁrst-hand experience with our techniques, reproduce ourresults, and verify our conclusions. In the operations research community,
usually with reference to proprietary data, these key ingredients of scientiﬁc
activity are often regarded as being of secondary importance. In contrast, we
would like to enable the reader to check claims which we do not explicitly
support by numerical results.
Our computational experiments are performed on a 700MHz Linux PC
with an execution time limit of three hours. We are provided with four practical problem instances from the Dutch Railways (Nederlandse Spoorwegen),ic97, ic98, ir98, and ar98. The names stem from the underlying network(InterCity, InterRegio, and AggloRegio) and the year. Variants of these
instances have been circulating in the community for years and, except ir98,
so far withstood a proven optimal solution, even with massive computational
power [2]. Only recently, Bixby optimally solved instance ic97 in a month’s
computation time using our model (CUT) (private communication), andweobtained an optimal solution to ic98 using the same model andcomparable computational eﬀorts. Since our comparison is intended to bequalitative in nature we do not present numerical results but illustrate ourexperiments in Figures 2 through 5. Horizontal solid lines indicate the optimal
objective function value, or for ar98, the best known upper and lower
bounds, respectively. All bounds are obtained using model (CUT).
In brief, as expected, there is no consistently best all purpose model.
However, when looking at upper and lower bounds separately, the picturechanges. Most strikingly, for three out of four instances our algorithm (FIX)
provides us with our by far best feasible solutions in about ﬁve minutes
computation time per instance. There is practically no improvement after
three hours of additional search with either model, (ILP) or (BLP). A bit
surprisingly, the integer solutions found by model (ILP) alone do not beneﬁt
from the smaller model size and the related quicker evaluation of branch-andbound nodes compared to (BLP). In what regards lower bounds, model(CUT) is signiﬁcantly superior except in one case; we give a cut statistic in
Table 1. (ILP) appears to be next best, again with one exception. Interestingly
enough, this one exception (with (FIX) unsuccessful, (CUT) and (ILP) inferior to (BLP), and the instance being optimally solvable within practical timebounds) is one and the same instance, namely ir98.
Results for model (ILP-c), both in terms of upper and lower bounds arepoor, and we refrain from plotting them. The structure of this seemingly
intuitive model is too heavily disturbed by the necessity to model the3

URL http://www.gamsworld.org/minlp/apps/blllop

216

Michael R. Bussieck et al.

Algorithm 1 Combining our models into the variable ﬁxing procedure (FIX)
Solve the LP relaxation of (ILP)
Let xILP ; y ILP denote an optimal fractional solution
Initialize variables xNLP ; y NLP of model (NLP)
for all r 2 R X
do¼
u  xILPxNLP
rr;u for all r 2 R
u2F
(P
ILP
u  yr;u=xNLP
if xNLP>0rr
yrNLP ¼ u2F
0otherwiseend forSolve the continuous relaxation of (NLP)
According to the optimal continuous solution, ﬁx variables in (ILP):
forall r 2 R, u 2 F do
 threshold (default: 105 ) thenif xNLPr
Add constraint xr;u ¼ 0 to (ILP)
end ifend forIntegrally solve the ﬁxed (ILP)
Drop again constraints xr;u ¼ 0 from (ILP)
Sequentially add valid inequalities:
Add (5) to (ILP)
if (4) violated thenAdd (4) to (ILP)
else if (6) violated thenAdd (6) to (ILP)
else if (7) violated thenAdd (7) to (ILP)
end ifSolve LP relaxation of (ILP), and remove all cuts which are not binding in the optimal solution
Warm start (ILP) with the integral solution foundIntegrally solve (ILP) by standard branch-and-bounddiscontinuous term in the objective function via (2). Note that these ﬁndings
stress the advantage of the other two linearizations. It should be mentionedthat Goossens et al. optimally solve a variant of instance ir98 within threeminutes of computation time using (BLP) as the basic model [8]. The eﬀectiveness of their branch-and-cut approach for (BLP) leads to the question
how well this algorithm performed with (ILP) as basic model. However, as
indicated in [8] the data instances these authors use are slightly diﬀerent fromours and serious comparisons cannot be drawn. From their experiments it
appears that problem speciﬁc branching rules, like branching on lines orcapacities are not signiﬁcantly superior over standard variable branching.
We made several experiments in order to identify the most eﬀective
combination of our models for the variable ﬁxing algorithm, cf. Table 2.
Admittedly, the beneﬁt of model (NLP) is marginal in terms of additional
solution quality. Almost equally good solutions are obtained by just solvinge.g., the relaxed (ILP) model, eliminate unpromising lines as above, and solve
the resulting smaller integer program. However, in our opinion, model (NLP)
most naturally reﬂects the practical background via (1) which may lead to aslightly better ‘‘guidance’’ of the variable ﬁxing. We consider it very important that – with or without (NLP) – we always obtain integer solutions of highquality; this points to the robustness of both approaches, and we value theadded ﬂexibility of using several models.

A fast algorithm for near cost optimal line plans

217

Table 1. Cut statistic and impact of valid inequalities on the linear programming relaxation ofmodel (ILP). The columns indicate the name of the instance, the respective number of violated
inequalities, and the respective objective function values of the LP relaxation, the LP relaxation
with violated inequalities added, the LP value in the root node of the branch-and-bound tree
(after CPLEX’ preprocessing), and the latter with violated inequalities added, in that orderInstance (4)

(5)

(6)

(7)

plain LPw/(4)–(7)(ILP) rootw/(4)–(7)

15
19
20
72

4
1
2
15

108
134
47
162

24
5
0
16

3845.474328.83
2230.545882.903915.544412.832291.076009.103875.454376.06
2286.286064.06

3922.634418.152299.176094.40

ic97
ic98
ir98
ar98

Table 2. Comparison of diﬀerent model combinations for our variable ﬁxing strategy. Headings‘init’, ‘guide’, and ‘ﬁxed’, respectively, refer to the relaxation which is used to ﬁnd an initialsolution for model (NLP), the relaxation for determining which lines to discard, and the partially
ﬁxed model to ﬁnd an integer feasible solution. We report the objective function value of the best
integer solution found after three minutes of the ‘ﬁxed’ run. A * indicates an optimal solution ofthe ‘ﬁxed’ model
init

guide

Fixed

ic97

ic98

ir98

ar98

(ILP)
–
(BLP)
–

(NLP)
(ILP)
(NLP)
(BLP)

(ILP)
(ILP)
(BLP)
(BLP)

*4088.52*4088.52
4371.664476.51*4521.674559.90
4811.934639.63*2385.84*2417.97*2378.33*2386.27

6336.446368.086326.64*6383.59

Fig. 2. Development of bounds for instance ic97

218

Michael R. Bussieck et al.

Fig. 3. Development of bounds for instance ic98

8 Conclusions
Nonlinear integer programming is quite often an adequate and intuitivemodeling approach for combinatorial optimization problems. Recently,
powerful software became available for actually solving these models, and wehave to re-think about such approaches. In this paper we support this claim,and demonstrate how to combine nonlinear techniques with ‘‘classical’’integer linear programming in order to successfully attack a very hardpractical combinatorial optimization problem. The heuristic we present isrobust in that it delivers an excellent integer solution for all our instances.
For better or for worse, the success or failure of an approach is linked tothe dramatic development of commercial solvers. Not seldom we see that we
can learn a lot from a tailored algorithm but an oﬀ-the-shelf product is muchsuperior in terms of computation time, and often quality. This is especiallytrue for linear and mixed integer programs. Solvers are much more sophisticated than those available when we set out with this research in 1997; in fact,
we experimented, among others, with CPLEX versions 2.0 through 8.0 and it
is not unlikely that future versions are competitive to our success with thevariable ﬁxing algorithm, or put its current usefulness in question altogether.
However, we would like to make the point that these solvers also improve
because of research like ours. Then again, in our opinion, the development ofindustry standard mixed integer nonlinear solvers is just at its beginning, andthe headway of using a nonlinear model may even be ampliﬁed in the future.
One must critically ask the question whether an intelligent use of diﬀerent
models and solvers is eligible to be termed an algorithm, even though thenotion is certainly technically appropriate. Our answer is aﬃrmative. The

A fast algorithm for near cost optimal line plans

Fig. 4. Development of bounds for instance ir98

Fig. 5. Development of bounds for instance ar98

219

220

Michael R. Bussieck et al.

outcome in terms of the building blocks model and solver looks simple but
beneath the surface a lot of algorithmic understanding is indispensable. On
the other hand, this simplicity at a higher level enables the practitioner toobtain reproducible and predictable results in ﬁnite time. We would like thenotion of a practical algorithm.
Besides fundamental investigations in the latter direction open research
avenues include the following speciﬁc modiﬁcations: Lines need not be
operated back and forth. This may lead to concatenations of several lines. An
additional complication in this setting is whether the resulting line plan isrobust as to the subsequent planning stage of train scheduling. Also, rollingstock does not need to be dedicated to speciﬁc lines. Column generation, an
entirely diﬀerent algorithmic approach, has been suggested for the maximization of the number of direct travelers in [2]. A continuation of our study
would investigate how suitable this approach is for the minimization ofoperational cost. We are conﬁdent that further computational and methodological progress will eventually lead to an integrated treatment of all thementioned networks.
References
1. Brook A, Kendrick D, Meeraus A (1992) GAMS: A user’s guide, Release 2.25. ScientiﬁcPress, San Francisco, CA
2. Bussieck MR (1998) Optimal lines in public rail transport. PhD thesis, Braunschweig
University of Technology, Available under http://www.biblio.tu-bs.de/ediss/data=19990121a=19990121a:pdf
3. Bussieck MR, Kreuzer P, Zimmermann UT (1997) Optimal lines for railway systems.
European J. Oper. Res. 96(1):54–63
4. Bussieck MR, Winter T, Zimmermann UT (1997) Discrete optimization in public railtransport. Math. Programming 79(1–3):415–444
5. Claessens MT (1994) De kost-lijnvoering. Master’s thesis, University of Amsterdam
6. Claessens MT, van Dijk NM, Zwaneveld PJ (1998) Cost optimal allocation of rail passenger
lines. European J. Oper. Res. 110(3):474–489
7. Fischetti M, Lodi A (2003) Local branching. Math. Programming 98(1–3):23–47
8. Goossens J-W, van Hoesel S, Kroon L (2001) A branch-and-cut approach for solving line
planning problems. METEOR Research Memorandum RM/01/016, University of Maastricht. Accepted to be published in Transportation Sci.
9. Grossman TA (2001) Reversing tradition. OR/MS Today 28(4):22–25
10. Kocis GR, Grossmann IE (1989) Computational experience with DICOPT solving MINLP
problems in process systems engineering. Computers and Chemical Engineering 13:307–315
11. Zwaneveld PJ (1997) Railway Planning—Routing of Trains and Allocation of Passenger
Lines. PhD thesis, Rotterdam School of Management