# Integer Programming Techniques > When the user wants to diagnose and fix slow MIP solves — branch-and-bound mechanics inside modern solvers, LP relaxation strength, MIP gap interpretation, formulation tightening, symmetry breaking, big-M versus indicator constraints, and presolve effects. Also use when the user mentions "branch and bound," "MIP gap," "LP relaxation," "symmetry breaking," "tighten formulation," "big-M," or when t… - **Type:** Skill - **Install:** `agentstack add skill-hajibabaie-combinatorial-optimization-skills-integer-programming-techniques` - **Verified:** Yes — security-reviewed for prompt injection and unsafe behavior - **Seller:** [hajibabaie](https://agentstack.voostack.com/s/hajibabaie) - **Installs:** 0 - **Category:** [Search](https://agentstack.voostack.com/c/search) - **Latest version:** 0.1.0 - **License:** MIT - **Upstream author:** [hajibabaie](https://github.com/hajibabaie) - **Source:** https://github.com/hajibabaie/combinatorial-optimization-skills/tree/main/skills/integer-programming-techniques ## Install ```sh agentstack add skill-hajibabaie-combinatorial-optimization-skills-integer-programming-techniques ``` Requires the [AgentStack CLI](https://agentstack.voostack.com/docs/cli). Works with Claude Code, Cursor, and any MCP-compatible agent. ## About # 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"): ```text 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). ``` ```python 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 $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. ```python 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: ```python # 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. ```python 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. ## Pricing - **Free** — Free ## Security capabilities Automated source analysis of v0.1.0 — what this tool can access: - **Network access:** no - **Filesystem access:** no - **Shell / process execution:** no - **Environment & secrets:** no - **Dynamic code execution:** no *"Yes" means the capability is present in the source — more access means more to trust, not that it is unsafe.* ## Versions - **0.1.0** — security scan: passed — Imported from the upstream source. ## Links - Listing page: https://agentstack.voostack.com/l/skill-hajibabaie-combinatorial-optimization-skills-integer-programming-techniques - Seller: https://agentstack.voostack.com/s/hajibabaie - Browse the marketplace: https://agentstack.voostack.com/browse --- Listed on AgentStack — the marketplace for AI agent skills and MCP servers. Every listing is security-reviewed. Creators keep 70%.