This post explores a simple example in the OR-Tools CP-SAT examples that shows how to use the cumulative constraint. I have never used this constraint, but I am wondering if it can be adapted to the bus scheduling examples to track the per-route drive time before needing to schedule a break. So I’m going to do my usual line-by-line examination of the example’s source code.
The cumulative variable example
In the OR-Tools source code under the CP-SAT code, there is a sample of how to use the cumulative constraint. I was really just searching for examples of using the interval constraint when I came across this program. I scanned through it, and had no idea how it worked or even what it was supposed to be doing, so I figured it would be a good candidate for a blog post write up.
What do the docs say
The doc string for the cumulative constraint call is as follows:
Adds Cumulative(intervals, demands, capacity).
This constraint enforces that:
for all t:
sum(demands[i]
if (start(intervals[i]) <= t < end(intervals[i])) and
(intervals[i] is present)) <= capacity
Args:
intervals: The list of intervals.
demands: The list of demands for each interval. Each demand must be >= 0.
Each demand can be a 1-var affine expression (a * x + b).
capacity: The maximum capacity of the cumulative constraint. It can be a
1-var affine expression (a * x + b).
So what does that mean? If there are some intervals, say N of them, and some demands, also N of them, then this constraint makes sure that the sum of the demands linked with the active intervals is less than the capacity expression.
Maybe another way to think about it is to take a slice of time t, and then
look at all of the intervals that are active at that time. Each interval has
some related demand. If you sum up the demands for all of the active intervals
that intersect with the slice of time t, then this constraint forces that sum
to be less than whatever the capacity expression is.
My own small example
So I have a problem where a warehouse has two loading docks, and so I don’t want to schedule more than two loading jobs at the same time. I’m using a multiple circuit constraint to constrain my vehicle trips, and I really don’t want to faff about with trying to work something into that. I suppose it might be possible to create two locations for the warehouse, and allow either (but not both) to be visited by the warehouse-bound movements. But I’d still have to work out an approach to ensuring that none of the visits to either of the warehouses are overlapping in time.
Instead, I’m trying out the cumulative constraint. All of the warehouse-based tasks are easily identified by their location (“warehouse”) and the task (“load” or “unload”). So I hacked out some code
warehouse_intervals = []
for id in node_interval_details_lookup.values():
if not (
id.event_type in [
EventType.load, EventType.unload
] and
id.location_type == LocationType.warehouse
):
continue
warehouse_intervals.append(id.interval)
loading_docks = Options.number_of_loading_docks
model.add_cumulative(
intervals = warehouse_intervals,
demands = [1 for _i in warehouse_intervals],
capacity = loading_docks
)
And as sure as Bob’s your uncle and Sadie’s your aunt, there are at most two simultaneous warehouse intervals in any valid solution. So easy!
The example program does something trickier
So much for my trivial application. The example program I linked to above gives better insight into some non-trivial things the cumulative constraint is supposed to do. Specifically, it shows how to use the constraint to enforce a minimum level of coverage. First, the data define the maximum and minimum loads for each two-hour interval in a day. (Note that the data are defined as strings so as to read them into Pandas dataframes.)
max_load_str: str = """
start_hour max_load
0 0
2 0
4 3
6 6
8 8
10 12
12 8
14 12
16 10
18 6
20 4
22 0
"""
min_load_str: str = """
start_hour min_load
0 0
2 0
4 0
6 0
8 3
10 3
12 1
14 3
16 3
18 1
20 1
22 0
"""
Then the tasks are defined as having a duration, a load, and a priority as follows:
tasks_str: str = """
name duration load priority
t1 60 3 2
t2 180 2 1
t3 240 5 3
t4 90 4 2
...
"""
The goal is to schedule the tasks such that at all times the minimum and maximum load constraints are met. If some loads must be dropped, then the idea is to keep those with the highest priority.
Then the program creates the variables it needs.
# Variables
starts = model.new_int_var_series(
name="starts",
lower_bounds=0,
upper_bounds=horizon - tasks_df.duration,
index=tasks_df.index,
)
performed = model.new_bool_var_series(name="performed", index=tasks_df.index)
intervals = model.new_optional_fixed_size_interval_var_series(
name="intervals",
index=tasks_df.index,
starts=starts,
sizes=tasks_df.duration,
are_present=performed,
)
The starts are the starting times of each task, performed is a list of booleans
whether at task is performed or not, and intervals is a list of optional
intervals (with start times set by starts and presence set by performed).
Set the maximum constraint
So far, to me at least, everything seems super obvious. But then the next chuck of code is a mystery:
time_period_max_intervals = model.new_fixed_size_interval_var_series(
name="time_period_max_intervals",
index=max_load_df.index,
starts=max_load_df.start_hour * minutes_per_hour,
sizes=minutes_per_hour * 2,
)
time_period_max_heights = max_load - max_load_df.max_load
# Cumulative constraint for the max profile.
model.add_cumulative(
intervals.to_list() + time_period_max_intervals.to_list(),
tasks_df.load.to_list() + time_period_max_heights.to_list(),
max_load,
)
First is the creation of the max intervals. That part is pretty clear—it is
just creating one fixed interval per minute. It is the next two blocks that confuse
me. The line for time_period_max_heights makes a variable that is the
difference between a load’s maximum value and the absolute maximum value of the
loads. So it is zero if the load is at the maximum. As I write this, I have no
idea why this is being done. So try it out and see, right?
So if you print out the variable time_period_max_heights it is a dataframe
that looks like this:
max heights time period is
0 12
1 12
2 9
3 6
4 4
5 0
6 4
7 0
8 2
9 6
10 8
11 12
The next block sets up the cumulative constraint. The intervals for the constraint are a combination of the intervals for the tasks, and the intervals for the timer period maximum values. This list looks like this:
intervals[t1](start = starts[t1], size = 60, end = starts[t1] + 60, is_present = performed[t1]),
intervals[t2](start = starts[t2], size = 180, end = starts[t2] + 180, is_present = performed[t2]),
intervals[t3](start = starts[t3], size = 240, end = starts[t3] + 240, is_present = performed[t3]),
intervals[t4](start = starts[t4], size = 90, end = starts[t4] + 90, is_present = performed[t4]),
intervals[t5](start = starts[t5], size = 120, end = starts[t5] + 120, is_present = performed[t5]),
intervals[t6](start = starts[t6], size = 300, end = starts[t6] + 300, is_present = performed[t6]),
intervals[t7](start = starts[t7], size = 120, end = starts[t7] + 120, is_present = performed[t7]),
intervals[t8](start = starts[t8], size = 100, end = starts[t8] + 100, is_present = performed[t8]),
intervals[t9](start = starts[t9], size = 110, end = starts[t9] + 110, is_present = performed[t9]),
intervals[t10](start = starts[t10], size = 300, end = starts[t10] + 300, is_present = performed[t10]),
intervals[t11](start = starts[t11], size = 90, end = starts[t11] + 90, is_present = performed[t11]),
intervals[t12](start = starts[t12], size = 120, end = starts[t12] + 120, is_present = performed[t12]),
intervals[t13](start = starts[t13], size = 250, end = starts[t13] + 250, is_present = performed[t13]),
intervals[t14](start = starts[t14], size = 120, end = starts[t14] + 120, is_present = performed[t14]),
intervals[t15](start = starts[t15], size = 40, end = starts[t15] + 40, is_present = performed[t15]),
intervals[t16](start = starts[t16], size = 70, end = starts[t16] + 70, is_present = performed[t16]),
intervals[t17](start = starts[t17], size = 90, end = starts[t17] + 90, is_present = performed[t17]),
intervals[t18](start = starts[t18], size = 40, end = starts[t18] + 40, is_present = performed[t18]),
intervals[t19](start = starts[t19], size = 120, end = starts[t19] + 120, is_present = performed[t19]),
intervals[t20](start = starts[t20], size = 60, end = starts[t20] + 60, is_present = performed[t20]),
intervals[t21](start = starts[t21], size = 180, end = starts[t21] + 180, is_present = performed[t21]),
intervals[t22](start = starts[t22], size = 240, end = starts[t22] + 240, is_present = performed[t22]),
intervals[t23](start = starts[t23], size = 90, end = starts[t23] + 90, is_present = performed[t23]),
intervals[t24](start = starts[t24], size = 120, end = starts[t24] + 120, is_present = performed[t24]),
intervals[t25](start = starts[t25], size = 300, end = starts[t25] + 300, is_present = performed[t25]),
intervals[t26](start = starts[t26], size = 120, end = starts[t26] + 120, is_present = performed[t26]),
intervals[t27](start = starts[t27], size = 100, end = starts[t27] + 100, is_present = performed[t27]),
intervals[t28](start = starts[t28], size = 110, end = starts[t28] + 110, is_present = performed[t28]),
intervals[t29](start = starts[t29], size = 300, end = starts[t29] + 300, is_present = performed[t29]),
intervals[t30](start = starts[t30], size = 90, end = starts[t30] + 90, is_present = performed[t30]),
time_period_max_intervals[0](start = 0, size = 120, end = 120),
time_period_max_intervals[1](start = 120, size = 120, end = 240),
time_period_max_intervals[2](start = 240, size = 120, end = 360),
time_period_max_intervals[3](start = 360, size = 120, end = 480),
time_period_max_intervals[4](start = 480, size = 120, end = 600),
time_period_max_intervals[5](start = 600, size = 120, end = 720),
time_period_max_intervals[6](start = 720, size = 120, end = 840),
time_period_max_intervals[7](start = 840, size = 120, end = 960),
time_period_max_intervals[8](start = 960, size = 120, end = 1080),
time_period_max_intervals[9](start = 1080, size = 120, end = 1200),
time_period_max_intervals[10](start = 1200, size = 120, end = 1320),
time_period_max_intervals[11](start = 1320, size = 120, end = 1440)]
I think I am beginning to get a clue what is going on. Rather than creating a
list of maximum values for each of the possible time slots (one per minute for
an entire day), the goal here is to use a single constant (max_load) as the
maximum capacity in the cumulative constraint. This list of intervals has all
of the actual task intervals we care about, plus a list of dummy fixed intervals
that eat up the excess between the two-hour time period’s allowed maximum and
the true maximum load.
Looking again a the add_cumulative call:
# Cumulative constraint for the max profile.
model.add_cumulative(
intervals.to_list() + time_period_max_intervals.to_list(),
tasks_df.load.to_list() + time_period_max_heights.to_list(),
max_load,
)
At any time, if all of the intervals are active, that demand would be much too
high, but even so, they are bumped still higher by the time period max heights.
So take time 8. At that time the maximum load is 8. In the complemented list
I copied above, the 4th slot in the list (which is for time 8) has a
time period max height of 4. The max_load is 12. So whatever loads are
active at time 8, plus a base level of 4, cannot be more than 12. Which is
another way of saying that at time 8, the max load will be 8, which is what we
want.
Similarly, at time 0, the max load is 0, and so the corresponding
time_period_max_height is 12. That means that none of the tasks can be
active at time 0, because doing so would exceed the limit of 12, given that we
have that fixed baseline of 12 units at time 0.
Set the minimum constraint
And finally, the last bit of trickery in this example program is how to set up a
required minimum load. Without having seen this example, I probably would have
no idea how to tackle a minimum load type constraint with the add_cumulative
call. I probably would have just gone through each possible time step and
written something like:
for i in timesteps:
model.add(timestep_staffing_level[i] > minimum_staffing_level[i])
That would require all sorts of contortions to link it back to the actual intervals being performed.
The example program uses the add_cumulative call by building complement
intervals. For each task, the task can be performed, or not, and if it is
performed, there is a period of time before and after when it is not performed.
When it is not being performed, it contributes nothing to meeting the minimum
load requirements.
The starting time for the task is the value start and the end is start + duration. If the task is performed, then it contributes zero to the minimum
required load from time 0 to the time start; it contributes its load value
from start to time start + duration; and finally contributes zero from time
start + duration to the horizon time.
However, remember that the cumulative constraint is a maximum limit. So we need to take the complement of the above values. If a task is not performed, it provides nothing to the minimum, or conversely, it contributes its entire load towards not satisfying the minimum. By looking at it from this complemented point of view, the solver can see which tasks would be most advantageous to select as active in order to meet the minimum demand requirement.
With that background out of the way, now let’s look at the actual code. Again, I’m writing this post live as I figure out what is going on, so although I’ll probably rewrite things, some code bits are a complete surprise to me.
First up, the program sets up the complemented intervals in which the tasks are not performed.
# Set up complemented intervals (from 0 to start, and from start + size to
# horizon).
prefix_intervals = model.new_optional_interval_var_series(
name="prefix_intervals",
index=tasks_df.index,
starts=0,
sizes=starts,
ends=starts,
are_present=performed,
)
suffix_intervals = model.new_optional_interval_var_series(
name="suffix_intervals",
index=tasks_df.index,
starts=starts + tasks_df.duration,
sizes=horizon - starts - tasks_df.duration,
ends=horizon,
are_present=performed,
)
No surprises there aside from one item. If a task is performed, these are the times when it is not active. If a task is not performed, then these intervals are also not performed.
Next up is creating the minimum load periods, just as with the maximum load, in which fixed intervals are created in blocks to represent the minimum load rules.
time_period_min_intervals = model.new_fixed_size_interval_var_series(
name="time_period_min_intervals",
index=min_load_df.index,
starts=min_load_df.start_hour * minutes_per_hour,
sizes=minutes_per_hour * 2,
)
time_period_min_heights = min_load_df.min_load
The min_load_df looks like:
min load df is
start_hour min_load
0 0 0
1 2 0
2 4 0
3 6 0
4 8 3
5 10 3
6 12 1
7 14 3
8 16 3
9 18 1
10 20 1
11 22 0
And time_period_min_heights is:
0 0
1 0
2 0
3 0
4 3
5 3
6 1
7 3
8 3
9 1
10 1
11 0
Next comes some crazy trickery. First, recall that we’re taking the complement
of things, so if all of the tasks are performed at exactly the same time, and
therefore have unperformed regions at the exact same times, then the magnitude
of the unperformed prefix and suffix period can be as high as the sum of all of
the active loads. So the variable sum_of_demands is calculated first.
# We take into account optional intervals. The actual capacity of the min load
# cumulative is the sum of all the active demands.
sum_of_demands = sum(tasks_df.load)
Next, a complement_capacity variable is created that can range from 0 to this
maximum value.
complement_capacity = model.new_int_var(0, sum_of_demands, "complement_capacity")
And now the tricky bit. A new constraint is added that is what above I said would require me to do lots of complicated contortions to get to work right:
model.add(complement_capacity == performed.dot(tasks_df.load))
What on earth is performed.dot() doing? It looks like matrix math cribbed from
Pandas, because I remember dot products being a thing back when I was getting my
undeserved A in linear algebra as an undergraduate. So opening up the Pandas
documentation for
Series.dot(), it says indeed that dot() is the dot
product of a series and “other”, where “other” can be another series or a
dataframe, as long as both share the same index.
Again, from my undeserved A in linear algebra, I recall that the dot product is the row (hold my fingers out horizontally) times a column (rotate my hand to the right) with all of the multiplied elements summed. Dot products reduce dimension, so two vectors becomes a scalar, two matrices become a vector, etc.
So here it is the row performed which is a sequence of booleans, times
the load of each task, summed. So easy enough, it is just the total load of the
actual tasks that are performed.
It is also easy enough to print the resulting constraint:
a = model.add(complement_capacity == performed.dot(tasks_df.load))
print(a.proto)
linear {
vars: 30
vars: 31
vars: 32
vars: 33
vars: 34
vars: 35
vars: 36
vars: 37
vars: 38
vars: 39
vars: 40
vars: 41
vars: 42
vars: 43
vars: 44
vars: 45
vars: 46
vars: 47
vars: 48
vars: 49
vars: 50
vars: 51
vars: 52
vars: 53
vars: 54
vars: 55
vars: 56
vars: 57
vars: 58
vars: 59
vars: 60
coeffs: -3
coeffs: -2
coeffs: -5
coeffs: -4
coeffs: -3
coeffs: -3
coeffs: -1
coeffs: -5
coeffs: -2
coeffs: -5
coeffs: -4
coeffs: -3
coeffs: -3
coeffs: -1
coeffs: -5
coeffs: -4
coeffs: -8
coeffs: -3
coeffs: -5
coeffs: -3
coeffs: -2
coeffs: -5
coeffs: -4
coeffs: -3
coeffs: -3
coeffs: -1
coeffs: -5
coeffs: -2
coeffs: -5
coeffs: -4
coeffs: 1
domain: 0
domain: 0
}
This constraint is the linear combination (read sum) of the variables times their coefficients, which must all add up to 0 because that’s the way the CP-SAT solver enforces the equality constraint—it subtracts the RHS of an equality from the LHS and makes them add up to zero.
Finally, the tricky code concludes with the cumulative constraint for the minimum itself:
# Cumulative constraint for the min profile.
model.add_cumulative(
prefix_intervals.to_list()
+ suffix_intervals.to_list()
+ time_period_min_intervals.to_list(),
tasks_df.load.to_list()
+ tasks_df.load.to_list()
+ time_period_min_heights.to_list(),
complement_capacity,
)
Even though I’ve been typing my thoughts, copying code, running code and looking at the intermediate results, I still look at this and think “but where’s the minimum”? That’s why I call this trickery. The minimum is there, it is just hidden a little bit.
First, the intervals of the cumulative constraint contain the prefix and suffix periods. So these are the times when a performed interval are not active. Skipping ahead to the demands part of the API call, these unperformed periods are paired with the actual loads of the performed tasks.
The the last chunk of intervals added to the cumulative constraint are the
time_period_min_intervals that were created above. These add the required
minimum load to the total time period. So what does that do, what does that
mean?
Take some time, say time 8 (so it is interesting). Suppose all of the tasks
start at time 10, so none of them are active at time 8. Well this means that
all of prefix intervals have a true value, and all of the loads are active at
time 8 in this complemented land. But we’ve got a problem. We also have the
minimum load value at hour 8 of 3 added to that list, and we have that cool
dot product constraint that says that the maximum value can only be all of the
active loads. Instead we have all of the active loads plus 3! So in order to
meet the constraint, the solver needs to drop 3 from the load at 8. It could
choose to make a task inactive (set the presence literal to False), but ooh,
then the cool dot product constraint would kick in and say, no, you have removed
three, but the goalposts also moved three, and so you’re still three over!
Instead what it must do is activate some task so that it is active at 8, and
inactive elsewhere. That would keep the presence bit set, keep the goalposts
where they are, but also remove 3 load units from the complemented
prefix_intervals.
All in all, this is a pretty cool example and there are lots of things to learn from it.