Solver Internals: BARON-Parity Features

Solver Internals: BARON-Parity Features#

This notebook documents the algorithmic features that close the gap between discopt and commercial global solvers like BARON [Tawarmalani and Sahinidis, 2005]. These are internal solver improvements that users benefit from automatically, but understanding them helps with tuning.

Topics covered:

  1. Convexity detection — automatic classification of expressions, constraints, and models

  2. Adaptive solver dispatch — skipping unnecessary relaxation overhead for convex problems

  3. Per-constraint OA cuts — outer approximation for individually convex constraints

  4. Reliability branching — pseudocost-guided variable selection in Branch & Bound

  5. OBBT with incumbent cutoff — tighter bounds using the best known solution

  6. FBBT with incumbent cutoff — fast constraint propagation with objective cutoff

  7. Specialized envelopes — tight relaxations for trig, signomial, and additional functions

import os

os.environ["JAX_PLATFORMS"] = "cpu"
os.environ["JAX_ENABLE_X64"] = "1"


import discopt.modeling as dm
from discopt._relax.convexity import classify_constraint, classify_expr, classify_model

1. Convexity Detection#

discopt includes a convexity detector that walks the expression DAG and classifies each subexpression using standard composition rules from convex analysis [Boyd and Vandenberghe, 2004]:

Rule

Result

constant, variable, parameter

AFFINE

convex + convex

CONVEX

concave + concave

CONCAVE

convex + concave

UNKNOWN

positive_const * convex

CONVEX

negative_const * convex

CONCAVE

exp(convex)

CONVEX (exp is convex & nondecreasing)

log(concave)

CONCAVE (log is concave & nondecreasing)

x^2 (affine base)

CONVEX

x^(2k) (even power, affine base)

CONVEX

x^p with p > 1, x >= 0

CONVEX

x^p with 0 < p < 1, x >= 0

CONCAVE

x * y (bilinear)

UNKNOWN

Results are cached on expression IDs for efficiency.

# Build some expressions and classify them
m = dm.Model("convexity_demo")
x = m.continuous("x", lb=0, ub=5)
y = m.continuous("y", lb=0, ub=5)

expressions = {
    "x": x,
    "x + y": x + y,
    "x^2": x**2,
    "x^2 + y^2": x**2 + y**2,
    "exp(x)": dm.exp(x),
    "log(x)": dm.log(x),
    "sqrt(x)": dm.sqrt(x),
    "exp(x^2)": dm.exp(x**2),
    "log(x + y)": dm.log(x + y),
    "x * y": x * y,
    "-x^2": -(x**2),
    "x^3 (nonneg)": x**3,
}

print(f"{'Expression':<20s} {'Curvature':<12s}")
print("-" * 32)
for name, expr in expressions.items():
    curv = classify_expr(expr, model=m)
    print(f"{name:<20s} {curv.value:<12s}")

Expression           Curvature   
--------------------------------
x                    affine      
x + y                affine      
x^2                  convex      
x^2 + y^2            convex      
exp(x)               convex      
log(x)               unknown     
sqrt(x)              concave     
exp(x^2)             convex      
log(x + y)           unknown     
x * y                unknown     
-x^2                 concave     
x^3 (nonneg)         convex

Constraint Convexity#

A constraint is convex when it defines a convex feasible set:

  • g(x) <= 0 is convex when g is convex

  • g(x) >= 0 is convex when g is concave (equivalently, -g is convex)

  • g(x) == 0 is convex only when g is affine

from discopt.modeling.core import Constraint

constraints = [
    ("x^2 + y^2 <= 10", Constraint(body=x**2 + y**2, sense="<=", rhs=10)),
    ("log(x) >= 1", Constraint(body=dm.log(x), sense=">=", rhs=1)),
    ("x + y == 3", Constraint(body=x + y, sense="==", rhs=3)),
    ("x * y <= 5", Constraint(body=x * y, sense="<=", rhs=5)),
    ("x^2 == 4", Constraint(body=x**2, sense="==", rhs=4)),
]

print(f"{'Constraint':<20s} {'Convex?':<8s}")
print("-" * 28)
for name, c in constraints:
    is_cvx = classify_constraint(c, model=m)
    print(f"{name:<20s} {str(is_cvx):<8s}")

Constraint           Convex? 
----------------------------
x^2 + y^2 <= 10      True    
log(x) >= 1          False   
x + y == 3           True    
x * y <= 5           False   
x^2 == 4             False

Model-Level Classification#

classify_model() checks whether the entire optimization problem is convex — meaning the objective is convex (for minimization) and all constraints are convex.

# A convex QP
m1 = dm.Model("convex_qp")
x = m1.continuous("x", lb=-5, ub=5)
y = m1.continuous("y", lb=-5, ub=5)
m1.minimize(x**2 + y**2)
m1.subject_to(x + y >= 1)

is_cvx, mask = classify_model(m1)
print(f"Convex QP: model_is_convex={is_cvx}, constraint_mask={mask}")

# A nonconvex problem (bilinear constraint)
m2 = dm.Model("nonconvex")
x = m2.continuous("x", lb=0, ub=5)
y = m2.continuous("y", lb=0, ub=5)
m2.minimize(x + y)
m2.subject_to(x * y >= 1)  # nonconvex
m2.subject_to(x + y <= 8)  # convex (affine)

is_cvx, mask = classify_model(m2)
print(f"Nonconvex:  model_is_convex={is_cvx}, constraint_mask={mask}")
print("  -> constraint 0 (x*y >= 1) is nonconvex, constraint 1 (x+y <= 8) is convex")

Convex QP: model_is_convex=True, constraint_mask=[True]
Nonconvex:  model_is_convex=False, constraint_mask=[False, True]
  -> constraint 0 (x*y >= 1) is nonconvex, constraint 1 (x+y <= 8) is convex

For a mixed-integer convex model this goes one step further: the model is routed to the MIP–NLP family and certified there, not in discopt’s spatial tree. The default method is LP/NLP-based branch and bound [Quesada and Grossmann, 1992] — single-tree outer approximation [Duran and Grossmann, 1986]: one MILP master (HiGHS) whose integer incumbents trigger fixed-integer NLP solves and new OA cuts. result.node_count is therefore the master’s node count, which is small (and is 0 when HiGHS’s presolve settles the master), and the log names the route rather than a backend. The cell below turns INFO logging on to show exactly that, and turns it back off afterwards — without the reset, every later cell in this notebook inherits it and buries its own output under the iteration log.

import logging
import time

# Enable logging to see the route the solver picks. Scoped to this cell: the
# handler and both levels are restored below, so the rest of the notebook is not
# buried under the OA iteration log.
logging.basicConfig(level=logging.INFO, format="%(name)s: %(message)s")
_root_level = logging.getLogger().level
_solver_logger = logging.getLogger("discopt.solver")
_oa_logger = logging.getLogger("discopt.solvers.oa")
_prev_levels = (_solver_logger.level, _oa_logger.level)
_solver_logger.setLevel(logging.INFO)
_oa_logger.setLevel(logging.INFO)

# Solve a convex MIQP — the solver detects convexity and skips relaxation overhead
m = dm.Model("convex_miqp")
x = m.continuous("x", lb=0, ub=10)
y = m.continuous("y", lb=0, ub=10)
z = m.binary("z")

m.minimize(x**2 + y**2 + z)
m.subject_to(x + y >= 2)
m.subject_to(x <= 8 * z)

t0 = time.perf_counter()
result = m.solve()
elapsed = time.perf_counter() - t0

print(f"\nStatus: {result.status}")
print(f"Objective: {result.objective:.6f}")
print(f"Time: {elapsed:.3f}s")
print(f"Nodes: {result.node_count}")
print(f"Route: {result.algorithm_route}")

# Restore, so the remaining cells show their own output and nothing else.
_solver_logger.setLevel(_prev_levels[0])
_oa_logger.setLevel(_prev_levels[1])
logging.getLogger().setLevel(logging.WARNING)

assert result.gap_certified, "a convex MIQP is certified by the MIP-NLP route"
assert "certified convex at the root" in (result.algorithm_route or ""), result.algorithm_route
assert "spatial" not in (result.algorithm_route or ""), result.algorithm_route

discopt.solver: Convex MINLP auto-route: mip-nlp/lp_nlp_bb: miqp certified convex at the root (DISCOPT_CONVEX_MINLP_ROUTE; master=highs)

discopt.solvers.oa: LP/NLP BB: 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.000000
Time: 0.514s
Nodes: 2
Route: mip-nlp/lp_nlp_bb: miqp certified convex at the root (DISCOPT_CONVEX_MINLP_ROUTE; master=highs)

3. Per-Constraint Outer Approximation#

Outer approximation (OA) cuts [Duran and Grossmann, 1986, Fletcher and Leyffer, 1994] are tangent hyperplanes that provide valid linear underestimators of convex constraints. In a nonconvex problem, some constraints may still be individually convex — OA cuts are valid for those constraints even though the overall problem is nonconvex.

discopt uses the convexity detector to identify which constraints are convex, then generates OA cuts only for those. This replaces the earlier approach that only generated OA for affine (linear) constraints.

# Mixed convex/nonconvex problem
m = dm.Model("per_constraint_oa")
x = m.continuous("x", lb=0.1, ub=5)
y = m.continuous("y", lb=0.1, ub=5)
z = m.binary("z")

m.minimize(x**2 + y + z)
m.subject_to(x**2 + y**2 <= 10)  # convex constraint -> OA valid
m.subject_to(x * y >= 1)  # nonconvex constraint -> no OA
m.subject_to(x + y <= 6 * z)  # affine -> OA valid

is_cvx, mask = classify_model(m)
print(f"Model is convex: {is_cvx}")
print(f"Per-constraint convexity: {mask}")
print("  constraint 0 (x^2+y^2 <= 10): convex =", mask[0])
print("  constraint 1 (x*y >= 1):       convex =", mask[1])
print("  constraint 2 (x+y <= 6z):      convex =", mask[2])
print()

# Solve with cutting planes — OA cuts generated for constraints 0 and 2
result = m.solve(cutting_planes=True, max_nodes=200)
print(f"Status: {result.status}, Objective: {result.objective:.4f}")

Model is convex: False
Per-constraint convexity: [True, False, True]
  constraint 0 (x^2+y^2 <= 10): convex = True
  constraint 1 (x*y >= 1):       convex = False
  constraint 2 (x+y <= 6z):      convex = True

Status: optimal, Objective: 2.8899

4. Reliability Branching#

In Branch & Bound, variable selection significantly impacts performance. discopt implements reliability branching [Achterberg et al., 2005], which combines:

  • Pseudocost branching: After each branch, record the objective improvement per unit of variable change. Maintain running averages (pseudocosts) per variable. Select the variable with the highest product score.

  • Strong branching: For variables with fewer than \(\eta_{\text{rel}}\) (default 8) observations, solve LP relaxations for both child nodes to get reliable bound estimates.

The hybrid approach starts with strong branching (expensive but informative) and transitions to pseudocost branching (cheap and accurate) as the tree matures.

How it works in discopt#

  1. The Rust TreeManager tracks pseudocosts and observation counts per integer variable

  2. score_candidates(solution) returns fractional variables with their pseudocost scores and observation counts

  3. For unreliable variables (obs < threshold), Python-side strong branching solves LP relaxations via HiGHS

  4. Branch hints are passed back to Rust via set_branch_hints()

# Reliability branching on a larger MINLP (built in, no configuration needed)
m = dm.Model("branching_demo")
xs = [m.continuous(f"x{i}", lb=0, ub=5) for i in range(4)]
ys = [m.binary(f"y{i}") for i in range(4)]

# Quadratic objective with indicator constraints
m.minimize(sum((x - 1.5) ** 2 for x in xs) + sum(2 * y for y in ys))
for i in range(4):
    m.subject_to(xs[i] <= 4 * ys[i])
m.subject_to(sum(xs) >= 3)
m.subject_to(sum(ys) <= 3)

result = m.solve(max_nodes=500)
print(f"Status: {result.status}, obj={result.objective:.4f}, nodes={result.node_count}")

Status: optimal, obj=8.2500, nodes=33

5. OBBT with Incumbent Cutoff#

Optimality-Based Bound Tightening (OBBT) [Gleixner et al., 2017] solves LP relaxations to find the tightest possible bounds on each variable:

\[\underline{x}_i = \min_{x} x_i \quad \text{s.t.} \quad Ax \leq b, \quad l \leq x \leq u\]
\[\overline{x}_i = \max_{x} x_i \quad \text{s.t.} \quad Ax \leq b, \quad l \leq x \leq u\]

When an incumbent solution with objective \(z^*\) is known, we can add the constraint \(c^T x \leq z^*\) to the OBBT LPs. This often dramatically tightens bounds because many variable values become infeasible when we require the objective to be at least as good as the incumbent.

discopt runs OBBT periodically during Branch & Bound whenever a new incumbent is found.

from discopt._relax.obbt import run_obbt

# Build a model with linear constraints
m = dm.Model("obbt_demo")
x = m.continuous("x", lb=0, ub=10)
y = m.continuous("y", lb=0, ub=10)

m.minimize(x + 2 * y)
m.subject_to(x + y >= 3)
m.subject_to(2 * x + y <= 12)

# Run OBBT without cutoff
result_no_cut = run_obbt(m)
print("OBBT without incumbent cutoff:")
print(f"  x bounds: [{result_no_cut.tightened_lb[0]:.4f}, {result_no_cut.tightened_ub[0]:.4f}]")
print(f"  y bounds: [{result_no_cut.tightened_lb[1]:.4f}, {result_no_cut.tightened_ub[1]:.4f}]")
print(f"  LPs solved: {result_no_cut.n_lp_solves}, bounds tightened: {result_no_cut.n_tightened}")
print()

# Run OBBT with incumbent cutoff z* = 8 (x + 2y <= 8)
result_cut = run_obbt(m, incumbent_cutoff=8.0)
print("OBBT with incumbent cutoff z* = 8:")
print(f"  x bounds: [{result_cut.tightened_lb[0]:.4f}, {result_cut.tightened_ub[0]:.4f}]")
print(f"  y bounds: [{result_cut.tightened_lb[1]:.4f}, {result_cut.tightened_ub[1]:.4f}]")
print(f"  LPs solved: {result_cut.n_lp_solves}, bounds tightened: {result_cut.n_tightened}")

OBBT without incumbent cutoff:
  x bounds: [0.0000, 6.0000]
  y bounds: [0.0000, 10.0000]
  LPs solved: 4, bounds tightened: 1

OBBT with incumbent cutoff z* = 8:
  x bounds: [0.0000, 6.0000]
  y bounds: [0.0000, 4.0000]
  LPs solved: 4, bounds tightened: 2

6. FBBT with Incumbent Cutoff#

Feasibility-Based Bound Tightening (FBBT) [Belotti et al., 2009] propagates bound information through the expression DAG without solving optimization problems, making it much cheaper than OBBT. When an incumbent solution with objective \(z^*\) is known, the cutoff constraint \(f(x) \leq z^*\) is added to the propagation, often triggering cascading domain reductions.

In discopt, FBBT is implemented in Rust (fbbt_with_cutoff) and exposed to Python via PyO3 bindings. It runs at every B&B node before the NLP relaxation solve, while OBBT is reserved for periodic invocation when a new incumbent is found. The combination provides complementary bound tightening: FBBT is fast but conservative; OBBT is expensive but finds the tightest LP-based bounds [Ryoo and Sahinidis, 1996].

import numpy as np
from discopt._rust import model_to_repr

# Build a simple model and run FBBT with and without cutoff
m = dm.Model("fbbt_demo")
x = m.continuous("x", lb=0, ub=100)
y = m.continuous("y", lb=0, ub=100)
m.minimize(x)
m.subject_to(x + y <= 10)

repr_ = model_to_repr(m)

# Without cutoff: FBBT propagates constraint x+y <= 10
lbs0, ubs0 = repr_.fbbt(max_iter=10, tol=1e-8)
lbs0, ubs0 = np.asarray(lbs0), np.asarray(ubs0)
print("FBBT without cutoff:")
print(f"  x: [0, {ubs0[0]:.0f}]  (was [0, 100])")
print(f"  y: [0, {ubs0[1]:.0f}]  (was [0, 100])")

# With cutoff z* = 7: adds x <= 7 to propagation
lbs1, ubs1 = repr_.fbbt_with_cutoff(max_iter=10, tol=1e-8, incumbent_bound=7.0)
lbs1, ubs1 = np.asarray(lbs1), np.asarray(ubs1)
print("\nFBBT with incumbent cutoff z* = 7:")
print(f"  x: [0, {ubs1[0]:.0f}]  (was [0, 100])")
print(f"  y: [0, {ubs1[1]:.0f}]  (was [0, 100])")
print(f"\nCutoff tightened x upper bound: {ubs0[0]:.0f} -> {ubs1[0]:.0f}")

FBBT without cutoff:
  x: [0, 10]  (was [0, 100])
  y: [0, 10]  (was [0, 100])

FBBT with incumbent cutoff z* = 7:
  x: [0, 7]  (was [0, 100])
  y: [0, 10]  (was [0, 100])

Cutoff tightened x upper bound: 10 -> 7

7. Univariate Envelopes#

Standard McCormick relaxations [McCormick, 1976] compose relaxations through the expression tree, which can be loose because an outer function sees the relaxation of its argument rather than the argument itself. The default relaxation engine (the numpy factorable engine in discopt._relax.uniform_relax) therefore relaxes each univariate atom — exp, log, x**p, sin, cos, tan, the inverse and hyperbolic functions, erf, x*log(x), softplus, sigmoid, … — with its own convex/concave envelope on the node box, choosing secants and tangents by the atom’s curvature on that box [Locatelli and Schoen, 2014]. Products of positive-variable powers are relaxed in log space [Maranas and Floudas, 1997].

The full inventory, with the soundness rules each envelope follows, is docs/design/relaxation-catalog.md.

Note

Earlier versions of this notebook demonstrated a separate JAX envelope library (relax_power_int, relax_sin_tight, …). It fed only the McCormick-NLP relaxation route, which was retired in #1873; its functions are gone.

The solver applies these automatically. A convex model needs none of them: below, discopt detects convexity and solves the model by outer approximation.

# A convex model: convexity is detected and the solve goes to outer approximation
m = dm.Model("envelope_integration")
x = m.continuous("x", lb=0.5, ub=4)
y = m.continuous("y", lb=0.5, ub=4)
z = m.binary("z")

m.minimize(x**3 + dm.exp(y) + z)
m.subject_to(x + y >= 2)
m.subject_to(x <= 3 * z)

result = m.solve(max_nodes=200)
print(f"Status: {result.status}")
print(f"Objective: {result.objective:.4f}")
print(f"x={result.x['x'].item():.4f}, y={result.x['y'].item():.4f}, z={result.x['z'].item():.0f}")

Status: optimal
Objective: 4.7137
x=0.9675, y=1.0325, z=1

8. PSD and SOC cuts for QCQP#

For quadratically constrained quadratic programs (QCQP), term-wise McCormick / RLT relaxations lift each product \(x_i x_j\) to an independent column \(X_{ij}\) but never enforce two structural facts that hold at every feasible point. discopt closes that gap with two on-demand separators that need no SDP or conic solver — both emit ordinary linear cuts [Belotti et al., 2013].

PSD (eigenvalue / moment) cuts#

The moment matrix

\[\begin{split}M(x, X) = \begin{bmatrix} 1 & x^\top \\ x & X \end{bmatrix}\end{split}\]

is positive semidefinite at every true point (there \(X = x x^\top\), so \(M = [1; x][1; x]^\top \succeq 0\)). The term-wise relaxation does not impose \(M \succeq 0\), and that missing condition is the dominant source of the QCQP gap. At a relaxation point \((x^*, X^*)\) whose moment submatrix has a negative eigenvalue \(\lambda_{\min}\) with eigenvector \(v\), the inequality

\[v^\top M(x, X)\, v \ge 0\]

is linear in \((x, X)\) and valid for the whole feasible region (it equals \((v_0 + v_{\text{rest}}^\top x)^2 \ge 0\) at any true point), while being violated at \((x^*, X^*)\). discopt separates it with a single dense eigh on the moment submatrix, over cliques of variables whose pairwise products and squares are all lifted — dense \(k\ge 3\) cliques capture multi-variable moment coupling that pairwise \(2\times2\) minors miss. Opt in with m.solve(psd_cuts=True).

SOC (second-order-cone) cuts#

A convex cone constraint \(\lVert y \rVert_2 \le t\) is not polyhedral, so an LP relaxation cannot represent it exactly. It can be outer-approximated to any accuracy by gradient cuts [Ben-Tal and Nemirovski, 2001]: at a point \((y^*, t^*)\) outside the cone (\(\lVert y^* \rVert > t^*\)), the unit vector \(u = y^*/\lVert y^* \rVert\) gives the valid inequality \(u^\top y \le t\). These tighten the convex-cone substructure that conic-representable terms expose (e.g. Euclidean-norm sqrt(xᵀQx) patterns), complementing the PSD cuts that target the nonconvex moment matrix.

A QCQP relaxation-strength scoreboard (discopt_benchmarks/benchmarks/qcqp_scoreboard.py) A/B’s these cut families against the baseline on registered nonconvex QCQP instances, tracking the root dual-bound gap, node count, and a hard correctness gate (the solved objective must match the known optimum). The AC optimal power flow problem is the capstone target for this machinery — its builder now lives in the standalone discopt-apps plugin (#431); see the AC-OPF page.

# A small nonconvex QCQP: minimize x + y over a nonconvex (bilinear) feasible set.
# The feasible region x*y >= 4 is nonconvex, so the moment matrix of the McCormick
# relaxation can be indefinite at LP points -> PSD cuts have something to separate.
m = dm.Model("qcqp_psd_demo")
x = m.continuous("x", lb=0, ub=4)
y = m.continuous("y", lb=0, ub=4)

m.minimize(x + y)
m.subject_to(x * y >= 4)  # nonconvex quadratic constraint
m.subject_to(x + 2 * y <= 8)

# Baseline (auto cut policy) vs. PSD eigenvalue cuts engaged. mccormick_bounds="lp"
# gives a valid global lower bound; subnlp_frequency=1 recovers incumbents from LP
# primals. Both must certify the same global optimum.
r_base = m.solve(mccormick_bounds="lp", subnlp_frequency=1, max_nodes=3000)
print("Baseline:")
print(f"  Status: {r_base.status}, Obj: {r_base.objective:.6f}, Nodes: {r_base.node_count}")

r_psd = m.solve(psd_cuts=True, mccormick_bounds="lp", subnlp_frequency=1, max_nodes=3000)
print("\nWith PSD (moment / eigenvalue) cuts:")
print(f"  Status: {r_psd.status}, Obj: {r_psd.objective:.6f}, Nodes: {r_psd.node_count}")

assert abs(r_base.objective - r_psd.objective) < 1e-4
print("\nSame certified optimum: PSD cuts tighten the relaxation without changing it.")

Baseline:
  Status: optimal, Obj: 4.000000, Nodes: 5

With PSD (moment / eigenvalue) cuts:
  Status: optimal, Obj: 4.000000, Nodes: 5

Same certified optimum: PSD cuts tighten the relaxation without changing it.

Objective-gating priority branching and best-estimate node selection#

Two search-order refinements help on problems where the dual bound stays pinned until the right integers are fixed [Achterberg et al., 2005].

Objective-gating priority branching (opt-in via the environment variable DISCOPT_OBJ_BRANCH_PRIORITY=1) branches first on the integer/binary variables that gate the objective’s nonlinear terms. On some models a product term is structurally zero — and the relaxation bound stuck at zero — until a whole set of selector binaries is jointly fixed. No single-variable fractionality or pseudocost score points at those variables, so most-fractional branching wanders. The priority set is the integers that either appear directly in a nonlinear term or are tied to a nonlinear-term variable by a two-variable linear equality (x_int = c·x_nl). It is branching-order metadata only — it never enters a bound or feasibility test, so it cannot affect soundness.

Best-estimate node selection (m.solve(strategy="best_estimate")) orders open nodes by a pseudocost-based estimate of the subtree optimum rather than by the dual bound alone. This surfaces good incumbents earlier (improving the cutoff for FBBT/OBBT), while strategy="best_first" remains the choice that drives the proof of optimality. The other strategies are "best_first" (default, lowest dual bound) and "depth_first".

import os

# An MINLP whose nonlinear objective terms are gated by selector binaries.
m = dm.Model("obj_gating_demo")
xs = [m.continuous(f"x{i}", lb=0, ub=2) for i in range(3)]
ys = [m.binary(f"y{i}") for i in range(3)]

m.minimize(sum(x * x for x in xs) + sum(ys))  # x_i^2 gated by x_i <= 2*y_i
for i in range(3):
    m.subject_to(xs[i] <= 2 * ys[i])
m.subject_to(sum(xs) >= 2)

# Best-estimate node selection (finds incumbents earlier).
r_be = m.solve(strategy="best_estimate", max_nodes=2000)
print(f"best_estimate: status={r_be.status}, obj={r_be.objective:.4f}, nodes={r_be.node_count}")

# Objective-gating priority branching is opt-in via an environment variable.
os.environ["DISCOPT_OBJ_BRANCH_PRIORITY"] = "1"
try:
    r_pri = m.solve(strategy="best_estimate", max_nodes=2000)
    print(
        f"obj-priority:  status={r_pri.status}, obj={r_pri.objective:.4f}, "
        f"nodes={r_pri.node_count}"
    )
finally:
    os.environ["DISCOPT_OBJ_BRANCH_PRIORITY"] = "0"

assert abs(r_be.objective - r_pri.objective) < 1e-4
print("Same optimum: both are search-order metadata only, never affecting soundness.")

best_estimate: status=optimal, obj=4.0000, nodes=13
obj-priority:  status=optimal, obj=4.0000, nodes=13
Same optimum: both are search-order metadata only, never affecting soundness.

Transcendental relaxation coverage#

The relaxation library has closed several gaps for transcendental and non-smooth functions [McCormick, 1976]. Inverse-trigonometric envelopes for asin, acos, and acosh are now covered, the previously non-smooth abs, min, and max relaxations are fixed, and a general-purpose LP relaxer engages for transcendental nonlinearity that no specialized envelope matches — so a model mixing these functions still gets a valid convex relaxation rather than falling back to a feasibility objective.

Summary: Solver Options for Tuning#

All these features are enabled automatically by default. Key parameters for tuning:

Parameter

Default

Description

cutting_planes

True

Enable RLT + per-constraint OA cuts

mccormick_bounds

"auto"

"auto", "lp" (valid global LB; recommended for nonconvex + bilinear, pair with subnlp_frequency=1), "none"

max_nodes

10000

Node limit for B&B tree

gap_tolerance

1e-4

Relative optimality gap tolerance

When to set mccormick_bounds="lp"#

The default "auto" resolves to "lp" for nonconvex models with an objective and to "none" for convex ones (whose NLP relaxation is already a valid bound). For pure-continuous nonconvex problems with bilinear/multilinear terms (pooling, polynomial NLP, QCQP), set mccormick_bounds="lp" and subnlp_frequency=1: this builds the full McCormick LP reformulation (a valid global LB), and the SubNLP heuristic turns each LP primal into a feasible incumbent. On a 22-instance MINLPLib bilinear+continuous sweep this delivered 13 strict wins (5 better objectives, 8 speedups to certified global optimum) over the default — at the cost of 2 instances that lose feasibility when the LP relaxation is loose. Stay with the default for convex or pure-integer models.

The convexity detector, adaptive dispatch, OBBT with incumbent cutoff, and the univariate envelopes all operate automatically without user configuration.

(mccormick_bounds="nlp", "midpoint" and partitions were removed: "midpoint" was not a valid bound (correctness issue C-18), and the McCormick-NLP route behind the other two was retired in #1873. Each is refused with a message saying why.)

Conflict analysis and no-good cuts#

When a partial assignment of binary variables is jointly infeasible, it can be excluded tree-wide with a no-good cut

\[ \sum_{j:\,a_j=1}(1-x_j) + \sum_{j:\,a_j=0} x_j \;\ge\; 1, \]

which removes exactly that assignment and nothing else. discopt.conflict uses feasibility-based bound tightening (FBBT) as the infeasibility oracle — FBBT never reports a feasible subproblem as infeasible, so every emitted cut is sound and never removes the global optimum. The detector keeps only minimal conflicts (it skips assignments already excluded by a shorter one).

import os
os.environ["JAX_PLATFORMS"] = "cpu"
os.environ["JAX_ENABLE_X64"] = "1"

import discopt.modeling.core as dm
from discopt.conflict import find_conflict_cuts, add_conflict_cuts

# A small packing model: at most one of three binaries may be 1.
m = dm.Model("conflict_demo")
x = [m.binary(f"x{i}") for i in range(3)]
m.maximize(x[0] + 2 * x[1] + 3 * x[2])
m.subject_to(x[0] + x[1] + x[2] <= 1)

# Each pair set to (1, 1) is jointly infeasible -> three minimal no-goods.
cuts = find_conflict_cuts(m, max_order=2)
print("minimal conflict cuts found:", len(cuts))

n_added = add_conflict_cuts(m, max_order=2)
res = m.solve()
print(f"added {n_added} no-good cuts; optimum still {float(res.objective):.1f}")

minimal conflict cuts found: 3

added 3 no-good cuts; optimum still 3.0