Mixed-integer programming (branch and bound)

Deep dive

How it works

Parts:

• presolve — shrinks the model before the tree starts
• LP solver — dual simplex, warm-started from the parent basis
• open-node list — nodes as bound changes plus the parent bound, ordered by the node-selection rule
• branching rule — picks the fractional variable to split
• incumbent and gap — best integer solution so far and best open bound
• cuts and heuristics — tighten the root and find incumbents early

1. Root: solve the LP relaxation with simplex after presolve. Infeasible means the MIP is infeasible. An integral solution means done.
2. Pick a node from the open list (see node selection below).
3. Branch: choose a fractional integer variable x_j = f; create children x_j ≤ ⌊f⌋ and x_j ≥ ⌈f⌉. A child stores only the bound changes along its path, not a model copy.
4. Solve the child LP by dual simplex from the parent's basis: a bound change keeps the basis dual feasible.
5. Prune when the LP is infeasible, when its value cannot beat the incumbent (integer objectives may round the bound), or when the solution is integral (then update the incumbent).
6. Stop when the open list is empty or |best bound − incumbent| / |incumbent| is below the MIP gap.

root LP ─► open list ─pick─► node LP ─infeasible / bound ≥ incumbent─► prune
              ▲                 │ fractional x_j = f
              └── 2 children ◄──┘ (x_j ≤ ⌊f⌋, x_j ≥ ⌈f⌉)        integral ─► new incumbent

Branching rules (which fractional variable):

• Most fractional: f closest to 0.5. Cheap, but often little better than random (Achterberg et al.).
• Pseudocost: per variable and direction, keep the average objective degradation per unit of fractionality seen in earlier branchings; pick the largest predicted product of the two children's gains (Bénichou et al.).
• Strong branching: tentatively solve both children for each candidate with a capped number of dual simplex iterations; the best score wins. Small trees, costly per node.
• Reliability branching: strong branching until a variable's pseudocosts have been observed a few times, then pseudocosts (Achterberg, Koch, Martin).

Node selection:

• Best-bound: always the open node with the best LP bound. Minimises explored nodes for a given branching rule, but the open list grows large and dives are rare, so incumbents arrive late.
• Depth-first (dive): follow a child down. Memory-light, warm starts are cheapest (child of the node just solved), and it finds incumbents early, but it can wander in a poor subtree.
• Hybrids: dive for a while, then jump to the best-bound node; best-estimate rules use pseudocosts to guess the best integer solution under a node (Linderoth and Savelsbergh).

Incumbents and primal heuristics: besides integral LP solutions, solvers run rounding, diving, the feasibility pump (Fischetti, Glover, Lodi) and neighbourhood searches such as RINS (Danna, Rothberg, Le Pape) from LP solutions. An early incumbent prunes more.

Cuts: at the root (and sparingly in the tree) tighten the relaxation, see cutting planes. Column generation instead adds columns at nodes (branch and price).

Code

Best-bound branch and bound for the 0/1 knapsack (max Σv·x, Σw·x ≤ cap). The node LP is the fractional greedy over free items (exact LP solution), a node holds only its fixings, and the branching variable is the single fractional item. Compiled; checked against brute force on 500 random instances (3-14 items), and a gap of 1 never lost more than 1.

use std::collections::BinaryHeap;

/// 0/1 knapsack, integer weights/values, items already sorted by value/weight descending.
/// LP relaxation = fractional greedy; node = bound changes only (fixings).
pub struct Item { pub w: i64, pub v: i64 }

#[derive(Clone, Copy, PartialEq)]
enum Fix { Free, In, Out }

struct Node { bound: i64, fix: Vec<Fix> } // best-bound order
impl PartialEq for Node { fn eq(&self, o: &Self) -> bool { self.bound == o.bound } }
impl Eq for Node {}
impl PartialOrd for Node { fn partial_cmp(&self, o: &Self) -> Option<std::cmp::Ordering> { Some(self.cmp(o)) } }
impl Ord for Node { fn cmp(&self, o: &Self) -> std::cmp::Ordering { self.bound.cmp(&o.bound) } }

/// Solve the node LP. Returns None if infeasible, else (floor of LP value, fractional item).
fn lp(items: &[Item], cap: i64, fix: &[Fix]) -> Option<(i64, Option<usize>)> {
    let (mut w, mut v) = (0, 0);
    for (i, f) in fix.iter().enumerate() {
        if *f == Fix::In { w += items[i].w; v += items[i].v; }
    }
    if w > cap { return None; }
    for (i, it) in items.iter().enumerate() {
        if fix[i] != Fix::Free { continue; }
        if w + it.w <= cap { w += it.w; v += it.v; continue; }
        let frac = (it.v * (cap - w)) as f64 / it.w as f64; // the one fractional variable
        return Some((v + frac.floor() as i64, Some(i)));
    }
    Some((v, None)) // LP solution is integral
}

/// Returns (best value, nodes explored). gap_abs: stop when best_bound - incumbent <= gap_abs.
pub fn branch_and_bound(items: &[Item], cap: i64, gap_abs: i64) -> (i64, usize) {
    let mut inc = 0; // incumbent (all-zero is feasible)
    let mut nodes = 0;
    let mut heap = BinaryHeap::new();
    let root = vec![Fix::Free; items.len()];
    if let Some((b, _)) = lp(items, cap, &root) { heap.push(Node { bound: b, fix: root }); }
    while let Some(node) = heap.pop() {
        if node.bound - inc <= gap_abs { break; } // heap top = global dual bound
        nodes += 1;
        let Some((b, frac)) = lp(items, cap, &node.fix) else { continue };
        if b <= inc { continue; }
        let Some(j) = frac else { inc = b; continue }; // integral: new incumbent
        for f in [Fix::In, Fix::Out] {
            let mut fix = node.fix.clone();
            fix[j] = f;
            if let Some((cb, _)) = lp(items, cap, &fix) {
                if cb > inc { heap.push(Node { bound: cb, fix }); }
            }
        }
    }
    (inc, nodes)
}

Complexity and constants

• Worst case exponential in the number of integer variables (2ⁿ leaves). In practice node count depends on relaxation strength far more than on n.
• Per node: one dual simplex re-solve (a few to tens of pivots, each touching the LU factors), plus branching scores, heuristics and bound propagation. Time is dominated by LP solves: the LP solve dominates.
• Measured by me (code above, release build, random instances, 3-14 items): average 6.7 nodes per instance. This only checks the logic; it says nothing about real MIPs.
• Memory: open nodes times bytes per node. Storing the full vector of fixings (as in the code) costs O(n) per node; solvers store only the bound changes versus the parent (a few words per node).

Efficiency

• Cost: one node is one warm-started dual simplex solve plus propagation. Bound: LP solve time (memory bandwidth and sparse index chasing in the LU solves), and for best-bound search, the memory of the open list.
• Levers:
  • Warm start — parent basis kept, so a child needs few pivots instead of a cold simplex solve — always.
  • Store bound changes only — a node is the parent pointer plus a handful of (var, bound) pairs; the basis is shared or recomputed — large trees.
  • Presolve and a tight formulation — fewer variables and a stronger relaxation shrink the tree itself, the biggest lever — before anything else (presolve, cuts).
  • Reliability branching — strong-branching cost only where pseudocosts are unknown — avoids the weakness of pure pseudocosts early in the tree — general MIPs.
  • Gap tolerance and early incumbent — stopping at a 1 % gap or with a heuristic solution can save most of the proof time — when a proof of exact optimality is not needed.
• Naive → standard → state of the art:
  • Depth-first, most-fractional branching, cold LP at every node · tree of up to millions of nodes, each a full solve.
  • Warm-started dual simplex, pseudocost branching, best-bound with dives, root cuts · standard in open-source MIP codes.
  • Reliability branching plus presolve, cuts, conflict analysis, restarts, primal heuristics and parallel tree search · what SCIP-class and commercial solvers do.

Implementation notes

• Tolerances: an LP value counts as integral within order of 1e-6 of an integer; set both an absolute and a relative gap. Check integral solutions against original constraints, not only the relaxation.
• Bound strengthening: with an integer-valued objective, prune at ⌈bound⌉ (max: ⌊bound⌋ as in the code).
• Memory vs search order: best-bound keeps every open node; depth-first keeps one path plus siblings. Solvers switch to depth-first or limited-discrepancy diving when the open list passes a memory limit.
• Parallel tree search: nodes are independent given the incumbent. Workers share the incumbent and the open list; the root is the bottleneck, so racing different settings on the root is common. Koch et al. describe scaling to very many cores in a distributed setting.
• Determinism: thread timing changes which node is solved first and so which incumbent appears first. Deterministic modes synchronise in batches (nodes processed in rounds, results merged in fixed order) at some loss of speed; budget heuristics in operation counts rather than wall-clock (see presolve).

Pitfalls

• Big-M constraints with large M: the relaxation is nearly useless and fractional solutions satisfy "x = 0.001 turns on the whole facility". Tighten M per constraint.
• Symmetry (identical vehicles, nurses, bins): thousands of equivalent subtrees; add symmetry-breaking constraints or orbital rules.
• Judging by the incumbent alone: a gap of 5 % at timeout means the optimum may be 5 % better, not that the answer is 5 % wrong.

Variants and alternatives

In production

SCIP, HiGHS, CBC, OR-Tools (via solvers and CP-SAT), and commercial Gurobi, CPLEX and Xpress run branch and cut. Typical applications are the goals above. Energy: day-ahead electricity market clearing at MISO (Carlson et al.). Logistics: vehicle routing. Workforce: nurse rostering. Benchmarking: MIPLIB 2017 collects instances from many application areas (Gleixner et al.).

Implementations

Sources

• A. H. Land, A. G. Doig, An automatic method of solving discrete programming problems, Econometrica 28(3), 1960. doi
• R. J. Dakin, A tree-search algorithm for mixed integer programming problems, The Computer Journal 8(3), 1965, pp. 250-255. doi
• M. Bénichou, J. M. Gauthier, P. Girodet, G. Hentges, G. Ribière, O. Vincent, Experiments in mixed-integer linear programming, Mathematical Programming 1, 1971. doi
• Jeff T. Linderoth, Martin W. P. Savelsbergh, A computational study of search strategies for mixed integer programming, INFORMS Journal on Computing 11(2), 1999. doi
• Tobias Achterberg, Thorsten Koch, Alexander Martin, Branching rules revisited, Operations Research Letters 33(1), 2005. doi
• Matteo Fischetti, Fred Glover, Andrea Lodi, The feasibility pump, Mathematical Programming 104, 2005. doi
• Emilie Danna, Edward Rothberg, Claude Le Pape, Exploring relaxation induced neighborhoods to improve MIP solutions, Mathematical Programming 102, 2005. doi
• Thorsten Koch, Ted Ralphs, Yuji Shinano, Could we use a million cores to solve an integer program?, Mathematical Methods of Operations Research 76, 2012. doi
• Ambros Gleixner et al., MIPLIB 2017: data-driven compilation of the 6th mixed-integer programming library, Math. Program. Comput. 13, 2021. doi
• Ksenia Bestuzheva et al., The SCIP Optimization Suite 8.0, 2021. arXiv