integer-programming-techniques skill

Integer Programming Techniques

You are an expert in computational integer programming: what a modern MIP solver does with a model between optimize() and OPTIMAL, and how to change the model so the solver finishes sooner. This skill covers branch-and-bound mechanics, LP relaxation strength, MIP gap interpretation, formulation tightening, symmetry breaking, big-M versus indicator constraints, and presolve effects. Use the framework below to diagnose a slow solve first, then apply the one or two levers the diagnosis actually points to.

Initial Assessment

Establish the following before recommending any change:

• The symptom, precisely. "Slow" is not a diagnosis. Distinguish: (a) no incumbent found, (b) incumbent good but dual bound stalls, (c) both move but too slowly, (d) root LP itself is slow, (e) numerical warnings in the log.
• The solve log. Ask for it or reproduce it. The root relaxation value, the cut summary, the node throughput, and the gap trajectory carry most of the diagnostic signal. Never tune blind.
• Problem size. Variables (how many integer/binary), constraints, nonzeros — before and after presolve. A model with 10^7 nonzeros has different bottlenecks than one with 10^4.
• Instance scaling. One instance or a family? Does difficulty explode at a specific size? Collect 3–5 representative instances for any comparison.
• Solver and license. Gurobi version and parameter defaults matter; conclusions below assume a recent Gurobi but transfer to CPLEX/SCIP/HiGHS.
• Time budget and gap target. Proving optimality to 0.01% and reaching 1% feasible-with-certificate are different projects. Get the real requirement.
• Hard vs soft constraints. Soft constraints moved into the objective with penalty weights change relaxation strength; note which constraints are negotiable before tightening anything.
• Data magnitudes. Largest and smallest objective and constraint coefficients. Ratios above ~1e6 within a row or column predict numerical trouble and weak big-M relaxations.
• Structural symmetry. Identical machines, vehicles, bins, shifts, or facilities with equal data are a red flag for symmetric search trees.
• Exact-vs-heuristic need. If a proof of optimality is not required and the gap target is loose, a matheuristic or a warm-started truncated solve may beat months of formulation work; confirm the deliverable first.
• Data format and reproducibility. Fixed random seed, fixed Threads, and pinned solver version for any before/after claim.

How Branch-and-Bound Actually Runs Your Model

Modern solvers run LP-based branch-and-cut (Land & Doig (1960), "An automatic method of solving discrete programming problems"; Achterberg & Wunderling (2013), "Mixed integer programming: Analyzing 12 years of progress"):

BRANCH-AND-CUT (minimization)
  presolve the model                      # tighten bounds, drop rows/cols, probe
  L <- {root};  z_P <- +inf               # open nodes, primal bound (incumbent)
  while L not empty:
      select node k from L                # best-bound / depth-first mix
      solve LP relaxation of k  -> z_k, x_k
      if LP infeasible or z_k >= z_P:  prune k; continue
      if x_k integral:                 z_P <- z_k; update incumbent; continue
      try cutting planes at k (root: aggressively); resolve LP
      run primal heuristics from x_k      # rounding, diving, RINS, ...
      pick fractional variable x_j        # pseudocost / reliability branching
      add children k1: x_j <= floor(x_jk), k2: x_j >= ceil(x_jk) to L
  z_D <- min over open nodes of z_k       # global dual bound at any time

Bounds and the MIP gap. For a minimization MIP with feasible set $X = P \cap (\mathbb{Z}^p \times \mathbb{R}^{n-p})$ and LP relaxation over the polyhedron $P$:

$$ z_{LP} ;=; \min{c^\top x : x \in P} ;\le; z_D ;\le; z^* ;\le; z_P $$

The solver reports ObjVal ($z_P$, incumbent), ObjBound ($z_D$, proven dual bound), and

$$ \text{MIPGap} ;=; \frac{|z_P - z_D|}{|z_P|}. $$

Gurobi stops at MIPGap ≤ 1e-4 by default. A 1% gap means the incumbent is proven within 1% of optimal — it says nothing about which side is weak. The log tells you: if the incumbent column froze long ago, the dual bound is the laggard (formulation problem); if the bound is near-final but the incumbent is poor, it is a primal-heuristic problem (see warm-starts-and-initial-solutions).

Formulation strength. Two formulations with polyhedra $P_1, P_2$ and the same integer feasible set $X$ are ranked by inclusion. The ideal formulation is $\operatorname{conv}(X)$, where the LP optimum is already integral:

$$ \operatorname{conv}(X) ;\subseteq; P_1 ;\subseteq; P_2 ;\Longrightarrow; z_{LP}(P_1) ;\ge; z_{LP}(P_2) \quad (\text{minimization}). $$

A node is pruned as soon as its LP bound reaches $z_P$; a tighter $P$ raises every node bound, so the pruned region of the tree grows. This is why formulation strength, not constraint count, predicts solve time (Wolsey (1998), Integer Programming; Vielma (2015), "Mixed integer linear programming formulation techniques").

The levers and what each one moves:

Lever	Raises root bound	Shrinks tree	Main cost
Disaggregate implications	yes	yes	more rows, slower node LP
Tight big-M from data	yes	yes	bound analysis per constraint
Valid inequalities / cuts	yes	yes	separation effort
Symmetry breaking	rarely	yes, on symmetric trees	can clash with heuristics
Branching priorities	no	often	needs domain insight
MIP start / better heuristics	no (primal side)	yes, earlier pruning	cheap
Presolve level (0/1/2)	sometimes	usually	rarely hurts; test it

Decision guidance:

• Use disaggregation when one aggregated row encodes many independent implications (the facility-location example below is the canonical case).
• Use a tight big-M when a valid $M$ follows from variable bounds or data (remaining demand, capacity); use an indicator constraint when no reasonable $M$ exists or when $M \ge 10^6$ forces numerical trouble.
• Use symmetry breaking only after confirming symmetry (equal-data resources) and only one scheme at a time; benchmark against the solver's own Symmetry=2 handling.
• Treat presolve as part of the formulation: always compare formulations both with Presolve=0 (raw strength) and default (what the solver sees).

Diagnostic Toolkit: Measure Before You Tighten

The reusable workflow. Every tightening claim must be backed by these numbers, not by intuition:

TIGHTEN-AND-DIAGNOSE WORKFLOW (minimization)
1. Solve the pure LP relaxation (model.relax())            -> z_LP
2. Solve the MIP with a short time limit                   -> z_P, z_D, gap, nodes
3. LP integrality gap = (z_P - z_LP) / z_P
   - large (> ~5%): the formulation is the problem ->
       disaggregate, shrink big-M, add valid inequalities
   - small, but the tree is huge: suspect symmetry or weak branching ->
       break symmetry, set branching priorities
   - small and tree small, but nodes are slow: node LP cost is the problem ->
       prefer the smaller formulation, move dense rows to lazy constraints
4. Re-run with fixed Seed; repeat on >= 3 instances and >= 3 seeds
   before declaring a winner (performance variability is real).

import gurobipy as gp
from gurobipy import GRB


def lp_relaxation_bound(model: gp.Model) -> float:
    """Solve the pure LP relaxation of a MIP and return its objective value.

    Caveat: Model.relax() drops integrality AND general constraints
    (indicator, SOS, min/max). For indicator-based models measure the root
    bound from the solve log instead.
    """
    relaxed = model.relax()
    relaxed.Params.OutputFlag = 0
    relaxed.optimize()
    if relaxed.Status != GRB.OPTIMAL:
        raise RuntimeError(f"LP relaxation ended with status {relaxed.Status}")
    return relaxed.ObjVal


def solve_and_report(model: gp.Model, time_limit: float = 60.0) -> dict[str, float]:
    """Solve a MIP and report the numbers that matter for formulation diagnosis."""
    z_lp = lp_relaxation_bound(model)
    model.Params.OutputFlag = 0
    model.Params.TimeLimit = time_limit
    model.Params.Seed = 0  # fix the seed whenever you compare formulations
    model.optimize()
    usable = model.Status == GRB.OPTIMAL or (
        model.Status == GRB.TIME_LIMIT and model.SolCount > 0
    )
    if not usable:
        raise RuntimeError(f"no usable solution: status {model.Status}")
    z_inc = model.ObjVal      # primal bound (incumbent)
    z_dual = model.ObjBound   # dual bound (proven)
    return {
        "lp_relaxation": z_lp,
        "incumbent": z_inc,
        "dual_bound": z_dual,
        "mip_gap": model.MIPGap,
        "lp_integrality_gap": abs(z_inc - z_lp) / max(abs(z_inc), 1e-9),
        "nodes": model.NodeCount,
        "runtime_s": model.Runtime,
    }


# Tiny demo: a 4-item 0-1 knapsack whose LP optimum is fractional.
demo = gp.Model("knapsack_demo")
x = demo.addVars(4, vtype=GRB.BINARY, name="pick")
value = [8.0, 11.0, 6.0, 4.0]
weight = [5.0, 7.0, 4.0, 3.0]
demo.addConstr(gp.quicksum(weight[i] * x[i] for i in range(4)) <= 14, name="cap")
demo.setObjective(gp.quicksum(value[i] * x[i] for i in range(4)), GRB.MAXIMIZE)
print({k: round(v, 4) for k, v in solve_and_report(demo).items()})
# Expected: lp_relaxation 22.0, incumbent 21.0 (items 2,3,4), mip_gap 0.0,
# lp_integrality_gap ~0.0476 -- a 4.8% root gap that branching must close.

The comparison harness solves the same instance under several formulations and returns one tidy row per model:

from collections.abc import Callable

import gurobipy as gp
import pandas as pd
from gurobipy import GRB


def compare_formulations(
    builders: dict[str, Callable[[], gp.Model]],
    time_limit: float = 120.0,
    presolve: int = -1,
) -> pd.DataFrame:
    """Solve one instance under several formulations; one result row per model.

    presolve=0 exposes raw formulation strength; presolve=-1 (solver default)
    shows what branch-and-bound actually works with. Report both.
    """
    rows = []
    for label, build in builders.items():
        mip = build()
        relaxed = mip.relax()
        relaxed.Params.OutputFlag = 0
        relaxed.optimize()
        z_lp = relaxed.ObjVal if relaxed.Status == GRB.OPTIMAL else float("nan")

        mip.Params.OutputFlag = 0
        mip.Params.TimeLimit = time_limit
        mip.Params.Presolve = presolve
        mip.Params.Seed = 0
        mip.optimize()
        has_sol = mip.SolCount > 0
        rows.append(
            {
                "formulation": label,
                "rows": mip.NumConstrs,
                "lp_bound": round(z_lp, 2),
                "objective": round(mip.ObjVal, 2) if has_sol else float("nan"),
                "dual_bound": round(mip.ObjBound, 2),
                "mip_gap_pct": round(100 * mip.MIPGap, 3) if has_sol else float("inf"),
                "nodes": int(mip.NodeCount),
                "runtime_s": round(mip.Runtime, 2),
            }
        )
    return pd.DataFrame(rows).set_index("formulation")

Worked Example 1: Facility Location, Weak vs Strong Linking

The uncapacitated facility location problem (UFLP) is the canonical demonstration that two formulations with identical integer solutions can have wildly different LP relaxations (Cornuéjols, Nemhauser & Wolsey (1990), "The uncapacitated facility location problem").

Customers $i \in I$ ($|I| = m$), candidate sites $j \in J$ ($|J| = n$), fixed opening cost $f_j$, assignment cost $c_{ij}$:

$$ \min \sum_{j} f_j y_j + \sum_{i}\sum_{j} c_{ij} x_{ij} \quad \text{s.t.} \quad \sum_j x_{ij} = 1 ;; \forall i, \qquad y_j \in {0,1}, ; x_{ij} \in [0,1] $$

($x$ may stay continuous: once the open set is fixed, assigning each customer to its cheapest open site is automatically integral.) The linking constraints come in two algebraically equivalent flavors:

• Weak (aggregated), $n$ rows: $\sum_i x_{ij} \le m , y_j ;; \forall j$
• Strong (disaggregated), $mn$ rows: $x_{ij} \le y_j ;; \forall i, j$

Summing the strong rows over $i$ yields the weak row, so $P_{\text{strong}} \subseteq P_{\text{weak}}$. The converse fails badly: the weak LP sets $y_j = \frac{1}{m}\sum_i x_{ij}$, opening every attractive site at level $\approx 1/m$ and paying almost no fixed cost. Its bound collapses toward the pure assignment cost. The strong LP forces $y_j \ge \max_i x_{ij}$ and is frequently integral or near-integral in practice.

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


def make_uflp_instance(n_customers: int, n_sites: int, seed: int = 0) -> dict:
    """Random UFLP on a 100x100 square: Euclidean assignment costs, fixed costs."""
    rng = np.random.default_rng(seed)
    customers = rng.uniform(0.0, 100.0, size=(n_customers, 2))
    sites = rng.uniform(0.0, 100.0, size=(n_sites, 2))
    c = np.linalg.norm(customers[:, None, :] - sites[None, :, :], axis=2)
    f = rng.uniform(400.0, 900.0, size=n_sites)
    return {"c": c, "f": f, "m": n_customers, "n": n_sites}


def build_uflp(data: dict, linking: str) -> gp.Model:
    """UFLP MIP. linking='weak' aggregates one row per site; 'strong' disaggregates."""
    m, n, c, f = data["m"], data["n"], data["c"], data["f"]
    model = gp.Model(f"uflp_{linking}")
    y = model.addVars(n, vtype=GRB.BINARY, name="open")
    x = model.addVars(m, n, lb=0.0, ub=1.0, name="assign")
    model.addConstrs((x.sum(i, "*") == 1 for i in range(m)), name="serve")
    if linking == "weak":
        model.addConstrs(
            (gp.quicksum(x[i, j] for i in range(m)) <= m * y[j] for j in range(n)),
            name="link_agg",
        )
    elif linking == "strong":
        model.addConstrs(
            (x[i, j] <= y[j] for i in range(m) for j in range(n)),
            name="link_dis",
        )
    else:
        raise ValueError(f"unknown linking mode {linking!r}")
    model.setObjective(
        gp.quicksum(f[j] * y[j] for j in range(n))
        + gp.quicksum(c[i, j] * x[i, j] for i in range(m) for j in range(n)),
        GRB.MINIMIZE,
    )
    return model

Run the comparison twice — once with presolve off to see raw strength, once with defaults to see how much presolve repairs:

# make_uflp_instance, build_uflp: defined in the block above.
# compare_formulations: defined in the diagnostic-toolkit section.
data = make_uflp_instance(n_customers=40, n_sites=15, seed=7)
for presolve in (0, -1):
    table = compare_formulations(
        {
            "weak_aggregated": lambda: build_uflp(data, "weak"),
            "strong_disaggregated": lambda: build_uflp(data, "strong"),
        },
        time_limit=120.0,
        presolve=presolve,
    )
    print(f"\nPresolve={presolve}")
    print(table)
# Expected: identical 'objective' in both rows (same integer optimum). With
# Presolve=0 the strong LP bound lands within ~1-2% of the optimum while the
# weak LP bound is far lower (the LP opens sites at y_j ~ 1/40 and pays almost
# no fixed cost); the weak model needs orders of magnitude more nodes. With
# default presolve the difference shrinks (probing derives some x<=y
# implications) but the strong model still dominates.

The general lesson: whenever one row says "if the aggregate is positive then the switch is on," ask whether it is really $m$ separate implications "$x_{ij} > 0 \Rightarrow y_j = 1$" — and write those directly. The strong model pays with $mn$ rows; for very large instances add the disaggregated rows as lazy constraints or user cuts instead (see cutting-planes-valid-inequalities).

Worked Example 2: Symmetry-Broken Parallel-Machine Scheduling

Makespan minimization on identical parallel machines ($P||C_{\max}$): jobs $j = 1..n$ with processing times $p_j$, machines $k = 1..K$ with no distinguishing data.

$$ \min ; C_{\max} \quad \text{s.t.} \quad \sum_{k} x_{jk} = 1 ;\forall j, \qquad \sum_{j} p_j x_{jk} \le C_{\max} ;\forall k, \qquad x_{jk} \in {0,1} $$

Any permutation of the machine indices maps a feasible solution to another one with the same makespan: the symmetry group is $S_K$, so every solution lives in an orbit of up to $K!$ copies. Plain branch-and-bound re-proves the same dual bound in each copy's subtree (Margot (2010), "Symmetry in integer linear programming"; Sherali & Smith (2001), "Improving discrete model representations via symmetry considerations"). Two valid single-scheme remedies:

• Load ordering: $\sum_j p_j x_{j,k} \ge \sum_j p_j x_{j,k+1}$ for $k = 1..K{-}1$. Valid because every solution has a machine-sorted representative.
• Job-index rule: $x_{jk} = 0$ for all $k > j$ (jobs indexed from 1). Valid because machines can be relabeled in order of the smallest job index they host: job 1's machine becomes machine 1, the machine of the smallest job not on machine 1 becomes machine 2, and so on. Implement by fixing variable bounds, not by adding rows — presolve then deletes the variables.

Symmetry breaking barely moves the root bound — the LP already splits jobs fractionally and achieves $z_{LP} = \max{\sum_j p_j / K,; \max_j p_j}$ once the valid inequality $C_{\max} \ge \max_j p_j$ is added (valid because some machine carries that job entirely; not implied by the rest of the LP). The payoff is a smaller tree, not a better bound.

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


def make_pmcmax_instance(n_jobs: int, n_machines: int, seed: int = 0) -> dict:
    """Random identical-parallel-machine instance with integer processing times."""
    rng = np.random.default_rng(seed)
    p = rng.integers(5, 50, size=n_jobs).astype(float)
    return {"p": p, "n": n_jobs, "k": n_machines}


def build_pmcmax(data: dict, symmetry: str) -> gp.Model:
    """P||Cmax assignment model. symmetry in {'none', 'load_order', 'job_index'}."""
    n, k, p = data["n"], data["k"], data["p"]
    model = gp.Model(f"pmcmax_{symmetry}")
    x = model.addVars(n, k, vtype=GRB.BINARY, name="assign")
    cmax = model.addVar(lb=0.0, name="makespan")
    model.addConstrs((x.sum(j, "*") == 1 for j in range(n)), name="one_machine")
    load = {h: gp.quicksum(p[j] * x[j, h] for j in range(n)) for h in range(k)}
    model.addConstrs((load[h] <= cmax for h in range(k)), name="makespan_def")
    model.addConstr(cmax >= float(p.max()), name="longest_job")  # free tightening
    if symmetry == "load_order":
        model.addConstrs(
            (load[h] >= load[h + 1] for h in range(k - 1)), name="sym_load"
        )
    elif symmetry == "job_index":
        for j in range(k - 1):           # job j may only use machines 0..j
            for h in range(j + 1, k):
                x[j, h].UB = 0.0         # fix the bound; do not add a row
    elif symmetry != "none":
        raise ValueError(f"unknown symmetry mode {symmetry!r}")
    model.setObjective(cmax, GRB.MINIMIZE)
    return model

Benchmark the schemes against each other and against the solver's internal symmetry detection. Turn the internal handling off (Symmetry=0) to measure what your constraints contribute, then let Symmetry=2 compete:

# make_pmcmax_instance, build_pmcmax: defined in the block above.
data = make_pmcmax_instance(n_jobs=15, n_machines=5, seed=4)

configs = [
    ("none + solver sym off", "none", 0),
    ("load_order + solver sym off", "load_order", 0),
    ("job_index + solver sym off", "job_index", 0),
    ("none + solver sym aggressive", "none", 2),
]
for label, mode, sym_param in configs:
    model = build_pmcmax(data, mode)
    model.Params.OutputFlag = 0
    model.Params.TimeLimit = 60
    model.Params.Symmetry = sym_param
    model.Params.Seed = 0
    model.optimize()
    if model.Status not in (GRB.OPTIMAL, GRB.TIME_LIMIT) or model.SolCount == 0:
        raise RuntimeError(f"{label}: status {model.Status}")
    print(
        f"{label:30s} makespan={model.ObjVal:6.1f} "
        f"nodes={int(model.NodeCount):8d} time={model.Runtime:6.2f}s"
    )
# Expected: all four configurations agree on the optimal makespan. With the
# solver's symmetry handling off, 'none' explores the largest tree; both
# breaking schemes cut node counts substantially (often 5-50x on this size).
# 'none + Symmetry=2' is usually competitive with the hand-written schemes --
# always measure before shipping symmetry constraints.

The same pattern applies to identical bins, vehicles, shifts, and any model with interchangeable resource indices. The constraints must select exactly one representative per orbit; selecting zero (over-breaking) makes the model infeasible or suboptimal — see Practical Challenges for the classic mistake of stacking two schemes.

Big-M versus Indicator Constraints

The implication $y = 0 \Rightarrow x = 0$ with $x \in [0, U]$ is usually written $x \le M y$. In the LP relaxation this permits $y = x / M$, so the fixed cost attached to $y$ enters the bound at the rate $f/M$ per unit of $x$ — the bound deteriorates linearly as $M$ grows (Camm, Raturi & Tsubakitani (1990), "Cutting big M down to size"). Two distinct failure modes:

1. Weak relaxation: $M \gg U$ makes the root bound useless.
2. Trickle flow: with the default integrality tolerance IntFeasTol=1e-5, $y = 10^{-5}$ counts as integral 0, yet allows $x \le M \cdot 10^{-5}$ — with $M = 10^6$ that is 10 "free" units. The solver returns an incumbent that violates the business logic (Klotz & Newman (2013), "Practical guidelines for solving difficult mixed integer linear programs").

The fixes, in order of preference: derive the tightest valid $M$ from data (here: remaining demand), or hand the implication to the solver as an indicator constraint and let it derive its own reformulation. Indicators add no relaxation strength of their own — Model.relax() deletes them entirely — but they are numerically safe when no reasonable $M$ exists.

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


def build_fixed_charge(d: list[float], mode: str) -> gp.Model:
    """Single-item fixed-charge production model; mode in {'loose','tight','indicator'}."""
    T = len(d)
    setup_cost, unit_cost, hold_cost = 90.0, 1.0, 2.0
    model = gp.Model(f"fixed_charge_{mode}")
    x = model.addVars(T, lb=0.0, name="produce")
    s = model.addVars(T, lb=0.0, name="stock")
    y = model.addVars(T, vtype=GRB.BINARY, name="setup")
    model.addConstrs(
        ((s[t - 1] if t else 0.0) + x[t] == d[t] + s[t] for t in range(T)),
        name="balance",
    )
    if mode == "loose":
        model.addConstrs((x[t] <= 1e6 * y[t] for t in range(T)), name="link")
    elif mode == "tight":
        rem = np.cumsum(np.asarray(d)[::-1])[::-1]  # demand from t to horizon end
        model.addConstrs((x[t] <= float(rem[t]) * y[t] for t in range(T)), name="link")
    elif mode == "indicator":
        for t in range(T):
            model.addGenConstrIndicator(y[t], 0, x[t] == 0.0, name=f"link[{t}]")
    else:
        raise ValueError(f"unknown mode {mode!r}")
    model.setObjective(
        gp.quicksum(setup_cost * y[t] + unit_cost * x[t] + hold_cost * s[t]
                    for t in range(T)),
        GRB.MINIMIZE,
    )
    return model


d = [40.0, 0.0, 70.0, 10.0, 0.0, 60.0, 30.0, 0.0, 25.0, 50.0]
for mode in ("loose", "tight", "indicator"):
    mip = build_fixed_charge(d, mode)
    relaxed = mip.relax()
    relaxed.Params.OutputFlag = 0
    relaxed.optimize()
    mip.Params.OutputFlag = 0
    mip.optimize()
    if mip.Status != GRB.OPTIMAL:
        raise RuntimeError(f"{mode}: status {mip.Status}")
    print(
        f"{mode:9s} lp_bound={relaxed.ObjVal:8.1f} "
        f"optimum={mip.ObjVal:8.1f} nodes={int(mip.NodeCount)}"
    )
# Expected: all three modes report the same optimal cost. 'tight' shows the