• Gurobi Staff

Hi Marcela,

I think the problem is with this set of constraints:

AllNodesVisited={i: quicksum(x[i,j,camion] for j in Nodes if i!=j for camion in Camiones) for i in Nodes} #(2.2)vrp.addConstrs(AllNodesVisited[i]==1 for i in Nodes) #(2.2)

which states that each node must be visited by exactly one truck. However, because each of the three trucks must pass through the depot, we should exclude the depot from these constraints:

AllNodesVisited={i: quicksum(x[i,j,camion] for j in Nodes if i!=j for camion in Camiones) for i in Nodes if i != 'DC'} #(2.2)vrp.addConstrs(AllNodesVisited[i]==1 for i in Nodes if i != 'DC') #(2.2)

I hope this helps!

Eli

Thanks, Eli... I'm gonna try out your suggestion.

The above code creates disconnected routes. For example, for C1,

Ruta[DC,Ret8,C1] 1
Ruta[Ret3,Ret4,C1] 1
Ruta[Ret4,Ret3,C1] 1
Ruta[Ret8,DC,C1] 1

This is because there is no constraint that says that any local loop is not permitted. Is there a different way to take care of these?

• Gurobi Staff

Hi Abir,

Yes, you will have to somehow disallow subtours like this. For the VRP, this can be done by adding lazy constraints inside of a callback. The lazy constraints would not be much different from those used in the TSP example.

Eli

Got it. I am assuming I have apply the lazy callback for each camione individually. Towards that I have modified the original functions in the following way:

# 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)
for cam in Camiones:
# print(t, nbTrucks)
selected = gp.tuplelist((i, j, k) for i, j, k in model._vars.keys()
if vals[i, j, k] > 0.5 and k == cam)
# print(selected)
# find the shortest cycle in the selected edge list
tour = subtour(selected)
if len(tour) < n:
# add subtour elimination constr. for every pair of cities in tour
model.cbLazy(gp.quicksum(model._vars[i, j, k]
for i, j, k in combinations(tour, 2))
<= len(tour)-1)

# Given a tuplelist of edges, find the shortest subtour
def subtour(edges):
unvisited = list(range(n))
cycle = range(n+1)  # initial length has 1 more city
# print(edges)
while unvisited:  # true if list is non-empty
thiscycle = []
neighbors = unvisited
while neighbors:
current = neighbors[0]
# print('current: ', current)
thiscycle.append(current)
unvisited.remove(current)
neighbors = [j for i, j, k in edges.select(current, '*')
if j in unvisited]
# print('neighbors: ', neighbors)
if len(cycle) > len(thiscycle):
cycle = thiscycle
# print(cycle)
return cycle
Does it look OK? I am not sure about n. Should it be simply the number of nodes including the warehouse? i.e.,  len(Nodes).

• Gurobi Staff

Hi Abir,

That is a start. You cannot use a fixed value of $$\texttt{n}$$ here, because the constraint you want to add depends on the length of the subtour you find.

In the $$\texttt{subtour}$$ function, you should extract the list of nodes visited by the truck and look for the shortest cycle among those that does not include the $$\texttt{DC}$$. The best way to do this is probably to not include the truck index when you construct the $$\texttt{selected}$$ list in the $$\texttt{subtourelim}$$ callback function (i.e., $$\texttt{selected}$$ should be a list of length-2 tuples). Then, $$\texttt{combinations(tour, 2)}$$ should contain only city-city pairs, so the $$\texttt{for}$$ loop in the lazy constraint should be $$\texttt{for i, j in combinations(tour, 2)}$$ (no $$\texttt{k}$$). The $$\texttt{k}$$-index in the constraint is really the local $$\texttt{cam}$$ variable. However, the lazy constraint you derive could apply to every truck.

By the way, there are lots of different ways to model the VRP and many different flavors of lazy constraints to add. I'm sure if you poke around the literature you could find some more efficient ways to solve it.

Thanks,

Eli

Hi Eli,

Here is how I modified the functions to enforce the subtour elimination:

I have a helper function to renumber the nodes.

def renumber(tuples):    numbers = []    for t in tuples:        numbers.append(t[0])        numbers.append(t[1])    numbers = list(set(numbers))    renum, renum2 = dict(), dict()    k = 0    for n in numbers:        renum[k] = n  # sorted keys        renum2[n] = k # original keys        k += 1    new_tuples = [(renum2[x], renum2[y]) for x,y in tuples]    return renum, renum2, new_tuples
def subtourelim(model, where):    if where == GRB.Callback.MIPSOL:        # make a list of edges selected in the solution        vals = model.cbGetSolution(model._vars)        for t in range(nbTrucks):            tuples = gp.tuplelist((i, j) for i, j, k in model._vars.keys()                                    if vals[i, j, k] > 0.5 and k == t)            d1, _, newtuples = renumber(tuples)            selected = gp.tuplelist(newtuples)            n = len(d1)            # find the shortest cycle in the selected edge list            tour = subtour(selected, n)            if len(tour) < n:                new_tour = [d1[x] for x in tour]  # map back to original node number                # add subtour elimination constr. for every pair of cities in tour                for t in range(nbTrucks):                    if len(new_tour) == 2:                        model.cbLazy(gp.quicksum(model._vars[i, j, t] + model._vars[j, i, t] for i, j in combinations(new_tour, 2)) <= len(tour)-1)                    else:                        model.cbLazy(gp.quicksum(model._vars[i, j, t] for i, j in combinations(new_tour, 2)) <= len(tour)-1)
def subtour(edges, n):    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]        # assuming 0 will be there for all vehicles and 0 means Depot        if len(cycle) > len(thiscycle) and 0 not in thiscycle:            cycle = thiscycle    return cycle
However, I am still not getting the subtours eliminated from the final solution although they are clearly captured by the above functions. Any suggestion would be extremely useful.

Thanks,

Abir

• Gurobi Staff

I don't see why you need the $$\texttt{renumber}$$ function here. Can't you calculate the subtour with the original values?

I think it would be easier to use a VRP formulation that doesn't explicitly include an index for each truck. Instead, the trucks are implicitly represented as cycles emanating from the depot location. With this approach, your variables are indexed like $$x_{ij}$$ instead of $$x_{ijk}$$.

It's not so hard to modify the tsp.py example to get a working version of this. Below is a $$\texttt{vrp.py}$$ example, made from the following changes to $$\texttt{tsp.py}$$:

• Use directed arcs instead of edges
• Modify flow balance constraints for all points (special constraints for depot, which is node $$0$$)
• Change $$\texttt{subtour}$$ function to only look for subtours over nodes not connected to the depot
import sysimport mathimport randomfrom itertools import permutationsimport gurobipy as gpfrom gurobipy import GRB# Callback - use lazy constraints to eliminate sub-toursdef subtourelim(model, where):    if where == GRB.Callback.MIPSOL:        # make a list of edges selected in the solution        vals = model.cbGetSolution(model._x)        selected = gp.tuplelist((i, j) for i, j in model._x.keys()                                if vals[i, j] > 0.5)        # find the shortest cycle in the selected edge list        tour = subtour(selected)        if len(tour) < n:            # add subtour elimination constr. for every pair of cities in tour            model.cbLazy(gp.quicksum(model._x[i, j]                                     for i, j in permutations(tour, 2))                         <= len(tour)-1)# Given a tuplelist of edges, find the shortest subtour not containing depotdef subtour(edges):    unvisited = list(range(1, n))    cycle = range(n+1)  # initial length has 1 more city    # First, remove all nodes connected to depot    depot_connected = [j for i, j in edges.select(0, '*')]    while depot_connected:        current = depot_connected.pop()        unvisited.remove(current)        neighbors = [j for i, j in edges.select(current, '*')                     if j in unvisited and j != 0]        depot_connected += neighbors    # Now, find subtour using tsp.py code    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(cycle) > len(thiscycle):            cycle = thiscycle    return cycle# Parse argumentif len(sys.argv) < 3:    print('Usage: vrp.py npoints ntrucks')    print('{:7}npoints includes depot'.format(""))    sys.exit(1)n = int(sys.argv[1])K = int(sys.argv[2])if K > n:    print("npoints must be at least as large as ntrucks")    sys.exit(1)# Create n random points# Depot is node 0, located at (50, 50)random.seed(1)points = [(50, 50)]points += [(random.randint(0, 100), random.randint(0, 100))           for i in range(n-1)]# Dictionary of Euclidean distance between each pair of pointsdist = {(i, j):        math.sqrt(sum((points[i][k]-points[j][k])**2 for k in range(2)))        for i in range(n) for j in range(n) if i != j}m = gp.Model()# Create variablesx = m.addVars(dist.keys(), obj=dist, vtype=GRB.BINARY, name='x')# Inbound and outbound flow is always 1, except for depotm.addConstrs(x.sum(i, '*') == 1 for i in range(1, n))m.addConstrs(x.sum('*', i) == 1 for i in range(1, n))# Depot has inbound and outbound flow equal to number of trucksm.addConstr(x.sum(0, '*') == K)m.addConstr(x.sum('*', 0) == K)# Optimize modelm._x = xm.Params.LazyConstraints = 1m.optimize(subtourelim)# Print optimal routesvals = m.getAttr('X', x)selected = gp.tuplelist((i, j) for i, j in vals.keys() if vals[i, j] > 0.5)for i, tup in enumerate(selected.select(0, '*')):    print("Route for truck {}: 0".format(i+1), end='')    neighbor = tup[1]    while neighbor:        print(" -> {}".format(neighbor), end='')        neighbor = selected.select(neighbor, '*')[0][1]    print(" -> 0")

To use it:

$python vrp.pyUsage: vrp.py npoints ntrucks npoints includes depot And an example using $$8$$ points (including depot) and $$3$$ trucks: $ python vrp.py 8 3Using license file /Users/towle/licenses/gurobi.licChanged value of parameter LazyConstraints to 1   Prev: 0  Min: 0  Max: 1  Default: 0Gurobi Optimizer version 9.0.2 build v9.0.2rc0 (mac64)Optimize a model with 16 rows, 56 columns and 112 nonzerosModel fingerprint: 0x7b99499cVariable types: 0 continuous, 56 integer (56 binary)Coefficient statistics:  Matrix range     [1e+00, 1e+00]  Objective range  [1e+01, 1e+02]  Bounds range     [1e+00, 1e+00]  RHS range        [1e+00, 3e+00]Presolve time: 0.00sPresolved: 16 rows, 56 columns, 112 nonzerosVariable types: 0 continuous, 56 integer (56 binary)Root relaxation: objective 3.323386e+02, 16 iterations, 0.00 seconds    Nodes    |    Current Node    |     Objective Bounds      |     Work Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time*    0     0               0     349.0383665  349.03837  0.00%     -    0sCutting planes:  Lazy constraints: 3Explored 0 nodes (17 simplex iterations) in 0.01 secondsThread count was 4 (of 4 available processors)Solution count 1: 349.038 Optimal solution found (tolerance 1.00e-04)Best objective 3.490383665458e+02, best bound 3.490383665458e+02, gap 0.0000%User-callback calls 44, time in user-callback 0.00 secRoute for truck 1: 0 -> 4 -> 1 -> 0Route for truck 2: 0 -> 5 -> 0Route for truck 3: 0 -> 6 -> 7 -> 2 -> 3 -> 0

You can play around with using different types of lazy constraints, strengthening the original formulation with some initial length-2 subtour elimination constraints, etc.

Hi Eli,

The reason I was going for an explicit index for vehicles is to add one more objective of minimizing the number of vehicles used. Also, for some reason the subtour function was not working with original values.

In any case, I have modified your code to have the capacity constraint. However, on the Google sample data the result is 9516, which is nowhere close to the minimum distance, i.e., 6208. The code is shown below:

import sysimport mathimport osimport jsonimport randomfrom itertools import permutationsimport gurobipy as gpfrom gurobipy import GRBdef make_edges(ll):    tuples = []    for ii in range(1, len(ll)):        tuples.append((ll[ii-1], ll[ii]))    return tuples# Callback - use lazy constraints to eliminate sub-toursdef subtourelim(model, where):    if where == GRB.Callback.MIPSOL:        # make a list of edges selected in the solution        vals = model.cbGetSolution(model._x)        selected = gp.tuplelist((i, j) for i, j in model._x.keys()                                if vals[i, j] > 0.5)        # add the capacity constraints        for i, tup in enumerate(selected.select(0, '*')):            capacity_k = 0            nodes_k = [0]            neighbor = tup[1]            while neighbor:                capacity_k += Q[neighbor]                nodes_k.append(neighbor)                neighbor = selected.select(neighbor, '*')[0][1]            if capacity_k > capacity:                model.cbLazy(gp.quicksum(Q[j]*model._x[i,j] for i, j in make_edges(nodes_k)) <= capacity)        # find the shortest cycle in the selected edge list        tour = subtour(selected)        if len(tour) < n:            # add subtour elimination constr. for every pair of cities in tour            model.cbLazy(gp.quicksum(model._x[i, j]                                     for i, j in permutations(tour, 2))                         <= len(tour)-1)# Given a tuplelist of edges, find the shortest subtour not containing depotdef subtour(edges):    unvisited = list(range(1, n))    cycle = range(n+1)  # initial length has 1 more city    # First, remove all nodes connected to depot    depot_connected = [j for i, j in edges.select(0, '*') if j != 0]    # print('Depot_connected:', depot_connected)    while depot_connected:        current = depot_connected.pop()        # print('Current:', current)        # print('Unvisited:', unvisited)        unvisited.remove(current)        neighbors = [j for i, j in edges.select(current, '*')                     if j in unvisited and j != 0]        depot_connected += neighbors    # Now, find subtour using tsp.py code    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(cycle) > len(thiscycle):            cycle = thiscycle    return cycledef create_data_model():    """Stores the data for the problem."""    data = {}    data['distance_matrix'] = [        [            0, 548, 776, 696, 582, 274, 502, 194, 308, 194, 536, 502, 388, 354,            468, 776, 662        ],        [            548, 0, 684, 308, 194, 502, 730, 354, 696, 742, 1084, 594, 480, 674,            1016, 868, 1210        ],        [            776, 684, 0, 992, 878, 502, 274, 810, 468, 742, 400, 1278, 1164,            1130, 788, 1552, 754        ],        [            696, 308, 992, 0, 114, 650, 878, 502, 844, 890, 1232, 514, 628, 822,            1164, 560, 1358        ],        [            582, 194, 878, 114, 0, 536, 764, 388, 730, 776, 1118, 400, 514, 708,            1050, 674, 1244        ],        [            274, 502, 502, 650, 536, 0, 228, 308, 194, 240, 582, 776, 662, 628,            514, 1050, 708        ],        [            502, 730, 274, 878, 764, 228, 0, 536, 194, 468, 354, 1004, 890, 856,            514, 1278, 480        ],        [            194, 354, 810, 502, 388, 308, 536, 0, 342, 388, 730, 468, 354, 320,            662, 742, 856        ],        [            308, 696, 468, 844, 730, 194, 194, 342, 0, 274, 388, 810, 696, 662,            320, 1084, 514        ],        [            194, 742, 742, 890, 776, 240, 468, 388, 274, 0, 342, 536, 422, 388,            274, 810, 468        ],        [            536, 1084, 400, 1232, 1118, 582, 354, 730, 388, 342, 0, 878, 764,            730, 388, 1152, 354        ],        [            502, 594, 1278, 514, 400, 776, 1004, 468, 810, 536, 878, 0, 114,            308, 650, 274, 844        ],        [            388, 480, 1164, 628, 514, 662, 890, 354, 696, 422, 764, 114, 0, 194,            536, 388, 730        ],        [            354, 674, 1130, 822, 708, 628, 856, 320, 662, 388, 730, 308, 194, 0,            342, 422, 536        ],        [            468, 1016, 788, 1164, 1050, 514, 514, 662, 320, 274, 388, 650, 536,            342, 0, 764, 194        ],        [            776, 868, 1552, 560, 674, 1050, 1278, 742, 1084, 810, 1152, 274,            388, 422, 764, 0, 798        ],        [            662, 1210, 754, 1358, 1244, 708, 480, 856, 514, 468, 354, 844, 730,            536, 194, 798, 0        ],    ]    data['num_vehicles'] = 4    data['depot'] = 0    data['demands'] = [0, 1, 1, 2, 4, 2, 4, 8, 8, 1, 2, 1, 2, 4, 4, 8, 8]    data['vehicle_capacities'] = [15, 15, 15, 15]    return datadata = create_data_model()distList = data['distance_matrix']Q = data['demands']capacities = data['vehicle_capacities']capacity = capacities[0]print(Q)n = len(distList)dist = {(i, j): distList[i][j] for i in range(n) for j in range(n)}K = data['num_vehicles']if K > n:    print("npoints must be at least as large as ntrucks")    sys.exit(1)print(f"Total {n-1} customers, {K} trucks with {capacity} capacity")m = gp.Model()# Create variablesx = m.addVars(dist.keys(), obj=dist, vtype=GRB.BINARY, name='e')# Inbound and outbound flow is always 1, except for depotm.addConstrs(x.sum(i, '*') == 1 for i in range(1, n))m.addConstrs(x.sum('*', i) == 1 for i in range(1, n))# all customers are servedAllNodesVisited={i: gp.quicksum(x[j,i] for j in range(n) if i!=j) for i in range(1,n)} # (2.2)m.addConstrs(AllNodesVisited[i]==1 for i in range(1,n))# Depot has inbound and outbound flow equal to number of trucksm.addConstr(x.sum(0, '*') <= K)m.addConstr(x.sum('*', 0) <= K)totaldist = gp.quicksum(x[i,j]*dist[i,j] for i in range(n) for j in range(n))m.setObjective(totaldist, GRB.MINIMIZE)# Optimize modelm._x = xm.Params.LazyConstraints = 1m.optimize(subtourelim)# Print optimal routesvals = m.getAttr('X', x)selected = gp.tuplelist((i, j) for i, j in vals.keys() if vals[i, j] > 0.5)total_distance = 0for i, tup in enumerate(selected.select(0, '*')):    capacity_k = 0    dist_k = 0    print("Route for truck {}: 0".format(i+1), end='')    neighbor = tup[1]    prev_node = 0    while neighbor:        capacity_k += Q[neighbor]        dist_k += dist[(prev_node, neighbor)]        print(" -> {} ({}, {})".format(neighbor, capacity_k, dist_k), end='')        prev_node = neighbor        neighbor = selected.select(neighbor, '*')[0][1]    dist_k += dist[(prev_node, 0)]    print(" -> 0 ({}, {})".format(capacity_k, dist_k))    total_distance += dist_kprint('Cumulative Distance of all Trucks: {}'.format(total_distance))

• Gurobi Staff

If you want a variable number of trucks, you can make $$\texttt{K}$$ a variable with a lower bound of $$1$$. Also, you don't need those $$\texttt{AllNodesVisited}$$ constraints. Those are equivalent to the inbound flow constraints directly above them.

I wouldn't try to handle the capacity constraints inside of a callback. That seems messy. You can just add capacity constraints directly to the model. In particular, define a new variable $$u \in \mathbb{R}^n_+$$ that is equal to the cumulative load of a truck that reaches node $$j$$. Then, add constraints that enforce $$u_j \geq u_i + Q_j$$ if $$x_{ij} = 1$$. If $$M$$ is your truck capacity, this can be done with the following constraints:

\begin{align*} u_j \geq u_i + Q_j - M(1 - x_{ij}) \quad \forall i \in N,\ j \in N \setminus \{ 0 \} \textrm{ s.t. } i \neq j.\end{align*}

• If $$x_{ij} = 1$$, then we have $$u_j \geq u_i + Q_j$$, as desired.
• If $$x_{ij} = 0$$, we instead have $$u_j \geq u_i + Q_j - M$$. Since $$u_i \leq M$$ and  $$u_j \geq Q_j$$, the constraint doesn't do anything.

Here is a quick implementation in Python:

# Track cumulative demand at each node; cannot exceed capacityu = m.addVars(n, ub=capacity, name='u')pairs = [(i, j) for i in range(n) for j in range(1, n) if i != j]m.addConstrs((u[j] >= Q[j] + u[i] - (1 - x[i, j]) * capacity             for (i, j) in pairs), 'demand')# Depot cumulative demand is always 0u[0].LB = 0u[0].UB = 0

Let's change the printing at the end to match the format in the OR-Tools page you linked:

# Print optimal routesvals = m.getAttr('X', x)selected = gp.tuplelist((i, j) for i, j in vals.keys() if vals[i, j] > 0.5)for i, tup in enumerate(selected.select(0, '*')):    print("\nRoute for truck {}:\n 0 Load(0)".format(i+1), end='')    neighbor = tup[1]    truck_dist = distList[0][neighbor]    truck_load = Q[neighbor]    while neighbor:        print(" -> {} Load({})".format(neighbor, truck_load), end='')        next_neighbor = selected.select(neighbor, '*')[0][1]        truck_dist += distList[neighbor][next_neighbor]        truck_load += Q[next_neighbor]        neighbor = next_neighbor    print(" -> 0 Load({})".format(truck_load))    print("Route distance: {}".format(truck_dist))    print("Route load: {}".format(truck_load))print("\nTotal distance for all routes: {}".format(m.ObjVal))

The results after solving with the OR-Tools example data:

...Optimal solution found (tolerance 1.00e-04)Best objective 6.208000000000e+03, best bound 6.208000000000e+03, gap 0.0000%User-callback calls 4104, time in user-callback 0.01 secRoute for truck 1: 0 Load(0) -> 7 Load(8) -> 1 Load(9) -> 4 Load(13) -> 3 Load(15) -> 0 Load(15)Route distance: 1552Route load: 15Route for truck 2: 0 Load(0) -> 8 Load(8) -> 2 Load(9) -> 6 Load(13) -> 5 Load(15) -> 0 Load(15)Route distance: 1552Route load: 15Route for truck 3: 0 Load(0) -> 12 Load(2) -> 11 Load(3) -> 15 Load(11) -> 13 Load(15) -> 0 Load(15)Route distance: 1552Route load: 15Route for truck 4: 0 Load(0) -> 14 Load(4) -> 16 Load(12) -> 10 Load(14) -> 9 Load(15) -> 0 Load(15)Route distance: 1552Route load: 15Total distance for all routes: 6208.0

Thank you so much Eli for the solution. It is working as expected.

Next, I am including time window constraints for CVRPTW problem. However, I have coded a slightly different approach of savings maximization (paper) instead of distance minimization. The function is shown below. I am not getting optimal result although the constraints are not violated. Would you please take a look and let me know if this is the right way of implementing multiple objectives?

def cvrp_3(data, params):

c = data['distance_matrix']
tau = data['time_matrix']
n = len(c)
q = data['demands']
d = {i: q[i] for i in range(1,n)}
capacities = data['vehicle_capacities']
Q = capacities[0]

dist = {(i, j): c[i][j] for i in range(n) for j in range(n)}
s = {(i, j): c[i][0] + c[0][j] - c[i][j] for i in range(1,n) for j in range(1,n)}
K = data['num_vehicles']

# Time window related
dur = {j : data['unloading_time'][j] for j in range(n)}
tStart = {j : data['time_windows'][j][0] for j in range(n)}
tEnd = {j : data['time_windows'][j][1] for j in range(n)}
tDue = {j : tEnd[j] for j in range(n)}
priority = {j : 1 for j in range(n)}

# print(max([data['time_windows'][j][0] for j in range(n)]), max([data['time_windows'][j][1] for j in range(n)]))

# create the model
m = gp.Model()

# Create variables
x = m.addVars(dist.keys(), obj=dist, vtype=GRB.BINARY, name='e')
y = m.addVars(d, vtype=GRB.CONTINUOUS, name='y')
# add truck number as a variable
p = m.addVar(vtype=GRB.INTEGER, obj=1.0, name="p", lb=1, ub=K)

# Start time of service
t = m.addVars(n, ub=3600, name="t")

# Lateness of service
z = m.addVars(n, name="z")

# Artificial variables to correct time window upper and lower limits
xa = m.addVars(n, name="xa")
xb = m.addVars(n, name="xb")

# Track cumulative demand at each node; cannot exceed capacity
# u = m.addVars(n, ub=capacities, name='u')
# pairs = [(i, j) for i in range(n) for j in range(1, n) if i != j]
# m.addConstrs((u[j] >= Q[j] + u[i] - (1 - x[i, j]) * Q for (i, j) in pairs), 'demand')

# Depot cumulative demand is always 0
# u[0].LB = 0
# u[0].UB = 0

m.update()

# Add the constraints
m.addConstr(gp.quicksum(x[0,j] for j in range(1,n)) == p)  # (9)

inflow={j:gp.quicksum(x[i,j] for i in range(n) if i!=j) for j in range(1,n)} # (10)
m.addConstrs(inflow[j]==1 for j in range(1,n))

outflow={i:gp.quicksum(x[i,j] for j in range(1,n) if j!=i) for i in range(1,n)}  # (11)
m.addConstrs(outflow[i]<=1 for i in range(1,n))

# Inbound and outbound flow is always 1, except for depot
m.addConstrs(x.sum(i, '*') == 1 for i in range(1, n))
m.addConstrs(x.sum('*', i) == 1 for i in range(1, n))

for i in range(1,n):
for j in range(1,n):
if i != j:

for i in range(1,n):

# Temporal constraints (8) for customer locations
M = {(i,j) : 3600 + dur[i] + dist[i, j] for i in range(1,n) for j in range(1,n)}
m.addConstrs((t[j] >= t[i] + dur[i] + dist[i, j] - M[i,j]*(1 - x[i,j]) for i in range(1,n) for j in range(1,n)), name="tempoCustomer")

# Temporal constraints (8) for depot locations
M = {j : 3600 + dist[0, j] for j in range(1,n)}
m.addConstrs((t[j] >= t[0] + dist[0, j] - M[j]*(1 - x.sum(0, j)) for j in range(1,n)), name="tempoDepot")

# Time window constraints (9 and 10)
m.addConstrs((t[j] + xa[j] >= tStart[j] for j in range(1,n)), name="timeWinA")
m.addConstrs((t[j] - xb[j] <= tEnd[j] for j in range(1,n)), name="timeWinB")

# Lateness constraint (11)
m.addConstrs((z[j] >= t[j] + dur[j] - tDue[j] for j in range(1,n)), name="lateness")

totalsaving = gp.quicksum(x[i,j]*s[i,j] for i in range(1,n) for j in range(1,n))
# Multi-objective
if params['multiObjective']:
m.ModelSense = GRB.MINIMIZE
m.setObjectiveN(-totalsaving, index=0, priority=3)
m.setObjectiveN(z.prod(priority) + gp.quicksum( 0.01 * 6100 * priority[j] * (xa[j] + xb[j]) for j in range(1,n)), index=1, priority=2)
m.setObjectiveN(p, index=2, priority=1)
else:
# single objective
m.setObjective(-totalsaving + z.prod(priority) + gp.quicksum( 0.01 * 6100 * priority[j] * (xa[j] + xb[j]) for j in range(1,n)) + p, GRB.MINIMIZE)

# Optimize model
m._x = x
m._y = y
m._t = t
# m.Params.LazyConstraints = 1
m.setParam("TimeLimit", params["TimeLimit"])
m.optimize()

# Print optimal routes
vals = m.getAttr('X', x)
selected = gp.tuplelist((i, j) for i, j in vals.keys() if vals[i, j] > 0.5)

# print(selected)
total_distance = 0
routes = []
customers_served = {i: 0 for i in range(1,n)}
for i, tup in enumerate(selected.select(0, '*')):
capacity_k = 0
dist_k = 0
print("\nRoute for truck {}: 0".format(i+1), end='')
neighbor = tup[1]
prev_node = 0
route_k = [0]
while neighbor:
route_k.append(neighbor)
capacity_k += d[neighbor]
dist_k += dist[(prev_node, neighbor)]
customers_served[neighbor] = i+1
print(" -> {} (c={}, d={}, t={:.2f})".format(neighbor, capacity_k, dist_k, t[neighbor].X), end='')
prev_node = neighbor
neighbor = selected.select(neighbor, '*')[0][1]
dist_k += dist[(prev_node, 0)]
route_k.append(0)
print(" -> 0 (c={}, d={})".format(capacity_k, dist_k))
total_distance += dist_k
routes.append(route_k)
print('\nCumulative Distance of all Trucks: {}\n'.format(total_distance))

count, count_late = 0, 0
for i in range(1,n):
if customers_served[i] == 0:
count += 1
else:
jobStr = "Customer-{}: Truck-{}, arrived at {:.2f}, Time-window - ({} - {})".format(i, customers_served[i], t[i].X, tStart[i], tEnd[i])
# print(f"Customer-{i}: Truck-{customers_served[i]}, arrived at {t[i].X}, Time-window - ({tStart[i]} to {tEnd[i]})")
if z[i].X > 1e-6:
jobStr += " {:.2f} minutes late.".format(z[i].X)
count_late += 1
if xa[i].X > 1e-6:
jobStr += " Start time corrected by {:.2f} minutes.".format(xa[i].X)
if xb[i].X > 1e-6:
jobStr += " End time corrected by {:.2f} minutes.".format(xb[i].X)
print(jobStr)
if count_late > 0:
print('\n{} customers are not served in time'.format(count_late))
if count > 0:
print('\n{} customers not served'.format(count))
print([i for i in range(1,n) if customers_served[i]==0])

return routes, total_distance, customers_served

• Gurobi Staff

Your approach for defining multiple objectives looks good to me. That said, there may be a problem in one of the objective expressions or in one of the constraints. You could try writing an LP file to disk with $$\texttt{m.write('cvrptw.lp')}$$. You can visually inspect this LP file to make sure that your objectives and constraints are written correctly.

Hey Eli, I love your tsp.py-like vrp.py programming. I think tsp.py should address the word "symmetric" somewhere to make it understandable. In the vrp.py case, you have already defined the model for an assymetric distance matrix. So, it is cool! I even think it can be displayed in Gurobi Examples.

One question I have though is that is there any easy to implement heuristic to get an upper bound for the vrp.py? Often, we would like to constrain the program by TimeLimit, and if we don't feed in a start, we may end up finding no solution within the TimeLimit. Are there any Gurobi settings that enforce finding a solution to VRP quickly?

Another question is that is it possible to have a distance constraint for each vehicle with this formulation or do we have to extend the indices with k?

• Gurobi Staff

I haven't studied the VRP much, but I'm sure there are many different heuristics in the literature. If you just want to tell Gurobi to focus more on finding feasible solutions, you could try setting MIPFocus=1, increasing the value of the Heuristics parameter, or enabling specific heuristics with NoRelHeurTime, ZeroObjNodes, MinRelNodes, and PumpPasses.

If all of the vehicles can travel the same maximum distance $$d_{\texttt{max}}$$, you could introduce a new variable for each node $$i$$ that is equal to the cumulative distance traveled by a vehicle when it reaches node $$i$$. Then, add constraints that this cumulative distance is no greater than $$d_{\texttt{max}}$$ less the distance from node $$i$$ to the depot.

If the vehicles have different distance constraints, it would probably be easiest to extend the formulation to include individual trucks, like you suggested.