Benders Decomposition

Benders Decomposition

You are an expert in exact decomposition methods for mixed-integer optimization, specifically Benders decomposition in its classic (iterative master/subproblem loop) and modern (branch-and-Benders-cut via lazy-constraint callbacks) forms, including the L-shaped method for two-stage stochastic programs. This skill covers the master–subproblem split, the derivation of optimality and feasibility cuts from LP duality, complete gurobipy implementations of both execution modes, and the acceleration techniques (Pareto-optimal cuts, stabilization, cut aggregation) that decide whether Benders converges in 20 iterations or 2,000. Use the framework below to verify the problem has the right structure, derive the cuts on paper, implement against the reusable loop or the callback pattern, and validate against the monolithic model.

Initial Assessment

Before decomposing anything, establish the following. Each answer changes the design.

• Identify the complicating variables. Which variables, once fixed, leave an easy remaining problem? Benders needs a clean split: integer/design variables y in the master, continuous recourse variables x in the subproblem. If no such split exists, Benders is the wrong tool.
• Subproblem class. Is the subproblem an LP for every fixed y? Classic Benders cuts come from LP duals. If the subproblem keeps integer variables, you need logic-based Benders or integer L-shaped cuts (see Advanced Techniques) — a different, weaker machinery.
• Subproblem separability. Does the subproblem split into independent blocks (per scenario, per customer, per period)? Separability is the main source of speedup and enables multi-cut formulations and parallel subproblem solves.
• Feasibility structure. Can the subproblem be infeasible for some master solutions? If yes, you need feasibility cuts (dual extreme rays / Farkas certificates). Check first whether a small master-side constraint (e.g., total capacity ≥ total demand) or penalized slack variables can give relatively complete recourse and remove feasibility cuts entirely.
• Why decompose at all? Estimate the extensive (monolithic) model size: variables = |y| + |x|·(blocks), constraints likewise. Modern solvers handle millions of nonzeros; decompose only when the monolith is too large, too slow, or the subproblem has special structure (closed-form duals, network structure) the solver cannot exploit.
• Execution mode. Classic loop (sequence of master MIPs) for prototyping, analysis, and cheap masters; branch-and-Benders-cut (one search tree, lazy cuts) for production runs where re-solving the master MIP from scratch each iteration is wasteful.
• Bound on the recourse term. What is a valid lower bound for eta (minimization)? Without one the first master solve is unbounded. Nonnegative recourse costs give eta ≥ 0; per-block bounds (e.g., cheapest assignment per customer) are tighter and free.
• Solver and license. gurobipy available? Callbacks, FarkasDual/UnbdRay, and lazy constraints are needed. With open-source solvers, SCIP supports Benders plugins; HiGHS supports only the classic loop.
• Scale and budget. Number of scenarios/blocks, master integer variables, target gap, wall-clock budget. These set single-cut vs multi-cut, stabilization needs, and whether subproblems must be parallelized.
• Validation baseline. Build the monolithic model on small instances first. Every Benders implementation must reproduce its optimum exactly before you trust it at scale.

Decomposition Anatomy

Benders decomposition (Benders 1962, "Partitioning procedures for solving mixed-variables programming problems") targets problems of the form

$$ \min_{y,\,x} \; f^\top y + q^\top x \quad \text{s.t.} \quad W x \ge h - T y, \qquad x \ge 0, \qquad y \in Y, $$

where y are the complicating variables (typically integer: open/close, build/buy, capacity levels) and, for fixed y, the remaining problem in x is an LP. Project x out:

$$ \min_{y \in Y} \; f^\top y + z(y), \qquad z(y) = \min \{\, q^\top x \;:\; W x \ge h - T y,\; x \ge 0 \,\}. $$

Dualize the inner LP. Its dual feasible region

$$ U = \{\, u \ge 0 \;:\; W^\top u \le q \,\} $$

does not depend on y — only the dual objective does:

$$ z(y) = \max \{\, u^\top (h - T y) \;:\; u \in U \,\}. $$

This is the single fact the whole method rests on. U is a fixed polyhedron with finitely many extreme points $u1, \dots, uP$ and extreme rays $r1, \dots, rR$, so for every y:

• if the subproblem is feasible, $z(y) = \max{p} \; up^\top (h - T y)$ (the max is attained at an extreme point);
• the subproblem is feasible iff $rk^\top (h - T y) \le 0$ for every extreme ray $rk$ — otherwise the dual is unbounded along some ray, which is exactly a Farkas certificate of primal infeasibility.

Substituting gives a master problem equivalent to the original:

$$ \min{y \in Y,\ \eta} \; f^\top y + \eta \quad \text{s.t.} \quad \underbrace{\eta \ge up^\top (h - T y)}{\text{optimality cuts, } p = 1..P}, \qquad \underbrace{rk^\top (h - T y) \le 0}_{\text{feasibility cuts, } k = 1..R}. $$

Nobody enumerates P and R. The algorithm keeps a relaxed master with a small subset of cuts and alternates:

1. Solve the relaxed master → candidate $(\hat{y}, \hat{\eta})$ and a lower bound $f^\top \hat{y} + \hat{\eta}$ (minimization).
2. Solve the dual subproblem at $\hat{y}$. Unbounded → add the feasibility cut from the ray. Optimal with value $z(\hat{y})$ → an upper bound $f^\top \hat{y} + z(\hat{y})$, and if $z(\hat{y}) > \hat{\eta} + \varepsilon$, add the optimality cut from the optimal extreme point.
3. Stop when $z(\hat{y}) \le \hat{\eta} + \varepsilon$: then UB ≤ LB, so $\hat{y}$ is optimal.

The lower bound is monotone nondecreasing (the master only gains constraints); the upper bound is not monotone — always track the incumbent best. Convergence is finite because each iteration produces an extreme point or ray not yet in the master, and there are finitely many.

Two execution modes

| Mode | How it runs | Use when | |---|---|---| | Classic iterative loop | Solve master MIP to optimality, add cuts, repeat | Prototyping; cheap or LP master; root-node cut warm-up; teaching the cut logic | | Branch-and-Benders-cut | One B&B tree on the master; separate cuts at each integer candidate via lazy-constraint callback | Production default; master MIP is expensive; avoids re-proving the same branching work every iteration |

The classic loop wastes effort: iteration k re-solves a master MIP that differs from iteration k−1 by one row. Branch-and-Benders-cut (also called one-tree Benders; see Fortz & Poss 2009 and the survey by Rahmaniani, Crainic, Gendreau & Rei 2017, "The Benders decomposition algorithm: A literature review") keeps a single tree and rejects integer candidates with lazy cuts. Worked Example 1 implements it.

When Benders pays off — and when it does not

Use Benders when at least one of these holds:

• The subproblem separates into many independent LPs (scenarios in stochastic programs, customers in facility location, commodities in network design). The extensive form is huge; the pieces are tiny.
• The subproblem has special structure: closed-form dual solutions (UFL — Worked Example 1), network flow (solvable by a combinatorial algorithm), or a structure destroyed by mixing with y.
• Memory: the monolithic model does not fit, but master + one block at a time does.

Avoid Benders when the subproblem keeps integer variables (classic duals do not exist), when T is dense so every cut is dense and the master degrades, or when the monolith solves in minutes anyway — a modern MIP solver with a strong formulation beats a naive decomposition embarrassingly often. Build the monolith first; it is both the baseline and the burden of proof.

The L-shaped method

A two-stage stochastic program with finite scenarios $s = 1..S$, probabilities $p_s$,

$$ \min{y \in Y} \; f^\top y + \sum{s} ps\, Qs(y), \qquad Qs(y) = \min \{\, qs^\top x \;:\; W x \ge hs - Ts y,\; x \ge 0 \,\}, $$

is exactly the Benders structure with a block-diagonal subproblem — one block per scenario. Benders applied to it is the L-shaped method (Van Slyke & Wets 1969). Two cut layouts:

• Single-cut: one variable theta approximates the whole expectation; each iteration adds one cut aggregating the probability-weighted duals of all scenarios.
• Multi-cut: one theta_s per scenario, one cut per violated scenario (Birge & Louveaux 1988). More master rows, far fewer iterations; the right default when S is moderate (≤ a few thousand).

Feasibility cuts disappear under relatively complete recourse — the subproblem is feasible for every master-feasible y. Design for it: add penalized shortfall variables to the second stage (Worked Example 2) instead of letting scenarios go infeasible.

Generic Benders Framework

The reusable artifact: a classic Benders loop for min f@y + q@x s.t. Wx >= h - Ty, x >= 0, y binary, with the dual subproblem built once (its feasible region never changes — only the objective moves with y), feasibility cuts from extreme rays, and the standard termination test. The worked examples follow the same anatomy; every block is self-contained and runs as-is.

BENDERS DECOMPOSITION (classic loop, minimization)
--------------------------------------------------
input: f, q, W, T, h; optional master-only constraints D y >= d
state: relaxed master over (y, eta); dual subproblem over u (built once)

1. master: min f@y + eta, s.t. D y >= d, eta >= eta_lb
2. repeat:
3.     solve master -> (y_hat, eta_hat); LB = 0
5.     if unbounded along ray r:
6.         add feasibility cut  (r@T) y >= r@h          # cuts y_hat off
7.     else, with optimal point u*, value z_hat:
8.         UB = u*@h
11. return incumbent; LB is monotone, UB is not -> track the best

"""Reusable classic Benders loop: binary master, LP subproblem via its dual."""
from __future__ import annotations

import math
from dataclasses import dataclass

import gurobipy as gp
import numpy as np
from gurobipy import GRB

@dataclass
class BendersData:
    """min f@y + q@x  s.t.  W x >= h - T y, x >= 0, y binary, D y >= d (master-only)."""

    f: np.ndarray            # (n_y,) master costs
    q: np.ndarray            # (n_x,) subproblem costs
    W: np.ndarray            # (m, n_x) subproblem matrix
    T: np.ndarray            # (m, n_y) linking matrix
    h: np.ndarray            # (m,) linking rhs
    D: np.ndarray | None = None   # optional master-only constraints D y >= d
    d: np.ndarray | None = None

class DualSubproblem:
    """Dual subproblem max u@(h - T y) s.t. W'u = 0 — built once.

    The feasible region is independent of y; only the objective changes.
    """

    def __init__(self, data: BendersData) -> None:
        self.data = data
        self.model = gp.Model("benders-dual-subproblem")
        self.model.Params.OutputFlag = 0
        self.model.Params.InfUnbdInfo = 1     # expose UnbdRay when unbounded
        self.model.Params.DualReductions = 0  # report UNBOUNDED, never INF_OR_UNBD
        self.model.Params.Method = 0          # primal simplex: returns extreme rays
        self.u = self.model.addMVar(data.W.shape[0], lb=0.0, name="u")
        self.model.addConstr(data.W.T @ self.u  tuple[str, float, np.ndarray]:
        """Return ('point', z(y), u*) or ('ray', inf, r) at the given master solution."""
        rhs = self.data.h - self.data.T @ y_val
        self.model.setObjective(rhs @ self.u, GRB.MAXIMIZE)
        self.model.optimize()
        if self.model.Status == GRB.OPTIMAL:
            return "point", self.model.ObjVal, self.u.X
        if self.model.Status == GRB.UNBOUNDED:
            ray = np.array(self.model.getAttr(GRB.Attr.UnbdRay, self.model.getVars()))
            return "ray", math.inf, ray
        raise RuntimeError(f"dual subproblem status {self.model.Status}: "
                           "dual infeasible means the primal recourse is unbounded")

def benders_solve(data: BendersData, eta_lb: float = 0.0, tol: float = 1e-6,
                  max_iters: int = 500) -> dict:
    """Classic Benders loop. eta_lb must be a valid lower bound on z(y)."""
    master = gp.Model("benders-master")
    master.Params.OutputFlag = 0
    y = master.addMVar(data.f.size, vtype=GRB.BINARY, name="y")
    eta = master.addVar(lb=eta_lb, name="eta")
    master.setObjective(data.f @ y + eta, GRB.MINIMIZE)
    if data.D is not None:
        master.addConstr(data.D @ y >= data.d, name="master_side")
    sub = DualSubproblem(data)

    best_ub, best_y, log = math.inf, None, []
    for it in range(1, max_iters + 1):
        master.optimize()
        if master.Status != GRB.OPTIMAL:
            raise RuntimeError(f"master status {master.Status}: feasibility cuts "
                               "may have emptied Y — check the original model")
        lb = master.ObjVal
        y_val = np.rint(y.X)              # clean integrality noise before cutting
        kind, z_hat, vec = sub.solve(y_val)
        if kind == "ray":                 # feasibility cut: (r@T) y >= r@h
            master.addConstr((vec @ data.T) @ y >= float(vec @ data.h),
                             name=f"feas_cut_{it}")
            log.append((it, lb, best_ub, "feasibility"))
            continue
        ub = float(data.f @ y_val) + z_hat
        if ub = float(vec @ data.h),
                         name=f"opt_cut_{it}")
    return {"y": best_y, "objective": best_ub, "lower_bound": lb, "log": log}

# --- tiny instance: open plants (capacity 8 each) to serve demand 10 --------
f = np.array([10.0, 14.0])                          # plant opening costs
q = np.array([1.0, 2.0])                            # unit production costs
W = np.array([[1.0, 1.0], [-1.0, 0.0], [0.0, -1.0]])  # demand row, capacity rows
h = np.array([10.0, 0.0, 0.0])
T = np.array([[0.0, 0.0], [8.0, 0.0], [0.0, 8.0]])  # -x_i >= -8 y_i
res = benders_solve(BendersData(f, q, W, T, h))
print(res["y"], round(res["objective"], 4), [step[3] for step in res["log"]])
# Expected: [1 1] 36.0 ['feasibility', 'optimality', 'optimality'] -- the first
# master tries y=(0,0); the ray cut 8*y0 + 8*y1 >= 10 forces both plants open
# (one plant covers only 8 = 1`, so the subproblem is always feasible — an instance of the general rule *move structural feasibility into the master when you can*.

```python
"""UFLP Benders separation in closed form: dual values and one cut per customer."""
from __future__ import annotations

import numpy as np

def ufl_cut(c_j: np.ndarray, y: np.ndarray) -> tuple[float, np.ndarray]:
    """Closed-form dual for one customer at binary y (>= 1 facility open).

    Returns (v, w) defining the optimality cut  eta_j >= v - w @ y,
    with v = cost of the cheapest open facility and w_i = max(0, v - c_ij).
    """
    v = float(c_j[y > 0.5].min())
    w = np.maximum(0.0, v - c_j)
    return v, w

# --- tiny instance: 3 facilities x 4 customers, facilities 0 and 2 open -----
c = np.array([[2.0, 7.0, 5.0, 4.0],
              [5.0, 3.0, 6.0, 8.0],
              [6.0, 4.0, 2.0, 3.0]])   # c[i, j]
y_hat = np.array([1, 0, 1])
for j in range(c.shape[1]):
    v, w = ufl_cut(c[:, j], y_hat)
    print(j, v, w)
# Expected: v = 2, 4, 2, 3 for customers 0..3; all w vectors are zero except
# customer 1, where w = [0, 1, 0]: its cut reads eta_1 >= 4 - 1*y_1, i.e. the
# bound drops to 3 if the (currently closed, cheaper) facility 1 opens.

The production implementation is branch-and-Benders-cut: one master tree over (y, eta), with the cuts above added lazily whenever Gurobi finds an integer candidate (MIPSOL). Params.LazyConstraints = 1 is mandatory — it tells the solver that not all constraints are present, disabling reductions that would otherwise cut off solutions only the callback knows are infeasible.

"""UFLP by branch-and-Benders-cut: lazy per

…

## 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.