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 ships with a pure-JAX interior-point method (IPM) as the default NLP backend. For the optional Ipopt backend:

# macOS
brew install ipopt
pip install cyipopt

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

import os

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

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)")

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

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 uses spatial Branch & Bound [Land and Doig, 1960] to find the global optimum.

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}")

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: 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}")

Status:    optimal
Objective: -0.000000  (expected: 0.0)
x = 2.700000  (expected: 2.7)
n = 3.000000  (expected: 3.0)
Nodes explored: 1

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 automatically differentiated by JAX.

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 pure NLP)

result.rust_time

Time spent in Rust backend

result.jax_time

Time spent in JAX (AD, relaxation evaluation)

result.python_time

Time spent in Python orchestration

 
# 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 ===")
print(f"Rust time:    {result.rust_time:.4f}s")
print(f"JAX time:     {result.jax_time:.4f}s")
print(f"Python time:  {result.python_time:.4f}s")
print()
print("=== Variable Values ===")
for name, val in result.x.items():
    print(f"  {name} = {val}")

=== Solve Result ===
Status:       optimal
Objective:    2.0000000000030003
Bound:        2.499999999999729
Gap:          0.0
Wall time:    0.0024s
Node count:   5

=== Layer Profiling ===
Rust time:    0.0000s
JAX time:     0.0016s
Python time:  0.0008s

=== Variable Values ===
  x = 1e-12
  y = 1e-12
  n = 2.000000000001

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.

from pathlib import Path

nl_dir = Path("../python/tests/data/minlplib")

instances = {
    "st_e13": {"opt": 2.0, "type": "NLP"},
    "nvs03": {"opt": 16.0, "type": "MINLP"},
}

for name, info in instances.items():
    nl_path = nl_dir / f"{name}.nl"
    if not nl_path.exists():
        print(f"{name}: .nl file not found, skipping")
        continue

    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"]) if result.objective is not None else float("inf")
    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"
    )

st_e13: .nl file not found, skipping
nvs03: .nl file not found, skipping

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

"ipm"

NLP backend: "ipm" (JAX IPM [Nocedal and Wright, 2006]), "ipopt" (Ipopt [Wächter and Biegler, 2006]), or "ripopt" (Rust)

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)

branching_policy

"fractional"

Variable selection: "fractional" or "gnn" [Gasse et al., 2019]

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")

Status:    optimal
Objective: 1.000000000066689
Gap:       0.0
Nodes:     3
Time:      0.002s

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.25 0.25 0.25 0.25]
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, GNN branching, 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