Travelling Salesman with Time Windows

Hello,

I am trying to modify the Milk Collection example and set contraints for farm 10 and 16 thus the distance droven by the driver doesn't excess the limits for respectively 10 and 20 when arriving at these farms. However, that seems that I am doing something wrong, because the cumulative sum (distance) for farm 10 and 16 are respectively 17.4 and 66.6.

Really, I am searching for an example of the solution of the travelling salesman with time windows which should be used for the business application at my work, thus I hope you can help me.

Thanks in advance.

import gurobipy as gp
from gurobipy import GRB
import pandas as pd
import math
from numpy import cumsum
from itertools import combinations,permutations,product

# Create a dictionary to capture the farm coordinates (10 miles) and DistanceMax.
farms, coordinates, DistanceMax  = gp.multidict({
    0: [(0,0), 10000],
    1: [(-3,3), 10000],
    2: [(1,11), 10000],
    3: [(4,7), 10000],
    4: [(-5,9), 10000],
    5: [(-5,-2), 10000],
    6: [(-4,-7), 10000],
    7: [(6,0), 10000],
    8: [(3,-6), 10000],
    9: [(-1,-3), 10000],
    10: [(0,-6), 10],
    11: [(6,4), 10000],
    12: [(2,5), 10000],
    13: [(-2,8), 10000],
    14: [(6,10), 10000],
    15: [(1,8), 10000],
    16: [(-3,1), 20],
    17: [(-6,5), 10000],
    18: [(2,9), 10000],
    19: [(-6,-5), 10000],
    20: [(5,-4), 10000]
})

# List of farm  including the depot
farms = [*range(0,21)]

# Compute pairwise distance matrix
# numpy linalg norm = euclidean n=2

def distance(city1, city2):
    c1 = coordinates[city1]
    c2 = coordinates[city2]
    diff = (c1[0]-c2[0], c1[1]-c2[1])
    return math.sqrt(diff[0]*diff[0]+diff[1]*diff[1])

dist = {(c1, c2): distance(c1, c2) for c1, c2 in product(farms, farms) if c1 != c2}

# Generate the edges list
list_ij = []
for i,j in dist.keys():
    tp = i,j
    list_ij.append(tp)
ij = gp.tuplelist(list_ij)

# Callback - use lazy constraints to eliminate sub-tours
def subtourelim(model, where):
    if where == GRB.Callback.MIPSOL:
        # make a list of edges selected in the solution
        vals = model.cbGetSolution(model._vars)
        # Selecting the farms visited during day 1
        selected = gp.tuplelist((i, j) for i, j in model._vars.keys() if (vals[i, j] > 0.5) )
        # find a cycle in the selected edge list of day 1
        nodes = set()
        for i,j in selected:
            nodes.add(i)
            nodes.add(j)
        tour = subtour(selected)

# Given a tuplelist of edges, find the shortest subtour
def subtour(edges):
    nodes = set()
    for i,j in edges:
        nodes.add(i)
        nodes.add(j)
    unvisited = list(nodes)
    cycle = list(nodes)
    while unvisited:
        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(thiscycle) <= len(cycle):
            cycle = thiscycle # New shortest subtour
    return cycle

# Create the model object m
m = gp.Model('MilkCollection.lp')

# Decision variables:
# Edge variables = 1, if farm 'i' is adjacent to farm 'j'
vars = m.addVars(ij, vtype=GRB.BINARY, name='x')

# Symmetry constraints.
for i,j in dist.keys():
    m.addConstr(vars[j, i] == vars[i, j], name='symmetry')  # edge in opposite direction

# Every day constraints: two edges incident to an every day farm
m.addConstrs((gp.quicksum(vars[tp] for tp in ij.select(i,'*') ) == 2 for i in farms  ), name='edgesConnection')

# Distance constraint
for j in farms:
    if j > 0:
        m.addConstr((gp.quicksum(dist[i, j] * vars[i, j] for i in farms if i != j) <= DistanceMax[j]), name='DistanceMax')

# Objective function: minimize total distance travel
m.setObjective((gp.quicksum(0.5*dist[i,j]*vars[i,j] for i,j in ij)), GRB.MINIMIZE)

# Optimize the model
m.Params.lazyConstraints = 1
m._vars = vars
m.optimize(subtourelim)

# Print optimal routes and distance traveled
print(f"The optimal distance traveled is: {10*round(m.objVal)} miles.")

vals = m.getAttr('X', vars)

edges = gp.tuplelist((i, j) for i, j in vals.keys() if (vals[i, j] > 0.5) )
nodes = set()

for i,j in edges:
    nodes.add(i)
    nodes.add(j)
unvisited = list(nodes)
while unvisited:
    thiscycle = []
    neighbors = unvisited
    while neighbors:
        current = neighbors[0]
        print(f"{current} -> ", end='')
        thiscycle.append(current)
        unvisited.remove(current)
        neighbors = [j for i,j in edges.select(current, '*') if j in unvisited]

# Retrieve solution
vals = m.getAttr('x', vars)
selected = gp.tuplelist((i, j) for i, j in vals.keys() if vals[i, j] > 0.5)

tour = subtour(selected)
assert len(tour) == len(farms)
FromLocation = tour
FinalLocation = tour[1:]+[0]
TourDataframe = pd.DataFrame({'FromLocation':FromLocation, 'FinalLocation':FinalLocation})

for index, row in TourDataframe.iterrows():
  distance = dist[row.FromLocation, row.FinalLocation]
  TourDataframe.at[index, 'Distance'] = distance

TourDataframe['CumulativeDistance'] = cumsum(TourDataframe.Distance)

print(TourDataframe)

Official comment

Simranjit Kaur

Gurobi Staff

June 23, 2025 02:58 Edited

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 why not try our AI Gurobot?.

Silke Horn

March 16, 2021 09:37

Hi Niels,

Thanks for posting your code. Your distance constraint seems to me as though you are summing up the edges adjacent to farm j when you actually want those that come before it. Is that correct?

At this point of your code, we don't know about any order of the edges though. I would suggest that the easiest way to include the distances in your code would be via additional lazy constraints. It shouldn't be too difficult to modify the subtour method so that it also sums up the distances driven until a certain node and compares this with the node's maximum allowed distance.

Also note that the original TSP example code is working with an undirected graph while for your application the direction of the edges/tour matters a lot. So you should also adapt the code to work on a directed graph.

Silke

March 16, 2021 10:47 Edited

I have given this some more thought and my previous suggestions might a bit too easy... While you can easily spot violations of the distance constraint in the lazy constraint code, it is not that easy to come up with a good constraint to cut off this solution. In fact, I would assume that the additional distance constraints may make it difficult to even find any feasible solution. So lazy constraints are probably not a wise approach here at all. On the contrary, I would expect that the solver would then just blindly walk through the entire branch-and-bound tree until it accidentally hits a solution that is not cut off by a lazy constraint. I hope that makes sense.

I guess you will need a different approach -- and a much bigger model/graph. For example, with discrete time steps, introduce one node for each combination of farm and time step, then model the path in this bigger graph.

Niels Strande

March 16, 2021 11:47 Edited

Hi Silke,

Thank you for the replies. I am very beginner at Gurobi and currently experimentering with Gurobi for eventually adoption at my work. Could you like to refer to some example in the Gurobi examples library which I can use as a start point for the modelling?

Niels

March 18, 2021 10:04

To get started with Gurobi, you can find many resources under the Resources tab on the Gurobi website. The code and modelling examples, demos, and webinars and especially anything related to scheduling could be interesting for you.

As a commercial user, you can also submit a support request if you have any specific questions.

Travelling Salesman with Time Windows

Comments

Didn't find what you were looking for?