Vehicle Routing Problem with Time Windows with Gurobi and Python

Hi 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 2-D 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)