I solved a Traveling Salesman Problem (TSP) by the cutting plane method, i.e. adding violated sub-tour constraints until the LP relaxation is a tour. Now, I am realizing that, in this symmetric TSP, driving route $1 \to 2 \to 3 \to 4$ is identical to driving $4 \to 3 \to 2 \to 1$.
I am wondering how I can add another constraint so that the symmetry is avoided. It should work by defining $x_{ij}$ variables only if $i < j$. However, I do not know how I can adapt my code to make it work.
I really appreciate your help!
Here is my python code:
# Load libraries
from itertools import product
from networkx import minimum_cut, DiGraph
from mip import Model, xsum, BINARY, OptimizationStatus, CutType
# This is the distance matrix
#
# 1 2 3 4 5 6 7 8 9 10
# 1 0 3 7 7 5 4 6 9 2 3
# 2 3 0 10 3 2 9 10 6 8 4
# 3 7 10 0 6 6 4 1 5 1 3
# 4 7 3 6 0 9 8 3 4 4 10
# 5 5 2 6 9 0 2 5 1 10 4
# 6 4 9 4 8 2 0 4 3 5 2
# 7 6 10 1 3 5 4 0 6 5 2
# 8 9 6 5 4 1 3 6 0 8 8
# 9 2 8 1 4 10 5 5 8 0 6
# 10 3 4 3 10 4 2 2 8 6 0
# ;
# List of locations
L = [x + 1 for x in range(10)]
# Set up distances between locations in a dictionary that uniquely identifies the nodes by a tuple
D = {
(1, 2): 3, (1, 3): 7, (1, 4): 7, (1, 5): 5, (1, 6): 4, (1, 7): 6, (1, 8): 9, (1, 9): 2, (1, 10): 3,
(2, 1): 3, (2, 3): 10, (2, 4): 3, (2, 5): 2, (2, 6): 9, (2, 7): 10, (2, 8): 6, (2, 9): 8, (2,10): 4,
(3, 1): 7, (3, 2): 10, (3, 4): 6, (3, 5): 6, (3, 6): 4, (3, 7): 1, (3, 8): 5, (3, 9): 1, (3, 10): 3,
(4, 1): 7, (4, 2): 3, (4, 3): 6, (4, 5): 9, (4, 6): 8, (4, 7): 3, (4, 8): 4, (4, 9): 4, (4, 10): 10,
(5, 1): 5, (5, 2): 2, (5, 3): 6, (5, 4): 9, (5, 6): 2, (5, 7): 5, (5, 8): 1, (5, 9): 10, (5, 10): 4,
(6, 1): 4, (6, 2): 9, (6, 3): 4, (6, 4): 8, (6, 5): 2, (6, 7): 4, (6, 8): 3, (6, 9): 5, (6, 10): 2,
(7, 1): 6, (7, 2): 10, (7, 3): 1, (7, 4): 3, (7, 5): 5, (7, 6): 4, (7, 8): 6, (7, 9): 5, (7, 10): 2,
(8, 1): 9, (8, 2): 6, (8, 3): 5, (8, 4): 4, (8, 5): 1, (8, 6): 3, (8, 7): 6, (8, 9): 8, (8, 10): 8,
(9, 1): 2, (9, 2): 8, (9, 3): 1, (9, 4): 4, (9, 5): 10, (9, 6): 5, (9, 7): 5, (9, 8):8, (9, 10): 6,
(10, 1): 3, (10, 2): 4, (10, 3): 3, (10, 4): 10, (10, 5): 4, (10, 6): 2, (10, 7): 2, (10, 8): 8, (10, 9): 6,
}
# Define movements from i to j and from j to i
Dout = {l: [d for d in D if d[0] == l] for l in L}
Din = {l: [d for d in D if d[1] == l] for l in L}
model = Model()
x = {d: model.add_var(name=f'x({d[0]},{d[1]})', var_type=BINARY) for d in D}
model.objective = xsum(d * x[l] for l, d in D.items())
# Add constraints: we want to leave and enter each location only once
for l in L:
model += xsum(x[l] for l in Dout[l]) == 1, f'out({l})'
model += xsum(x[l] for l in Din[l]) == 1, f'in({l})'
# Solve the relaxed model (without subtour elimination)
model.optimize(relax=True)
# Yields the same solution as GAMS :) {VW}
relaxed_solution = {y: x[y].x for y in x if x[y].x != 0}
# Solve TSP by the cutting plane method
violated_inequalities = True
while violated_inequalities:
model.optimize(relax=True)
# print("status: {} objective value : {}".format(m.status, m.objective_value))
# Base class for directed graphs. A DiGraph stores nodes and edges with data.
# Edges of the graph G are expected to have an attribute capacity that indicates how much flow the edge can support.
G = DiGraph()
for d in D:
G.add_edge(d[0], d[1], capacity=x[d].x)
# Separation routine:
# Choose an arbitrary location ๐ โ {1, ... , ๐}
# Solve a maximum flow problem from ๐ to each location ๐ โ {1, ... , ๐}\{๐} with edge capacities ๐ฅาง๐๐
# If the optimal objective value of one of these maximum flow problems is < 1, a subtour constraint is violated.
# The set ๐ of the violated subtour constraint is one half of the corresponding minimum cut.
# A cut divides the node set into two disjoint sets ๐ and ๐.
# The cutโs capacity is the sum of capacities of all edges between ๐ and ๐.
violated_inequalities = False
for (n1, n2) in [(i, j) for (i, j) in product(L, L) if i != j]:
# The capacity of a minimum capacity cut is equal to the flow value of a maximum flow.
cut_capacity, (S, T) = minimum_cut(G, n1, n2)
if cut_capacity < 1:
model += (xsum(x[l] for l in D if (l[0] in S and l[1] in S)) <= len(S) - 1)
violated_inequalities = True
final_solution = {y: x[y].x for y in x if x[y].x != 0}
# The tour is complete
final_solution
x[i,j] == x[j,i]
fori != j
. $\endgroup$