discopt Advanced Features#
This notebook showcases the advanced capabilities that distinguish discopt from traditional MINLP solvers:
Differentiable optimization — Compute sensitivities of optimal solutions to parameters
Where JAX still lives — autodiff and differentiable layers, now that the solve path is Rust
Piecewise McCormick relaxations — Tighter convex underestimators for nonlinear terms
Cutting planes — RLT and outer approximation cuts
Multi-start heuristics — Feasibility pump and diverse starting points
CUTEst interface — Loading standard NLP test problems
Prerequisites: familiar with the basics from notebooks/quickstart.ipynb.
import os
os.environ["JAX_PLATFORMS"] = "cpu"
os.environ["JAX_ENABLE_X64"] = "1"
import discopt.modeling as dm
import numpy as np
print("discopt loaded successfully")
discopt loaded successfully
1. Differentiable Optimization#
discopt can differentiate the optimal objective value with respect to model parameters. This is useful for sensitivity analysis, bilevel optimization, and end-to-end learning.
Level 1: Envelope Theorem#
For a parametric NLP \(\min_x f(x; p)\) subject to \(g(x; p) \leq 0\), the sensitivity of the optimal value is given by the envelope theorem [Fiacco, 1983]:
where \(\lambda^*\) are the optimal dual variables (Lagrange multipliers) returned by Ipopt. No additional linear system solve is needed.
Example: how does optimal cost change with the price parameter?#
from discopt._relax.differentiable import differentiable_solve
m = dm.Model("parametric_cost")
x = m.continuous("x", lb=0, ub=5)
y = m.continuous("y", lb=0, ub=5)
price = m.parameter("price", value=3.0)
m.minimize(price * x + y**2)
m.subject_to(x + y >= 2)
result = differentiable_solve(m)
print(f"Status: {result.status}")
print(f"Objective: {result.objective:.6f}")
print(f"x* = {result.x['x'].item():.6f}")
print(f"y* = {result.x['y'].item():.6f}")
print(f"d(obj*)/d(price) = {result.gradient(price):.6f}")
print()
print("Interpretation: increasing price by 1 changes optimal cost by ~d(obj*)/d(price)")
Status: optimal
Objective: 3.750000
x* = 0.500000
y* = 1.500000
d(obj*)/d(price) = 0.500000
Interpretation: increasing price by 1 changes optimal cost by ~d(obj*)/d(price)
Tracing sensitivity across a range of parameter values#
Let’s sweep the price parameter and see how the optimal objective changes.
prices = np.linspace(0.5, 5.0, 10)
objectives = []
sensitivities = []
for p_val in prices:
price.value = np.float64(p_val)
res = differentiable_solve(m)
objectives.append(res.objective)
sensitivities.append(res.gradient(price))
print(f"{'Price':>8s} {'Obj*':>10s} {'dObj*/dPrice':>14s}")
print("-" * 34)
for p, o, s in zip(prices, objectives, sensitivities):
print(f"{p:8.2f} {o:10.4f} {s:14.6f}")
Price Obj* dObj*/dPrice
----------------------------------
0.50 0.9375 1.750000
1.00 1.7500 1.500000
1.50 2.4375 1.250000
2.00 3.0000 1.000000
2.50 3.4375 0.750000
3.00 3.7500 0.500000
3.50 3.9375 0.250000
4.00 4.0000 0.000047
4.50 4.0000 -0.000000
5.00 4.0000 -0.000000
Level 3: Implicit Differentiation via KKT System#
For more accurate sensitivities (including dx*/dp, not just dobj*/dp), use L3 implicit differentiation through the KKT conditions [Agrawal et al., 2019, Amos and Kolter, 2017].
The L3 result also caches the KKT factorization for sIPOPT-style sensitivity queries [Pirnay et al., 2012] — see the dedicated Sensitivity Analysis notebook for approximate_resolve, dual_sensitivity, reduced_hessian, and sensitivity methods.
from discopt._relax.differentiable import differentiable_solve_l3
# Reset parameter
price.value = np.float64(3.0)
result_l3 = differentiable_solve_l3(m)
print(f"Status: {result_l3.status}")
print(f"Objective: {result_l3.objective:.6f}")
print()
print("L1 gradient (envelope theorem):")
print(f" d(obj*)/d(price) = {result_l3.gradient(price):.6f}")
print()
print("L3 gradient (implicit differentiation):")
print(f" d(obj*)/d(price) = {result_l3.implicit_gradient(price):.6f}")
print()
sm = result_l3.sensitivity_matrix()
if sm is not None:
print(f"Full dx*/dp matrix shape: {sm.shape}")
print(f" dx*/d(price) = {sm[:, 0]}")
Status: optimal
Objective: 3.750000
L1 gradient (envelope theorem):
d(obj*)/d(price) = 0.500000
L3 gradient (implicit differentiation):
d(obj*)/d(price) = 0.500000
Full dx*/dp matrix shape: (2, 1)
dx*/d(price) = [-0.5 0.5]
JAX-native differentiable solve (custom_jvp)#
For embedding optimization inside a larger JAX computation, use _make_jax_differentiable_solve which returns a function compatible with jax.grad.
import jax
from discopt._relax.differentiable import _flatten_params, _make_jax_differentiable_solve
# Build a differentiable solve function
price.value = np.float64(3.0)
solve_fn = _make_jax_differentiable_solve(m)
# Get the flat parameter vector
p_flat = _flatten_params(m)
print(f"Parameter vector: {p_flat}")
# Evaluate and differentiate with jax.grad
obj_star = solve_fn(p_flat)
grad_fn = jax.grad(solve_fn)
grad_val = grad_fn(p_flat)
print(f"Optimal objective via JAX: {float(obj_star):.6f}")
print(f"Gradient via jax.grad: {float(grad_val[0]):.6f}")
Parameter vector: [3.]
Optimal objective via JAX: 3.750000
Gradient via jax.grad: 0.500000
2. Where JAX Still Lives#
discopt’s solve path is Rust, not JAX. The NLP engine is POUNCE (a pure-Rust
Ipopt port), LP/MILP run on the in-house Rust simplex, and the branch-and-bound
tree is managed in Rust. The pure-JAX interior-point method that used to serve as
the NLP backend has been retired; nlp_solver="ipm" is kept only as a
back-compat alias that resolves to POUNCE on the NLP/MINLP path.
What JAX still does:
Differentiable layers —
_make_jax_differentiable_solveabove, which wraps a completed solve in acustom_jvpso the solution can be differentiated inside a larger JAX computation. The derivative information comes from the KKT system at the solved point, not from tracing the solver.Autodiff for model derivatives — gradients, Jacobians and Hessians of the user’s expressions, where the JAX evaluator is the one in use.
Because the solve itself no longer runs in JAX, setting JAX_PLATFORMS=gpu
accelerates only those JAX-backed pieces; it does not move the branch-and-bound
search or the NLP solves onto a GPU.
os.environ["JAX_PLATFORMS"] = "gpu" # affects the JAX autodiff/layer path only
# The `nlp_solver="ipm"` alias and `batch_size` still work; `batch_size` is the
# number of B&B nodes exported from the Rust tree per iteration.
import time
m = dm.Model("batch_demo")
x = m.continuous("x", lb=0, ub=5)
y = m.continuous("y", lb=0, ub=5)
z = m.binary("z")
m.minimize(x**2 + y**2 + z)
m.subject_to(x + y >= 1)
t0 = time.perf_counter()
result_ipm = m.solve(nlp_solver="ipm", batch_size=8)
t_ipm = time.perf_counter() - t0
print('nlp_solver="ipm" (resolves to POUNCE):')
print(f" Status: {result_ipm.status}")
print(f" Objective: {result_ipm.objective}")
print(f" Time: {t_ipm:.3f}s")
print(f" Nodes: {result_ipm.node_count}")
nlp_solver="ipm" (resolves to POUNCE):
Status: optimal
Objective: 0.5000000000348728
Time: 0.037s
Nodes: 1
3. Partitioned McCormick Relaxations (AMP)#
Standard McCormick relaxations [McCormick, 1976] can be loose for highly nonlinear functions. Piecewise McCormick [Bergamini et al., 2005, Castro, 2015] divides a variable’s domain into partitions and relaxes each piece separately, which tightens the relaxation where the partition is fine.
In discopt this lives in the Adaptive Multivariate Partitioning solver, solver="amp" [Nagarajan et al., 2019]. AMP refines the partition only where the incumbent says it matters, rather than splitting every domain uniformly; see the AMP notebook for its options.
Note
solve() used to accept partitions=k. That option configured a separate
McCormick-NLP relaxation route, which was retired in #1873: measured on the in-repo
corpus, its bound never changed a result. partitions and the related options are now
refused with a message saying so. The default spatial branch-and-bound already uses the
LP-form McCormick relaxation for nonconvex models.
# The default spatial B&B vs. AMP's adaptive partitioning on a nonconvex MINLP
m = dm.Model("mccormick_demo")
x = m.continuous("x", lb=0.1, ub=5)
y = m.continuous("y", lb=0.1, ub=5)
z = m.binary("z")
m.minimize(dm.exp(x) * y + z)
m.subject_to(x + y >= 2)
m.subject_to(x * y <= 5 * z + 1)
result_bb = m.solve(time_limit=60)
print("Default spatial B&B (LP-form McCormick):")
print(f" Status: {result_bb.status}, Obj: {result_bb.objective:.6f}, Bound: {result_bb.bound:.6f}")
result_amp = m.solve(solver="amp", time_limit=60)
print("\nAMP (adaptive McCormick partitioning):")
print(f" Status: {result_amp.status}, Obj: {result_amp.objective:.6f}, Bound: {result_amp.bound:.6f}")
# Both certify the same global optimum.
assert result_bb.status == result_amp.status == "optimal"
assert abs(result_bb.objective - result_amp.objective) < 1e-5
AMP: direct OA skips 1 of 2 constraint rows not certified convex: rows=[(1, 'nonconvex_quadratic_alpha_bb_candidate')]
Default spatial B&B (LP-form McCormick):
Status: optimal, Obj: 0.668589, Bound: 0.668589
AMP (adaptive McCormick partitioning):
Status: optimal, Obj: 0.668589, Bound: 0.668589
4. Cutting Planes#
discopt supports two types of cutting planes to tighten relaxations:
RLT (Reformulation-Linearization Technique) [Sherali and Adams, 1990]: McCormick linearization envelopes for bilinear terms \(x_i \cdot x_j\)
OA (Outer Approximation) [Duran and Grossmann, 1986, Fletcher and Leyffer, 1994]: gradient-based tangent hyperplanes at NLP relaxation solutions
Enable with cutting_planes=True. Cuts are accumulated across B&B iterations.
m = dm.Model("cutting_planes_demo")
x1 = m.continuous("x1", lb=0, ub=5)
x2 = m.continuous("x2", lb=0, ub=5)
y = m.binary("y")
m.minimize(x1**2 + x2**2 + y)
m.subject_to(x1 * x2 >= 1) # bilinear term -> RLT cuts apply
m.subject_to(x1 + x2 <= 5 * y)
# Without cutting planes
r1 = m.solve(cutting_planes=False, max_nodes=200)
print("Without cutting planes:")
print(f" Status: {r1.status}, Obj: {r1.objective}, Nodes: {r1.node_count}")
# With cutting planes
r2 = m.solve(cutting_planes=True, max_nodes=200)
print("\nWith cutting planes (RLT + OA):")
print(f" Status: {r2.status}, Obj: {r2.objective}, Nodes: {r2.node_count}")
Without cutting planes:
Status: optimal, Obj: 3.000000009090619, Nodes: 5
With cutting planes (RLT + OA):
Status: optimal, Obj: 3.000000009090619, Nodes: 5
4b. RLT as a first-class solve option#
The Reformulation-Linearization Technique [Sherali and Adams, 1990] is now a first-class solver choice rather than something bundled implicitly under cutting_planes. Level-1 RLT works by multiplying pairs of valid inequalities and linearizing the resulting products in lifted bilinear variables \(X_{ij}\approx x_i x_j\):
A bound factor \((x_j - l_j)\ge 0\) or \((u_j - x_j)\ge 0\) times a linear constraint factor \((b - a^\top x)\ge 0\) gives the product \((b - a^\top x)(x_j - l_j)\ge 0\). Each operand is non-negative on the feasible region, so the product is too — and replacing every \(x_i x_j\) by its lifted column \(X_{ij}\) turns it into a linear valid inequality in \((x, X)\).
Because \(X = x x^\top\) at any true point, the cut is exact there and never removes a feasible point — it only tightens the relaxation toward the convex hull. discopt separates only violated RLT rows per node (the “targeted” part), reusing the same lifted-column map as the PSD cuts rather than enumerating the full Sherali-Adams hierarchy.
Enable it directly with m.solve(rlt_cuts=True). With the default rlt="auto", a structure-gated policy already decides per-node RLT-vs-PSD separation (RLT on QCQP with linear constraints, PSD on box-QP); passing rlt_cuts=True forces RLT-cut separation on regardless.
cuts="auto" is the solver default [Belotti et al., 2013]: it balances bound-tightening against node reduction by selecting at most one cut family per relaxation rather than stacking them. Use cuts="manual" to disable the auto policy and control the psd_cuts / rlt_cuts flags yourself.
# A small nonconvex QCQP: bilinear objective over boxes with linear constraints.
m = dm.Model("rlt_qcqp")
x = m.continuous("x", lb=0, ub=4)
y = m.continuous("y", lb=0, ub=4)
m.minimize(x * y - 3 * x - 2 * y) # bilinear term x*y -> RLT cuts apply
m.subject_to(x + y <= 5)
m.subject_to(x - y <= 2)
# Solve WITHOUT RLT cuts. mccormick_bounds="lp" gives a valid global lower bound
# for this nonconvex bilinear model; subnlp_frequency=1 turns each LP primal into
# a feasible incumbent.
r_off = m.solve(rlt_cuts=False, mccormick_bounds="lp", subnlp_frequency=1, max_nodes=2000)
print("Without RLT cuts:")
print(f" Status: {r_off.status}, Obj: {r_off.objective:.6f}, Nodes: {r_off.node_count}")
# Solve WITH first-class RLT cuts (per-node targeted constraint x bound products).
r_on = m.solve(rlt_cuts=True, mccormick_bounds="lp", subnlp_frequency=1, max_nodes=2000)
print("\nWith RLT cuts (rlt_cuts=True):")
print(f" Status: {r_on.status}, Obj: {r_on.objective:.6f}, Nodes: {r_on.node_count}")
# Both must reach the same certified global optimum.
assert abs(r_off.objective - r_on.objective) < 1e-4
print("\nSame global optimum reached either way (RLT only tightens the relaxation).")
Without RLT cuts:
Status: optimal, Obj: -8.250000, Nodes: 5
With RLT cuts (rlt_cuts=True):
Status: optimal, Obj: -8.250000, Nodes: 5
Same global optimum reached either way (RLT only tightens the relaxation).
5. Multi-Start Heuristics#
For difficult MINLP problems, finding a good initial feasible solution is critical. discopt provides two heuristic strategies:
Multi-start NLP#
Launches NLP solves from diverse starting points generated by stratified random sampling. Tracks the best feasible (and integer-feasible) solution found.
Feasibility Pump#
The feasibility pump [Bonami et al., 2008, Fischetti et al., 2005] iteratively rounds and re-solves to find integer-feasible solutions:
Round integer variables to the nearest integer
Fix integer variables and re-solve the NLP for continuous variables
If feasible, return. Otherwise perturb and retry.
from discopt._relax.primal_heuristics import MultiStartNLP
m = dm.Model("multistart_demo")
x = m.continuous("x", lb=-5, ub=5)
y = m.continuous("y", lb=-5, ub=5)
# Rosenbrock: notoriously hard for single-start solvers due to banana valley
m.minimize((1 - x) ** 2 + 100 * (y - x**2) ** 2)
ms = MultiStartNLP(m, n_starts=16, seed=42)
ms_result = ms.solve()
print("Multi-start results:")
print(f" Starts launched: {ms_result.n_starts}")
print(f" Feasible solutions: {ms_result.n_feasible}")
print(f" Best objective: {ms_result.best_objective:.8f} (expected: 0.0)")
if ms_result.all_objectives:
print(f" Worst among feasible: {max(ms_result.all_objectives):.4f}")
print(f" Best among feasible: {min(ms_result.all_objectives):.8f}")
Multi-start results:
Starts launched: 16
Feasible solutions: 16
Best objective: 0.00000000 (expected: 0.0)
Worst among feasible: 0.0000
Best among feasible: 0.00000000
from discopt._relax.nlp_evaluator import NLPEvaluator
from discopt._relax.primal_heuristics import feasibility_pump
from discopt.solvers.nlp_ipopt import solve_nlp
# Build an MINLP and get the NLP relaxation solution
m2 = dm.Model("pump_demo")
x = m2.continuous("x", lb=0, ub=5)
n = m2.integer("n", lb=0, ub=3)
m2.minimize((x - 1.7) ** 2 + (n - 2) ** 2)
m2.subject_to(x + n <= 4)
# Solve the NLP relaxation (ignoring integrality)
evaluator = NLPEvaluator(m2)
lb, ub = evaluator.variable_bounds
x0 = 0.5 * (np.clip(lb, -10, 10) + np.clip(ub, -10, 10))
nlp_result = solve_nlp(evaluator, x0, options={"print_level": 0})
print(f"NLP relaxation: x={nlp_result.x[0]:.4f}, n={nlp_result.x[1]:.4f} (fractional!)")
# Run feasibility pump to get an integer-feasible solution
x_feas = feasibility_pump(m2, nlp_result.x)
if x_feas is not None:
print(f"Feasibility pump: x={x_feas[0]:.4f}, n={x_feas[1]:.4f} (integer-feasible)")
else:
print("Feasibility pump did not find an integer-feasible solution")
******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
Ipopt is released as open source code under the Eclipse Public License (EPL).
For more information visit https://github.com/coin-or/Ipopt
******************************************************************************
NLP relaxation: x=1.7000, n=2.0000 (fractional!)
Feasibility pump: x=1.7000, n=2.0000 (integer-feasible)
6. CUTEst Interface#
discopt can load and solve problems from the CUTEst [Gould et al., 2015] collection, the standard test set for nonlinear optimization with over 1500 problems.
Requirements: pycutest is an optional dependency.
pip install discopt[cutest]
You also need gfortran and the CUTEst/SIFDecode libraries. See the pycutest documentation for setup.
API overview#
from discopt.interfaces.cutest import CUTEstProblem, list_cutest_problems
# List available problems matching criteria
names = list_cutest_problems(constraints="U", max_n=10) # unconstrained, <= 10 vars
# Load a specific problem
prob = CUTEstProblem("ROSENBR")
print(prob.info) # metadata: n, m, classification
print(prob.x0) # starting point
print(prob.bl, prob.bu) # variable bounds
# Create an evaluator for use with discopt's solvers
evaluator = prob.to_evaluator()
obj = evaluator.evaluate_objective(prob.x0)
grad = evaluator.evaluate_gradient(prob.x0)
NLPEvaluatorFromCUTEst presents the same interface as NLPEvaluator, so a CUTEst problem can be handed to any of discopt’s NLP backends (POUNCE or cyipopt).
# Demonstrate the CUTEst interface (only runs if pycutest + CUTEst are installed)
try:
from discopt.interfaces.cutest import CUTEstProblem
prob = CUTEstProblem("ROSENBR")
print(f"Problem: {prob.name}")
print(f" Variables: {prob.n}")
print(f" Constraints: {prob.m}")
print(f" x0 = {prob.x0}")
print(f" Classification: {prob.info.classification}")
evaluator = prob.to_evaluator()
obj = evaluator.evaluate_objective(prob.x0)
grad = evaluator.evaluate_gradient(prob.x0)
print(f" f(x0) = {obj:.6f}")
print(f" ||grad|| = {np.linalg.norm(grad):.6f}")
except ImportError:
print("pycutest is not installed. Install with: pip install discopt[cutest]")
print("Skipping CUTEst demo.")
except RuntimeError as e:
print(f"CUTEst libraries not found: {e}")
print("Install CUTEst and set the CUTEST environment variable.")
print("See: https://jfowkes.github.io/pycutest/")
pycutest is not installed. Install with: pip install discopt[cutest]
Skipping CUTEst demo.
Summary of Advanced Solver Options#
All advanced features are controlled through m.solve() parameters or dedicated API functions:
Feature |
How to Use |
|---|---|
Differentiable solving (L1) |
|
Differentiable solving (L3) |
|
sIPOPT sensitivity queries |
|
JAX-native differentiation |
|
NLP backend (default) |
|
NLP backend (cyipopt) |
|
Warm-started Rust simplex MILP B&B |
|
B&B nodes exported per iteration |
|
Adaptive McCormick partitioning |
|
Cutting planes |
|
Multi-start |
|
Feasibility pump |
|
CUTEst problems |
|
nlp_solver="ipm" / "sparse_ipm" are accepted as back-compat aliases and resolve
to POUNCE on the NLP/MINLP path; the pure-JAX IPM they used to name is retired.
See the Sensitivity Analysis notebook for detailed examples of the sIPOPT-style features [Pirnay et al., 2012].