Integer Programming and LP, in depth: LP relaxations, branch and bound, cutting planes and models a solver can finish

The branch-and-bound tree in one picture

Linear programs and why integers change everything

Write a linear program as: maximise c^Tx subject to Ax ≤ b and x ≥ 0. Each constraint is a half-space, and the feasible region is a convex polyhedron. If an optimum exists, one occurs at a vertex. The simplex method walks between adjacent vertices; interior-point methods, which are provably polynomial, cut through the middle. Both handle millions of variables in practice.

Every linear program also has a dual whose optimal value equals the primal's whenever either has an optimum. Any feasible dual solution bounds the primal objective, which is how solvers prove optimality, and the optimal dual values, the shadow prices, say which capacity is worth buying.

Now require x to be integer. The feasible set becomes the lattice points inside the polyhedron, and the optimum need not be near the LP vertex. Rounding can produce infeasible points, or feasible points far from optimal. Binary variables can encode satisfiability, so integer programming is NP-hard.

The LP relaxation and the integrality gap

The LP relaxation of an integer program is the same model with the integrality requirements dropped. Because it allows more solutions, its optimum is at least as good as the integer optimum: for a maximisation, the relaxation's value is an upper bound. The difference between the two is the integrality gap, and much of the art of integer programming is writing models where that gap is small, because the solver's search effort grows with it.

Some problems have no gap at all. If the constraint matrix is totally unimodular, meaning every square submatrix has determinant 0, 1 or -1, and the right-hand side is integer, every vertex of the polyhedron is integer, so the plain LP returns an integer optimum. Network flow constraints and bipartite matching constraints have this property, which is why maximum flow with integer capacities always has an integer optimal flow and why bipartite matching can be solved as a linear program. When you can make a model look like a network, do.

The LP relaxation of 0/1 knapsack is fractional knapsack, which a greedy algorithm solves, and it is exactly the bound used by the branch-and-bound solver in the knapsack deep dive.

Branch and bound

Branch and bound searches integer solutions as a tree. Solve the root relaxation; if the solution is integer, it is optimal. Otherwise pick a fractional variable, say x_i = 3.4, and create two children with x_i ≤ 3 and x_i ≥ 4. Every integer solution lies in exactly one child, the fractional point in neither. The best integer solution so far is the incumbent, and any node whose bound cannot beat it is pruned with its whole subtree.

The search is complete when no open node has a bound better than the incumbent. The best bound among open nodes and the incumbent together give a gap, which is what a solver reports while it runs: an incumbent of 980 with a best bound of 1,000 means the current answer is within 2% of optimal, even if the solver never finishes. The minimal implementation below uses SciPy's linprog for the relaxations, branches on the most fractional variable and searches depth first, which finds incumbents early so that pruning starts sooner.

import math
import numpy as np
from scipy.optimize import linprog

def branch_and_bound(c, A, b, bounds, tol=1e-6):
    """Maximise c @ x subject to A @ x <= b, bounds, x integer."""
    best_val, best_x = -math.inf, None
    stack = [list(bounds)]                      # each node is a list of (lo, hi)
    nodes = 0
    while stack:
        node = stack.pop()                      # depth first finds incumbents early
        nodes += 1
        res = linprog(-np.asarray(c), A_ub=A, b_ub=b, bounds=node, method="highs")
        if res.status != 0:
            continue                            # infeasible node: prune
        val = -res.fun
        if val <= best_val + tol:
            continue                            # bound no better than incumbent: prune
        frac = [i for i, v in enumerate(res.x) if abs(v - round(v)) > tol]
        if not frac:
            best_val, best_x = val, np.round(res.x)   # integer: new incumbent
            continue
        i = min(frac, key=lambda k: abs(res.x[k] - math.floor(res.x[k]) - 0.5))  # most fractional
        lo, hi = node[i]
        down, up = list(node), list(node)
        down[i] = (lo, math.floor(res.x[i]))
        up[i] = (math.ceil(res.x[i]), hi)
        stack += [down, up]
    return best_val, best_x, nodes

c = [5, 4]
A = [[6, 4], [1, 2]]
b = [24, 6]
print(branch_and_bound(c, A, b, [(0, None), (0, None)]))

Real solvers add smarter branching, warm-started child LPs, primal heuristics and presolve, but the skeleton is the right mental model for reading their logs.

Worked example: solving a two-variable IP by hand

Maximise 5x₁ + 4x₂ subject to 6x₁ + 4x₂ ≤ 24, x₁ + 2x₂ ≤ 6, with x₁, x₂ non-negative integers. Think of two products sharing a machine with 24 hours and a material supply of 6 units.

Root. The relaxation's optimum is where both constraints are tight: solving 6x₁ + 4x₂ = 24 and x₁ + 2x₂ = 6 gives x = (3, 1.5) and z = 21. So no integer solution beats 21. Branch on x₂, the fractional variable.

Right child, x₂ ≥ 2. With the material constraint tight, x₁ = 6 - 2x₂ and z = 30 - 6x₂, which falls as x₂ grows, so the optimum is x = (2, 2) with z = 18. It is integer: the first incumbent is 18.

Left child, x₂ ≤ 1. Per hour of machine time, x₂ earns 4/4 = 1 and x₁ earns 5/6, so the relaxation sets x₂ = 1 and spends the remaining 20 hours on x₁ = 3.33, giving z = 20.67. That beats 18, so branch on x₁. With x₁ ≤ 3 the optimum is (3, 1) with z = 19, a better incumbent. With x₁ ≥ 4 the machine constraint forces x₂ = 0 and z = 20, integer and better again. No open nodes remain, so x = (4, 0) with z = 20 is optimal.

Rounding the LP optimum (3, 1.5) gives (3, 2), which is infeasible, or (3, 1), worth only 19. And once the bound was 20.67, finding 20 ended the search, because integer coefficients give integer objective values and none lies between 20 and 20.67.

Cutting planes

Branching splits the region; cutting planes shrink it. A cut is an inequality every integer feasible point satisfies but the current fractional optimum violates, so adding it tightens the relaxation without losing any integer solution.

The simplest is the Chvátal-Gomory cut: add non-negative multiples of constraints, round down every left-hand coefficient (valid since x ≥ 0), then round down the right-hand side (valid since the left side is now an integer). In the example, take one eighth of the machine constraint and one quarter of the material constraint: (6/8 + 1/4)x₁ + (4/8 + 2/4)x₂ ≤ 3 + 1.5, which is x₁ + x₂ ≤ 4.5. Integer points satisfy x₁ + x₂ ≤ 4. The LP optimum (3, 1.5) sums to 4.5 and is cut off.

Re-solve with the cut. The vertex where x₁ + x₂ = 4 meets the machine constraint is (4, 0) with z = 20, and the vertex where it meets the material constraint is (2, 2) with z = 18. The relaxation's optimum is now (4, 0), which is integer, so one cut solved the problem with no branching. That rarely happens on real models, but modern solvers combine cuts with branching, an approach called branch and cut, and manage the trade-off between tighter bounds and slower LPs themselves.

Modelling patterns that decide solve time

The formulation matters more than the solver. Two models of the same problem can have the same integer solutions and very different relaxations, and the one with the tighter relaxation can solve orders of magnitude faster. These patterns cover most models.

• Yes-or-no decisions. A binary y_j for each option: open warehouse j, assign job i to machine j, include item j.
• Fixed charge. A cost F paid if production x > 0: write x ≤ U·y with y binary and add F·y to the cost. U must be a true upper bound on x, and should be the smallest one you can justify.
• Either-or constraints. Either a^Tx ≤ b or d^Tx ≤ e: add a binary z and write a^Tx ≤ b + M·z and d^Tx ≤ e + M·(1 - z). This is the big-M pattern, and M is the danger: a huge M makes the relaxation nearly useless and invites numerical trouble. Many solvers accept indicator constraints that state the logic directly; prefer them where available.
• Products of binaries. w = x·y for binaries becomes w ≤ x, w ≤ y, w ≥ x + y - 1, w ≥ 0.
• Disaggregate. For facility location, x_ij ≤ y_j per customer is a much tighter relaxation than Σ_i x_ij ≤ n·y_j.
• Symmetry. Ten identical trucks give 10! equivalent solutions; require truck k to be used only if truck k-1 is.

Using a real solver

For real work, use a solver. SciPy's milp wraps HiGHS and takes arrays: an objective, a LinearConstraint, Bounds and an integrality array (0 continuous, 1 integer). It minimises, so negate to maximise.

import numpy as np
from scipy.optimize import milp, LinearConstraint, Bounds

c = -np.array([5, 4])                    # milp minimises, so negate to maximise
A = np.array([[6, 4], [1, 2]])
cons = LinearConstraint(A, lb=-np.inf, ub=[24, 6])
res = milp(c, constraints=cons, integrality=np.ones(2),
           bounds=Bounds(0, np.inf),
           options={"time_limit": 60, "mip_rel_gap": 1e-4})

print(res.status, res.message)           # 0 means optimal
print(res.x, -res.fun)                   # x = (4, 0), objective 20
print(res.mip_node_count, res.mip_gap, res.mip_dual_bound)

Read three things from every result. status says whether the answer is proven optimal (0), stopped at a limit (1), infeasible (2) or unbounded (3). mip_gap and mip_dual_bound say how far from proven optimal a limited run stopped. mip_node_count says how hard the search was; a large count on a small model signals a weak formulation. The default mip_rel_gap of 10^-4 stops at a proven 0.01% gap; raise it when a 1% answer now beats an exact one later. For larger models, modelling layers such as PuLP, Pyomo or OR-Tools let you name constraints and switch between open-source and commercial solvers.

Failure modes

• Huge big-M values. The relaxation becomes weak, the search explodes, and with M near 10⁶ times the other coefficients a solver's integrality tolerance can accept a binary of 0.000001 as zero while it still switches on a constraint. Compute the smallest valid M per constraint.
• Reading 0.9999999 as 0. Integer variables come back within a feasibility tolerance of an integer; round them before using them as indices or counts.
• Infeasible with no explanation. Ask the solver for an irreducible infeasible subsystem where it supports one, or relax constraints with penalised slack variables to find which ones conflict.
• Optimal for the wrong problem. A missing constraint lets the solver exploit the hole; check solutions with independent code.
• Unbounded time. Always set a time limit and a gap, and decide in advance what to do with an unproven answer.

Trade-offs

Integer programming gives a provable bound on distance from optimal and absorbs side constraints as extra rows. It pays with unpredictable run times: a model that solves in seconds can stall after a small data change. Flows, matchings and dynamic programs are faster when the problem fits them exactly. Constraint programming solvers such as CP-SAT often do better on scheduling with weak LP relaxations. Metaheuristics scale further but give no bound. A common pattern combines them: a heuristic supplies a starting solution and the integer program improves it under a time limit.

What to do next

1. Write your problem as an LP first and solve the relaxation; its value bounds what any integer solution can achieve.
2. Check whether the constraints form a network or bipartite structure, and if so solve the LP or a flow algorithm directly.
3. Model with binaries for decisions, keep every big-M as small as you can justify, and use indicator constraints when your solver supports them.
4. Solve with SciPy milp or a modelling layer, with a time limit and a gap, and log status, gap, bound and node count.
5. If the node count is large, tighten the formulation: disaggregate constraints, break symmetry, add valid bounds.
6. Validate each solution with independent code that checks the business rules, and round integer variables before use.
7. Compare against a simple heuristic on real data so you know what the solver's extra minutes are buying.