Constraint programming for routing

Read first

• Post 3: The TSP family. The time-window version of our eight stops: B must be reached by minute 60 and F by minute 140, one distance unit takes one minute. Of the 5,040 directed orders, 144 are feasible, and the best is A-B-G-F-E-H-D-C-A at 359.
• Post 5: Branch and bound. Search a tree of partial solutions depth-first, keep the best complete one as an upper bound, and prune any branch whose lower bound can’t beat it.

The question this post answers: can we write a routing problem down as rules, and let the rules themselves throw away impossible routes before we ever search them?

Intuition

Most real routing requests are not “find the shortest tour”. They are a list of rules: this customer only before 10:00, that one needs the refrigerated van, this driver must be back by 17:00. Constraint programming (CP) is built for exactly that. You state the rules, and the solver uses them to rule things out as it goes.

We’ll use the time-window instance from Post 3. Here are the eight stops again, with the depot at A.

Distances are in units of 100 m, rounded, and driving one unit takes one minute. Darker cells are longer legs.

A CP model has three parts:

• Variables: the decisions. Here, “which stop comes next” at each step of the route.
• Domains: the values each variable may still take. At the start, the next stop can be any of the seven unvisited stops.
• Constraints: the rules. Visit every stop exactly once, and reach B and F before their deadlines.

The key move is propagation: using a constraint to delete values from a domain before searching them. The reasoning we need is one line long. Straight-line distances obey the triangle inequality, and our rounded matrix still does (I checked every triple of stops), so from wherever the van is now, no route reaches B sooner than driving there directly. If even the direct drive is too late, every route through this point is dead.

Level 1: which stop comes first? For each candidate first stop, the table below shows when the van gets there and the earliest it could then reach each deadline stop.

Six of the seven values are gone, so the first stop is forced: it has to be B. No route was built to learn this.

Level 2: after A-B (minute 41), which stop comes second?

Level 3, on the G branch. After A-B-G the clock reads 86. The cheapest detour before F is through E, which arrives at F at 86 + 41 + 41 = 168, too late. So F is forced next, at minute 128.

That accounts for every feasible order from Post 3. Either the route starts A-B-F and the other five stops go in any order (5! = 120 orders), or it starts A-B-G-F and the other four go in any order (4! = 24). That is 144, and propagation found the structure after looking at three levels.

The diagram below shows the search tree that propagation leaves behind. Each pruned box stands for a whole subtree that is never visited.

flowchart TD
  R["A, minute 0"] --> B["A-B, minute 41"]
  R --> X1["6 other first stops: pruned, B late"]
  B --> G["A-B-G, minute 86"]
  B --> F["A-B-F, minute 127"]
  B --> X2["C, D, E, H second: pruned, F late"]
  G --> GF["A-B-G-F, minute 128"]
  G --> X3["anything else third: pruned, F late"]
  F --> L1["any order of C, D, E, G, H: 120 routes"]
  GF --> L2["any order of C, D, E, H: 24 routes"]

Feasible is not optimal, so the search still has to choose among the 144 routes. It does that the way Post 5 did, with a bound: the cost so far plus the direct drive home is a lower bound on any completion, and a branch whose bound already reaches the best route found is cut. (Here time and distance are the same number, because the windows only have deadlines and nobody waits.)

Step through the diagram below to follow the search. The lengths it shows are unrounded kilometres, a few tenths off the matrix units.

The optimal route, priced leg by leg:

B is reached at minute 41 and F at minute 128, both on time.

The formal idea

Constraint satisfaction. A constraint satisfaction problem has variables x₁, …, xₖ, a domain Dᵢ of allowed values for each xᵢ, and a set of constraints, each of which says which combinations of values are allowed. A solution gives every variable a value from its domain that satisfies every constraint. Add an objective to minimize and it becomes a constraint optimization problem.

A routing model. One standard model gives every stop i a successor variable next(i), the stop visited right after i, and a time variable Tᵢ, the arrival time at i. With travel time t(i, j) and deadline lᵢ:

next(i) ∈ {all stops} ∖ {i}                    for every stop i
all next(i) different, forming one cycle       (one tour, no subtours)
T(next(i)) ≥ Tᵢ + t(i, next(i))                for every stop i except the last before the depot
Tᵢ ≤ lᵢ                                        for every stop with a deadline
minimize  Σᵢ d(i, next(i))

My Python version builds the route from the depot forwards, which fixes next(A), then next of that stop, and so on. It is the same model, searched in route order.

Propagation. A propagator for a constraint removes values from domains that cannot appear in any solution, given the domains of the other variables. Propagation runs every propagator until nothing changes. The rule I used above is a bound propagator on the time variables: by the triangle inequality, Tₛ ≥ T_here + t(here, s) for any stop s not yet visited, so if T_here + t(here, s) > lₛ, the current partial route has no feasible completion. Checking that rule on the stops still to come, right after each choice, is called forward checking (Haralick and Elliott, 1980). Stronger propagators exist. The best known is for “all different”: Régin (1994) showed how to remove every value that cannot be part of any all-different assignment, using bipartite matching.

Search. Propagation alone rarely decides everything. CP alternates: propagate, pick a variable, try a value, propagate again, and back up when a domain empties. For optimization it adds the constraint “cost < best found so far” after each improvement, which is branch and bound from Post 5 expressed as one more constraint.

The diagram below shows that loop.

flowchart TD
  P["Propagate all constraints"] --> E{"Any domain empty?"}
  E -- yes --> BK["Back up to the last choice"]
  E -- no --> D{"Every variable fixed?"}
  D -- yes --> S["New best route: require a lower cost from now on"]
  D -- no --> C["Choose a variable and a value"]
  C --> P
  S --> BK
  BK --> P

How OR-Tools uses CP for routing. Google’s OR-Tools ships a routing library built on its original CP solver. The model is close to the one above: one “next” variable per location (NextVar), and dimensions for quantities that accumulate along a route, such as time or load, each with a cumulative variable at every location (CumulVar) whose bounds the solver propagates. Time windows and capacities are constraints on those cumulative variables.

Where it differs from my toy is the search. OR-Tools does not try to prove optimality on real-size problems. It builds a first solution with a fast rule, such as PATH_CHEAPEST_ARC, which extends each route by its cheapest arc, then improves it with local search moves, optionally guided by a metaheuristic such as GUIDED_LOCAL_SEARCH, until a time limit. CP’s job inside that loop is feasibility: every candidate move is checked by propagation before it is accepted.

flowchart LR
  M["Model: next variables, dimensions, windows"] --> F["First solution: fast construction rule"]
  F --> L["Local search move"]
  L --> P{"Propagation: still feasible?"}
  P -- no --> L
  P -- yes --> A{"Better, or accepted by metaheuristic?"}
  A --> L
  A -. time limit .-> R["Return best route found"]

Shaw (1998) was one of the first to combine the two: his large neighbourhood search removes a group of related visits from the routes and reinserts them with a CP tree search. We’ll meet it again in Part 6.

Common misconceptions

“Constraint programming doesn’t search; it just solves the constraints.” It searches, depth-first, with backtracking. Propagation shrinks the tree, it rarely removes it. On our instance, propagation reduced 13,700 search nodes to 408, and the solver still needed a cost bound and 197 nodes to prove 359 optimal.

“More constraints make the problem harder.” For CP, tight constraints usually make it easier, because each one gives propagation something to cut with. The two deadlines here took the search from 13,700 nodes to 408. Without any time windows, the same cost-bound search on the plain eight-stop TSP needs 3,127 nodes to prove 315 optimal, against 197 with the windows. (Deciding whether any feasible route exists can still be hard in general: Savelsbergh (1985) showed it is NP-complete for the TSPTW.)

“Propagation finds the best route.” Propagation only removes what is impossible. All 144 routes it left were feasible, and only one was optimal. Finding it took search plus a bound.

“OR-Tools’ routing solver returns the optimal route.” Its own documentation says it “sometimes returns solutions that are good, but not optimal.” For realistic sizes it returns the best route its local search found within the time limit, with no proof of optimality. CP inside it guarantees feasibility, not optimality.

Try it in Python

This script is the whole model: the variable is the next stop, its domain is the unvisited stops, and the four runs switch on one more piece of reasoning each.

import math

STOPS = {"A": (2, 2), "B": (6, 1), "C": (9, 3), "D": (8, 7),
         "E": (5, 9), "F": (1, 8), "G": (4, 5), "H": (10, 9)}  # km; A is the depot
LATE = {"B": 60, "F": 140}  # latest arrival in minutes; one distance unit takes one minute

def d(a, b):  # straight-line distance in units of 100 m, rounded
    (x1, y1), (x2, y2) = STOPS[a], STOPS[b]
    return int(10 * math.hypot(x1 - x2, y1 - y2) + 0.5)

def on_time(route):  # does a finished route meet every deadline?
    time = 0
    for a, b in zip(route, route[1:]):
        time += d(a, b)
        if b in LATE and time > LATE[b]:
            return False
    return True

def solve(check_arrival, forward_check, cost_bound):
    """Build the route one stop at a time. At each level the variable is
    'which stop comes next' and its domain is the set of unvisited stops."""
    nodes, complete, best = 0, 0, (math.inf, None)

    def search(route, time):
        nonlocal nodes, complete, best
        nodes += 1
        here = route[-1]
        unvisited = [s for s in STOPS if s not in route]
        if check_arrival and here in LATE and time > LATE[here]:
            return                     # this stop's deadline is already missed
        if forward_check and any(time + d(here, s) > LATE[s] for s in LATE if s in unvisited):
            return                     # a stop still to come can no longer make its deadline
        if cost_bound and time + d(here, "A") >= best[0]:
            return                     # even driving straight home can't beat the best route
        if not unvisited:
            complete += 1
            total = time + d(here, "A")
            if on_time(route) and total < best[0]:
                best = (total, "-".join(route + ["A"]))
            return
        for s in sorted(unvisited, key=lambda s: d(here, s)):  # nearest stop first
            search(route + [s], time + d(here, s))

    search(["A"], 0)
    return nodes, complete, best

for name, flags in [("Generate and test", (False, False, False)),
                    ("Check each arrival", (True, False, False)),
                    ("+ forward checking", (True, True, False)),
                    ("+ cost bound", (True, True, True))]:
    nodes, complete, (cost, route) = solve(*flags)
    print(f"{name:<19} nodes {nodes:>6}   complete routes {complete:>5}   best {cost}  {route}")

It prints:

Each row is the same search with one more piece of reasoning. “Generate and test” builds all 5,040 orders and checks the deadlines at the end; its 13,700 nodes are every partial route of every length, plus the root. Checking each arrival as it happens stops a branch the moment a deadline is missed. Forward checking stops it earlier, the moment a deadline will be missed, which is the reasoning from the tables above. The cost bound turns the search into branch and bound: it builds only two complete routes, 365 and then 359, and spends the rest of its 197 nodes proving that nothing beats 359.

Exercise

Tighten F’s deadline from minute 140 to minute 127 and keep B’s at 60. Without a computer, use the two tables from the Intuition section to decide which stops can now come second, and how many feasible orders remain. Then tighten F to minute 126. What does propagation tell you, and how many nodes do you expect the search to need before it can say so?

Hint: the level-2 table already lists the earliest possible arrival at F for every second stop. Then run the script with LATE = {"B": 60, "F": 127} and check your count against “complete routes” in the “Check each arrival” row.

What this method can’t do yet

Constraint programming did well here because the deadlines were tight and the instance was tiny. Take them away and the picture changes. With no time windows there is nothing for propagation to cut, and the only pruning left is the cost bound. I ran that cost-bound search, nearest stop first, on the first n cities of berlin52:

Thirteen cities already took about 9 seconds in plain Python on my laptop, and each extra city multiplies the work several times over. Real CP solvers propagate far more cleverly than my one-line rule, but the shape is the same: proving optimality grows exponentially, and a dispatcher planning hundreds of stops every morning can’t wait for a proof. Daily operations need a good route in seconds, and the proof can go. That is why OR-Tools runs CP inside a local search rather than a full tree search, and it is the problem Post 9, “Construction heuristics”, starts to solve: build a reasonable route fast, with no search at all.

References

• Rossi, F., van Beek, P., and Walsh, T., eds. (2006). Handbook of Constraint Programming. Elsevier.
• Mackworth, A. K. (1977). “Consistency in Networks of Relations.” Artificial Intelligence 8(1): 99–118.
• Haralick, R. M., and Elliott, G. L. (1980). “Increasing Tree Search Efficiency for Constraint Satisfaction Problems.” Artificial Intelligence 14(3): 263–313.
• Régin, J.-C. (1994). “A Filtering Algorithm for Constraints of Difference in CSPs.” Proceedings of AAAI-94, 362–367.
• Savelsbergh, M. W. P. (1985). “Local Search in Routing Problems with Time Windows.” Annals of Operations Research 4: 285–305.
• Pesant, G., Gendreau, M., Potvin, J.-Y., and Rousseau, J.-M. (1998). “An Exact Constraint Logic Programming Algorithm for the Traveling Salesman Problem with Time Windows.” Transportation Science 32(1): 12–29.
• Shaw, P. (1998). “Using Constraint Programming and Local Search Methods to Solve Vehicle Routing Problems.” In Principles and Practice of Constraint Programming (CP98), LNCS 1520, Springer, 417–431.
• Kilby, P., and Shaw, P. (2006). “Vehicle Routing.” Chapter 23 in Rossi, van Beek and Walsh (eds.), Handbook of Constraint Programming, Elsevier, 801–836.
• Google. OR-Tools documentation: “Vehicle Routing” and “Routing Options.” https://developers.google.com/optimization/routing

See the whole series · ← All writing