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('Last-mile 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.")
0
-
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 [6e-02, 6e+00] Bounds range [1e+00, 1e+00] RHS range [3e-01, 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 Relax-and-lift: 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 -
-
-
Gurobi Staff
Hi Thuy Anh Pham,
Vehicle routing problems are in general very difficult combinatorial optimization problems. State-of-the-art solution approaches involve a highly-tuned and sophisticated branch-price-and-cut 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 Big-M 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 general-purpose 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 branch-and-bound 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 Big-M constraints.
- Additionally, adding connectivity or capacity cuts dynamically in callbacks (branch-and-cut 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