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 JAX-compiled 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 os

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

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 via JAX automatic differentiation), 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 JAX IPM 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: 1.0016s

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 JAX IPM#

On larger problems, Ipopt’s sparse linear algebra gives it an advantage over discopt’s dense JAX IPM. Let us compare both backends on an \(n = 50\) variable problem with nonlinear constraints.

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
/var/folders/gq/k1kgbl7n539_4dl1md8x3jt80000gn/T/ipykernel_45480/2292490133.py:6: RuntimeWarning: divide by zero encountered in matmul
  q_diag = np.diag(Q_half.T @ Q_half) + 0.1  # diagonal of Q^T Q + 0.1 I
/var/folders/gq/k1kgbl7n539_4dl1md8x3jt80000gn/T/ipykernel_45480/2292490133.py:6: RuntimeWarning: overflow encountered in matmul
  q_diag = np.diag(Q_half.T @ Q_half) + 0.1  # diagonal of Q^T Q + 0.1 I
/var/folders/gq/k1kgbl7n539_4dl1md8x3jt80000gn/T/ipykernel_45480/2292490133.py:6: RuntimeWarning: invalid value encountered in matmul
  q_diag = np.diag(Q_half.T @ Q_half) + 0.1  # diagonal of Q^T Q + 0.1 I
# Solve with Ipopt
t0 = time.perf_counter()
res_ipopt = m_large.solve(nlp_solver="ipopt")
t_ipopt = time.perf_counter() - t0

# Solve with JAX IPM
t0 = time.perf_counter()
res_ipm = m_large.solve(nlp_solver="ipm")
t_ipm = time.perf_counter() - t0

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"{'JAX IPM':<12} {str(res_ipm.status):<12} {res_ipm.objective:>12.6f} {t_ipm:>10.4f}")
print(f"\nObjective difference: {abs(res_ipopt.objective - res_ipm.objective):.2e}")
******************************************************************************
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.9254
JAX IPM      optimal         -6.629880     0.0066

Objective difference: 0.00e+00

Both backends should find the same (or very close) optimal objective. Ipopt’s sparse factorization typically wins on larger problems, while the JAX IPM can be faster on small, dense problems and benefits from JIT warm-start when solving many similar problems in sequence.

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: 1.6423s

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

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 JAX-traceable functions by the NLPEvaluator. These functions are JIT-compiled for efficient evaluation.

  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. JAX computes exact gradients, Jacobians, and Hessians via reverse-mode AD. This is more accurate than finite differences and avoids the \(O(n)\) cost of forward-mode for gradients.

  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 JAX Hessian and flattens it to match hessianstructure().

from discopt._jax.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="ipm" (the JAX IPM) when:

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

  • You need JIT compilation benefits for solving many similar problems.

  • You want a pure-Python/JAX stack with no C library dependency.

  • You need differentiable optimization (backpropagating through the solve).

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      iteration_limit   0.00011320     0.5391
ipm        iteration_limit   0.00011320     0.4480

Ipopt solution: x1 = 0.989367, x2 = 0.978810

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, ripopt, cyipopt) and provides decision criteria for choosing the right one.