Interfacing with Gurobi: A Travelling Salesman Solver¶
You can see this app running online at: Travelling Salesman Solver App Online
What this app does:
- Solves the travelling salesman problem for up to 30 locations
This app aims to demonstrate:
- How to interface to Gurobi using Gurobi’s python interface
- How to render a route in KML
Gurobi worked example
- The code that solves the optimisation problem within this app is taken from an online Gurobi example and is used with permission.
Setup Instructions¶
Before you can run this app, you will need to install Gurobi, and obtain a trial license. Details can be found on the Gurobi website
Next, use the app name 'tropofy_gurobi_travelling_salesman'
to quickstart as in Running and Debugging Tropofy Apps
Full Source¶
"""
Authors: www.tropofy.com and www.gurobi.com
Copyright 2015 Tropofy Pty Ltd, all rights reserved.
Copyright 2013, Gurobi Optimization, Inc.
This source file (where not indicated as under the copyright of Gurobi)
is part of Tropofy and governed by the Tropofy terms of service
available at: http://www.tropofy.com/terms_of_service.html
Parts of the formulation provided by Gurobi have been modified.
The original example is in the Gurobi installation in the example file dietmodel.py
Used with permission.
This source file is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
or FITNESS FOR A PARTICULAR PURPOSE. See the license files for details.
"""
import os
from sqlalchemy.types import Integer, Text, Float
from sqlalchemy.schema import Column, ForeignKeyConstraint, UniqueConstraint
from sqlalchemy.orm import relationship, backref
from simplekml import Kml, Style, IconStyle, Icon, LineStyle
import random
from math import radians, cos, sin, asin, sqrt
import pkg_resources
import gurobipy
from tropofy.app import AppWithDataSets, Step, StepGroup
from tropofy.widgets import ExecuteFunction, SimpleGrid, KMLMap
from tropofy.database.tropofy_orm import DataSetMixin
class TravellingSalesManStop(DataSetMixin):
name = Column(Text, nullable=False)
latitude = Column(Float, nullable=False)
longitude = Column(Float, nullable=False)
ordinal = Column(Integer, nullable=False)
def __init__(self, name, latitude, longitude, ordinal):
self.name = name
self.latitude = latitude
self.longitude = longitude
self.ordinal = ordinal
@classmethod
def get_table_args(cls):
return (UniqueConstraint('data_set_id', 'name'),)
class KMLMapForStops(KMLMap):
def get_kml(self, app_session):
kml = Kml()
stop_style = Style(iconstyle=IconStyle(scale=0.8, icon=Icon(href='https://maps.google.com/mapfiles/kml/paddle/blu-circle-lv.png')))
stops = kml.newfolder(name="Stops")
for s in [stops.newpoint(name=s.name, coords=[(s.longitude, s.latitude)]) for s in app_session.data_set.query(TravellingSalesManStop).all()]:
s.style = stop_style
return kml.kml()
class KMLMapForRoute(KMLMap):
def get_kml(self, app_session):
kml = Kml()
def LongLat(l):
return (l.longitude, l.latitude)
mylocstyle = Style(iconstyle=IconStyle(scale=0.8, icon=Icon(href='https://maps.google.com/mapfiles/kml/paddle/blu-circle-lv.png')))
LocsFolder = kml.newfolder(name="Locations")
for p in [LocsFolder.newpoint(name=loc.name, coords=[LongLat(loc)]) for loc in app_session.data_set.query(TravellingSalesManStop).all()]:
p.style = mylocstyle
mylinestyle = Style(linestyle=LineStyle(color='FFFF00FF', width=4))
LinesFolder = kml.newfolder(name="Route")
route = app_session.data_set.query(TravellingSalesManStop).all()
route.sort(key=lambda x: x.ordinal)
for i, _ in enumerate(route):
if i > 0:
line = LinesFolder.newlinestring(name='line', coords=[LongLat(route[i]), LongLat(route[i - 1])])
else:
line = LinesFolder.newlinestring(name='line', coords=[LongLat(route[-1]), LongLat(route[i])])
line.style = mylinestyle
return kml.kml()
class ExecuteGurobiSolver(ExecuteFunction):
def get_button_text(self, app_session):
return "Solve Travelling Salesman Problem"
def execute_function(self, app_session):
if len(app_session.data_set.query(TravellingSalesManStop).all()) > 30:
app_session.task_manager.send_progress_message("More than 30 stops exceeds to license limits of the free trial version of Gurobi")
else:
formulate_and_solve_travelling_salesman_problem(app_session)
class GurobiRosteringApp(AppWithDataSets):
def get_name(self):
return 'Travelling Salesman Solver'
def get_examples(self):
return {"Demo 30 Stop random data set": load_random_data}
def get_static_content_path(self, app_session):
return pkg_resources.resource_filename('te_gurobi_travelling_salesman', 'static')
def get_gui(self):
step_group1 = StepGroup(name='Enter your stops')
step_group1.add_step(Step(name='Enter your stops', widgets=[SimpleGrid(TravellingSalesManStop)]))
step_group1.add_step(Step(name='Review the stops', widgets=[KMLMapForStops()]))
step_group2 = StepGroup(name='Solve')
step_group2.add_step(Step(name='Solve travelling salesman problem using Gurobi', widgets=[ExecuteGurobiSolver()]))
step_group3 = StepGroup(name='View the Solution')
step_group3.add_step(Step(
name='Optimal Route',
widgets=[
{"widget": SimpleGrid(TravellingSalesManStop), "cols": 6},
{"widget": KMLMapForRoute(), "cols": 6}
]
))
return [step_group1, step_group2, step_group3]
def get_icon_url(self):
return "/{}/static/{}/travelling_salesman.png".format(
self.url_name,
self.get_app_version(),
)
def load_random_data(app_session):
# Extents roughly cover central Chicago
top = 41.936060
bottom = 41.833714
left = -87.801123
right = -87.652781
locations = []
for i in range(0, 30):
locations.append(TravellingSalesManStop(str(i), random.uniform(bottom, top), random.uniform(right, left), 0))
app_session.data_set.add_all(locations)
def haversine(lon1, lat1, lon2, lat2):
"""
Calculate the great circle distance between two points on the earth (specified in decimal degrees)
From http://stackoverflow.com/questions/4913349/haversine-formula-in-python-bearing-and-distance-between-two-gps-points
"""
# convert decimal degrees to radians
lon1, lat1, lon2, lat2 = map(radians, [lon1, lat1, lon2, lat2])
# haversine formula
dlon = lon2 - lon1
dlat = lat2 - lat1
a = sin(dlat / 2) ** 2 + cos(lat1) * cos(lat2) * sin(dlon / 2) ** 2
c = 2 * asin(sqrt(a))
km = 6367 * c
return km
def formulate_and_solve_travelling_salesman_problem(app_session):
# Copyright 2013, Gurobi Optimization, Inc.
# Solve a traveling salesman problem on a randomly generated set of
# points using lazy constraints. The base MIP model only includes
# 'degree-2' constraints, requiring each node to have exactly
# two incident edges. Solutions to this model may contain subtours -
# tours that don't visit every city. The lazy constraint callback
# adds new constraints to cut them off.
# Callback - use lazy constraints to eliminate sub-tours
def subtourelim(model, where):
if where == gurobipy.GRB.callback.MIPSOL:
selected = []
# make a list of edges selected in the solution
for i in range(n):
sol = model.cbGetSolution([model._vars[i, j] for j in range(n)])
selected += [(i, j) for j in range(n) if sol[j] > 0.5]
# find the shortest cycle in the selected edge list
tour = subtour(selected)
if len(tour) < n:
# add a subtour elimination constraint
expr = 0
for i in range(len(tour)):
for j in range(i + 1, len(tour)):
expr += model._vars[tour[i], tour[j]]
model.cbLazy(expr <= len(tour) - 1)
# Euclidean distance between two points
def distance(points, i, j):
return haversine(points[i].longitude, points[i].latitude, points[j].longitude, points[j].latitude)
# Given a list of edges, finds the shortest subtour
def subtour(edges):
visited = [False] * n
cycles = []
lengths = []
selected = [[] for i in range(n)]
for x, y in edges:
selected[x].append(y)
while True:
current = visited.index(False)
thiscycle = [current]
while True:
visited[current] = True
neighbors = [x for x in selected[current] if not visited[x]]
if len(neighbors) == 0:
break
current = neighbors[0]
thiscycle.append(current)
cycles.append(thiscycle)
lengths.append(len(thiscycle))
if sum(lengths) == n:
break
return cycles[lengths.index(min(lengths))]
points = app_session.data_set.query(TravellingSalesManStop).all()
n = len(points)
m = gurobipy.Model()
# Create variables
variables = {}
for i in range(n):
for j in range(i + 1):
variables[i, j] = m.addVar(obj=distance(points, i, j), vtype=gurobipy.GRB.BINARY, name='e' + str(i) + '_' + str(j))
variables[j, i] = variables[i, j]
m.update()
# Add degree-2 constraint, and forbid loops
for i in range(n):
m.addConstr(gurobipy.quicksum(variables[i, j] for j in range(n)) == 2)
variables[i, i].ub = 0
m.update()
# Optimize model
m._vars = variables
m.params.LazyConstraints = 1
m.optimize(subtourelim)
solution = m.getAttr('x', variables)
selected = [(i, j) for i in range(n) for j in range(n) if solution[i, j] > 0.5]
assert len(subtour(selected)) == n
app_session.task_manager.send_progress_message('Optimal tour: %s' % subtour(selected))
app_session.task_manager.send_progress_message('Optimal distance: %s' % m.objVal)
counter = 0
for l in subtour(selected):
points[l].ordinal = counter
counter += 1