Integer programming formulations

Read first

• Post 1: What is the TSP? A tour visits every stop once and returns to the start, and its cost is the sum of its legs, read from the distance matrix d(i, j). Our running example is eight stops, A to H, whose optimal tour A-B-C-D-H-E-F-G-A costs 315 (in units of 100 m).
• Primer C: linear and integer programming A linear program (LP) has variables, linear constraints and a linear objective. An integer program (IP) also requires some variables to be whole numbers. Dropping that requirement gives the LP relaxation, whose optimum is a lower bound on the IP’s optimum, because it searches a bigger set.

The question this post answers: how do you write the TSP as an integer program, and how close does its LP relaxation get to the true optimum?

Intuition

In Post 5 I built lower bounds by hand from structures that only make sense for the TSP: 1-trees, assignments, Lagrangian penalties. Integer programming is the general-purpose alternative. You describe the problem once as variables and linear constraints, and the whole machinery of LP bounds, branching and cutting applies to it without any TSP-specific code. Most exact routing methods today, including the strongest TSP code, start from such a formulation. The question is which one, because the TSP can be written in very different ways, and they are not equally good.

Variables. For every ordered pair of stops, one yes/no decision: x(i, j) = 1 if the tour drives straight from i to j, and 0 otherwise. Eight stops give 8 × 7 = 56 variables. The cost of a solution is the sum of d(i, j) over the legs that are switched on.

Degree constraints. Every stop is left exactly once and entered exactly once. That is two constraints per stop, 16 in all. For stop B, for example: x(B, A) + x(B, C) + … + x(B, H) = 1 for the way out, and x(A, B) + x(C, B) + … + x(H, B) = 1 for the way in.

Is that enough? Let me solve the model with just these constraints. A 0/1 solution that leaves and enters every stop once is simply a successor for every stop, so I can list all of them (14,833 for eight stops) and take the cheapest. In the matrix below, a solution is one outlined cell in every row (where each stop goes next) and one in every column (where each stop is entered from).

Every row and every column has exactly one outlined cell, so the degree constraints are satisfied. But each stop simply drives to its nearest partner and straight back. The table below lists the four pairs.

These four little loops are subtours: closed cycles that do not visit every stop. They cost 282, which is less than the optimal tour, because they are not a tour at all. This model is the assignment problem, and its optimum is the assignment bound from Post 5.

Subtour elimination. To forbid the loop B→C→B, say that inside the set {B, C} at most one leg may be used: x(B, C) + x(C, B) ≤ 1. In general, for any set S of stops, at most |S| − 1 legs may run between stops of S. A real tour always passes this test: inside S its legs form paths, never a closed cycle. A subtour that stays inside S fails it, because it uses |S| legs there.

I added the four constraints for the four pairs and solved again. The table below shows both rounds.

Round 2 is a single tour, and it is the optimal one from Post 1. Step through the diagram below to see both rounds on the map.

The diagram draws unrounded kilometres, so its totals differ by a few tenths from the matrix costs in the tables.

Notice what made this easy: I saw the subtours, so I knew which constraints to add. Four were enough. But the formulation itself does not know in advance which subtours a solver will stumble into, so it has to list a constraint for every set S that could hold one. For eight stops that is 238 constraints, and the count doubles with every extra stop.

The formal idea

Number the stops 1 to n, with stop 1 as the depot. For every ordered pair (i, j) with i ≠ j, let d(i, j) be the cost of the leg and x(i, j) ∈ {0, 1} the decision to drive it. The Dantzig-Fulkerson-Johnson (DFJ) formulation comes from the 1954 paper that solved a 49-city tour of the United States: one city in each of the 48 states, plus Washington, D.C. They actually solved a 42-city version with seven East Coast cities removed. Its optimal tour used the leg from Washington to Boston, and the shortest road between those two runs through the seven missing cities, so it also solved the 49-city problem. The formulation is:

minimize    Σ  d(i,j) · x(i,j)                over all i ≠ j

subject to  Σ  x(i,j) = 1                     for every stop i      (leave i once)
            j

            Σ  x(i,j) = 1                     for every stop j      (enter j once)
            i

            Σ    Σ   x(i,j) ≤ |S| − 1          for every set S with 2 ≤ |S| ≤ n − 2
           i∈S  j∈S                                                  (no subtour inside S)

            x(i,j) ∈ {0, 1}

The sizes follow directly. There are n(n − 1) variables and 2n degree constraints. The subtour constraints run over all subsets of size 2 to n − 2, which is 2ⁿ − 2n − 2 of them. Sets of size 1 and n − 1 are left out because the degree constraints already take care of them. A set and its complement give the same restriction, so in practice half of them can go, but half of an exponential is still exponential.

The diagram below shows how the pieces fit together, and the two roads from the assignment problem to a correct formulation.

flowchart TD
  X["One 0/1 variable per leg"] --> DEG["Degree constraints: leave once, enter once"]
  DEG --> AP["Assignment problem: subtours allowed"]
  AP --> Q{"How to forbid subtours?"}
  Q --> DFJ["DFJ: fewer legs than stops inside every set S"]
  Q --> MTZ["MTZ: order variables u, one row per pair of stops"]
  DFJ --> D2["2ⁿ − 2n − 2 rows, strong LP bound"]
  MTZ --> M2["(n − 1)(n − 2) rows, weak LP bound"]

MTZ: a compact alternative. Miller, Tucker and Zemlin (1960) avoided the exponential count with extra variables. Give every stop except the depot a number u(i) between 1 and n − 1, its position in the tour. Then, for every pair of non-depot stops i ≠ j, add:

u(i) − u(j) + (n − 1) · x(i,j) ≤ n − 2

If x(i, j) = 1, this says u(j) ≥ u(i) + 1: positions increase along every leg. If x(i, j) = 0, it says u(i) − u(j) ≤ n − 2, which is always true for numbers between 1 and n − 1, so it does nothing. Now take any subtour that avoids the depot and walk around it: every leg raises u by at least 1, and you arrive back where you started with a larger u than you began with. That is impossible, so the subtour is infeasible. The real tour escapes the contradiction because it passes through the depot, which has no u. The diagram below runs the argument on the loop B→C→B from round 1.

flowchart LR
  B["B, position u(B)"] -->|"leg B→C forces u(C) ≥ u(B) + 1"| C["C, position u(C)"]
  C -->|"leg C→B forces u(B) ≥ u(C) + 1"| B
  C --> X["Together: u(B) ≥ u(B) + 2, impossible"]

MTZ needs only (n − 1)(n − 2) of these rows and n − 1 extra variables. For berlin52 that is 2,550 rows, against 4.5 × 10¹⁵ for DFJ. The table below compares the two sizes.

LP relaxation and the gap. Both formulations describe exactly the same set of tours when x must be 0 or 1. They differ once you relax x to any value between 0 and 1. The LP relaxation can be solved fast, and its optimum is a lower bound on the optimal tour. Its gap is how far below the optimum it lands: (optimum − LP value) / optimum. A small gap means the branch-and-bound search of Post 5 prunes early; a large gap means it explores a huge tree.

I computed these LP relaxations with an LP solver (SciPy’s interface to HiGHS); the script at the end of the post does not need one. On our eight stops:

And on berlin52, where the optimum is 7,542:

MTZ barely improves on the assignment bound. The reason is in its constraint: with a fractional x(i, j) = 0.5 and n = 52, it only demands u(j) ≥ u(i) − 24.5, which the positions can satisfy without any real ordering. DFJ’s LP, on the other hand, lands exactly on the optimum here. That is not a coincidence of formulation: Held and Karp (1970) showed that this LP bound equals the Lagrangian 1-tree bound from Post 5, so the two posts arrive at the same number from different directions. On these instances the bound is perfect. On harder ones it can sit below the optimum, but it is still the bound that every serious TSP code builds on.

The symmetric version. For a symmetric instance you can halve the variables: one variable per edge {i, j} instead of one per direction, and “exactly two edges at every stop” instead of “leave once, enter once”. Two-stop loops then become impossible, and on our eight stops the degree constraints alone already return the optimal tour. On berlin52 they do not: the degree-only LP gives 7,163 and breaks into four subtours of 33, 10, 6 and 3 cities. Adding the four subtour constraints for exactly those sets lifts the LP to 7,542, and the solution is a single tour.

Common misconceptions

“Leave every stop once and enter every stop once, and you have a tour.” You have an assignment, which can be any collection of closed loops. On our eight stops the cheapest one is four two-stop loops costing 282. Without subtour elimination the model is a different problem.

“MTZ is better because it is smaller.” It is smaller: 2,550 rows against 4.5 × 10¹⁵ for berlin52. But its LP relaxation is 16.2% below the optimum, almost as weak as having no subtour constraints at all, so a branch-and-bound search on it has very little to prune with. Solvers spend their time where the bound is weak, and MTZ’s is weak. Strengthened versions exist, such as the lifted constraints of Desrochers and Laporte (1991). They narrow the gap but stay well short of DFJ. On ten asymmetric TSPLIB instances, Öncan, Altınel and Laporte (2009) report an average LP gap of 7.4% for MTZ, 4.7% for the lifted version, and 1.1% for a formulation whose LP equals DFJ’s.

“The LP relaxation gives you a tour.” It gives you a number and a solution that may use legs at fractional values like 0.5, or split into subtours. The number is a lower bound. Only when the LP solution happens to be a whole-number single tour, as with DFJ on berlin52, is it also the answer.

“To use DFJ, you must write down all of its constraints.” On our stops, 4 of the 238 subtour constraints were enough, and on berlin52, 4 of more than 4.5 × 10¹⁵. You only need the ones a candidate solution actually violates. The hard part is finding those automatically, for solutions that are fractional, without listing the rest. That is the next post.

Try it in Python

This script needs no solver. It prints the size of both formulations, then solves the eight-stop model exactly the way I did by hand above: list every solution of the degree constraints, take the cheapest, and if it has subtours, add a subtour constraint for each one and repeat. It finishes by checking that the optimal tour, with u = position in the tour, satisfies every MTZ row.

import math
from itertools import permutations

STOPS = {"A": (2, 2), "B": (6, 1), "C": (9, 3), "D": (8, 7),
         "E": (5, 9), "F": (1, 8), "G": (4, 5), "H": (10, 9)}  # km; A is the depot
N = list(STOPS)

def d(a, b):  # straight-line distance in units of 100 m, rounded
    (x1, y1), (x2, y2) = STOPS[a], STOPS[b]
    return int(10 * math.hypot(x1 - x2, y1 - y2) + 0.5)

# 1. Size of the two formulations (directed: one variable x_ij per ordered pair)
print("n   variables  degree rows  DFJ subtour rows       MTZ rows")
for n in (8, 15, 52):
    print(f"{n:<3} {n * (n - 1):>9,}  {2 * n:>11}  {2**n - 2 * n - 2:>21,}  {(n - 1) * (n - 2):>8,}")

# 2. Every 0/1 solution of the degree constraints is a successor map: succ[i] = the j with x_ij = 1.
#    So "solve the IP with degree constraints only" = try every permutation with no fixed point.
def cycles(succ):
    out, seen = [], set()
    for s in N:
        if s not in seen:
            c, x = [s], succ[s]
            seen.add(s)
            while x != s:
                c.append(x)
                seen.add(x)
                x = succ[x]
            out.append("".join(c))
    return out

def satisfies(succ, cuts):  # DFJ: inside any subset S, at most |S| - 1 arcs
    return all(sum(succ[i] in S for i in S) <= len(S) - 1 for S in cuts)

solutions = sorted(
    ((sum(d(i, succ[i]) for i in N), succ)
     for succ in (dict(zip(N, p)) for p in permutations(N))
     if all(succ[i] != i for i in N)),
    key=lambda t: t[0])

cuts = []
for rnd in range(1, 10):
    cost, succ = next((c, s) for c, s in solutions if satisfies(s, cuts))
    subtours = cycles(succ)
    print(f"round {rnd}: cost {cost}, cycles {subtours}")
    if len(subtours) == 1:
        break
    cuts += [set(c) for c in subtours]  # forbid exactly the subtours we just saw
print(f"subtour constraints used: {len(cuts)} of {2**len(N) - 2 * len(N) - 2}")

# 3. MTZ: with u = position in the tour, every constraint u_i - u_j + (n-1) x_ij <= n-2 holds
tour = "ABCDHEFG"
u = {s: k for k, s in enumerate(tour)}  # A gets 0 and is left out of the MTZ rows
n = len(N)
ok = all(u[i] - u[j] + (n - 1) * (tour[(tour.index(i) + 1) % n] == j) <= n - 2
         for i in N[1:] for j in N[1:] if i != j)
print("optimal tour satisfies all MTZ rows:", ok)

It prints the sizes, then the two rounds:

Two things to notice. The loop in part 2 is the whole idea of DFJ in a few lines: subtour constraints only matter when a solution breaks them, and the four that round 1 broke were all it took. And part 2 only works because eight stops have 14,833 successor maps; an IP solver gets the same answers without listing them, which is what the LP relaxation and branch and bound are for.

Exercise

Take the six stops from Posts 4 and 5 (A, B, C, D, G and H). Post 5 found that the cheapest solution of the degree constraints is three two-stop loops, A⇄G, B⇄C and D⇄H, costing 200.

1. Write the three DFJ subtour constraints that forbid those loops.
2. Check that the optimal tour from Posts 4 and 5, A-B-C-H-D-G-A, satisfies all three.
3. How many DFJ subtour rows and how many MTZ rows does the full six-stop model have?
4. Add your three constraints and find the cheapest solution again. Is it a tour?

Hint: for step 4, delete E and F from STOPS in the script, set tour = "ABCHDG", and run it. For step 3, use the formulas in the formal idea with n = 6.

What this method can’t do yet

The DFJ formulation is the right one: on both of our instances its LP relaxation lands exactly on the optimum. But it is written over every subset of stops, and that is 4.5 × 10¹⁵ constraints for berlin52, far more than any solver can be handed. MTZ fits in 2,550 rows, but its LP bound is barely better than no subtour constraints at all. On the eight stops I got around this by looking at the subtours and adding only the four constraints they broke. A solver needs to do that automatically, for LP solutions that are fractional rather than neat little loops, and it needs to combine it with branching when the LP alone does not give a tour. Post 7, “Branch-and-cut and Concorde”, shows how: add violated constraints only when needed, and how Concorde uses that idea to solve very large TSPs exactly.

References

• Dantzig, G., Fulkerson, R., and Johnson, S. (1954). “Solution of a Large-Scale Traveling-Salesman Problem.” Journal of the Operations Research Society of America 2(4): 393–410.
• Miller, C. E., Tucker, A. W., and Zemlin, R. A. (1960). “Integer Programming Formulation of Traveling Salesman Problems.” Journal of the ACM 7(4): 326–329.
• Held, M., and Karp, R. M. (1970). “The Traveling-Salesman Problem and Minimum Spanning Trees.” Operations Research 18(6): 1138–1162.
• Desrochers, M., and Laporte, G. (1991). “Improvements and Extensions to the Miller-Tucker-Zemlin Subtour Elimination Constraints.” Operations Research Letters 10(1): 27–36.
• Padberg, M., and Sung, T.-Y. (1991). “An Analytical Comparison of Different Formulations of the Travelling Salesman Problem.” Mathematical Programming 52(1–3): 315–357.
• Öncan, T., Altınel, İ. K., and Laporte, G. (2009). “A Comparative Analysis of Several Asymmetric Traveling Salesman Problem Formulations.” Computers & Operations Research 36(3): 637–654.
• Wolsey, L. A. (1998). Integer Programming. Wiley.
• Applegate, D. L., Bixby, R. E., Chvátal, V., and Cook, W. J. (2006). The Traveling Salesman Problem: A Computational Study. Princeton University Press.
• Reinelt, G. (1991). “TSPLIB—A Traveling Salesman Problem Library.” ORSA Journal on Computing 3(4): 376–384.

See the whole series · ← All writing