Speeding up the bus shift scheduling example

James Marca

Overview

• Circuit constraint example
• Bus shift scheduling
• Circuit constraint motivation
• Multiple circuit constraint for VRP
• Faster version of bus shift sheduling

Circuit constraint example

• simple example program from OR Tools
• Make a Hamiltonian graph
• Dropping nodes
• Why do it?

Simple example program

def main():
    """Entry point of the program."""
    num_nodes = len(DISTANCE_MATRIX)
    all_nodes = range(num_nodes)
    print("Num nodes =", num_nodes)

    # Model.
    model: cp_model.CpModel = cp_model.CpModel()

    obj_vars: list[cp_model.LiteralT] = []
    obj_coeffs: list[int] = []

    # Create the circuit constraint.
    arcs: list[cp_model.ArcT] = []
    arc_literals: dict[tuple[int, int], cp_model.LiteralT] = {}
    for i in all_nodes:
        for j in all_nodes:
            if i == j:
                continue

            lit = model.new_bool_var("%i follows %i" % (j, i))
            arcs.append((i, j, lit))
            arc_literals[i, j] = lit

            obj_vars.append(lit)
            obj_coeffs.append(DISTANCE_MATRIX[i][j])

    model.add_circuit(arcs)

    # Minimize weighted sum of arcs. Because this s
    model.minimize(sum(obj_vars[i] * obj_coeffs[i] for i in range(len(obj_vars))))

Example graph

Maybe a TSP tour

Make a Hamiltonian cycle

• The circuit constraint makes a Hamiltonian cycle.
• A cycle in a graph is a closed trail that only repeats the first and last vertices.
• A Hamiltonian cycle is a cycle that visits all the vertices of G.

Another code example

arcs: list[cp_model.ArcT] = []
for i in all_tasks:
    # if node i is first.
    start_lit = model.new_bool_var(f"start_{i}")
    arcs.append((0, i + 1, start_lit))

    # if node i is last.
    end_lit = model.new_bool_var(f"end_{i}")
    arcs.append((i + 1, 0, end_lit))

    for j in all_tasks:
        if i == j:
            continue
        lit = model.new_bool_var(f"arc_{i}_to_{j}")
        arcs.append((i + 1, j + 1, lit))

But what if we need to drop nodes

• cheat by adding self loops
• If a node is not visited in the Hamiltonian cycle, then it must loop on itself
• This makes the graph a Hamiltonian graph (as far as I can remember)

In code

arcs: list[cp_model.ArcT] = []
for i in all_tasks:
    # if node i is first.
    start_lit = model.new_bool_var(f"start_{i}")
    arcs.append((0, i + 1, start_lit))

    # if node i can be skipped, add a self loop
    skip_lit = model.new_bool_var(f"skip_{i}")
    arcs.append((i + 1, i + 1, skip_lit))
    
    # if node i is last.
    end_lit = model.new_bool_var(f"end_{i}")
    arcs.append((i + 1, 0, end_lit))

    for j in all_tasks:
        if i == j:
            continue
        lit = model.new_bool_var(f"arc_{i}_to_{j}")
        arcs.append((i + 1, j + 1, lit))

Why do it?

• The circuit constraint gives the solver knowledge of precedence
• The key is keeping track of the arc literals
• If a literal is true, then you know the node that came before
• If the literal from node 0 (source) is true, then you know the node is first
• If the literal to node 0 (sink) is true, then you know the node is last

The TSP example

for i in all_nodes:
    for j in all_nodes:
        if i == j:
            continue

        lit = model.new_bool_var("%i follows %i" % (j, i))
        arcs.append((i, j, lit))
        arc_literals[i, j] = lit

        obj_vars.append(lit)
        obj_coeffs.append(DISTANCE_MATRIX[i][j])

model.add_circuit(arcs)

# Minimize weighted sum of arcs.
model.minimize(sum(obj_vars[i] * obj_coeffs[i] for i in range(len(obj_vars))))
        

Result

• If an arc literal is true, the arc cost is in the objective
• If an arc is false, the arc cost is zeroed out
• Minimizing the cost therefore gets the shortest TSP solution
• This is not at all a fast way to solve the TSP
• The routing solver is better, and special purpose solvers are better still

Another example: Bus shift scheduling

• Assign bus shifts to drivers
• Minimize number of drivers
• Minimize excess break time
• Respect minimum break rules
• Assign a minimum amount of work
• Respect maximum driving and working limits

Two of these requirements are difficult to handle

• Minimum break rule: drivers cannot drive more than 4 hours without taking a 30 minute break
• Minimum work rule: all drivers must work at least 4 hours (start of first shift to end of last shift)
• Both of these require the solver to do some backtracking from good solutions

For each driver and each shift, we store:

• the total driving time including this shift
• the acrued driving time since the last 30 minute break

Special arcs have the following effect:

• ‘from source to shift’ sets the starting time and accumulate the first shift
• ‘from shift to end’ sets the ending time, and fill 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

The code

# Per (driver, node) info (driving time, performed,
# driving time since break)
total_driving = {}
no_break_driving = {}
performed = {}
starting_shifts = {}

# Per driver info (start, end, driving times, is working)
start_times = []
end_times = []
driving_times = []
working_drivers = []
working_times = []

# Weighted objective
delay_literals = []
delay_weights = []

# Used to propagate more between drivers
shared_incoming_literals = collections.defaultdict(list)
shared_outgoing_literals = collections.defaultdict(list)

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)
    )

    incoming_literals = collections.defaultdict(list)
    outgoing_literals = collections.defaultdict(list)
    outgoing_source_literals = []
    incoming_sink_literals = []

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

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

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

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

        # Node not performed
        #    - set both driving times to 0
        #    - add a looping arc on the node
        model.add(total_driving[d, s] == 0).only_enforce_if(~performed[d, s])
        model.add(no_break_driving[d, s] == 0).only_enforce_if(~performed[d, s])
        incoming_literals[s].append(~performed[d, s])
        outgoing_literals[s].append(~performed[d, s])
        # negated adding to the shared lists, because, globally, each node will
        # have one incoming literal, and one outgoing literal.

        # 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]
        )

        for o in range(num_shifts):
            other = shifts[o]
            delay = other[3] - shift[4]
            if delay < min_delay_between_shifts:
                continue
            lit = model.new_bool_var("%i from %i to %i" % (d, s, o))

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

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

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

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

    if minimize_drivers:
        # Driver is not working.
        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)
    else:
        # Working time constraints
        model.add(working_times[d] >= min_working_time)

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

# Each shift is covered.
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])

The add_circuit call is not used!

for d in range(num_drivers):
    ...
    # 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)

# Each shift is covered.
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])

Running the program

• tiny instance is fast
• small instance takes 40 to 50 seconds
• medium instance does not work
• large instance does not work

Does adding the proper circuit constraint speed this up?

• Nope

Does recognizing that this is really a multiple circuit speed this up?

• Nope

My solution

• A multiple circuit to determine arcs from break to break
• A hack to link up multiple break loops to a single vehicle
• Intervals to enforce no overlap constraints

Running my version

• tiny instance is fast
• small instance is fast, but slow to converge
• medium instance works, takes about 300s or more
• large instance does not work (laptop RAM issue)

The end for now

Circuit constraint motivation

Multiple circuit constraint for VRP

Faster version of bus shift sheduling