Generalized Benders Decomposition

Generalized Benders Decomposition#

Generalized Benders Decomposition (GBD) [Geoffrion, 1972] solves a MINLP by alternating between two much smaller problems: an NLP with the integers fixed, and a MILP master that decides the next integer assignment. It is the nonlinear generalization of classical Benders [Benders, 1962], which requires the subproblem to be an LP — see the Benders tutorial for that case.

This is the answer to a question that comes up often: is there such a thing as a “local” MINLP solver? Yes — GBD and Outer Approximation are both decomposition schemes that never build a spatial branch & bound tree. On a convex MINLP they converge finitely to the global optimum; on a non-convex one they are local-solution heuristics with no certificate.

The two-stage picture#

Split the variables into a first stage \(y\) (the complicating ones — typically the integers) and a recourse stage \(x\) (continuous):

\[\min_{y \in Y,\; x} \; f(y) + g(y, x) \quad \text{s.t.}\quad h(y, x) \le 0\]

Fixing \(y = \hat{y}\) leaves a pure NLP in \(x\). Its optimal value defines the value function \(v(\hat y) = \min_x \{\, g(\hat y, x) : h(\hat y, x) \le 0 \,\}\), and the original problem collapses to \(\min_{y \in Y} f(y) + v(y)\) — an optimization over the integers alone, if only we knew \(v\).

We don’t. But when the model is convex, \(v\) is convex in \(y\), and the Lagrangian dual of the recourse NLP hands us a supporting hyperplane of \(v\) at \(\hat y\) for free: with multipliers \(\mu \ge 0\) from the recourse solve,

\[v(y) \;\ge\; L(\hat x, \hat y, \mu) + \nabla_y L(\hat x, \hat y, \mu)^{\!\top}(y - \hat y).\]

Accumulating these gives a piecewise-linear under-estimator of \(v\), which is exactly the MILP master:

\[\min_{y \in Y,\; \eta} \; f(y) + \eta \quad \text{s.t.}\quad \eta \ge \text{(all cuts so far)}.\]

The master value is a lower bound (it under-estimates \(v\)); each recourse solve gives a feasible point, hence an upper bound. The loop stops when they meet.

The algorithm, written out#

The mechanics are easiest to see with no solver in the loop at all. Take

\[\min_{y \in \{0,1\}^2,\; x \in [0,5]} \; 4y_1 + 3y_2 + x^2 \quad \text{s.t.}\quad x + 2y_1 + y_2 \ge 3,\]

whose recourse value and multiplier are both closed form:

\[v(y) = \big(3 - 2y_1 - y_2\big)_+^2, \qquad \mu(y) = 2\big(3 - 2y_1 - y_2\big)_+,\]

so the cut generated at \(\hat y\) is \(\eta \ge v(\hat y) - \mu(\hat y)\big[2(y_1 - \hat y_1) + (y_2 - \hat y_2)\big]\). With only four candidate \(y\) vectors we can solve the master by enumeration and watch the bounds close.

import itertools

import numpy as np

C = np.array([4.0, 3.0])   # first-stage cost
A = np.array([2.0, 1.0])   # coupling coefficients
D = 3.0                    # demand
YS = list(itertools.product([0, 1], repeat=2))


def value(y):
    '''Recourse optimal value with the first stage fixed.'''
    return max(0.0, D - A @ np.asarray(y, float)) ** 2


def multiplier(y):
    '''Optimal multiplier of the coupling row (the Lagrangian dual solution).'''
    return 2.0 * max(0.0, D - A @ np.asarray(y, float))


cuts = []                       # each entry: (const, grad) meaning eta >= const + grad . y
y_hat, ub, lb = (0, 0), np.inf, -np.inf

print(f"{'iter':>4}  {'y_hat':>8}  {'v(y_hat)':>9}  {'mu':>6}  {'LB':>9}  {'UB':>9}")
for it in range(1, 11):
    v, mu = value(y_hat), multiplier(y_hat)
    ub = min(ub, C @ np.array(y_hat, float) + v)          # incumbent from the recourse solve
    cuts.append((v + mu * (A @ np.array(y_hat, float)), -mu * A))

    # MILP master, solved here by enumerating the four y vectors.
    lb, y_next = min(
        (C @ np.array(y, float) + max([0.0] + [c + g @ np.array(y, float) for c, g in cuts]), y)
        for y in YS
    )
    print(f"{it:>4}  {str(y_hat):>8}  {v:>9.4f}  {mu:>6.3f}  {lb:>9.6f}  {ub:>9.6f}")
    if ub - lb <= 1e-9:
        break
    y_hat = y_next

print(f"\nconverged in {it} iterations: optimum {ub} at y={y_hat}")

iter     y_hat   v(y_hat)      mu         LB         UB
   1    (0, 0)     9.0000   6.000   4.000000   9.000000
   2    (1, 0)     1.0000   2.000   5.000000   5.000000

converged in 2 iterations: optimum 5.0 at y=(1, 0)

Two iterations, two cuts, and the master’s lower bound rises from 4 to 5 while the incumbent falls from 9 to 5. Note what the master never saw: the nonlinearity. All it ever handles is \(f(y)\) plus a few linear inequalities in \(y\) and \(\eta\).

solve_gbd does the same thing with real NLP solves and real multipliers.

import discopt.modeling as dm
from discopt.decomposition.benders import solve_gbd

m = dm.Model("toy")
y = m.binary("y", shape=(2,))
x = m.continuous("x", lb=0.0, ub=5.0)
m.first_stage(y)                          # declare the complicating variables
m.minimize(4 * y[0] + 3 * y[1] + x * x)
m.subject_to(x + 2 * y[0] + y[1] >= D)

r = solve_gbd(m, time_limit=60)
print(f"status={r.status}  objective={r.objective:.6f}  bound={r.bound:.6f}  "
      f"certified={r.gap_certified}")

status=optimal  objective=5.000000  bound=5.000000  certified=True

A worked example: capacity expansion with congestion#

A more realistic shape. Three candidate units, each with a fixed build cost \(f_i\) and a capacity \(\bar{u}_i\); the operating cost is convex quadratic in throughput (a congestion or efficiency-loss term), which is what makes the recourse an NLP rather than an LP:

\[\begin{split} \begin{aligned} \min \quad & \sum_i f_i y_i + \sum_i \big(a_i x_i + b_i x_i^2\big) \\ \text{s.t.}\quad & \sum_i x_i \ge D, \qquad x_i \le \bar{u}_i\, y_i, \qquad y_i \in \{0,1\},\; x_i \ge 0 \end{aligned} \end{split}\]

The binaries are the complicating variables — with \(y\) fixed this is a small convex QP.

FIXED = [3.5, 2.0, 4.5]     # build cost
LIN = [1.0, 2.0, 0.6]       # linear operating cost
QUAD = [0.30, 0.10, 0.45]   # congestion coefficient
CAP = [4.0, 5.0, 3.0]
DEMAND = 7.0


def capacity_expansion():
    m = dm.Model("capacity")
    y = m.binary("y", shape=(3,))
    x = m.continuous("x", shape=(3,), lb=0.0, ub=5.0)
    m.first_stage(y)
    m.minimize(
        sum(FIXED[i] * y[i] for i in range(3))
        + sum(LIN[i] * x[i] + QUAD[i] * x[i] * x[i] for i in range(3))
    )
    m.subject_to(sum(x[i] for i in range(3)) >= DEMAND)
    for i in range(3):
        m.subject_to(x[i] <= CAP[i] * y[i])
    return m


r_gbd = solve_gbd(capacity_expansion(), time_limit=120)
print(f"GBD: status={r_gbd.status}  objective={r_gbd.objective:.4f}  "
      f"bound={r_gbd.bound:.4f}  certified={r_gbd.gap_certified}")
print("built units:", [round(float(v)) for v in r_gbd.x["y"]])

GBD: status=optimal  objective=20.8000  bound=20.8000  certified=True
built units: [1, 1, 0]

Three ways to ask for it#

solve_gbd is the direct entry point. solve_benders inspects the model and dispatches to GBD when the subproblem is nonlinear (and to classical Benders when it is an LP), and model.solve(decomposition="benders") routes through the same dispatcher.

from discopt.decomposition.benders import solve_benders

r_dispatch = solve_benders(capacity_expansion(), time_limit=120)
r_kwarg = capacity_expansion().solve(decomposition="benders", time_limit=120)

print(f"solve_benders(...)                    -> {r_dispatch.objective:.4f}")
print(f"model.solve(decomposition='benders')  -> {r_kwarg.objective:.4f}")

solve_benders(...)                    -> 20.8000
model.solve(decomposition='benders')  -> 20.8000

Declaring the structure#

m.first_stage(y) annotates the complicating variables. If you omit it, detect_decomposition infers a structure — by default the integers become the first stage, which is the right guess for most MINLPs but not for every model.

from discopt import detect_decomposition

structure = detect_decomposition(capacity_expansion())
print("first-stage (complicating) variables:", structure.complicating_vars)
print("inferred from:", structure.source)

first-stage (complicating) variables: ['y']
inferred from: annotated

GBD vs. Outer Approximation vs. spatial B&B#

All three solve this model to the same optimum, by very different routes. Grossmann [2002] establishes the relationship between the first two: the OA master problem is at least as tight as the GBD master, which is at least as tight as the continuous relaxation,

\[z_{\text{OA}} \;\ge\; z_{\text{GBD}} \;\ge\; z_{\text{LP}},\]

because a single GBD cut is a particular non-negative aggregation of the per-constraint OA linearizations. GBD therefore usually needs more iterations than OA — the trade-off is that its master grows by exactly one row per iteration regardless of how many nonlinear constraints the model has, whereas OA adds one per constraint. GBD wins when the master must stay small; OA wins almost everywhere else.

The three agree on \(20.8\) below. They agree only to solver tolerance, not to the last bit: GBD stops at the default relative gap of \(10^{-4}\), and its incumbent carries the recourse NLP’s own feasibility tolerance, so the reported bound and incumbent can differ in the eighth digit — including in the direction that looks like an inverted certificate. Read the pair as “equal to within gap_tolerance”.

r_oa = capacity_expansion().solve(solver="mip-nlp", mip_nlp_method="oa", time_limit=120)
# solver="bb" forces spatial branch & bound. A bare solve() is re-routed to the
# MIP-NLP family on a convex MINLP like this one, which would make the "spatial
# B&B" row another outer-approximation run (it would report node_count = 0).
r_bb = capacity_expansion().solve(solver="bb", time_limit=120)
assert getattr(r_bb, "mip_nlp_trace", None) is None, "expected the spatial path"

rows = [
    ("GBD", r_gbd, None),
    ("OA", r_oa, r_oa.mip_nlp_trace["summary"]),
    ("spatial B&B", r_bb, None),
]
print(f"{'method':>12}  {'objective':>11}  {'bound':>11}  {'time (s)':>9}  {'nodes':>6}  {'masters':>8}")
for name, res, summary in rows:
    masters = "-" if summary is None else summary["mip_count"]
    print(f"{name:>12}  {res.objective:>11.4f}  {res.bound:>11.4f}  "
          f"{res.wall_time:>9.3f}  {res.node_count:>6}  {str(masters):>8}")

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

      method    objective        bound   time (s)   nodes   masters
         GBD      20.8000      20.8000      0.074       0         -
          OA      20.8000      20.8000      0.025       0         4
 spatial B&B      20.8000      20.8000      0.015       9         -

When GBD is the right tool#

Situation	Use
Convex MINLP, moderate number of nonlinear constraints	OA — tighter master, fewer iterations
Convex MINLP where the master must stay small (many nonlinear rows, few first-stage vars)	GBD
Two-stage stochastic program with linear recourse	Classical Benders [Benders, 1962]
Non-convex MINLP, global certificate required	Spatial B&B — `model.solve()`

The Decomposition Advisor will make this call for you; it ranks OA above GBD automatically on a convex model.

Note

model.solve() reaches spatial B&B on a non-convex model, as the last row says. On a convex MINLP a bare solve() is re-routed to the MIP–NLP family instead, and comes back with node_count = 0 and a mip_nlp_trace. To force the spatial path for a like-for-like comparison — as the cell above does — pass solver="bb" explicitly.

Caveats#

Warning

Convexity is what makes the bound rigorous. The cut is a supporting hyperplane of \(v\) only when \(v\) is convex. On a non-convex model solve_gbd runs in heuristic mode: it returns a feasible point and withholds the dual bound (bound is None). A None bound is the honest answer, not a failure — it means “no certificate”, and you should reach for spatial B&B if you need one.

Warning

Degenerate recourse. If the recourse NLP has an empty interior at some \(\hat y\) (Slater’s condition fails), no finite multiplier exists, and the Lagrangian cut is tight at \(\hat y\) but vacuous everywhere else. The master then re-proposes points it has already cut and the lower bound stalls — GBD exits with status="iteration_limit" and a weak-but-valid bound, and logs a warning saying exactly this. When the first stage is all 0/1, discopt adds a multiplier-free integer L-shaped cut [Laporte and Louveaux, 1993] alongside the Lagrangian one, which restores progress in most such cases.

Two structural limits of the current implementation: master-only constraints (those touching no recourse variable) must be linear, and infeasible recourse is handled by a no-good cut, which requires an all-0/1 first stage. A model with general-integer first-stage variables therefore needs relatively complete recourse — every first-stage assignment must leave the recourse problem feasible.

Generalized Benders Decomposition

Contents

Generalized Benders Decomposition#

The two-stage picture#

The algorithm, written out#

A worked example: capacity expansion with congestion#

Three ways to ask for it#

Declaring the structure#

GBD vs. Outer Approximation vs. spatial B&B#

When GBD is the right tool#

Caveats#

Further reading#