In my previous post, I looked at how slack variables were used to
determine the maximum cost value of a node on a route, and how that
value could be assigned as a “cost” for the entire route. The
program,
vrp_node_max.py.
was developed by Mizux and is really clever. However, as I mentioned
in my previous post, the final note in the original forum thread from
Guy was a note that the approach did not seem to work in practice.
With this post, I want to explore what that means and how to fix it.
A failing case
The first step with every debugging and development exercise is to try to reproduce the problem. To do that, I modified the original program to increase the numbers of high-cost nodes, and to reduce the numbers of vehicles, until I could see similar output to what Guy described. He said:
The structure @MIZUX suggested runs and gives reasonable results. (a demo case, small scale > 16 Nodes)
However, very often the result could be easily improved by a human, for example:
With 4 Nodes with low-cost to visit plus 4 Nodes with high-cost to visit, the solver mixes Nodes in each vehicle and the result is not optimal:
- V0 > high, high, high, low
- V1> low, low, low, high
The clear assumption here is that the best solution would be to keep all the high value nodes together, and keep all the low value nodes together. In the above solution, both vehicles pay the high cost. In the preferred solution, only one vehicle would pay the high cost, while the other one could pay the low cost as it has only low cost nodes to visit.
So my goal is to reproduce that by assigning approximately equal numbers of high cost and low cost nodes, and see if the solver mixes them. To do this, I set up a little loop to assign the state, low or high, based on the node’s distance from the depot. In this case, low=8 and high=42.
def initialize_trip_values(data):
# rules:
# if distance from depot is less than 500, value of node is 8;
# if greater, value is 42
for i, d in enumerate(data['distance_matrix'][0]):
if i == 0:
data['value'].append(0)
elif d < 500:
data['value'].append(8)
else:
data['value'].append(42)
Most of the other settings are the same as before. However, in order to make sure I was testing the right things, I commented out the API call that tries to balance travel between vehicles:
# distance_dimension.SetGlobalSpanCostCoefficient(10)
I also initially set the weight of the soft lower bound to 1, so as to proceed methodically.
# Should force all others dim_two slack var to zero...
for v in range(manager.GetNumberOfVehicles()):
end = routing.End(v)
dim_two.SetCumulVarSoftUpperBound(end, 0, 1)
Running with 3 vehicles gives the following solution.
Objective: 4788
Route for vehicle 0:
N:0 value:0, one:(0,0), two:(0,0)
-> N:0 one:0 two:0
Distance of the route: 0m
Route for vehicle 1:
N:0 value:0, one:(0,0), two:(0,0)
-> N:7 value:8, one:(0,8), two:(0,0)
-> N:1 value:42, one:(8,42), two:(0,0)
-> N:4 value:42, one:(50,42), two:(0,0)
-> N:3 value:42, one:(92,42), two:(0,0)
-> N:15 value:42, one:(134,42), two:(0,0)
-> N:11 value:42, one:(176,42), two:(0,0)
-> N:12 value:8, one:(218,42), two:(0,0)
-> N:13 value:8, one:(260,42), two:(0,42)
-> N:0 one:302 two:42
Distance of the route: 2352m
Route for vehicle 2:
N:0 value:0, one:(0,0), two:(0,0)
-> N:5 value:8, one:(0,8), two:(0,0)
-> N:8 value:8, one:(8,8), two:(0,0)
-> N:6 value:42, one:(16,42), two:(0,0)
-> N:2 value:42, one:(58,42), two:(0,0)
-> N:10 value:42, one:(100,42), two:(0,0)
-> N:16 value:42, one:(142,42), two:(0,0)
-> N:14 value:8, one:(184,42), two:(0,0)
-> N:9 value:8, one:(226,42), two:(0,42)
-> N:0 one:268 two:42
Distance of the route: 2352m
Maximum of the route distances: 2352m
So here the high and low values are mixed in the two routes, as was
reported by Guy. One vehicle is unused, and both route distances are
well under the allowed maximum distance of 3500. The solver thinks
this is optimal. The
objective function of 4788 is indeed the sum of the route distances
and the final dim_two values (2352+2352+48+48). It is hard to
tell by inspection alone, but perhaps combining all of the value 8
nodes into one route will exceed the maximum allowed distance. To be
safe, I boosted the maximum distance to 10000. That produced the
following result:
Objective: 4426
Route for vehicle 0:
N:0 value:0, one:(0,0), two:(0,0)
-> N:0 one:0 two:0
Distance of the route: 0m
Route for vehicle 1:
N:0 value:0, one:(0,0), two:(0,0)
-> N:0 one:0 two:0
Distance of the route: 0m
Route for vehicle 2:
N:0 value:0, one:(0,0), two:(0,0)
-> N:9 value:8, one:(0,8), two:(0,0)
-> N:5 value:8, one:(8,8), two:(0,0)
-> N:8 value:8, one:(16,8), two:(0,0)
-> N:6 value:42, one:(24,42), two:(0,0)
-> N:2 value:42, one:(66,42), two:(0,0)
-> N:10 value:42, one:(108,42), two:(0,0)
-> N:16 value:42, one:(150,42), two:(0,0)
-> N:14 value:8, one:(192,42), two:(0,0)
-> N:13 value:8, one:(234,42), two:(0,0)
-> N:12 value:8, one:(276,42), two:(0,0)
-> N:11 value:42, one:(318,42), two:(0,0)
-> N:15 value:42, one:(360,42), two:(0,0)
-> N:3 value:42, one:(402,42), two:(0,0)
-> N:4 value:42, one:(444,42), two:(0,0)
-> N:1 value:42, one:(486,42), two:(0,0)
-> N:7 value:8, one:(528,42), two:(0,42)
-> N:0 one:570 two:42
Distance of the route: 4384m
This solution has a better objective function as the path is shorter, but it lumps all the nodes into one route.
The next run keeps the maximum distance at 10k, and increased the weight of the soft lower bound to 10000, so that a value of 42 will produce 420,000, which will definitely dominate the distance (4384) in the objective function. That produced the same exact solution, but with a new objective value of 424,384.
In a way this makes sense. If two vehicles are used, then one of them will have at least 80,000 additional points added to the objective function, while using just one vehicle saves that cost. So back to reducing the maximum distance to force two vehicles to be used, setting the maximum to 4000. At that limit, the all-on-one solution is denied, as the total route distance of 4384m will exceed 4000. I expected that the optimal solution would be to put most of the nodes into one vehicle, and then put a few value 8 nodes onto the second one. It turns out this was not true:
Objective: 844704
Route for vehicle 0:
N:0 value:0, one:(0,0), two:(0,0)
-> N:0 one:0 two:0
Distance of the route: 0m
Route for vehicle 1:
N:0 value:0, one:(0,0), two:(0,0)
-> N:7 value:8, one:(0,8), two:(0,0)
-> N:1 value:42, one:(8,42), two:(0,0)
-> N:4 value:42, one:(50,42), two:(0,0)
-> N:3 value:42, one:(92,42), two:(0,0)
-> N:15 value:42, one:(134,42), two:(0,0)
-> N:11 value:42, one:(176,42), two:(0,0)
-> N:12 value:8, one:(218,42), two:(0,0)
-> N:13 value:8, one:(260,42), two:(0,42)
-> N:0 one:302 two:42
Distance of the route: 2352m
Route for vehicle 2:
N:0 value:0, one:(0,0), two:(0,0)
-> N:9 value:8, one:(0,8), two:(0,0)
-> N:14 value:8, one:(8,8), two:(0,0)
-> N:16 value:42, one:(16,42), two:(0,0)
-> N:10 value:42, one:(58,42), two:(0,0)
-> N:2 value:42, one:(100,42), two:(0,0)
-> N:6 value:42, one:(142,42), two:(0,0)
-> N:8 value:8, one:(184,42), two:(0,0)
-> N:5 value:8, one:(226,42), two:(0,42)
-> N:0 one:268 two:42
Distance of the route: 2352m
This is the same as the very first solution! A better objective could definitely be obtained by using the all-on-one solution, and just removing one value 8 node to the second vehicle. In that case, the second vehicle would pay just 80,000 for the value 8 node, giving an objective on the order of 500,000, rather than 840,000.
I think I have reproduced the issue that Guy observed, as something is
very wrong here. I suspect that the issue is the second order effects
of accumulating the maximum value node via the slack values is not
working well with the search heuristics used by the solver. My guess
is that the solver expects that each node is more or less independent
of the others, in that picking up a node and moving it to a different
route or to a different place in a route will not affect the stored
values of any other nodes. However, in this case, it does affect
other nodes. Further, the effect is not consistent and deterministic.
Moving a single value 42 node from one route to another can have no
effect if there are other value 42 nodes remaining in the first route.
On the other hand, moving the very last value 42 node has a drastic
effect, in that all of the nodes remaining on the route will drop to
dim_two cost of 8. The solver has no easy way to figure this out,
so I suspect that even though using the slack to carry along the
maximum value is extremely clever, in the end it is useless.
A different approach
Since writing the previous post on this topic, I’ve been thinking that the best way forward would be to switch to using the CP-SAT solver, and treating this more like a knapsack problem of some sort. But a few days ago, yet another question floated by on the mailing list or in the GitHub discussions asking how to achieve time-dependent travel times. My thought has always been that the best way to do that is to use the usual, tried and true, ever-since the 50’s approach of using time of day periods. It is well known that traffic varies by time of day, but until very recently (thanks to GPS on cellphones), it has been next to impossible to estimate instantaneous travel times on arbitrary links by time of day. Instead, traffic engineers used time periods and weights. The time periods are tuned to the area, but typically are over night (low volume, free flow speeds on all links), morning peak (short sharp congestion effects), midday (everything is busy), and afternoon peak (flatter than the morning peak, and longer-tailed). This same approach can be used with OR-Tools, by duplicating every node into the 4 different time of day periods, and then limiting access based on the value of the dimension. You can only arrive at a morning peak node if the CumulVar of time is inside the morning peak, and links to and from the morning peak nodes use the slower morning peak travel times. And so on.
I didn’t reply to the question because I already wrote that up at least once on the mailing list and so the original poster either didn’t search the list first, or else did not like that solution. But thinking about it, I realized the same approach can be used for this problem. The basic idea is that if a low value node is included in a route with at least one high value node, then the low value node becomes a high value node. So if all low value nodes are duplicated into the high value set with a suitable disjunction so that only one of the duplicated nodes is visited, and if things are set up so that only high value nodes can reach other, and only low value nodes can reach each other, then the solver will not have to have magical powers to solve the problem. Instead we can adhere to the usual convention that a node’s incremental addition to the totals of the dimensions on the route is independent from all other nodes.
Grouping and duplicating nodes by value
I hacked out a quick implementation that is bloated and messy, but works okay. My ulterior motive for writing this post is really to go through my own code, clean it up, and explain the moving parts in case I need to use it for a real problem some day. With that in mind, the first place to start is in the dimension callbacks, which will dictate what needs to happen when I duplicate nodes.
Again the idea behind duplicating nodes is that only one of the many duplicates should be visited. However, the solver thinks each node is unique, so the callback needs to understand that the nodes passed to it are dummy nodes, and how to translate those dummies back to the original node value.
The distance callback function
For the distance callback function, handling the duplicate nodes looks like this:
def distance_callback(from_index, to_index):
"""Returns the distance between the two nodes."""
# Convert from routing variable Index to distance matrix NodeIndex.
from_dummy_node = manager.IndexToNode(from_index)
to_dummy_node = manager.IndexToNode(to_index)
# convert the dummy nodes to the nodes
from_node = data['node_from_dummynode_lookup'][from_dummy_node]
to_node = data['node_from_dummynode_lookup'][to_dummy_node]
# look up the distance between nodes
distance = data['distance_matrix'][from_node][to_node]
# if the from_dummy and to_dummy are not in the same value
# group, return infinity
# the depot can get to any node
if from_node == 0 or to_node == 0:
return distance
from_value_group = data['value_group_from_dummy_node'][from_dummy_node]
to_value_group = data['value_group_from_dummy_node'][to_dummy_node]
if from_value_group == to_value_group:
return distance
return distance * 10000
All fairly straightforward. The solver passes indices (why I’ll never know), which are then converted immediately into the correct dummy nodes (why doesn’t the solver just pass nodes?). Then those dummy nodes need to be converted back to the original node (so that means I need to keep track of the conversion using a lookup dictionary). The real nodes are used to compute the distance. In order to allow movements only between dummy nodes with the same value, I have to look up the corresponding value of each dummy node (which is yet another lookup dictionary that needs to be generated).
The value dimension callback function
Rather than the prior solutions combination of two dimension, here I just use one dimension to track the value of the nodes on the route. The callback function returns the value if both the from node and to node are in the same value group, or else it returns infinity. Trips from the depot have a cost of 0, and trips returning to the depot have the value of the from node.
def value_callback(from_index, to_index):
"""Returns the value of the node group."""
# Convert from routing variable Index to distance matrix NodeIndex.
from_dummy_node = manager.IndexToNode(from_index)
to_dummy_node = manager.IndexToNode(to_index)
# if coming from the depot, return 0
if from_dummy_node == 0:
return 0
# if returning to the depot, return the from node value
if to_dummy_node == 0:
return data['value_group_from_dummy_node'][from_dummy_node]
# compare the from and to values
from_value_group = data['value_group_from_dummy_node'][from_dummy_node]
to_value_group = data['value_group_from_dummy_node'][to_dummy_node]
# if they are the same, return that value
if from_value_group == to_value_group:
return from_value_group
# if they are not in the same group, return infinity
return data['maximum_value'] * 1_000_000
This function is similar to the distance function. Again, a lookup dictionary is needed to find the value of the dummy node. No other lookup dictionary seems to be needed.
Prevent double visits using disjunctions
In order to prevent visiting the same real node twice, the disjunction API call is used with a list of nodes. The call will by default visit at most one of the nodes in the list. If none of the nodes are visited, then the penalty value is added to the objective function. The disjunction call has to iterate over all of the real nodes, to make sure that each node is only visited once, as follows:
# disjunctions for each value grouped node
for node in range(1, len(data['value'])):
duplicate_nodes = [manager.NodeToIndex(n) for n in data['dummy_nodes_lookup'][node]]
routing.AddDisjunction(
duplicate_nodes, 10000000000
)
The depot node (0 in this case) is not optional, so it is skipped in
the loop. For each real node, the list of dummy nodes is loaded, and
then passed to the AddDisjunction API call, with an arbitrary, large
penalty to ensure that the highest priority is given to serving all
nodes given the constraints. By default, the maximum cardinality of
the disjunction is 1, so the third parameter in the API call is
skipped. This exposes a need for a final lookup dictionary, one that
obtains all duplicate nodes given the real node.
Generate the duplicated nodes
The duplicate nodes are generated using the following idea. First, cheap nodes (value 8, etc), can be “promoted” to be expensive nodes, but expensive nodes cannot be demoted to become cheap nodes. It follows that the cheapest nodes will be duplicated the most times, and the most expensive nodes will only be “duplicated” once. In the test case with just two values possible (8 and 42), the value 8 nodes will be duplicated two times. My function is as follows:
def create_dummy_nodes_for_value_groups(data):
# create the dictionaries we need in callback fns, disjunction call
data['node_from_dummynode_lookup'] = {}
data['value_group_from_dummy_node'] = {}
data['dummy_nodes_lookup'] = {}
data['value_groups'] = {}
maximum_value = 0
# pick off unique values, create groups
for i, v in enumerate(data['value']):
# create a list to hold dummy nodes for this node i
data['dummy_nodes_lookup'][i] = []
if i == 0:
# skip the depot node
continue
if v in data['value_groups']:
# append this node to the list
data['value_groups'][v].append(i)
else:
# create the list for this value
data['value_groups'][v] = [i]
if v > maximum_value:
# update the running maximum
maximum_value = v
data['maximum_value'] = maximum_value
# duplicate nodes as needed.
#
# Iterate over each node, value pair, and generate a duplicate
# node as needed for each of the relevant value groups...those
# that are equal to or greater than the value
dummy_node = -1 # initialize such that first node is node 0
for i, v in enumerate(data['value']):
if v == maximum_value or i == 0:
dummy_node += 1
# cannot boost max cost nodes (or depot) to a higher cost group
data['dummy_nodes_lookup'][i] = [dummy_node]
data['node_from_dummynode_lookup'][dummy_node] = i
data['value_group_from_dummy_node'][dummy_node] = v
data['dummy_nodes'].append(dummy_node)
# skip to the next node
continue
# iterate over all of the value groups
for boost_value in data['value_groups'].keys():
if boost_value < v:
# cannot demote from value v to a lesser "boost" value
# so skip to the next value group
continue
# Still here? Can reach this node, either by boosting or by regular
dummy_node += 1
data['node_from_dummynode_lookup'][dummy_node] = i
data['dummy_nodes_lookup'][i].append(dummy_node)
data['value_group_from_dummy_node'][dummy_node] = boost_value
data['dummy_nodes'].append(dummy_node)
This function can probably be broken up into two. The first half has a loop that finds each of the potential values. In the demo case there are just two values, but I want to test with more so I wrote this more general approach.
The second part of the function loops over each of the real nodes and generates all of the dummy nodes that are required. The first case is if the node is either the depot node or has the maximum value. In the depot case, depots are never evaluated on the value dimension, but I do need to keep track of the correspondence between real node and dummy node. Similarly, the nodes that are already in the highest value group only get one “duplicate” node, but I still need to keep track of this value.
In the more general case, the value of the node is not already at the
maximum. In that case, the code iterates over all of the discovered
value groups. If the value of the group is less than the node value,
then the entry is skipped. If the value of the group (boost_value)
is equal to or greater than the value of the node, then a new dummy
node is created, and all of the relevant dictionaries and lists are
updated as needed.
Creating the dimensions
Putting it all together, the distance dimension and the new value dimension need to be created. The distance dimension is mostly the same as the prior version. The value dimension is defined as follows:
routing.AddDimension(
value_callback_index,
0, # no slack
data['maximum_value'] * num_nodes, # essentially infinity
True, # start cumul to zero
'Value')
value_dimension = routing.GetDimensionOrDie('Value')
Here I do not want the capacity of the dimension to be a limiting factor for normal nodes, so I set it to be equal to the maximum value times the total number of nodes. Recall in the value callback function I return “infinity” as the maximum value times one million. As long as we are not processing on the order of a million nodes, this is okay.
The final step is to tell the solver to minimize the value dimension
at the end of each route. This is more or less the same as what the
original code did for dim_two.
# Should force minimizing the cost of each route
for v in range(manager.GetNumberOfVehicles()):
end = routing.End(v)
value_dimension.SetCumulVarSoftUpperBound(end, 0, 1)
Running the revised code
The first run uses a weight of 1 in the soft upper bound for the value dimension. The resulting routes are almost identical to the previous solution.
Objective: 5376
Route for vehicle 0:
N:0 base value:0, real value:0, one:(0)
-> N:13 base value:8, real value:42, one:(0)
-> N:12 base value:8, real value:42, one:(42)
-> N:11 base value:42, real value:42, one:(84)
-> N:15 base value:42, real value:42, one:(126)
-> N:3 base value:42, real value:42, one:(168)
-> N:4 base value:42, real value:42, one:(210)
-> N:1 base value:42, real value:42, one:(252)
-> N:7 base value:8, real value:42, one:(294)
-> N:0 one:336
Distance of the route: 2352m
Route for vehicle 1:
N:0 base value:0, real value:0, one:(0)
-> N:9 base value:8, real value:42, one:(0)
-> N:14 base value:8, real value:42, one:(42)
-> N:16 base value:42, real value:42, one:(84)
-> N:10 base value:42, real value:42, one:(126)
-> N:2 base value:42, real value:42, one:(168)
-> N:6 base value:42, real value:42, one:(210)
-> N:8 base value:8, real value:42, one:(252)
-> N:5 base value:8, real value:42, one:(294)
-> N:0 one:336
Distance of the route: 2352m
Route for vehicle 2:
N:0 base value:0, real value:0, one:(0)
-> N:0 one:0
Distance of the route: 0m
Maximum of the route distances: 2352m
Each of the two assigned vehicles has a route distance of 2352m, and each route has a mix of low value and high value nodes. However, if I boost the weight assigned to the soft lower bound to be 10,000 (so that it dominates the trip distance in the objective function) as follows:
for v in range(manager.GetNumberOfVehicles()):
end = routing.End(v)
value_dimension.SetCumulVarSoftUpperBound(end, 0, 10_000)
then I see the following result:
Objective: 4346416
Route for vehicle 0:
N:0 base value:0, real value:0, one:(0)
-> N:11 base value:42, real value:42, one:(0)
-> N:15 base value:42, real value:42, one:(42)
-> N:3 base value:42, real value:42, one:(84)
-> N:4 base value:42, real value:42, one:(126)
-> N:1 base value:42, real value:42, one:(168)
-> N:0 one:210
Distance of the route: 2192m
Route for vehicle 1:
N:0 base value:0, real value:0, one:(0)
-> N:16 base value:42, real value:42, one:(0)
-> N:10 base value:42, real value:42, one:(42)
-> N:2 base value:42, real value:42, one:(84)
-> N:6 base value:42, real value:42, one:(126)
-> N:0 one:168
Distance of the route: 2192m
Route for vehicle 2:
N:0 base value:0, real value:0, one:(0)
-> N:9 base value:8, real value:8, one:(0)
-> N:5 base value:8, real value:8, one:(8)
-> N:8 base value:8, real value:8, one:(16)
-> N:14 base value:8, real value:8, one:(24)
-> N:13 base value:8, real value:8, one:(32)
-> N:12 base value:8, real value:8, one:(40)
-> N:7 base value:8, real value:8, one:(48)
-> N:0 one:56
Distance of the route: 2032m
Maximum of the route distances: 2192m
At last, this is the desired output. The first two vehicles share the value 42 nodes, and the last vehicle serves all of the value 8 nodes. The solver can easily recognize that serving the higher value versions of the value 8 nodes is more expensive to the objective function than directly serving the lower value nodes in a low-value-specific route. The total trip distance is of course greater than the previous mixed solution, but that is okay because we really prefer to minimize the value of the nodes served.
As a final test, I changed the value assignment step to divide the distance by 100, as follows:
for i, d in enumerate(data['distance_matrix'][0]):
if i == 0:
data['value'].append(0)
elif d < 500:
data['value'].append(int(d/100)+1) # (8)
else:
data['value'].append(int(d/100)+1) # (42)
The distances from the depot span values from 194 though 776, so the “values” should range from 2 through 8. The solver should try to minimize the value assigned to each route, which it seems to do successfully:
Objective: 968400
Route for vehicle 0:
N:0 base value:0, real value:0, one:(0)
-> N:15 base value:8, real value:8, one:(0)
-> N:2 base value:8, real value:8, one:(8)
-> N:0 one:16
Distance of the route: 3104m
Route for vehicle 1:
N:0 base value:0, real value:0, one:(0)
-> N:11 base value:6, real value:7, one:(0)
-> N:4 base value:6, real value:7, one:(7)
-> N:3 base value:7, real value:7, one:(14)
-> N:1 base value:6, real value:7, one:(21)
-> N:6 base value:6, real value:7, one:(28)
-> N:10 base value:6, real value:7, one:(35)
-> N:16 base value:7, real value:7, one:(42)
-> N:14 base value:5, real value:7, one:(49)
-> N:0 one:56
Distance of the route: 3424m
Route for vehicle 2:
N:0 base value:0, real value:0, one:(0)
-> N:7 base value:2, real value:4, one:(0)
-> N:12 base value:4, real value:4, one:(4)
-> N:13 base value:4, real value:4, one:(8)
-> N:9 base value:2, real value:4, one:(12)
-> N:8 base value:4, real value:4, one:(16)
-> N:5 base value:3, real value:4, one:(20)
-> N:0 one:24
Distance of the route: 1872m
Maximum of the route distances: 3424m
The solver splits off the two most expensive nodes with a value of 8 into their own route, assigns the most nodes to a route with a maximum value of 7, and then uses the third vehicle to serve nodes with values of 2, 3, and 4.
There are only a two questions left in my mind. First is how scalable this solution is. It seems to run and converge right away, but as everybody knows who has used OR-Tools on real data, large data sets can always pose problems. Second, I wonder about how this solver formulation will perform if nodes must be dropped. I suspect that the current disjunction value is too high, and that the solver will end up generating routes that are too expensive in order to serve more nodes. If the disjunction value is dropped below the soft bound weight, then the solver might end up dropping too many nodes. I think the best approach will be to carefully assess the value of serving the node versus the cost of keeping it in, but that would take real data to do properly.
The full program
The final program is pasted below.
#!/usr/bin/env python3
# Copyright 2010-2021 Google LLC
# modifications related to the Value dimension Copyright 2022 Activimetrics LLC
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
#
#
# [START program]
"""Vehicles Routing Problem (VRP).
Each route as an associated objective cost equal to the max node value along the
road multiply by a constant factor (4200)
"""
# [START import]
from ortools.constraint_solver import routing_enums_pb2
from ortools.constraint_solver import pywrapcp
# [END import]
# [START data_model]
def create_data_model():
"""Stores the data for the problem."""
data = {}
data['distance_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
],
]
data['value'] = []
data['num_vehicles'] = 3
data['depot'] = 0
initialize_trip_values(data)
create_dummy_nodes_for_value_groups(data)
return data
# [END data_model]
def initialize_trip_values(data):
# rule: if distance from depot is less than 500, value of node is 8; if greater, value is 42
for i, d in enumerate(data['distance_matrix'][0]):
if i == 0:
data['value'].append(0)
elif d < 500:
data['value'].append(8) #int(d/100)+1)
else:
data['value'].append(42) #int(d/100)+1)
assert len(data['distance_matrix']) == len(data['value'])
assert len(data['distance_matrix'][0]) == len(data['value'])
def create_dummy_nodes_for_value_groups(data):
# create the dictionaries we need in callback fns, disjunction call
data['node_from_dummynode_lookup'] = {}
data['value_group_from_dummy_node'] = {}
data['dummy_nodes_lookup'] = {}
data['value_groups'] = {}
maximum_value = 0
# pick off unique values, create groups
for i, v in enumerate(data['value']):
# create a list to hold dummy nodes for this node i
data['dummy_nodes_lookup'][i] = []
if i == 0:
# skip the depot node
continue
if v in data['value_groups']:
# append this node to the list
data['value_groups'][v].append(i)
else:
# create the list for this value
data['value_groups'][v] = [i]
if v > maximum_value:
# update the running maximum
maximum_value = v
data['maximum_value'] = maximum_value
# duplicate nodes as needed.
#
# Iterate over each node, value pair, and generate a duplicate
# node as needed for each of the relevant value groups...those
# that are equal to or greater than the value
data['dummy_nodes'] = []
dummy_node = -1 # initialize such that first node is node 0
for i, v in enumerate(data['value']):
if v == maximum_value or i == 0:
dummy_node += 1
# cannot boost max cost nodes (or depot) to a higher cost group
data['dummy_nodes_lookup'][i] = [dummy_node]
data['node_from_dummynode_lookup'][dummy_node] = i
data['value_group_from_dummy_node'][dummy_node] = v
data['dummy_nodes'].append(dummy_node)
# skip to the next node
continue
# iterate over all of the value groups
for boost_value in data['value_groups'].keys():
if boost_value < v:
# cannot demote from value v to a lesser "boost" value
# so skip to the next value group
continue
# Still here? Can reach this node, either by boosting or by regular
dummy_node += 1
data['node_from_dummynode_lookup'][dummy_node] = i
data['dummy_nodes_lookup'][i].append(dummy_node)
data['value_group_from_dummy_node'][dummy_node] = boost_value
data['dummy_nodes'].append(dummy_node)
print(f"dummy_node is {dummy_node}, len dummy nodes is {len(data['dummy_nodes'])}, and max value of dummy nodes is {max(data['dummy_nodes'])}")
# [START solution_printer]
def print_solution(data, manager, routing, solution):
"""Prints solution on console."""
print(f'Objective: {solution.ObjectiveValue()}')
max_route_distance = 0
dim_one = routing.GetDimensionOrDie('Value')
for vehicle_id in range(data['num_vehicles']):
index = routing.Start(vehicle_id)
plan_output = 'Route for vehicle {}:\n '.format(vehicle_id)
route_distance = 0
while not routing.IsEnd(index):
one_var = dim_one.CumulVar(index)
# one_slack_var = dim_one.SlackVar(index)
dummy_node = manager.IndexToNode(index)
node = data['node_from_dummynode_lookup'][dummy_node]
plan_output += f" N:{node} base value:{data['value'][node]}, real value:{data['value_group_from_dummy_node'][dummy_node]}, one:({solution.Value(one_var)}) \n -> "
index = solution.Value(routing.NextVar(index))
next_dummy_node = manager.IndexToNode(index)
next_node = data['node_from_dummynode_lookup'][next_dummy_node]
route_distance += data['distance_matrix'][node][next_node]
one_var = dim_one.CumulVar(index)
plan_output += f"N:{manager.IndexToNode(index)} one:{solution.Value(one_var)}\n"
plan_output += 'Distance of the route: {}m\n'.format(route_distance)
print(plan_output)
max_route_distance = max(route_distance, max_route_distance)
print('Maximum of the route distances: {}m'.format(max_route_distance))
# [END solution_printer]
def main():
"""Solve the CVRP problem."""
# Instantiate the data problem.
# [START data]
data = create_data_model()
# [END data]
# Create the routing index manager.
num_nodes = len(data['dummy_nodes'])
# [START index_manager]
manager = pywrapcp.RoutingIndexManager(
num_nodes,
data['num_vehicles'],
data['depot'])
# [END index_manager]
# Create Routing Model.
# [START routing_model]
model_parameters = pywrapcp.DefaultRoutingModelParameters()
model_parameters.max_callback_cache_size = 2 * num_nodes * num_nodes
# Create Routing Model.
routing = pywrapcp.RoutingModel(manager, model_parameters)
# [END routing_model]
# Create and register a transit callback.
# [START transit_callback]
def distance_callback(from_index, to_index):
"""Returns the distance between the two nodes."""
# Convert from routing variable Index to distance matrix NodeIndex.
from_dummy_node = manager.IndexToNode(from_index)
to_dummy_node = manager.IndexToNode(to_index)
# convert the dummy nodes to the nodes
from_node = data['node_from_dummynode_lookup'][from_dummy_node]
to_node = data['node_from_dummynode_lookup'][to_dummy_node]
# look up the distance between nodes
distance = data['distance_matrix'][from_node][to_node]
# if the from_dummy and to_dummy are not in the same value
# group, return infinity
# the depot can get to any node
if from_node == 0 or to_node == 0:
return distance
from_value_group = data['value_group_from_dummy_node'][from_dummy_node]
to_value_group = data['value_group_from_dummy_node'][to_dummy_node]
if from_value_group == to_value_group:
return distance
return distance * 10000
transit_callback_index = routing.RegisterTransitCallback(distance_callback)
# [END transit_callback]
# Define cost of each arc.
# [START arc_cost]
routing.SetArcCostEvaluatorOfAllVehicles(transit_callback_index)
# [END arc_cost]
# Add Distance constraint.
# [START distance_constraint]
dimension_name = 'Distance'
routing.AddDimension(
transit_callback_index,
0, # no slack
3_500, # vehicle maximum travel distance
True, # start cumul to zero
dimension_name)
distance_dimension = routing.GetDimensionOrDie(dimension_name)
# distance_dimension.SetGlobalSpanCostCoefficient(10)
# [END distance_constraint]
# cost of a link is based on the highest cost zone reached by trip.
# disjunctions for each value grouped node
for node in range(1, len(data['value'])):
duplicate_nodes = [manager.NodeToIndex(n) for n in data['dummy_nodes_lookup'][node]]
routing.AddDisjunction(
duplicate_nodes, 10000000000
)
# Max Node value Constraint.
# Create and register a transit callback.
# [START transit_callback]
def value_callback(from_index, to_index):
"""Returns the value of the node group."""
# Convert from routing variable Index to distance matrix NodeIndex.
from_dummy_node = manager.IndexToNode(from_index)
to_dummy_node = manager.IndexToNode(to_index)
# if coming from the depot, return 0
if from_dummy_node == 0:
return 0
# if returning to the depot, return the from node value
if to_dummy_node == 0:
return data['value_group_from_dummy_node'][from_dummy_node]
# compare the from and to values
from_value_group = data['value_group_from_dummy_node'][from_dummy_node]
to_value_group = data['value_group_from_dummy_node'][to_dummy_node]
# if they are the same, return that value
if from_value_group == to_value_group:
return from_value_group
# if they are not in the same group, return infinity
return data['maximum_value'] * 1_000_000
value_callback_index = routing.RegisterTransitCallback(value_callback)
# [END transit_callback]
# assert 0
# Dimension One will be used to compute the max node value up to the node in
# the route and store the result in the SlackVar of the node.
routing.AddDimension(
value_callback_index,
0, # no slack
data['maximum_value'] * num_nodes, # essentially infinity
True, # start cumul to zero
'Value')
value_dimension = routing.GetDimensionOrDie('Value')
# Should force minimizing the cost of each route
for v in range(manager.GetNumberOfVehicles()):
end = routing.End(v)
value_dimension.SetCumulVarSoftUpperBound(end, 0, 1)
# try 10_000 to prefer lower value routes over shorter routes
# value_dimension.SetCumulVarSoftUpperBound(end, 0, 10_000)
# Setting first solution heuristic.
# [START parameters]
search_parameters = pywrapcp.DefaultRoutingSearchParameters()
search_parameters.first_solution_strategy = (
routing_enums_pb2.FirstSolutionStrategy.PATH_CHEAPEST_ARC)
search_parameters.local_search_metaheuristic = (
routing_enums_pb2.LocalSearchMetaheuristic.GUIDED_LOCAL_SEARCH)
search_parameters.log_search = False
search_parameters.time_limit.FromSeconds(2)
# [END parameters]
# Solve the problem.
# [START solve]
solution = routing.SolveWithParameters(search_parameters)
# [END solve]
# Print solution on console.
# [START print_solution]
if solution:
print_solution(data, manager, routing, solution)
else:
print('No solution found !')
# [END print_solution]
if __name__ == '__main__':
main()
# [END program]