My foray into discrete finite automata from the last few blogs has led to exploring more of the examples in the OR-Tools source code. Last Friday I was comparing two different approaches to the bus scheduling problem. One ran in a reasonable time on two easy problems but didn’t finish on the medium or large problems, and the other was infeasible on the easy problem, and then took forever on the middle problem. So I’m kicking out one or two (or six or seven?) blog posts to look at these programs. This first post in the series will introduce both programs, and then look in detail at bus_driver_scheduling_sat.py.

Two bus scheduling problems

In the OR-Tools source code, there is a big list of examples for python here. I was looking through the examples at different scheduling approaches to see if I could learn any new tricks. Business is slow right now, so I may as well improve my skills, right?

I came across two programs for bus driver scheduling. I’m not sure what I was expecting, but when I looked at the programs, neither of them did what I would have expected.

bus_driver_scheduling_flow_sat

The first one I looked at was called bus_driver_scheduling_flow_sat. I opened that one first because I’d never seen using a flow-based model to optimize scheduling. I had no idea how that might be done. The comment in the code states:

# We are going to build a flow from a the start of the day to the end
# of the day.
#
# Along the path, we will accumulate driving time, accrued time since the
# last break, and total working time.

The code that follows that comment looks like the standard approach to how one creates a circuit constraint. A network of nodes is defined, one for each shift, and then arcs are created from all nodes to all other nodes. The arcs are active if the node at the tail of the arc (the source) is followed by the node at the head of the arc (the sink). Further, there are the same number of unconnected sources and sinks as there are drivers. This is a list of directed acyclic paths through all of the nodes, one path per driver.

I come from an engineering background (literally, my undergraduate degree from HMC says “Engineering”), not a programming background. So to my mind, flow means the movement of a fluid, not a network of nodes and arcs. I expected to see analogues of flow velocity, pressure, and density, and at the very least to have the arcs have some property that relates to capacity. In this code, flow is used in the computer science sense. The capacity is just 1 vehicle per link, and a limited number of links per route.
So I guess, if you squint and look sideways, this sort of describes “flow” through the shifts. But there isn’t any notion of min-cut/max-flow, or bottlenecks, or anything like that.

Moving on to the code itself, the program creates a network of all shifts to all shifts, with arcs between them allowing just one active inbound and outbound arc per node (shift). It then manually enforces a directed acyclic graph constraint, with one path per driver out of the source and into the sink. It does not call or use the OR-Tools circuit constraint.

# Create dag constraint.
for shift in range(num_shifts):
    model.Add(sum(outgoing_literals[shift]) == 1)
    model.Add(sum(incoming_literals[shift]) == 1)

# Num drivers
num_drivers = model.NewIntVar(min_num_drivers, min_num_drivers * 3, "num_drivers")
model.Add(sum(incoming_sink_literals) == num_drivers)
model.Add(sum(outgoing_source_literals) == num_drivers)

The real problem with this example though is that it doesn’t appear to work properly. I ran it on my laptop using the smallest data set (option 1, 50 shifts) and it was infeasible.

$ python bus_driver_scheduling_flow_sat.py --instance 1
----------- first pass: minimize the number of drivers
Bus driver scheduling
  num shifts = 50
  total driving time = 2355 minutes
  min num drivers = 5
  min start time = 318
  max end time = 1248

Starting CP-SAT solver v9.15.6755

...

Preloading model.
#Bound   0.05s best:inf   next:[5,15]     initial_domain
#Model   0.05s var:599/599 constraints:1337/1337

Starting search at 0.05s with 22 workers.

...

#Done    0.08s reduced_costs [loading]
...
LRAT_status: NA
CpSolverResponse summary:
status: INFEASIBLE
objective: NA
best_bound: NA

The second problem instance (with 200 shifts) generated some solutions, but then got stuck hunting, as there is no time limit set in the code:

$ python bus_driver_scheduling_flow_sat.py --instance 2
----------- first pass: minimize the number of drivers
Bus driver scheduling
  num shifts = 200
  total driving time = 7793 minutes
  min num drivers = 15
  min start time = 270
  max end time = 1512

Starting CP-SAT solver v9.15.6755

...

#Model   1.48s var:6534/6566 constraints:17880/17956
#Bound   1.61s best:inf   next:[16,45]    reduced_costs
#1       2.03s best:34    next:[16,33]    fs_random_quick_restart
#2       2.06s best:33    next:[16,32]    quick_restart_no_lp
#3       2.09s best:31    next:[16,30]    quick_restart_no_lp
#4       2.11s best:30    next:[16,29]    quick_restart_no_lp
#5       4.59s best:29    next:[16,28]    graph_arc_lns (d=4.36e-01 s=153 t=0.10 p=0.50 stall=0 h=base)

...

It showed promise at first, but then it just got stuck.

Finally, I ran the program on problem instance three, with 1356 shifts to schedule. I ran it, then I made tea and unloaded the dishwasher.

$ python bus_driver_scheduling_flow_sat.py --instance 3
----------- first pass: minimize the number of drivers
Bus driver scheduling
  num shifts = 1356
  total driving time = 55483 minutes
  min num drivers = 103
  min start time = 258
  max end time = 1526

Starting CP-SAT solver v9.15.6755

...

Preloading model.
#Bound  31.97s best:inf   next:[103,309]  initial_domain
#Model  32.24s var:270200/270200 constraints:798726/798726

Starting search at 32.55s with 22 workers.

...


#Model  96.52s var:269678/270200 constraints:797846/798726
#1     234.38s best:167   next:[103,166]  fs_random_no_lp
#2     249.17s best:166   next:[103,165]  rnd_cst_lns (d=5.00e-01 s=1800 t=0.10 p=0.00 stall=0 h=base) [hint]
#3     267.41s best:165   next:[103,164]  graph_cst_lns (d=5.00e-01 s=1813 t=0.10 p=0.00 stall=0 h=base) [hint]
#4     286.67s best:155   next:[103,154]  quick_restart_no_lp
...
#12    854.02s best:147   next:[103,146]  lb_relax_lns_bool (d=6.47e-02 s=8263 t=0.50 p=0.38 stall=0 h=base) [hint]
#13    855.20s best:146   next:[103,145]  lb_relax_lns_bool (d=6.47e-02 s=8262 t=0.50 p=0.38 stall=0 h=base)

So something is wrong in the flow example code and I was tempted to take a deeper look, but I moved on to the sat program first. Maybe after writing these posts I can fix it and submit a patch!

bus_driver_scheduling_sat

The other example for bus scheduling is the program bus_driver_scheduling_sat. This one appears at first glance to be similar to the “flow” version. This program also creates a network from all of the shifts to all other shifts, using literals to capture whether an arc is active or not, and tallying up the incoming arcs and outgoing arcs at each node and making sure there is one of each. But unlike the “flow” program, in this formulation the all to all generation of nodes and arcs is done for each possible driver. This seems like it is going to be terribly slow with lots of work for the solver to do. I switched over to my terminal to run the program, fully expecting it to be even slower than the “flow” version and to be able to go and do some laundry while it ran. In fact, this program runs lickety split on the embedded “tiny” data, and in okay time on the “small” data set.

Line by line analysis

I will eventually come back to the flow version of the bus solver in some future post. For the remainder of this post I will just concentrate on the non-flow version. This bus driver scheduling program is set up for two passes of the same function over the data. The first pass has a goal of minimizing the numbers of drivers, and the second pass has a goal to minimize the total “working” time over all of the drivers.

def bus_driver_scheduling(minimize_drivers: bool, max_num_drivers: int) -> int:

    ...

def main(_):
    """Optimize the bus driver allocation in two passes."""
    print("----------- first pass: minimize the number of drivers")
    num_drivers = bus_driver_scheduling(True, -1)
    if num_drivers == -1:
        print("no solution found, skipping the final step")
    else:
        print("----------- second pass: minimize the sum of working times")
        bus_driver_scheduling(False, num_drivers)

In the first pass we have no idea how many drivers might be needed. So the best approach is to guess high, and let the solver work its way back to the minimum number required. The input data contains only shifts, not the number of drivers, but it does contain parameters describing how much work each driver can do. A small bit of manual calculation is used to generate an estimate of too many drivers.

total_driving_time = sum(shift[5] for shift in shifts)
min_num_drivers = int(math.ceil(total_driving_time * 1.0 / max_driving_time))
num_drivers = 2 * min_num_drivers 

The first step is to compute the total driving time from serving all of the shifts. Next, because there is a hard maximum driving time per driver, the absolute minimum number of drivers is equal to if each driver can drive his or her maximum amount of time. Of course this is unlikely to be possible, as the shifts all overlap each other in erratic ways. Instead, each driver will likely drive less than that maximum amount. So here, knowing how the data looks, the programmer decided that double the minimum number of drivers is a comfortable value for “too many” drivers. If that guess is too low, then the solver run will be infeasible. If the guess is too high, then we will waste some time eliminating spare drivers. The decision to use 2 * min_num_drivers is entirely dependent upon the real data. If the range of possible inputs are known to be highly erratic, then this factor might need to be larger. Or perhaps the minimum expected driver time might be used, just to make sure that this initial guess at the maximum number of drivers results in a feasible solution. I find that it is usually faster to drop extra drivers etc inside of the solver, than it is to set up a solver run, have it fail, try again, etc. etc.

With the estimated drivers in hand (which hopefully has more drivers than needed!), the program next does the dance over each shift to create each driver’s network. The comment in the code describes it best:

# For each driver and each shift, we store:
#   - the total driving time including this shift
#   - the accrued driving time since the last 30 minute break
# Special arcs have the following effect:
#   - 'from source to shift' sets the starting time and accumulates the first
#      shift
#   - 'from shift to end' sets the ending time, and fills the driving_times
#      variable
# Arcs between two shifts have the following impact
#   - add the duration of the shift to the total driving time
#   - reset the accumulated driving time if the distance between the two
#     shifts is more than 30 minutes, add the duration of the shift to the
#     accumulated driving time since the last break otherwise

So there will be a lot going on. The arcs between shifts are used to carry state. Because the shifts occur at fixed times, an active arc has a known starting time, and a known ending time. Thus the arcs allow us not only to link up shifts, but also to fix what time things are happening, and therefore create longer periods between shifts that allow drivers to take their required breaks. Let’s break down the code chunk by chunk.

Stepping through the optimization function for bus_driver_scheduling_sat.py

First comes the initialization. Because this is just a demo program, all of these variables are hard-coded. In a real problem they would be parameters or input values.

def bus_driver_scheduling(minimize_drivers: bool, max_num_drivers: int) -> int:
    ...
    num_shifts = len(shifts)

    # All durations are in minutes.
    max_driving_time = 540  # 8 hours.
    max_driving_time_without_pauses = 240  # 4 hours
    min_pause_after_4h = 30
    min_delay_between_shifts = 2
    max_working_time = 720
    min_working_time = 390  # 6.5 hours
    setup_time = 10
    cleanup_time = 15

    # Computed data.
    total_driving_time = sum(shift[5] for shift in shifts)
    min_num_drivers = int(math.ceil(total_driving_time * 1.0 / max_driving_time))
    num_drivers = 2 * min_num_drivers if minimize_drivers else max_num_drivers
    min_start_time = min(shift[3] for shift in shifts)
    max_end_time = max(shift[4] for shift in shifts)

Note that the variable num_drivers is dependent upon whether the function is being called with minimize_drivers set to True or False.

num_drivers = 2 * min_num_drivers if minimize_drivers else max_num_drivers

If it is set to false, then the assumption is that the input variable max_num_drivers has already been discovered in a previous pass, so it is used directly as the number of drivers.

Next the code creates a handful of containers to hold variables, and then starts the loop over the drivers. The first thing that is done with each driver is to create the driver’s time variables–the start time, end time, driving time, and working time–which will be determined by the solver run.

for d in range(num_drivers):
    start_times.append(
        model.new_int_var(min_start_time - setup_time, max_end_time, "start_%i" % d)
    )
    end_times.append(
        model.new_int_var(min_start_time, max_end_time + cleanup_time, "end_%i" % d)
    )
    driving_times.append(model.new_int_var(0, max_driving_time, "driving_%i" % d))
    working_times.append(
        model.new_int_var(0, max_working_time, "working_times_%i" % d)
    )

Then, for each driver, the shift variables are created. These are the total driving by the driver up to and including the shift, the driving up to and including the shift since the last break period, and a boolean indicating if this driver performs the shift.

    # Create all the shift variables before iterating on the transitions
    # between these shifts.
    for s in range(num_shifts):
        total_driving[d, s] = model.new_int_var(
            0, max_driving_time, "dr_%i_%i" % (d, s)
        )
        no_break_driving[d, s] = model.new_int_var(
            0, max_driving_time_without_pauses, "mdr_%i_%i" % (d, s)
        )
        performed[d, s] = model.new_bool_var("performed_%i_%i" % (d, s))

These variables are important, so at the risk of repeating myself, I’ll go over each one. First, for this driver and for each of the possible shifts, an IntVar called total_driving is created that can range from 0 to the max_driving_time allowed. This variable will contain the “driving state” of the driver if and when he or she completes shift s. This variable holds the state of the current driver’s total shift assignment up to this shift. All of the shifts up to shift s as well as the driving time of shift s will be included. Because the maximum value is equal to the maximum allowed driving time, this makes sure that no assignment can result in a drive time that exceeds the maximum allowed.

Similarly, the second variable no_break_driving is another IntVar that accumulates how long the driver has been driving. However, instead of the total time driving, it only adds up driving time since the last break up to and including shift s. It ranges from 0 to the value of max_driving_time_without_pauses. As with total_driving, this variable ensures that at no time can the driver drive more than the maximum allowed drive time without a break before going on another break (which resets the counter).

Finally, a simple boolean literal is created to track whether or not shift s is performed by driver d. All shifts must be performed, but for any reasonably sized problem, no single driver can perform all shifts. The nice thing about the list performed[d, s] is that you can set two constraints. If you hold the shift constant and sum up over all drivers, then you can require that the sum be exactly 1 to make sure one and only one driver performs the shift. If you hold the driver constant, then you have a boolean sequence of ones and zeros over all shifts that tells you how many shifts the driver performs, and lets you multiply by other features of the shift like the drive time or elapsed time to compute aggregate variables for each driver.

With all of these variables in hand, the double loop over the shifts can begin.

Sidebar: A network creates an ordering

As I was writing this up I realized that when I first started working on circuit constraints, it took me a long time to grasp the concept of why it is useful. But now, several years on, I am writing this up and completely skipping any discussion of why we need the network and the circuit constraint.

The big idea behind a network of nodes is to establish the order of things.

In a network, each driver has to start from some node (the source), end at some other node (the sink), and then can, theoretically, reach every node in one long sequence of arcs between those nodes. The arcs dictate a sequence. If the driver does nothing, there is still an arc from the source to the sink, and still the notion that the driver is first at the source node, and then at the sink node.

If you don’t care about the order of things, then it is much faster and easier to just skip creating a network and invoking the circuit constraint rules. In this case, with timed shifts, you can just use intervals and a no-overlap constraint. If you just cared about total driving time (which really is just filling up bags with items), you don’t even need intervals. You can just keep track of the drive time of each shift. But in this case, the shifts have a fixed starting and ending time, and so the valid assignment has to respect those values. So an interval works best.

In addition, because the shifts are timed, you can also generate the ordering of the nodes. However, and this is important, the order of nodes is not something that can be computed inside of the solver run. You can really only figure it out after the solver is finished. This means that the order of shifts is not a variable that the solver logic can use. If you use intervals, you can only know inside the solver logic if one shift comes before another. You can’t really know if a shift is the immediate predecessor of another shift unless their ending and starting times match up exactly. Further, you can’t know cumulative values, like how much driving has been done up to this point.

If we only had a requirement that drivers work a minimum number of hours, and do not drive or work more than a maximum number of hours, then using optional intervals that do not overlap would be fine. You can even handle some minimum time between shifts by adding a small buffer to the shift start and end times. In this problem, however, we also have a requirement that drivers must take breaks after a certain amount of driving. This constraint is more difficult to satisfy if you only know the shifts, but not their order. If you try to insert breaks after the shift sequences are assigned, then you will probably create infeasible assignments.

To make this obvious, suppose we have an unordered list of shifts for a driver that meets the aggregate constraints (max driving time, max work day, etc), and also incorporates rules fixed setup and cleanup times, minimum times between shifts, and the fact that shifts cannot overlap. With that list of shifts, we can sort them into the right order by just looking at the start and end times of each shift. Suppose a shift sequence looks like this:

shift   2 : 05:40 - 05:56   16   16
shift   3 : 06:06 - 06:51   45   61
shift   7 : 06:59 - 08:07   68  129
shift  13 : 08:11 - 09:41   90  219
shift  20 : 09:57 - 10:20   23  242
shift  23 : 11:00 - 11:10   10  252
shift  24 : 11:45 - 12:24   39  291

Now we’ll try to add in break times. Starting from the beginning of the assignment, summing up the incremental driving times, we need to insert a break after 240 minutes (4 hours) of work (not of elapsed clock time). So here the schedule starts with shift 2 at 5:40, and the incremental driving times hit 242 minutes at the end of shift 20. So between shift 13 and shift 20, we need to add a 30 minute break, from 9:41 to 10:11. But of course, that means the driver cannot start shift 20 on time, and the assignment is invalid.

shift   2 : 05:40 - 05:56   16   16
shift   3 : 06:06 - 06:51   45   61
shift   7 : 06:59 - 08:07   68  129
shift  13 : 08:11 - 09:41   90  219
**break** : 09:41 - 10:11
shift  20 : 09:57 - 10:20   23   23  *invalid!*
shift  23 : 11:00 - 11:10   10   33
shift  24 : 11:45 - 12:24   39   72

So the solver needs to know this, so that it doesn’t waste time creating invalid schedules. And that means that we need to set up a network that gives us the complete order of nodes.

You might just think to use a floating “break” interval that is 30 minutes long and that can occur anytime. I’ve seen some problem formulations, for example assigning tasks to hourly workers, where the break is specified as something like “a half hour break between 11AM and 1PM”. But that won’t work here. If a driver cannot drive more than 4 hours in a row without a 30 minute break, then that rule must be followed. In this example problem, we have a hard requirement that the driver cannot drive more than 240 minutes without a 30 minute break. That means we need to track the driving time as we go, and slot in the 30 minute break before we hit the 240 minutes of driving. To do that, we need to know the order of shifts inside the solver as we compute the best assignment of shifts.

That is the whole point of creating this network of “shifts to shifts” for each driver. It is expensive both computationally and in terms of RAM, but this network allows us to create variables that know the current order of shifts that are assigned to a driver. Without that network and without that knowledge, we can run faster and use less memory, but we’re going to generate invalid schedules, so what’s the point?

Starts and ends: From source to shift, shift to sink

The first question to answer is which node is visited first, and the second question is which node is visited last. Then we want to know, for each node, if it doesn’t move to the sink, then which node does it transition to? So, with those questions in mind, we are ready to begin the double loop over all of the nodes.

First the outer loop, and the definition of the current shift s and its duration.

for s in range(num_shifts):
    shift = shifts[s]
    duration = shift[5]

So the first thing that must be done is to determine whether this node is visited first. If it is, then there will be an active link from the starting source node to this shift for driver d. So that needs a literal to be created:

# Arc from source to shift.
#    - set the start time of the driver
#    - increase driving time and driving time since break
source_lit = model.new_bool_var("%i from source to %i" % (d, s))
outgoing_source_literals.append(source_lit)
incoming_literals[s].append(source_lit)
shared_incoming_literals[s].append(source_lit)

Here the source_lit is created and then is tracked by adding it to various containers. The list outgoing_source_literals holds all of the literals from the source node to somewhere else. In this case, if source_lit is true, then the driver moves from the source node to the shift s. Next, the source_lit is added to the list incoming_literals for the node s. This list tracks every link that moves from any other shift node and the source node into the node for s. If the lit is true, then the movement is performed. Because the driver can only perform any given shift at most once, only one of the literals in the list incoming_literals[s] can be true. If all of them are false, then the shift s is not performed by this driver. If one is true, then that driver performs shift s, and further you can know which shift (or source node) preceded shift s in this driver’s path.

Finally, the source_lit literal is added to the special list shared_incoming_literals[s]. This list tracks all incoming nodes for all drivers to shift s. This is how you can state that only one driver can perform a shift s by saying that at most one of the literals in the list can be true. Further, if you required that each shift must be performed, then you can set the constraint that exactly one of these literals must be true. That would mean that one and only one driver can perform shift s.

Next, the literal is used to set some truths about the start and end time variables.

model.add(start_times[d] == shift[3] - setup_time).only_enforce_if(
    source_lit
)
model.add(total_driving[d, s] == duration).only_enforce_if(source_lit)
model.add(no_break_driving[d, s] == duration).only_enforce_if(source_lit)
starting_shifts[d, s] = source_lit

The first constrain says that if the source lit is true (the first shift performed by the driver is shift s), then the start time of the driver must necessarily be early enough to perform the shift. In this case, that means that the start time of the driver must be equal to the start time of the shift minus the fixed setup time for each shift. This setup time might cover things like loading any required items for the shift, making a log entry, fueling up the vehicle, and so on.

Next, the variable for total_driving for driver d at node s is set equal to the duration of node s, but only if the source_lit is true. Again, if this literal is true, then s is the first node in the driver’s job, and so the total driving duration, which starts at zero by definition, can be set to the duration of the shift.

Similarly for the variable no_break_driving, the value is also set to the shift duration if the source_lit is true.

Finally the source_lit is added to the global list of starting_shifts for each driver d and each shift s. This list allows constraints to be placed on the starting shifts of each driver. For example, if you require that the sum of all of the literals in this list must be equal to the number of drivers, that would mean every driver must have at least one shift performed. If however you add to this list a literal that says the driver goes straight from the start node to the sink node, then that would allow some drivers to do nothing.

Next a similar set of instructions are coded up for transition from the shift s to the sink node, or the end of the driver’s tasks. This code is exactly like the above code, except it captures the movement out of the shift s to the end of the route, and sets constraints on the shift end time.

# Arc from shift to sink
#    - set the end time of the driver
#    - set the driving times of the driver
sink_lit = model.new_bool_var("%i from %i to sink" % (d, s))
outgoing_literals[s].append(sink_lit)
shared_outgoing_literals[s].append(sink_lit)
incoming_sink_literals.append(sink_lit)
model.add(end_times[d] == shift[4] + cleanup_time).only_enforce_if(sink_lit)
model.add(driving_times[d] == total_driving[d, s]).only_enforce_if(sink_lit)

As before, first a literal sink_lit is created that is true if the shift s is the last one performed by driver d. This literal is added to the list of outgoing_literals from s (at most one can be true), the shared_outgoing_literals for s (again at most one can be true), and the list of global list of incoming_sink_literals over all drivers. Finally, some truths about time are imposed. If you know that sink_lit is true and shift s is the last shift performed by d, then you know that the end time for driver d’s day must be the end time of the shift plus the fixed post-shift clean up time. Similarly, you know that the final total driving_time for driver d must be equal to the total_driving time for d at node s, because if s is the very last shift performed, then no additional driving time can be added on. Unlike with the start to s rules, there is no need to worry about the no_break_driving variable, because if shift s is last, then there is no longer any need to track when a break must be inserted into the list of shifts to perform.

Finally, before beginning the second loop over the shifts, you know two things that are true if this shift is performed, regardless of whether it is first, last, or somewhere in the middle.

# Node performed:
#    - add upper bound on start_time
#    - add lower bound on end_times
model.add(start_times[d] <= shift[3] - setup_time).only_enforce_if(
    performed[d, s]
)
model.add(end_times[d] >= shift[4] + cleanup_time).only_enforce_if(
    performed[d, s]
)

First, if the shift s is performed by driver d, then it must be true that the value of start_times for the driver d is less than or equal to the start time of the shift (shift[3]) less the required setup_time. In other words, if the driver performs a shift, then it must be true that the driver starts early enough to perform the shift! Second, you also know that the value of end_times for the driver d is greater than or equal to the ending time of the shift (shift[4]) plus the required cleanup_time. Again, if the shift is performed, then the driver must end his or her day after the shift is performed. These constraints are true and enforced if and only if the shift is performed. Importantly, while these constraints are redundant (we already dictate the exact equality for the start and end times when setting the first and last shifts), these conditional constraints are always true, and so help the solver to reach the right conclusions as it works through the variables and constraints on its hunt for the optimal assignment.

The inner loop over shifts

And now we come to the inner loop. This inner loop handles moving from some shift s to some other shift o. Conceptually, it is the same as handling moving from the source to s, or from s to the sink. You make a literal to represent the arc, and if the literal is true, then the arc is active. The literal is an exit literal from s, and an entrance literal to o. But the code does some other things too, so I’ll step through it bit by bit.

for o in range(num_shifts):
    other = shifts[o]
    delay = other[3] - shift[4]
    if delay < min_delay_between_shifts:
        continue

First up, as with any double loop, you want to have some early quit conditions right away. There is a rule that there must be a minimum delay between shifts, as well as an implicit rule that the end of shift s must come before the beginning of shift o. In either case, if the delay is less than the minimum delay (either because it is negative or too small), then we quit and move on to the next node.

If we haven’t quit, then it is time to create the literal representing movement from s to o, and save it in the various containers. One thing I learned from looking through the OR-Tools examples is that slotting away the same thing in multiple places is generally more convenient than creating one giant container and doing fancy, R-style indexing over that container. The solver run is going to eat up time and space, so there isn’t much point in minimizing repeated copies or multiple passes over the data in the model construction code.

lit = model.new_bool_var("%i from %i to %i" % (d, s, o))

# add arc
outgoing_literals[s].append(lit)
shared_outgoing_literals[s].append(lit)
incoming_literals[o].append(lit)
shared_incoming_literals[o].append(lit)

This next block increases the driving and no break driving time variables.

# Increase driving time
model.add(
    total_driving[d, o] == total_driving[d, s] + other[5]
).only_enforce_if(lit)

# Increase no_break_driving or reset it to 0 depending on the delay
if delay >= min_pause_after_4h:
    model.add(no_break_driving[d, o] == other[5]).only_enforce_if(lit)
else:
    model.add(
        no_break_driving[d, o] == no_break_driving[d, s] + other[5]
    ).only_enforce_if(lit)

The the increase of total_driving is just like what was done earlier for source to s and s to sink If the arc from s to o is active, then that means you should add the driving time in.

The code for the no_break_driving is somewhat different. First the delay is checked. If the time time between the end of shift s and the start of shift o is large enough, then the delay can count as an official break. Recall that the variable no_break_driving can only reach a maximum of max_driving_time_without_pauses:

no_break_driving[d, s] = model.new_int_var(
    0, max_driving_time_without_pauses, "mdr_%i_%i" % (d, s)
)

That means as the solver computes possible feasible solutions, it keeps track of how large no_break_driving is getting. If the break between the two shifts is large enough (larger than min_pause_after_4h), the the no_break_driving is reset, dropping down to just the drive time from shift o. If the break is not large enough, the drive time from shift o is added to the cumulative value of no_break_driving. As the cumulative value increases closer and closer the limit, the only feasible solutions will become those with an active arc to a node with a “large enough” break time.

There is one more bit of code inside the handler for shift o from shift s. This is to keep track of the delay between shifts:

# Cost part
delay_literals.append(lit)
delay_weights.append(delay)

If the link is active from s to o, then the literal lit will be true, or 1. In that case, the 1 multiplied by the corresponding delay weight value of delay will be positive (the value of delay, of course). If the arc is not used, and o does not follow s, then lit times delay will be zero. This means that you can use these values in the objective function to represent the total delay in the overall shift assignments. My minimizing the weighted sum of all of the literals, you in turn minimize the total delay of the assignment.

That ends the double loop over the shifts. Before finish up with the current driver d, a few more things need to be done. First the total working time for the driver is defined as the value of the end time minus the start time for that driver:

model.add(working_times[d] == end_times[d] - start_times[d])

Next, the code has to handle how to minimize the total number of drivers. Recall there are two ways to invoke this function. The first pass has minimize_drivers set to true and uses too many drivers. The second pass, with the minimum number of drivers in hand, attempts to optimize the schedule. So this next chunk code handles the minimize drivers case.

working = model.new_bool_var("working_%i" % d)
model.add(start_times[d] == min_start_time).only_enforce_if(~working)
model.add(end_times[d] == min_start_time).only_enforce_if(~working)
model.add(driving_times[d] == 0).only_enforce_if(~working)
working_drivers.append(working)
outgoing_source_literals.append(~working)
incoming_sink_literals.append(~working)
# Conditional working time constraints
model.add(working_times[d] >= min_working_time).only_enforce_if(working)
model.add(working_times[d] == 0).only_enforce_if(~working)

In the above code, a boolean is created to indicate if the driver is working or not. That indicator variable is used to define correct values for other variables. If the driver is not working (using the shorthand ~working), then the driver’s start_times[d] and end_times[d] are both set equal to min_start_time, and the total driving_times[d] is fixed at zero. In short, this boolean is the literal for an arc that goes directly from the source node to the sink node. Therefore, just like any other arc, it is added to the appropriate list of outgoing and incoming literals: outgoing_source_literals for the source node, and incoming_sink_literals for the sink node.

Finally, if the literal working is true, then the solver is required to make sure that the working_time[d] for the driver is greater than or equal to the min_working_time requirement. If the driver is not working (~working is true), then the working time of the driver is set to 0.

If the code is not trying to minimize the driver count, if minimize_drivers is false, then the above block is skipped and just this very simple code is done:

# Working time constraints
model.add(working_times[d] >= min_working_time)

Thus once we have determined the minimum number of drivers that can hit all of the shifts, then we no longer allow drivers to do nothing (there is no arc directly from the source to the sink). However, we still need to set the requirement that the working_times[d] for the driver is at least the required min_working_time.

Finally, before moving on to the next driver, the code does what it says is creating the circuit constraint.

# Create circuit constraint.
model.add_exactly_one(outgoing_source_literals)
for s in range(num_shifts):
    model.add_exactly_one(outgoing_literals[s])
    model.add_exactly_one(incoming_literals[s])
model.add_exactly_one(incoming_sink_literals)

These are just the constraints for this driver, not the global ones. The first requirement says that the driver, starting from the source, can only move to one node. That node can be directly to the sink (no work), or to one of the shift nodes. Then, for each shift node, only one of the incoming and outgoing nodes can be true. This means that no shift node can be visited twice. But what if the shift is not visited by the driver? After all, no driver can possibly visit all of the shifts, so if so how can the per-shift constraints be the add_exactly_one type rather than at most one? Well, inside of the shift setup, there is a literal called performed[d, s] that is true if the shift s is performed by driver d. Then there is a little bit of code I didn’t talk about:

incoming_literals[s].append(~performed[d, s])
outgoing_literals[s].append(~performed[d, s])

This is essentially a self-loop. If the driver is not assigned the shift, then performed[d, s] is false, and ~performed[d, s] is true. This is added to the incoming_literals[s] as well as the outgoing_literals[s] for the driver. So if the shift is not performed, then there is a loop that goes out from s and then back into s, isolating it from the rest of the route but still maintaining the add_exactly_one condition.

Global constraints and prep for solving

The final part of the function covers overall constraints on the model across all drivers and shifts. First, every shift must be performed, and only once:

for s in range(num_shifts):
    model.add_exactly_one(performed[d, s] for d in range(num_drivers))
    # Globally, each node has one incoming and one outgoing literal
    model.add_exactly_one(shared_incoming_literals[s])
    model.add_exactly_one(shared_outgoing_literals[s])

Next, a small hack is added for “symmetry breaking”. The problem is that all of the drivers are identical, and so the solver treats each driver as equally likely to drive the first shift. So what it will do is build out a complete candidate solution with every driver picking up the first shift just in case it turns out that one is better than the others. This is expensive and slow. It is much better to just pick a few drivers randomly to be preferred.

# Symmetry breaking

# The first 3 shifts must be performed by 3 different drivers.
# Let's assign them to the first 3 drivers in sequence
model.add(starting_shifts[0, 0] == 1)
model.add(starting_shifts[1, 1] == 1)
model.add(starting_shifts[2, 2] == 1)

The code is taking advantage of the fact that it is known that the first three shifts overlap with each other. In a general problem, this could probably be automated, by just checking how many shifts overlap the first shift. In fact, I think I’m going to add that right now because the large data set has lots of overlaps. My hack looks like this:

shift0 = shifts[0]
needed_unique_drivers = 1
for s in range(1, num_shifts):
    shift = shifts[s]
    # does it overlap shift0?
    delay = shift[3] - shift0[4]
    if delay >= min_delay_between_shifts:
        # done
        break
    # the first driver is still on shift0, cannot cover this shift
    needed_unique_drivers += 1

# Let's assign them to the first N drivers in sequence
for d in range(needed_unique_drivers):
    print(f"symmetry breaking: linking driver {d} with shift {d}")
    model.add(starting_shifts[d, d] == 1)

I reused the delay code. Maybe I should make it a function call. Anyway, it doesn’t make much difference. Problem instance 2 uses four drivers, not three, but it still takes an age on my laptop.

Back to the program analysis, there is another bit of code for symmetry breaking. If the mode of the solver run is to minimize the drivers, then it is easier for the solver if all of these unused drivers are pushed to the back. This lets the solver know to prefer using earlier drivers than later drivers.

if minimize_drivers:
    # Push non working drivers to the end
    for d in range(num_drivers - 1):
        model.add_implication(~working_drivers[d], ~working_drivers[d + 1])

This constraint works as follows. If driver d is working, then the implication is false, and there is no impact on driver d + 1. If driver d is not working, then ~working_drivers[d] is true (remember, ~ is the negation operator for booleans), and so the implication holds, and driver d + 1 also must not be working.

Another hack to speed things up is the insertion of some redundant constraints.

# Redundant constraints: sum of driving times = sum of shift driving times
model.add(cp_model.LinearExpr.sum(driving_times) == total_driving_time)
if not minimize_drivers:
    model.add(
        cp_model.LinearExpr.sum(working_times)
        == total_driving_time
        + num_drivers * (setup_time + cleanup_time)
        + cp_model.LinearExpr.weighted_sum(delay_literals, delay_weights)
    )

There aren’t any hard rules for when to use these, but they do tend to help. Here, we know before the solver runs that since all of the shifts must be served, the summation of all the assigned driving times must be equal to the total_driving_time. This stops the solver from having to figure that out. The other rule comes into play when drivers are not being minimized, because we know the number of drivers in advance.

The final bit of logic before creating and calling the solver is to set up the optimization goal.

if minimize_drivers:
    # minimize the number of working drivers
    model.minimize(cp_model.LinearExpr.sum(working_drivers))
else:
    # minimize the sum of delays between tasks, which in turns minimize the
    # sum of working times as the total driving time is fixed
    model.minimize(cp_model.LinearExpr.weighted_sum(delay_literals, delay_weights))

If the mode is to minimize drivers, then the optimization is to minimize the sum of all working_drivers. If the goal is not to minimize drivers, then it must be to minimize the total delay between tasks, which is a good measure of how efficient the assignment is. The assignment minimizes drivers waiting for their next shift as much as possible, allowing for the required breaks and the minimum pause between shifts.

Does it work?

The code works okay on the tiny (1.8s) and small (28s) instances, but not well at all on the medium (no result) and large (no result) instances! This surprised me. As a rule of thumb, Google’s OR-Tools team is good about including high quality examples. Perhaps my laptop isn’t fast enough, but I’m also wondering if maybe this is just older code. Looking at the commit logs, I see that the major work on these programs was some 5 years ago or more. I know that the circuit constraint got much better about a year ago, so maybe they just haven’t gotten around to updating this example.

I did some work on integrating the circuit constraint call, but I’ve been sitting on this post long enough. That work will be the subject of my next post, along with the false starts and steadily worsening results I was getting because I was thinking about it wrong. But that is for next time.