Certified-global MINLP with AMP

Certified-global MINLP with AMP#

discopt’s default solver path is fast and handles many nonconvex problems well in practice, but for nonconvex MINLPs (bilinear, signomial, concave objectives, trig, etc.) it cannot in general issue a certificate that the returned solution is globally optimal — it could be stuck in a local basin.

The Adaptive Multivariate Partitioning (AMP) solver path closes that gap. AMP is a piecewise-relaxation method in the lineage of Alpine.jl [Nagarajan et al., 2019], BARON, and Couenne. At every iteration it maintains a valid lower bound $LB_k$ and upper bound $UB_k$ on the global optimum:

\[LB_k \;\le\; \mathtt{global\_opt} \;\le\; UB_k\]

and refines a partition of the nonconvex variables until $\frac{UB_k - LB_k}{|UB_k|} \le \mathtt{rel\_gap}$. At termination the returned result.gap is a hard guarantee on suboptimality.

How AMP works#

Each iteration $k$:

MILP relaxation — Build a piecewise McCormick [McCormick, 1976] (or convex-hull / SOS2) relaxation of every nonconvex term, augmented with $\alpha$-BB cuts [Adjiman et al., 1998] for nonconvex quadratics. Solve with HiGHS → valid lower bound $LB_k$ and a relaxed candidate point.
NLP subproblem — Fix integers and solve the original NLP at the candidate (Ipopt). Any feasible solution gives a valid upper bound $UB_k$ and an incumbent.
Refine — Pick the variable / interval with the largest relaxation gap and split it. Add outer-approximation cuts from the incumbent. Loop.

Setup#

import discopt.modeling as dm

print("discopt loaded")

discopt loaded

Example 1 — A nonconvex separable QP#

Consider the concave-quadratic problem

\[\begin{split} \begin{aligned} \min_{x \in \mathbb{R}^3} \quad & \sum_{i=1}^{3} -(x_i - c_i)^2 \\ \text{s.t.} \quad & -1 \;\le\; \sum_i x_i \;\le\; 3 \\ \quad & x_i \in [-2,\, 2] \end{aligned} \end{split}\]

with offsets $c = (-1,\, 0.5,\, 1.5)$. The objective is concave (the negative-of-square terms each open downward), so the minimum sits at a vertex of the feasible polytope — a setting where a generic NLP solver is not guaranteed to find the global optimum.

Default solver:

c = [-1.0, 0.5, 1.5]


def build_concave_qp():
    m = dm.Model("concave_qp")
    xs = [m.continuous(f"x{i}", lb=-2.0, ub=2.0) for i in range(3)]
    m.subject_to(sum(xs) >= -1.0)
    m.subject_to(sum(xs) <= 3.0)
    m.minimize(sum(-((xs[i] - c[i]) ** 2) for i in range(3)))
    return m, xs


m, xs = build_concave_qp()
result_default = m.solve()
print("Default solver:")
print(f"  status    = {result_default.status}")
print(f"  objective = {result_default.objective:.6f}")
print(
    f"  x         = ({result_default.x['x0']:.4f}, "
    f"{result_default.x['x1']:.4f}, {result_default.x['x2']:.4f})"
)

Default solver:
  status    = optimal
  objective = -23.500000
  x         = (2.0000, 2.0000, -2.0000)

Now solve with AMP — note the explicit gap field on the result:

m, xs = build_concave_qp()
result_amp = m.solve(solver="amp", rel_gap=1e-4)
print("AMP (certified global):")
print(f"  status    = {result_amp.status}")
print(f"  objective = {result_amp.objective:.6f}")
print(f"  gap       = {result_amp.gap:.2e}     <-- certified suboptimality bound")
print(f"  MIPs      = {result_amp.mip_count}")
print(f"  wall time = {result_amp.wall_time:.2f}s")
print(
    f"  x         = ({result_amp.x['x0']:.4f}, {result_amp.x['x1']:.4f}, {result_amp.x['x2']:.4f})"
)

AMP (certified global):
  status    = optimal
  objective = -23.500000
  gap       = 6.81e-05     <-- certified suboptimality bound
  MIPs      = 3
  wall time = 0.05s
  x         = (2.0000, 2.0000, -2.0000)

Both solvers find the same vertex, but AMP additionally returns a gap that proves the answer is within rel_gap of the true global optimum. That certificate is what AMP buys you.

Example 2 — A nonconvex MINLP#

Now add binary on/off variables. We pick 2 of 3 “facilities” and choose continuous loadings on each picked facility to minimize a concave operating cost subject to a minimum total throughput:

\[\begin{split} \begin{aligned} \min_{x,z} \quad & \sum_{i=1}^{3} -(x_i - c_i)^2 \\ \text{s.t.} \quad & \sum_i x_i \;\ge\; 2 \\ & \sum_i z_i \;=\; 2 \\ & x_i \le 4\, z_i, \quad x_i \in [0,4], \quad z_i \in \{0,1\} \end{aligned} \end{split}\]

c = [1.0, 2.5, 3.5]

m = dm.Model("facility")
xs = [m.continuous(f"x{i}", lb=0.0, ub=4.0) for i in range(3)]
zs = [m.binary(f"z{i}") for i in range(3)]

m.subject_to(sum(xs) >= 2.0)
m.subject_to(sum(zs) == 2)
for i in range(3):
    m.subject_to(xs[i] <= 4.0 * zs[i])

m.minimize(sum(-((xs[i] - c[i]) ** 2) for i in range(3)))

result = m.solve(solver="amp", rel_gap=1e-4, max_iter=100)
print(f"status    = {result.status}")
print(f"objective = {result.objective:.6f}")
print(f"gap       = {result.gap:.2e}")
print(f"MIPs      = {result.mip_count}")
print(f"wall time = {result.wall_time:.2f}s")
for i in range(3):
    z = int(round(float(result.x[f"z{i}"])))
    x = float(result.x[f"x{i}"])
    print(f"  z{i}={z:d}  x{i}={x:.4f}")

status    = optimal
objective = -27.500000
gap       = 0.00e+00
MIPs      = 1
wall time = 0.30s
  z0=1  x0=4.0000
  z1=0  x1=-0.0000
  z2=1  x2=-0.0000

Tuning AMP#

All AMP-specific options are passed as keyword arguments to Model.solve(solver="amp", ...). The most useful:

Option	Default	What it does
`rel_gap`	`1e-4`	Stop when $\tfrac{UB - LB}{
`abs_tol`	`1e-6`	Absolute optimality tolerance
`max_iter`	`100`	Hard cap on partition-refinement iterations
`n_init_partitions`	`4`	Initial partitions per discretized variable
`convhull_formulation`	`"disaggregated"`	`"sos2"` and `"facet"` give tighter relaxations at higher MILP cost
`convhull_ebd`	`False`	Logarithmic Gray-code embedded SOS2 binaries (fewer integer vars per partition)
`convhull_ebd_encoding`	`"gray"`	`"gray"` is the only SOS2-compatible choice for >2 partitions
`presolve_bt`	`True`	OBBT + FBBT bound tightening before the first MILP
`presolve_bt_algo`	`"obbt"`	`"obbt"`, `"fbbt"`, or `"both"`
`obbt_at_root`	`True`	Strengthen variable bounds at the root using LP-based OBBT
`obbt_with_cutoff`	`True`	Use the current incumbent as a cutoff during OBBT
`alphabb_cutoff_obbt`	`True`	Add $\alpha$-BB OA cuts during cutoff OBBT
`partition_method`	`"adaptive"`	How to choose the next refinement variable / interval
`disc_var_pick`	`"max_cover"`	Variable-selection heuristic for partitioning
`use_start_as_incumbent`	`False`	Trust a user-supplied `initial_solution` as a starting incumbent
`milp_time_limit`	`None`	Per-MILP wall-clock cap (seconds)
`milp_gap_tolerance`	`1e-4`	MILP-solver gap stop

Comparing convex-hull formulations#

The default "disaggregated" McCormick is the simplest; "sos2" builds the convex hull of each piecewise segment using SOS2 lambda variables (tighter but more binaries); "sos2" with convhull_ebd=True swaps the per-interval binaries for a logarithmic Gray-code encoding.

configs = [
    ("disaggregated", False),
    ("sos2", False),
    ("sos2", True),
]
for formulation, ebd in configs:
    m, xs = build_concave_qp()
    result = m.solve(
        solver="amp",
        rel_gap=1e-4,
        convhull_formulation=formulation,
        convhull_ebd=ebd,
    )
    label = f"{formulation}{' + embedded' if ebd else ''}"
    print(
        f"{label:25s} obj={result.objective:.6f}  gap={result.gap:.2e}  "
        f"MIPs={result.mip_count}  time={result.wall_time:.2f}s"
    )

disaggregated             obj=-43.500000  gap=3.68e-05  MIPs=3  time=0.04s
sos2                      obj=-43.500000  gap=3.68e-05  MIPs=3  time=0.04s
sos2 + embedded           obj=-43.500000  gap=3.68e-05  MIPs=3  time=0.04s

Watching convergence with `iteration_callback`#

Pass a callback that receives {iteration, lower_bound, upper_bound} after each iteration to inspect the bound trajectory:

trace = []


def cb(info):
    trace.append((info["iteration"], info["lower_bound"], info["upper_bound"]))


m, xs = build_concave_qp()
_ = m.solve(solver="amp", rel_gap=1e-4, iteration_callback=cb)

print(f"{'iter':>4} {'LB':>14} {'UB':>14} {'gap':>10}")
for it, lb, ub in trace:
    gap = float("inf") if abs(ub) < 1e-12 else abs(ub - lb) / max(1.0, abs(ub))
    print(f"{it:>4} {lb:>14.6f} {ub:>14.6f} {gap:>10.2e}")

iter             LB             UB        gap
   -44.500000     -43.500000   2.30e-02
   -43.540000     -43.500000   9.20e-04
   -43.501600     -43.500000   3.68e-05

When not to use AMP#

Convex MINLPs. The default solver (outer approximation / convex B&B) is faster and equally certified for convex problems — AMP’s relaxation overhead is wasted.
Pure LP / MILP / convex QP. Use the dedicated solvers; AMP adds no value.
Very large nonconvex problems where a fast feasible solution matters more than a certificate. A local NLP from a good warm start may return a usable point faster, at the cost of no guarantee.
Problems with operator combinations AMP cannot factorize cleanly. AMP requires the nonlinearities to decompose into bilinear / trilinear / monomial / supported-univariate (exp, log, sqrt, sin, cos, tan, abs) terms. Rich polynomial objectives like $(x^2-1)^2 + xy$ that contain mixed repeated-factor products may fall back to a feasibility relaxation and lose the certificate.

Use AMP when you have a small-to-medium nonconvex (MI)NLP, the structure is factorable, and you need a certified-global answer or a verified bound on suboptimality.

References#

AMP’s algorithm follows the Alpine.jl design [Nagarajan et al., 2019]; the piecewise-McCormick relaxations originate with McCormick [1976]; the $\alpha$-BB cuts for nonconvex quadratics follow Adjiman et al. [1998].