hajibabaie 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
yin the master, continuous recourse variablesxin 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 giveeta ≥ 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 $u_1, \dots, u_P$ and extreme rays $r_1, \dots, r_R$, so for every y:
- if the subproblem is feasible, $z(y) = \max_{p} ; u_p^\top (h - T y)$ (the max is attained at an extreme point);
- the subproblem is feasible iff $r_k^\top (h - T y) \le 0$ for every extreme ray $r_k$ — 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 u_p^\top (h - T y)}{\text{optimality cuts, } p = 1..P}, \qquad \underbrace{r_k^\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:
- Solve the relaxed master → candidate $(\hat{y}, \hat{\eta})$ and a lower bound $f^\top \hat{y} + \hat{\eta}$ (minimization).
- 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.
- 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} p_s, Q_s(y), \qquad Q_s(y) = \min {, q_s^\top x ;:; W x \ge h_s - T_s 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
thetaapproximates the whole expectation; each iteration adds one cut aggregating the probability-weighted duals of all scenarios. - Multi-cut: one
theta_sper scenario, one cut per violated scenario (Birge & Louveaux 1988). More master rows, far fewer iterations; the right default whenSis 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 <- master objective
4. solve dual subproblem: max u@(h - T y_hat) s.t. W'u <= q, u >= 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 <- min(UB, f@y_hat + z_hat)
9. if z_hat <= eta_hat + eps: STOP (y_hat optimal)
10. add optimality cut eta + (u*@T) y >= 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 <= q, 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 <= data.q, name="dual_feasibility")
def solve(self, y_val: np.ndarray) -> 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 < best_ub:
best_ub, best_y = ub, y_val.astype(int)
log.append((it, lb, best_ub, "optimality"))
if z_hat <= eta.X + tol * max(1.0, abs(z_hat)):
break # eta matches the true recourse: optimal
master.addConstr(eta + (vec @ data.T) @ y >= 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 < 10); convergence at iteration 3 with cost 36.
Two details repay attention. First, the dual subproblem is constructed once and re-optimized with a new objective each call — simplex warm starts make repeat solves nearly free, and the construction makes the theory visible: extreme points and rays come from one fixed polyhedron. Second, eta needs a valid lower bound (0 here, because q ≥ 0); without it the first master solve is unbounded and Gurobi raises an error before the loop ever produces a cut.
When the subproblem is solved in its primal form instead (often more natural — Worked Example 2), the same information comes from Gurobi attributes: Constr.Pi gives the optimal dual point, and Constr.FarkasDual (with InfUnbdInfo=1, on an infeasible primal) gives the ray. The cut algebra is identical.
Worked Example 1: Uncapacitated Facility Location via Benders Callbacks
The uncapacitated facility location problem (UFLP): open facilities $i \in F$ at fixed cost $f_i$, serve every customer $j \in C$ from an open facility at cost $c_{ij}$:
$$ \min ; \sum_i f_i y_i + \sum_j \sum_i c_{ij} x_{ij} \quad \text{s.t.} \quad \sum_i x_{ij} = 1 ;; \forall j, \qquad x_{ij} \le y_i, \qquad x \ge 0, ; y \in {0,1}^{|F|}. $$
For fixed $\hat{y} \ne 0$ the x-LP separates by customer, and each block has a closed-form optimum: assign customer j to the cheapest open facility, $z_j(\hat{y}) = \min { c_{ij} : \hat{y}i = 1 }$. The dual of customer j's block (dual variable $v_j$ for the assignment row, $w{ij} \ge 0$ for each linking row) is
$$ \max ; v_j - \sum_i w_{ij} \hat{y}i \quad \text{s.t.} \quad v_j - w{ij} \le c_{ij} ;; \forall i, \qquad w_{ij} \ge 0, $$
with the closed-form optimal extreme point $v_j = z_j(\hat{y})$, $w_{ij} = \max(0,, v_j - c_{ij})$. The disaggregated (per-customer) optimality cut is therefore
$$ \eta_j ;\ge; z_j(\hat{y}) - \sum_i \bigl(z_j(\hat{y}) - c_{ij}\bigr)^+ , y_i , $$
valid for every y by weak duality and tight at $\hat{y}$. Read it as: the bound starts at the current best cost and drops by $(z_j - c_{ij})^+$ for each cheaper facility i that opens. No LP is ever solved — the entire separation is a sort-free min and a maximum over a cost column. This closed-form, disaggregated scheme is what makes Benders the method of choice for very large UFLP (Fischetti, Ljubić & Sinnl 2017, "Redesigning Benders decomposition for large-scale facility location"). Feasibility cuts are avoided by construction: the master keeps sum(y) >= 1, so the subproblem is always feasible — an instance of the general rule move structural feasibility into the master when you can.
"""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-customer optimality cuts in Gurobi."""
from __future__ import annotations
import gurobipy as gp
import numpy as np
from gurobipy import GRB
def generate_ufl(n_fac: int, n_cust: int, seed: int) -> tuple[np.ndarray, np.ndarray]:
"""Euclidean UFLP instance in the unit square; returns (fixed costs, c[i, j])."""
rng = np.random.default_rng(seed)
fac, cust = rng.random((n_fac, 2)), rng.random((n_cust, 2))
c = 10.0 * np.linalg.norm(fac[:, None, :] - cust[None, :, :], axis=2)
return rng.uniform(5.0, 15.0, size=n_fac), c
def solve_ufl_benders(f: np.ndarray, c: np.ndarray,
tol: float = 1e-6) -> tuple[float, np.ndarray, int]:
"""Master over (y, eta_j); closed-form cuts separated in a MIPSOL callback."""
n_fac, n_cust = c.shape
model = gp.Model("ufl-benders")
model.Params.OutputFlag = 0
model.Params.LazyConstraints = 1 # mandatory for cbLazy
y = model.addVars(n_fac, vtype=GRB.BINARY, name="y")
eta = model.addVars(n_cust, lb=c.min(axis=0).tolist(), name="eta")
model.setObjective(gp.quicksum(f[i] * y[i] for i in range(n_fac))
+ gp.quicksum(eta.values()), GRB.MINIMIZE)
model.addConstr(gp.quicksum(y.values()) >= 1, name="open_at_least_one")
n_cuts = 0
def separate(model: gp.Model, where: int) -> None:
"""Lazily add every violated per-customer cut at an integer candidate."""
nonlocal n_cuts
if where != GRB.Callback.MIPSOL:
return
y_val = np.array(model.cbGetSolution([y[i] for i in range(n_fac)]))
eta_val = model.cbGetSolution([eta[j] for j in range(n_cust)])
open_mask = y_val > 0.5 # candidate values are not exact
for j in range(n_cust):
v = float(c[open_mask, j].min())
if eta_val[j] >= v - tol:
continue # this customer's eta is correct
w = np.maximum(0.0, v - c[:, j])
model.cbLazy(eta[j] >= v - gp.quicksum(
float(w[i]) * y[i] for i in range(n_fac) if w[i] > 0.0))
n_cuts += 1
model.optimize(separate)
assert model.Status == GRB.OPTIMAL and model.SolCount > 0
y_opt = np.array([y[i].X for i in range(n_fac)]).round().astype(int)
return model.ObjVal, y_opt, n_cuts
f, c = generate_ufl(n_fac=15, n_cust=40, seed=7)
obj, y_opt, n_cuts = solve_ufl_benders(f, c)
print(round(obj, 4), int(y_opt.sum()), n_cuts > 0)
# Expected: 109.0666 4 True -- four facilities open; the objective matches the
# monolithic MIP below to full precision because both are exact.
Validation is non-negotiable: a buggy separation routine does not crash — it silently returns a wrong "optimum", because with LazyConstraints=1 the solver trusts the callback to reject every bad integer point. Reproduce the monolithic optimum on the same seeded instance before scaling up.
"""Monolithic UFLP MIP: the correctness baseline for the Benders solver."""
from __future__ import annotations
import gurobipy as gp
import numpy as np
from gurobipy import GRB
def generate_ufl(n_fac: int, n_cust: int, seed: int) -> tuple[np.ndarray, np.ndarray]:
"""Identical generator and seed => identical instance to the Benders block."""
rng = np.random.default_rng(seed)
fac, cust = rng.random((n_fac, 2)), rng.random((n_cust, 2))
c = 10.0 * np.linalg.norm(fac[:, None, :] - cust[None, :, :], axis=2)
return rng.uniform(5.0, 15.0, size=n_fac), c
def solve_ufl_direct(f: np.ndarray, c: np.ndarray) -> float:
"""Strong (disaggregated) formulation with x[i, j] <= y[i]."""
n_fac, n_cust = c.shape
model = gp.Model("ufl-direct")
model.Params.OutputFlag = 0
y = model.addVars(n_fac, vtype=GRB.BINARY, name="y")
x = model.addVars(n_fac, n_cust, lb=0.0, ub=1.0, name="x")
model.setObjective(
gp.quicksum(f[i] * y[i] for i in range(n_fac))
+ gp.quicksum(c[i, j] * x[i, j] for i in range(n_fac) for j in range(n_cust)),
GRB.MINIMIZE)
model.addConstrs((gp.quicksum(x[i, j] for i in range(n_fac)) == 1
for j in range(n_cust)), name="assign")
model.addConstrs((x[i, j] <= y[i]
for i in range(n_fac) for j in range(n_cust)), name="link")
model.optimize()
assert model.Status == GRB.OPTIMAL
return model.ObjVal
f, c = generate_ufl(n_fac=15, n_cust=40, seed=7)
print(round(solve_ufl_direct(f, c), 4))
# Expected: 109.0666 -- identical to the Benders objective on the same seed.
Worked Example 2: Two-Stage Stochastic Capacity Planning (L-shaped Method)
Capacity planning under demand uncertainty: choose integer capacity units $y_i \in {0, \dots, u^{\max}}$ at plant i (unit cost $f_i$, unit size $\kappa$) before demand is known; after scenario s reveals demands $d_{js}$, ship $x_{ij} \ge 0$ at cost $c_{ij}$ and absorb any shortfall $g_j \ge 0$ at penalty $\rho$:
$$ \min_y ; \sum_i f_i y_i + \sum_s p_s, Q_s(y), \qquad Q_s(y) = \min_{x, g \ge 0} ; \sum_{ij} c_{ij} x_{ij} + \rho \sum_j g_j $$
$$ \text{s.t.} \quad \sum_j x_{ij} \le \kappa, y_i ;; \forall i ;;(\mu_{is} \le 0), \qquad \sum_i x_{ij} + g_j \ge d_{js} ;; \forall j ;;(\pi_{js} \ge 0). $$
The shortfall variables give relatively complete recourse — every scenario LP is feasible for every y, so the L-shaped method needs no feasibility cuts. From LP duality, any optimal dual pair $(\pi_s, \mu_s)$ at $\hat{y}$ satisfies $Q_s(\hat{y}) = \pi_s^\top d_s + \sum_i \mu_{is}, \kappa, \hat{y}_i$ and stays dual-feasible for every other y, giving the multi-cut optimality cut
$$ \theta_s ;\ge; \pi_s^\top d_s + \sum_i \mu_{is}, \kappa, y_i . $$
First the deterministic-equivalent (extensive form) model — the validation baseline and the thing the L-shaped method avoids building at scale: with S scenarios it has $S(|F||C| + |C|)$ second-stage variables, while the decomposition holds one scenario block at a time.
"""Two-stage capacity planning: extensive form (deterministic equivalent)."""
from __future__ import annotations
import gurobipy as gp
import numpy as np
from gurobipy import GRB
def generate_capacity_instance(m: int, n: int, n_scen: int, seed: int) -> dict:
"""m plants x n customers; uniform scenario demands; shortfall penalty rho."""
rng = np.random.default_rng(seed)
return {"fixed": rng.uniform(40.0, 80.0, size=m), "kappa": 10.0, "u_max": 6,
"cost": rng.uniform(1.0, 6.0, size=(m, n)),
"demand": rng.uniform(5.0, 25.0, size=(n_scen, n)),
"prob": np.full(n_scen, 1.0 / n_scen), "rho": 50.0}
def solve_extensive_form(inst: dict) -> tuple[float, np.ndarray]:
"""One monolithic MIP with a copy of the second stage per scenario."""
m, n = inst["cost"].shape
n_scen = inst["demand"].shape[0]
model = gp.Model("capacity-extensive")
model.Params.OutputFlag = 0
y = model.addVars(m, vtype=GRB.INTEGER, ub=inst["u_max"], name="y")
x = model.addVars(m, n, n_scen, lb=0.0, name="x")
g = model.addVars(n, n_scen, lb=0.0, name="g")
model.setObjective(
gp.quicksum(inst["fixed"][i] * y[i] for i in range(m))
+ gp.quicksum(inst["prob"][s] * (
gp.quicksum(inst["cost"][i, j] * x[i, j, s]
for i in range(m) for j in range(n))
+ inst["rho"] * gp.quicksum(g[j, s] for j in range(n)))
for s in range(n_scen)),
GRB.MINIMIZE)
model.addConstrs((gp.quicksum(x[i, j, s] for j in range(n))
<= inst["kappa"] * y[i]
for i in range(m) for s in range(n_scen)), name="capacity")
model.addConstrs((gp.quicksum(x[i, j, s] for i in range(m)) + g[j, s]
>= inst["demand"][s, j]
for j in range(n) for s in range(n_scen)), name="demand")
model.optimize()
assert model.Status == GRB.OPTIMAL
return model.ObjVal, np.array([round(y[i].X) for i in range(m)])
inst = generate_capacity_instance(m=4, n=6, n_scen=20, seed=3)
obj, y_ef = solve_extensive_form(inst)
print(round(obj, 4), y_ef)
# Expected: 773.5573 [6 4 1 0] -- eleven capacity units over three plants;
# the L-shaped method below reproduces this value exactly.
The L-shaped implementation uses the primal subproblem pattern: one scenario LP built once, re-targeted per (y, scenario) by updating constraint right-hand sides (Constr.RHS), with duals read from Constr.Pi. RHS updates preserve the basis, so each re-solve is a warm-started dual simplex — orders of magnitude cheaper than rebuilding.
"""Multi-cut L-shaped method for two-stage stochastic capacity planning."""
from __future__ import annotations
import math
import gurobipy as gp
import numpy as np
from gurobipy import GRB
def generate_capacity_instance(m: int, n: int, n_scen: int, seed: int) -> dict:
"""Identical generator and seed => identical instance to the extensive form."""
rng = np.random.default_rng(seed)
return {"fixed": rng.uniform(40.0, 80.0, size=m), "kappa": 10.0, "u_max": 6,
"cost": rng.uniform(1.0, 6.0, size=(m, n)),
"demand": rng.uniform(5.0, 25.0, size=(n_scen, n)),
"prob": np.full(n_scen, 1.0 / n_scen), "rho": 50.0}
class ScenarioSubproblem:
"""Second-stage LP built once; only constraint RHS change per (y, scenario)."""
def __init__(self, inst: dict) -> None:
m, n = inst["cost"].shape
self.inst = inst
self.model = gp.Model("second-stage")
self.model.Params.OutputFlag = 0
x = self.model.addVars(m, n, lb=0.0, name="x")
g = self.model.addVars(n, lb=0.0, name="g")
self.model.setObjective(
gp.quicksum(inst["cost"][i, j] * x[i, j] for i in range(m) for j in range(n))
+ inst["rho"] * gp.quicksum(g.values()), GRB.MINIMIZE)
self.cap = [self.model.addConstr(
gp.quicksum(x[i, j] for j in range(n)) <= 0.0, name=f"cap_{i}")
for i in range(m)]
self.dem = [self.model.addConstr(
gp.quicksum(x[i, j] for i in range(m)) + g[j] >= 0.0, name=f"dem_{j}")
for j in range(n)]
def solve(self, y_val: np.ndarray, demand: np.ndarray) -> tuple[float, np.ndarray, np.ndarray]:
"""Return (Q_s(y), pi over demand rows, mu over capacity rows)."""
for i, con in enumerate(self.cap):
con.RHS = self.inst["kappa"] * y_val[i]
for j, con in enumerate(self.dem):
con.RHS = demand[j]
self.model.optimize()
assert self.model.Status == GRB.OPTIMAL # complete recourse: never infeasible
return (self.model.ObjVal,
np.array([con.Pi for con in self.dem]),
np.array([con.Pi for con in self.cap]))
def solve_lshaped(inst: dict, tol: float = 1e-6, max_iters: int = 100) -> dict:
"""Multi-cut L-shaped loop; integer master re-solved each iteration."""
m, n = inst["cost"].shape
n_scen = inst["demand"].shape[0]
master = gp.Model("lshaped-master")
master.Params.OutputFlag = 0
y = master.addVars(m, vtype=GRB.INTEGER, ub=inst["u_max"], name="y")
theta = master.addVars(n_scen, lb=0.0, name="theta") # recourse cost >= 0
master.setObjective(
gp.quicksum(inst["fixed"][i] * y[i] for i in range(m))
+ gp.quicksum(inst["prob"][s] * theta[s] for s in range(n_scen)),
GRB.MINIMIZE)
sub = ScenarioSubproblem(inst)
best_ub, best_y, log = math.inf, None, []
for it in range(1, max_iters + 1):
master.optimize()
lb = master.ObjVal
y_val = np.array([round(y[i].X) for i in range(m)], dtype=float)
expected, n_cuts = 0.0, 0
for s in range(n_scen):
q_s, pi, mu = sub.solve(y_val, inst["demand"][s])
expected += inst["prob"][s] * q_s
if theta[s].X < q_s - tol * max(1.0, abs(q_s)):
master.addConstr(
theta[s] >= float(pi @ inst["demand"][s]) + gp.quicksum(
float(mu[i]) * inst["kappa"] * y[i] for i in range(m)),
name=f"opt_cut_s{s}_it{it}")
n_cuts += 1
ub = float(inst["fixed"] @ y_val) + expected
if ub < best_ub:
best_ub, best_y = ub, y_val.astype(int)
log.append((it, round(lb, 4), round(best_ub, 4), n_cuts))
if n_cuts == 0:
break # every theta_s is exact: optimal
return {"objective": best_ub, "y": best_y, "log": log}
inst = generate_capacity_instance(m=4, n=6, n_scen=20, seed=3)
res = solve_lshaped(inst)
print(round(res["objective"], 4), res["y"], len(res["log"]))
for row in res["log"]:
print(row)
# Expected: 773.5573 [6 4 1 0] in 6 iterations -- identical objective and y to
# the extensive form; the log shows LB rising monotonically to meet UB, with
# all 20 scenarios cutting in early rounds (20, 20, 20 cuts), then 18, 7, 0.
Two contrasts with Worked Example 1 are deliberate. The classic loop here re-solves the master MIP per iteration — affordable because the master is tiny and iterations are few; switch to the callback pattern of Example 1 when it is not. And the subproblem is solved in primal form with Pi duals — the right pattern when the primal is the natural model and complete recourse rules out rays. For scenario generation, sample-average approximation, and EVPI/VSS analysis around this model, see stochastic-optimization.
Advanced Techniques
Pareto-optimal cuts (Magnanti–Wong)
Degenerate subproblems have many optimal dual solutions, and the one the simplex method happens to return often yields a weak cut. Magnanti & Wong (1981, "Accelerating Benders decomposition: algorithmic enhancement and model selection criteria") select, among all duals optimal at $\hat{y}$, one maximizing the cut value at a core point $y^0$ (a point in the relative interior of conv(Y)): solve a second LP over the optimal-dual face. The resulting cut is Pareto-optimal — no other cut from the same subproblem dominates it.
"""Magnanti-Wong Pareto-optimal cut: re-optimize over the optimal-dual face."""
from __future__ import annotations
import gurobipy as gp
import numpy as np
from gurobipy import GRB
def pareto_optimal_cut(W: np.ndarray, T: np.ndarray, h: np.ndarray, q: np.ndarray,
y_hat: np.ndarray, z_hat: float,
y_core: np.ndarray) -> np.ndarray:
"""Among duals optimal at y_hat (value z_hat), maximize depth at y_core."""
model = gp.Model("magnanti-wong")
model.Params.OutputFlag = 0
u = model.addMVar(W.shape[0], lb=0.0, name="u")
model.addConstr(W.T @ u <= q, name="dual_feasibility")
model.addConstr((h - T @ y_hat) @ u == z_hat, name="optimal_face")
model.setObjective((h - T @ y_core) @ u, GRB.MAXIMIZE)
model.optimize()
assert model.Status == GRB.OPTIMAL
return u.X
# --- plant-opening instance with demand 8: y_hat=(1,0) is already feasible --
W = np.array([[1.0, 1.0], [-1.0, 0.0], [0.0, -1.0]])
T = np.array([[0.0, 0.0], [8.0, 0.0], [0.0, 8.0]])
h, q = np.array([8.0, 0.0, 0.0]), np.array([1.0, 2.0])
u_mw = pareto_optimal_cut(W, T, h, q, y_hat=np.array([1.0, 0.0]), z_hat=8.0,
y_core=np.array([0.5, 0.5]))
print(u_mw, float(u_mw @ h), (u_mw @ T))
# Expected: one maximizer is u = [2, 1, 0], giving the cut eta >= 16 - 8*y_0,
# which dominates the naive tight cut eta >= 8 everywhere with y_0 < 1: the
# subproblem at y_hat is degenerate, and MW picks the deepest dual.
Practical notes: the face constraint can make the MW problem unbounded when the optimal face is — add a normalization (bound u's components or its sum) or use the simpler $\varepsilon$-perturbation variant of Papadakos (2008). A fixed core point such as 0.5 * ones works for binary masters; updating it as a running average of incumbents is common.
Single-cut, multi-cut, and partial aggregation
Multi-cut (one theta_s per block) carries more information per iteration and typically cuts iteration counts by 5–10× at the price of a larger master; single-cut keeps the master lean but learns one aggregate hyperplane per round. Between them, partial aggregation clusters scenarios into K groups with one theta_k each — for thousands of scenarios, K in the tens is a strong default (Trukhanov, Ntaimo & Schaefer 2010, adaptive multicut). Choose by master-solve share of runtime: master-bound → aggregate more; iteration-bound → disaggregate.
Stabilization and in-out separation
The master's argmin jumps wildly between distant y in early iterations ("bang-bang" behavior), producing shallow cuts and the long flat tail typical of naive Benders. In-out separation (Ben-Ameur & Neto 2007; used to great effect in Fischetti, Ljubić & Sinnl 2017) separates not at the master optimum $\hat{y}$ but at an interpolation $\lambda \hat{y} + (1-\lambda) \tilde{y}$ toward a stabilizing point $\tilde{y}$ (e.g., the incumbent or a core point), with $\lambda \approx 0.5$ and a fallback to $\hat{y}$ when no cut is violated. Alternatives with the same intent: trust-region constraints on y (limit Hamming distance per iteration), level-method stabilization, and warm-starting eta/theta with cuts harvested from the LP relaxation before integrality is enforced.
Cuts at fractional points
Branch-and-Benders-cut as implemented in Worked Example 1 separates only at integer candidates (MIPSOL). Separating also at fractional LP solutions (MIPNODE + cbCut, as user cuts) tightens the tree bound — for UFLP the closed-form separation works unchanged at fractional y by replacing the open-facility minimum with the LP value $z_j(\hat{y}) = \min_i$ over a sorted scan. Standard recipe: run a classic-loop warm-up on the LP relaxation until the root bound stalls, load those cuts into the MIP, then go single-tree with MIPSOL lazy cuts plus throttled MIPNODE cuts (e.g., only at the root). See gurobi-advanced-features for callback plumbing details.
Integer or nonlinear subproblems
Classic cuts need LP duality. When the subproblem keeps discrete decisions: logic-based Benders (Hooker & Ottosson 2003) derives problem-specific cuts from an inference dual — standard in planning-plus-scheduling splits where the subproblem is a CP feasibility check; combinatorial Benders cuts (Codato & Fischetti 2006) exclude minimal infeasible assignments via no-good cuts sum(1 - y_i over ones) + sum(y_i over zeros) >= 1; the integer L-shaped method (Laporte & Louveaux 1993) handles integer recourse with optimality cuts built from a global lower bound L. All are much weaker per cut than LP-dual cuts — expect many more iterations and design the master to carry as much structure as possible. For convex nonlinear subproblems, generalized Benders (Geoffrion 1972) replaces LP duals with convex duality.
Practical Challenges
The first master solve is unbounded. eta (or theta_s) has no lower bound until the first optimality cut arrives. Always set a valid bound: 0 when recourse costs are nonnegative, per-block bounds like eta_j >= min_i c[i, j] (Worked Example 1), or the LP-relaxation recourse value. A tight initial bound also speeds up the first iterations measurably.
The subproblem reports status INF_OR_UNBD and no ray or dual is available. Presolve's dual reductions blur the infeasible/unbounded distinction, and barrier solves do not produce vertex certificates. Set DualReductions=0 plus InfUnbdInfo=1, and force simplex (Method=0 or 1) on subproblems whose rays or FarkasDual you intend to read.
Cuts with wrong signs silently break everything. A sign error in dualization makes cuts that are invalid (cut off the optimum) or vacuous (never bind). Discipline: derive the dual on paper for your exact inequality directions; then assert at runtime that every added cut is violated by the current master point and satisfied by a known feasible solution; and reproduce the monolithic optimum on small seeded instances before scaling. Gurobi's Pi convention — duals of <= rows are ≤ 0 in minimization, duals of >= rows are ≥ 0 — means you can use Pi values directly in the formulas of Worked Example 2 without manual sign flips.
Convergence stalls with a long flat tail. The textbook loop often closes 95% of the gap in a few iterations and crawls for hundreds more. In order of impact: switch to multi-cut, add Magnanti–Wong or in-out stabilization, keep a partial copy of linking constraints in the master (partial Benders retains some x-aggregates in the master to strengthen its relaxation), and harvest LP-relaxation cuts before going integer.
The callback returns a wrong optimum without any error. With LazyConstraints=1, correctness rests entirely
…(truncated)