Global optimization: what discopt can and can’t certify

Global optimization: what discopt can and can’t certify#

discopt’s branch-and-bound engine is a deterministic global optimizer: when it returns status="optimal" it also returns a certificate — a valid lower bound on the true optimum together with a feasible incumbent — so the reported solution is provably within the optimality gap of the global optimum [Belotti et al., 2013, Tawarmalani and Sahinidis, 2005]. This is fundamentally different from a local NLP solve, which only finds a stationary point and can get stuck in a poor local minimum on a nonconvex problem.

That guarantee is not free, and it does not apply to every model you can write. This guide explains which problems are amenable to certified global solving in discopt, which are clearly not, and what to do about the gap between the two. It complements rather than repeats the mechanics covered elsewhere: the AMP global solver, convexity detection, bound tightening, and the relaxation internals.

Did I actually get a certified global optimum?#

Every SolveResult carries the evidence. You do not have to take status on faith — inspect the bound and the gap:

result = model.solve()

result.status          # "optimal" only when the gap is closed within tolerance
result.objective       # the incumbent (best feasible objective found)
result.bound           # the dual bound (valid LB for min / UB for max)
result.gap             # |incumbent - bound| / |incumbent|
result.gap_certified   # True when `bound` is a *sound* global bound

A result is a certified global optimum when status == "optimal", gap_certified is True, and gap is at or below your tolerance. If instead you see status == "feasible", a non-None gap that never closed, or gap_certified is False, then discopt found a feasible point but could not prove it is global — treat it as a (possibly very good) heuristic solution, not a certificate.

Important

gap_certified is the single most important field to check on a nonconvex model. A small gap with gap_certified=False means the bound itself is not trustworthy as a global guarantee.

What makes a problem amenable#

Certified global solving in discopt rests on two requirements. Meet both and the spatial branch-and-bound will converge with a sound certificate; violate either and the guarantee weakens or disappears.

1. Finite variable bounds#

The relaxations that produce the dual bound are built over the variable box. A nonconvex term such as a product \(x\,y\) is replaced by its McCormick envelope [McCormick, 1976], whose tightness depends entirely on the finite bounds \([x^L,x^U]\times[y^L,y^U]\). Without finite bounds there is no valid envelope and no valid global bound. discopt warns when declared bounds are missing or enormous (it treats \(\pm 10^{20}\) as “unbounded”), and a model left that way will typically return status="feasible" with no certified gap rather than status="optimal".

Give every variable that participates in a nonconvex term a realistic finite range:

x = m.continuous("x", lb=0, ub=1000)     # bounded → envelopes are valid & tight
theta = m.continuous("theta", lb=-3.15, ub=3.15)

Bounds need not be tight by hand — discopt’s FBBT/OBBT presolve will contract them automatically — but they must be finite.

2. Factorable structure from supported intrinsics#

discopt builds relaxations by decomposing each expression into elementary operations it knows how to under- and over-estimate — factorable programming [Tawarmalani and Sahinidis, 2005]. The relaxation library covers a broad menu of intrinsics, each with a convex/concave envelope or αBB underestimator [Adjiman et al., 1998]:

  • Bilinear, trilinear, and general multilinear products (McCormick / RLT [Sherali and Adams, 1990]);

  • Polynomials and quadratics, including QCQP (PSD/SOC cuts);

  • Fractional terms — reciprocals, ratios of products, fractional powers;

  • Transcendentals — exp, log, the trig family, asin/acos/acosh, erf; and

  • Saturating activations — sigmoid, tanh, softplus (used by neural-network embeddings).

If your model is composed of these building blocks over bounded variables, it is amenable. Integer and binary variables are fully supported — they are branched just like spatial branching on continuous nonconvex terms — as is generalized disjunctive programming once reformulated.

Amenable problem classes (worked elsewhere)#

Class

discopt support

See

Bilinear / QCQP

McCormick + RLT/PSD/SOC cuts

Solver Internals: BARON-Parity Features

Pooling / blending

pq-formulation builder

The pooling problem and the pq-formulation (discopt.pooling)

Signomial / posynomial

factorable lift; GP log-transform

Geometric programming

Polynomial / concave

piecewise-McCormick, αBB

Certified-global MINLP with AMP

MINLP (general nonconvex)

spatial B&B; AMP for hard cases

Mixed-Integer Nonlinear Programming (MINLP)

Embedded ML models

big-M / full-space formulations

Neural Network Embedding

What is clearly not amenable#

The following do not yield a certified global optimum. In most cases discopt will still return a feasible point (from its local solves and primal heuristics), but the global guarantee is gone — so do not read result.objective as global.

  • Unbounded (or effectively unbounded) variables in nonconvex terms. No finite box ⇒ no valid envelope ⇒ no sound bound. This is the most common cause of a status="feasible", gap_certified=False result. Fix: add finite bounds.

  • Black-box / user-defined functions (udf, custom, external, CustomCall). These call arbitrary Python and have no algebraic relaxation, so the bounding step cannot reason about them. They are evaluated pointwise during local NLP solves only — useful for modeling convenience, but they reduce the run to local optimization with no global certificate. If the body traces through the reduced-space MCBox type you do get a certificate — see Reduced-Space Global Optimization of Hidden-Function Models. If it does not, solver="direct" gives you a systematic global search over the box in place of a single local solve (Black-box global search with DIRECT (solver="direct")) — still with no certificate, but a much better answer on a multimodal objective.

    The degradation is whole-model, not per-term. One opaque body anywhere in the objective or a constraint routes the entire solve to the local NLP path, so the algebraic structure in the rest of the model stops being exploited too. If only a small part of your model is genuinely opaque, it is worth asking whether that part can be written algebraically (or traced through MCBox) before accepting the loss for everything else.

    solver="amp" refuses such a model outright rather than returning something: it certifies by linearizing a partitioned relaxation, and there is nothing to linearize. That is a deliberate refusal, not a limitation to work around.

  • External simulators, table lookups, or interpolated data used as objective or constraint terms. Same reason as black-box functions: discopt cannot construct a valid under/over-estimator for something it cannot see symbolically. External Functions: Grey-Box Models with Your Own Derivatives shows how to embed one with dm.external(fn, jac=..., hess=...) — you supply the derivatives, the block becomes twice differentiable and solves on the NLP path, and it is status="feasible" with no dual bound, by contract. Supplying derivatives buys you a solvable block, not a boundable one; a certificate would additionally need a valid interval or Lipschitz enclosure of the block over a box, which discopt cannot obtain from a function it can only evaluate pointwise and would have no way to verify if you asserted it.

  • Operators with no implemented envelope, or genuinely non-factorable expressions. If an intrinsic is outside the supported menu above, or an expression cannot be decomposed into supported elementary operations, no sound relaxation is available.

  • Discontinuous or pathologically non-smooth terms (step functions, floor, expressions with poles inside the variable box, regions where the function is undefined / returns NaN). Convex relaxation theory assumes the function is defined and finite on the box.

  • Very high-dimensional nonconvexity. This is a tractability limit, not a correctness one: the certificate remains sound, but spatial branch-and-bound scales poorly with the number of nonconvex terms and integer variables. A model with hundreds of independent bilinear terms or a deep combinatorial structure may not close its gap within any reasonable time budget. Expect status="feasible" at the time limit with a non-zero gap.

Warning

A black-box udf does not raise an error — it silently turns a global solve into a local one. If you need a certificate, keep the model algebraic and built from the supported intrinsics. If you cannot, solver="direct" at least searches the box instead of taking the first local minimum it finds — but it does not restore the certificate, and nothing else will either.

A quick amenability checklist#

Before expecting a certified global optimum, confirm:

  1. Bounds. Does every variable in a nonlinear term have finite lb/ub? (No bounds warning at solve time.)

  2. Algebraic. Is every nonlinearity written with discopt intrinsics — no udf/custom, no external calls, no table lookups?

  3. Supported. Are all operators in the relaxation menu (products, polynomials, fractional, the transcendental/activation list)?

  4. Smooth & defined. Is every function finite and free of poles/NaN across the whole box (no floor/step inside nonconvex terms)?

  5. Sized sensibly. Is the count of nonconvex terms and integers small enough to close the gap in your time budget?

If all five are “yes,” solve() should return status="optimal" with gap_certified=True. If one is “no,” see the remedies below.

Remediation#

  • Add or tighten bounds. Provide finite lb/ub; let FBBT/OBBT presolve contract them further.

  • Reformulate to expose structure. Use the pooling pq-formulation for blending, the geometric-programming log-transform for posynomials, or GDP for logical/disjunctive structure — each turns a hard or unsupported form into one discopt relaxes well.

  • Replace black boxes with algebra. If a udf wraps something expressible in intrinsics, write it out; if it wraps a trained model, embed it as algebraic NN constraints instead.

  • Escalate to AMP for hard nonconvex MINLPs. The Adaptive Multivariate Partitioning solver [Nagarajan et al., 2019] (model.solve(solver="amp")) refines a piecewise relaxation where the gap is largest and is often the better choice when the default spatial branch-and-bound stalls. Tighter cuts (rlt_cuts=True, PSD/SOC separators) can also close stubborn gaps — see Solver Internals: BARON-Parity Features.

  • Register a composite the relaxer should treat as one atom. Factorable relaxation decomposes an expression term by term, which throws away cancellation between terms. dm.register_function(name, lower) declares a named univariate composite whose definition is built from dm.* primitives; the relaxation layer then envelopes the composite as a single atom instead of relaxing its expansion:

    import discopt.modeling as dm

rk = dm.register_function(
    "redlich_kister",
    lambda x: x * (1 - x) * (L0 + L1 * (2 * x - 1)),
    description="binary excess Gibbs energy",
)
m.constraint(g == rk(x))          # relaxed as one atom, not four products

    The name must not collide with a built-in operator — register_function raises rather than silently re-defining an operator’s envelope process-wide. Re-binding an existing name needs replace=True, and expressions already built keep the old body and lose the atom envelope. dm.registered_names() lists what is registered. This is the mechanism plugins use to hand the relaxer structure they know and generic factorable reformulation cannot recover (#1248).

  • Accept a local solution deliberately. If global is genuinely out of reach, read result.objective as a high-quality heuristic value and report it as such — but check gap_certified so you know that is what you have.

Engine choice at a glance#

Situation

Recommended call

Convex, or bounded factorable nonconvex (default)

model.solve()

Hard nonconvex MINLP, default B&B stalls

model.solve(solver="amp")

Hard nonconvex MINLP with Gurobi available for MILP masters

model.solve(solver="amp", milp_solver="gurobi")

Cross-check against, or hand the whole model to, SCIP (reported, not certified)

model.solve(solver="scip") (pip install discopt[scip])

Stubborn bilinear/QCQP gap

model.solve(rlt_cuts=True) (+ PSD/SOC cuts)

Only a feasible/local point is needed quickly

model.solve(time_limit=...), read the incumbent

solver="gurobi" is a direct matrix backend for LP/MILP/QP/QCP-family models. It intentionally does not translate arbitrary nonlinear expression DAGs into Gurobi nonlinear constraints. For general NLP/MINLP, use discopt’s global algorithms and select Gurobi only where it is a matrix subsolver, for example solver="amp", milp_solver="gurobi".

solver="scip" hands the whole model, nonlinear expressions included, to the open-source global solver SCIP through PySCIPOpt, whose wheel bundles SCIP, so no license or separate install is needed. The expression DAG is translated in memory, not through an .nl file. Every node is either translated exactly or refused by name: SCIP has no atan, erf, min/max or opaque dm.custom; x**y and xlogx are only translated where their argument is provably positive; and a fractional power or sqrt whose base may be negative gets the row base >= 0, because discopt evaluates those to NaN there.

SCIP’s claims are reported, never certified. A SCIP optimal comes back as status="feasible", and a SCIP infeasible or unbounded comes back as status="unknown", all with gap_certified=False. SCIP’s own status is in result.termination, and its dual bound is kept as information with bound_source="external_solver" and bound_valid=False. This is measured, not cautious: on MINLPLib ex1252a, SCIP’s presolve certified 131452.97 with a bound above the true optimum 128893.74, on a model whose reference solution is feasible. discopt cannot check a bound it did not prove.

The point is checked: discopt re-verifies it on the declared model and recomputes its objective. An answer that rests on SCIP’s numerical limits is reported without a bound. That means a point at a singularity the box admits (1/x near 0, tan near a pole, log near 0), or a value past SCIP’s numerics/hugeval. The usual cause is an unbounded problem: SCIP stops at x = -1e-9 on min 1/x and calls -1e9 optimal.

Pass SCIP parameters as scip_options={"limits/nodes": 1000, ...}, and use {"display/verblevel": 4} to see SCIP’s log. The route bypasses discopt’s own search, so options that steer it (threads, partitions, strategy, …) are reported as ignored.

The throughline: discopt gives a sound global certificate exactly when the model is bounded and algebraically factorable from supported intrinsics. Keep your formulation inside that envelope and you get a proof; step outside it and you get a point.