Integer Programming Techniques

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  z_k, x_k
      if LP infeasible or z_k >= z_P:  prune k; continue
      if x_k integral:                 z_P = ceil(x_jk) to L
  z_D  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))  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 $fj$, assignment cost $c{ij}$:

$$ \min \sum{j} fj yj + \sum{i}\sum{j} c{ij} x{ij} \quad \text{s.t.} \quad \sumj x{ij} = 1 \;\; \forall i, \qquad yj \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: $\sumi x{ij} \le m \, y_j \;\; \forall j$
• Strong (disaggregated), $mn$ rows: $x{ij} \le yj \;\; \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 $yj = \frac{1}{m}\sumi 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 $yj \ge \maxi 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))  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.

```python
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] = 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:
    """Sin

…

## Source & license

This open-source skill is cataloged on AgentStack and links to its original source — we do not rehost the code.

- **Author:** [hajibabaie](https://github.com/hajibabaie)
- **Source:** [hajibabaie/combinatorial-optimization-skills](https://github.com/hajibabaie/combinatorial-optimization-skills)
- **License:** MIT

Install and usage instructions live in the source repository linked above.