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 `IntVar`

s.

```
/// 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.