Modifying tsp.py for generic ordering of cities/nodes.
AnsweredIt seems as if this file only handles ordered nodes starting from 0, like these N = [0,1,2,3,...,n]. However, I got the cities N = [1,0,4,8] where 0 is the depot. The program runs, but it constructs a tour with a single city and that city is 2, which isn't even in my set of cities N. I don't get why this happens. Below is a completely reproducible code snippet.
I also have one more question. Does the optimum result of running tsp.py also take into consideration the return to depot? It seems as if when the tour is printed, 0->2->3... it does not end with a 0 again. Ofcourse one can just add a 0 when printing but how do I change it so that the cost of goin from last city back to start is counted in the minimum objective?
------Code------
#######################################################################
# Fix the ordering of the tour to minimize cost via TSP
#######################################################################
from itertools import combinations
import gurobipy as gp
from gurobipy import GRB
# Callback - use lazy constraints to eliminate sub-tours
def subtourelim(model, where):
if where == GRB.Callback.MIPSOL:
vals = model.cbGetSolution(model._vars)
# find the shortest cycle in the selected edge list
tour = subtour(vals)
if len(tour) < n:
# add subtour elimination constr. for every pair of cities in tour
model.cbLazy(gp.quicksum(model._vars[i, j]
fori, jincombinations(tour, 2))
<= len(tour)-1)
# Given a tuplelist of edges, find the shortest subtour
def subtour(vals):
# make a list of edges selected in the solution
edges = gp.tuplelist((i, j) for i, j in vals.keys()
if vals[i, j] > 0.5)
unvisited = list(range(n))
cycle = range(n+1) # initial length has 1 more city
while unvisited: # true if list is non-empty
thiscycle = []
neighbors = unvisited
while neighbors:
current = neighbors[0]
thiscycle.append(current)
unvisited.remove(current)
neighbors = [j for i, j in edges.select(current, '*')
if j in unvisited]
if len(cycle) > len(thiscycle):
cycle = thiscycle
returncycle
N = [1,0,4,8] # The cities
n = len(N) # number of cities
# distances between cities. Note that for example (1,4) has been removed
# since dist(4,1) = dist(1,4) etc.
dist = {(4, 1): 63.34824385884742,
(8, 1): 26.1725046566048,
(1, 0): 79.42921376924235,
(8, 4): 39.395431207184416,
(4, 0): 119.20570456148481,
(8, 0): 85.47514258543241}
m = gp.Model()
# Create variables
x = m.addVars(dist.keys(), obj=dist, vtype=GRB.BINARY, name='x')
for i,j in x.keys():
x[j, i] = x[i, j] # edge in opposite direction
# Add degree-2 constraint
m.addConstrs(x.sum(i, '*') == 2 for i in N)
# Optimize model
m._vars = x
m.Params.LazyConstraints = 1
m.optimize(subtourelim)
vals = m.getAttr('x', x)
tour = subtour(vals)
#assert len(tour) == n
print('')
print('Optimal tour: %s' % str(tour))
print('Optimal cost: %g' % m.ObjVal)
print('')
-----Log-----
Gurobi Optimizer version 9.5.2 build v9.5.2rc0 (win64)
Thread count: 4 physical cores, 8 logical processors, using up to 8 threads
Optimize a model with 4 rows, 6 columns and 12 nonzeros
Model fingerprint: 0x304fd742
Variable types: 0 continuous, 6 integer (6 binary)
Coefficient statistics:
Matrix range [1e+00, 1e+00]
Objective range [3e+01, 1e+02]
Bounds range [1e+00, 1e+00]
RHS range [2e+00, 2e+00]
Presolve removed 3 rows and 3 columns
Presolve time: 0.00s
Presolved: 1 rows, 3 columns, 3 nonzeros
Variable types: 0 continuous, 3 integer (3 binary)
[0]
Found heuristic solution: objective 294.2015957
Root relaxation: objective 2.642029e+02, 1 iterations, 0.00 seconds (0.00 work units)
[2]
Nodes | Current Node | Objective Bounds | Work
Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
* 0 0 0 264.2028542 264.20285 0.00% - 0s
Explored 1 nodes (1 simplex iterations) in 0.02 seconds (0.00 work units)
Thread count was 8 (of 8 available processors)
Solution count 2: 264.203 294.202
Optimal solution found (tolerance 1.00e-04)
Best objective 2.642028541945e+02, best bound 2.642028541945e+02, gap 0.0000%
User-callback calls 188, time in user-callback 0.00 sec
Optimal tour: [2]
Optimal cost: 264.203
-
Official comment
-
Gurobi Staff
This post is more than three years old. Some information may not be up to date. For current information, please check the Gurobi Documentation or Knowledge Base. If you need more help, please create a new post in the community forum, or try Gurobot, our chatbot interface offering instant, expert-level support. -
-
-
I am not sure if I really understand the issue here. Can't you just rename your cities and order them as is expected by the tsp.py example?
N = [1,0,4,8] # The cities
n = len(N) # number of cities
namemap = {}
for i in range(n):
namemap[i] = N[i]
dist = {(0, 1): 79.42921376924235,
(0, 2): 63.34824385884742,
(0, 3): 26.1725046566048,
(1, 2): 119.20570456148481,
(1, 3): 85.47514258543241,
(2, 3): 39.395431207184416}
m = gp.Model()
# Create variables
x = m.addVars(dist.keys(), obj=dist, vtype=GRB.BINARY, name='x')
for i,j in x.keys():
x[j, i] = x[i, j] # edge in opposite direction
# Add degree-2 constraint
m.addConstrs(x.sum(i, '*') == 2 for i in range(n))
# Optimize model
m._vars = x
m.Params.LazyConstraints = 1
m.optimize(subtourelim)
vals = m.getAttr('x', x)
tour = subtour(vals)
#assert len(tour) == n
print('')
for i in range(len(tour)):
tour[i] = namemap[i]
print('Optimal tour: %s' % str(tour))
print('Optimal cost: %g' % m.ObjVal)
print('')Does the optimum result of running tsp.py also take into consideration the return to depot?
Yes, the final tour is a cycle. Just the output does not make it clear.
I also don't really understand why a depot matters for tsp? You are searching for a cycle visiting all cities, so any city can be a depot no?
Best regards,
Jaromił0 -
-
My understanding of the TSP is that given a set of cities N = [1,2,...,n], it returns the order in which the cities should be visited such as to minimize the travel distance. In the example above my cities are, unordered, N = [1,0,4,8]. I want to know, starting from 0, how can I visit 1,4, 8 and return to zero while minimizing the distance.
Running the code you provided just returns [1,0,4,8] which doesn't tell me anything about the order. I assume, ,since this is a cycle, I can rearrange the output to be [0,4,8,1,0] but I get the same regardless of N.
Obviously the output [1,0,4,8] could be the optimal route but I get the same regardless of problem size, even If I try n = 50, the tsp.py never changes the ordering but just returns N.
0 -
-
-
Running the code you provided just returns [1,0,4,8] which doesn't tell me anything about the order. I assume, ,since this is a cycle, I can rearrange the output to be [0,4,8,1,0] but I get the same regardless of N.
Correct, it looks like [1,0,4,8,1] is just the optimal solution in your case.
If I try n = 50, the tsp.py never changes the ordering but just returns N.
Could you post an example where the final solution is actually not optimal? When I run the tsp.py example without your changes, the optimal route is very different from just [0,1,2,3,...].
0 -
-
Here is a completely reproducible example. Note that instead of hardcoding dist, I just calculated all distances between 0,...,max(N) in dist and then create a new dictionary dists that is passed to the model, determined by the contents of N. In the code below, no matter what N is, that particular order is always printed as output from tsp.py. The heuristic can't be that good that it accidentally constructs optimal N each time.
Just do describe my full problem: I'm using a construction heuristic to determine feasible routes for a cVRP. Once each vehicle k has been assigned all it's customers such that their combined demand does not exceed the maximum load capacity of the vehicle, I decide the optimal order for each vehicle using the TSP. This means for example that if I have ALL the customers A = [1,2,3,4,5,6,7,8,9,10,11,12] (0 = depot) and the heuristic assigns
vehicle_1 = [12,8,9,4]
vehicle_1 = [2,7,3,11]
vehicle_1 = [10,1,6,5]
I'll myself just insert the depot 0 in each of the above using vehicle_k.insert(0,0). Now, I want this TSP program to output the optimal route for each vehicle as vehicle_1_opt =[0,...,0] and the optimal order should be in between. In the example below, my N represents only one vehicle. Once I fix the case with one vehicle then it's trivial to just loop/iteratively solve the TSP for each vehicle and sum their costs. I hope my explanation makes sense.
from itertools import combinations, permutations
import gurobipy as gp
from gurobipy import GRB
# Callback - use lazy constraints to eliminate sub-tours
def subtourelim(model, where):
if where == GRB.Callback.MIPSOL:
vals = model.cbGetSolution(model._vars)
# find the shortest cycle in the selected edge list
tour = subtour(vals)
if len(tour) < n:
# add subtour elimination constr. for every pair of cities in tour
model.cbLazy(gp.quicksum(model._vars[i, j]
for i, j in combinations(tour, 2))
<= len(tour)-1)
# Given a tuplelist of edges, find the shortest subtour
def subtour(vals):
# make a list of edges selected in the solution
edges = gp.tuplelist((i, j) for i, j in vals.keys()
if vals[i, j] > 0.5)
unvisited = list(range(n))
cycle = range(n+1) # initial length has 1 more city
while unvisited: # true if list is non-empty
thiscycle = []
neighbors = unvisited
while neighbors:
current = neighbors[0]
thiscycle.append(current)
unvisited.remove(current)
neighbors = [j for i, j in edges.select(current, '*')
if j in unvisited]
if len(cycle) > len(thiscycle):
cycle = thiscycle
return cycle
random.seed(1)
N = [12,8,9,4] # the cities for one vehicle to visit
n = len(N) # number of cities
# Create n random points
points = [(0, 0)]
points += [(random.randint(0, 100), random.randint(0, 100)) for i in range(max(N))]
# Dictionary of ALL distances
dist = {(i, j):
math.sqrt(sum((points[i][k]-points[j][k])**2 for k in range(2)))
for i in range(max(N)+1) for j in range(max(N)+1)}
namemap = {}
for i in range(n):
namemap[i] = N[i]
# Dictionary of distances to be passed to model
dists = {}
for i in list(combinations(namemap.keys(),2)):
dists[i] = dist[(namemap[i[0]], namemap[i[1]])]
m = gp.Model()
# Create variables
x = m.addVars(dists.keys(), obj=dists, vtype=GRB.BINARY, name='x')
for i,j in x.keys():
x[j, i] = x[i, j] # edge in opposite direction
# Add degree-2 constraint
m.addConstrs(x.sum(i, '*') == 2 for i in range(n))
# Optimize model
m._vars = x
m.Params.LazyConstraints = 1
m.optimize(subtourelim)
vals = m.getAttr('x', x)
tour = subtour(vals)
assert len(tour) == n
print('')
for i in range(len(tour)):
tour[i] = namemap[i]
print('Original tour: %s' % str(N))
print('Optimal tour: %s' % str(tour))
print('Optimal cost: %g' % m.ObjVal)
print('')0 -
-
-
It looks like the mapping part is not fully correct. It should read
for i in range(len(tour)):
tour[i] = namemap[tour[i]]Best regards,
Jaromił1 -
Post is closed for comments.
Comments
6 comments