Vehicle Routing Problem on Python using Gurobi
AnsweredHi,
I have written python code to solve a Vehicle Routing Problem Pickup and Delivery with Time Windows using Gurobi. It takes, however, really long time to solve and I still received no solution.
I based it on below mathematical model:
Objective function is to minimize the total travel distance.
Could you please help if something is wrong with the code?
import gurobipy as gp
from gurobipy import GRB
import pandas as pd
import numpy as np
# Read input data from Excel file
input_df = pd.read_excel('C:/Users/49172/Downloads/Vytal/TAP01.xlsx')
# Read distance matrix from file
distance_df = pd.read_excel('C:/Users/49172/Downloads/Vytal/distance_matrix_1803.xlsx')
# Read travel time matrix from file
traveltime_df = pd.read_excel('C:/Users/49172/Downloads/Vytal/traveltime_matrix_2003.xlsx')
# Read binary matrix from file
binary_df = pd.read_excel('C:/Users/49172/Downloads/Vytal/binary_matrix_1803.xlsx')
# Define the model
model = gp.Model('Lastmile Delivery and Pickup Routing Problem')
# Set of vehicles
V = range(1, 8)
#V = ['mo', 'tu', 'we', 'th', 'fr', 'sa', 'su']
# Set of customers
n = len(input_df)  1 #Subtract 1 for column name #Depot is indexed in files as both 0 and n
C = input_df["Store_ID"].tolist()
M = 99999
Q = 500
# Start time for each customer
st = dict(zip(C, input_df["StartTime"].tolist()))
# End time for each customer
et = dict(zip(C, input_df["EndTime"].tolist()))
# Pickup demand of each customer
pu = dict(zip(C, input_df["Pickup"].tolist()))
# Delivery demand of each customer
de = dict(zip(C, input_df["Delivery"].tolist()))
# Distance between each pair of customers
d = {(i, j): distance_df.loc[i, j] for i in C for j in C}
# Travel time between each pair of customers
tv = {(i, j): traveltime_df.loc[i, j] for i in C for j in C}
# Binary matrix indicating if a customer can be visited by a particular type of vehicle
b = {(i, v): binary_df.loc[i, v] for i in C for v in V}
# Create decision variables
x = {} # Binary variable for whether a route includes an arc
for v in V:
for i in C:
for j in C:
x[i,j,v] = model.addVar(vtype=gp.GRB.BINARY, name=f"x_{i}_{j}_{v}")
y = {} # Integer variable for load of products leaving depot
for v in V:
y[0, v] = model.addVar(vtype=GRB.INTEGER, name=f"y_0_{v}")
z = {} # Integer variable for load of products after servicing j
for j in range(1,n+1):
z[j] = model.addVar(vtype=gp.GRB.INTEGER, name=f"z_{j}")
t = {} # Continuous variable for arrival time at the location
for v in V:
for i in C:
t[i, v] = model.addVar(vtype=gp.GRB.CONTINUOUS, name=f"t_{i}_{v}")
# Define objective function
model.setObjective(gp.quicksum(d[i, j] * x[i, j, v] for i in C for j in C if i !=j for v in V), GRB.MINIMIZE)
# Define constraints
# Each customer must be visited exactly once
for j in C:
model.addConstr(gp.quicksum(x[i, j, v] for v in V for i in C if i != j) == 1, name='firstConst')
# Vehicle leaves node that it enters
for v in V:
for h in range(1,n+1):
model.addConstr(gp.quicksum(x[i,h,v] for i in C) == gp.quicksum(x[h,j,v] for j in C), name= '2nd_b')
# Each vehicle must start and end at the depot
for v in V:
model.addConstr(gp.quicksum(x[0, j, v] for j in range(1, n+1)) == gp.quicksum(x[i, 0, v] for i in range(1, n+1)), name=f"depot_{v}")
for v in V:
model.addConstr(y[0, v] == gp.quicksum(de[j]* x[i, j, v] for i in C for j in range(1,n+1) if j!=i), name='4Const')
for v in V:
for j in range(1,n+1):
model.addConstr(z[j] >= y[0, v]  de[j] + pu[j]  M* (1  x[0, j, v]), name='5Const')
for i in range(1,n+1):
for j in range(1,n+1):
model.addConstr(z[j] >= z[i]  de[j] + pu[j]  M*(1  gp.quicksum(x[i, j, v] for v in V)), name='6Const')
for v in V:
model.addConstr(y[0,v] <= Q, name='7Const')
for v in V:
for j in range(1,n+1):
model.addConstr(z[j] <= Q + M*(1  gp.quicksum(x[i, j, v] for i in C)), name='8Const')
for i in C:
for j in C:
model.addConstr(t[i,v] + 1/144 + tv[i, j]  M*(1  gp.quicksum(x[i, j, v] for v in V)) <= t[j, v], name='9Const')
for i in C:
for v in V:
model.addConstr(t[i, v] >= st[i], name='10Const')
model.addConstr(t[i, v] <= et[i], name='11Const')
model.setParam(GRB.Param.TimeLimit, 100.0)
# Solve the model
model.optimize()
if model.status == GRB.INFEASIBLE:
model.computeIIS()
model.write("infeasible.lp")
print('\nThe following constraint(s) cannot be satisfied:')
for c in model.getConstrs():
if c.IISConstr:
print('%s' % c.constrName)
# Check if the model is solved to optimality
if model.status == gp.GRB.OPTIMAL:
# Print the optimal objective value
print(f"Optimal objective value: {model.objVal}")
# Print the solution for each variable
for v in model.getVars():
print(f"{v.varName} = {v.x}")
else:
print("No solution found.")

I have set the time limit to 100 seconds and the report is following:
Set parameter TimeLimit to value 100 Gurobi Optimizer version 9.5.0 build v9.5.0rc5 (win64) Thread count: 4 physical cores, 8 logical processors, using up to 8 threads Optimize a model with 3301 rows, 7893 columns and 55427 nonzeros Model fingerprint: 0xa97ff95d Variable types: 231 continuous, 7662 integer (7623 binary) Coefficient statistics: Matrix range [1e+00, 1e+05] Objective range [6e02, 6e+00] Bounds range [1e+00, 1e+00] RHS range [3e01, 1e+05] Presolve removed 624 rows and 744 columns Presolve time: 0.16s Presolved: 2677 rows, 7149 columns, 51719 nonzeros Variable types: 33 continuous, 7116 integer (7077 binary) Root relaxation: objective 1.631622e+01, 647 iterations, 0.06 seconds (0.03 work units) Nodes  Current Node  Objective Bounds  Work Expl Unexpl  Obj Depth IntInf  Incumbent BestBd Gap  It/Node Time 0 0 16.31622 0 71  16.31622   0s 0 0 18.02170 0 94  18.02170   0s 0 0 18.03463 0 103  18.03463   1s 0 0 18.74425 0 122  18.74425   1s 0 0 18.74926 0 129  18.74926   1s 0 0 18.75114 0 130  18.75114   1s 0 0 18.75749 0 137  18.75749   1s 0 0 18.89854 0 111  18.89854   2s 0 0 18.99712 0 92  18.99712   2s 0 0 18.99712 0 84  18.99712   2s 0 0 19.08739 0 76  19.08739   2s 0 0 19.09536 0 76  19.09536   2s 0 0 19.09536 0 73  19.09536   2s 0 0 19.09536 0 76  19.09536   2s 0 0 19.10368 0 83  19.10368   3s 0 0 19.10368 0 60  19.10368   3s 0 2 19.12233 0 58  19.12233   4s 11 16 19.12463 4 68  19.12233  69.1 5s 456 444 infeasible 70  19.12233  45.3 10s 1080 955 25.30838 89 60  19.12233  44.0 15s 1086 959 19.41407 15 101  19.15521  43.7 20s 1090 962 19.65402 26 140  19.23144  43.6 25s 1097 966 19.49847 14 84  19.38201  43.3 30s 1105 972 32.45293 88 104  19.50556  43.0 35s 1112 976 33.33750 74 172  19.54711  42.7 40s 1119 981 31.49252 132 96  19.54958  42.4 45s 1128 987 19.55443 11 194  19.55443  42.1 50s 1136 992 19.55739 18 126  19.55739  41.8 55s 1145 999 20.27492 24 117  19.55739  74.4 61s 1150 1003 35.98030 139 168  19.55739  74.0 65s 1155 1006 19.97654 17 120  19.57099  73.7 70s 1161 1010 28.11875 104 180  19.58055  73.3 76s 1170 1016 19.58919 25 185  19.58919  72.8 80s 1177 1021 34.82520 97 116  19.60531  72.3 88s 1178 1024 19.60531 19 112  19.60531  50.0 93s 1184 1030 19.66622 21 146  19.66622  50.8 96s 1251 1079 19.80732 28 101  19.78275  58.5 100s Cutting planes: Learned: 4 Gomory: 30 Implied bound: 6 Clique: 1 MIR: 12 StrongCG: 1 Flow cover: 72 Zero half: 26 RLT: 103 Relaxandlift: 10 Explored 1259 nodes (124060 simplex iterations) in 100.04 seconds (55.85 work units) Thread count was 8 (of 8 available processors) Solution count 0 Time limit reached Best objective , best bound 1.978274782307e+01, gap  No solution found.
0 
Hi Thuy Anh Pham,
Vehicle routing problems are in general very difficult combinatorial optimization problems. Stateoftheart solution approaches involve a highlytuned and sophisticated branchpriceandcut algorithm. You can read more details for example in the famous book:
Toth, Paolo, and Daniele Vigo, eds. Vehicle routing: problems, methods, and applications. Society for industrial and applied mathematics, 2014.
There are plenty of different ways to model the same vehicle routing problem, in particular to deal with subtour elimination and resource constraints like capacity or time. Your model is the classical BigM model and this is known to provide very weak dual bounds.
Additionally, the combinatorial structure of the problem makes it difficult to find feasible solutions by modern generalpurpose MIP solvers like Gurobi (but this applies to others as well).
There are different ways to at least improve the chance to find good solutions for this model: You could create a solution with some simple heuristic and hand it over to Gurobi as MIP start. Then, the builtin improvement heuristics might be able to find better solutions more easily.
 You could try the optional NoRel heuristic in Gurobi that often helps for highly combinatorial problems like scheduling, routing, or other graph problems. You can enable it with parameter NoRelHeurTime defining the time in seconds spent for this heuristic before the actual branchandbound part starts.
Additionally, I would recommend to work on your formulation:
 You could push the limits further by using stronger modeling techniques like commodity flows to avoid the BigM constraints.
 Additionally, adding connectivity or capacity cuts dynamically in callbacks (branchandcut approach) usually help to improve performance.
 If the fleet is homogeneous, you can get rid of the vehicle index k in your variables (by using appropriate modeling techniques from the literature).
Best regards,
Mario1
Please sign in to leave a comment.
Comments
2 comments