Presolve

Presolve#

Before any branch-and-bound or relaxation work, discopt runs a presolve orchestrator that walks a registered list of passes over the model to a fixed point. Each pass is a small, focused transformation that either tightens variable bounds, eliminates variables, drops constraints, rewrites constraint bodies, or detects structure. The orchestrator iterates passes under a global iteration / time / work budget until no pass reports progress.

The motivating numbers are large: classical MIP presolve already delivers a ~10.7× shifted geomean speedup on standard benchmarks [Achterberg et al., 2020], and for nonlinear and mixed-integer-nonlinear problems the literature documents 17–80% root-LP gap reductions and comparable wall-time savings [Belotti et al., 2009, Puranik and Sahinidis, 2017]. Every preprocessing dollar amortizes over the entire B&B tree, so leverage on MINLP is structurally larger than on MIP.

This notebook walks through the orchestrator’s API and demonstrates the new structural passes that ship in discopt-core’s presolve module.

Architecture#

The presolve module lives at crates/discopt-core/src/presolve/. Each pass implements the PresolvePass trait and emits a PresolveDelta describing what it changed. The orchestrator drives them to fixed point and records the chronological delta log. Python passes (e.g. JAX-backed reverse-mode AD interval tightening, NN-embedded MINLP presolve) plug into the same protocol via discopt._relax.presolve.

All shipped passes are deterministic by construction — no RNG, no unsorted hash iteration. The crate’s determinism harness (crates/discopt-core/tests/presolve_determinism.rs) verifies that byte-identical inputs produce byte-identical delta sequences.

Calling the orchestrator#

From Python, the orchestrator is reachable in two ways:

  1. Implicitly through Model.solve(presolve=True) (the default), which routes through discopt._relax.presolve_pipeline.run_root_presolve.

  2. Explicitly via PyModelRepr.presolve(passes=[...]) for users who want to inspect or customise the pass schedule.

Let’s build a small MINLP and look at what presolve does to it.

import discopt as do
from discopt._rust import model_to_repr

m = do.Model("presolve_demo")
x = m.continuous("x", lb=0.0, ub=10.0)
y = m.continuous("y", lb=0.0, ub=10.0)
b = m.binary("b")

# Two constraints with intentional structure for presolve to find:
#  - the equality fixes y in terms of x (factorable elimination, C2)
#  - the big-M row admits coefficient strengthening (C4)
m.subject_to(y - x * x == 0.0)
m.subject_to(5.0 * b + x <= 5.0)
m.minimize(x + b)

repr_ = model_to_repr(m)
print(
    f"before: {repr_.n_var_blocks} var blocks, "
    f"{repr_.n_constraints} constraints, "
    f"{repr_.n_nodes} arena nodes"
)

before: 3 var blocks, 2 constraints, 12 arena nodes

new_repr, stats = repr_.presolve()
print(f"after:  {new_repr.n_var_blocks} var blocks, {new_repr.n_constraints} constraints")
print(f"terminated_by: {stats['terminated_by']}")
print(f"iterations:    {stats['iterations']}")
for d in stats["deltas"]:
    progress = (
        int(d["bounds_tightened"]) > 0
        or len(d["vars_fixed"]) > 0
        or len(d["constraints_removed"]) > 0
        or len(d["constraints_rewritten"]) > 0
    )
    if progress:
        print(
            f"  {d['pass_name']:30s} "
            f"removed={list(d['constraints_removed'])} "
            f"rewritten={list(d['constraints_rewritten'])} "
            f"tightened={int(d['bounds_tightened'])}"
        )

after:  3 var blocks, 2 constraints
terminated_by: NoProgress
iterations:    2
  fbbt                           removed=[] rewritten=[] tightened=1

The default schedule runs eliminate → factorable_elim → simplify → coefficient_strengthening → fbbt → probing. Each pass ships as a first-class PresolvePass in crates/discopt-core/src/presolve/passes.rs and is dispatched by name from Python.

Coefficient strengthening (C4)#

Following Savelsbergh [1994], for a row \(\sum_j a_j x_j \le b\) with binary \(x_k\) and positive \(a_k\), define \(U_{-k} = \sum_{j\ne k} \max(a_j \ell_j, a_j u_j)\) — the largest LHS contribution from the non-\(k\) terms. If the slack \(b - U_{-k}\) is strictly between \(0\) and \(a_k\), we can simultaneously shift both the coefficient and the right-hand side by that slack: replace \(a_k \leftarrow a_k - \text{slack}\) and \(b \leftarrow b - \text{slack}\). The integer corners of the row are unchanged; the LP relaxation strictly tightens.

m = do.Model("c4")
b = m.binary("b")
y = m.continuous("y", lb=0.0, ub=2.0)
m.subject_to(5.0 * b + y <= 5.0)  # b=1 ⇒ y ≤ 0; b=0 ⇒ y ≤ 2 (slack 3)
m.minimize(b + y)
repr_ = model_to_repr(m)
_, stats = repr_.presolve(passes=["coefficient_strengthening"])
for d in stats["deltas"]:
    if d["pass_name"] == "coefficient_strengthening":
        print("row indices rewritten:", list(d["constraints_rewritten"]))

row indices rewritten: [0]
row indices rewritten: []

The original row 5b + y ≤ 5 becomes 2b + y ≤ 2 — same integer feasible set, strictly tighter LP relaxation.

Factorable-expression elimination (C2)#

When an equality uniquely determines a variable that appears nowhere else — even when the right-hand side is nonlinear in other variables — the equation can be dropped after a bounded-range safety check:

\[ v - g(x_1, \dots, x_n) = 0 \quad\Longrightarrow\quad v = g(x_1, \dots, x_n) \]

The pass uses interval forward-propagation (the same machinery FBBT uses) to verify that the derived range of \(v\) stays inside \(v\)’s declared box for every feasible assignment of the other variables; otherwise it abstains. This generalises the M10 singleton-equality elimination pass to nonlinear right-hand sides.

m = do.Model("c2")
v = m.continuous("v", lb=-100.0, ub=100.0)
x = m.continuous("x", lb=0.0, ub=2.0)
y = m.continuous("y", lb=0.0, ub=3.0)
m.subject_to(v - x * y == 0.0)  # v determined by x*y; v not in objective
m.minimize(x + y)
repr_ = model_to_repr(m)
new_repr, stats = repr_.presolve(passes=["factorable_elim"])
print(f"constraints before: {repr_.n_constraints}, after: {new_repr.n_constraints}")

constraints before: 1, after: 0

Symmetry detection (D4)#

Variable permutation symmetries inflate the search tree by a factor of the orbit size [Margot, 2010]. The detector compares each scalar variable’s column signature — its participation across constraints and the objective — and groups variables with identical fingerprints into candidate orbits. Detection is sound for fully-linear models; for nonlinear models the kernel emits orbits as candidates the caller can choose to break (or to leave alone, if structural verification is unavailable).

m = do.Model("symmetry")
x1 = m.binary("x1")
x2 = m.binary("x2")
x3 = m.binary("x3")
m.subject_to(x1 + x2 + x3 <= 2.0)
m.minimize(x1 + x2 + x3)
repr_ = model_to_repr(m)
result = repr_.detect_symmetries()
for orbit in result["orbits"]:
    print("orbit:", orbit["vars"])

orbit: [0, 1, 2]

Persistent in-tree FBBT (B3)#

Root presolve only catches what’s visible at the root. Inside the B&B tree, every branching decision shifts variable bounds, and FBBT can tighten other variables in response — but only if it runs at the node. The in-tree FBBT pass [Belotti et al., 2009] re-runs FBBT on the node’s local bounds, gated by a depth-stride schedule so the cost amortises over the tree.

import numpy as np

m = do.Model("in_tree_demo")
x = m.continuous("x", lb=0.0, ub=10.0)
y = m.continuous("y", lb=0.0, ub=10.0)
m.subject_to(x + y <= 5.0)
m.minimize(x + y)
repr_ = model_to_repr(m)

# Imagine the B&B branched on x ∈ [3, 10] at this node. FBBT propagates
# x + y ≤ 5 to derive y ≤ 2 — a tightening that persists down the subtree.
delta = repr_.in_tree_presolve(
    np.array([3.0, 0.0]),
    np.array([10.0, 10.0]),
    node_depth=0,
    depth_stride=1,
)
print("ran:        ", delta["ran"])
print("infeasible: ", delta["infeasible"])
print("tightened:  ", int(delta["bounds_tightened"]))
print("new ub:     ", list(delta["ub"]))

ran:         True
infeasible:  False
tightened:   2
new ub:      [np.float64(5.0), np.float64(2.0)]

Pass in_tree_presolve_stride > 0 to Model.solve(...) to enable this pass during the actual B&B search. The default is 0 (off) to keep behaviour conservative; values of 4–8 are typical, trading a small per-stride FBBT cost for tighter LP relaxations downstream.

Separability detection (D5)#

When the model splits into disjoint variable blocks linked by neither constraints nor objective, each block can be solved independently. Detection is a union-find over the constraint hypergraph; the result is reported as a list of variable blocks that the relaxation compiler can dispatch in parallel.

from discopt._relax.presolve import SeparabilityPass

m = do.Model("separable")
a = m.continuous("a", lb=0.0, ub=1.0)
b = m.continuous("b", lb=0.0, ub=1.0)
c = m.continuous("c", lb=0.0, ub=1.0)
d = m.continuous("d", lb=0.0, ub=1.0)
m.subject_to(a + b <= 1.5)
m.subject_to(c + d <= 1.5)
m.minimize(a + b + c + d)
repr_ = model_to_repr(m)
p = SeparabilityPass(m)
delta = p.run(repr_)
print("separable:", delta["separable"])
print("blocks:   ", [sorted(b) for b in delta["separable_blocks"]])

separable: False
blocks:    [['a', 'b', 'c', 'd']]

Periodic-variable bound reduction#

Angular variables are a recurring source of unbounded boxes in nonlinear models: a voltage angle in AC optimal power flow [Carpentier, 1962, Molzahn and Hiskens, 2019], or any variable that enters the model only as sin(theta) / cos(theta), is mathematically free, yet spatial branch-and-bound and the local sub-NLP cannot converge over an infinite domain.

The periodic-variable bound reduction pass (PeriodicVariableBoundRule in discopt._relax.nonlinear_bound_tightening) detects variables that appear solely as the bare argument of a periodic function and restricts each one to a single 2*pi window. The reduction is lossless: because the variable reaches the rest of the model only through (sin theta, cos theta), any 2*pi-wide window reproduces every (sin, cos) value the original box could take. A free variable is anchored to [-pi, pi]; a one-sided-finite box is closed to one period from its finite anchor; a > 2*pi finite box is shrunk to its first period. A variable used anywhere outside sin/cos, or inside a non-bare argument such as cos(2*theta), is conservatively skipped, so the pass only ever shrinks a box it can prove periodic-only.

Below, theta is unbounded on both sides and appears only as sin(theta); the pass anchors it to [-pi, pi], which is exactly what lets the solve converge.

import numpy as np

from discopt._relax.nonlinear_bound_tightening import (
    build_flat_variable_metadata,
    tighten_nonlinear_bounds,
)

m = do.Model("periodic")
theta = m.continuous("theta")          # free on both sides: [-inf, +inf]
x = m.continuous("x", lb=-2.0, ub=2.0)
m.minimize(x)
m.subject_to(x - do.sin(theta) >= 0.7)  # theta enters ONLY as sin(theta)

# Apply the nonlinear bound rules to the (unbounded) angular variable's box.
meta = build_flat_variable_metadata(m)
lb = np.array([-np.inf, -2.0])
ub = np.array([np.inf, 2.0])
new_lb, new_ub, stats = tighten_nonlinear_bounds(m, lb, ub)
print(f"theta before: [{lb[0]}, {ub[0]}]")
print(f"theta after : [{new_lb[0]:.4f}, {new_ub[0]:.4f}]   (one period, width {new_ub[0] - new_lb[0]:.4f} = 2*pi)")
print(f"applied rules: {stats.applied_rules}")

theta before: [-inf, inf]
theta after : [-3.1416, 3.1416]   (one period, width 6.2832 = 2*pi)
applied rules: ('periodic_variable_bound',)

With theta anchored to a single period, the otherwise-free angular variable is finitely bounded and the global solve converges (the minimum of x is min sin(theta) + 0.7 = -0.3).

result = m.solve()
print(f"status: {result.status}, objective: {float(result.objective):.4f}")

status: optimal, objective: -0.3000

/Users/jkitchin/projects/discopt/python/discopt/debug/__init__.py:98: UserWarning: Variables with very large or infinite declared bounds: theta (lb=-9.999e+19 [default], ub=9.999e+19 [default]). Nonlinear tightening can adjust 2 bounds via periodic_variable_bound. A bound marked [default] is the box applied to a column you declared with no bounds: 9.999e+19, which is FINITE (it sits just below the 1e+20 infinity sentinel), so the solve can return a certified 'optimal' sitting on that corner rather than 'unbounded'. NLP solvers may fail (NaN, iteration_limit) when bounds exceed ~1e15. Add tighter explicit bounds, e.g. m.continuous('x', lb=0, ub=1000).
  return fn(*args, **kwargs)

What’s next#

The roadmap (in crates/discopt-core/src/presolve/ROADMAP.md) documents the four-phase rollout (P1 foundations, P2 quick wins, P3 structural detection, P4 advanced) and ties each pass to its literature source. Items beyond detection — e.g. orbital branching for D4, automatic symmetry-breaking constraint emission, decomposition dispatch on D5 — are tracked there for future work.