SIAM REVIEW
c
2003 Society for Industrial and Applied Mathematics
Vol. 45, No. 1, pp. 116–123
Teaching Integer Programming
Formulations Using the
Traveling Salesman Problem
∗
G
´
abor Pataki
†
Abstract. We designed a simple computational exercise to compare weak and strong integer pro-
gramming formulations of the traveling salesman problem. Using commercial IP software,
and a short (60 line long) MATLAB code, students can optimally solve instances with up
to 70 cities in a few minutes by adding cuts from the stronger formulation to the weaker,
but simpler one.
Key words. integer programming, traveling salesman problem, subtour elimination constraints, cut-
ting planes
AMS subject classifications. 90C10, 90C11, 90C27, 90C57, 90-01, 97D40
PII. S0036144502368551
1. Introduction.
1.1. Integer Programs and Their Formulations. An integer program (IP) is an
optimization problem of the form
min c
T
x
s.t. x ∈ X,
where X = Z
n
∩{x ∈ R
n
|Ax ≤ b }, with some matrix A, vector b, and the symbols
Z
n
and R
n
denoting the set of integer and real n-vectors, respectively. The polyhe-
dron {x ∈ R
n
|Ax ≤ b } is a formulation of the set X. An IP has many formulations;
for instance, in Figure 1 the solid and dashed lines enclose the same integer points;
i.e., the corresponding polyhedra are both formulations of the same set.
Given two such polyhedra, P and Q (i.e., P ∩ Z
n
= Q ∩ Z
n
), we say that P is a
better,orstronger, formulation, if P Q. (The definitions here are from [9].)
One of the most important skills that a practitioner of integer programming must
acquire is that of designing a strong formulation for a particular problem. The main
component of all commercial IP solvers is a branch-and-bound algorithm that uses
the linear programming (LP) relaxation of an IP (i.e., the problem obtained from an
IP by discarding the integrality constraints). Hence, a stronger formulation usually
can be solved with fewer branch-and-bound nodes; if the tighter LPs are not “much”
more time-consuming, this translates into a smaller overall solution time. Even if the
∗
Received by the editors May 2, 2001; accepted for publication (in revised form) May 28, 2002;
published electronically February 3, 2003. Part of this work was done while the author was affiliated
with the Department of IE/OR at Columbia University, and was supported by NSF grant DMS
95-27124.
http://www.siam.org/journals/sirev/45-1/36855.html
†
Department of Operations Research, University of North Carolina at Chapel Hill, 212 Smith
Building, Chapel Hill, NC 27599-3180 (gabor@unc.edu).
116
Downloaded 02/05/21 to 136.56.14.111. Redistribution subject to SIAM license or copyright; see https://epubs.siam.org/page/terms
TEACHING INTEGER PROGRAMMING USING THE TSP 117
Fig. 1 Two formulations of the same set.
problem cannot be solved to optimality within the available time, a strong formulation
provides a good bound on the optimal value of the problem. Hence it can also serve
as a counterpoint to an effective heuristic, by proving that a solution provided by the
latter is close enough to being optimal.
1.2. Exercises to Compare Formulations in Practice. A classroom lecture on
comparing weak and strong formulations is best accompanied by an assignment asking
students to do a computational comparison. Such a comparison can simply consist
of feeding a strong and a weak formulation of the same problem to an IP solver and
comparing the number of branch-and-bound nodes and times required to solve them
to optimality. A problem is well suited for the purpose of a comparison if
(1) it is relevant in practice, interesting, and easy to understand;
(2) the advantages of the strong formulation are not immediately apparent, for
some of the following reasons:
• the strong formulation uses many more variables and/or constraints;
• the weak one can accommodate more versatile objective functions;
(3) the weak and strong formulations are easy to generate.
Several such problems are (see, e.g., [9])
• facility location problems (with their weak, i.e., the aggregated, and the
strong, i.e., the disaggregated formulation);
• knapsack problems (their usual formulation can easily be strengthened by
cover inequalities);
• lot-sizing problems (their usual formulation can be strengthened by various
inequalities).
For many problems, however, the advances in IP software disguise the advantage of
providing a stronger formulation to solvers. Most IP solvers now incorporate auto-
matic reformulation techniques that can substantially strengthen a weak formulation.
Such techniques include disaggregation (i.e., replacing the constraint
m
i=1
x
i
≤ my
on the 0–1 variables x
i
and y with the inequalities x
i
≤ y (i =1,...,m)), generating
many of the inequalities that one would add to the knapsack, and lot-sizing problems
(e.g., covers, and flow-covers). Hence, in the case of the models listed above,
• frequently there is no significant difference in the solution times, when feeding
a weak or a strong formulation to the solver;
Downloaded 02/05/21 to 136.56.14.111. Redistribution subject to SIAM license or copyright; see https://epubs.siam.org/page/terms
118 G
´
ABOR PATAKI
• sometimes the solution times are even better for the weaker formulation. The
reason is that IP solvers strengthen the weak formulation quite intelligently.
For instance, when disaggregating, they only add those x
i
≤ y inequalities,
which are violated by the current LP solution, thus not enlarging the LP
relaxation unnecessarily;
• or, the strong one must be quite sophisticated and hence not too useful at
the beginning of a course, when the importance of strong formulations should
already “sink in.” Such an example is strengthening the lot-sizing problem
with path-inequalities ([9, p. 223]).
Of course, one could always ask the students to turn off the preprocessor, or the
cut-generator of the solver; this would hardly make a convincing exercise, though.
A problem that well fits criteria (1) and (2) is the traveling salesman problem
(TSP): given n cities, and pairwise travel costs between them, find the shortest
tour, i.e., a directed cycle containing all cities. (In fact, students invariably find
the TSP more fascinating than the somewhat mundane knapsack, location, or lot-
sizing problems.) The formulation that serves as a benchmark for all others is the
subtour formulation; there are also many weaker ones. Three articles that survey
and analytically compare some of them are [1], [5], and [8]. The weaker formulations
all use extra variables; i.e., the set of tours X is written as X = Z
n
∩ P , where
P = {x |Bx + Du ≤ f for some u ∈ R
m
}.NowP is the projection of the poly-
hedron {(x, u) |Bx + Du ≤ f }, hence it is also a polyhedron. Thus we can write
P = {x |Gx ≤ h } for some matrix G and vector h. The latter representation, how-
ever, may contain many more and/or more complicated inequalities than the first.
Designing an experiment that allows students to do a computational comparison
is not quite obvious though. This note gives a simple guide to comparing the subtour
formulation to a weak one due to Miller, Tucker, and Zemlin [6] on relatively small, but
far from trivial instances. It utilizes a primitive “cutting plane algorithm” of about
60 lines of MATLAB code that, used with a commercial IP solver, allows students to
do the comparison and solve real world instances with up to 70 nodes to optimality.
2. The Formulations.
2.1. Problem Statement and the Node Constraints. In the complete directed
graph D =(N,A), with arc-costs c
ij
, we seek the tour (a directed cycle that contains
all n cities) of minimal length. Define the variables
x
ij
=
1 if arc (i, j) is in the tour,
0 otherwise.
The feasible solutions of the integer program (in fact, there is always an integer
optimal solution to the LP relaxation)
(2.1)
min
i,j
c
ij
x
ij
s.t.
i
x
ij
=1 ∀i,
j
x
ji
=1 ∀i,
0 ≤ x
ij
≤ 1,x
ij
integer,
may contain several directed cycles, called subtours. The constraints of (2.1) are
called the assignment constraints.
Downloaded 02/05/21 to 136.56.14.111. Redistribution subject to SIAM license or copyright; see https://epubs.siam.org/page/terms
TEACHING INTEGER PROGRAMMING USING THE TSP 119
2.2. The MTZ Formulation. To exclude subtours, one can use extra variables
u
i
(i =1,...,n), and the constraints
(2.2)
u
1
=1,
2 ≤ u
i
≤ n ∀i =1,
u
i
− u
j
+1≤ (n −1)(1 −x
ij
) ∀i =1, ∀j =1.
We call the last inequality in (2.2) an arc-constraint. The formulation consisting of
(2.1) and (2.2) is called the Miller–Tucker–Zemlin (MTZ) formulation of the TSP; see
[6]. It indeed excludes subtours, as (1) the arc-constraint for (i, j) forces u
j
≥ u
i
+1,
when x
ij
= 1; (2) if a feasible solution of (2.1)–(2.2) contained more than one subtour,
then at least one of these would not contain node 1, and along this subtour the u
i
values would have to increase to infinity. This argument, with the bounds on the u
i
variables, also implies that the only feasible value of u
i
is the position of node i in the
tour. The advantages of the MTZ formulation are
• its small size (we need only n extra variables and roughly n
2
/2 extra con-
straints),
• if it is preferable to visit, say, city i early in the tour, one can easily model
this by adding a term −αu
i
with some α>0 to the objective.
2.3. The Subtour Formulation. The other way to exclude subtours is by adding
to (2.1) the family of subtour (or subtour elimination) constraints
i∈S,j∈S
x
ij
≤|S|−1(S V, |S| > 1),(2.3)
to obtain the subtour formulation of the TSP consisting of (2.1) and (2.3). (As the
subtour inequality for N \ S is a linear combination of the inequality for S and of
the assignment constraints, it is enough to use the subtour inequalities with S having
size at most n/2.) It does not have the advantages of the MTZ formulation. The
disadvantage of its exponential size is mitigated by the fact that not all subtour in-
equalities must be put into the formulation from the start. They can be generated
as needed by a separation algorithm: one can start with the formulation (2.1), then
generate subtour inequalities that are violated by the current LP solution. The sep-
aration algorithm for subtour constraints is based on network flow techniques; for
further details, see [4].
3. The Comparison.
3.1. The Strength of the Two Formulations. Solving a reasonably large (with
at least, say, 50 cities) problem to optimality is only possible using the subtour for-
mulation; at least, we are not aware of any published computational studies that use
the pure MTZ formulation. There is an analytical explanation of this fact: the LP
relaxation of the MTZ formulation is much weaker. Since it contains the extra u
i
variables, for the comparison we need to eliminate them by taking appropriate linear
combinations of the inequalities in (2.1)–(2.2). One way of doing this is by summing
the arc-inequalities along a directed cycle. If the arc set of the cycle is denoted by C,
then the result is
(i,j)∈C
x
ij
≤
1 −
1
n − 1
|C|.(3.1)
On the other hand, if we denote the node set of the cycle by N(C), then the subtour
Downloaded 02/05/21 to 136.56.14.111. Redistribution subject to SIAM license or copyright; see https://epubs.siam.org/page/terms
120 G
´
ABOR PATAKI
formulation contains the obviously stronger inequality
i,j∈N(C)
x
ij
≤|C|−1.(3.2)
The intuitive fact that adding the arc-inequalities along a cycle is the only “essential”
way to eliminate the u
i
variables can be made precise. Balas [1] and Padberg and
Sung [8] showed that all nondominated inequalities in the projection of the MTZ
formulation onto the space of the x variables are of the form (3.1).
3.2. How to Combine Them?. An efficient implementation of a separation rou-
tine for subtour inequalities is beyond the ability of most practitioners and under-
graduate/MS students. A collection of excellent public-domain routines for the TSP,
called Concorde, is available, with such routines among them; see [3]. However, most
problems encountered in practice are not pure TSPs, only “TSP-like.” One example
is the vehicle routing problem, in which there are k “salesmen” (i.e., trucks), and each
of the cities must be visited by one to make a given amount of delivery, while obeying
capacity restrictions. For TSP-like problems the TSP formulations are usually not
hard to generalize, but the TSP separation routines cannot be used.
One can combine the MTZ and subtour formulations to obtain the ease of use of
the first and some of the strength of the second. Violated subtour inequalities are easy
to identify if a solution satisfies the assignment constraints and is 0–1: one can take
S to be the union of several subtours. Such an approach is likely to work for TSP-like
problems as well. For small problems, one can identify such subtour inequalities even
by inspection: see, for instance, the article [7] for such an approach to a variant of the
TSP arising in aircraft mission planning. We used the following “cutting plane algo-
rithm” to solve several TSP instances to optimality, with maxrounds an integer be-
tween 0 and 2, depending on how much we wanted to strengthen the MTZ formulation.
1. Let the IP formulation consist of the assignment constraints only, and k =1.
2. while k ≤ maxrounds
2a. Solve the IP over the current formulation. Assume that the optimal
solution consists of r subtours S
1
,...,S
r
.
2b. If r = 1, STOP; the solution is optimal to the TSP. Otherwise, add
to the formulation at most 1000 subtour constraints, in which S is the
union of several S
i
sets and |S|≤n/2. Set k = k +1.
end while
3. Add the arc constraints to the formulation, and solve the IP to optimality.
We used a simple implementation of step 2b: we output the (i, j) pairs corresponding
to the nonzero x
ij
values into a file, then ran a MATLAB routine to generate the
required constraints. These were then pasted into the formulation. For example,
suppose that the number of cities is 10, and the arcs corresponding to the variables
at 1 are as shown in Figure 2. The nontrivial parts (i.e., excluding input, output, and
obvious initializations) of the MATLAB code are given in Figure 3. The code works
as follows:
• The first part constructs the list of subtours from the list of arcs. In the
example, it will create the list {1, 4, 6}, {2, 3, 8}, {5, 7}, {9, 10}.
• The second part constructs all appropriate subsets of nodes in the variable S,
and the index-pairs corresponding to arcs within S, in the variable arcs
in S.
The code works with all IP solvers that have a “reasonable” user interface: all
we need to do is output the index-pairs corresponding to the variables at 1 and paste
Downloaded 02/05/21 to 136.56.14.111. Redistribution subject to SIAM license or copyright; see https://epubs.siam.org/page/terms