discopt Quickstart Tutorial

discopt Quickstart Tutorial#

This notebook walks you through the basics of formulating and solving optimization problems with discopt. discopt solves problems using spatial Branch and Bound [Belotti et al., 2013, Land and Doig, 1960] with McCormick relaxations [McCormick, 1976].

By the end you will know how to:

  1. Build a simple model with continuous, binary, and integer variables

  2. Add linear, nonlinear, and equality constraints

  3. Use mathematical functions (exp, log, sqrt, sin, cos)

  4. Inspect the solution (status, objective, variable values, timing)

  5. Load benchmark problems from .nl files

  6. Control solver behavior with options

Installation#

pip install discopt

discopt’s default NLP engine is POUNCE (pounce-solver), a pure-Rust port of the Ipopt interior-point method [Wächter and Biegler, 2006]. It is a core dependency, so nothing extra is needed. For the optional Ipopt backend, used here to cross-check POUNCE:

# macOS
brew install ipopt
pip install cyipopt

# Linux: install libipopt-dev (or build from source), then pip install cyipopt

import discopt.modeling as dm
import numpy as np

print("discopt loaded successfully")

discopt loaded successfully

1. Your First Model: Unconstrained Quadratic#

The simplest possible optimization problem. Minimize a quadratic with two variables:

\[\min_{x,y} \; (x - 3)^2 + (y + 1)^2, \quad x, y \in [-10, 10]\]

The minimum is obviously at \(x^* = 3\), \(y^* = -1\) with \(f^* = 0\).

m = dm.Model("first_model")
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()

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

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

The pattern is always the same:

  1. Create a Model

  2. Add variables with m.continuous(), m.binary(), or m.integer()

  3. Set the objective with m.minimize() or m.maximize()

  4. Optionally add constraints with m.subject_to()

  5. Call m.solve()

2. Adding Constraints#

Inequality constraints#

\[\min_{x,y} \; x^2 + y^2 \quad \text{s.t.} \quad x + y \geq 1, \quad x, y \in [-5, 5]\]

Without the constraint the minimum would be at the origin. The constraint pushes the solution to \(x^* = y^* = 0.5\) with \(f^* = 0.5\).

m = dm.Model("inequality_example")
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)

result = m.solve()
print(f"Status:    {result.status}")
print(f"Objective: {result.objective:.6f}  (expected: 0.5)")
print(f"x = {result.x['x']:.6f}  (expected: 0.5)")
print(f"y = {result.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)

Equality constraints#

Use == for equality constraints:

\[\min_{x,y} \; x^2 + y^2 \quad \text{s.t.} \quad x \cdot y = 1, \quad x, y \in [0.1, 5]\]

Known optimal: \(x^* = y^* = 1\), \(f^* = 2\).

m = dm.Model("equality_example")
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)

result = m.solve()
print(f"Status:    {result.status}")
print(f"Objective: {result.objective:.6f}  (expected: 2.0)")
print(f"x = {result.x['x']:.6f}  (expected: 1.0)")
print(f"y = {result.x['y']:.6f}  (expected: 1.0)")

Status:    optimal
Objective: 2.000000  (expected: 2.0)
x = 1.000000  (expected: 1.0)
y = 1.000000  (expected: 1.0)

Multiple constraints#

You can add as many constraints as needed:

\[\min_{x,y} \; (x - 2)^2 + (y - 2)^2 \quad \text{s.t.} \quad x + y \leq 3, \; x \geq 0, \; y \geq 0\]
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)

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

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

3. Variable Types#

discopt supports three variable types:

Method

Type

Domain

m.continuous("x", lb, ub)

Continuous

\([\text{lb}, \text{ub}]\)

m.binary("y")

Binary

\(\{0, 1\}\)

m.integer("n", lb, ub)

Integer

\(\{\text{lb}, \ldots, \text{ub}\}\)

When a model contains binary or integer variables, discopt solves it to global optimality. Spatial Branch & Bound [Land and Doig, 1960] is the general engine, but it is not the only route: at solve entry discopt classifies the model, and a convex MINLP (such as the convex MIQP below) is handed to the MIP–NLP family — outer approximation [Duran and Grossmann, 1986] — which certifies the optimum without ever opening a B&B node. That is why result.node_count reads 0 in the next two cells; the cells print result.algorithm_route so the route is visible rather than inferred; the first of them also solves a second time with solver="bb" to show the spatial B&B arm reaching the same optimum on the same model.

Binary variables#

\[\min_{x, y_1, y_2} \; (x - 3)^2 + 2y_1 + 3y_2 \quad \text{s.t.} \quad x \leq 5y_1, \; x \leq 4y_2, \; x \in [0, 5]\]

The binary variables act as “switches” that must be turned on for \(x\) to take a nonzero value.

m = dm.Model("binary_example")
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)

result = m.solve()
print(f"Status:    {result.status}")
print(f"Objective: {result.objective:.6f}  (expected: 5.0)")
print(f"x  = {result.x['x']:.6f}  (expected: 3.0)")
print(f"y1 = {result.x['y1']:.6f}  (expected: 1.0)")
print(f"y2 = {result.x['y2']:.6f}  (expected: 1.0)")
print(f"Nodes explored: {result.node_count}")
print(f"Route:          {result.algorithm_route}")

# The same model down the spatial-B&B arm, so "0 nodes" above is read as a route
# choice and not as "B&B did nothing". Both must agree on the optimum.
bb = m.solve(solver="bb")
print(f"\nsolver='bb':    obj={bb.objective:.6f}, nodes={bb.node_count}")
assert result.node_count == 0, "the convex MIQP route should not open a B&B node"
assert bb.node_count > 0, "the solver='bb' arm must actually branch"
assert abs(result.objective - bb.objective) < 1e-6, "the two routes must agree"

OA: generating OA cuts for 0 of 2 rows (2 convex; the rest are either nonconvex or already exact in the master)

Status:    optimal
Objective: 5.000000  (expected: 5.0)
x  = 3.000000  (expected: 3.0)
y1 = 1.000000  (expected: 1.0)
y2 = 1.000000  (expected: 1.0)
Nodes explored: 0
Route:          mip-nlp/oa: miqp certified convex at the root (DISCOPT_CONVEX_MINLP_ROUTE)

solver='bb':    obj=5.000000, nodes=5

Integer variables#

\[\min_{x, n} \; (x - 2.7)^2 + (n - 3)^2 \quad \text{s.t.} \quad x + n \leq 6, \; x \in [0, 5], \; n \in \{0, 1, \ldots, 5\}\]
m = dm.Model("integer_example")
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)

result = m.solve()
print(f"Status:    {result.status}")
print(f"Objective: {result.objective:.6f}  (expected: 0.0)")
print(f"x = {result.x['x']:.6f}  (expected: 2.7)")
print(f"n = {result.x['n']:.6f}  (expected: 3.0)")
print(f"Nodes explored: {result.node_count}")
print(f"Route:          {result.algorithm_route}")
assert abs(result.objective) < 1e-6, "this model's optimum is exactly 0"

OA: generating OA cuts for 0 of 1 rows (1 convex; the rest are either nonconvex or already exact in the master)

Status:    optimal
Objective: 0.000000  (expected: 0.0)
x = 2.700000  (expected: 2.7)
n = 3.000000  (expected: 3.0)
Nodes explored: 0
Route:          mip-nlp/oa: miqp certified convex at the root (DISCOPT_CONVEX_MINLP_ROUTE)

4. Mathematical Functions#

discopt provides these nonlinear functions in the dm namespace:

Function

Description

dm.exp(x)

Exponential \(e^x\)

dm.log(x)

Natural logarithm \(\ln(x)\)

dm.log2(x)

Base-2 logarithm

dm.log10(x)

Base-10 logarithm

dm.sqrt(x)

Square root \(\sqrt{x}\)

dm.sin(x)

Sine

dm.cos(x)

Cosine

dm.tan(x)

Tangent

dm.abs(x)

Absolute value

All of these are differentiated by the POUNCE AD tape in the Rust core. JAX is not on the default solve path — it is an optional dependency of the differentiable-solve and learned-relaxation subsystems only, and a plain solve() imports no jax module.

Example: exponential objective#

\[\min_{x,y} \; e^x + y^2 \quad \text{s.t.} \quad x + y \geq 1, \quad x, y \in [-2, 2]\]
m = dm.Model("exp_example")
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)

result = m.solve()
print(f"Status:    {result.status}")
print(f"Objective: {result.objective:.6f}  (expected: ~1.8395)")
print(f"x = {result.x['x']:.6f}")
print(f"y = {result.x['y']:.6f}")

Status:    optimal
Objective: 1.839484  (expected: ~1.8395)
x = 0.314923
y = 0.685077

Example: logarithmic objective#

\[\min_{x,y} \; -\log(x) - \log(y) \quad \text{s.t.} \quad x + y \leq 1, \quad x, y \in [0.01, 0.99]\]

Known optimal: \(x^* = y^* = 0.5\), \(f^* = 2\ln(2) \approx 1.3863\).

import math

m = dm.Model("log_example")
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)

expected = 2 * math.log(2)
result = m.solve()
print(f"Status:    {result.status}")
print(f"Objective: {result.objective:.6f}  (expected: {expected:.6f})")
print(f"x = {result.x['x']:.6f}  (expected: 0.5)")
print(f"y = {result.x['y']:.6f}  (expected: 0.5)")

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

5. Inspecting Results#

The SolveResult object returned by m.solve() contains:

Attribute

Description

result.status

"optimal", "feasible", "infeasible", "time_limit", "node_limit"

result.objective

Optimal objective value (float or None)

result.x

Solution dictionary: {var_name: value}

result.bound

Best dual bound (lower bound for minimization)

result.gap

Relative optimality gap

result.wall_time

Total solve time in seconds

result.node_count

Number of Branch & Bound nodes explored (0 for a pure NLP, and 0 on the MIP–NLP routes)

result.algorithm_route

Which route solved it, or None for plain spatial B&B

result.rust_time

Time spent in the Rust backend

result.python_time

Time spent in Python orchestration

result.pounce_time

Time inside POUNCE — a subset of rust_time, not a peer of it

result.jax_time

Time inside JAX — 0.0 on any default solve, since JAX is not on that path

Warning

rust_time and python_time partition the wall clock; pounce_time and jax_time are subsets of them. Subtracting a subset from the wall clock double-counts, so wall_time - jax_time - rust_time is not “the rest” — it is meaningless.

# Solve a small MINLP and inspect all result fields
m = dm.Model("inspect_demo")
x = m.continuous("x", lb=0, ub=5)
y = m.binary("y")
n = m.integer("n", lb=0, ub=3)

m.minimize(x**2 + 2 * y + n)
m.subject_to(x + n >= 2)
m.subject_to(x <= 4 * y)

result = m.solve()

print("=== Solve Result ===")
print(f"Status:       {result.status}")
print(f"Objective:    {result.objective}")
print(f"Bound:        {result.bound}")
print(f"Gap:          {result.gap}")
print(f"Wall time:    {result.wall_time:.4f}s")
print(f"Node count:   {result.node_count}")
print()
print("=== Layer Profiling ===")
# rust + python partition the wall clock; pounce/jax are subsets of that total.
print(f"Rust time:    {result.rust_time:.4f}s")
print(f"Python time:  {result.python_time:.4f}s")
print(f"  of which POUNCE: {result.pounce_time:.4f}s")
print(f"  of which JAX:    {result.jax_time:.4f}s  (0.0 -- JAX is not on the solve path)")
partition = result.rust_time + result.python_time
assert abs(partition - result.wall_time) < 0.05 * max(result.wall_time, 1e-3) + 1e-3, (
    f"rust+python ({partition:.4f}s) should account for the wall clock "
    f"({result.wall_time:.4f}s)"
)
print()
print("=== Variable Values ===")
for name, val in result.x.items():
    print(f"  {name} = {val}")

OA: generating OA cuts for 0 of 2 rows (2 convex; the rest are either nonconvex or already exact in the master)

=== Solve Result ===
Status:       optimal
Objective:    2.0
Bound:        2.0
Gap:          0.0
Wall time:    0.0075s
Node count:   0

=== Layer Profiling ===
Rust time:    0.0073s
Python time:  0.0001s
  of which POUNCE: 0.0073s
  of which JAX:    0.0000s  (0.0 -- JAX is not on the solve path)

=== Variable Values ===
  x = 3.334936609471964e-16
  y = 0.0
  n = 2.0

6. Loading .nl Files#

discopt can load standard AMPL .nl format files, commonly used for benchmark libraries like MINLPLib [Bussieck et al., 2003]. The Rust backend parses these files directly.

model = dm.from_nl("path/to/problem.nl")
result = model.solve()

Here we load two instances from the included MINLPLib test data and check each against its reference optimum. A missing corpus raises rather than printing skipping: a demo that quietly does nothing reads as a pass.

from pathlib import Path

# The notebook is executed from its own directory (jupyter-book) but may also be run
# from the repo root, so resolve against both rather than guessing one. The previous
# path, "../python/tests/data/minlplib", resolved to docs/python/... from here --- it
# never existed, and the exists() guard turned that into a silent "skipping".
NL_DIRS = [
    Path("../../python/tests/data/minlplib"),  # from docs/notebooks/
    Path("python/tests/data/minlplib"),  # from the repo root
]
nl_dir = next((d for d in NL_DIRS if d.is_dir()), None)
if nl_dir is None:
    raise FileNotFoundError(
        "MINLPLib test data not found in any of: " + ", ".join(str(d) for d in NL_DIRS)
    )

# Reference optima: nvs03 from python/tests/data/known_optima.toml, st_e13 from the
# in-repo reference in python/tests/test_monomial_lp_bound.py.
instances = {
    "st_e13": {"opt": 2.0, "type": "NLP"},
    "nvs03": {"opt": 16.0, "type": "MINLP"},
}

checked = 0
for name, info in instances.items():
    nl_path = nl_dir / f"{name}.nl"
    if not nl_path.exists():
        raise FileNotFoundError(nl_path)

    model = dm.from_nl(str(nl_path))
    # .nl models work with any backend; here we use Ipopt for demonstration.
    result = model.solve(time_limit=30, nlp_solver="ipopt")

    error = abs(result.objective - info["opt"])
    print(
        f"{name} ({info['type']}): status={result.status}, "
        f"obj={result.objective:.6f}, expected={info['opt']:.6f}, "
        f"error={error:.2e}, time={result.wall_time:.2f}s"
    )
    assert result.status == "optimal", f"{name}: {result.status}"
    assert error < 1e-6, f"{name}: {result.objective} != {info['opt']}"
    checked += 1

assert checked == len(instances), f"only {checked} instance(s) were checked"
print(f"\n{checked} instance(s) solved and checked against their reference optima.")

OA: generating OA cuts for 1 of 2 rows (2 convex; the rest are either nonconvex or already exact in the master)

st_e13 (NLP): status=optimal, obj=2.000000, expected=2.000000, error=2.52e-09, time=0.01s

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

nvs03 (MINLP): status=optimal, obj=16.000000, expected=16.000000, error=0.00e+00, time=0.01s

2 instance(s) solved and checked against their reference optima.

7. Solver Options#

The solve() method accepts several options to control solver behavior:

Parameter

Default

Description

time_limit

3600

Wall-clock time limit in seconds

max_nodes

100000

Maximum Branch & Bound nodes

gap_tolerance

1e-4

Relative optimality gap for termination

nlp_solver

"pounce"

NLP backend: "pounce" (the Rust Ipopt port [Wächter and Biegler, 2006]), "ipopt" (cyipopt), or "simplex". "ipm"/"sparse_ipm" are back-compat aliases that resolve to POUNCE [Nocedal and Wright, 2006] — the pure-JAX IPM they once named is retired

strategy

"best_first"

Node selection: "best_first" or "depth_first"

batch_size

16

Number of nodes solved per B&B iteration

partitions

0

Piecewise McCormick partitions [Bergamini et al., 2005, Castro, 2015] (0=standard relaxation)

cutting_planes

False

Enable OA [Duran and Grossmann, 1986, Fletcher and Leyffer, 1994] / RLT [Sherali and Adams, 1990] cut generation

 
# Solve with a tight time limit and node limit
m = dm.Model("options_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 + 3 * y)
m.subject_to(x1 + x2 >= 1)
m.subject_to(x1 <= 4 * y)

result = m.solve(
    time_limit=10,
    max_nodes=50,
    gap_tolerance=1e-6,
)

print(f"Status:    {result.status}")
print(f"Objective: {result.objective}")
print(f"Gap:       {result.gap}")
print(f"Nodes:     {result.node_count}")
print(f"Time:      {result.wall_time:.3f}s")

OA: generating OA cuts for 0 of 2 rows (2 convex; the rest are either nonconvex or already exact in the master)

Status:    optimal
Objective: 1.000000000001402
Gap:       0.0
Nodes:     0
Time:      0.010s

8. Array Variables#

For problems with many variables of the same type, use the shape parameter to create array variables. Use dm.sum() for aggregation.

# Simple portfolio-like problem: minimize variance subject to return target
n_assets = 4
returns = np.array([0.12, 0.10, 0.07, 0.03])
target_return = 0.08

m = dm.Model("portfolio")
w = m.continuous("w", shape=(n_assets,), lb=0, ub=1)

# Minimize sum of squared weights (simplified variance proxy)
m.minimize(dm.sum(lambda i: w[i] ** 2, over=range(n_assets)))

# Weights sum to 1
m.subject_to(dm.sum(w) == 1)

# Minimum return constraint
m.subject_to(dm.sum(lambda i: returns[i] * w[i], over=range(n_assets)) >= target_return)

result = m.solve()
print(f"Status: {result.status}")
print(f"Objective (variance proxy): {result.objective:.6f}")
print(f"Weights: {result.x['w']}")
actual_return = np.dot(returns, result.x["w"])
print(f"Portfolio return: {actual_return:.4f} (target: {target_return})")

Status: optimal
Objective (variance proxy): 0.250000
Weights: [0.25002565 0.25001283 0.24999359 0.24996793]
Portfolio return: 0.0800 (target: 0.08)

9. Model Summary#

You can inspect a model before solving with print(m) or m.summary().

m = dm.Model("demo")
x = m.continuous("x", shape=(3,), lb=0, ub=10)
y = m.binary("y", shape=(2,))
n = m.integer("n", lb=0, ub=100)
m.minimize(dm.sum(x) + 10 * dm.sum(y) + n)
m.subject_to(x[0] + x[1] >= 1)
m.subject_to(dm.exp(x[2]) <= 5)

print(m.summary())

Model: demo
  Variables: 6 (3 continuous, 3 integer/binary)
  Constraints: 2
  Objective: minimize ((sum(x[3]) + (10 * sum(y[2]))) + n)
  Parameters: 0

Next Steps#

  • More examples: See notebooks/minlp_examples.ipynb for a comprehensive set of NLP and MINLP problems

  • Advanced features: See notebooks/advanced_features.ipynb for differentiable optimization, GPU acceleration, piecewise McCormick relaxations, and more

  • Benchmarks: Run python discopt_benchmarks/run_benchmarks.py --suite smoke to benchmark discopt against other solvers

  • API reference: import discopt.modeling as dm; help(dm.Model) for full API documentation