ripopt: Rust Interior Point Method

ripopt: Rust Interior Point Method#

This tutorial introduces ripopt, a native Rust implementation of an interior-point method (IPM) for nonlinear programming. ripopt is compiled to a shared library via maturin and exposed to Python through PyO3 bindings.

By the end of this notebook you will be able to:

  1. Solve NLPs using the nlp_solver="ripopt" backend.

  2. Understand the architecture: Rust solver with Python evaluation callbacks.

  3. Configure ripopt solver options.

  4. Compare ripopt performance against the pure-JAX IPM and Ipopt backends.

  5. Use ripopt as the NLP subsolver inside Branch-and-Bound for MINLPs.

When to choose ripopt:

  • Deployment without C dependencies. Unlike Ipopt (which requires BLAS, LAPACK, and optionally HSL), ripopt is a single compiled .so with no external C/Fortran dependencies.

  • Single-threaded performance. Rust’s zero-cost abstractions and lack of garbage collection yield competitive single-threaded speed.

  • No Python GIL concerns. The Rust solver releases the GIL during computation, making it safe to use in multi-threaded Python applications.

References: The IPM algorithm is based on the primal-dual interior-point framework described in [Nocedal and Wright, 2006], Chapter 19.

import os

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

import time

import discopt.modeling as dm
import numpy as np

Architecture#

ripopt follows a callback-based design common to NLP solvers like Ipopt [Wächter and Biegler, 2006]. The architecture has four layers:

Python (discopt.modeling)
  |  formulates Model, builds NLPEvaluator
  v
Python (discopt.solvers.nlp_ripopt)
  |  maps evaluator callbacks to ripopt API
  v
PyO3 bindings (discopt._rust.solve_ripopt)
  |  wraps Python evaluator as ripopt::NlpProblem
  v
Rust solver (ripopt crate)
  |  primal-dual IPM with line search
  |  calls back into Python for f(x), grad(x), Hessian(x)

Key points:

  • The Rust solver owns the IPM iteration loop, KKT system assembly, and line search logic.

  • Objective, gradient, constraint, Jacobian, and Hessian evaluations are performed by the Python NLPEvaluator via GIL-acquiring callbacks.

  • No external C libraries (BLAS, LAPACK, HSL) are required — all linear algebra is implemented in Rust.

  • Variable and constraint bounds are passed once at initialization; the Rust solver caches them as owned Vec<f64> data.

Simple NLP Examples#

Let us start with two small problems to verify that ripopt produces correct results.

# Unconstrained quadratic: min (x - 3)^2 + (y + 1)^2
# Known optimal: x* = 3, y* = -1, f* = 0

m = dm.Model("unconstrained_quad")
x = m.continuous("x", lb=-10, ub=10)
y = m.continuous("y", lb=-10, ub=10)

m.minimize((x - 3) ** 2 + (y + 1) ** 2)

result = m.solve(nlp_solver="ripopt")

print(f"Status:    {result.status}")
print(f"Objective: {result.objective:.6f}")
print(f"x = {result.x['x']:.6f}  (expected 3.0)")
print(f"y = {result.x['y']:.6f}  (expected -1.0)")
print(f"Time:      {result.wall_time:.4f}s")

Status:    optimal
Objective: 0.000000
x = 3.000000  (expected 3.0)
y = -1.000000  (expected -1.0)
Time:      0.4493s

# Constrained: min x^2 + y^2  s.t. x + y >= 1
# Known optimal: x* = y* = 0.5, f* = 0.5

m2 = dm.Model("constrained_simple")
x = m2.continuous("x", lb=-5, ub=5)
y = m2.continuous("y", lb=-5, ub=5)

m2.minimize(x**2 + y**2)
m2.subject_to(x + y >= 1.0, name="sum_lb")

result2 = m2.solve(nlp_solver="ripopt")

print(f"Status:    {result2.status}")
print(f"Objective: {result2.objective:.6f}  (expected 0.5)")
print(f"x = {result2.x['x']:.6f}  (expected 0.5)")
print(f"y = {result2.x['y']:.6f}  (expected 0.5)")

Status:    optimal
Objective: 0.500000  (expected 0.5)
x = 0.500000  (expected 0.5)
y = 0.500000  (expected 0.5)

Constrained Rosenbrock#

The Rosenbrock function is a classic nonlinear test problem with a narrow curved valley. We add a disk constraint to make it more challenging:

\[ \min_{x,y} \quad (1 - x)^2 + 100\,(y - x^2)^2 \qquad \text{s.t.} \quad x^2 + y^2 \le 2, \quad x \in [-2, 2],\; y \in [-2, 2]. \]

The unconstrained optimum is \((1, 1)\) with \(f^* = 0\). The disk constraint \(x^2 + y^2 \le 2\) is inactive at that point, so the constrained optimum should also be \((1, 1)\).

mr = dm.Model("rosenbrock")
x = mr.continuous("x", lb=-2, ub=2)
y = mr.continuous("y", lb=-2, ub=2)

mr.minimize((1 - x) ** 2 + 100 * (y - x**2) ** 2)
mr.subject_to(x**2 + y**2 <= 2.0, name="disk")

rr = mr.solve(nlp_solver="ripopt")

print(f"Status:    {rr.status}")
print(f"Objective: {rr.objective:.8f}  (expected ~0)")
print(f"x = {rr.x['x']:.6f}  (expected 1.0)")
print(f"y = {rr.x['y']:.6f}  (expected 1.0)")
print(f"Disk constraint: x^2 + y^2 = {rr.x['x'] ** 2 + rr.x['y'] ** 2:.4f} <= 2.0")

Status:    optimal
Objective: 0.00000000  (expected ~0)
x = 0.999998  (expected 1.0)
y = 0.999996  (expected 1.0)
Disk constraint: x^2 + y^2 = 2.0000 <= 2.0

Configuration#

ripopt inherits solver options from the solve() call. The most commonly used options are:

Option

Default

Description

time_limit

300

Maximum wall-clock seconds

max_iter

3000

Maximum IPM iterations

tol

1e-8

Optimality tolerance

print_level

0

Verbosity (0 = silent, 1+ = verbose)

These are passed through m.solve() keyword arguments.

# Solve with a time limit and iteration cap
mc = dm.Model("config_demo")
x = mc.continuous("x", lb=-10, ub=10)
y = mc.continuous("y", lb=-10, ub=10)
mc.minimize((x - 2) ** 2 + (y - 3) ** 2)

rc = mc.solve(nlp_solver="ripopt", time_limit=10)
print(f"Status: {rc.status}, Objective: {rc.objective:.6f}")
print(f"x = {rc.x['x']:.4f}, y = {rc.x['y']:.4f}")
print(f"Wall time: {rc.wall_time:.4f}s (limit was 10s)")

Status: optimal, Objective: 0.000000
x = 2.0000, y = 3.0000
Wall time: 0.0139s (limit was 10s)

Performance Comparison#

We compare three NLP solver backends on five small problems:

  • "ipm" — pure-JAX interior-point method (Mehrotra predictor-corrector)

  • "ripopt" — Rust IPM via PyO3

  • "ipopt" — Ipopt via cyipopt (requires BLAS/LAPACK)

Each problem is solved once per backend. Note that the first JAX IPM call includes JIT compilation overhead; subsequent calls of the same size reuse the compiled kernel.

The interior-point framework is described in [Wright, 1997].

def make_unconstrained_quad(n=2):
    """min sum_i (x_i - i)^2"""
    m = dm.Model(f"quad_{n}")
    x = m.continuous("x", shape=(n,), lb=-10, ub=10)
    m.minimize(dm.sum(lambda i: (x[i] - i) ** 2, over=range(n)))
    return m


def make_constrained(n=5):
    """min sum x_i^2  s.t. sum x_i >= n/2"""
    m = dm.Model(f"constr_{n}")
    x = m.continuous("x", shape=(n,), lb=-5, ub=5)
    m.minimize(dm.sum(lambda i: x[i] ** 2, over=range(n)))
    m.subject_to(dm.sum(x) >= n / 2.0)
    return m


def make_rosenbrock():
    """min (1-x)^2 + 100*(y-x^2)^2"""
    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


def make_sum_of_squares(n=10):
    """min sum_i (x_i - sin(i))^2  s.t. sum x_i <= n"""
    m = dm.Model(f"sos_{n}")
    x = m.continuous("x", shape=(n,), lb=-5, ub=5)
    targets = [float(np.sin(i)) for i in range(n)]
    m.minimize(dm.sum(lambda i: (x[i] - targets[i]) ** 2, over=range(n)))
    m.subject_to(dm.sum(x) <= float(n))
    return m


def make_portfolio(n=8):
    """min w' Sigma w  s.t. sum w = 1, mu'w >= 0.08"""
    np.random.seed(42)
    L = np.random.randn(n, n) * 0.05
    Sigma = L @ L.T + 0.01 * np.eye(n)
    mu = np.linspace(0.03, 0.15, n)
    m = dm.Model(f"portfolio_{n}")
    w = m.continuous("w", shape=(n,), lb=0, ub=1)
    m.minimize(
        dm.sum(
            lambda i: dm.sum(
                lambda j: Sigma[i, j] * w[i] * w[j],
                over=range(n),
            ),
            over=range(n),
        )
    )
    m.subject_to(dm.sum(w) == 1.0)
    m.subject_to(dm.sum(lambda i: mu[i] * w[i], over=range(n)) >= 0.08)
    return m


problems = [
    ("Quadratic (n=2)", make_unconstrained_quad(2)),
    ("Constrained (n=5)", make_constrained(5)),
    ("Rosenbrock (n=2)", make_rosenbrock()),
    ("Sum-of-squares (n=10)", make_sum_of_squares(10)),
    ("Portfolio (n=8)", make_portfolio(8)),
]

print(f"Prepared {len(problems)} test problems.")

Prepared 5 test problems.

backends = ["ipm", "ripopt", "ipopt"]
results_table = []

for name, model in problems:
    row = {"Problem": name}
    for backend in backends:
        try:
            t0 = time.perf_counter()
            r = model.solve(nlp_solver=backend)
            elapsed = time.perf_counter() - t0
            row[f"{backend}_time"] = f"{elapsed:.4f}s"
            row[f"{backend}_obj"] = f"{r.objective:.6f}"
            row[f"{backend}_status"] = r.status
        except Exception as e:
            row[f"{backend}_time"] = "error"
            row[f"{backend}_obj"] = str(e)[:30]
            row[f"{backend}_status"] = "error"
    results_table.append(row)

# Display results
header = f"{'Problem':<25s} | {'ipm':>10s} | {'ripopt':>10s} | {'ipopt':>10s}"
print(header)
print("-" * len(header))
for row in results_table:
    print(
        f"{row['Problem']:<25s} | "
        f"{row['ipm_time']:>10s} | "
        f"{row['ripopt_time']:>10s} | "
        f"{row['ipopt_time']:>10s}"
    )

******************************************************************************
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
******************************************************************************

Problem                   |        ipm |     ripopt |      ipopt
----------------------------------------------------------------
Quadratic (n=2)           |    0.0129s |    0.0117s |    0.0120s
Constrained (n=5)         |    0.0248s |    0.0007s |    0.0005s
Rosenbrock (n=2)          |    0.0708s |    0.0013s |    0.0039s
Sum-of-squares (n=10)     |    0.8051s |    0.0568s |    0.0584s
Portfolio (n=8)           |    0.0274s |    0.0011s |    0.0009s

# Verify all backends agree on objective values
print(f"{'Problem':<25s} | {'ipm obj':>12s} | {'ripopt obj':>12s} | {'ipopt obj':>12s}")
print("-" * 70)
for row in results_table:
    print(
        f"{row['Problem']:<25s} | "
        f"{row['ipm_obj']:>12s} | "
        f"{row['ripopt_obj']:>12s} | "
        f"{row['ipopt_obj']:>12s}"
    )

Problem                   |      ipm obj |   ripopt obj |    ipopt obj
----------------------------------------------------------------------
Quadratic (n=2)           |     0.000000 |     0.000000 |     0.000000
Constrained (n=5)         |     1.250000 |     1.250000 |     1.250000
Rosenbrock (n=2)          |     0.000000 |     0.000000 |     0.000000
Sum-of-squares (n=10)     |     0.000000 |     0.000000 |     0.000000
Portfolio (n=8)           |     0.002064 |     0.002064 |     0.002064

Notes on timing:

  • The pure-JAX IPM ("ipm") includes JIT compilation time on the first call. Subsequent solves of the same problem size are significantly faster.

  • ripopt times include the overhead of PyO3 callback crossings (Python to Rust and back for each evaluation).

  • Ipopt typically has the most mature linear algebra (via HSL or MUMPS), which gives it an advantage on larger problems.

Using ripopt in Branch-and-Bound#

When solving a mixed-integer nonlinear program (MINLP), discopt uses Branch-and-Bound (B&B). At each tree node, a continuous NLP relaxation is solved. By setting nlp_solver="ripopt", every node relaxation is solved with the Rust IPM instead of the default backend.

# Small MINLP: min x^2 + y  s.t. x + y >= 2, y in {0, 1}
mm = dm.Model("minlp_ripopt")
x = mm.continuous("x", lb=0, ub=5)
y = mm.binary("y")

mm.minimize(x**2 + y)
mm.subject_to(x + y >= 2.0, name="coupling")

rm = mm.solve(nlp_solver="ripopt")

print(f"Status:     {rm.status}")
print(f"Objective:  {rm.objective:.6f}")
print(f"x = {rm.x['x']:.6f}")
print(f"y = {rm.x['y']:.0f}")
print(f"Node count: {rm.node_count}")
print(f"Gap:        {rm.gap:.2e}")
print(f"Wall time:  {rm.wall_time:.4f}s")

Status:     optimal
Objective:  2.000000
x = 1.000000
y = 1
Node count: 1
Gap:        0.00e+00
Wall time:  0.1589s

# Slightly larger MINLP with 3 binary variables
ml = dm.Model("minlp_larger")
x = ml.continuous("x", shape=(3,), lb=0, ub=5)
z = ml.binary("z", shape=(3,))

# Nonlinear objective
ml.minimize(dm.sum(lambda i: (x[i] - float(i + 1)) ** 2, over=range(3)) + dm.sum(z))

# Linking constraints: x[i] <= 5 * z[i]  (x active only if z selected)
for i in range(3):
    ml.subject_to(x[i] <= 5.0 * z[i], name=f"link_{i}")

# Budget: at most 2 active
ml.subject_to(dm.sum(z) <= 2, name="budget")

rl = ml.solve(nlp_solver="ripopt")

print(f"Status:     {rl.status}")
print(f"Objective:  {rl.objective:.6f}")
print(f"x = {np.round(rl.x['x'], 4)}")
print(f"z = {np.round(rl.x['z'])}")
print(f"Node count: {rl.node_count}")
print(f"Wall time:  {rl.wall_time:.4f}s")

Status:     optimal
Objective:  3.000000
x = [0. 2. 3.]
z = [0. 1. 1.]
Node count: 9
Wall time:  1.4046s

When to Choose ripopt#

Scenario

Recommended backend

Deployment in containers / minimal dependencies

ripopt — no C toolchain needed

Batch solving many small NLPs

ipm — JIT-compiled, amortized overhead

Large-scale NLPs (n > 1000)

ipopt — sparse linear algebra (HSL/MUMPS)

Differentiable optimization / sensitivity

ipm — native JAX autodiff

Single-threaded, no-GIL performance

ripopt — Rust execution

General-purpose robustness

ipopt — decades of development

Drawbacks of ripopt:

  • Dense linear algebra only — no sparse KKT factorization yet, so performance degrades on large problems.

  • Evaluation overhead from PyO3 callbacks (each f(x), grad(x), etc. crosses the Python/Rust boundary).

  • Fewer algorithmic refinements than Ipopt (e.g., no adaptive barrier strategy, no inertia-free regularization).

Exercise#

Formulate a 10-variable constrained NLP and compare ripopt to the default JAX IPM backend. The problem:

\[ \min_{x} \quad \sum_{i=0}^{9} (x_i - t_i)^4 \qquad \text{s.t.} \quad \sum_{i=0}^{9} x_i = 5, \quad x_i \in [-3, 3]. \]

where \(t_i = \sin(i)\).

n = 10
targets = [float(np.sin(i)) for i in range(n)]

# YOUR CODE HERE
# 1. Create a Model
# 2. Add continuous variables x with shape (n,), lb=-3, ub=3
# 3. Minimize sum_i (x[i] - targets[i])^4
# 4. Add equality constraint: sum(x) == 5
# 5. Solve with nlp_solver="ripopt" and time it
# 6. Solve with nlp_solver="ipm" and time it
# 7. Print both objectives and times

Solution#

Click to reveal
def build_quartic_model():
    m = dm.Model("quartic")
    x = m.continuous("x", shape=(n,), lb=-3, ub=3)
    m.minimize(dm.sum(lambda i: (x[i] - targets[i]) ** 4, over=range(n)))
    m.subject_to(dm.sum(x) == 5.0, name="sum_eq")
    return m


# Solve with ripopt
m_rip = build_quartic_model()
t0 = time.perf_counter()
r_rip = m_rip.solve(nlp_solver="ripopt")
t_rip = time.perf_counter() - t0

# Solve with JAX IPM
m_ipm = build_quartic_model()
t0 = time.perf_counter()
r_ipm = m_ipm.solve(nlp_solver="ipm")
t_ipm = time.perf_counter() - t0

print(f"{'Backend':<10s} | {'Status':<10s} | {'Objective':>12s} | {'Time':>8s}")
print("-" * 48)
print(f"{'ripopt':<10s} | {r_rip.status:<10s} | {r_rip.objective:12.6f} | {t_rip:8.4f}s")
print(f"{'ipm':<10s} | {r_ipm.status:<10s} | {r_ipm.objective:12.6f} | {t_ipm:8.4f}s")
print(f"\nObjective difference: {abs(r_rip.objective - r_ipm.objective):.2e}")

Backend    | Status     |    Objective |     Time
------------------------------------------------
ripopt     | optimal    |     0.085947 |   0.1447s
ipm        | optimal    |     0.085947 |   0.1476s

Objective difference: 9.71e-17

Summary#

In this tutorial you learned to:

  • Select the ripopt backend with m.solve(nlp_solver="ripopt").

  • Understand the Rust IPM architecture: compiled .so, PyO3 bindings, Python evaluation callbacks.

  • Configure solver options (time limit, tolerance, verbosity).

  • Compare ripopt against the pure-JAX IPM and Ipopt backends on small-to-medium NLPs.

  • Use ripopt as the NLP subsolver inside Branch-and-Bound for MINLPs.

ripopt is the right choice when you need a dependency-free NLP solver for deployment, or when you want single-threaded Rust performance without the GIL. For large-scale problems or when sparse linear algebra is critical, Ipopt remains the stronger option [Wächter and Biegler, 2006]. For batch solving or differentiable optimization, the pure-JAX IPM is preferred [Nocedal and Wright, 2006].

See also: