Benders Decomposition

Benders Decomposition#

Benders decomposition [Benders, 1962] solves block-angular / two-stage problems by splitting the variables into complicating (first-stage / master) variables and recourse (subproblem) variables. For a fixed master point the recourse subproblem is an independent program; its optimal dual (or KKT multipliers) yields a cut that under-estimates the recourse cost, and its infeasibility yields a cut that excludes the offending master point. Iterating the master/subproblem loop drives the master’s lower bound up to meet the incumbent’s upper bound.

discopt’s solve_benders implements classical Benders for linear-recourse problems (linear constraints and objective, integer variables first-stage so the recourse is a continuous LP) and Generalized Benders Decomposition [Geoffrion, 1972] for a convex nonlinear recourse subproblem — the solver detects nonlinearity and dispatches to GBD automatically. Because each cut is a valid global under-estimator (LP weak duality, or — for GBD — the envelope theorem on the convex value function), the master objective is a rigorous lower bound whenever the model is convex; the solver never reports an unsound bound. For a broad survey of the method and its modern variants, see Rahmaniani et al. [2017].

import discopt.modeling as dm

A two-stage facility problem#

Open a facility (\(y \in \{0, 1\}\), first stage) and then ship from it (\(x_1, x_2 \geq 0\), recourse) to meet a demand of 3 units. Shipping is only possible from an open facility (\(x_i \le 5y\)).

\[\min_{y, x}\; 2y + x_1 + x_2 \quad\text{s.t.}\quad x_1 + x_2 \ge 3,\; x_1 \le 5y,\; x_2 \le 5y.\]

If \(y = 0\) the recourse is infeasible (no shipping can meet demand), so Benders adds a feasibility cut forcing \(y = 1\); then the recourse cost is 3 and the optimum is \(2 + 3 = 5\).

m = dm.Model("facility")
y = m.binary("y")
x1 = m.continuous("x1", lb=0, ub=10)
x2 = m.continuous("x2", lb=0, ub=10)
m.minimize(2 * y + x1 + x2)
m.subject_to(x1 + x2 >= 3)
m.subject_to(x1 <= 5 * y)
m.subject_to(x2 <= 5 * y)

result = m.solve(decomposition="benders")
print("status   :", result.status)
print("objective:", round(result.objective, 4))
print("bound    :", round(result.bound, 4))
print("y =", result.x["y"], " x1 =", result.x["x1"], " x2 =", result.x["x2"])

Benders iteration 1 added no separating cut (degenerate or missing duals); stall=1.

status   : optimal
objective: 5.0
bound    : 5.0
y = 1.0  x1 = 3.0  x2 = 0.0

The reported bound equals the objective: the optimality gap is certified. This matches the monolithic solve:

mono = m.solve()
print("monolithic objective:", round(mono.objective, 4))

monolithic objective: 5.0

Declaring the decomposition structure#

By default the complicating variables are the integer variables — the canonical Benders split for two-stage mixed-integer programs. For a purely continuous two-stage model, or to override the split, annotate the model directly:

  • m.first_stage(u) / m.second_stage(v) — mark complicating vs recourse variables.

  • m.mark_coupling(constraint) — mark a linking constraint (used by Lagrangian relaxation).

discopt.detect_decomposition(model) resolves these annotations and reports the structure it will use.

import discopt

lp = dm.Model("continuous_two_stage")
u = lp.continuous("u", lb=0, ub=10)
v = lp.continuous("v", lb=0, ub=10)
lp.first_stage(u)  # u is the complicating variable
lp.minimize(u + 2 * v)
lp.subject_to(u + v >= 4)

print(discopt.detect_decomposition(lp).summary())
print("objective:", round(lp.solve(decomposition="benders").objective, 4))

Benders iteration 1 added no separating cut (degenerate or missing duals); stall=1.

DecompositionStructure (annotated)
  blocks: 1 (separable=False)
  complicating vars: 1
  coupling constraints: 0
objective: 4.0

When to use Benders#

Benders shines when fixing a few complicating variables makes the remaining problem an easy, possibly decomposable, subproblem — two-stage stochastic programs, network/facility design, and capacity-expansion models. The implementation handles:

  • classical Benders: linear constraints and objective, with all integer variables first-stage (continuous LP recourse), and

  • Generalized Benders [Geoffrion, 1972]: a convex nonlinear recourse subproblem (cuts from the recourse NLP’s KKT multipliers), dispatched automatically when the model is nonlinear.

The reported lower bound is rigorous when the model is convex. For a nonconvex model GBD still returns a heuristic incumbent but reports bound=None (the value function is no longer guaranteed convex, so no valid global cut can be certified); use m.solve() (spatial branch & bound) for a guaranteed bound there, or m.solve(solver="mip-nlp", mip_nlp_method="oa") (Outer Approximation) [Duran and Grossmann, 1986] for convex MINLP.

GBD gets a tutorial of its own — see Generalized Benders Decomposition, which derives the Lagrangian cut, runs the loop by hand, and compares GBD against OA and spatial B&B on the same model.

gbd = dm.Model("facility_quadratic")
yb = gbd.binary("y")
xs = gbd.continuous("x", lb=0, ub=5)
gbd.first_stage(yb)
gbd.minimize(2 * yb + xs**2)  # convex quadratic recourse cost
gbd.subject_to(xs >= 2)  # must ship >= 2 units
gbd.subject_to(xs <= 5 * yb)  # only from an open facility

res = gbd.solve(decomposition="benders")  # nonlinear -> Generalized Benders
print("status   :", res.status)
print("objective:", round(res.objective, 4))
print("bound    :", round(res.bound, 4))
print("monolithic:", round(gbd.solve().objective, 4))

OA: generating OA cuts for 0 of 2 rows (2 convex; the rest are either nonconvex or already exact in the master)

status   : optimal
objective: 6.0
bound    : 6.0
monolithic: 6.0