Cutting Planes Valid Inequalities

Cutting Planes and Valid Inequalities

You are an expert in polyhedral methods for integer programming: deriving valid inequalities, implementing separation routines, and wiring them into a branch-and-cut solver as user cuts or lazy constraints. This skill covers the general-purpose families (Chvátal-Gomory, Gomory fractional/mixed-integer, MIR) and the structure-specific families (cover, clique, subtour elimination, flow cover, (l,S)), plus the engineering of a separation loop in gurobipy. Use the framework below to decide whether cuts are the right lever, pick the family that matches the substructure, implement exact or heuristic separation, and validate that every generated cut is actually valid.

Initial Assessment

Before deriving or coding any cut, establish the following. Each answer changes the plan.

• Why is the formulation weak? Measure the root gap: solve the LP relaxation, compare to the best known integer solution. Cuts attack a weak relaxation. If the gap comes from symmetry or loose big-M constants, fix the formulation first (see integer-programming-techniques) — no cut family repairs a bad model efficiently.
• Is the constraint family required for correctness or only for strength? Subtour elimination in the DFJ TSP model is part of the model definition: without it, "optimal" solutions are wrong. That is a lazy constraint. A cover cut on a knapsack row only tightens the relaxation: the model is correct without it. That is a user cut. This single distinction drives the callback type, the parameters, and the failure modes.
• What substructure does the model contain? Knapsack rows (sums of bounded variables against a capacity) suggest cover and lifted cover cuts. Pairwise conflicts between binaries suggest clique cuts. Single mixed rows suggest MIR. Connectivity requirements suggest cut-set inequalities. No recognizable structure suggests letting Gomory/MIR machinery inside the solver do the work.
• Has the solver's own cut engine been given a chance? Gurobi separates Gomory, MIR, cover, clique, flow-cover, zero-half, and more by default. Try Params.Cuts = 2 (aggressive) and read the cut summary in the log before writing custom separation. Custom cuts pay off when you know structure the solver cannot see.
• Exact or heuristic separation? Exact separation for covers and cliques is itself NP-hard (a knapsack and a max-weight clique, respectively), but pseudo-polynomial DP or greedy heuristics are standard. Subtour separation is polynomial (min cut). Decide what the time budget per node allows.
• Problem size and cut volume. Estimate how many cuts a round can produce. Thousands of dense cuts can slow each LP reoptimization more than the bound improvement is worth.
• Solver and license. The callback recipes here use gurobipy (cbCut, cbLazy). PySCIPOpt and python-mip/CBC offer equivalents; CP-SAT does not expose user cuts.
• Data integrality. Cover-cut DP separation and pure Gomory arguments assume integer row data. Scale or round the data, or switch to MIR, when coefficients are fractional.
• Validation plan. Every separation bug produces a cut that slices off the optimum, and the solver will not warn you. Keep at least one known feasible (ideally optimal) solution of a small instance and assert that every generated cut is satisfied by it.
• Reporting needs. For research output, record bound-versus-round traces and percentage of root gap closed per cut family, not just final runtimes.

Cut Families and the Separation Problem

Let $X = {x \in \mathbb{Z}^n_+ : Ax \le b}$ (or a mixed-integer set) and let $P = {x \ge 0 : Ax \le b}$ be its LP relaxation. An inequality $\pi^\top x \le \pi_0$ is valid for $X$ if every point of $X$ satisfies it. A cutting plane at the fractional point $x^$ is a valid inequality with $\pi^\top x^ > \pi_0$: it separates $x^$ from $\mathrm{conv}(X)$. The separation problem for a family $\mathcal{F}$ of valid inequalities is: given $x^$, return a member of $\mathcal{F}$ violated by $x^*$ or certify that none exists. By the ellipsoid-method equivalence of separation and optimization (Grötschel, Lovász & Schrijver 1981, "The ellipsoid method and its consequences in combinatorial optimization"), optimizing over the closure of a family is polynomially equivalent to separating it — which is why NP-hard separation problems (covers, cliques) are attacked heuristically in practice.

The master template for general-purpose cuts is the Chvátal-Gomory rounding step: for any multiplier vector $u \ge 0$,

$$ \sum_{j=1}^{n} \lfloor u^\top A_j \rfloor , x_j ;\le; \lfloor u^\top b \rfloor $$

is valid for $X$ (Chvátal 1973, "Edmonds polytopes and a hierarchy of combinatorial problems"). Every family below is either a CG cut for a clever $u$, or its mixed-integer generalization.

The catalog

Key formulas. The Gomory fractional cut comes from a tableau row of an optimal LP basis for a pure integer program: if basic variable $x_{B(i)}$ has fractional value $\bar b_i$ in row $x_{B(i)} + \sum_{j \in N} \bar a_{ij} x_j = \bar b_i$, then with $f_{ij} = \bar a_{ij} - \lfloor \bar a_{ij} \rfloor$ and $f_{i0} = \bar b_i - \lfloor \bar b_i \rfloor$,

$$ \sum_{j \in N} f_{ij}, x_j ;\ge; f_{i0} $$

is valid and cuts off the current LP optimum (Gomory 1958, "Outline of an algorithm for integer solutions to linear programs"). The GMI cut strengthens this for mixed problems: with $f_j = f_{ij}$, $f_0 = f_{i0}$,

$$ \sum_{\substack{j \in N_I:, f_j \le f_0}} \frac{f_j}{f_0} x_j

• \sum_{\substack{j \in N_I:, f_j > f_0}} \frac{1 - f_j}{1 - f_0} x_j
• \sum_{\substack{j \in N_C:, \bar a_{ij} > 0}} \frac{\bar a_{ij}}{f_0} x_j

• \sum_{\substack{j \in N_C:, \bar a_{ij} < 0}} \frac{\bar a_{ij}}{1 - f_0} x_j ;\ge; 1 . $$

The MIR inequality is the single-row engine behind GMI: for the set ${(x, s) \in \mathbb{Z}^n_+ \times \mathbb{R}_+ : a^\top x \le b + s}$ with $f = b - \lfloor b \rfloor > 0$ and $f_j = a_j - \lfloor a_j \rfloor$,

$$ \sum_j \Big( \lfloor a_j \rfloor + \frac{\max(f_j - f,, 0)}{1 - f} \Big) x_j ;\le; \lfloor b \rfloor + \frac{s}{1-f} . $$

Marchand & Wolsey (2001, "Aggregation and mixed integer rounding to solve MIPs") showed that aggregating rows and applying MIR reproduces most cut families solvers use; GMI is exactly MIR applied to a tableau row.

User cuts versus lazy constraints

When cuts help — and when they do not

Cuts help when the root gap is large and traceable to a known substructure; when the constraint family is exponential so it must be generated on demand (connectivity, subtours, combinatorial Benders); and when the same model family is solved repeatedly, so separation engineering amortizes. Cuts do not help when the gap comes from symmetry or weak big-M (reformulate instead), when the LP is already the bottleneck and extra dense rows make every reoptimization slower, or when instances are small enough that branching alone finishes in seconds. Measure with one number: percentage of root gap closed, $(z_{\text{cut}} - z_{\text{LP}}) / (z_{\text{IP}} - z_{\text{LP}})$, and watch for tailing off — the typical pattern is strong improvement for a few rounds, then stalling (branching beats round twelve of cutting; this is why branch-and-cut won over pure cutting planes, Padberg & Rinaldi 1991, "A branch-and-cut algorithm for the resolution of large-scale symmetric traveling salesman problems").

Cutting-Plane Framework and General-Purpose Cuts

The reusable artifact is a root-node separation loop: solve the LP relaxation, call a list of separation routines on the fractional point, add the violated cuts, repeat until no violation or tailing off. The same loop runs inside a solver via callbacks; building it standalone first makes separation routines testable in isolation.

ROOT CUTTING-PLANE LOOP
-----------------------
input: MIP model M, separation routines S1..Sk
1. R <- LP relaxation of M
2. repeat up to max_rounds:
3.     solve R; record bound; x* <- optimal point
4.     cuts <- union of Si(x*) for all i, keep violation > eps
5.     if no cuts: stop                      # x* is in all known closures
6.     sort by violation, keep the top K     # cut management
7.     add chosen cuts to R and to M         # M goes to branch-and-cut after
8. return strengthened M, bound trace

"""Reusable root cutting-plane loop: solve the LP relaxation, separate, add, repeat."""
from __future__ import annotations

from dataclasses import dataclass, field
from typing import Callable

import gurobipy as gp
from gurobipy import GRB

# A cut is (coefficients keyed by variable name, sense "<=" or ">=", rhs).
Cut = tuple[dict[str, float], str, float]
Separator = Callable[[dict[str, float]], list[Cut]]


@dataclass
class CutLoopReport:
    """Round-by-round trace of the root cutting-plane loop."""

    bounds: list[float] = field(default_factory=list)
    cuts_per_round: list[int] = field(default_factory=list)

    @property
    def total_cuts(self) -> int:
        return sum(self.cuts_per_round)


def _violation(cut: Cut, point: dict[str, float]) -> float:
    """Amount by which `point` violates `cut` (positive means violated)."""
    coeffs, sense, rhs = cut
    lhs = sum(c * point[name] for name, c in coeffs.items())
    return lhs - rhs if sense == "<=" else rhs - lhs


def root_cutting_plane_loop(
    model: gp.Model,
    separators: list[Separator],
    max_rounds: int = 50,
    min_violation: float = 1e-6,
    max_cuts_per_round: int = 100,
) -> CutLoopReport:
    """Tighten the LP relaxation of a MIP by iterated separation.

    Every accepted cut is added to the working relaxation and to the original
    MIP (named gcut_*), so the strengthened model can then go to the MIP solve.
    """
    report = CutLoopReport()
    model.update()
    relax = model.relax()
    relax.Params.OutputFlag = 0
    cut_id = 0
    for _ in range(max_rounds):
        relax.optimize()
        if relax.Status != GRB.OPTIMAL:
            break
        report.bounds.append(relax.ObjVal)
        point = {v.VarName: v.X for v in relax.getVars()}
        violated = [
            cut
            for sep in separators
            for cut in sep(point)
            if _violation(cut, point) > min_violation
        ]
        if not violated:
            break
        violated.sort(key=lambda cut: -_violation(cut, point))
        for coeffs, sense, rhs in violated[:max_cuts_per_round]:
            for target in (relax, model):
                expr = gp.quicksum(
                    c * target.getVarByName(name) for name, c in coeffs.items()
                )
                constr = expr <= rhs if sense == "<=" else expr >= rhs
                target.addConstr(constr, name=f"gcut_{cut_id}")
            cut_id += 1
        report.cuts_per_round.append(min(len(violated), max_cuts_per_round))
        relax.update()
        model.update()
    return report


# Tiny instance: max 4a + 3b s.t. 2a + 2b <= 3, a, b binary.
# LP optimum (1, 0.5) with value 5.5; the cover cut a + b <= 1 closes the gap.
def _cover_separator(point: dict[str, float]) -> list[Cut]:
    """Hand-written separator for the single cover inequality a + b <= 1."""
    if point["a"] + point["b"] > 1.0 + 1e-6:
        return [({"a": 1.0, "b": 1.0}, "<=", 1.0)]
    return []


demo = gp.Model("demo")
demo.Params.OutputFlag = 0
a = demo.addVar(vtype=GRB.BINARY, name="a")
b = demo.addVar(vtype=GRB.BINARY, name="b")
demo.addConstr(2 * a + 2 * b <= 3, name="capacity")
demo.setObjective(4 * a + 3 * b, GRB.MAXIMIZE)
trace = root_cutting_plane_loop(demo, [_cover_separator])
print(trace.bounds, trace.total_cuts)
# Expected: bounds [5.5, 4.0] and 1 cut -- one cover cut moves the root
# bound from 5.5 to the integer optimum 4.0 before any branching.

MIR as a callable formula

MIR needs no search: given a $\le$ row over nonnegative integer variables with a continuous slack, the cut is a closed-form rounding. The craft is in which row you feed it — aggregated rows, rows with complemented variables, and rows divided by a chosen coefficient all produce different MIR cuts.

"""Mixed-integer rounding (MIR) cut for a single <= row with a continuous slack."""
import math

import numpy as np


def mir_cut(a: np.ndarray, b: float) -> tuple[np.ndarray, float, float] | None:
    """MIR inequality for X = {x in Z^n_+, s in R_+ : a @ x <= b + s}.

    Returns (alpha, beta, gamma) encoding alpha @ x <= beta + gamma * s, with
    alpha_j = floor(a_j) + max(f_j - f, 0) / (1 - f), beta = floor(b),
    gamma = 1 / (1 - f), where f = b - floor(b) and f_j = a_j - floor(a_j).
    Returns None when b is integral (the MIR inequality reduces to the base row).
    """
    f = b - math.floor(b)
    if f < 1e-9 or f > 1 - 1e-9:
        return None
    f_j = a - np.floor(a)
    alpha = np.floor(a) + np.maximum(f_j - f, 0.0) / (1.0 - f)
    return alpha, float(math.floor(b)), 1.0 / (1.0 - f)


# Tiny instance: the set {x integer >= 0, s >= 0 : x <= 2.5 + s}.
alpha, beta, gamma = mir_cut(np.array([1.0]), 2.5)
print(alpha, beta, gamma)
# Expected: alpha=[1.], beta=2.0, gamma=2.0 -- the MIR cut x <= 2 + 2s.
# It cuts off (x, s) = (2.5, 0) while keeping every mixed-integer point:
# at s = 0 it forces x <= 2, and at s = 0.5 the original bound returns.

Gomory fractional cuts from the optimal basis

The implementation below recovers the optimal simplex tableau ($B^{-1}[A \mid I]$ and $B^{-1}b$) from Gurobi's basis statuses and generates one fractional cut per fractional basic row, expressed back in the structural variables by substituting each slack $s_i = b_i - A_i x$. The slack of every original row and of every Gomory cut is integral at integer points, which is what keeps later rounds valid.

"""Gomory fractional cuts recovered from the optimal LP basis (pure integer program)."""
from __future__ import annotations

import numpy as np
import gurobipy as gp
from gurobipy import GRB

TOL = 1e-7


def _frac(values: np.ndarray) -> np.ndarray:
    """Fractional part, snapping entries within TOL of an integer to zero."""
    f = values - np.floor(values)
    f[f < TOL] = 0.0
    f[f > 1.0 - TOL] = 0.0
    return f


def solve_ip_by_gomory(
    c: np.ndarray, A: np.ndarray, b: np.ndarray, max_rounds: int = 100
) -> tuple[np.ndarray, float, int]:
    """Solve max c @ x s.t. A @ x <= b, x >= 0 integer, by Gomory fractional cuts only.

    Requires integer A and b. The slack of every original row and of every
    Gomory cut is integral at integer points, so each fractional basic row
    of the optimal tableau yields a valid cut. Returns (x, objective, n_cuts).
    """
    A_cur = A.astype(float)
    b_cur = b.astype(float)
    n = c.size
    n_cuts = 0
    for _ in range(max_rounds):
        m_rows = A_cur.shape[0]
        lp = gp.Model("gomory-lp")
        lp.Params.OutputFlag = 0
        lp.Params.Method = 0  # primal simplex: we need a basic optimal solution
        x = lp.addMVar(n, lb=0.0, name="x")
        lp.addConstr(A_cur @ x <= b_cur, name="row")
        lp.setObjective(c @ x, GRB.MAXIMIZE)
        lp.optimize()
        if lp.Status != GRB.OPTIMAL:
            raise RuntimeError(f"LP relaxation not optimal (status {lp.Status})")
        x_val = x.X
        if np.all(np.abs(x_val - np.round(x_val)) < 1e-6):
            return np.round(x_val), float(lp.ObjVal), n_cuts
        # Recover the optimal tableau: rows of B^-1 [A | I] and B^-1 b.
        vbasis = np.array(lp.getAttr("VBasis", lp.getVars()))
        cbasis = np.array(lp.getAttr("CBasis", lp.getConstrs()))
        basic = np.concatenate(
            [np.where(vbasis == GRB.BASIC)[0], n + np.where(cbasis == GRB.BASIC)[0]]
        )
        A_ext = np.hstack([A_cur, np.eye(m_rows)])
        basis_matrix = A_ext[:, basic]
        tableau = np.linalg.solve(basis_matrix, A_ext)
        rhs = np.linalg.solve(basis_matrix, b_cur)
        new_rows: list[np.ndarray] = []
        new_rhs: list[float] = []
        for i in range(m_rows):
            f0 = rhs[i] - np.floor(rhs[i])
            if f0 < 1e-6 or f0 > 1.0 - 1e-6:
                continue
            f = _frac(tableau[i])
            # Cut in (x, s) space: f_x @ x + f_s @ s >= f0; substitute s = b - A x.
            alpha = f[:n] - f[n:] @ A_cur
            beta = f0 - float(f[n:] @ b_cur)
            new_rows.append(-alpha)  # store in <= form
            new_rhs.append(-beta)
        if not new_rows:
            raise RuntimeError("fractional optimum but no cut generated (numerics)")
        A_cur = np.vstack([A_cur] + new_rows)
        b_cur = np.concatenate([b_cur, np.array(new_rhs)])
        n_cuts += len(new_rows)
    raise RuntimeError("round limit reached before integrality")


# Tiny instance: max 5x1 + 4x2 s.t. 6x1 + 4x2 <= 24, x1 + 2x2 <= 6, x integer.
c = np.array([5.0, 4.0])
A = np.array([[6.0, 4.0], [1.0, 2.0]])
b = np.array([24.0, 6.0])
x_opt, obj, cuts_used = solve_ip_by_gomory(c, A, b)
print(x_opt, obj, cuts_used)
# Expected: x = [4, 0] with objective 20.0 after 3 Gomory cuts. The LP
# relaxation was 21.0 at (3, 1.5); cuts close the gap with no branching,
# which is Gomory's 1958 finiteness result in action.

Pure Gomory loops are finite in theory but numerically fragile in long runs — fractional residues accumulate until "cuts" become invalid (Zanette, Fischetti & Balas 2011, "Lexicography and degeneracy: can a pure cutting plane algorithm work?"). Solvers therefore apply a few Gomory passes and hand over to branching; do the same.

Worked Example: Knapsack Cover Cuts with Exact Separation and Lifting

Take a knapsack row $\sum_{j \in N} a_j x_j \le b$ with $x$ binary and integer data. A set $C \subseteq N$ is a cover if $\sum_{j \in C} a_j > b$: all of $C$ cannot be packed, so

$$ \sum_{j \in C} x_j ;\le; |C| - 1 $$

is valid. A cover is minimal if dropping any item destroys the cover property; only minimal covers can be facet-inducing on the restricted polytope, and lifting extends them toward facets of the full polytope. These cuts powered the first large-scale 0-1 successes (Crowder, Johnson & Padberg 1983, "Solving large-scale zero-one linear programming problems"); modern computational practice is surveyed in Gu, Nemhauser & Savelsbergh (1998, "Lifted cover inequalities for 0-1 integer programs").

Exact separation. At the fractional point $x^$, the most violated cover minimizes $\sum_{j \in C}(1 - x^_j)$ over covers; a violated cover exists iff

$$ \zeta = \min \Big{ \sum_j (1 - x^*_j), z_j ;:; \sum_j a_j z_j \ge b + 1,; z \in {0,1}^n \Big} < 1 . $$

This is itself a knapsack — NP-hard in general, solved below by a DP over the capped weight in $O(nb)$ time.

Sequential lifting. For $j \notin C$, the exact lifting coefficient is $\alpha_j = (|C| - 1) - \max{\sum_i \alpha_i x_i : \sum_{i} a_i x_i \le b - a_j}$, maximizing over the variables already in the inequality — again a small knapsack. Lifting order changes the result; descending $x^*_j$ is the standard choice because it strengthens the cut where the LP point puts weight.

"""Exact cover separation for one knapsack row, with sequential lifting."""
from __future__ import annotations

import numpy as np


def separate_cover(
    a: np.ndarray, b: int, x_frac: np.ndarray
) -> tuple[list[int], float] | None:
    """Most violated cover for sum(a_j x_j) <= b at the fractional point x_frac.

    Solves min sum((1 - x*_j) z_j) s.t. sum(a_j z_j) >= b + 1 exactly by a DP
    over capped weight (needs integer a, b; O(n*b) time). Returns the cover as
    a minimal-cover index list plus its violation, or None if none is violated.
    """
    n = int(a.size)
    target = int(b) + 1
    cost = 1.0 - x_frac
    INF = float("inf")
    dp = np.full((n + 1, target + 1), INF)
    dp[0, 0] = 0.0
    take = np.zeros((n + 1, target + 1), dtype=bool)
    prev_w = np.zeros((n + 1, target + 1), dtype=int)
    for j in range(1, n + 1):
        dp[j] = dp[j - 1]
        prev_w[j] = np.arange(target + 1)
        for w in range(target + 1):
            if dp[j - 1, w] == INF:
                continue
            w2 = min(w + int(a[j - 1]), target)
            if dp[j - 1, w] + cost[j - 1] < dp[j, w2]:
                dp[j, w2] = dp[j - 1, w] + cost[j - 1]
                take[j, w2] = True
                prev_w[j, w2] = w
    if dp[n, target] >= 1.0 - 1e-9:
        return None
    cover: list[int] = []
    w = target
    for j in range(n, 0, -1):
        if take[j, w]:
            cover.append(j - 1)
        w = prev_w[j, w]
    # Make the cover minimal; dropping j raises violation by 1 - x*_j >= 0.
    weight = int(a[cover].sum())
    for j in sorted(cover, key=lambda i: x_frac[i]):
        if weight - int(a[j]) >= target:
            cover.remove(j)
            weight -= int(a[j])
    violation = float(x_frac[cover].sum()) - (len(cover) - 1)
    return sorted(cover), violation


def _knapsack_max(profit: np.ndarray, weight: np.ndarray, cap: int) -> int:
    """0-1 knapsack by DP over capacity; returns the max total profit."""
    dp = np.zeros(cap + 1, dtype=int)
    for p, w in zip(profit.astype(int), weight.astype(int)):
        if w <= cap:
            dp[w:] = np.maximum(dp[w:], dp[: cap + 1 - w] + p)
    return int(dp[cap])


def lift_cover(
    cover: list[int], a: np.ndarray, b: int, x_frac: np.ndarray
) -> dict[int, int]:
    """Sequential up-lifting of variables outside a minimal cover.

    For each outside j (in descending x*_j order), the exact lifting
    coefficient is alpha_j = (|C| - 1) - max{current LHS : weight <= b - a_j},
    a small knapsack solved by DP. Returns {index: coefficient} for all
    variables with a nonzero coefficient.
    """
    coeffs = {j: 1 for j in cover}
    rhs = len(cover) - 1
    outside = sorted(
        (j for j in range(int(a.size)) if j not in coeffs),
        key=lambda j: -x_frac[j],
    )
    for j in outside:
        cap = int(b) - int(a[j])
        if cap < 0:  # x_j = 1 is infeasible on its own; any coefficient <= rhs works
            coeffs[j] = rhs
            continue
        items = list(coeffs)
        best = _knapsack_max(
            np.array([coeffs[i] for i in items]),
            np.array([int(a[i]) for i in items]),
            cap,
        )
        if rhs - best > 0:
            coeffs[j] = rhs - best
    return coeffs


# Tiny instance (Wolsey 1998, Integer Programming, ch. 9 flavor):
# 79x1 + 53x2 + 53x3 + 45x4 + 45x5 <= 178, x binary.
a = np.array([79, 53, 53, 45, 45])
b = 178
x_frac = np.array([0.0, 1.0, 1.0, 0.75, 0.5])
cover, violation = separate_cover(a, b, x_frac)
print(cover, round(violation, 4))
print(lift_cover(cover, a, b, x_frac))
# Expected: cover [1, 2, 3, 4] (i.e., x2 + x3 + x4 + x5 <= 3) with violation
# 0.25; lifting then gives x1 the coefficient 1, so the final lifted cover
# inequality is x1 + x2 + x3 + x4 + x5 <= 3.

Cover cuts inside a Gurobi user-cut callback

In the tree, separation runs at MIPNODE on the node relaxation. The callback below uses the standard greedy heuristic (items enter in increasing $(1 - x^*_j)/a_j$ order — the LP greedy for the separation knapsack); swap in the exact DP above when the heuristic dries up. PreCrush = 1 is mandatory so Gurobi can translate cuts onto the presolved model.

"""Cover cuts as Gurobi user cuts on a 0-1 multidimensional knapsack."""
from __future__ import annotations

import numpy as np
import gurobipy as gp
from gurobipy import GRB


def greedy_cover(a: np.ndarray, b: float, x_frac: np.ndarray) -> list[int] | None:
    """Heuristic separation: fill a cover cheapest-violation-first, then minimalize.

    Items enter in increasing (1 - x*_j) / a_j order (the LP-relaxation greedy
    for the exact separation knapsack); integer row data assumed. Returns a
    violated minimal cover or None.
    """
    order = np.argsort((1.0 - x_frac) / a)
    cover: list[int] = []
    weight = 0.0
    for j in order:
        cover.append(int(j))
        weight += float(a[j])
        if weight > b + 0.5:
            break
    else:
        return None  # all items together do not exceed b: no cover exists
    for j in sorted(cover, key=lambda i: x_frac[i]):
        if weight - float(a[j]) > b + 0.5:
            cover.remove(j)
            weight -= float(a[j])
    violation = float(x_frac[cover].sum()) - (len(cover) - 1)
    return sorted(cover) if violation > 1e-6 else None


def cover_cut_callback(model: gp.Model, where: int) -> None:
    """Separate cover cuts on every knapsack row at the root relaxation."""
    if where != GRB.Callback.MIPNODE:
        return
    if model.cbGet(GRB.Callback.MIPNODE_STATUS) != GRB.OPTIMAL:
        return
    if model.cbGet(GRB.Callback.MIPNODE_NODCNT) > 0:
        return  # root only; drop this test to separate throughout the tree
    x_frac = np.array(model.cbGetNodeRel(model._xvars))
    for i in range(model._A.shape[0]):
        cover = greedy_cover(model._A[i], float(model._b[i]), x_frac)
        if cover is not None:
            model.cbCut(
                gp.quicksum(model._xvars[j] for j in cover) <= len(cover) - 1
            )
            model._n_cover_cuts += 1


def solve_mkp_with_cover_cuts(
    profits: np.ndarray, A: np.ndarray, b: np.ndarray
) -> tuple[float, int]:
    """Solve max p @ x, A @ x <= b, x binary, separating our own cover cuts."""
    model = gp.Model("mkp-cover")
    model.Params.OutputFlag = 0
    model.Params.PreCrush = 1  # mandatory for user cuts: map cuts to presolved model
    model.Params.Cuts = 0      # disable Gurobi's own cuts to isolate the effect
    n = profits.size
    x = model.addVars(n, vtype=GRB.BINARY, name="x")
    model.addConstrs(
        (gp.quicksum(float(A[i, j]) * x[j] for j in range(n)) <= float(b[i])
         for i in range(A.shape[0])),
        name="knapsack",
    )
    model.setObjective(
        gp.quicksum(float(profits[j]) * x[j] for j in range(n)), GRB.MAXIMIZE
    )
    model._xvars = [x[j] for j in range(n)]
    model._A, model._b = A, b
    model._n_cover_cuts = 0
    model.optimize(cover_cut_callback)
    if model.Status != GRB.OPTIMAL:
        raise RuntimeError(f"unexpected status {model.Status}")
    return model.ObjVal, model._n_cover_cuts


# Tiny instance: 30 items, 5 knapsack rows, capacities at 35% of row weight.
rng = np.random.default_rng(7)
profits = rng.integers(10, 50, size=30)
A = rng.integers(5, 30, size=(5, 30))
b = (0.35 * A.sum(axis=1)).astype(int)
obj, n_cuts = solve_mkp_with_cover_cuts(profits, A, b)
print(obj, n_cuts)
# Expected: objective 445.0 with about 10 cover cuts separated at the root
# (seed 7; the exact cut count varies slightly across Gurobi versions).

Worked Example: TSP with Lazy Subtour Elimination

The Dantzig-Fulkerson-Johnson formulation of the symmetric TSP on $n$ cities (Dantzig, Fulkerson & Johnson 1954, "Solution of a large-scale traveling-salesman problem") uses one binary $x_e$ per edge:

$$ \min \sum_{e} d_e x_e \quad \text{s.t.} \quad x(\delta(i)) = 2 ;; \forall i, \qquad x(\delta(S)) \ge 2 ;; \forall, \emptyset \ne S \subsetneq V, \qquad x \in {0,1}^E , $$

where $\delta(S)$ is the set of edges with exactly one endpoint in $S$. The cut-set form $x(\delta(S)) \ge 2$ is equivalent, under the degree constraints, to the subtour form $\sum_{e \in E(S)} x_e \le |S| - 1$. There are exponentially many such constraints, so they are lazy: the model starts with degree constraints only (its LP is the 2-matching relaxation), and every candidate integer solution is checked for connectivity. Integer separation is just connected components, $O(n + |E^+|)$; this beats the compact MTZ formulation by orders of magnitude on nontrivial instances (see traveling-salesman-problem for the comparison and for heuristics).

"""Symmetric TSP, DFJ formulation: lazy subtour-elimination cuts in a callback."""
from __future__ import annotations

import itertools

import numpy as np
import gurobipy as gp
from gurobipy import GRB


def euclidean_instance(n: int, seed: int) -> np.ndarray:
    """Random points in the unit square; returns the symmetric distance matrix."""
    rng = np.random.default_rng(seed)
    pts = rng.random((n, 2))
    diff = pts[:, None, :] - pts[None, :, :]
    return np.sqrt((diff**2).sum(axis=2))


def components(edges: list[tuple[int, int]], n: int) -> list[list[int]]:
    """Connected components of the support graph of the selected edges."""
    adj: list[list[int]] = [[] for _ in range(n)]
    for i, j in edges:
        adj[i].append(j)
        adj[j].append(i)
    seen = [False] * n
    comps: list[list[int]] = []
    for start in range(n):
        if seen[start]:
            continue
        stack, comp = [start], []
        seen[start] = True
        while stack:
            u = stack.pop()
            comp.append(u)
            for v in adj[u]:
                if not seen[v]:
                    seen[v] = True
                    stack.append(v)
        comps.append(comp)
    return comps


def solve_tsp(dist: np.ndarray, time_limit: float = 60.0) -> tuple[list[int], float]:
    """Exact symmetric TSP: degree constraints up front, DFJ cuts added lazily."""
    n = dist.shape[0]
    model = gp.Model("tsp-dfj")
    model.Params.OutputFlag = 0
    model.Params.LazyConstraints = 1  # mandatory whenever a callback calls cbLazy
    model.Params.TimeLimit = time_limit
    x = model.addVars(
        itertools.combinations(range(n), 2), vtype=GRB.BINARY, name="x"
    )
    model.addConstrs(
        (
            gp.quicksum(x[min(i, j), max(i, j)] for j in range(n) if j != i) == 2
            for i in range(n)
        ),
        name="degree",
    )
    model.setObjective(
        gp.quicksum(dist[i, j] * x[i, j] for i, j in x), GRB.MINIMIZE
    )
    model._n_sec = 0

    def subtour_callback(m: gp.Model, where: int) -> None:
        """Reject every integer point whose support graph is disconnected."""
        if where != GRB.Callback.MIPSOL:
            return
        vals = m.cbGetSolution(x)
        chosen = [(i, j) for (i, j), v in vals.items() if v > 0.5]
        comps = components(chosen, n)
        if len(comps) == 1:
            return  # a single Hamiltonian cycle: accept
        for comp in comps:
            m.cbLazy(
                gp.quicksum(
                    x[i, j] for i, j in itertools.combinations(sorted(comp), 2)
                )
                <= len(comp) - 1
            )
            m._n_sec += 1

    model.optimize(subtour_callback)
    if model.SolCount == 0:
        raise RuntimeError(f"no tour found (status {model.Status})")
    chosen = [(i, j) for (i, j), v in x.items() if v.X > 0.5]
    # Independent validation: degree 2 everywhere and a single component.
    degree = [0] * n
    for i, j in chosen:
        degree[i] += 1
        degree[j] += 1
    assert all(d == 2 for d in degree), "degree constraint violated"
    assert len(components(chosen, n)) == 1, "subtour survived the callback"
    adj = {i: [] for i in range(n)}
    for i, j in chosen:
        adj[i].append(j)
        adj[j].append(i)
    tour, prev = [0], -1
    while len(tour) < n:
        nxt = adj[tour[-1]][0] if adj[tour[-1]][0] != prev else adj[tour[-1]][1]
        prev = tour[-1]
        tour.append(nxt)
    length = sum(dist[tour[k], tour[(k + 1) % n]] for k in range(n))
    print(f"SEC added: {model._n_sec}, status: {model.Status}")
    return tour, length


dist = euclidean_instance(n=30, seed=42)
tour, length = solve_tsp(dist)
print(len(tour), round(length, 4))
# Expected: optimal tour over 30 cities with length 4.5698 (seed 42), after
# roughly two dozen lazy subtour cuts (the exact count varies with Gurobi
# version and heuristics; status 2 = OPTIMAL either way).

Three implementation rules for lazy constraints. First, the callback must reject every violating integer point, including those proposed by Gurobi's heuristics — never assume MIPSOL points come from your own branching. Second, add one cut per component (the complement set gives the same cut-set, so skip it or dedupe). Third, do not read .X inside the callback; only cbGetSolution/cbGetNodeRel are valid there.

Advanced Techniques

Fractional subtour separation by minimum cut

Waiting for integer points wastes information: a fractional LP point can already violate $x(\delta(S)) \ge 2$, and separation is polynomial — a global minimum cut on the support graph with capacities $x^*_e$. Any violated $S$ separates node 0 from some node $t$, so $n - 1$ max-flow calls suffice; Padberg & Rinaldi (1990, "An efficient algorithm for the minimum capacity cut problem") add shrinking to make this fast at scale. Add the resulting sets at MIPNODE via cbCut (they are user cuts there — the lazy MIPSOL check still guarantees correctness).

"""Fractional subtour separation: min cuts on the LP support graph (no solver needed)."""
from __future__ import annotations

import numpy as np


def max_flow_min_cut(cap: np.ndarray, s: int, t: int) -> tuple[float, np.ndarray]:
    """Edmonds-Karp max flow on a dense capacity matrix.

    Returns (flow value, boolean mask of the source side of a min s-t cut).
    O(V * E^2) worst case; fine for LP support graphs, which are sparse.
    """
    n = cap.shape[0]
    residual = cap.astype(float).copy()
    flow = 0.0
    while True:
        parent = np.full(n, -1)
        parent[s] = s
        queue = [s]
        while queue and parent[t] == -1:
            u = queue.pop(0)
            for v in np.where(residual[u] > 1e-9)[0]:
                if parent[v] == -1:
                    parent[v] = u
                    queue.append(int(v))
        if parent[t] == -1:
            break
        bottleneck = float("inf")
        v = t
        while v != s:
            u = int(parent[v])
            bottleneck = min(bottleneck, residual[u, v])
            v = u
        v = t
        while v != s:
            u = int(parent[v])
            residual[u, v] -= bottleneck
            residual[v, u] += bottleneck
            v = u
        flow += bottleneck
    side = np.zeros(n, dtype=bool)
    side[s] = True
    stack = [s]
    while stack:
        u = stack.pop()
        for v in np.where((residual[u] > 1e-9) & ~side)[0]:
            side[v] = True
            stack.append(int(v))
    return flow, side


def violated_subtour_sets(
    x_frac: dict[tuple[int, int], float], n: int, tol: float = 1e-6
) -> list[set[int]]:
    """All distinct sets S with x*(delta(S)) < 2 - tol, via n-1 max-flow calls.

    Exact separation for the DFJ cut-set inequalities x(delta(S)) >= 2: any
    violated S separates node 0 from some t, so min cuts (0 -> t) for every
    t suffice (Padberg & Rinaldi 1990 refine this with graph shrinking).
    """
    cap = np.zeros((n, n))
    for (i, j), v in x_frac.items():
        if v > tol:
            cap[i, j] = cap[j, i] = v
    found: list[set[int]] = []
    seen: set[frozenset[int]] = set()
    for t in range(1, n):
        flow, side = max_flow_min_cut(cap, 0, t)
        if flow < 2.0 - tol:
            s_set = frozenset(int(v) for v in np.where(side)[0])
            small = (
                s_set
                if len(s_set) <= n - len(s_set)
                else frozenset(range(n)) - s_set
            )
            if 2 <= len(small) <= n - 2 and small not in seen:
                seen.add(small)
                found.append(set(small))
    return found


# Tiny instance: 6 nodes, the LP point puts two disjoint triangles at value 1.
# Degree constraints hold (every node has incident weight 2), but the point
# is far outside the TSP polytope.
x_lp = {
    (0, 1): 1.0, (1, 2): 1.0, (0, 2): 1.0,
    (3, 4): 1.0, (4, 5): 1.0, (3, 5): 1.0,
}
print(violated_subtour_sets(x_lp, n=6))
# Expected: [{0, 1, 2}] -- one violated set (its complement gives the same
# cut), yielding the user cut x(delta({0,1,2})) >= 2.

Clique cuts from conflict graphs

When the model implies pairwise conflicts $x_i + x_j \le 1$ (from set-packing rows, from probing on big-M rows, or from one-machine overlap), any clique $Q$ in the conflict graph gives $\sum_{j \in Q} x_j \le 1$ — strictly dominating the edge inequalities it merges. Exact separation is max-weight clique (NP-hard); the greedy below grows a clique in descending-$x^*$ order, which is the standard cheap heuristic.

"""Clique cut separation on a conflict graph by greedy weighted clique growing."""
from __future__ import annotations

import numpy as np


def separate_clique(
    conflicts: np.ndarray, x_frac: np.ndarray, tol: float = 1e-6
) -> list[int] | None:
    """Grow a clique greedily in descending-x* order; report it if violated.

    conflicts[i, j] = True means x_i + x_j <= 1 is implied by the model.
    Returns a clique Q with sum(x*[Q]) > 1 + tol -- the violated clique cut
    sum_{j in Q} x_j <= 1 -- or None. Exact separation is NP-hard.
    """
    order = np.argsort(-x_frac, kind="stable")
    clique: list[int] = []
    for v in order:
        if x_frac[v] <= tol:
            break
        if all(conflicts[v, u] for u in clique):
            clique.append(int(v))
    weight = float(x_frac[clique].sum())
    return clique if weight > 1.0 + tol else None


# Tiny instance: conflicts (0,1), (0,2), (1,2), (2,3); x* = (.5, .5, .5, .3, 0).
n = 5
conflicts = np.zeros((n, n), dtype=bool)
for i, j in [(0, 1), (0, 2), (1, 2), (2, 3)]:
    conflicts[i, j] = conflicts[j, i] = True
x_frac = np.array([0.5, 0.5, 0.5, 0.3, 0.0])
print(separate_clique(conflicts, x_frac))
# Expected: [0, 1, 2] -- every pair sums to exactly 1.0, so the point obeys
# all edge inequalities, yet the triangle has weight 1.5 > 1: the clique cut
# x0 + x1 + x2 <= 1 is violated by 0.5.

(l,S) inequalities for lot-sizing

For uncapacitated lot-sizing with production $x_t$, setup $y_t$, demand $d_t$, and $d_{tl} = \sum_{u=t}^{l} d_u$: for every $l$ and every $S \subseteq L = {1, \dots, l}$,

$$ \sum_{t \in L \setminus S} x_t ;+; \sum_{t \in S} d_{tl}, y_t ;\ge; d_{1l} . $$

Demand through $l$ is met either by production in periods outside $S$ or by periods in $S$, each capped at $d_{tl}$ once its setup is paid. Separation is exact and $O(n^2)$: for

…(truncated)