Need Help with Vehicle Routing Problem with Time Windows (VRPTW) in Python
AnsweredHello,
I've been working on developing a solution for a Vehicle Routing Problem with Time Windows (VRPTW). I found a code example on GitHub which I adapted to work in Python. However, we're encountering issues when we try to run the code with our own distance matrix, tested with a 10x10 matrix – it's not working.
We're using Gurobi for optimization, which should theoretically be able to solve the problem. Additionally, we think there might be some constraints in our model that potentially render the problem infeasible, but we're struggling to pinpoint the exact issue.
Here's a brief overview of our task:
- We have 90 locations, but are currently testing with 10 to simplify.
- Arrival at the first stop (a central depot) is at 10:00 (600 minutes) - this is an estimate.
- Return to the second stop (end location) must be by 12:30 - also an estimate.
- Our matrix indicates time in minutes between locations. Currently, the route can be completed with 9 vehicles, but we aim to reduce this to 8 or fewer.
- Each point must be visited exactly once (departures from location 1 = number of vehicles, arrivals at location 2 = number of vehicles) as all routes must start at location 1 and end at location 2.
- All other locations must be visited exactly once with precisely one departure (arrival at a location followed by departure).
- Departure from each stop is 5 minutes after arrival (service time)
It's worth noting that the code we're basing this on has a capacity constraint, which we don't have in our problem. We've based our work on the following code from GitHub:
https://github.com/ruthmair/vrp/blob/main/vrptw.ipynb
Hoping for some advice, guidelines, or assistance with the code. We're a bit out of our depth here, but any help would be greatly appreciated if you can spare the time. 😊
Here are our current code:
import gurobipy as gp
from gurobipy import GRB
# INSTANCE CONFIGURATION
numCustomers = 9 # 9 customers plus the depot
maxNumVehicles = 6 # Increased number of vehicles
vehicleCapacity = 40 # Increased vehicle capacity
timeHorizon = 900 # Extended time horizon to 3:00 pm (900 minutes)
# Time matrix (depot + 9 customers)
time_matrix = [
[0.00, 17.00, 60.80, 38.08, 23.88, 31.58, 50.05, 50.17, 50.17, 28.20],
[16.82, 0.00, 55.80, 35.98, 21.78, 30.42, 48.83, 48.05, 48.05, 29.15],
[60.65, 55.40, 0.00, 42.52, 47.17, 36.97, 19.80, 19.78, 19.78, 36.50],
[36.30, 36.67, 43.40, 0.00, 19.92, 9.67, 27.57, 27.68, 27.68, 9.27],
[21.35, 21.73, 46.75, 18.95, 0.00, 13.40, 32.10, 31.03, 31.03, 11.83],
[30.13, 31.22, 37.95, 9.35, 14.47, 0.00, 22.12, 22.23, 22.23, 4.08],
[48.47, 48.83, 19.80, 27.45, 32.10, 21.88, 0.00, 0.97, 0.97, 21.43],
[48.73, 49.10, 19.75, 27.72, 32.35, 22.15, 1.08, 0.00, 0.00, 23.08],
[48.73, 49.10, 19.75, 27.72, 32.35, 22.15, 1.08, 0.00, 0.00, 23.08],
[28.20, 28.58, 38.80, 11.00, 11.83, 4.08, 22.97, 23.08, 23.08, 0.00]
]
# Convert the time matrix to a dictionary format suitable for Gurobi
travelTimes = {(i, j): time_matrix[i][j] for i in range(numCustomers + 1) for j in range(numCustomers + 1)}
# Reduced Service time to 2 minutes at each stop
serviceTimes = {i: 2 for i in range(numCustomers + 1)}
# More flexible uniform time window for all customers
timeWindows = {i: (600, timeHorizon) for i in range(numCustomers + 1)}
# Create model
model = gp.Model("VRPTW")
locations = [i for i in range(numCustomers + 1)]
connections = [(i, j) for i in locations for j in locations if i != j]
# Decision variable for routing
x = model.addVars(connections, vtype=GRB.BINARY, name="x")
model.setObjective(gp.quicksum(travelTimes[i, j] * x[i, j] for i, j in connections), GRB.MINIMIZE)
# Constraints
model.addConstrs((x.sum("*", j) == 1 for j in locations if j != 0), name="incoming")
model.addConstrs((x.sum(i, "*") == 1 for i in locations if i != 0), name="outgoing")
model.addConstr(x.sum(0, "*") <= maxNumVehicles, name="maxNumVehicles")
# Time Constraints - Adjusted for more flexibility
y = model.addVars(locations, name="y")
for i in locations:
y[i].LB = timeWindows[i][0]
y[i].UB = timeWindows[i][1]
model.addConstrs((y[i] + serviceTimes[i] + travelTimes[i, j] <= y[j] + (timeHorizon - y[j]) * (1 - x[i, j])
for i, j in connections if i != j), name="timeBigM")
# Solve the model
model.params.Threads = 4
model.optimize()
# Check for infeasibility
if model.status == GRB.INFEASIBLE:
print('Model is infeasible')
model.computeIIS()
model.write("model.ilp")
else:
# Solution Retrieval
usedConnections = [(i, j) for (i, j) in x.keys() if x[i, j].X > 0.5]
nextCustomer = {}
for (i, j) in usedConnections:
if i == 0:
if 0 not in nextCustomer:
nextCustomer[0] = []
nextCustomer[0].append(j)
else:
nextCustomer[i] = j
print(f"Solution contains {len(nextCustomer[0])} routes:")
routeNumber = 0
visitedCustomers = [False] * (numCustomers + 1)
for firstCustomer in nextCustomer[0]:
print(f"Route {routeNumber}: 0 -> ", end="")
vehicleLoad = 0
time = travelTimes[0, firstCustomer]
violatedTimeWindows = False
currentCustomer = firstCustomer
while currentCustomer != 0:
print(f"{currentCustomer} (L:{vehicleLoad}, T:{time}) -> ", end="")
visitedCustomers[currentCustomer] = True
vehicleLoad += 1 # Assuming each customer has a demand of 1 unit
time = max(time, timeWindows[currentCustomer][0])
if time > timeWindows[currentCustomer][1]:
violatedTimeWindows = True
time += (serviceTimes[currentCustomer] + travelTimes[currentCustomer, nextCustomer[currentCustomer]])
currentCustomer = nextCustomer[currentCustomer]
print(f"0 (L:{vehicleLoad}/{vehicleCapacity}, T:{time})")
if vehicleLoad > vehicleCapacity:
print("Vehicle capacity is exceeded!")
if violatedTimeWindows:
print("Time windows are violated!")
routeNumber += 1
print("Unvisited customers: ", end="")
for c in range(1, numCustomers + 1):
if not visitedCustomers[c]:
print(f"{c}, ", end="")
# Free resources
model.dispose()
gp.disposeDefaultEnv()
-
-
Gurobi Staff
Hi Kasoer,
Your timeBigM-constraints are also defined if j is the depot, this makes the model infeasible.
Cheers,
Marika1 -
-
Hi Marika,
Thank you very much for your help.
We tried modifying the code as you suggested. It then worked and gave us a result for our 10x10 matrix
Here is the new code:import gurobipy as gp
from gurobipy import GRB
# INSTANCE CONFIGURATION
numCustomers = 9 # 9 customers plus the depot
maxNumVehicles = 6# Increased number of vehicles
vehicleCapacity = 100 # Increased vehicle capacity
timeHorizon = 900 # Extended time horizon to 3:00 pm (900 minutes)
# Time matrix (depot 89 customers)
time_matrix = [
[0.00, 17.00, 60.80, 38.08, 23.88, 31.58, 50.05, 50.17, 50.17, 28.20],
[16.82, 0.00, 55.80, 35.98, 21.78, 30.42, 48.83, 48.05, 48.05, 29.15],
[60.65, 55.40, 0.00, 42.52, 47.17, 36.97, 19.80, 19.78, 19.78, 36.50],
[36.30, 36.67, 43.40, 0.00, 19.92, 9.67, 27.57, 27.68, 27.68, 9.27],
[21.35, 21.73, 46.75, 18.95, 0.00, 13.40, 32.10, 31.03, 31.03, 11.83],
[30.13, 31.22, 37.95, 9.35, 14.47, 0.00, 22.12, 22.23, 22.23, 4.08],
[48.47, 48.83, 19.80, 27.45, 32.10, 21.88, 0.00, 0.97, 0.97, 21.43],
[48.73, 49.10, 19.75, 27.72, 32.35, 22.15, 1.08, 0.00, 0.00, 23.08],
[48.73, 49.10, 19.75, 27.72, 32.35, 22.15, 1.08, 0.00, 0.00, 23.08],
[28.20, 28.58, 38.80, 11.00, 11.83, 4.08, 22.97, 23.08, 23.08, 0.00]
]
# Convert the time matrix to a dictionary format suitable for Gurobi
travelTimes = {(i, j): time_matrix[i][j] for i in range(numCustomers + 1) for j in range(numCustomers + 1)}
# Reduced Service time to 2 minutes at each stop
serviceTimes = {i: 2 for i in range(numCustomers + 1)}
# More flexible uniform time window for all customers
timeWindows = {i: (600, timeHorizon) for i in range(numCustomers + 1)}
# Create model
model = gp.Model("VRPTW")
locations = [i for i in range(numCustomers + 1)]
connections = [(i, j) for i in locations for j in locations if i != j]
# Decision variable for routing
x = model.addVars(connections, vtype=GRB.BINARY, name="x")
model.setObjective(gp.quicksum(travelTimes[i, j] * x[i, j] for i, j in connections), GRB.MINIMIZE)
# Constraints
model.addConstrs((x.sum("*", j) == 1 for j in locations if j != 0), name="incoming")
model.addConstrs((x.sum(i, "*") == 1 for i in locations if i != 0), name="outgoing")
model.addConstr(x.sum(0, "*") <= maxNumVehicles, name="maxNumVehicles")
# Time Constraints - Adjusted for more flexibility and excluding depot as a destination
y = model.addVars(locations, name="y")
for i in locations:
y[i].LB = timeWindows[i][0]
y[i].UB = timeWindows[i][1]
for i, j in connections:
if j != 0: # Exclude the depot as a destination in these constraints
model.addConstr(y[i] + serviceTimes[i] + travelTimes[i, j] <= y[j] + (timeHorizon - y[j]) * (1 - x[i, j]),
name=f"timeBigM_{i}_{j}")
# Solve the model
model.params.Threads = 4
model.optimize()
# Check for infeasibility
if model.status == GRB.INFEASIBLE:
print('Model is infeasible')
model.computeIIS()
model.write("model.ilp")
else:
print("\nOptimal Solution:")
print("Total travel time:", model.objVal)
# Solution Retrieval
routes = {}
for v in range(maxNumVehicles):
route = [0] # Start from depot
while True:
next_stop = None
for j in locations:
if j != route[-1] and x[route[-1], j].X > 0.5:
next_stop = j
break
if next_stop is None or next_stop == 0:
break # No more stops or back to depot
route.append(next_stop)
if len(route) > 1: # Valid route has more than just the depot
routes[v] = route
# Print the routes
for v, route in routes.items():
print(f"Vehicle {v + 1} route: {' -> '.join(map(str, route))} -> 0")
# Free resources
model.dispose()
gp.disposeDefaultEnv()The result looked like this:
Optimal Solution:
Total travel time: 167.34
Vehicle 1 route: 0 -> 1 -> 2 -> 8 -> 7 -> 6 -> 5 -> 3 -> 9 -> 4 -> 0
Vehicle 2 route: 0 -> 1 -> 2 -> 8 -> 7 -> 6 -> 5 -> 3 -> 9 -> 4 -> 0
Vehicle 3 route: 0 -> 1 -> 2 -> 8 -> 7 -> 6 -> 5 -> 3 -> 9 -> 4 -> 0
Vehicle 4 route: 0 -> 1 -> 2 -> 8 -> 7 -> 6 -> 5 -> 3 -> 9 -> 4 -> 0
Vehicle 5 route: 0 -> 1 -> 2 -> 8 -> 7 -> 6 -> 5 -> 3 -> 9 -> 4 -> 0
Vehicle 6 route: 0 -> 1 -> 2 -> 8 -> 7 -> 6 -> 5 -> 3 -> 9 -> 4 -> 0
Freeing default Gurobi environment
I did not expect the output to look like this. Although I believe it is the correct solution, there is an issue: the routes are displayed as being identical. This is problematic, especially as I have tried extending the code to accommodate my 90x90 matrix. For one instance, the solution appears infeasible, and it seems the code does not remember the locations visited by the first vehicle, which I believe contributes to the infeasibility. I attempted to make the constraints more lenient by allowing the use of 90 vehicles, but it still returned as infeasible.
What I need it to do is remember the locations visited by vehicle one, as subsequent vehicles should not/cannot visit these locations. Each location, apart from the depot, should be visited exactly once.
Do you have any suggestions?0 -
-
-
Hi Kasoer,
I see two issues.
- The printing of the routes is not correct. In your example, only one vehicle is needed. So, there is no route for vehicles 2 to 6. Here is an idea of how to fix the code
# Solution Retrieval
routes = {}
numVeh = 0
for start in range(1, numCustomers+1):
route = [0] # Start from depot
if x[0, start].X < 0.5:
continue
next_stop = start
route.append(start)
while True:
next_stop = None
for j in locations:
if j != route[-1] and x[route[-1], j].X > 0.5:
next_stop = j
break
if next_stop is None or next_stop == 0:
break # No more stops or back to depot
route.append(next_stop)
if len(route) > 1: # Valid route has more than just the depot
routes[numVeh] = route
numVeh += 1
- The timeBigM constraints look strange. I think it should be
y[i] + serviceTimes[i] + travelTimes[i, j] <= y[j] + (timeHorizon) * (1 - x[i, j])
Best regards,
Marika0 -
Please sign in to leave a comment.
Comments
3 comments