VRP model is infeasible
Answered
Good day,
I'm struggling with solving a VRP through MIP. I followed the general VRP model (except for time window constraints) as it is shown in El-Sherbeny & Nasser A (2010)*, but the result is infeasible. What am I doing wrong? I appreciate any help! (below is my code).
*Model from El-Sherbeny, Nasser A. "Vehicle routing with time windows: An overview of exact, heuristic and metaheuristic methods." Journal of King Saud University-Science 22, no. 3 (2010): 123-131.
-------------------------
# Mathematical model for VRP
vrp=Model()
#SETS
Nodes=["DC","Ret1","Ret2","Ret3","Ret4","Ret5","Ret6","Ret7","Ret8","Ret9","Ret10"]
Camiones=["C1","C2","C3"]
#PARAMETERS
capacity=200
demand={"DC":0,
"Ret1":22,
"Ret2":18,
"Ret3":53,
"Ret4":54,
"Ret5":31,
"Ret6":80,
"Ret7":78,
"Ret8":84,
"Ret9":85,
"Ret10":42}
arcs, Dist=multidict({
('DC','Ret1'):6.0828,('DC','Ret2'):3.000000,('DC','Ret3'):7.2801,('DC','Ret4'):8.6023,('DC','Ret5'):6.3246,('DC','Ret6'):7.0711,('DC','Ret7'):10.8167,('DC','Ret8'):2.8284,('DC','Ret9'):5.0000,('DC','Ret10'):8.6023,
('Ret1','DC'):6.0828,('Ret1','Ret2'):9.0554,('Ret1','Ret3'):10.00,('Ret1','Ret4'):6.0828,('Ret1','Ret5'):12.0416,('Ret1','Ret6'):12.53,('Ret1','Ret7'):15.6205,('Ret1','Ret8'):4.1231,('Ret1','Ret9'):5.831,('Ret1','Ret10'):8.0623,
('Ret2','DC'):3.0000,('Ret2','Ret1'):9.0554,('Ret2','Ret3'):7.0711,('Ret2','Ret4'):10.6301,('Ret2','Ret5'):3.6056,('Ret2','Ret6'):5.3852,('Ret2','Ret7'):9.4868,('Ret2','Ret8'):5.3852,('Ret2','Ret9'):7.2111,('Ret2','Ret10'):10.6301,
('Ret3','DC'):7.2801,('Ret3','Ret1'):10.000,('Ret3','Ret2'):7.0711,('Ret3','Ret4'):7.0000,('Ret3','Ret5'):6.4031,('Ret3','Ret6'):12.3693,('Ret3','Ret7'):16.4924,('Ret3','Ret8'):6.4031,('Ret3','Ret9'):12.083,('Ret3','Ret10'):15.6525,
('Ret4','DC'):8.6023,('Ret4','Ret1'):6.0828,('Ret4','Ret2'):10.6301,('Ret4','Ret3'):7.000,('Ret4','Ret5'):12.083,('Ret4','Ret6'):15.6205,('Ret4','Ret7'):19.4165,('Ret4','Ret8'):5.831,('Ret4','Ret9'):11.1803,('Ret4','Ret10'):14,
('Ret5','DC'):6.3246,('Ret5','Ret1'):12.0416,('Ret5','Ret2'):3.6056,('Ret5','Ret3'):6.4031,('Ret5','Ret4'):12.083,('Ret5','Ret6'):7.0711,('Ret5','Ret7'):11,('Ret5','Ret8'):8,('Ret5','Ret9'):10.8167,('Ret5','Ret10'):14.2127,
('Ret6','DC'):7.0711,('Ret6','Ret1'):12.5300,('Ret6','Ret2'):5.3852,('Ret6','Ret3'):12.3693,('Ret6','Ret4'):15.6205,('Ret6','Ret5'):7.0711,('Ret6','Ret7'):4.1231,('Ret6','Ret8'):9.8995,('Ret6','Ret9'):8.0623,('Ret6','Ret10'):10.198,
('Ret7','DC'):10.8167,('Ret7','Ret1'):15.6205,('Ret7','Ret2'):9.4868,('Ret7','Ret3'):16.4924,('Ret7','Ret4'):19.4165,('Ret7','Ret5'):11,('Ret7','Ret6'):4.1231,('Ret7','Ret8'):13.6015,('Ret7','Ret9'):10.2956,('Ret7','Ret10'):11.1803,
('Ret8','DC'):2.8284,('Ret8','Ret1'):4.1231,('Ret8','Ret2'):5.3852,('Ret8','Ret3'):6.4031,('Ret8','Ret4'):5.8310,('Ret8','Ret5'):8.000,('Ret8','Ret6'):9.8995,('Ret8','Ret7'):13.6015,('Ret8','Ret9'):6.0828,('Ret8','Ret10'):9.4868,
('Ret9','DC'):5.0000,('Ret9','Ret1'):5.8310,('Ret9','Ret2'):7.2111,('Ret9','Ret3'):12.083,('Ret9','Ret4'):11.1803,('Ret9','Ret5'):10.8167,('Ret9','Ret6'):8.0623,('Ret9','Ret7'):10.2956,('Ret9','Ret8'):6.0828,('Ret9','Ret10'):3.6056,
('Ret10','DC'):8.6023,('Ret10','Ret1'):8.0623,('Ret10','Ret2'):10.6301,('Ret10','Ret3'):15.6525,('Ret10','Ret4'):14,('Ret10','Ret5'):14.2127,('Ret10','Ret6'):10.198,('Ret10','Ret7'):11.1803,('Ret10','Ret8'):9.4868,('Ret10','Ret9'):3.6056
})
#VARIABLES
x=vrp.addVars(arcs, Camiones,vtype=GRB.BINARY,name="Ruta")
#OBJECTIVE FUNCTION
TotalDist=quicksum(x[i,j,camion]*Dist[i,j]for i,j in arcs for camion in Camiones)
vrp.setObjective(TotalDist,GRB.MINIMIZE)
#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)
TruckCapacity={camion: quicksum(demand[i]*x[i,j,camion] for i,j in arcs) for camion in Camiones} #(2.3)
InFlow={(h,camion):quicksum(x[i,h,camion] for i in Nodes if h!=i) for h in Nodes for camion in Camiones} #(2.5)
OutFlow={(h,camion):quicksum(x[h,j,camion] for j in Nodes if h!=j ) for h in Nodes for camion in Camiones} #(2.5)
DepositOut={camion:quicksum(x['DC',j,camion] for j in Nodes if j!='DC') for camion in Camiones} #(2.4)
DepositIn={camion:quicksum(x[i,'DC',camion] for i in Nodes if i!='DC') for camion in Camiones} #(2.6)
vrp.addConstrs(AllNodesVisited[i]==1 for i in Nodes) #(2.2)
vrp.addConstrs(TruckCapacity[camion]<=capacity for camion in Camiones) #(2.3)
vrp.addConstrs(InFlow[h,camion]==OutFlow[h,camion] for h in Nodes for camion in Camiones) #(2.5)
vrp.addConstrs(DepositOut[camion]==1 for camion in Camiones) #(2.4)
vrp.addConstrs(DepositIn[camion]==1 for camion in Camiones) #(2.6)
#SOLVE
vrp.optimize()
vrp.printAttr('X')
-----------------
Thanks.
-
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
0 -
Thanks, Eli... I'm gonna try out your suggestion.
0 -
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] 1This 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?
0 -
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
0 -
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-toursdef subtourelim(model, where):if where == GRB.Callback.MIPSOL:# make a list of edges selected in the solutionvals = 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 listtour = subtour(selected)if len(tour) < n:# add subtour elimination constr. for every pair of cities in tourmodel.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 subtourdef 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-emptythiscycle = []neighbors = unvisitedwhile 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 cycleDoes it look OK? I am not sure about n. Should it be simply the number of nodes including the warehouse? i.e., len(Nodes).0 -
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
0 -
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_tuplesdef 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):
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.
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 cycleThanks,
Abir
0 -
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 sys
import math
import random
from itertools import 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:
# 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 depot
def 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 argument
if 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 points
dist = {(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 variables
x = m.addVars(dist.keys(), obj=dist, vtype=GRB.BINARY, name='x')
# 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))
# Depot has inbound and outbound flow equal to number of trucks
m.addConstr(x.sum(0, '*') == K)
m.addConstr(x.sum('*', 0) == K)
# Optimize model
m._x = x
m.Params.LazyConstraints = 1
m.optimize(subtourelim)
# 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)
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.py
Usage: vrp.py npoints ntrucks
npoints includes depotAnd an example using \( 8 \) points (including depot) and \( 3 \) trucks:
$ python vrp.py 8 3
Using license file /Users/towle/licenses/gurobi.lic
Changed value of parameter LazyConstraints to 1
Prev: 0 Min: 0 Max: 1 Default: 0
Gurobi Optimizer version 9.0.2 build v9.0.2rc0 (mac64)
Optimize a model with 16 rows, 56 columns and 112 nonzeros
Model fingerprint: 0x7b99499c
Variable 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.00s
Presolved: 16 rows, 56 columns, 112 nonzeros
Variable 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% - 0s
Cutting planes:
Lazy constraints: 3
Explored 0 nodes (17 simplex iterations) in 0.01 seconds
Thread 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 sec
Route for truck 1: 0 -> 4 -> 1 -> 0
Route for truck 2: 0 -> 5 -> 0
Route for truck 3: 0 -> 6 -> 7 -> 2 -> 3 -> 0You can play around with using different types of lazy constraints, strengthening the original formulation with some initial length-2 subtour elimination constraints, etc.
1 -
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 sys
import math
import os
import json
import random
from itertools import permutations
import gurobipy as gp
from gurobipy import GRB
def 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-tours
def 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 depot
def 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 cycle
def 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 data
data = 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 variables
x = m.addVars(dist.keys(), obj=dist, vtype=GRB.BINARY, name='e')
# 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))
# all customers are served
AllNodesVisited={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 trucks
m.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 model
m._x = x
m.Params.LazyConstraints = 1
m.optimize(subtourelim)
# 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)
total_distance = 0
for 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_k
print('Cumulative Distance of all Trucks: {}'.format(total_distance))0 -
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 capacity
u = 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 0
u[0].LB = 0
u[0].UB = 0Let's change the printing at the end to match the format in the OR-Tools page you linked:
# 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)
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 sec
Route for truck 1:
0 Load(0) -> 7 Load(8) -> 1 Load(9) -> 4 Load(13) -> 3 Load(15) -> 0 Load(15)
Route distance: 1552
Route load: 15
Route for truck 2:
0 Load(0) -> 8 Load(8) -> 2 Load(9) -> 6 Load(13) -> 5 Load(15) -> 0 Load(15)
Route distance: 1552
Route load: 15
Route for truck 3:
0 Load(0) -> 12 Load(2) -> 11 Load(3) -> 15 Load(11) -> 13 Load(15) -> 0 Load(15)
Route distance: 1552
Route load: 15
Route for truck 4:
0 Load(0) -> 14 Load(4) -> 16 Load(12) -> 10 Load(14) -> 9 Load(15) -> 0 Load(15)
Route distance: 1552
Route load: 15
Total distance for all routes: 6208.00 -
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 relateddur = {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 modelm = gp.Model()
# Create variablesx = 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 variablep = m.addVar(vtype=GRB.INTEGER, obj=1.0, name="p", lb=1, ub=K)
# Start time of servicet = m.addVars(n, ub=3600, name="t")# Lateness of servicez = m.addVars(n, name="z")# Artificial variables to correct time window upper and lower limitsxa = 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 constraintsm.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 depotm.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:m.addConstr(y[i]+d[j]*x[i,j]-Q*(1-x[i,j])<=y[j])
for i in range(1,n):m.addConstr(y[i] >= d[i])m.addConstr(y[i] <= Q)
# Temporal constraints (8) for customer locationsM = {(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 locationsM = {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-objectiveif params['multiObjective']:m.ModelSense = GRB.MINIMIZEm.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 objectivem.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 modelm._x = xm._y = ym._t = t# m.Params.LazyConstraints = 1m.setParam("TimeLimit", params["TimeLimit"])m.optimize()
# Print optimal routesvals = 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 = 0routes = []customers_served = {i: 0 for i in range(1,n)}for i, tup in enumerate(selected.select(0, '*')):capacity_k = 0dist_k = 0print("\nRoute for truck {}: 0".format(i+1), end='')neighbor = tup[1]prev_node = 0route_k = [0]while neighbor:route_k.append(neighbor)capacity_k += d[neighbor]dist_k += dist[(prev_node, neighbor)]customers_served[neighbor] = i+1print(" -> {} (c={}, d={}, t={:.2f})".format(neighbor, capacity_k, dist_k, t[neighbor].X), end='')prev_node = neighborneighbor = 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_kroutes.append(route_k)print('\nCumulative Distance of all Trucks: {}\n'.format(total_distance))
count, count_late = 0, 0for i in range(1,n):if customers_served[i] == 0:count += 1else: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 += 1if 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_served0 -
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.
0 -
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?
0 -
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.
0
Please sign in to leave a comment.
Comments
14 comments