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:

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