cyipopt: The Ipopt Solver Backend

cyipopt: The Ipopt Solver Backend#

This tutorial covers the cyipopt solver backend in discopt — a Python binding for the industry-standard Ipopt (Interior Point OPTimizer) solver. By the end, you will be able to:

  1. Select the Ipopt backend via nlp_solver="ipopt" and understand when it is the right choice.

  2. Configure Ipopt options (iteration limits, tolerances, verbosity).

  3. Interpret Ipopt’s iteration output.

  4. Solve .nl files from MINLPLib using the Ipopt backend.

  5. Understand how discopt wraps its evaluator callbacks for cyipopt.

Prerequisites: familiarity with the discopt modeling API (see the quickstart tutorial) and basic knowledge of nonlinear programming.

Ipopt implements a primal-dual interior-point algorithm with a filter line-search strategy [Wächter and Biegler, 2006].

import time

import discopt.modeling as dm
import numpy as np

Ipopt Background#

Ipopt (Interior Point OPTimizer) is an open-source solver for large-scale nonlinear programming, first released in 2002 as part of the COIN-OR project. It implements a primal-dual interior-point method with a filter line-search mechanism that accepts trial steps when they either reduce the objective or reduce constraint violation — avoiding the need to tune a penalty parameter [Wächter and Biegler, 2006].

Key characteristics:

  • Sparse linear algebra. Ipopt factors the KKT system using sparse symmetric indefinite solvers (MA27, MA57, MUMPS). This makes it highly efficient for problems with hundreds or thousands of variables.

  • Robust convergence. Decades of testing across chemical engineering, optimal control, and operations research have made Ipopt one of the most reliable NLP solvers available.

  • Second-order information. Ipopt uses exact Hessians of the Lagrangian when provided (which discopt supplies from the POUNCE AD tape in the Rust core — JAX is not on the default solve path), or falls back to L-BFGS approximations.

The foundational reference is [Nocedal and Wright, 2006], Chapter 19, which covers interior-point methods for nonlinear programming in detail.

Installation#

The cyipopt package provides Python bindings for Ipopt. It requires the Ipopt C library to be installed separately.

macOS (Homebrew):

brew install ipopt
pip install cyipopt

Ubuntu / Debian:

sudo apt-get install coinor-libipopt-dev
pip install cyipopt

Conda (all platforms):

conda install -c conda-forge cyipopt

You can verify the installation by importing cyipopt:

import cyipopt
print(cyipopt.__version__)

Basic NLP Solve#

We start with a medium-sized NLP: 10 variables, a quadratic objective, and coupling constraints that link consecutive variables. The key API is m.solve(nlp_solver="ipopt"), which routes the NLP through cyipopt instead of the default POUNCE backend.

n = 10
m = dm.Model("coupled_quadratic")
x = m.continuous("x", shape=(n,), lb=-5.0, ub=10.0)

# Objective: sum of (x[i] - i)^2
m.minimize(dm.sum(lambda i: (x[i] - float(i)) ** 2, over=range(n)))

# Coupling constraints: x[i] + x[i+1] <= 2*i + 3
for i in range(n - 1):
    m.subject_to(x[i] + x[i + 1] <= 2.0 * i + 3.0, name=f"couple_{i}")

# Sum constraint: total <= 30
m.subject_to(dm.sum(lambda i: x[i], over=range(n)) <= 30.0, name="budget")

print(m.summary())

Model: coupled_quadratic
  Variables: 10 (10 continuous, 0 integer/binary)
  Constraints: 10
  Objective: minimize Σ[10 terms]
  Parameters: 0

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

print(f"Status   : {result.status}")
print(f"Objective: {result.objective:.6f}")
print(f"Wall time: {result.wall_time:.4f}s")
print("\nSolution:")
for i in range(n):
    print(f"  x[{i}] = {result.x['x'][i]:>8.4f}  (target = {i})")

Status   : optimal
Objective: 22.500000
Wall time: 0.0342s

Solution:
  x[0] =  -1.5000  (target = 0)
  x[1] =  -0.5000  (target = 1)
  x[2] =   0.5000  (target = 2)
  x[3] =   1.5000  (target = 3)
  x[4] =   2.5000  (target = 4)
  x[5] =   3.5000  (target = 5)
  x[6] =   4.5000  (target = 6)
  x[7] =   5.5000  (target = 7)
  x[8] =   6.5000  (target = 8)
  x[9] =   7.5000  (target = 9)

Ipopt converges quickly on this smooth, convex problem. The coupling constraints prevent some variables from reaching their targets exactly, which is why the objective is greater than zero.

Understanding Ipopt’s Iteration Output#

Setting print_level to 5 (the Ipopt default) produces a detailed iteration log. This is invaluable for diagnosing convergence issues.

result_verbose = m.solve(
    nlp_solver="ipopt",
    ipopt_options={"print_level": 5, "max_iter": 50},
)
print(f"\nFinal status: {result_verbose.status}")
print(f"Objective   : {result_verbose.objective:.6f}")

Final status: optimal
Objective   : 22.500000

Reading the Iteration Log#

Each row in the Ipopt iteration log contains:

Column

Meaning

iter

Iteration number

objective

Current objective value

inf_pr

Primal infeasibility (max constraint violation)

inf_du

Dual infeasibility (max KKT violation)

lg(mu)

Log10 of the barrier parameter

`

lg(rg)

Log10 of the inertia regularization (if any)

alpha_du

Dual step size

alpha_pr

Primal step size

ls

Number of line-search backtracks

What to look for:

  • inf_pr and inf_du should both converge toward zero.

  • lg(mu) should decrease steadily (the barrier vanishes at optimality).

  • If ls is consistently large, the line-search is struggling — consider scaling the problem or tightening bounds.

  • The primal-dual interior-point theory behind these quantities is covered in [Wright, 1997].

Larger NLP: Ipopt vs POUNCE#

Let us compare both backends on an \(n = 50\) variable problem with nonlinear constraints.

Warning

nlp_solver="ipm" is a back-compat alias, not a third engine. The pure-JAX interior-point method it once named has been retired; on the NLP/MINLP path it resolves silently to POUNCE, discopt’s pure-Rust port of Ipopt. So this comparison is Ipopt vs POUNCE — two interior-point implementations of the same method — and not “JAX vs Ipopt”. The cell asserts that "ipm" and "pounce" return the identical objective, which is what an alias means.

np.random.seed(42)
n_large = 50

# Random positive-definite quadratic
Q_half = np.random.randn(n_large, n_large) * 0.1
q_diag = np.diag(Q_half.T @ Q_half) + 0.1  # diagonal of Q^T Q + 0.1 I

m_large = dm.Model("large_nlp")
x = m_large.continuous("x", shape=(n_large,), lb=-2.0, ub=2.0)

# Quadratic objective: sum of q_i * x_i^2 - c_i * x_i
c_vec = np.random.randn(n_large) * 0.5
m_large.minimize(dm.sum(lambda i: q_diag[i] * x[i] ** 2 - c_vec[i] * x[i], over=range(n_large)))

# Coupling constraints: x[i]^2 + x[i+1]^2 <= 3.0
for i in range(0, n_large - 1, 3):
    m_large.subject_to(
        x[i] ** 2 + x[min(i + 1, n_large - 1)] ** 2 <= 3.0,
        name=f"quad_couple_{i}",
    )

print(m_large.summary())

Model: large_nlp
  Variables: 50 (50 continuous, 0 integer/binary)
  Constraints: 17
  Objective: minimize Σ[50 terms]
  Parameters: 0

# Solve with Ipopt
t0 = time.perf_counter()
res_ipopt = m_large.solve(nlp_solver="ipopt")
t_ipopt = time.perf_counter() - t0

# Solve with POUNCE. `nlp_solver="ipm"` is an alias that resolves to it; both
# spellings are run so the alias is demonstrated rather than asserted in prose.
t0 = time.perf_counter()
res_pounce = m_large.solve(nlp_solver="pounce")
t_pounce = time.perf_counter() - t0
res_ipm = m_large.solve(nlp_solver="ipm")
assert res_ipm.objective == res_pounce.objective, (
    'nlp_solver="ipm" is meant to resolve to POUNCE: '
    f"{res_ipm.objective} != {res_pounce.objective}"
)

print(f"{'Backend':<12} {'Status':<12} {'Objective':>12} {'Time (s)':>10}")
print("-" * 48)
print(f"{'Ipopt':<12} {str(res_ipopt.status):<12} {res_ipopt.objective:>12.6f} {t_ipopt:>10.4f}")
print(f"{'POUNCE':<12} {str(res_pounce.status):<12} {res_pounce.objective:>12.6f} {t_pounce:>10.4f}")
print(f"\nObjective difference: {abs(res_ipopt.objective - res_pounce.objective):.2e}")
print('nlp_solver="ipm" resolved to the same engine: '
      f"{res_ipm.objective == res_pounce.objective}")

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

Backend      Status          Objective   Time (s)
------------------------------------------------
Ipopt        optimal         -6.629880     0.0280
POUNCE       optimal         -6.629880     0.0114

Objective difference: 0.00e+00
nlp_solver="ipm" resolved to the same engine: True

Both backends find the same optimal objective. They are two implementations of the same interior-point method [Wächter and Biegler, 2006]: cyipopt calls the reference C++ Ipopt, POUNCE is discopt’s pure-Rust port of it. The timings above are a single un-replicated pair on one problem and should not be read as a benchmark — see Benchmarks by Problem Class for measured comparisons with a spread.

Solving MINLPLib .nl Files#

The .nl file format (from AMPL) is a standard interchange format for nonlinear programs. discopt can load .nl files via dm.from_nl() and solve them with any backend, including Ipopt [Bussieck et al., 2003].

from pathlib import Path

nl_path = Path("../../python/tests/data/minlplib/st_e13.nl")
model_nl = dm.from_nl(str(nl_path))

result_nl = model_nl.solve(nlp_solver="ipopt", time_limit=30)

print("Problem : st_e13")
print(f"Status  : {result_nl.status}")
print(f"Objective: {result_nl.objective:.6f}")
print(f"Wall time: {result_nl.wall_time:.4f}s")

Problem : st_e13
Status  : optimal
Objective: 2.000000
Wall time: 0.0117s

The from_nl() function parses the .nl file into a discopt expression DAG using a Rust-based parser (discopt-core). From there, the same evaluator pipeline builds the objective, gradient, Hessian, constraint, and Jacobian callbacks that cyipopt needs — there is no separate .nl code path.

Callback Structure: How discopt Wraps cyipopt#

When you call m.solve(nlp_solver="ipopt"), the following pipeline executes:

  1. DAG compilation. The modeling API’s expression DAG is compiled into a Rust AD tape by TapeNLPEvaluator (discopt._tape_nlp_evaluator). A model the tape cannot represent falls back to the JAX-backed NLPEvaluator, which is where JAX still earns its place — but that is the exception. Measured on the model in this notebook: an nlp_solver="ipopt" solve loads zero jax modules and the evaluator is a TapeNLPEvaluator.

  2. Callback adapter. An internal _IpoptCallbacks class maps the NLPEvaluator methods to the interface that cyipopt.Problem expects:

    • objective(x) — scalar objective value

    • gradient(x) — gradient of the objective

    • constraints(x) — vector of constraint body values

    • jacobian(x) — flattened constraint Jacobian

    • hessian(x, lagrange, obj_factor) — lower triangle of the Hessian of the Lagrangian

  3. Automatic differentiation. The tape computes exact gradients, Jacobians, and Hessians by reverse-mode AD. This is more accurate than finite differences and avoids the \(O(n)\) cost of forward-mode for gradients. What reaches Ipopt is the same mathematics either way — the change since this section was written is which engine produces it.

  4. Sparsity. Currently, discopt provides dense Jacobian and Hessian arrays to cyipopt (with dense sparsity structures). For very large problems, the sparse Jacobian and sparse KKT modules can be used to exploit structure.

Note. Ipopt expects the Hessian in CSR sparse lower-triangular form. The callback adapter extracts the lower triangle from the dense Hessian the evaluator returns and flattens it to match hessianstructure().

from discopt._relax.nlp_evaluator import NLPEvaluator

# Build the evaluator for our first model
m_demo = dm.Model("demo")
x_demo = m_demo.continuous("x", shape=(3,), lb=-5, ub=5)
m_demo.minimize(x_demo[0] ** 2 + x_demo[1] ** 2 + x_demo[2] ** 2)
m_demo.subject_to(x_demo[0] + x_demo[1] + x_demo[2] >= 1.0, name="sum_lb")

ev = NLPEvaluator(m_demo)
x_test = np.array([1.0, 0.5, 0.2])

print(f"Objective   : {ev.evaluate_objective(x_test):.4f}")
print(f"Gradient    : {ev.evaluate_gradient(x_test)}")
print(f"Constraints : {ev.evaluate_constraints(x_test)}")
print(f"Jacobian    :\n{ev.evaluate_jacobian(x_test)}")
print(f"Hessian     :\n{ev.evaluate_hessian(x_test)}")

Objective   : 1.2900
Gradient    : [2.  1.  0.4]
Constraints : [-0.7]
Jacobian    :
[[-1. -1. -1.]]
Hessian     :
[[2. 0. 0.]
 [0. 2. 0.]
 [0. 0. 2.]]

When cyipopt Is the Best Choice#

Use nlp_solver="ipopt" when:

  • Problem size exceeds ~200 variables. Ipopt’s sparse linear algebra (MA27/MUMPS) scales much better than dense factorization.

  • You need battle-tested robustness. Ipopt has been validated on thousands of benchmark instances over two decades.

  • You need specific Ipopt options — warm-starting, custom linear solvers, adaptive barrier strategies, etc.

  • Working with .nl files from AMPL, Pyomo, or other modeling languages.

  • Constrained problems with many inequality constraints. Ipopt’s filter line-search handles constraint activation robustly.

Use nlp_solver="pounce" (the default) when:

  • The problem is small to medium (\(n < 200\)) and dense.

  • You want no external toolchain — POUNCE is pure Rust and a core dependency, so there is no libipopt, no BLAS/LAPACK, and nothing to build.

  • You are deploying somewhere a C/Fortran stack is inconvenient.

There is no third option to weigh: nlp_solver="ipm" and "sparse_ipm" are back-compat aliases for POUNCE (see the warning above), so choosing one of them is choosing this row. Differentiable optimization is a separate axis, not a backend choice — see Sensitivity Analysis and Advanced Features.

Ipopt Options Reference#

Pass options via ipopt_options={...} in m.solve(). Below are the most commonly used options:

Option

Type

Default

Description

print_level

int

0 (in discopt)

Verbosity: 0 = silent, 5 = default Ipopt, 12 = debug

max_iter

int

3000

Maximum number of iterations

tol

float

1e-8

Convergence tolerance (KKT error)

acceptable_tol

float

1e-6

Acceptable convergence tolerance

max_wall_time

float

1e6

Maximum wall-clock time (seconds)

mu_strategy

str

"monotone"

Barrier update: "monotone" or "adaptive"

linear_solver

str

"mumps"

Linear solver: "mumps", "ma27", "ma57", "pardiso"

warm_start_init_point

str

"no"

Use "yes" to warm-start from previous solution

hessian_approximation

str

"exact"

Use "limited-memory" for L-BFGS Hessian

The full list of Ipopt options is documented in the Ipopt documentation.

Tip. discopt sets print_level=0 by default to keep output clean. Set it to 5 for debugging, or use the ipopt_options parameter to override any default.

Exercise: Constrained Rosenbrock#

The Rosenbrock function is a classic test for nonlinear optimizers:

\[ \min_{x_1, x_2} \quad (1 - x_1)^2 + 100\,(x_2 - x_1^2)^2 \]

subject to: $\( x_1^2 + x_2^2 \le 2, \qquad x_1 + x_2 \ge 0.5 \)$

Task: Formulate this in discopt and solve with both "ipopt" and "ipm". Compare the objective values and wall times.

# YOUR CODE HERE
# 1. Create a Model("rosenbrock")
# 2. Add 2 continuous variables x1, x2 with lb=-1.5, ub=1.5
# 3. Set the Rosenbrock objective
# 4. Add the two constraints
# 5. Solve with nlp_solver="ipopt" and nlp_solver="ipm"
# 6. Print both results

Solution#

Click to reveal
mr = dm.Model("rosenbrock")
x1 = mr.continuous("x1", lb=-1.5, ub=1.5)
x2 = mr.continuous("x2", lb=-1.5, ub=1.5)

mr.minimize((1.0 - x1) ** 2 + 100.0 * (x2 - x1**2) ** 2)
mr.subject_to(x1**2 + x2**2 <= 2.0, name="circle")
mr.subject_to(x1 + x2 >= 0.5, name="halfplane")

results = {}
for backend in ["ipopt", "ipm"]:
    t0 = time.perf_counter()
    res = mr.solve(nlp_solver=backend)
    elapsed = time.perf_counter() - t0
    results[backend] = (res, elapsed)

print(f"{'Backend':<10} {'Status':<12} {'Objective':>12} {'Time (s)':>10}")
print("-" * 46)
for backend, (res, elapsed) in results.items():
    print(f"{backend:<10} {str(res.status):<12} {res.objective:>12.8f} {elapsed:>10.4f}")

# Both should find the constrained optimum near x1=1, x2=1
res_ipopt = results["ipopt"][0]
print(f"\nIpopt solution: x1 = {res_ipopt.x['x1']:.6f}, x2 = {res_ipopt.x['x2']:.6f}")

Backend    Status          Objective   Time (s)
----------------------------------------------
ipopt      optimal        0.00000002     0.6465
ipm        optimal        0.00000002     0.6215

Ipopt solution: x1 = 0.999862, x2 = 0.999724

Summary#

In this tutorial you learned to:

  • Select the Ipopt backend via m.solve(nlp_solver="ipopt").

  • Configure Ipopt using ipopt_options for verbosity, iteration limits, tolerances, and solver strategy.

  • Interpret the iteration log — tracking primal/dual infeasibility, barrier parameter, and step sizes.

  • Solve .nl files from MINLPLib using the Ipopt backend.

  • Understand the callback pipeline — how discopt compiles JAX functions into the interface that cyipopt expects.

Ipopt is the recommended backend for large-scale NLPs (\(n > 200\)) and when robust convergence on difficult constrained problems is required [Wächter and Biegler, 2006].

Next: Solver Selection Guide compares all three backends (JAX IPM, POUNCE, cyipopt) and provides decision criteria for choosing the right one.