Outer Approximation Solver

Outer Approximation Solver#

This tutorial introduces discopt’s general-purpose Outer Approximation (OA) solver for Mixed-Integer Nonlinear Programming (MINLP). OA decomposes a MINLP into alternating NLP subproblems and MILP master problems, and is especially effective for convex MINLP problems with many binary/integer variables [Duran and Grossmann, 1986].

Selecting the solver#

OA is one of the MIP–NLP decomposition methods, selected with solver="mip-nlp":

result = model.solve(solver="mip-nlp", mip_nlp_method="oa")

Note

The old spelling gdp_method="oa" is deprecated and no longer selects OA. It now emits a DeprecationWarning and is interpreted as gdp_method="big-m", a GDP reformulation — so a solve written that way silently runs spatial branch & bound. The gdp_method= and solver="mip-nlp" axes are orthogonal: gdp_method says how disjunctions are reformulated, solver says which algorithm solves the result. (gdp_method="loa" is a different thing again — native logic-based OA, which linearizes inside disjunctions; see the GDP tutorial.)

When to use OA vs. Branch & Bound#

Criterion

OA

Spatial B&B

Convex MINLP

Preferred — finite convergence to global optimum

Valid but often slower

Non-convex MINLP

Local optimum only — cuts may be invalid

Global optimum with relaxation guarantees

Many binary variables

Strong — MILP solver handles combinatorics

Can struggle with large B&B tree

Few integers, hard NLP

Overkill — NLP dominates

Preferred

Algorithm Overview#

The OA algorithm [Duran and Grossmann, 1986, Fletcher and Leyffer, 1994] works as follows:

  1. Initialize: Solve continuous NLP relaxation → first linearization point

  2. Loop:

    • Solve MILP master (linear constraints + accumulated OA cuts) → lower bound

    • Fix integers at master solution, solve NLP subproblem → upper bound

    • Generate OA cuts (tangent hyperplanes) at NLP solution

    • If NLP infeasible: add feasibility cuts [Fletcher and Leyffer, 1994] or no-good cuts

  3. Terminate when gap converges or master becomes infeasible

Because the master accumulates supporting hyperplanes of a convex feasible set, its optimum is a valid lower bound, and the fixed-integer NLP gives a valid upper bound. On a convex model the two meet in finitely many iterations.

import discopt.modeling as dm

Example 1: Simple Convex MINLP#

\[\min_{x,y} \quad x_1^2 + x_2^2 + x_3\]
\[\text{s.t.} \quad x_1 + x_2 \geq 1, \quad x_1^2 + x_2 \leq 3, \quad x_3 \in \{0,1\}\]

The optimal solution is \(x_1 = x_2 = 0.5\), \(x_3 = 0\), with objective 0.5.

from discopt.modeling.examples import example_simple_minlp

m = example_simple_minlp()
result = m.solve(solver="mip-nlp", mip_nlp_method="oa", time_limit=60)

print(f"Status: {result.status}")
print(f"Objective: {result.objective:.4f}")
print(f"Solution: {result.x}")
print(f"Wall time: {result.wall_time:.2f}s")

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

Model: textbook
  Variables: 3 (2 continuous, 1 integer/binary)
  Constraints: 2
  Objective: minimize (((x1 ** 2) + (x2 ** 2)) + x3)
  Parameters: 0
Status: optimal
Objective: 0.5000
Solution: {'x1': array(0.5), 'x2': array(0.5), 'x3': array(0.)}
Wall time: 0.01s

Watching the loop: a problem that actually iterates#

The model above is so small that OA closes the gap on its first master/subproblem pair — which teaches nothing about the alternation. synthes1, the first process-synthesis test problem of Duran and Grossmann [1986], is the standard textbook instance: three binaries choosing between process alternatives, three continuous flows, and logarithmic (concave, hence convex when negated in a minimization) yield terms.

\[\begin{split} \begin{aligned} \min \quad & 5y_1 + 6y_2 + 8y_3 + 10x_1 - 7x_3 - 18\ln(x_2+1) - 19.2\ln(x_1-x_2+1) + 10\\ \text{s.t.}\quad & 0.8\ln(x_2+1) + 0.96\ln(x_1-x_2+1) - 0.8x_3 \geq 0\\ & \ln(x_2+1) + 1.2\ln(x_1-x_2+1) - x_3 - 2y_3 \geq -2\\ & x_2 \le x_1, \quad x_2 \le 2y_1, \quad x_1 - x_2 \le 2y_2, \quad y_1 + y_2 \le 1 \end{aligned} \end{split}\]

The reported optimum is \(6.00976\).

def synthes1():
    '''Duran & Grossmann (1986), test problem 1 — a convex MINLP.'''
    m = dm.Model("synthes1")
    x1 = m.continuous("x1", lb=0, ub=2)
    x2 = m.continuous("x2", lb=0, ub=2)
    x3 = m.continuous("x3", lb=0, ub=1)
    y1, y2, y3 = m.binary("y1"), m.binary("y2"), m.binary("y3")
    m.minimize(
        5 * y1 + 6 * y2 + 8 * y3 + 10 * x1 - 7 * x3
        - 18 * dm.log(x2 + 1) - 19.2 * dm.log(x1 - x2 + 1) + 10
    )
    m.subject_to(0.8 * dm.log(x2 + 1) + 0.96 * dm.log(x1 - x2 + 1) - 0.8 * x3 >= 0)
    m.subject_to(dm.log(x2 + 1) + 1.2 * dm.log(x1 - x2 + 1) - x3 - 2 * y3 >= -2)
    m.subject_to(x2 - x1 <= 0)
    m.subject_to(x2 - 2 * y1 <= 0)
    m.subject_to(x1 - x2 - 2 * y2 <= 0)
    m.subject_to(y1 + y2 <= 1)
    return m


result_oa = synthes1().solve(solver="mip-nlp", mip_nlp_method="oa", time_limit=60)
print(f"Status: {result_oa.status}, objective: {result_oa.objective:.6f}")
print(f"Dual bound: {result_oa.bound:.6f}, certified: {result_oa.gap_certified}")

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

Status: optimal, objective: 6.009759
Dual bound: 6.009759, certified: True

The iteration trace#

Every MIP–NLP solve attaches a structured mip_nlp_trace to the result: a summary dict of counters and an iterations list with the lower bound, upper bound, gap and cuts added at each pass. This is the algorithm made visible — watch the bounds squeeze together from both sides.

trace = result_oa.mip_nlp_trace
summary = trace["summary"]

print(f"{'iter':>4}  {'lower bound':>12}  {'upper bound':>12}  {'gap':>10}  {'cuts':>5}")
for it in trace["iterations"]:
    ub = it["ub"] if it["ub"] is not None else float("inf")
    print(
        f"{it['index']:>4}  {it['lb']:>12.6f}  {ub:>12.6f}  "
        f"{it['gap']:>10.2e}  {it['cuts_added']:>5}"
    )

print(f"\nterminated because: {trace['termination_reason']}")
print(f"MILP masters solved:     {summary['mip_count']}")
print(f"NLP subproblems solved:  {summary['nlp_subproblem_count']}")
print(f"OA cuts generated:       {summary['cut_count']}")

iter   lower bound   upper bound         gap   cuts
   0      1.411895      7.092732    8.01e-01      3
   1      1.687599      6.009759    7.19e-01      3
   2      5.312482      6.009759    1.16e-01      3
   3      6.009759      6.009759    0.00e+00      3

terminated because: gap
MILP masters solved:     4
NLP subproblems solved:  4
OA cuts generated:       15

Example 2: Quadratic Objective with Binary Activation#

A common pattern in process design: a quadratic cost with binary decisions that activate or deactivate options.

\[\min \quad (x - 3)^2 + 2y_1 + 3y_2\]
\[\text{s.t.} \quad y_1 + y_2 \leq 1, \quad x \leq 2y_1 + 4y_2, \quad x \in [0, 5],\; y_i \in \{0,1\}\]
m = dm.Model("activation")
x = m.continuous("x", lb=0, ub=5)
y1 = m.binary("y1")
y2 = m.binary("y2")

m.minimize((x - 3) ** 2 + 2 * y1 + 3 * y2)
m.subject_to(y1 + y2 <= 1)
m.subject_to(x <= 2 * y1 + 4 * y2)

result = m.solve(solver="mip-nlp", mip_nlp_method="oa", time_limit=60)
print(f"Status: {result.status}, Objective: {result.objective:.4f}")
print(f"x={result.x['x']:.3f}, y1={result.x['y1']:.0f}, y2={result.x['y2']:.0f}")

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: 3.0000
x=3.000, y1=0, y2=1

Comparing OA with Branch & Bound#

Both find the same optimum on a convex MINLP — they differ in how they get there. OA spends its effort in a MILP master and a handful of NLPs and never builds a search tree (node_count is 0); spatial B&B builds a tree and reports nodes.

Important

To compare against spatial B&B you must ask for it explicitly with solver="bb". A bare model.solve() on synthes1 does not run spatial B&B: discopt detects that the model is a convex MINLP and routes it to the MIP–NLP family, so the result comes back with node_count = 0 and a mip_nlp_trace whose selected_strategy is "oa". Comparing solve(solver="mip-nlp", ...) against a bare solve() therefore compares OA with itself.

At synthes1’s size the two are within noise of each other on wall time — three binaries is not enough combinatorics for either approach to pull ahead. The separation appears as the integer count grows: OA’s work per iteration is one MILP, while B&B’s tree can grow exponentially in the binaries.

# solver="bb" forces spatial branch & bound. A bare solve() would be re-routed to
# the MIP-NLP family for this convex model -- see the note above.
result_bb = synthes1().solve(solver="bb", time_limit=60)

assert getattr(result_bb, "mip_nlp_trace", None) is None, (
    "expected the spatial path, but a MIP-NLP trace came back"
)

print(f"OA:  obj={result_oa.objective:.6f}  bound={result_oa.bound:.6f}  "
      f"time={result_oa.wall_time:.3f}s  nodes={result_oa.node_count}  "
      f"masters={result_oa.mip_nlp_trace['summary']['mip_count']}")
print(f"B&B: obj={result_bb.objective:.6f}  bound={result_bb.bound:.6f}  "
      f"time={result_bb.wall_time:.3f}s  nodes={result_bb.node_count}")

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

OA:  obj=6.009759  bound=6.009759  time=0.033s  nodes=0  masters=4
B&B: obj=6.009759  bound=6.009759  time=0.096s  nodes=5

Extended Cutting Plane (ECP) Mode#

The Extended Cutting Plane method [Westerlund and Pettersson, 1995] avoids NLP subproblem solves entirely. Instead, it generates cuts at the MILP master solution for violated nonlinear constraints. This is simpler but converges more slowly — the counters below make the trade-off concrete: ECP solves zero NLP subproblems, and pays for that with more master solves.

result_ecp = synthes1().solve(solver="mip-nlp", mip_nlp_method="ecp", time_limit=60)
ecp = result_ecp.mip_nlp_trace["summary"]
oa = result_oa.mip_nlp_trace["summary"]

print(f"ECP: status={result_ecp.status}, objective={result_ecp.objective:.6f}")
print(f"{'':>6}{'masters':>9}{'NLPs':>7}{'cuts':>7}")
print(f"{'OA':>6}{oa['mip_count']:>9}{oa['nlp_subproblem_count']:>7}{oa['cut_count']:>7}")
print(f"{'ECP':>6}{ecp['mip_count']:>9}{ecp['nlp_subproblem_count']:>7}{ecp['cut_count']:>7}")

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

ECP: status=optimal, objective=6.009759
        masters   NLPs   cuts
    OA        4      4     15
   ECP        8      0     16

Equality Relaxation#

For problems with nonlinear equality constraints, the OA linearizations can make the MILP master infeasible. The Equality Relaxation strategy [Viswanathan and Grossmann, 1990] relaxes these to inequalities in the OA cuts, restoring master feasibility.

Warning

A nonlinear equality makes the feasible set non-convex, so the convex-case convergence guarantee no longer applies. Equality relaxation (and the augmented penalty that usually accompanies it) is a robustness heuristic: it keeps the master solvable, but it does not restore the guarantee. Use spatial B&B — the default model.solve() — when you need a global certificate.

m = dm.Model("eq_relax")
x = m.continuous("x", lb=0, ub=2)
y = m.binary("y")
m.minimize(x**2 + y)
m.subject_to(x**2 - y == 0)  # nonlinear equality

result_er = m.solve(
    solver="mip-nlp",
    mip_nlp_method="oa",
    equality_relaxation=True,
    time_limit=60,
)
print(f"Status: {result_er.status}, Objective: {result_er.objective:.6f}")

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

Status: optimal, Objective: 0.000000

Other MIP–NLP methods#

mip_nlp_method selects among several members of the same family:

Method

What it does

"oa"

Outer approximation: MILP master + fixed-integer NLP [Duran and Grossmann, 1986]

"ecp"

Extended cutting plane: no NLP subproblems [Westerlund and Pettersson, 1995]

"goa"

Global OA — a convex under-estimator is built first, so the cuts stay valid on a non-convex model

"fp"

Feasibility pump: chases a feasible point rather than an optimal one

"lp_nlp_bb"

Single-tree LP/NLP branch & bound [Quesada and Grossmann, 1992] — OA cuts injected as lazy constraints into one MILP tree instead of re-solving the master from scratch. Requires milp_solver="gurobi".

Level-method regularization [Kronqvist et al., 2020], which stabilizes the integer selection between iterations, is available as add_regularization=True with strength level_coef.

Configuration Options#

Parameter

Default

Description

solver="mip-nlp"

—

Select the MIP–NLP decomposition family

mip_nlp_method

"oa"

"oa", "ecp", "goa", "fp", "lp_nlp_bb"

time_limit

3600

Wall-clock time limit (seconds)

gap_tolerance

1e-4

Relative optimality gap

max_iterations

100

Maximum OA iterations

equality_relaxation

False

Relax nonlinear equalities (heuristic — see the warning above)

feasibility_cuts

True

Gradient-based feasibility cuts

add_no_good_cuts

False

Exclude visited integer assignments

add_regularization

False

Level-method regularization [Kronqvist et al., 2020]

level_coef

—

Regularization strength

init_strategy

"rNLP"

"rNLP", "initial_binary", "max_binary", "fp"

nlp_solver

"pounce"

NLP backend: "pounce", "ipm", "ipopt", "cyipopt", "sparse_ipm"

milp_solver

auto

MILP backend; "gurobi" is required by lp_nlp_bb and solution_pool

All of these are passed straight through model.solve(...) as keyword arguments.

Limitations#

  • Non-convex problems: OA only guarantees local optimality, because a linearization of a non-convex constraint can cut off feasible points. For a global certificate on a non-convex MINLP, use the default spatial B&B solver (model.solve()), or mip_nlp_method="goa", which linearizes a convex under-estimator instead.

  • Nonlinear equalities are non-convex by construction; see the equality-relaxation warning above.

  • Single-tree mode needs a commercial MILP backend: lp_nlp_bb is implemented but relies on Gurobi’s lazy-constraint callbacks (milp_solver="gurobi"). Without it, OA runs in multi-tree mode and re-solves the master each iteration.

For a comprehensive review and comparison of convex MINLP solvers, see [Kronqvist et al., 2019]. The Bonmin solver [Bonami et al., 2008] provides multiple algorithmic variants including OA, NLP-BB, and hybrid approaches. A related decomposition, Generalized Benders [Geoffrion, 1972], cuts with Lagrangian duals of the NLP rather than per-constraint linearizations; its master is smaller but weaker than OA’s [Grossmann, 2002] — see the GBD tutorial.