Vehicle Routing Problem with Time Windows with Gurobi and Python
AnsweredHi Team Gurobi!
I'm trying to model VRPTW formulation in Python with Gurobi. Any suggestions as to why the model is still infeasible? Would love to have some feedback, especially the subtour elimination part
Consider 3 vehicles of 130 capacity each.
the data will be imported using data.txt:
The datasets for this problem consists of multiple rows; each row defines a delivery task
with the following information:
 LOC_ID: identification number for a location; the first row (with DEMAND = 0)
represent the depot
 XCOORD: x coordinates of the location
 YCOORD: y coordinates of the location
 DEMAND: parcel volume to be delivered at each location (for the depot: zero)
 READYTIME: earliest time for delivery (for the depot: opening time)
 DUETIME: latest time for delivery (for the depot: closing time)
 SERVICETIME: time needed for a delivery at the location
Here is our current code:
import gurobipy as gp
from gurobipy import GRB
import numpy as np
import pandas as pd
file_path = "C:/Users/ahmad/Downloads/data_small.txt"
# Read the data from the text file into a pandas DataFrame
data_df = pd.read_csv(file_path, sep='\s+', header=None)
# Define the column names
data_df.columns = ["LOC_ID", "XCOORD", "YCOORD", "DEMAND", "READYTIME", "DUETIME", "SERVICETIME"]
# Extract the data into lists
locations = [(x, y) for x, y in zip(data_df["XCOORD"], data_df["YCOORD"])]
demands = data_df["DEMAND"].tolist()
ready_times = data_df["READYTIME"].tolist()
due_times = data_df["DUETIME"].tolist()
service_times = data_df["SERVICETIME"].tolist()
num_locations = len(locations)
num_vehicles = 3
vehicle_capacity = 130
depot = 0
# Create a 2D array of distances between each pair of locations
distances = np.hypot(
np.subtract.outer([x for x, y in locations], [x for x, y in locations]),
np.subtract.outer([y for x, y in locations], [y for x, y in locations])
)
# Initialize the model
m = gp.Model('CVRPTW')
# Add variables
x = m.addVars(num_locations, num_locations, num_vehicles, vtype=GRB.BINARY, name="x")
w = m.addVars(num_locations, num_vehicles, vtype=GRB.CONTINUOUS, name="w")
u = m.addVars(num_locations, num_vehicles, vtype=GRB.CONTINUOUS, name="u")
# Set objective: minimize the total distance traveled
m.setObjective(gp.quicksum(distances[i, j] * x[i, j, k] for i in range(num_locations) for j in range(num_locations) for k in range(num_vehicles)), GRB.MINIMIZE)
# Add constraints
# Each location is visited exactly once by exactly one vehicle
for j in range(1, num_locations):
m.addConstr(gp.quicksum(x[i, j, k] for i in range(num_locations) for k in range(num_vehicles) if i != j) == 1)
# Vehicle capacity constraints
for k in range(num_vehicles):
m.addConstr(gp.quicksum(demands[j] * x[i, j, k] for i in range(num_locations) for j in range(num_locations) if i != j) <= vehicle_capacity)
# Time window constraints
for j in range(num_locations):
for k in range(num_vehicles):
m.addConstr(w[j, k] >= ready_times[j])
m.addConstr(w[j, k] <= due_times[j])
# Service start time consistency
for i in range(num_locations):
for j in range(num_locations):
if i != j:
for k in range(num_vehicles):
# When vehicle k travels from i to j, we need to enforce the service time at i and the travel time to j
m.addConstr((w[i, k] + service_times[i] + distances[i, j]) * x[i, j, k] <= w[j, k],
name=f'time_{i}_{j}_{k}')
# Subtour elimination and vehicle load consistency
for i in range(1, num_locations):
for j in range(1, num_locations):
if i != j:
for k in range(num_vehicles):
m.addConstr(u[i, k] + demands[j] <= u[j, k] + (1  x[i, j, k]) * vehicle_capacity)
# Each vehicle leaves and enters the depot
for k in range(num_vehicles):
m.addConstr(gp.quicksum(x[depot, j, k] for j in range(1, num_locations)) == 1)
m.addConstr(gp.quicksum(x[i, depot, k] for i in range(1, num_locations)) == 1)
# Flow conservation at each customer location
for k in range(num_vehicles):
for j in range(1, num_locations): # Exclude the depot
m.addConstr(
gp.quicksum(x[i, j, k] for i in range(num_locations) if i != j) ==
gp.quicksum(x[j, h, k] for h in range(num_locations) if h != j),
name=f'flow_conservation_vehicle{k}_location{j}'
)
# Ensure that the vehicle load at the depot starts at zero (u[0, k] == 0)
# and the vehicle load at each customer location i does not exceed its capacity (u[i, k] <= vehicle_capacity)
for k in range(num_vehicles):
m.addConstr(u[depot, k] == 0, name=f'vehicle_load_depot_{k}')
for i in range(1, num_locations):
m.addConstr(u[i, k] >= demands[i], name=f'vehicle_load_min_{i}_{k}')
m.addConstr(u[i, k] <= vehicle_capacity, name=f'vehicle_load_max_{i}_{k}')
# Solve the model
m.optimize()
# Extract the solution
routes = {}
if m.status == GRB.OPTIMAL:
print('Optimal solution found with objective value:', m.objVal)
for k in range(num_vehicles):
routes[k] = []
for i in range(num_locations):
for j in range(num_locations):
if x[i, j, k].X > 0.5: # Threshold for binary variables
routes[k].append((i, j))
print(f'Vehicle {k}\'s route: {routes[k]}')
else:
print('No optimal solution found. Status code:', m.status)

Hi Ahmad,
This could simply be a case of the number of vehicles, and their capacity, cannot satisfy the demand and the time windows, as opposed to something wrong in your formulation.
The first step would be to review our article here:
How do I determine why my model is infeasible?
In an ideal scenario you will have a small IIS which helps you discover the root cause of the infeasibility. If you suspect there is something wrong with your model it would help to try and construct a smaller problem with a known solution, for which the solver also declares infeasibility.
 Riley
0
Please sign in to leave a comment.
Comments
1 comment