DM559/DM545 Assignment 1.2¶
Part of this script is taken from the very nice tutorial by Peter Norvig: http://nbviewer.jupyter.org/url/norvig.com/ipython/TSP.ipynb
Consider the Traveling Salesperson Problem:
Given a set of cities and the distances between each pair of cities, what is the shortest possible tour that visits each city exactly once, and returns to the starting city?
Representing an instance of the problem¶
The only thing that matters about cities is the distance between them. We will ignore the fully general TSP where distances can be defined in any arbitrary way and concentrate on an important special case, the Euclidean TSP, where the distance between any two cities is the Euclidean distance, the straight-line distance between points in a two-dimensional plane. So a city can be represented by a two-dimensional point: a pair of $x$ and $y$ coordinates in the Cartesian plane. Python already has the two-dimensional point as a built-in numeric data type, but in a non-obvious way: as complex numbers, which inhabit the two-dimensional (real $\times$ imaginary) plane. We will use those but you will not have to care about this aspect as all functions to handle data are made available to you.
We will work with three predetermined, historical instances, dantzig42.dat berlin52.dat bier127.dat, and with randomly generated instances.
You can download the three instances from these links: dantzig42.dat,
berlin52.dat,
bier127.dat.
In the file tsputil.py you will find the functions to read the instances from the files and to generate instances at random.
The constructor function City, creates a city object, so that City(300, 0) creates a city with x-coordinate of 300 and y-coordinate of 0.
Then, distance(A, B) will be a function that uses the $x$ and $y$ coordinates to compute the distance between cities A and B.
Let's import the files:
from tsputil import *
from gurobipy import *
from collections import OrderedDict
%run ../Models/TSP/python_pub/tsputil.py
%matplotlib inline
We can then generate an instance of random points with the function Cities(n,seed). The function returns a frozen set because these are the input data and cannot be modified. We plot the instance with plot_situation. When you generate an instance make sure that you use a seed number different from the one of other groups working at this project.
ran_points = list(Cities(n=20,seed=35))
plot_situation(ran_points)
Alternatively, we can read the dantiz42.dat instance which represents locations in USA.
dantzig42 = read_instance("dantzig42.dat")
plot_situation(dantzig42)
Dantzig, Fulkerson and Johnson (DFJ) Formulation¶
Consider the following formulation of the symmetric traveling salesman problem due to Dantzig, Fulkerson and Johnson, 1954 (DFJ) that we have been discussing during the course. Let $V=\{0..n-1\}$ be the set of nodes and $E$ the set of edges. Let $E(S)$ be the set of edges induced by the subset of vertices $S$ and $\delta(v)$ the set of edges in the cut $(v,V\setminus\{v\})$. (We will assume that an edge between two nodes $i$ and $j$ is present in $E$ in the form $ij$ if $j>i$ and $ji$ if $ j < i $.) $$\begin{aligned} \text{(TSPIP)}\quad \min\; & \sum c_{ij} x_{ij} \\ \text{s.t.}\;&\sum_{ij \in \delta(i)} x_{ij}+\sum_{ji \in \delta(i)}x_{ji}=2 \text{ for all } i \in V\\ \label{subtour}&\sum_{ij \in E(S)} x_{ij} \leq |S|-1 \text{ for all } \emptyset \subset S\subset V, 2 \leq |S| \leq n-1\\ &x_{ij} \in \{0,1\} \text{ for all } {ij} \in E\end{aligned}$$
We can generate all subsets of the set of 20 randomly generated cities as follows:
from itertools import chain, combinations
def powerset(iterable):
"powerset([1,2,3]) --> () (1,) (2,) (3,) (1,2) (1,3) (2,3) (1,2,3)"
s = list(iterable)
return chain.from_iterable(combinations(s, r) for r in range(len(s)+1))
subtours = list(powerset(range(len(ran_points))))
# The first element of the list is the empty set and the last element is the full set, hence we remove them.
subtours = subtours[1:(len(subtours)-1)]
print(len(subtours))
import sys
print(sys.getsizeof(subtours)/1024/1024," GB")
Task 1¶
Implement the DFJ formulation in the model below and solve your random instance:
def solve_tsp(points, subtours=[], VAR_TYPE=GRB.INTEGER, silent=True):
points=list(points)
V = range(len(points))
E = [(i,j) for i in V for j in V if i<j] # complete graph
E = tuplelist(E)
m = Model("TSP0")
if silent:
m.setParam(GRB.Param.OutputFlag,0)
m.setParam(GRB.param.Presolve, 0) # no preprocessing
m.setParam(GRB.param.Method, 0)
m.setParam(GRB.param.MIPGap,1e-7)
######### BEGIN: Write here your model for Task 1
x = OrderedDict()
for (i,j) in E:
x[i,j]=m.addVar(lb=0.0, ub=1.0, obj=distance(points[i],points[j]), vtype=VAR_TYPE, name="x[%d,%d]" % (i,j))
## Objective
m.modelSense = GRB.MINIMIZE
## Constraints
for v in V:
m.addConstr(quicksum(x[v,i] for v,i in E.select(v,'*')) +
quicksum(x[i,v] for i,v in E.select('*',v)) == 2,"flow_balance_%d" % v)
for S in subtours:
m.addConstr(quicksum(x[i,j] for i in S for j in S if i<j)<=len(S)-1,"tour_elim_%d" % sum(S))
######### END
m.optimize()
m.write("tsplp.lp")
if m.status == GRB.status.OPTIMAL:
print('Solve TSP: the optimal objective value is %g' % m.objVal)
m.write("tsplp.sol") # write the solution
return {(i,j) : x[i,j].x for i,j in x}
else:
print("Something wrong in solve_tsplp")
exit(0)
If we try to solve the small instance on 20 cities with this model, ie:
solve_tsp(ran_points, subtours)
we get out of memory (you can be more lucky but if so try increasing your instance with a few more cities).
solve_tsp(ran_points,subtours,silent=False)
Changed value of parameter Presolve to 0
Prev: -1 Min: -1 Max: 2 Default: -1
Changed value of parameter Method to 0
Prev: -1 Min: -1 Max: 4 Default: -1
Changed value of parameter MIPGap to 1e-07
Prev: 0.0001 Min: 0.0 Max: 1e+100 Default: 0.0001
Optimize a model with 1048594 rows, 190 columns and 49807550 nonzeros
Variable types: 0 continuous, 190 integer (0 binary)
Coefficient statistics:
Matrix range [1e+00, 1e+00]
Objective range [3e+01, 9e+02]
Bounds range [1e+00, 1e+00]
RHS range [1e+00, 2e+01]
Explored 0 nodes (0 simplex iterations) in 0.34 seconds
Thread count was 1 (of 8 available processors)
Solution count 0
Pool objective bound 0
Best objective -, best bound -, gap -
---------------------------------------------------------------------------
GurobiError Traceback (most recent call last)
<ipython-input-7-8784813f6174> in <module>()
----> 1 solve_tsp(ran_points,subtours)
<ipython-input-6-8bb38fbd53a5> in solve_tsp(points, subtours)
27 ######### END
28
---> 29 m.optimize()
30 m.write("tsplp.lp")
31
model.pxi in gurobipy.Model.optimize (../../src/python/gurobipy.c:54012)()
GurobiError: Out of memory
</span>
A Cutting Plane Approach¶
The number of subtour elimination constraints in the DFJ formulation can be huge ($2^n$) and even though we can remove half of those due to symmetries, there are still exponentially many. A way to deal with large families of constraints is to introduce them only when needed. This idea is similar to the cutting plane algorithm that we saw in the intro classes. We start by solving the relaxed version of (TSPIP) obtained by removing the integrality constraints and the exponentially many subtour elimination constraints. We dub this relaxation TSPLP$_0$. Let $x^*$ be the optimal solution of TSPLP$_0$.
Task 2¶
Implement the new function solve_tsp(points, subtours=[]) that solves the relaxed version TSPLP$_0$ and solve your random instance with 20 points. Inspect the plot of the solution obtained.
tsplp0 = solve_tsp(ran_points, [])
plot_situation(ran_points, tsplp0)
Describe what you observe. Are all variables integer? Is the matrix TUM? Why? Is the solution a tour? Is the solution matching your expectations?
Relaxing only the subtour elimination constraints in solve_tsp we obtain:
tsplp0 = solve_tsp(ran_points, [], GRB.INTEGER, silent=False)
plot_situation(ran_points, tsplp0)
Relaxing also the integrality constraints `vtype=GRB.INTEGER` is more convenient since the solution will be faster:
tsplp0 = solve_tsp(ran_points, [], VAR_TYPE=GRB.CONTINUOUS, silent=False)
plot_situation(ran_points, tsplp0)
You should pay attention to the two different outputs of gurobi. In particular, in the latter output with continuous variables you find the log of the iterations of the primal simplex, while in the output of the model with interger variables you find the log of the branch and bound process in which the 48 iterations of the LP relaxation at the root is summarized in one single row. In the linear relaxation the variables are not all integer. The matrix is indeed not TUM because although we have only two terms different from zero in each column of the matrix, they are both positive and the graph is not bipartite, hence we are not able to separate the rows in the fashion described by the theorem. The solution is not a tour, since, as expected, there are subtours and fractional variables.
An alternative directed graph formulation. Note that if we had used a directed representation of the problem, by introducing for each edge two oriented arcs then we would have had a TUM matrix, although the number of arcs would have been now double. Let $A=\{ij\mid i\in V,j\in V,i\neq j\}$ and let $\delta^+(i)$ and $\delta^-(i)$ the set of arcs, respectively, outgoing and ingoing $i$: $$\begin{aligned} \text{(TSPIP)}\quad \min\; & \sum_{ij \in A} c_{ij} x_{ij} \\ \text{s.t.}\;&\sum_{ij \in \delta^+(i)} x_{ij}=1 \text{ for all } i \in V\\ &\sum_{ji \in \delta^-(i)} x_{ji}=1 \text{ for all } i \in V\\ &\sum_{ij \in A(S)} x_{ij} \leq |S|-1 \text{ for all } \emptyset \subset S\subset V, 2 \leq |S| \leq n-1\\ &x_{ij} \in \{0,1\} \text{ for all } {ij} \in E \end{aligned}$$
def solve_tsp_directed(points, subtours=[], VAR_TYPE=GRB.INTEGER, silent=True):
points=list(points)
V = range(len(points))
A = [(i,j) for i in V for j in V if i!=j] # complete graph
A = tuplelist(A)
m = Model("TSP0")
if silent:
m.setParam(GRB.Param.OutputFlag,0)
m.setParam(GRB.param.Presolve, 0) # no preprocessing
m.setParam(GRB.param.Method, 0)
m.setParam(GRB.param.MIPGap,1e-7)
######### BEGIN: Write here your model for Task 1
x = OrderedDict()
for (i,j) in A:
x[i,j]=m.addVar(lb=0.0, ub=1.0, obj=distance(points[i],points[j]), vtype=VAR_TYPE, name="x[%d,%d]" % (i,j))
m.update()
## Objective
m.modelSense = GRB.MINIMIZE
## Constraints
for v in V:
m.addConstr(quicksum(x[i,v] for i,v in A.select('*',v)) == 1,"incoming_flow_balance_%d" % v)
m.addConstr(quicksum(x[v,i] for v,i in A.select(v,'*')) == 1,"outgoing_flow_balance_%d" % v)
for S in subtours:
m.addConstr(quicksum(x[i,j] for i in S for j in S if i!=j)<=len(S)-1,"tour_elim_%d" % sum(S))
######### END
m.optimize()
m.write("tsplp.lp")
if m.status == GRB.status.OPTIMAL:
print('Solve TSP: the optimal objective value is %g' % m.objVal)
m.write("tsplp.sol") # write the solution
return {(i,j) : x[i,j].x for i,j in x}
else:
print("Something wrong in solve_tsplp")
exit(0)
tsplp0 = solve_tsp_directed(ran_points, [], GRB.CONTINUOUS, False)
plot_situation(ran_points, tsplp0)
As we see, an integer solution is found without need of a branch and bound (see the output of gurobi). However, a side effect of considering a directed graph formulation is the fact that subtours of lenght 2 are not eliminated and we actually get 4 of them in the case above.
Task 3¶
Detect by visual inspection some subtour inequalities to add to TSPLP0. List those subtours in the subtour argument of the function solve_tsp and solve the new problem TSPLP$_1$=TSPLP$_0 \cup c$, where $c$ is the constraint added.
Report the visualization produced by plot_situation on the new solution. Is it a tour? If not iterate this procedure until you cannot find anymore subtours. Show your last solution and comment on it.
tsplp0 = solve_tsp(ran_points, [], GRB.CONTINUOUS)
plot_situation(ran_points, tsplp0)
tsplp1 = solve_tsp(ran_points, [[2,7,9], [0,13,5], [3,6,16,14,1,11,10]], GRB.CONTINUOUS)
plot_situation(ran_points, tsplp1)
tsplp1 = solve_tsp(ran_points, [[2,7,9], [0,13,5], [3,6,16,14,1,11,10],
[3,6,16,14,1,11,18,10],[2,7,9,17]], GRB.CONTINUOUS)
plot_situation(ran_points, tsplp1)
tsplp1 = solve_tsp(ran_points, [[2,7,9], [0,13,5], [3,6,16,14,1,11,10],
[3,6,16,14,1,11,18,10],[2,7,9,17],
[6,16,14,1,12,8,18,11,10,3],[2,7,9,0,13,17],
[2,7,9,0,5,13,17],[11,14,1],[4,19,15,12,8,18,10,3,6,16]], GRB.CONTINUOUS)
plot_situation(ran_points, tsplp1)
tsplp1 = solve_tsp(ran_points, [[2,7,9], [0,13,5], [3,6,16,14,1,11,10],
[3,6,16,14,1,11,18,10],[2,7,9,17],
[6,16,14,1,12,8,18,11,10,3],[2,7,9,0,13,17],
[2,7,9,0,5,13,17],[11,14,1],[4,19,15,12,8,18,10,3,6,16],
[2,7,9,0,13,5,4,17],[19,15,12,8,18,11,10,3,6,16,14,1]], GRB.CONTINUOUS)
plot_situation(ran_points, tsplp1)
tsplp1 = solve_tsp(ran_points, [[2,7,9], [0,13,5], [3,6,16,14,1,11,10],
[3,6,16,14,1,11,18,10],[2,7,9,17],
[6,16,14,1,12,8,18,11,10,3],[2,7,9,0,13,17],
[2,7,9,0,5,13,17],[11,14,1],[4,19,15,12,8,18,10,3,6,16],
[2,7,9,0,13,5,4,17],[19,15,12,8,18,11,10,3,6,16,14,1],
[6,3,10]], GRB.CONTINUOUS)
plot_situation(ran_points, tsplp1)
We might finally find a tour continuing in this way and that would be an optimal solution since it would satisfy all DFJ constraints. However, in the general case we could have no subtour constraints unsatisfied but a solution not yet integer. This is for example the case for the instance `dantzig42` as we will see below.
Task 4¶
We can automatize the manual process described above. Let $V'=V\setminus\{0\}$ and $E'=E\setminus\{\delta(0)\}$. A set $S\subseteq V'$ that contains a vertex $k$ forms a subtour in $x^*$ if $$\sum_{e=ij \in E'(S)} x_e^* > |S \setminus \{k\} |=|S|-1.$$
We can find the subtour elimination constraint that is most violated by $x^*$ by solving the following separation problem for $k=\{1..n-1\}$. Let’s represent the unknown set $S \subseteq V'$ by binary variables $z$ with $z_i=1$ if $i \in S$ then the separation problem for $k=\{1..n-1\}$ is
$$ (SEP) \quad \xi_k=\max \left\{ \sum_{e=ij\in E':i < j} x^*_e z_i z_j - \sum_{i \in V' \setminus \{k\}} z_i : z \in \mathbb{B}^n, z_k=1\right\} $$
If $ \xi_k > 0 $ for some $k$ then the solution $z^*$ of SEP gives a subtour $S$, containing $k$, whose corresponding subtour elimination constraint in TSPIP is violated by $x^*$. We insert this constraint in TSPLP$_0$ and solve the new problem TSPLP$_1$.
Rewrite the (SEP) problem for some $k$ as an integer linear program and report the mathematical model in your answers.
Solution
We introduce the variables $y_{e}$ with $y_{e}=1$ if $z_i=z_j=1$ for $e=ij \in E'$ with $i<j$, and we obtain the integer program
$$\begin{aligned} \text{(SEPILP)}\quad \xi_k=\max\; & \sum_{e \in E'} x_e^* y_e-\sum_{i \in V'\setminus\{k\}} z_i \\ \label{sep1}\text{s.t.}\;&y_e \leq z_i \text{ for } e=ij \in E'\\ \label{sep2}&y_e \leq z_j \text{ for } e=ij \in E'\\ \label{sep3}&y_e \geq z_i+z_j-1 \text{ for } e=ij \in E'\\ &z_k=1\\ &y \in \{0,1\}^{|E'|}, z\in \{0,1\}^{|V'|}\end{aligned}$$ End Solution </span>
Task 5¶
Implement your model in the function solve_separation(points, x_star, k) below and solve the LP relaxation of the (SEP) problem to separate the optimal solution $x^*$ of TSPLP$_0$.
def solve_separation(points, x_star, k, VAR_TYPE=GRB.CONTINUOUS, silent=True):
points=list(points)
V = range(len(points))
Vprime = range(1,len(points))
E = [(i,j) for i in V for j in V if i<j]
Eprime = [(i,j) for i in Vprime for j in Vprime if i<j]
E = tuplelist(E)
Eprime = tuplelist(Eprime)
m = Model("SEP")
if silent:
m.setParam(GRB.Param.OutputFlag,0)
m.setParam(GRB.Param.Presolve, 0)
m.setParam(GRB.Param.Method, 0)
m.setParam(GRB.Param.MIPGap,1e-7)
######### BEGIN: Write here your model for Task 4
y={}
for i,j in Eprime:
y[i,j]=m.addVar(lb=0.0,ub=1.0,obj=x_star[i,j],vtype=VAR_TYPE,name="y[%d,%d]" % (i,j))
z={}
for v in Vprime:
if v!=k:
z[v]=m.addVar(lb=0.0,ub=1.0,obj=-1.0,vtype=VAR_TYPE,name="z[%d]" % (v))
else:
z[v]=m.addVar(lb=0.0,ub=1.0,obj=0.0,vtype=VAR_TYPE,name="z[%d]" % (v))
m.update()
# Objective
m.modelSense = GRB.MAXIMIZE
# Constraints
for i,j in Eprime:
m.addConstr(y[i,j]<=z[i])
m.addConstr(y[i,j]<=z[j])
m.addConstr(y[i,j]>=z[i]+z[j]-1)
m.addConstr(z[k]==1)
######### END
m.optimize()
m.write("sep.lp")
if m.status == GRB.status.OPTIMAL:
m.write("sep.sol") # write the solution
subtour = list(filter(lambda i: z[i].x>=0.99,z))
return m.objVal, subtour
else:
print("Something wrong in solve_tsplp")
exit(0)
Try your implementation of solve_separation on your previous TSPLP$_0$ solution. Do you get the expected answer?
Solution
tsplp0 = solve_tsp(ran_points, [], GRB.CONTINUOUS)
plot_situation(ran_points, tsplp0)
S_star = solve_separation(ran_points, tsplp0, 12, VAR_TYPE=GRB.CONTINUOUS)
print(S_star)
This means that by solving for the vertex k=12 we find the set S that contains 12 and violates the most the subtour elimination constraint. We should now solve for all k.
Task 6¶
The following procedure cutting_plane_alg implements the cutting plane algorithm that uses your implementation of solve_tsp and solve_separation.
Finish the implementation of the function. That is: add the condition for the
presence of violated inequalities.
Solution
def cutting_plane_alg(points, silent=True):
Vprime = range(1,len(points))
subtours = []
found = True
while found:
lpsol = solve_tsp(points,subtours, GRB.CONTINUOUS)
if not silent:
plot_situation(points, lpsol)
found = False
tmp_subtours = []
best_val = float('-inf')
for k in Vprime:
value, subtour = solve_separation(points,lpsol,k,GRB.CONTINUOUS)
if not silent:
print('Separation problem solved for k=%d, solution value %g' % (k,value))
best_val = value if value > best_val else best_val
######### BEGIN: write here the condition. Include a tollerance
if value > 0.0001:
######### END
found = True
tmp_subtours += [subtour]
subtours += tmp_subtours
if not silent:
print('*'*60)
print("Subtours found: ",tmp_subtours,"\nwith best value : ",best_val)
print('*'*60)
plot_situation(points, lpsol)
We added the condition `value>0`.
Run the cutting plane algorithm thus implemented until termination. Report the plot of the last solution found by solving the TSPLP$_n$ problem, where $n$ is the number of iterations in which violated inequalities were found. Is the solution optimal?
Solution
cutting_plane_alg(ran_points, False)
We find the same optimal solution as eariler.
End Solution
Task 7 (Optional)¶
Reflecting upon the process implemented you may wonder whether we are gaining anything in efficiency. After all we circumvented having to solve the original integer programming formulation of the TSP by having to solve many separation problems, which are in principle all integer programming problems and hence possibly hard to solve. It turns out that the SEP problems are actually solvable efficiently. In task 5, you were given the tip to implement the LP relaxation of the SEP problem. Why are the solutions to SEP always integer? It turns out that the SEP corresponds to one problem treated during the course. Which one? In efficient implementations the SEP problem is solved by specialized algorithms for this problem. [Hint: It may help to consider the other way to formulate subtour elimination constraints, that is, cut-set constraints.]
Solution
Solutions to the LP relaxation of (SEP) are always integer because the matrix is TUM. (Since $x^*\geq 0$ then we wish $y_e$ as large as possible. Thus by - $y_e=\min(z_i,z_j)$ and $\min(z_i,z_j)\geq z_i+z_j-1$ is always satisfied and can be dropped. We are left with a matrix with a block that contains only two consecutive ones on each column (for the variables $y$) concatenated with two identity matrices (for the variables $z$). Such a matrix can be shown to be TUM (consecutive ones can be rearranged and a partition satisfying the sufficiency conditions found + concatation with identity matrices preserve TUM properties).
The separation problem can be solved by running a max flow algorithm from any vertex $k$ as source to any other vertex as target and $x^*$ as capacity on the arcs oriented from source to tank. This can be better seen from the dual of the (SEP) problem that gives an arc-incidence matrix. Another way to see this is the following. Consider the alternative cut set constraints in place of the subtour elimination constraints : $$\sum_{e\in \delta(S)} x_e\geq 2, \forall \emptyset \subset S \subset V$$ which imposes that every solution must cross every demarcation line separating the set of cities into two nonempty parts at least twice. Let $N_{kt}(x)$ be the network with arcs obtained from $E$ oriented from a source vertex $k\in V$ to an arbitrary fixed tank vertex $t\in V$ and capacity $x^*$ on these arcs. Any nonempty proper subset $S$ of $V$ that generates a cut $N_{kt}(x)$ with $k$ in $S$ and $t$ in $V\setminus S$ with capacity less than 2 corresponds to a set $S$ violating the cut separation constraint. One such set, if it exists, can be found by solving a global minimum cut problem in $N_{kt}(x)$ from any vertex $k$ to an arbitrary fixed vertex $t$ (the choice does not matter). If we use the push-relabel algorithm the complexity of each single max flow-min cut problem is $O(n^2m)$ and since we need to solve it $n-1$ times for all $k$ the total complexity of the separation problem becomes $O(n^4m)$. The procedure is implemented in Concorde in a very efficient way @AppBixChvCoo06 [page 170].
The problem to solve is the global min cut problem, which can be solved solving $n-1$ min cut problems. The MIP formulation of the min cut problem is the dual of the dual the maximum flow problem. Hence for $k=\{2..n\}$ and an arbitrarily chosen $t$: $$\begin{aligned} \min\;& \alpha_{ij}x^*_{ij}\\ &u_i-u_j+\alpha_{ij}\geq 0 \quad \forall ij \in E'\\ &u_k-u_t\geq 1\\ &u_k=1\\ &u_i\in \{0,1\}\quad \forall i \in V'\\ &\alpha_{ij}\geq 0 \quad \forall ij \in E'\\\end{aligned}$$
End Solution
Task 8 (Optional)¶
In task 6, you run the cutting plane algorithm until no cutting
plane can be found anymore. Is this cutting plane algorithm enough
to solve your instance to optimality? What about the general case on
any arbitrary instance? Write your considerations and conclusions.
[Hint: try the whole procedure on the instances dantzig42.dat and berlin52.dat to see what happens in the general case. Report the
final plots on these instance.]
Solution
As argued before, the solution found to the random instance is optimal.
For the instance with 20 vertices it would have been possible to
generate all cuts at the beginning and solve the IP once. For the
instance berlin52.dat, the number of constraints would be simply
too large to enumerate. We need therefore to use the dynamic
procedure illustrated in this assignment. If we apply the procedure to
the berlin52.dat example we finish again with an integer solution
that is a tour and hence it is optimal. For the instance dantzig42
however the last solution, TSPLP$_5$ has no subtour elimination
constraint to add, that is $\xi\leq 0$ for all $k$. Nevertheless
the solution is not integer. See Figure. This is
due to the fact that the formulation of the TSP given at the
beginning of Task 2 is not tight, that is, it is not the description
of the convex hull. We do not know the formulation of the convex
hull for the TSP problem. What one could do is to proceed by
applying further tighter constraints such as the comb constraints.
Alternatively one could reinsert the integrality constraints in the
problem TSP$_5$ and continue the cutting plane procedure with branch and bound from there.
cutting_plane_alg(dantzig42,True)
Task 10 (Optional)¶
Implement the Miller, Tucker and Zemlin (MTZ) formulation
from Ex. 7 Sheet 7. Compare the results of the DFJ (that is the name
of our previous TSP formulation from the authors Dantzig,
Fulkerson, Johnson) and MTZ formulations first on your instance on
20 nodes and then on the instances berlin52.dat and dantzig42.
In both formulations keep the integrality constraints of the main
variables so that an optimal solution can be found. Comment the
results in your answers: which is the fastest formulation to solve
the problem? How many branch and bound nodes are needed in the MTZ
formulation to find the optimal solution?
Compare then the results on the instance berlin52.dat. Comment the
results in your answers: can the MTZ solve the problem? Which of the
two formulations provides the best LP relaxation?
Compare the DFJ and MTZ formulations on the basis of the
quality of their linear relaxations on the instance dantzig42.
Already an instance on 127 vertices is challenging for
our implementation. Try your cutting plane procedure on the instance
bier127.dat, which asks to find the shortest tour among the 127
biergardens in Augsburg, Bavaria. You will have to comment the line
for plotting at each iteration, since it is inefficient. Moreover it
might be convenient not to wait for all $k$s to be evaluated but to
insert a cut as soon as one is found.
Experiment ideas to solve the instance. For example, try combining the cutting plane procedure with integrality constraints in the DFJ formulation or combine the two formulations DFJ and MTZ.
The following formulation (C. E. Miller, A. W. Tucker, and R. A. Zemlin, “Integer programming formulations and traveling salesman problems,” J. ACM, 7 (1960), pp. 326–329) which also eliminates subtours has only a polynomial number of constraints. Let $x_{ij}$ be the variables on the edges with the same meaning as for undirected graphs described above and let $u_{i}$, $i=0,1,\ldots,n-1$ be the number of cities visited already when we arrive at city $i$. $$\begin{aligned} \min & \sum_{(ij) \in A} c_{ij} x_{ij}\\ &\sum_{i:i\neq j} x_{ij}=1 & \forall j=0,\ldots,n-1 \\ &\sum_{j:i\neq j} x_{ij}=1 &\forall i=0,\ldots,n-1 \\ &u_i - u_j + n x_{ij} \leq n - 1, & \forall i, j = 1,2,\ldots, n-1, i\neq j\\ &x_{ij}\in \mathbb{B} &\forall i,j,i\neq j\\ &u_i\in \mathbb{R} & \forall i = 1,\ldots,n-1 \end{aligned}$$
Note that the third constraint can be equivalently written as: $$u_i +1 \leq u_j + n (1-x_{ij}), \qquad i, j = 1, 2,\ldots, n-1, i\neq j $$
Let's fix a vertex, say $0$, to be the home base and for each other vertex $i$ let $u_i$ be an arbitrary real number. The $u_i$ variables play a role similar to node potentials in a network and the inequalities involving them serve to eliminate tours that do not begin and end at city $0$ and tours that visit more than $n-1$ cities. Consider a feasible solution to the TSP: it can be expressed as a permutation of cities, $\pi=(\pi_1=0,\pi_2,\ldots\pi_n,\pi_1=0)$. Hence $x_{\pi_1,\pi_2}=1$. Unless $\pi_2=0$ then there exists $\pi_3$ such that $x_{\pi_2\pi_3}=1$. We proceed in this way until some $\pi_i=0$, which if the solution is feasible happens only at the end of the permutation. The constraint forces $u_{\pi_j}\geq u_{\pi_i}+1$ when $x_{\pi_i,\pi_j}=1$, except when $\pi_{i+1}=0$. Hence an assignment of $\vec u$ can be found when the solution is feasible. On the contrary if there were subtours or a node was visited more than once, then an assignment for the variables $u_i$ that satisfies $u_{\pi_j}\geq u_{\pi_i}+1$ when $x_{\pi_i,\pi_j}=1$ could be found (the constraint would be violated where the subtour closes).
To be continued