I am new to using pywrapcp library to solve vehicle routing problems, but I am well acquainted with cpsat solver. So, for a vehicle routing problem, I thought to put to use cpsat solver rather than the routing specific library because of the additional learning overhead.
I looked at the below vehicle routing problem with pickup and delivery, already solved in or tools webpage: https://developers.google.com/optimization/routing/pickup_delivery
I expressed the exact problem (same data) in cpsat solver, and got an optimal solution. Now, I wanted to check if my formulation / code is correct, so I decided to match the solution of cpsat with pywrapcp library, but found that with pywrapcp we are getting a heuristic solution (albeit a good one). With cpsat, I have got an exact solution.
So the question is how can we get an optimal solution with pywrapcp ?
pywrapcp reports the below solution:
Route for vehicle 0:
0 -> 13 -> 15 -> 11 -> 12 -> 0
Distance of the route: 1552m
Route for vehicle 1:
0 -> 5 -> 2 -> 10 -> 16 -> 14 -> 9 -> 0
Distance of the route: 2192m
Route for vehicle 2:
0 -> 4 -> 3 -> 0
Distance of the route: 1392m
Route for vehicle 3:
0 -> 7 -> 1 -> 6 -> 8 -> 0
Distance of the route: 1780m
So the maximum distance travelled by a vehicle = 2192
With cpsat formulation, I get the following:
# vehicle : [node_list, distance]
{0: [[0, 4, 3, 15, 11, 0], 2032],
1: [[0, 10, 2, 8, 7, 0], 1940],
2: [[0, 13, 16, 14, 12, 0], 2008],
3: [[0, 9, 5, 6, 1, 0], 1940]}
maximum distance travelled by a vehicle = 2032
So, cpsat reports the optimal solution which is lower than pywrapcp, but how to confirm this with pywrapcp ?
Below is my formulation with cpsat solver. The code it works fine, but any improvements are welcome.
from ortools.sat.python import cp_model as cp
locations = range(17)
data = [
# fmt: off
[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],
# fmt: on
]
distance_matrix = {}
for i in locations:
for j in locations:
distance_matrix[(i, j)] = data[i][j]
num_vehicles = 4
pickups_deliveries = [
[1, 6],
[2, 10],
[4, 3],
[5, 9],
[7, 8],
[15, 11],
[13, 12],
[16, 14],
]
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))
route_length = {}
for k in range(num_vehicles):
route_length[k] = model.NewIntVar(0, 10000, "")
for k in range(num_vehicles):
model.Add(
route_length[k]
== sum(dv[i, j, k] * distance_matrix[i, j] for i in locations for j in locations)
)
# =============================================================================
# Objective function
# =============================================================================
# minimize the length of the longest single route among all vehicles
max_route_length = model.NewIntVar(0, 10000, "max_route_length")
model.AddMaxEquality(max_route_length, [route_length[k] for k in route_length])
model.Minimize(max_route_length)
# =============================================================================
# CONSTRAINTS :
# =============================================================================
# 1. Vehicle leaves node that it enters
for j in locations:
for k in range(num_vehicles):
model.Add(
sum(dv[i, j, k] for i in locations) == sum(dv[j, i, k] for i in locations)
)
# 2. Ensure that every node is entered once
for j in locations[1:]:
model.Add(sum(dv[i, j, k] for k in range(num_vehicles) for i in locations) == 1)
# 3. Every vehicle leaves the depot
for k in range(num_vehicles):
model.Add(sum(dv[0, j, k] for j in locations[1:]) == 1)
# 4. Capacity constraint
# No capacity constraint here
# 5. no travel from a node to itself
for k in range(num_vehicles):
for i in locations:
model.Add(dv[i, i, k] == 0)
# 5. time windows constraints
# these are created to enforce precedences amongst the nodes
time_wind = {}
for i in locations:
time_wind[i] = model.NewIntVar(0, 100, "start_node_%i" % i)
model.Add(time_wind[0] == 0)
# =============================================================================
# # sub-tour elimination
# =============================================================================
# list of nodes that a vehicle visits
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
# whether a node was visited by a vehicle or not
for k in range(num_vehicles):
for i in locations[1:]:
model.AddMaxEquality(
dv_vehc_node_bin[k][i],
[dv[i, j, k] for j in locations] + [dv[j, i, k] for j in locations],
)
# if [1, 6] is one of the precedence
# then both the nodes have to be visited by the same vehicle
# so if a vehicle visits 1, it has to visit 6, if not 1 then not 6
for p in pickups_deliveries:
for k in range(num_vehicles):
model.AddImplication(dv_vehc_node_bin[k][p[0]], dv_vehc_node_bin[k][p[1]])
model.AddImplication(dv_vehc_node_bin[k][p[1]], dv_vehc_node_bin[k][p[0]])
model.AddImplication(
dv_vehc_node_bin[k][p[0]].Not(), dv_vehc_node_bin[k][p[1]].Not()
)
model.AddImplication(
dv_vehc_node_bin[k][p[1]].Not(), dv_vehc_node_bin[k][p[0]].Not()
)
# enforcing the precedences
for p in pickups_deliveries:
model.Add(time_wind[p[0]] < time_wind[p[1]])
# eliminating sub-tours using AddCircuit constraint
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):
lit = model.NewBoolVar("vehc_%i_node_%i" % (k, i))
model.AddImplication(dv_vehc_node_bin[k][i], lit.Not())
model.AddImplication(dv_vehc_node_bin[k][i].Not(), lit)
arcs_list.append([i, i, lit])
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])
solver = cp.CpSolver()
solver.parameters.num_search_workers = 8
status = solver.Solve(model)
print(solver.ResponseStats())
# output the solution
path_dict = {}
for k in range(num_vehicles):
path = [y for y, v in distance_matrix.items() if solver.Value(dv[y[0], y[1], k]) > 0]
distance = 0
for p in path:
distance += distance_matrix[p[0], p[1]]
i = 0
lst = [0]
nxt = 10000
while nxt != 0:
nxt = 0 if i == 0 else nxt
nxt = [j for i, j in path if i == nxt][0]
lst.append(nxt)
i += 1
path_dict[k] = [lst, distance]