More fun with Slack: Understanding vrp_node_max.py 2022-03-22

A few weeks ago while working on a project I came across a very interesting program in the OR-Tools code base. The origin of the code started from a question from Guy on the OR-Tools forum. When the question was originally posted last summer, I didn’t have the opportunity to dig into the proposed result, but now I do, and I discovered it is a very interesting (ab)use of the slack dimension. Since I haven’t written anything here in ages, it seemed like a good idea to write up a blog post on the code and my understanding of how it works.

The original question

In the original thread on the forum, Guy asked about how to use the vehicle routing solver for a problem in which the cost of each route was based on the most expensive node reached. In his words, he was “looking for a ‘peak-detector’, that reflects only the cost of the highest cost Node.” The difficulty is that in a routing problem the usual approach is to create dimensions that collect something at each node visited. For example, “time” at a location can be represented by summing the service time at the previous node, plus the time required to get from the previous node to the current node. The main constraint is that the value of these dimension can only be dependent on the current node and the node immediately prior. The value of the dimension at the node should be independent of the value at all other locations. This is what allows the solver to add and drop nodes in a path without having to worry about the underlying objective function changing in unexpected ways.

The clever hack that figures out the maximum cost node

The solution suggested by Mizux from the OR-Tools team was to use the slack value at a node to store the value of interest. This is not done in the dimension callback, but rather by setting up constraints prior to running the solver. This works because the value at a node is determined by the following relationship (one that I’ve talked about previously):

Given:
 i = current node
 j = next node
 cumul(x) = quantity at node x
Then:
 cumul(j) = cumul(i) + transit(i, j) + slack(i)

In Guy’s problem, the cost of reaching some node i is fixed, based on its original distance from the depot. Therefore it is not necessary to know anything about the other nodes on the path to know the cost of reaching that node.

The creative solution hit upon by Mizux is to set up a recursive relationship between every pair of nodes (i, j), and to use the slack variable at i to store the maximum cost of either i or j. Recursive algorithms are cool, but they used to give me a headache when I first studied them in school. Conceptually, if you start from the vehicle start, then the first location is the most expensive node reached so far. The trouble is that each node can potentially be served by any of the vehicles’ start nodes as a potential first location. The trick that Mizux’s code uses is to multiply each value by a boolean variable that indicates if the vehicle start is actually the preceding node. Each potential value is then stored in an array. The array can include at most one non-zero value (because only one node can be the actual predecessor node) and so summing the array will always give the maximum cost of either the predecessor or the current node.

So that’s the key insight, plus a few more clever hacks that make it work. Next I’ll step through the code to explain how it all works.

Dim One is for the recursion

The code reviewed here can be found in its entirety in the GitHub repository here. The beginning follows the usual pattern of the example files—setting up the data, defining the solution printer function, etc. The interesting stuff gets going around line 214, with the creation of dim_one.

routing.AddConstantDimensionWithSlack(
    0,  # transit 0
    42 * 16,  # capacity: be able to store PEAK*ROUTE_LENGTH in worst case
    42,  # slack_max: to be able to store peak in slack
    True,  #  Fix StartCumulToZero not really matter here
    'One')
dim_one = routing.GetDimensionOrDie('One')

In a real project, I would of course change the hard coded numbers to variables, something like

routing.AddConstantDimensionWithSlack(
    0,  # transit 0
    peak_value * num_nodes,  # capacity: be able to store PEAK*ROUTE_LENGTH in worst case
    peak_value,  # slack_max: to be able to store peak in slack
    True,  #  Fix StartCumulToZero not really matter here
    'One')

While the inline comments are helpful, I want to go over some important points here. First of all, dim_one is a “constant” dimension. A constant dimension increases by a constant value each time a vehicle visits a node. For example, I often use a constant dimension with a value of 1 to count the nodes visited on a route.

In this case, the constant value added to the dimension at each node is zero. This does not mean that the dimension will sum to zero. Rather it just means that moving from node to node in and of itself does not add any value to the dimension.

Second, the capacity of the dimension is set quite high. This is to ensure that the dimension is never forced to avoid a node because its value is too high. In the worst case, the first node visited is the most expensive node, and then all of the nodes are on a single route. In that case, each slack variable will be set to the maximum cost value, and that value will be summed up for each visit to each node. So the worst case upper bound for the dimension is the maximum value of the data['valuie'] array (in this case, 42), times the total number of nodes (in this case, 16).

The recursive relationship uses dim_one, so I will focus on that for now and ignore dim_two, even though they are defined at the same time in the code.

Initialize the start nodes

Every recursive algorithm needs to first start by setting up the boundary conditions. In this case, the first location of each trip is a vehicle start node. At those nodes, the slack value is undefined if you follow the algorithm outlined above, because there is no predecessor node to use in the max()expression. Therefore the slack value is fixed to the value assigned to node zero.

# force depot Slack to be value since we don't have any predecessor...
for v in range(manager.GetNumberOfVehicles()):
    start = routing.Start(v)
    dim_one.SlackVar(start).SetValue(data['value'][0])

The data['value'] array is equal to zero at the starting location, but it could be anything.

With the initial conditions set up, the next step is to define the recursive relationship itself. This is defined for every node other than 0, because the 0 node represents the vehicle start node or depot value, which has already been taken care of.

Each node must consider all possible prior nodes. This is made a little bit more complicated because the 0 node is the depot, which is really one unique node per vehicle. So the set of possible prior nodes is defined as the combination of all of the vehicle start nodes, and all of the regular nodes from 1 to N.

For every possible prior node, the test is the same. First, a zero-one variable is created that is 1 if the possible prior is the predecessor node of the current node, and 0 if it is not. Then the potential value is computed as the maximum of that potential prior and the current node. The zero-one boolean is then multiplied by the value, to create a variable that is zero if the possible predecessor is not actually the predecessor. This code looks like this:

cond = routing.NextVar(previous_index) == index
value = solver.Max(dim_one.SlackVar(previous_index),
                   data['value'][node])
test.append((cond * value).Var())

The variable cond uses a trick that is rarely seen in the example programs, but that is quite useful. The routing solver keeps track of all possible next nodes for every node. In a generic problem, that creates an all-to-all matrix of boolean values. Each node i can only transition to one other node, and that value at each iteration of the solver is contained in routing.NextVar(i). This can be used to create a boolean expression. The Python bindings for OR-Tools are nice here in that they hide the details behind some syntactic sugar. In C++, the same call would be something like:

IntVar* cond = solver->MakeIsEqualCstVar(routing.NextVar(previous_index), index);

The last thing that is unusual compared to the other example programs is the line

test.append((cond * value).Var())

Looking at this from the inside out, the most important part is the combination of (cond * value). This multiplies the value by zero or one and makes the magic happen. But then right after then close parenthesis, the function Var() is called. Again, this is not usually done in the example programs. What is happening is hidden a bit behind the syntactic sugar of the Python API. Behind the scenes, the multiplication of two IntVar types creates an IntExpr. An IntVar is a subset of the IntExpr class. Quoting the source code:

The class IntVar is a subset of IntExpr. In addition to the IntExpr protocol, it offers persistence, removing values from the domains, and a finer model for events.

Every IntExpr can be turned into an IntVar by calling the virtual method Var(). If it is not allowed, the default return value is 0, but that is really only an issue when calling from C++.

The reason the Var() method must be called is that later the code sums up all of the expressions. The API has three methods available for summing up IntExpr values:

  // ----- Integer Expressions -----

  /// left + right.
  IntExpr* MakeSum(IntExpr* const left, IntExpr* const right);
  /// expr + value.
  IntExpr* MakeSum(IntExpr* const expr, int64_t value);
  /// sum of all vars.
  IntExpr* MakeSum(const std::vector<IntVar*>& vars);

In this case we are summing up a vector of values, and the only call available is the one that has a vector of IntVars.

  /// sum of all vars.
  IntExpr* MakeSum(const std::vector<IntVar*>& vars);

So the call to Var() is required because later we want to call Sum() over the array test.

The full loop first examines the vehicle start nodes as the predecessor node, and then all the regular nodes, as follows:

# Step by step relation
# Slack(N) = max( Slack(N-1) , value(N) )
solver = routing.solver()
for node in range(1, 17):
    index = manager.NodeToIndex(node)
    test = []
    for v in range(manager.GetNumberOfVehicles()):
        previous_index = routing.Start(v)
        cond = routing.NextVar(previous_index) == index
        value = solver.Max(dim_one.SlackVar(previous_index),
                           data['value'][node])
        test.append((cond * value).Var())
    for previous in range(1, 17):
        previous_index = manager.NodeToIndex(previous)
        cond = routing.NextVar(previous_index) == index
        value = solver.Max(dim_one.SlackVar(previous_index),
                           data['value'][node])
        test.append((cond * value).Var())
    solver.Add(solver.Sum(test) == dim_one.SlackVar(index))

The last line in the routine creates a constraint for each regular node. The constraint is that the slack value at the current node is equal to the sum of the test array. That summation sums up a bunch of zeros, and one non-zero value, as only one of the previous_index values can actually be the prior node. One could also use the Max() API call, as in:

solver.Add(solver.Max(test) == dim_one.SlackVar(index))

This is because the maximum value of zeros and one non-zero value is the same as the sum of that array.

This is the end of dim_one code. To summarize, dim_one’s slack value at each node is equal to the maximum of either the slack value at the prior node, or the data value at the current node. By recursion, this means that at each point in the path, the current index’s slack value holds the highest data value seen so far in the path. If the first node visited has the highest value, then that value is carried along through the sequence of slack values. If the last node visited has the highest value, then that final node’s slack will contain the maximum data value.

The slack value by itself is not useful to the optimization. The optimizer API calls only deal with the CumulVar() values, not the SlackVar() values. Because the slack at each point could equal the maximum value, the final CumulVar() of dim_one for each path could be ridiculously large. So in order to actually use the slack variable at the final node which holds the quantity we want to use, we need to create a second dimension, dim_two.

Dim Two stores the net cost of the route

The second dimension is created and initialized just like the first dimension.

# Dimension Two will be used to store the max node value in the route end node
# CumulVar so we can use it as an objective cost.
routing.AddConstantDimensionWithSlack(
    0,  # transit 0
    42 * 16,  # capacity: be able to have PEAK value in CumulVar(End)
    42,  # slack_max: to be able to store peak in slack
    True,  #  Fix StartCumulToZero YES here
    'Two')
dim_two = routing.GetDimensionOrDie('Two')

   ...

# force depot Slack to be value since we don't have any predecessor...
for v in range(manager.GetNumberOfVehicles()):
    start = routing.Start(v)
    ...
    dim_two.SlackVar(start).SetValue(data['value'][0])

The point of the dimension is to use the last slack value of dim_one as the actual value of the last node of dim_two, thereby setting the cost of the route. Once again, instead of directly assigning a value to the CumulVar of the last node, the value will be assigned to the slack of the second to last node (last real delivery), which will set the final CumulVar value.

The only complication is determining the actual last delivery node for each route. This is accomplished with the following loop:

for node in range(1, 17):
    index = manager.NodeToIndex(node)
    values = []
    for v in range(manager.GetNumberOfVehicles()):
        next_index = routing.End(v)
        cond = routing.NextVar(index) == next_index
        value = dim_one.SlackVar(index)
        values.append((cond * value).Var())
    solver.Add(solver.Sum(values) == dim_two.SlackVar(index))

This loop iterates over each of the actual delivery nodes. For each node, a second loop is performed over each of the vehicles. If the node is the last node visited on the route, then that means, by definition, that the next node will be the End(v) location for vehicle v. That is what the zero-one condition variable is here.

next_index = routing.End(v)
cond = routing.NextVar(index) == next_index

If the condition is true (cond == 1), then that means the slack variable we care about is the one for the current index. If the condition is false, cond==0, and the slack value of the current index is meaningless. As before, the slack variable is multiplied by the condition, and the expression’s value is added to the array of values.

values.append((cond * value).Var())

Finally, after looping over all of the available vehicles, the actual slack value of the index in question is assigned to be the sum of the values array.

solver.Add(solver.Sum(values) == dim_two.SlackVar(index))

This is almost exactly like the assignment of the SlackVar in the recursive loop, with one important difference. If the node in question is not the last node for any of the routes, then it will always be the case that the SlackVar(index) will be zero, because all of the cond values will be zero. That means that the value of dim_two.CumulVar() at every node in every route in the solution will be zero. The only non-zero value will be the final, vehicle end node. That is, only:

dim_two.CumulVar(routing.End(v)) > 0

The final piece of the puzzle is to insert the value of dim_two into the objective function. This is done by using the API call:

# 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, 4200)

The penalty 4200 was used without any notes, but I did notice that it was added in a later revision according to the git log. I assume that a value of 1 was determined to be too small.

Final notes

The last entry in the forum thread on this topic from the original poster Guy, mentioned that an actual solver based on this idea did not perform so well. 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

I was going to explore this issue in this post, but I’ve delayed pushing this out for too long now. Instead I will leave that for a future post. I also want to explore how this can be implemented using the CP-SAT solver.