I am solving a simple capacitated VRP with CP-SAT solver with 6 vehicles and 17 nodes. The issue that is happening is that the solution ends up using all 6 vehicles, when clearly a smaller number of vehicles could have sufficed (considering each node demand, and sufficient vehicle capacity). How can I use the right number of vehicles that would minimize the objective function
Please find below the code:
from ortools.sat.python import cp_model as cp
locations = range(17)
dist_matrix = [
[
0, 548, 776, 696, 582, 274, 502, 194, 308, 194, 536, 502, 388, 354,
468, 776, 662
],
[
548, 0, 684, 308, 194, 502, 730, 354, 696, 742, 1084, 594, 480, 674,
1016, 868, 1210
],
[
776, 684, 0, 992, 878, 502, 274, 810, 468, 742, 400, 1278, 1164,
1130, 788, 1552, 754
],
[
696, 308, 992, 0, 114, 650, 878, 502, 844, 890, 1232, 514, 628, 822,
1164, 560, 1358
],
[
582, 194, 878, 114, 0, 536, 764, 388, 730, 776, 1118, 400, 514, 708,
1050, 674, 1244
],
[
274, 502, 502, 650, 536, 0, 228, 308, 194, 240, 582, 776, 662, 628,
514, 1050, 708
],
[
502, 730, 274, 878, 764, 228, 0, 536, 194, 468, 354, 1004, 890, 856,
514, 1278, 480
],
[
194, 354, 810, 502, 388, 308, 536, 0, 342, 388, 730, 468, 354, 320,
662, 742, 856
],
[
308, 696, 468, 844, 730, 194, 194, 342, 0, 274, 388, 810, 696, 662,
320, 1084, 514
],
[
194, 742, 742, 890, 776, 240, 468, 388, 274, 0, 342, 536, 422, 388,
274, 810, 468
],
[
536, 1084, 400, 1232, 1118, 582, 354, 730, 388, 342, 0, 878, 764,
730, 388, 1152, 354
],
[
502, 594, 1278, 514, 400, 776, 1004, 468, 810, 536, 878, 0, 114,
308, 650, 274, 844
],
[
388, 480, 1164, 628, 514, 662, 890, 354, 696, 422, 764, 114, 0, 194,
536, 388, 730
],
[
354, 674, 1130, 822, 708, 628, 856, 320, 662, 388, 730, 308, 194, 0,
342, 422, 536
],
[
468, 1016, 788, 1164, 1050, 514, 514, 662, 320, 274, 388, 650, 536,
342, 0, 764, 194
],
[
776, 868, 1552, 560, 674, 1050, 1278, 742, 1084, 810, 1152, 274,
388, 422, 764, 0, 798
],
[
662, 1210, 754, 1358, 1244, 708, 480, 856, 514, 468, 354, 844, 730,
536, 194, 798, 0
],
]
distance_matrix = {}
for i in locations:
for j in locations:
distance_matrix[(i, j)] = dist_matrix[i][j]
num_vehicles = 6
demands = [0, 1, 1, 2, 4, 2, 4, 8, 8, 1, 2, 1, 2, 4, 4, 8, 8]
vehicle_capacity = 30
depot = 0
model = cp.CpModel()
# Build decision variables :
dv = {}
for i in locations:
for j in locations:
for k in range(num_vehicles):
dv[(i, j, k)] = model.NewBoolVar("from_%i_to_%i_vehc_%i" % (i, j, k))
# Objective function
total_distance_route = sum(dv[i, j, k] * distance_matrix[i, j] for i in locations for j in locations for k in range(num_vehicles))
model.Minimize(total_distance_route)
# Capacity constraint
for k in range(num_vehicles):
model.Add(sum(dv[i, j, k] * demands[j] for i in locations for j in locations[1:]) <= vehicle_capacity)
# no travel from a node to itself
for k in range(num_vehicles):
model.Add(dv[0, 0, k] == 0)
# which nodes were covered in a tour by a vehicle
dv_vehc_node_bin = {}
for k in range(num_vehicles):
lst = {}
for i in locations[1:]:
lst[i] = model.NewBoolVar("vehc_%i_node_%i" % (k, i))
dv_vehc_node_bin[k] = lst
# constraint: each vertex is in exactly one tour
for v in locations[1:]:
model.add(sum(dv_vehc_node_bin[i][v] for i in range(num_vehicles)) == 1)
# sub-tour elimination
arcs_dict = {}
for k in range(num_vehicles):
arcs_list = []
for i in locations:
for j in locations:
if (i == j) & (i == 0):
arcs_list.append([0, 0, dv[0, 0, k]])
elif (i == j) & (i > 0):
arcs_list.append([i, i, ~dv_vehc_node_bin[k][i]])
else:
arcs_list.append([i, j, dv[(i, j, k)]])
arcs_dict[k] = arcs_list
for a in arcs_dict:
model.AddCircuit(arcs_dict[a])
# solve the model
solver = cp.CpSolver()
solver.parameters.num_search_workers = 8
solver.parameters.log_search_progress = True
solver.log_callback = print
status = solver.Solve(model)
The solution I get is as follows. No sub-tours, but using all available vehicles.
0
[(0, 12), (1, 0), (3, 4), (4, 1), (11, 15), (12, 11), (15, 3)]
1
[(0, 14), (2, 6), (6, 8), (8, 0), (10, 2), (14, 16), (16, 10)]
2
[(0, 13), (13, 0)]
3
[(0, 5), (5, 0)]
4
[(0, 7), (7, 0)]
5
[(0, 9), (9, 0)]
I am solving a simple capacitated VRP with CP-SAT solver with 6 vehicles and 17 nodes. The issue that is happening is that the solution ends up using all 6 vehicles, when clearly a smaller number of vehicles could have sufficed (considering each node demand, and sufficient vehicle capacity). How can I use the right number of vehicles that would minimize the objective function
Please find below the code:
from ortools.sat.python import cp_model as cp
locations = range(17)
dist_matrix = [
[
0, 548, 776, 696, 582, 274, 502, 194, 308, 194, 536, 502, 388, 354,
468, 776, 662
],
[
548, 0, 684, 308, 194, 502, 730, 354, 696, 742, 1084, 594, 480, 674,
1016, 868, 1210
],
[
776, 684, 0, 992, 878, 502, 274, 810, 468, 742, 400, 1278, 1164,
1130, 788, 1552, 754
],
[
696, 308, 992, 0, 114, 650, 878, 502, 844, 890, 1232, 514, 628, 822,
1164, 560, 1358
],
[
582, 194, 878, 114, 0, 536, 764, 388, 730, 776, 1118, 400, 514, 708,
1050, 674, 1244
],
[
274, 502, 502, 650, 536, 0, 228, 308, 194, 240, 582, 776, 662, 628,
514, 1050, 708
],
[
502, 730, 274, 878, 764, 228, 0, 536, 194, 468, 354, 1004, 890, 856,
514, 1278, 480
],
[
194, 354, 810, 502, 388, 308, 536, 0, 342, 388, 730, 468, 354, 320,
662, 742, 856
],
[
308, 696, 468, 844, 730, 194, 194, 342, 0, 274, 388, 810, 696, 662,
320, 1084, 514
],
[
194, 742, 742, 890, 776, 240, 468, 388, 274, 0, 342, 536, 422, 388,
274, 810, 468
],
[
536, 1084, 400, 1232, 1118, 582, 354, 730, 388, 342, 0, 878, 764,
730, 388, 1152, 354
],
[
502, 594, 1278, 514, 400, 776, 1004, 468, 810, 536, 878, 0, 114,
308, 650, 274, 844
],
[
388, 480, 1164, 628, 514, 662, 890, 354, 696, 422, 764, 114, 0, 194,
536, 388, 730
],
[
354, 674, 1130, 822, 708, 628, 856, 320, 662, 388, 730, 308, 194, 0,
342, 422, 536
],
[
468, 1016, 788, 1164, 1050, 514, 514, 662, 320, 274, 388, 650, 536,
342, 0, 764, 194
],
[
776, 868, 1552, 560, 674, 1050, 1278, 742, 1084, 810, 1152, 274,
388, 422, 764, 0, 798
],
[
662, 1210, 754, 1358, 1244, 708, 480, 856, 514, 468, 354, 844, 730,
536, 194, 798, 0
],
]
distance_matrix = {}
for i in locations:
for j in locations:
distance_matrix[(i, j)] = dist_matrix[i][j]
num_vehicles = 6
demands = [0, 1, 1, 2, 4, 2, 4, 8, 8, 1, 2, 1, 2, 4, 4, 8, 8]
vehicle_capacity = 30
depot = 0
model = cp.CpModel()
# Build decision variables :
dv = {}
for i in locations:
for j in locations:
for k in range(num_vehicles):
dv[(i, j, k)] = model.NewBoolVar("from_%i_to_%i_vehc_%i" % (i, j, k))
# Objective function
total_distance_route = sum(dv[i, j, k] * distance_matrix[i, j] for i in locations for j in locations for k in range(num_vehicles))
model.Minimize(total_distance_route)
# Capacity constraint
for k in range(num_vehicles):
model.Add(sum(dv[i, j, k] * demands[j] for i in locations for j in locations[1:]) <= vehicle_capacity)
# no travel from a node to itself
for k in range(num_vehicles):
model.Add(dv[0, 0, k] == 0)
# which nodes were covered in a tour by a vehicle
dv_vehc_node_bin = {}
for k in range(num_vehicles):
lst = {}
for i in locations[1:]:
lst[i] = model.NewBoolVar("vehc_%i_node_%i" % (k, i))
dv_vehc_node_bin[k] = lst
# constraint: each vertex is in exactly one tour
for v in locations[1:]:
model.add(sum(dv_vehc_node_bin[i][v] for i in range(num_vehicles)) == 1)
# sub-tour elimination
arcs_dict = {}
for k in range(num_vehicles):
arcs_list = []
for i in locations:
for j in locations:
if (i == j) & (i == 0):
arcs_list.append([0, 0, dv[0, 0, k]])
elif (i == j) & (i > 0):
arcs_list.append([i, i, ~dv_vehc_node_bin[k][i]])
else:
arcs_list.append([i, j, dv[(i, j, k)]])
arcs_dict[k] = arcs_list
for a in arcs_dict:
model.AddCircuit(arcs_dict[a])
# solve the model
solver = cp.CpSolver()
solver.parameters.num_search_workers = 8
solver.parameters.log_search_progress = True
solver.log_callback = print
status = solver.Solve(model)
The solution I get is as follows. No sub-tours, but using all available vehicles.
0
[(0, 12), (1, 0), (3, 4), (4, 1), (11, 15), (12, 11), (15, 3)]
1
[(0, 14), (2, 6), (6, 8), (8, 0), (10, 2), (14, 16), (16, 10)]
2
[(0, 13), (13, 0)]
3
[(0, 5), (5, 0)]
4
[(0, 7), (7, 0)]
5
[(0, 9), (9, 0)]
Thanks for the relevant comments. I understood where the problem was. At one place, I have stated:
# no travel from a node to itself
for k in range(num_vehicles):
model.Add(dv[0, 0, k] == 0)
What this was doing was making sure that every vehicle is leaving the depot (node 0), which was not desirable.
Now, with the above correction (removing the above code block), and minimising the number of vehicle used (to check if the logic is working), I get the correct result. In addition, I have added symmetry breaking constraints - since vehicles are indistinguishable.
Below is the revised code for capacitated VRP with CP-SAT solver. I like the flexibility that CP-SAT solver gives to the user.
from ortools.sat.python import cp_model as cp
locations = range(17)
dist_matrix = [
[
0, 548, 776, 696, 582, 274, 502, 194, 308, 194, 536, 502, 388, 354,
468, 776, 662
],
[
548, 0, 684, 308, 194, 502, 730, 354, 696, 742, 1084, 594, 480, 674,
1016, 868, 1210
],
[
776, 684, 0, 992, 878, 502, 274, 810, 468, 742, 400, 1278, 1164,
1130, 788, 1552, 754
],
[
696, 308, 992, 0, 114, 650, 878, 502, 844, 890, 1232, 514, 628, 822,
1164, 560, 1358
],
[
582, 194, 878, 114, 0, 536, 764, 388, 730, 776, 1118, 400, 514, 708,
1050, 674, 1244
],
[
274, 502, 502, 650, 536, 0, 228, 308, 194, 240, 582, 776, 662, 628,
514, 1050, 708
],
[
502, 730, 274, 878, 764, 228, 0, 536, 194, 468, 354, 1004, 890, 856,
514, 1278, 480
],
[
194, 354, 810, 502, 388, 308, 536, 0, 342, 388, 730, 468, 354, 320,
662, 742, 856
],
[
308, 696, 468, 844, 730, 194, 194, 342, 0, 274, 388, 810, 696, 662,
320, 1084, 514
],
[
194, 742, 742, 890, 776, 240, 468, 388, 274, 0, 342, 536, 422, 388,
274, 810, 468
],
[
536, 1084, 400, 1232, 1118, 582, 354, 730, 388, 342, 0, 878, 764,
730, 388, 1152, 354
],
[
502, 594, 1278, 514, 400, 776, 1004, 468, 810, 536, 878, 0, 114,
308, 650, 274, 844
],
[
388, 480, 1164, 628, 514, 662, 890, 354, 696, 422, 764, 114, 0, 194,
536, 388, 730
],
[
354, 674, 1130, 822, 708, 628, 856, 320, 662, 388, 730, 308, 194, 0,
342, 422, 536
],
[
468, 1016, 788, 1164, 1050, 514, 514, 662, 320, 274, 388, 650, 536,
342, 0, 764, 194
],
[
776, 868, 1552, 560, 674, 1050, 1278, 742, 1084, 810, 1152, 274,
388, 422, 764, 0, 798
],
[
662, 1210, 754, 1358, 1244, 708, 480, 856, 514, 468, 354, 844, 730,
536, 194, 798, 0
],
]
distance_matrix = {}
for i in locations:
for j in locations:
distance_matrix[(i, j)] = dist_matrix[i][j]
num_vehicles = 6
demands = [0, 1, 1, 2, 4, 2, 4, 8, 8, 1, 2, 1, 2, 4, 4, 8, 8]
vehicle_capacity = 30
depot = 0
model = cp.CpModel()
# Build decision variables :
dv = {}
for i in locations:
for j in locations:
for k in range(num_vehicles):
dv[(i, j, k)] = model.NewBoolVar("from_%i_to_%i_vehc_%i" % (i, j, k))
# Objective function
total_distance_route = sum(dv[i, j, k] * distance_matrix[i, j] for i in locations for j in locations for k in range(num_vehicles))
# Capacity constraint
for k in range(num_vehicles):
model.Add(sum(dv[i, j, k] * demands[j] for i in locations for j in locations[1:]) <= vehicle_capacity)
# # no travel from a node to itself
# for k in range(num_vehicles):
# model.Add(dv[0, 0, k] == 0)
# which nodes were covered in a tour by a vehicle
dv_vehc_node_bin = {}
for k in range(num_vehicles):
lst = {}
for i in locations[1:]:
lst[i] = model.NewBoolVar("vehc_%i_node_%i" % (k, i))
dv_vehc_node_bin[k] = lst
# each node (apart from depot) is in exactly one tour
for v in locations[1:]:
model.Add(sum(dv_vehc_node_bin[i][v] for i in range(num_vehicles)) == 1)
# sub-tour elimination
arcs_dict = {}
for k in range(num_vehicles):
arcs_list = []
for i in locations:
for j in locations:
if (i == j) & (i == 0):
arcs_list.append([0, 0, dv[0, 0, k]])
elif (i == j) & (i > 0):
arcs_list.append([i, i, ~dv_vehc_node_bin[k][i]])
else:
arcs_list.append([i, j, dv[(i, j, k)]])
arcs_dict[k] = arcs_list
for a in arcs_dict:
model.AddCircuit(arcs_dict[a])
# variable to identify which vehicle was used
dv_is_vehicle_used = {i : model.NewBoolVar("") for i in range(num_vehicles)}
# symmetry breaking constraint
for i in range(num_vehicles-1):
model.Add( dv_is_vehicle_used[i] >= dv_is_vehicle_used[i+1])
for i in dv_is_vehicle_used:
model.AddMaxEquality(dv_is_vehicle_used[i], [dv_vehc_node_bin[i][n] for n in locations[1:]])
# objective function
model.Minimize(sum(dv_is_vehicle_used[i] for i in range(num_vehicles)))
# model.Minimize(total_distance_route)
# solve the model
solver = cp.CpSolver()
solver.parameters.num_search_workers = 8
solver.parameters.log_search_progress = True
solver.log_callback = print
status = solver.Solve(model)
Optimal vehicle route looks like below
for k in range(num_vehicles):
print(k)
print([y for y, v in distance_matrix.items() if solver.Value(dv[y[0], y[1], k]) > 0])
print("\n")
0
[(0, 14), (2, 16), (4, 11), (6, 4), (11, 0), (14, 2), (15, 6), (16, 15)]
1
[(0, 7), (1, 9), (3, 1), (5, 10), (7, 8), (8, 12), (9, 5), (10, 0), (12, 13), (13, 3)]
2
[(0, 0)]
3
[(0, 0)]
4
[(0, 0)]
5
[(0, 0)]
Only 2 out of 6 vehicles were used.
total_distance_route. Would using fewer vehicles result in a better value of that function? If not then there is no incentive for the solver to use fewer vehicles. You could add a penalty for each vehicle used to your objective function. – Daniel Junglas Commented Jan 16 at 12:56