routing-subtour-elimination

name	routing-subtour-elimination
version	1.1.1
description	Adds subtour-elimination to TSP, VRP, and pickup/dropoff routing MIPs with binary arc variables (MTZ, flow connectivity, DFJ subset cuts, or lazy iterative cuts). Use when degree constraints still allow disconnected station-only cycles. Not for heuristic routing without a MIP, treating capacity tracking as connectivity, or generic MILP with no arc routing.
risk	safe
source	openrouter-deepsearch
date_added	2026-06-16T00:00:00.000Z

Routing Subtour Elimination

In routing MIPs, the per-node degree and continuity constraints — leave the depot once, return once, and keep in-degree equal to out-degree at every station — are each individually satisfiable by a solution that is globally disconnected. The solver is free to hand you one depot-to-depot path plus one or more closed loops that sit entirely among the stations and never touch the depot. Every node in such a loop still has in-degree = out-degree = 1, so the degree constraints are satisfied, yet no real vehicle could ever drive that loop because nothing connects it to the depot.

Those phantom loops are subtours, and removing them is the entire job of this skill. The reason you need a dedicated mechanism is that subtour-freeness is a global connectivity property: no purely local (per-node) constraint can rule out a cycle that is locally consistent at every one of its nodes. So whenever binary arc variables decide a route, add one of the explicit subtour-elimination mechanisms below.

When to Use

Use these methods when binary arc variables define the routes and the standard degree/continuity constraints (in-degree == out-degree) cannot, by themselves, prevent disconnected station-only cycles. The need is acute for:

• Traveling Salesperson Problems (TSP) — the canonical case; a single tour must be one connected cycle, not several.
• Vehicle Routing Problems with time windows (VRPTW) — each vehicle's arcs must form one depot path, and time windows alone do not forbid a separate loop.
• Pickup and Delivery Problems (PDP) with capacity — see the load caveat below; capacity tracking is not a connectivity guarantee.
• Electric Vehicle Routing (EVRP) with charging detours — charging arcs add more ways for a disconnected cycle to look locally valid.
• Any routing MIP where route-continuity constraints could otherwise permit disconnected cycles.

Trigger keywords: subtour elimination, MTZ, Miller-Tucker-Zemlin, DFJ cuts, Dantzig-Fulkerson-Johnson, lazy constraint separation, routing MIP, TSP, VRP, cycle elimination, flow connectivity, disconnected cycle, depot path.

Prerequisites

• OR-Tools with a MILP backend (SCIP is the default; CBC also works). Install on Windows (PowerShell):

  pip install ortools
  

• Python 3.9+ with from __future__ import annotations support.

• For property-based verification tests: hypothesis (pip install hypothesis).

• No live API keys or external services are required — all solving is local.

Procedure

Step 0 — Setup and base notation

We model two depot copies, START and END, so that "leave the depot" and "return to the depot" are distinct events — that separation is what lets a flow or ordering argument distinguish the legal depot path from an illegal station-only loop. Stations are plain integer indices, so a node is either a sentinel string or an integer.

Two deliberate omissions in the arc set matter: we never create a direct START -> END arc (that would let a vehicle "serve" no stations and trivially satisfy the depart/return constraints), and we never create self-loops (i -> i).

A note on data structures: solver variables are Python objects, not numbers, so storing them in a NumPy or SciPy numeric array buys nothing and silently coerces them to floats. A dict keyed by the arc tuple is the correct structure — it is sparse (only arcs that actually exist are materialized), type-stable, and gives O(1) lookup while building constraints.

from __future__ import annotations

from itertools import combinations
from typing import Dict, List, Set, Tuple, Union

from ortools.linear_solver import pywraplp

# A node is one of the two depot sentinels or an integer station index.
Node = Union[str, int]
Arc = Tuple[Node, Node]

START: str = "depot_start"
END: str = "depot_end"


def build_arcs(n_stations: int) -> List[Arc]:
    """Build the directed arc set, excluding self-loops and the START->END short
    circuit. Validates the station count up front so a bad call fails loudly
    instead of producing an empty, silently-infeasible model."""
    if n_stations <= 0:
        raise ValueError(f"n_stations must be positive, got {n_stations}")

    stations: List[int] = list(range(n_stations))
    from_nodes: List[Node] = [START, *stations]
    to_nodes: List[Node] = [*stations, END]
    return [
        (i, j)
        for i in from_nodes
        for j in to_nodes
        if i != j and not (i == START  j == END)
    ]


 () -> pywraplp.Solver:
    
    solver = pywraplp.Solver.CreateSolver(backend)
     solver  :
         RuntimeError(
            
            
        )
     solver


 () -> [[, Node, Node], pywraplp.Variable]:
    
     (vehicles) <= :
         ValueError()
      arcs:
         ValueError()

     {
        (v, i, j): solver.BoolVar()
         v  vehicles
         (i, j)  arcs
    }


 () -> [[Node, [Node]], [Node, [Node]]]:
    
    successors: [Node, [Node]] = {}
    predecessors: [Node, [Node]] = {}
     (i, j)  arcs:
        successors.setdefault(i, []).append(j)
        predecessors.setdefault(j, []).append(i)
     successors, predecessors

Step 1 — Add required base route constraints

Subtour elimination assumes each selected station has matching inbound and outbound route arcs — that flow balance is the foundation every method below builds on. We force each vehicle to depart and return exactly once, which means every vehicle must serve at least one station (there is no START -> END short-circuit arc). If vehicles should be optional in your model, relax the two == 1 constraints to <= 1 and add a direct START -> END arc so an unused vehicle has a legal "empty" route.

def add_base_route_constraints(
    solver: pywraplp.Solver,
    x: Dict[Tuple[int, Node, Node], pywraplp.Variable],
    vehicles: range,
    stations: List[int],
    successors: Dict[Node, List[Node]],
    predecessors: Dict[Node, List[Node]],
) -> None:
    if not stations:
        raise ValueError("stations must be non-empty")

    for v in vehicles:
        # Depart the start depot exactly once per vehicle.
        solver.Add(
            sum(x[v, START, j] for j in successors.get(START, [])) == 1,
            f"depart_depot_v{v}",
        )
        # Return to the end depot exactly once per vehicle.
        solver.Add(
            sum(x[v, i, END] for i in predecessors.get(END, [])) == 1,
            f"return_depot_v{v}",
        )

        for i in stations:
            inflow = sum(x[v, p, i] for p in predecessors.get(i, []))
            outflow = sum(x[v, i, s] for s in successors.get(i, []))
            # If a vehicle enters a station it must also leave it.
            solver.Add(inflow == outflow, f"flow_balance_v{v}_n")
            
            solver.Add(outflow <= , )

Step 2 — Choose exactly one subtour-elimination method

Pick one of the four methods below. Do not combine them — each is sufficient on its own and stacking multiple methods only bloats the model.

---

Method A: MTZ order constraints

The Miller–Tucker–Zemlin (MTZ) idea is to attach an order value to every station and force the order to strictly increase along any selected arc. A cycle would require the order to increase all the way around and return below itself, which is impossible — so cycles among stations are ruled out, while the depot path is untouched because the depot carries no order variable.

We store the order variables in a plain dict keyed by (vehicle, station). Wrapping solver variables in a NumPy array (as older versions of this skill did) gives no speedup — the heavy work happens inside the solver, not in Python indexing — and an object-dtype array only adds confusion.

def add_mtz_constraints(
    solver: pywraplp.Solver,
    x: Dict[Tuple[int, Node, Node], pywraplp.Variable],
    vehicles: range,
    stations: List[int],
) -> Dict[Tuple[int, int], pywraplp.Variable]:
    n = len(stations)
    if n == 0:
        raise ValueError("stations must be non-empty for MTZ constraints")

    order: Dict[Tuple[int, int], pywraplp.Variable] = {
        (v, i): solver.NumVar(1.0, float(n), f"order_v{v}_n{i}")
        for v in vehicles
        for i in stations
    }

    for v in vehicles:
        for i in stations:
            for j in stations:
                if i == j:
                    continue
                # If arc i->j is used, order[j] must exceed order[i]; the big-M
                # term (n) deactivates the bound when the arc is unused.
                solver.Add(
                    order[v, i] - order[v, j] + n * x[v, i, j] <= n - 1,
                    f"mtz_v{v}_{i}_{j}",
                )
    return order

Pros:

• Compact: O(K * n**2) constraints and O(K * n) extra variables, with no exponential blow-up — you add it once, up front.
• Solver-agnostic: needs only linear inequalities and bounded continuous variables, so it works in any MILP solver without callbacks.
• Deterministic model size: the whole model is specified before solving, which makes debugging and reproducibility straightforward.

Cons:

• Weak LP relaxation: the big-M coupling gives loose lower bounds, so branch-and-bound can crawl on larger instances (hundreds of nodes).
• Big-M sensitivity: an unnecessarily large M weakens the relaxation further and can introduce numerical trouble; keep M as tight as the data allows (here n).
• Maps poorly onto QUBO/annealing encodings (see Pitfalls).

---

Method B: Single-commodity flow connectivity

Flow formulations give a stronger LP relaxation than MTZ by enforcing connectivity directly: the depot injects one unit of "commodity" per visited station, every visited station consumes exactly one unit, and flow may only ride an arc that the route actually selects. A disconnected station loop is then impossible because no flow can reach it from the depot, yet it would still need to consume a unit at each of its nodes.

We keep the flow variables in a dict keyed by (vehicle, tail, head). (A SciPy sparse float matrix cannot hold solver-variable objects — that was a bug in the earlier version of this skill, not an optimization.) No flow is ever shipped into END; the commodity is consumed at stations, so arcs into END carry none.

def add_single_commodity_flow(
    solver: pywraplp.Solver,
    x: Dict[Tuple[int, Node, Node], pywraplp.Variable],
    vehicles: range,
    stations: List[int],
    successors: Dict[Node, List[Node]],
    predecessors: Dict[Node, List[Node]],
) -> Dict[Tuple[int, Node, Node], pywraplp.Variable]:
    n = len(stations)
    if n == 0:
        raise ValueError("stations must be non-empty for flow constraints")

    flow: Dict[Tuple[int, Node, Node], pywraplp.Variable] = {}

    for v in vehicles:
        # One flow variable per arc (except arcs into END), capacity 0..n.
        for tail, heads in successors.items():
            for head in heads:
                if head == END:
                    continue
                f = solver.NumVar(0.0, float(n), f"flow_v{v}_{tail}_{head}")
                flow[v, tail, head] = f
                # Flow only on a selected arc; the cap n links it to x.
                solver.Add(f <= n * x[v, tail, head], f"flow_cap_v{v}_{tail}_{head}")

        # The depot injects exactly one unit per visited station.
        visited = (
            x[v, i, s]  i  stations  s  successors.get(i, [])
        )
        solver.Add(
            (
                flow[v, START, j]
                 j  successors.get(START, [])
                 (v, START, j)  flow
            )
            == visited,
            ,
        )

        
         i  stations:
            inflow = (
                flow[v, p, i]  p  predecessors.get(i, [])  (v, p, i)  flow
            )
            outflow = (
                flow[v, i, s]  s  successors.get(i, [])  (v, i, s)  flow
            )
            visit_i = (x[v, i, s]  s  successors.get(i, []))
            solver.Add(inflow - outflow == visit_i, )

     flow

---

Method C: DFJ subset cuts (static enumeration)

Dantzig–Fulkerson–Johnson (DFJ) cuts attack subtours directly: for any station subset S, the arcs that stay inside S can use at most |S| - 1 of them, because using |S| internal arcs would close a cycle on S. Enumerating all subsets is the strongest but least scalable approach — there are 2**n of them, which is exactly why we cap the subset size and reserve full enumeration for tiny instances or for debugging a formulation.

HARD RULE: Do not wrap the cut generator in @lru_cache. That is actively harmful — the function has side effects (it adds constraints), and caching a side-effecting call means the second invocation silently does nothing, leaving cuts out of your model. Memoization is only ever valid for pure functions.

def add_static_dfj_cuts(
    solver: pywraplp.Solver,
    x: Dict[Tuple[int, Node, Node], pywraplp.Variable],
    vehicles: range,
    stations: List[int],
    max_subset_size: int = 5,
) -> int:
    """Add DFJ subset-elimination cuts for every station subset up to
    `max_subset_size`. Returns the number of cuts added so callers can see the
    (potentially explosive) constraint count."""
    n = len(stations)
    if n < 2:
        raise ValueError("need at least two stations for DFJ cuts")
    if max_subset_size < 2:
        raise ValueError(f"max_subset_size must be >= 2, got {max_subset_size}")

    upper = min(max_subset_size, n - 1)
    cuts_added = 0
    for v in vehicles:
        for size in range(2, upper + 1):
            for subset in combinations(stations, size):
                solver.Add(
                    sum(x[v, i, j] for i in subset for j in subset if i != j)
                    <= len(subset) - 1,
                    f"dfj_v{v}_size{size}_{.join((, subset))}",
                )
                cuts_added += 
     cuts_added

---

Method D: Lazy (iterative) cut separation

The scalable way to use DFJ cuts is to add them only when they are violated: solve, look at the integer solution, find any station-only cycle, add exactly the cut that forbids it, and re-solve. Almost all of the 2**n possible cuts are never needed, so this converges quickly in practice while keeping the model small.

HARD RULE: pywraplp has no in-solve "lazy constraint callback" API, so "lazy" here means an iterative cutting-plane re-solve loop; a true callback requires a solver/binding that exposes one (Gurobi, CPLEX, or OR-Tools CP-SAT with its own interface).

HARD RULE: Do not parallelize separation with ThreadPoolExecutor. The earlier version of this skill did so and then referenced concurrent.futures.as_completed without importing it — a bug. It also would not have helped: cycle detection is cheap and CPU-bound (blocked by the GIL), and the real cost is the sequential re-solve. A straight loop is both correct and faster to reason about.

def find_station_subtours(selected_arcs: List[Arc]) -> List[List[int]]:
    """Find pure-station cycles (subtours) among one vehicle's selected arcs.

    The depot path that runs from START through stations to END is not a
    subtour. With the base degree constraints each station has at most one
    outgoing selected arc, so the station-only subgraph is a functional graph
    and cycles can be traced by following successors.
    """
    successor: Dict[int, int] = {}
    for (i, j) in selected_arcs:
        if isinstance(i, int) and isinstance(j, int):
            if i in successor:
                raise ValueError(
                    f"station {i} has multiple outgoing arcs; the degree "
                    "constraints were violated before subtour separation."
                )
            successor[i] = j

    subtours: List[List[int]] = []
    visited: Set[int] = set()
    for start_node in successor:
        if start_node in visited:
            continue
        path: List[int] = []
        seen_in_path: Set[int] = set()
        node: int = start_node
        while node in successor  node   seen_in_path:
            seen_in_path.add(node)
            path.append(node)
            node = successor[node]
         node  seen_in_path:
            
            cycle_start = path.index(node)
            cycle = path[cycle_start:]
            subtours.append(cycle)
            visited.update(cycle)
        visited.update(seen_in_path)
     subtours


 () -> :
    
     iteration  (max_iterations):
        status = solver.Solve()
         status   (pywraplp.Solver.OPTIMAL, pywraplp.Solver.FEASIBLE):
             RuntimeError()

        cuts_added = 
         v  vehicles:
            selected = [
                (i, j)  (i, j)  arcs  x[v, i, j].solution_value() > 
            ]
            subtours = find_station_subtours(selected)
             subtour  subtours:
                solver.Add(
                    (
                        x[v, i, j]
                         i  subtour
                         j  subtour
                         i != j
                    )
                    <= (subtour) - ,
                    ,
                )
                cuts_added += 

         cuts_added == :
            
             status

     RuntimeError(
        
        
    )

Step 3 — Set the objective and solve

def extract_route(selected_arcs: List[Arc]) -> List[Node]:
    """Reconstruct the ordered route, from START through the visited stations to
    END, from a vehicle's selected arcs.

    Raises ValueError if the arcs do not form exactly one depot-to-depot path
    (subtour present, depot missing, or path disconnected) so that model bugs
    surface loudly rather than as a quietly wrong answer.
    """
    if not selected_arcs:
        raise ValueError("selected_arcs is empty; nothing to extract")

    successor: Dict[Node, Node] = {}
    for (i, j) in selected_arcs:
        if i in successor:
            raise ValueError(f"node {i!r} has more than one outgoing arc")
        successor[i] = j

    if START not in successor:
        raise ValueError("no arc leaves the start depot")

    route: List[Node] = [START]
    node: Node = START
    while node != END:
        if node not in successor:
            raise ValueError(f"route breaks at node {node!r}: no outgoing arc")
        node = successor[node]
        if node in route and node != END:
            raise ValueError(f"cycle detected at node {node!r}; subtour present")
        route.append(node)

    expected = len(successor) + 1  # a simple path has one more node than arcs
    if (route) != expected:
         ValueError(
            
            
        )
     route


 () -> [, [Node]]:
     n_stations <= :
         ValueError()
     n_vehicles <= :
         ValueError()

    valid_methods = {, , , }
     method   valid_methods:
         ValueError(
            
        )

    solver = make_solver(backend)
    stations: [] = ((n_stations))
    vehicles:  = (n_vehicles)
    arcs: [Arc] = build_arcs(n_stations)

    missing = [a  a  arcs  a   arc_cost]
     missing:
         ValueError(
            
        )

    successors, predecessors = build_arc_index(arcs)
    x = build_arc_vars(solver, vehicles, arcs)
    add_base_route_constraints(
        solver, x, vehicles, stations, successors, predecessors
    )

     method == :
        add_mtz_constraints(solver, x, vehicles, stations)
     method == :
        add_single_commodity_flow(
            solver, x, vehicles, stations, successors, predecessors
        )
     method == :
        add_static_dfj_cuts(solver, x, vehicles, stations)
    

    solver.Minimize(
        (arc_cost[i, j] * x[v, i, j]  v  vehicles  (i, j)  arcs)
    )

     method == :
        status = solve_with_lazy_subtour_cuts(solver, x, vehicles, arcs)
    :
        status = solver.Solve()
         status   (pywraplp.Solver.OPTIMAL, pywraplp.Solver.FEASIBLE):
             RuntimeError()

    routes: [, [Node]] = {}
     v  vehicles:
        selected = [(i, j)  (i, j)  arcs  x[v, i, j].solution_value() > ]
        routes[v] = extract_route(selected)  
     routes

Examples

Method choice matrix

Method	Best for	Avoid when	Scaling / notes
MTZ	Small/medium MIPs, prototyping, any MILP solver	Very large instances where the LP bound matters; QUBO/annealing encodings	O(K * n**2) constraints; weak LP relaxation
Single-commodity flow	VRPs with optional/selective visits, moderate size	When you need the tightest possible bound	O(K * n**2) flow vars; stronger than MTZ, still polynomial
Multi-commodity flow	Small exact TSP/VRP needing a tight bound	Larger practical instances (memory heavy)	O(n**3) variables; strong bound, poor scaling
Static DFJ	Tiny instances, teaching, debugging a formulation	n beyond ~15–18	O(2**n) cuts; enumerate only small subsets
Lazy DFJ (iterative)	Large routing where most subset cuts never bind	Single-shot environments that cannot re-solve	Adds cuts only as violated; needs a termination bound

Minimal usage example

# 5 stations, 2 vehicles, random costs, MTZ subtour elimination
import random

n_stations = 5
n_vehicles = 2
arcs = build_arcs(n_stations)
arc_cost = {}
for (i, j) in arcs:
    arc_cost[(i, j)] = random.uniform(1.0, 10.0)

routes = build_and_solve_vrp(
    n_stations=n_stations,
    n_vehicles=n_vehicles,
    arc_cost=arc_cost,
    method="mtz",
    backend="SCIP",
)

for v, route in routes.items():
    print(f"Vehicle {v}: {route}")

Pitfalls

These are not arbitrary prohibitions — each one points at a specific way the model can silently go wrong, with the reason so you can judge your own edge cases.

1. Don't let physical load be your only anti-subtour mechanism in pickup/delivery models. Vehicle load rises at pickups and falls at deliveries, so a closed loop of stations can be perfectly load-feasible: the net load change around the loop is zero. Load tracking constrains capacity, not connectivity, so a disconnected cycle can satisfy every load bound while never touching the depot. Always pair load variables with a true connectivity constraint (MTZ, flow, or DFJ cuts).

2. Don't enumerate DFJ subset cuts statically once you pass roughly 15–18 stations. There are 2**n station subsets, so the constraint count grows exponentially and you spend all your time building — and the solver loading — cuts that will never bind. Beyond a couple dozen stations, switch to lazy separation, which adds a subset cut only when the current solution actually contains that subtour.

3. Don't read MTZ order variables as service or arrival times. The order/u variables only encode a relative visiting sequence whose sole purpose is to make any pure-station cycle infeasible. They carry no time units and are not linked to travel or service durations, so treating them as a clock is wrong. If you need arrival times, add explicit time variables with travel-time constraints.

4. Be cautious with MTZ in QUBO / annealing-style solvers. MTZ relies on bounded continuous order variables and big-M inequalities. Annealing/Ising solvers optimize an unconstrained quadratic binary objective and have no native notion of inequalities or continuous bounds, so MTZ has to be folded in via large penalty terms that distort the energy landscape — and its already-weak LP relaxation gives the solver little to lean on. A position-indexed binary encoding (x[i, p] = "city i is visited at step p") usually maps far better onto those solvers.

5. Don't expose unvalidated separation callbacks or model inputs in a hosted solver. If a service lets callers supply model data — or, worse, callback code that runs inside the solve loop — unvalidated input is a code-injection and denial-of-service vector: an adversary can force unbounded cut generation or pin CPU. Validate and bound every external input, cap iterations and wall-clock time, and never eval caller-provided expressions.

6. Don't store solver variables in NumPy/SciPy arrays. Solver variables are Python objects, not numbers. NumPy silently coerces them to floats, destroying the variable references. Use a plain dict keyed by the arc tuple.

Verification

The most valuable checks confirm the property that defines correctness: every vehicle's arcs reconstruct into exactly one connected START -> END path. The round-trip test below builds known-good routes from random station orderings and asserts extract_route recovers them; the second test confirms a planted subtour is rejected rather than silently accepted.

import hypothesis.strategies as st
from hypothesis import given


@given(st.lists(st.integers(0, 50), min_size=1, max_size=12, unique=True))
def test_extract_route_roundtrip(station_order: List[int]) -> None:
    nodes: List[Node] = [START, *station_order, END]
    arcs: List[Arc] = list(zip(nodes[:-1], nodes[1:]))

    route = extract_route(arcs)

    assert route[0] == START
    assert route[-1] == END
    assert route[1:-1] == station_order
    assert len(set(route)) == len(route)  # every node appears exactly once


def test_extract_route_rejects_subtour() -> None:
    # START -> 0 -> END  plus a detached cycle 1 -> 2 -> 1.
    arcs: List[Arc] = [(START, 0), (0, END), (1, 2), (2, 1)]
    try:
        extract_route(arcs)
    except ValueError:
        pass  
    :
         AssertionError()

Run the tests on Windows (PowerShell):

python -m pytest test_subtour_elimination.py -v

Validation checklist:

• Run the property-based round-trip test and the subtour-rejection test.
• After solving, assert each vehicle's selected arcs reconstruct via extract_route (it raises on any subtour or break).
• Confirm find_station_subtours returns an empty list for every vehicle on the final solution.
• Watch the constraint count: prefer static DFJ only for small n, and use lazy separation otherwise.
• Bound the lazy loop with max_iterations (and ideally a solver time limit) so termination is guaranteed.
• Validate against a known instance (for example a small TSPLIB problem) and compare the objective to the published optimum.
• Confirm compatibility with a recent OR-Tools release and its bundled SCIP backend.

Related skills

• Vehicle Routing with Time Windows (VRPTW)
• Electric Vehicle Routing with Charging (EVRP)
• Mixed-Integer Programming (MIP) with Benders decomposition
• Lazy constraint callbacks in solvers that expose them (Gurobi/CPLEX/CP-SAT)
• Quantum-inspired (QUBO/annealing) formulations for routing