Logo
Published on

Aspects of OR - Integer Programming (Part 2)

Authors
  • avatar
    Name
    Till Heller
    Twitter

In Part 1 we saw that integer programs are much harder than LPs, and that naive rounding of the LP relaxation does not work. We need a systematic way to search the space of integer solutions while avoiding exhaustive enumeration. The answer is branch and bound — an algorithm that combines divide-and-conquer with LP relaxation bounds to explore only a fraction of the solution space.

Branch and bound (B&B) was introduced by Land and Doig in 1960 and remains the backbone of every modern IP solver. In this post, we explain the algorithm, walk through a complete example that demonstrates all three types of pruning, and discuss the design choices that make B&B efficient.

The idea

Branch and bound works by building a search tree. At the root, we solve the LP relaxation. If the solution happens to be integer, we are done. If not, we pick a fractional variable and branch: create two child nodes that force this variable to be either rounded down or rounded up. We then solve the LP relaxation at each child.

The key insight is that we do not need to explore every node. The LP relaxation at each node provides a bound on the best IP solution achievable in that subtree. If this bound is worse than the best integer solution found so far (the incumbent), we can prune the entire subtree — no integer solution in it can improve the incumbent. This is what makes B&B practical: good bounds let us discard large parts of the search space without ever looking at them.

The algorithm

Branch and bound for a maximization IP works as follows:

Initialize: Set the incumbent value to -\infty. Create a list of active nodes containing only the root.

Step 1 — Select. Choose an active node from the list.

Step 2 — Relax. Solve the LP relaxation at this node.

Step 3 — Prune. Check if the node can be pruned by one of three criteria:

  • Pruning by bound: The LP optimal value \leq incumbent value. No integer solution in this subtree can beat what we already have.
  • Pruning by infeasibility: The LP relaxation is infeasible. The branching constraints made the subproblem impossible.
  • Pruning by integrality: The LP solution is integer. Update the incumbent if this solution is better.

Step 4 — Branch. If the node is not pruned, select a fractional variable xjx_j with LP value xj=fx_j^* = f (where fZf \notin \mathbb{Z}). Create two child nodes:

  • Left child: add constraint xjfx_j \leq \lfloor f \rfloor.
  • Right child: add constraint xjfx_j \geq \lceil f \rceil.

Add both children to the active node list.

Step 5 — Terminate. If the active node list is empty, the incumbent is optimal. Otherwise, go to Step 1.

Worked example: knapsack

Let us trace B&B on the knapsack instance from Part 1. This small example demonstrates all three pruning types — integrality, bound, and infeasibility.

max  16x1+22x2+12x3+8x4s.t.6x1+8x2+5x3+4x413,xj{0,1}.\max \; 16 x_1 + 22 x_2 + 12 x_3 + 8 x_4 \quad \text{s.t.} \quad 6 x_1 + 8 x_2 + 5 x_3 + 4 x_4 \leq 13, \quad x_j \in \{0, 1\}.

ItemWeightValueValue/Weight
16162.67
28222.75
35122.40
4482.00

Capacity W=13W = 13.

Node 0 (root)

Solve the LP relaxation (allow 0xj10 \leq x_j \leq 1). Sorting by value-to-weight ratio: item 2 (2.75), item 1 (2.67), item 3 (2.40), item 4 (2.00). Greedily: take x2=1x_2 = 1 (weight 8), then x1=5/60.833x_1 = 5/6 \approx 0.833 (weight 5, filling the remaining capacity).

LP solution: x1=5/6,  x2=1,  x3=0,  x4=0x_1 = 5/6, \; x_2 = 1, \; x_3 = 0, \; x_4 = 0. LP value =5/6×16+22=35.33= 5/6 \times 16 + 22 = 35.33.

Fractional variable: x1=5/6x_1 = 5/6. Branch on x1x_1.

Incumbent: none. Upper bound: 35.33.

Node 1: x1=0x_1 = 0

Subproblem: max  22x2+12x3+8x4\max \; 22 x_2 + 12 x_3 + 8 x_4 s.t. 8x2+5x3+4x4138 x_2 + 5 x_3 + 4 x_4 \leq 13, xj[0,1]x_j \in [0,1].

LP solution: x2=1x_2 = 1 (weight 8), x3=1x_3 = 1 (weight 5). Total weight 13, LP value =22+12=34= 22 + 12 = 34.

This solution is integer! Pruned by integrality. Update incumbent to 34.

Node 2: x1=1x_1 = 1

Subproblem: add item 1 (weight 6), leaving capacity 136=713 - 6 = 7. Maximize 16+22x2+12x3+8x416 + 22 x_2 + 12 x_3 + 8 x_4 s.t. 8x2+5x3+4x478 x_2 + 5 x_3 + 4 x_4 \leq 7, xj[0,1]x_j \in [0,1].

LP solution: x2=7/8=0.875x_2 = 7/8 = 0.875 (weight 7). LP value =16+0.875×22=35.25= 16 + 0.875 \times 22 = 35.25.

LP bound 35.25 > incumbent 34, so we cannot prune yet. Fractional variable: x2=0.875x_2 = 0.875. Branch on x2x_2.

Node 3: x1=1,x2=0x_1 = 1, x_2 = 0

Remaining capacity: 136=713 - 6 = 7. Maximize 16+12x3+8x416 + 12 x_3 + 8 x_4 s.t. 5x3+4x475 x_3 + 4 x_4 \leq 7, xj[0,1]x_j \in [0,1].

LP solution: x3=1x_3 = 1 (weight 5), x4=0.5x_4 = 0.5 (weight 2). LP value =16+12+4=32= 16 + 12 + 4 = 32.

LP bound 32 \leq incumbent 34. Pruned by bound. No integer solution in this subtree can beat 34 — we discard it entirely without exploring further.

Node 4: x1=1,x2=1x_1 = 1, x_2 = 1

Remaining capacity: 1368=1<013 - 6 - 8 = -1 < 0. Items 1 and 2 together weigh 14, exceeding the capacity of 13.

Pruned by infeasibility. This subproblem has no feasible solution.

Termination

All nodes are pruned. The optimal solution is x1=0,x2=1,x3=1x_1 = 0, x_2 = 1, x_3 = 1 with value 34.

Here is the search tree:

                  Node 0: LP=35.33
                 /              \
        x₁=0  /                \ x₁=1
             /                  \
      Node 1: LP=34          Node 2: LP=35.25
      INTEGER/            \
      (incumbent=34)    x₂=0 /              \ x₂=1
                             /                \
                      Node 3: LP=32      Node 4: INFEASIBLE
                      PRUNED BY BOUND    PRUNED

The tree had 5 nodes. Exhaustive enumeration would check 24=162^4 = 16 solutions. For larger problems, the savings from pruning are dramatic — modern solvers routinely prune 99.9%+ of the tree.

Notice what happened: we found a good incumbent early (Node 1, value 34) which then allowed us to prune Node 3 by bound. This is the engine of branch and bound — a good incumbent powers aggressive pruning.

Node selection strategies

The order in which we explore nodes significantly affects performance:

Best-first (best bound): Always explore the node with the best LP bound. This improves the global upper bound quickly, which is useful for proving optimality. However, it may take long to find good integer solutions.

Depth-first: Always explore the most recently created node. This dives deep into the tree quickly, finding integer solutions early (which helps pruning) and using less memory. The downside is that it may explore many suboptimal subtrees before reaching the optimum.

Best-first with diving: A hybrid. Mostly use best-first, but periodically dive deep from a promising node to find feasible solutions. This is close to what modern solvers do.

Variable selection (branching rules)

Which fractional variable should we branch on? This choice has a major impact on the tree size.

Most fractional: Branch on the variable closest to 0.5 — the most "undecided" variable. Simple but often not the most effective.

Strong branching: For each candidate variable, tentatively solve the LP relaxations of both children. Pick the variable that gives the biggest bound improvement. This produces small trees but is expensive per node.

Reliability branching: Use strong branching for the first few occurrences of each variable, then switch to a cheaper heuristic based on observed improvements. This is the default in most modern solvers.

Pseudocost branching: Track the historical improvement per unit change for each variable and use these "pseudocosts" to estimate branching quality. Cheap per node once initialized.

Bounding and the integrality gap

The effectiveness of B&B depends directly on the quality of the LP bounds. The tighter the LP relaxation, the more nodes can be pruned.

The integrality gap is the difference between the LP optimal value and the IP optimal value. For our knapsack: LP bound =35.33= 35.33, IP optimum =34= 34, gap =(35.3334)/343.9%= (35.33 - 34)/34 \approx 3.9\%.

In practice, solvers report the MIP gap: the relative difference between the best known integer solution (lower bound for maximization) and the best LP bound (upper bound). The solve is complete when the MIP gap reaches zero (or a user-specified tolerance).

Implementation in Python

Here is a branch-and-bound implementation for binary programs that prints every step of the search, so you can see the algorithm in action:

from scipy.optimize import linprog
import numpy as np

def branch_and_bound(c, A_ub, b_ub, n, verbose=True):
    """
    Solve max c^T x s.t. A_ub x <= b_ub, x in {0,1}^n.
    Returns (optimal_value, optimal_x).
    """
    best_val = -np.inf
    best_x = None
    node_id = 0

    # Stack of subproblems: (lower_bounds, upper_bounds, description)
    stack = [(np.zeros(n), np.ones(n), 'root')]

    while stack:
        lb, ub, desc = stack.pop()  # depth-first
        current_id = node_id
        node_id += 1

        # Solve LP relaxation (linprog minimizes, so negate c)
        bounds = list(zip(lb, ub))
        res = linprog(-c, A_ub=A_ub, b_ub=b_ub, bounds=bounds, method='highs')

        if not res.success:
            if verbose:
                print(f'  Node {current_id} ({desc}): INFEASIBLE — pruned')
            continue

        lp_val = -res.fun
        x = res.x

        # Check integrality
        is_integer = np.all(np.abs(x - np.round(x)) < 1e-6)

        if lp_val <= best_val:
            if verbose:
                xstr = ', '.join(f'x{j+1}={x[j]:.3f}' for j in range(n))
                print(f'  Node {current_id} ({desc}): LP={lp_val:.2f}'
                      f' [{xstr}] — PRUNED BY BOUND (incumbent={best_val:.0f})')
            continue

        if is_integer:
            best_val = lp_val
            best_x = np.round(x)
            if verbose:
                xstr = ', '.join(f'x{j+1}={x[j]:.3f}' for j in range(n))
                print(f'  Node {current_id} ({desc}): LP={lp_val:.2f}'
                      f' [{xstr}] — INTEGER, incumbent={best_val:.0f}')
            continue

        # Branch on most fractional variable
        frac_vals = np.minimum(x - np.floor(x), np.ceil(x) - x)
        frac_vals[np.abs(x - np.round(x)) < 1e-6] = 0  # exclude integer vars
        branch_var = np.argmax(frac_vals)

        if verbose:
            xstr = ', '.join(f'x{j+1}={x[j]:.3f}' for j in range(n))
            print(f'  Node {current_id} ({desc}): LP={lp_val:.2f}'
                  f' [{xstr}] — branch on x{branch_var+1}')

        # Right child first (so left is popped first = depth-first left)
        new_lb = lb.copy()
        new_lb[branch_var] = np.ceil(x[branch_var])
        stack.append((new_lb, ub.copy(), f'x{branch_var+1}=1'))

        new_ub = ub.copy()
        new_ub[branch_var] = np.floor(x[branch_var])
        stack.append((lb.copy(), new_ub, f'x{branch_var+1}=0'))

    return best_val, best_x

Applying this to our knapsack:

c = np.array([16.0, 22.0, 12.0, 8.0])
A_ub = np.array([[6.0, 8.0, 5.0, 4.0]])
b_ub = np.array([13.0])

val, x = branch_and_bound(c, A_ub, b_ub, n=4)
print(f'\nOptimal value: {val:.0f}')
print(f'Solution: {x}')

Output:

  Node 0 (root): LP=35.33 [x1=0.833, x2=1.000, x3=0.000, x4=0.000] — branch on x1
  Node 1 (x1=0): LP=34.00 [x1=0.000, x2=1.000, x3=1.000, x4=0.000]INTEGER, incumbent=34
  Node 2 (x1=1): LP=35.25 [x1=1.000, x2=0.875, x3=0.000, x4=0.000] — branch on x2
  Node 3 (x2=0): LP=32.00 [x1=1.000, x2=0.000, x3=1.000, x4=0.500]PRUNED BY BOUND (incumbent=34)
  Node 4 (x2=1): INFEASIBLE — pruned

Optimal value: 34
Solution: [0. 1. 1. 0.]

This matches our manual trace exactly. The implementation is purely educational — production solvers add warm-starting, sophisticated branching rules, cutting planes (Part 3), heuristics, and parallelism.

Looking ahead

Branch and bound alone can solve many IPs, but the LP relaxation bounds are sometimes too loose, leading to enormous search trees. In Part 3, we introduce cutting planes — valid inequalities that tighten the LP relaxation by cutting off fractional solutions without removing any integer point. In Part 4, we combine both into branch-and-cut, the state of the art for solving IPs.