NLP Backend Comparison: ripopt vs IPM vs Ipopt#
discopt supports three NLP backends for solving continuous relaxations:
Backend |
Flag |
Description |
|---|---|---|
ripopt |
|
Rust interior-point method via PyO3 bindings |
discopt IPM |
|
Pure-JAX primal-dual interior point method [Nocedal and Wright, 2006] (default) |
Ipopt |
|
Industry-standard NLP solver [Wächter and Biegler, 2006] via cyipopt |
This notebook compares them on a representative set of problems, measuring:
Accuracy: objective value agreement with known optima
Speed: wall-clock solve time
Robustness: convergence status
import os
os.environ["JAX_PLATFORMS"] = "cpu"
os.environ["JAX_ENABLE_X64"] = "1"
import time
import discopt.modeling as dm
import numpy as np
print("discopt loaded successfully")
discopt loaded successfully
Test Problems#
We define a mix of problem types built with dm.Model() and solve with all three backends.
All three backends support both dm.Model()-built models and .nl file models. The .nl evaluator uses the Rust expression evaluator with finite differences, which is compatible with all backends.
def make_unconstrained_quadratic():
"""Simple unconstrained quadratic. Optimal: f*=0 at (3, -1)."""
m = dm.Model("quadratic")
x = m.continuous("x", lb=-10, ub=10)
y = m.continuous("y", lb=-10, ub=10)
m.minimize((x - 3) ** 2 + (y + 1) ** 2)
return m, 0.0
def make_rosenbrock():
"""Rosenbrock function. Optimal: f*=0 at (1, 1)."""
m = dm.Model("rosenbrock")
x = m.continuous("x", lb=-5, ub=5)
y = m.continuous("y", lb=-5, ub=5)
m.minimize((1 - x) ** 2 + 100 * (y - x**2) ** 2)
return m, 0.0
def make_constrained_quadratic():
"""Constrained quadratic. Optimal: f*=0.5 at (0.5, 0.5)."""
m = dm.Model("constrained_quad")
x = m.continuous("x", lb=-5, ub=5)
y = m.continuous("y", lb=-5, ub=5)
m.minimize(x**2 + y**2)
m.subject_to(x + y >= 1)
return m, 0.5
def make_nonlinear_equality():
"""Nonlinear equality constraint. Optimal: f*=2 at (1, 1)."""
m = dm.Model("nl_equality")
x = m.continuous("x", lb=0.1, ub=5)
y = m.continuous("y", lb=0.1, ub=5)
m.minimize(x**2 + y**2)
m.subject_to(x * y == 1)
return m, 2.0
def make_exp_nlp():
"""Exponential objective with constraint. Optimal: f*~1.8395."""
from scipy.optimize import fsolve
# Solve KKT: exp(x)=2y, x+y=1 => exp(x)=2(1-x)
x_star = fsolve(lambda x: np.exp(x) - 2 * (1 - x), 0.0)[0]
y_star = 1 - x_star
opt = np.exp(x_star) + y_star**2
m = dm.Model("exp_nlp")
x = m.continuous("x", lb=-2, ub=2)
y = m.continuous("y", lb=-2, ub=2)
m.minimize(dm.exp(x) + y**2)
m.subject_to(x + y >= 1)
return m, float(opt)
def make_log_nlp():
"""Logarithmic objective. Optimal: f*=2*ln(2)~1.3863 at (0.5, 0.5)."""
import math
m = dm.Model("log_nlp")
x = m.continuous("x", lb=0.01, ub=0.99)
y = m.continuous("y", lb=0.01, ub=0.99)
m.minimize(-dm.log(x) - dm.log(y))
m.subject_to(x + y <= 1)
return m, 2 * math.log(2)
def make_multi_constraint():
"""Multiple constraints. Optimal: f*=0.5 at (1.5, 1.5)."""
m = dm.Model("multi_constraint")
x = m.continuous("x", lb=0, ub=5)
y = m.continuous("y", lb=0, ub=5)
m.minimize((x - 2) ** 2 + (y - 2) ** 2)
m.subject_to(x + y <= 3)
return m, 0.5
def make_bilinear():
"""Bilinear constraints. Optimal: f*=2 at (1, 1)."""
m = dm.Model("bilinear")
x = m.continuous("x", lb=0.1, ub=5)
y = m.continuous("y", lb=0.1, ub=5)
m.minimize(x**2 + y**2)
m.subject_to(x * y >= 1)
m.subject_to(x + y <= 4)
return m, 2.0
def make_trig_nlp():
"""Trigonometric objective. Optimal: f*~-1.0 at x~pi/2."""
m = dm.Model("trig_nlp")
x = m.continuous("x", lb=0.1, ub=3.0)
y = m.continuous("y", lb=0.1, ub=3.0)
m.minimize(-dm.sin(x) * dm.cos(y))
m.subject_to(x + y <= 4)
return m, -1.0 # approximate
problems = [
("Unconstrained quadratic", make_unconstrained_quadratic),
("Rosenbrock", make_rosenbrock),
("Constrained quadratic", make_constrained_quadratic),
("Nonlinear equality", make_nonlinear_equality),
("Exponential NLP", make_exp_nlp),
("Logarithmic NLP", make_log_nlp),
("Multiple constraints", make_multi_constraint),
("Bilinear constraints", make_bilinear),
("Trigonometric NLP", make_trig_nlp),
]
print(f"Defined {len(problems)} test problems")
Defined 9 test problems
Run Comparison#
For each problem, solve with all three backends and record objective, time, and status.
backends = ["ripopt", "ipm", "ipopt"]
results = []
for name, make_fn in problems:
model, expected = make_fn()
row = {"problem": name, "expected": expected}
for backend in backends:
m_copy, _ = make_fn() # fresh model each time
try:
t0 = time.perf_counter()
result = m_copy.solve(nlp_solver=backend, time_limit=30)
elapsed = time.perf_counter() - t0
row[f"{backend}_obj"] = result.objective
row[f"{backend}_time"] = elapsed
row[f"{backend}_status"] = result.status
if result.objective is not None:
row[f"{backend}_error"] = abs(result.objective - expected)
else:
row[f"{backend}_error"] = None
except Exception as e:
row[f"{backend}_obj"] = None
row[f"{backend}_time"] = None
row[f"{backend}_status"] = f"ERROR: {e}"
row[f"{backend}_error"] = None
results.append(row)
print(f"Solved {len(results)} problems with {len(backends)} backends")
******************************************************************************
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
******************************************************************************
Solved 9 problems with 3 backends
Results: Accuracy Comparison#
header = (
f"{'Problem':<26s} {'Expected':>10s} "
f"{'ripopt obj':>12s} {'ripopt err':>10s} "
f"{'IPM obj':>12s} {'IPM err':>10s} "
f"{'Ipopt obj':>12s} {'Ipopt err':>10s}"
)
print(header)
print("-" * len(header))
for r in results:
cols = []
cols.append(f"{r['problem']:<26s} {r['expected']:>10.6f}")
for b in backends:
obj = r.get(f"{b}_obj")
err = r.get(f"{b}_error")
obj_s = f"{obj:.6f}" if obj is not None else "FAIL"
err_s = f"{err:.2e}" if err is not None else "N/A"
cols.append(f"{obj_s:>12s} {err_s:>10s}")
print(" ".join(cols))
Problem Expected ripopt obj ripopt err IPM obj IPM err Ipopt obj Ipopt err
-------------------------------------------------------------------------------------------------------------
Unconstrained quadratic 0.000000 0.000000 2.31e-14 0.000000 2.31e-14 0.000000 2.31e-14
Rosenbrock 0.000000 0.000000 1.05e-25 0.000000 9.83e-20 0.000000 9.83e-20
Constrained quadratic 0.500000 0.500000 1.11e-16 0.500000 1.11e-16 0.500000 1.11e-16
Nonlinear equality 2.000000 2.000000 9.50e-14 2.000000 0.00e+00 2.000000 0.00e+00
Exponential NLP 1.839484 1.839484 1.38e-14 1.839484 1.12e-08 1.839484 1.12e-08
Logarithmic NLP 1.386294 1.386294 3.69e-14 1.386294 1.75e-08 1.386294 1.75e-08
Multiple constraints 0.500000 0.500000 0.00e+00 0.500000 0.00e+00 0.500000 0.00e+00
Bilinear constraints 2.000000 2.000000 1.00e-11 2.000000 1.75e-08 2.000000 1.75e-08
Trigonometric NLP -1.000000 -0.995004 5.00e-03 -0.995004 5.00e-03 -0.995004 5.00e-03
Results: Speed Comparison#
header = f"{'Problem':<26s} {'ripopt':>10s} {'IPM':>10s} {'Ipopt':>10s} {'Fastest':>8s}"
print(header)
print("-" * len(header))
totals = {b: 0.0 for b in backends}
for r in results:
times = {}
for b in backends:
t = r.get(f"{b}_time")
times[b] = t
if t is not None:
totals[b] += t
time_strs = []
for b in backends:
t = times[b]
time_strs.append(f"{t:.4f}s" if t is not None else "FAIL")
valid = {b: t for b, t in times.items() if t is not None}
fastest = min(valid, key=valid.get) if valid else "N/A"
print(f"{r['problem']:<26s} " + " ".join(f"{s:>10s}" for s in time_strs) + f" {fastest:>8s}")
print()
for b in backends:
print(f" {b} total: {totals[b]:.3f}s")
Problem ripopt IPM Ipopt Fastest
--------------------------------------------------------------------
Unconstrained quadratic 0.6490s 0.0127s 0.0114s ipopt
Rosenbrock 0.0362s 0.0724s 0.0604s ripopt
Constrained quadratic 0.0218s 0.0006s 0.0005s ipopt
Nonlinear equality 0.0791s 0.0763s 0.0806s ipm
Exponential NLP 0.0815s 0.0834s 0.0852s ripopt
Logarithmic NLP 0.0857s 0.0843s 0.0762s ipopt
Multiple constraints 0.0940s 0.0149s 0.0137s ipopt
Bilinear constraints 0.1315s 0.0862s 0.0869s ipm
Trigonometric NLP 0.0765s 0.0798s 0.0769s ripopt
ripopt total: 1.255s
ipm total: 0.511s
ipopt total: 0.492s
Results: Status Summary#
header = f"{'Problem':<26s} {'ripopt':<18s} {'IPM':<18s} {'Ipopt':<18s}"
print(header)
print("-" * len(header))
for r in results:
cols = [f"{r['problem']:<26s}"]
for b in backends:
s = r.get(f"{b}_status", "N/A")
cols.append(f"{s:<18s}")
print(" ".join(cols))
Problem ripopt IPM Ipopt
-----------------------------------------------------------------------------------
Unconstrained quadratic optimal optimal optimal
Rosenbrock optimal optimal optimal
Constrained quadratic optimal optimal optimal
Nonlinear equality optimal optimal optimal
Exponential NLP optimal optimal optimal
Logarithmic NLP optimal optimal optimal
Multiple constraints optimal optimal optimal
Bilinear constraints optimal optimal optimal
Trigonometric NLP optimal optimal optimal
MINLP Comparison (Branch & Bound)#
For MINLP problems, the NLP backend is used at each Branch & Bound [Land and Doig, 1960] node. Let’s compare on a few small MINLPs.
def make_simple_minlp():
"""Simple MINLP. Optimal: f*=0.5."""
m = dm.Model("simple_minlp")
x1 = m.continuous("x1", lb=0, ub=5)
x2 = m.continuous("x2", lb=0, ub=5)
x3 = m.binary("x3")
m.minimize(x1**2 + x2**2 + x3)
m.subject_to(x1 + x2 >= 1)
m.subject_to(x1**2 + x2 <= 3)
return m, 0.5
def make_binary_knapsack():
"""Binary knapsack-like MINLP. Optimal: f*=5.0."""
m = dm.Model("binary_knapsack")
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(x <= 5 * y1)
m.subject_to(x <= 4 * y2)
return m, 5.0
def make_integer_quadratic():
"""Integer quadratic. Optimal: f*=0.0."""
m = dm.Model("integer_quad")
x = m.continuous("x", lb=0, ub=5)
n = m.integer("n", lb=0, ub=5)
m.minimize((x - 2.7) ** 2 + (n - 3) ** 2)
m.subject_to(x + n <= 6)
return m, 0.0
def make_exp_minlp():
"""Exponential MINLP. Binary switches continuous region."""
m = dm.Model("exp_minlp")
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)
return m, None # optimal unknown a priori
minlp_problems = [
("Simple MINLP", make_simple_minlp),
("Binary knapsack", make_binary_knapsack),
("Integer quadratic", make_integer_quadratic),
("Exp MINLP", make_exp_minlp),
]
minlp_results = []
for name, make_fn in minlp_problems:
row = {"problem": name}
_, expected = make_fn()
row["expected"] = expected
for backend in backends:
m_copy, _ = make_fn()
try:
t0 = time.perf_counter()
result = m_copy.solve(nlp_solver=backend, time_limit=30, max_nodes=500)
elapsed = time.perf_counter() - t0
row[f"{backend}_obj"] = result.objective
row[f"{backend}_time"] = elapsed
row[f"{backend}_status"] = result.status
row[f"{backend}_nodes"] = result.node_count
except Exception as e:
row[f"{backend}_obj"] = None
row[f"{backend}_time"] = None
row[f"{backend}_status"] = f"ERROR: {e}"
row[f"{backend}_nodes"] = None
minlp_results.append(row)
# Print MINLP results
header = f"{'Problem':<20s} " + " ".join(
f"{b + ' obj':>12s} {b + ' time':>10s} {b + ' nodes':>8s}" for b in backends
)
print(header)
print("-" * len(header))
for r in minlp_results:
cols = [f"{r['problem']:<20s}"]
for b in backends:
obj = r.get(f"{b}_obj")
t = r.get(f"{b}_time")
n = r.get(f"{b}_nodes")
obj_s = f"{obj:.4f}" if obj is not None else "FAIL"
t_s = f"{t:.3f}s" if t is not None else "N/A"
n_s = str(n) if n is not None else "N/A"
cols.append(f"{obj_s:>12s} {t_s:>10s} {n_s:>8s}")
print(" ".join(cols))
Problem ripopt obj ripopt time ripopt nodes ipm obj ipm time ipm nodes ipopt obj ipopt time ipopt nodes
--------------------------------------------------------------------------------------------------------------------------------
Simple MINLP 0.5000 0.206s 1 0.5000 1.205s 1 0.5000 0.205s 1
Binary knapsack 5.0000 0.940s 5 5.0000 0.045s 5 5.0000 0.043s 5
Integer quadratic -0.0000 0.171s 1 -0.0000 0.033s 1 -0.0000 0.031s 1
Exp MINLP 0.6686 2.224s 1 0.6686 2.164s 1 0.6686 1.316s 1
Summary#
Key takeaways:
ripopt (
nlp_solver="ripopt") is the project’s own Rust interior-point solver, integrated via PyO3 bindings. It provides the Hessian of the full Lagrangian for better convergence with constraints.discopt IPM (
nlp_solver="ipm", default) is a pure-JAX implementation that supportsjax.vmapfor batch solving at B&B nodes and can run on GPU.Ipopt (
nlp_solver="ipopt") is the industry-standard NLP solver, more mature and robust. Required for differentiable solving.All three backends should produce the same optimal objective values (to solver tolerance).
The IPM backend’s advantage is batch solving on GPU — this notebook uses CPU only.
ripopt avoids the cyipopt/Ipopt dependency and provides a fully self-contained solver stack.