Modeling Guide: Expressive vs. Fast Construction

Modeling Guide: Expressive vs. Fast Construction#

discopt provides two complementary APIs for building optimization models:

Expressive API

Fast API

Syntax

Operator overloading (x + y <= 3)

Matrix/vector (A, b, c, Q)

Best for

Prototyping, nonlinear, readable models

Large-scale LP/QP/MILP/MIQP

Construction speed

O(n) Python objects per constraint

Single Rust call for all constraints

Nonlinear support

Full (exp, log, sin, bilinear, …)

Linear and quadratic only

Hybrid

Yes — mix both in one model

Yes — mix both in one model

Both APIs feed into the same solver pipeline: automatic problem classification, algebraic coefficient extraction, pure-JAX interior-point method, and Branch & Bound.

This notebook covers:

  1. The Expressive API — syntax, functions, and examples

  2. The Fast API — bulk construction for large models

  3. Hybrid mode — combining both in one model

  4. Performance comparison — when the fast API matters

  5. Decision guide — which API to choose

import os

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

import time

import discopt.modeling as dm
import numpy as np
import scipy.sparse as sp

Part 1: The Expressive API#

The expressive API uses Python operator overloading to build an expression DAG (directed acyclic graph). Every arithmetic operation creates a node; the DAG is later compiled to JAX for autodiff and to Rust for structure detection.

Pattern#

m = dm.Model("name")
x = m.continuous("x", lb=0, ub=10)     # 1. Variables
m.minimize(x**2 + 1)                    # 2. Objective
m.subject_to(x >= 0.5)                  # 3. Constraints
result = m.solve()                       # 4. Solve

1.1 Variables#

Three variable types are available:

Method

Domain

Example

m.continuous(name, lb, ub, shape)

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

Flow rates, temperatures

m.binary(name, shape)

\(\{0, 1\}\)

On/off decisions

m.integer(name, lb, ub, shape)

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

Batch counts

 
m = dm.Model("variable_demo")

# Scalar variables
T = m.continuous("T", lb=300, ub=600)  # temperature (K)
on = m.binary("on")  # reactor on/off
n_batches = m.integer("n_batches", lb=0, ub=20)

# Array variables (shape tuple)
x = m.continuous("x", shape=(5,), lb=0, ub=1)  # 5-element vector
y = m.binary("y", shape=(3, 4))  # 3x4 binary matrix

print(f"Scalar T: shape={T.shape}, bounds=[{T.lb}, {T.ub}]")
print(f"Array x:  shape={x.shape}, {x.size} elements")
print(f"Matrix y: shape={y.shape}, {y.size} elements")

Scalar T: shape=(), bounds=[300.0, 600.0]
Array x:  shape=(5,), 5 elements
Matrix y: shape=(3, 4), 12 elements

1.2 Building Expressions#

Arithmetic on variables creates expression nodes in the DAG. All standard Python operators are supported, plus mathematical functions from the dm namespace.

m = dm.Model("expression_demo")
x = m.continuous("x", shape=(3,), lb=0, ub=5)

# Arithmetic operators
linear = 2 * x[0] + 3 * x[1] - x[2]  # linear combination
quadratic = x[0] ** 2 + x[1] * x[2]  # quadratic / bilinear
rational = x[0] / (1 + x[1])  # rational expression

# Mathematical functions
nonlinear = dm.exp(-x[0]) + dm.log(1 + x[1]) + dm.sqrt(x[2])
trig = dm.sin(x[0]) * dm.cos(x[1])

# Aggregations
total = dm.sum(x)  # sum of array
indexed = dm.sum(lambda i: x[i] ** 2, over=range(3))  # indexed sum

print(f"linear:    {linear}")
print(f"quadratic: {quadratic}")
print(f"nonlinear: {nonlinear}")

linear:    (((2 * x[3][0]) + (3 * x[3][1])) - x[3][2])
quadratic: ((x[3][0] ** 2) + (x[3][1] * x[3][2]))
nonlinear: ((exp(neg(x[3][0])) + log((1 + x[3][1]))) + sqrt(x[3][2]))

1.3 Constraints#

Comparisons on expressions produce Constraint objects:

m.subject_to(expr1 <= expr2)     # inequality
m.subject_to(expr1 >= expr2)     # inequality
m.subject_to(expr1 == expr2)     # equality

Constraints can be named for diagnostics:

m = dm.Model("constraint_demo")
x = m.continuous("x", shape=(3,), lb=0, ub=10)

m.minimize(dm.sum(x))

# Named constraints
m.subject_to(x[0] + x[1] <= 5, name="capacity")
m.subject_to(x[0] * x[1] >= 1, name="quality")
m.subject_to(dm.sum(x) == 7, name="balance")

# Loop-based constraints
for i in range(3):
    m.subject_to(x[i] >= 0.5, name=f"min_flow_{i}")

print(m.summary())

Model: constraint_demo
  Variables: 3 (3 continuous, 0 integer/binary)
  Constraints: 6
  Objective: minimize sum(x[3])
  Parameters: 0

1.4 Example: Geometric Programming (via log transform)#

Minimize the surface area of a box with fixed volume \(V = 1\):

\[ \min_{w, d, h} \; 2(wd + wh + dh) \quad \text{s.t.} \quad wdh = 1, \; w, d, h \in [0.1, 10] \]

The optimal solution is a cube: \(w^* = d^* = h^* = 1\), surface area \(= 6\).

m = dm.Model("box")
w = m.continuous("w", lb=0.1, ub=10)
d = m.continuous("d", lb=0.1, ub=10)
h = m.continuous("h", lb=0.1, ub=10)

m.minimize(2 * (w * d + w * h + d * h))
m.subject_to(w * d * h == 1)  # volume constraint

result = m.solve()
print(f"Status: {result.status}")
print(f"Surface area: {result.objective:.4f}  (expected: 6.0)")
print(f"w={result.x['w']:.4f}, d={result.x['d']:.4f}, h={result.x['h']:.4f}")

******************************************************************************
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
Surface area: 6.0000  (expected: 6.0)
w=1.0000, d=1.0000, h=1.0000

1.5 Example: Reactor Design (nonlinear MINLP)#

Select a reactor configuration and optimize temperature to maximize yield. The reaction rate follows Arrhenius kinetics:

\[ \max_{T, y} \; y \cdot k_0 \, e^{-E_a / (RT)} \quad \text{s.t.} \quad 300 \leq T \leq 600, \; y \in \{0, 1\} \]

where \(y\) is an on/off indicator for the reactor.

k0 = 1e6  # pre-exponential factor
Ea = 5000.0  # activation energy / R (in K)

m = dm.Model("reactor")
T = m.continuous("T", lb=300, ub=600)
y = m.binary("y")  # reactor on/off

# Maximize rate (= minimize negative rate)
rate = k0 * dm.exp(-Ea / T)
m.minimize(-(y * rate))

# Reactor must be on for temperature to matter
m.subject_to(T <= 300 + 300 * y)  # T = 300 when y=0

result = m.solve()
print(f"Status: {result.status}")
print(f"Reactor on: {result.x['y']:.0f}")
print(f"Temperature: {result.x['T']:.1f} K")
print(f"Rate: {-result.objective:.2f}")

Status: optimal
Reactor on: 1
Temperature: 600.0 K
Rate: 240.37

1.6 Example: Portfolio Optimization (QP via expressions)#

Minimize portfolio variance subject to a target return:

\[ \min_w \; w^\top \Sigma\, w \quad \text{s.t.} \quad \mu^\top w \geq r_{\min}, \; \mathbf{1}^\top w = 1, \; w \geq 0 \]
n_assets = 4
mu = np.array([0.12, 0.10, 0.07, 0.03])  # expected returns
# Covariance matrix
Sigma = np.array(
    [
        [0.04, 0.006, 0.002, 0.0],
        [0.006, 0.025, 0.004, 0.001],
        [0.002, 0.004, 0.01, 0.002],
        [0.0, 0.001, 0.002, 0.005],
    ]
)
r_min = 0.08

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

# Variance = w' Sigma w (expressed element-wise)
variance = dm.sum(
    lambda i: dm.sum(lambda j: Sigma[i, j] * w[i] * w[j], over=range(n_assets)),
    over=range(n_assets),
)

m.minimize(variance)
m.subject_to(dm.sum(lambda i: mu[i] * w[i], over=range(n_assets)) >= r_min)
m.subject_to(dm.sum(w) == 1)

result = m.solve()
print(f"Status: {result.status}")
print(f"Variance: {result.objective:.6f}")
print(f"Weights: {np.round(result.x['w'], 4)}")
print(f"Return: {np.dot(mu, result.x['w']):.4f}  (target >= {r_min})")

Status: optimal
Variance: 0.006602
Weights: [0.2065 0.2187 0.4027 0.1721]
Return: 0.0800  (target >= 0.08)

Part 2: The Fast Construction API#

For large-scale LP and QP models, the expressive API becomes a bottleneck: creating thousands of Python expression objects is slow. The fast API bypasses expression objects entirely, building constraints directly into the Rust expression arena via a single PyO3 call.

Methods#

Method

Purpose

m.add_linear_constraints(A, x, sense, b)

Bulk linear constraints: \(Ax \leq b\)

m.add_linear_objective(c, x)

Linear objective: \(\min\; c^\top x\)

m.add_quadratic_objective(Q, c, x)

Quadratic objective: \(\min\; \tfrac{1}{2} x^\top Q x + c^\top x\)

Key points:

  • A can be a scipy sparse matrix (any format) or dense numpy array — automatically converted to CSR

  • sense is "<=", "==", or ">="

  • b can be a vector or a scalar (broadcast to all rows)

  • Q must be symmetric; uses the convention \(f = \tfrac{1}{2} x^\top Q x + c^\top x\)

2.1 Example: Large LP (Transportation)#

Ship goods from 50 warehouses to 100 customers. This creates 150 constraints that we add in a single call.

n_wh, n_cust = 50, 100
n_vars = n_wh * n_cust  # 5000 decision variables

np.random.seed(42)
supply = np.random.uniform(50, 200, n_wh)
demand = np.random.uniform(10, 50, n_cust)
# Scale demand to match total supply
demand *= supply.sum() / demand.sum()
cost = np.random.uniform(1, 20, (n_wh, n_cust))

m = dm.Model("transport_fast")
x = m.continuous("x", shape=(n_vars,), lb=0)

# --- Build constraint matrix ---
# Supply: sum_j x[i,j] <= supply[i]  for each warehouse i
# Demand: sum_i x[i,j] >= demand[j]  for each customer j
rows_supply, cols_supply, data_supply = [], [], []
for i in range(n_wh):
    for j in range(n_cust):
        rows_supply.append(i)
        cols_supply.append(i * n_cust + j)
        data_supply.append(1.0)

A_supply = sp.csr_matrix((data_supply, (rows_supply, cols_supply)), shape=(n_wh, n_vars))

rows_demand, cols_demand, data_demand = [], [], []
for j in range(n_cust):
    for i in range(n_wh):
        rows_demand.append(j)
        cols_demand.append(i * n_cust + j)
        data_demand.append(1.0)

A_demand = sp.csr_matrix((data_demand, (rows_demand, cols_demand)), shape=(n_cust, n_vars))

# Add all constraints in two bulk calls
t0 = time.perf_counter()
m.add_linear_constraints(A_supply, x, "<=", supply, name="supply")
m.add_linear_constraints(A_demand, x, ">=", demand, name="demand")
m.add_linear_objective(cost.flatten(), x)
t_build = time.perf_counter() - t0

print(f"Model: {n_vars} variables, {n_wh + n_cust} constraints")
print(f"Build time: {t_build * 1000:.1f} ms")

result = m.solve()
print(f"Status: {result.status}")
print(f"Total cost: {result.objective:.2f}")
print(f"Solve time: {result.wall_time:.3f}s")

Model: 5000 variables, 150 constraints
Build time: 0.6 ms

/Users/jkitchin/Dropbox/projects/discopt/python/discopt/modeling/core.py:1903: UserWarning: Variables with very large or infinite bounds: x[0] (lb=0, ub=1e+20), x[1] (lb=0, ub=1e+20), x[2] (lb=0, ub=1e+20), x[3] (lb=0, ub=1e+20), x[4] (lb=0, ub=1e+20). NLP solvers may fail (NaN, iteration_limit) when bounds exceed ~1e15. Add tighter explicit bounds, e.g. m.continuous('x', lb=0, ub=1000).
  result = solve_model(

Status: optimal
Total cost: 8761.51
Solve time: 1.532s

2.2 Example: Ridge Regression as QP#

Given data matrix \(X \in \mathbb{R}^{m \times n}\) and targets \(y \in \mathbb{R}^m\), ridge regression solves:

\[ \min_\beta \; \tfrac{1}{2} \|X\beta - y\|^2 + \tfrac{\lambda}{2} \|\beta\|^2 = \min_\beta \; \tfrac{1}{2} \beta^\top (X^\top X + \lambda I) \beta - (X^\top y)^\top \beta + \text{const} \]

This is a QP with \(Q = X^\top X + \lambda I\) and \(c = -X^\top y\).

# Generate synthetic regression data
np.random.seed(0)
m_samples, n_features = 200, 20
X_data = np.random.randn(m_samples, n_features)
beta_true = np.random.randn(n_features)
y_data = X_data @ beta_true + 0.5 * np.random.randn(m_samples)

lam = 1.0  # regularization strength

# QP coefficients
Q = sp.csr_matrix(X_data.T @ X_data + lam * np.eye(n_features))
c_vec = -(X_data.T @ y_data)

# Build and solve
mod = dm.Model("ridge")
beta = mod.continuous("beta", shape=(n_features,), lb=-100, ub=100)
mod.add_quadratic_objective(Q, c_vec, beta)

result = mod.solve()

# Compare with closed-form solution
beta_closed = np.linalg.solve(X_data.T @ X_data + lam * np.eye(n_features), X_data.T @ y_data)

print(f"Status: {result.status}")
print(f"Max |discopt - closed_form|: {np.max(np.abs(result.x['beta'] - beta_closed)):.2e}")
print(f"Solve time: {result.wall_time:.3f}s")

Status: optimal
Max |discopt - closed_form|: 1.27e-09
Solve time: 0.002s

/var/folders/gq/k1kgbl7n539_4dl1md8x3jt80000gn/T/ipykernel_44057/2005627179.py:6: RuntimeWarning: divide by zero encountered in matmul
  y_data = X_data @ beta_true + 0.5 * np.random.randn(m_samples)
/var/folders/gq/k1kgbl7n539_4dl1md8x3jt80000gn/T/ipykernel_44057/2005627179.py:6: RuntimeWarning: overflow encountered in matmul
  y_data = X_data @ beta_true + 0.5 * np.random.randn(m_samples)
/var/folders/gq/k1kgbl7n539_4dl1md8x3jt80000gn/T/ipykernel_44057/2005627179.py:6: RuntimeWarning: invalid value encountered in matmul
  y_data = X_data @ beta_true + 0.5 * np.random.randn(m_samples)
/var/folders/gq/k1kgbl7n539_4dl1md8x3jt80000gn/T/ipykernel_44057/2005627179.py:11: RuntimeWarning: divide by zero encountered in matmul
  Q = sp.csr_matrix(X_data.T @ X_data + lam * np.eye(n_features))
/var/folders/gq/k1kgbl7n539_4dl1md8x3jt80000gn/T/ipykernel_44057/2005627179.py:11: RuntimeWarning: overflow encountered in matmul
  Q = sp.csr_matrix(X_data.T @ X_data + lam * np.eye(n_features))
/var/folders/gq/k1kgbl7n539_4dl1md8x3jt80000gn/T/ipykernel_44057/2005627179.py:11: RuntimeWarning: invalid value encountered in matmul
  Q = sp.csr_matrix(X_data.T @ X_data + lam * np.eye(n_features))
/var/folders/gq/k1kgbl7n539_4dl1md8x3jt80000gn/T/ipykernel_44057/2005627179.py:12: RuntimeWarning: divide by zero encountered in matmul
  c_vec = -(X_data.T @ y_data)
/var/folders/gq/k1kgbl7n539_4dl1md8x3jt80000gn/T/ipykernel_44057/2005627179.py:12: RuntimeWarning: overflow encountered in matmul
  c_vec = -(X_data.T @ y_data)
/var/folders/gq/k1kgbl7n539_4dl1md8x3jt80000gn/T/ipykernel_44057/2005627179.py:12: RuntimeWarning: invalid value encountered in matmul
  c_vec = -(X_data.T @ y_data)
/var/folders/gq/k1kgbl7n539_4dl1md8x3jt80000gn/T/ipykernel_44057/2005627179.py:22: RuntimeWarning: divide by zero encountered in matmul
  beta_closed = np.linalg.solve(X_data.T @ X_data + lam * np.eye(n_features), X_data.T @ y_data)
/var/folders/gq/k1kgbl7n539_4dl1md8x3jt80000gn/T/ipykernel_44057/2005627179.py:22: RuntimeWarning: overflow encountered in matmul
  beta_closed = np.linalg.solve(X_data.T @ X_data + lam * np.eye(n_features), X_data.T @ y_data)
/var/folders/gq/k1kgbl7n539_4dl1md8x3jt80000gn/T/ipykernel_44057/2005627179.py:22: RuntimeWarning: invalid value encountered in matmul
  beta_closed = np.linalg.solve(X_data.T @ X_data + lam * np.eye(n_features), X_data.T @ y_data)

2.3 Example: Facility Location (MILP)#

Choose \(p\) facility locations from \(n\) candidates to minimize total weighted distance to \(m\) customers.

\[ \min_{x, y} \; \sum_{i,j} d_{ij}\, x_{ij} \quad \text{s.t.} \quad \sum_j x_{ij} = 1 \; \forall i, \; x_{ij} \leq y_j \; \forall i,j, \; \sum_j y_j = p, \; x_{ij} \geq 0, \; y_j \in \{0,1\} \]
n_facilities = 10
n_customers = 30
p = 3  # number to open

np.random.seed(7)
# Random locations
fac_loc = np.random.rand(n_facilities, 2) * 100
cust_loc = np.random.rand(n_customers, 2) * 100
dist = np.sqrt(((cust_loc[:, None, :] - fac_loc[None, :, :]) ** 2).sum(axis=2))

n_x = n_customers * n_facilities  # assignment variables
n_y = n_facilities  # open/close variables
n_total = n_x + n_y

mod = dm.Model("facility")
# x_ij: fraction of customer i served by facility j (continuous)
x = mod.continuous("x", shape=(n_x,), lb=0, ub=1)
# y_j: whether facility j is open (binary)
y = mod.binary("y", shape=(n_y,))

# Objective: minimize total distance (linear in x)
mod.add_linear_objective(dist.flatten(), x)

# Assignment constraints: sum_j x_ij = 1 for each customer
rows_a, cols_a, data_a = [], [], []
for i in range(n_customers):
    for j in range(n_facilities):
        rows_a.append(i)
        cols_a.append(i * n_facilities + j)
        data_a.append(1.0)
A_assign = sp.csr_matrix((data_a, (rows_a, cols_a)), shape=(n_customers, n_x))
mod.add_linear_constraints(A_assign, x, "==", 1.0, name="assign")

# Linking: x_ij <= y_j (via expressions, since they involve two variable groups)
for i in range(n_customers):
    for j in range(n_facilities):
        mod.subject_to(x[i * n_facilities + j] <= y[j])

# Cardinality: sum_j y_j = p
mod.subject_to(dm.sum(y) == p, name="cardinality")

result = mod.solve(time_limit=30)
print(f"Status: {result.status}")
print(f"Total distance: {result.objective:.2f}")
print(f"Nodes: {result.node_count}, Time: {result.wall_time:.2f}s")

opened = [j for j in range(n_facilities) if result.x["y"][j] > 0.5]
print(f"Opened facilities: {opened}")

Status: optimal
Total distance: 708.58
Nodes: 1, Time: 0.05s
Opened facilities: [1, 2, 4]

Part 3: Hybrid Mode#

The real power of discopt is that you can mix both APIs in a single model. Use the fast API for the bulk linear/quadratic structure, and the expressive API for the few nonlinear constraints.

Example: Blending with Nonlinear Quality#

A refinery has 500 linear flow-balance constraints (perfect for the fast API) plus a few nonlinear quality blending constraints (need expressions).

n_streams = 50

mod = dm.Model("hybrid_blend")
flow = mod.continuous("flow", shape=(n_streams,), lb=0, ub=1000)

# --- Bulk linear constraints: flow balance ---
# Random sparse flow network
np.random.seed(99)
A_balance = sp.random(40, n_streams, density=0.15, format="csr", random_state=99)
b_balance = np.random.uniform(10, 100, 40)
mod.add_linear_constraints(A_balance, flow, "<=", b_balance, name="balance")

# Linear objective
cost_vec = np.random.uniform(1, 10, n_streams)
mod.add_linear_objective(cost_vec, flow)

# --- A few nonlinear constraints via expressions ---
# Nonlinear blending: quality is a ratio
mod.subject_to(
    (3.0 * flow[0] + 5.0 * flow[1]) / (flow[0] + flow[1] + 0.01) >= 4.0,
    name="quality_blend",
)
# Exponential decay constraint
mod.subject_to(
    dm.exp(-flow[2] / 100) <= 0.5,
    name="conversion",
)

print(mod.summary())

result = mod.solve(time_limit=10)
print(f"\nStatus: {result.status}")
print(f"Objective: {result.objective:.4f}")
print(f"Solve time: {result.wall_time:.3f}s")

Model: hybrid_blend
  Variables: 50 (50 continuous, 0 integer/binary)
  Constraints: 2
  Objective: minimize 0
  Parameters: 0

Status: optimal
Objective: 0.0000
Solve time: 0.115s

Part 4: Performance Comparison#

The fast API is designed for model construction speed. Both APIs produce models that solve at the same speed, but building the model object differs dramatically for large instances.

sizes = [100, 500, 1000, 5000]
results_table = []

for n_con in sizes:
    n_var = 50
    np.random.seed(42)
    A = sp.random(n_con, n_var, density=0.1, format="csr", random_state=42)
    b = np.ones(n_con)
    c = np.ones(n_var)

    # --- Fast API ---
    t0 = time.perf_counter()
    m_fast = dm.Model("fast")
    x_fast = m_fast.continuous("x", shape=(n_var,), lb=0, ub=100)
    m_fast.add_linear_constraints(A, x_fast, "<=", b)
    m_fast.add_linear_objective(c, x_fast)
    t_fast = time.perf_counter() - t0

    # --- Expressive API ---
    A_dense = A.toarray()
    t0 = time.perf_counter()
    m_expr = dm.Model("expr")
    x_expr = m_expr.continuous("x", shape=(n_var,), lb=0, ub=100)
    for row in range(n_con):
        lhs = dm.sum(lambda j: A_dense[row, j] * x_expr[j], over=range(n_var))
        m_expr.subject_to(lhs <= b[row])
    m_expr.minimize(dm.sum(lambda j: c[j] * x_expr[j], over=range(n_var)))
    t_expr = time.perf_counter() - t0

    speedup = t_expr / t_fast if t_fast > 0 else float("inf")
    results_table.append((n_con, t_fast, t_expr, speedup))

print(f"{'Constraints':>12s} {'Fast (ms)':>10s} {'Expressive (ms)':>16s} {'Speedup':>8s}")
print("-" * 50)
for n_con, t_fast, t_expr, speedup in results_table:
    print(f"{n_con:>12d} {t_fast * 1000:>10.1f} {t_expr * 1000:>16.1f} {speedup:>7.0f}x")

 Constraints  Fast (ms)  Expressive (ms)  Speedup
--------------------------------------------------
         100        0.1              6.5      77x
         500        0.1             63.4     521x
        1000        0.2            144.2     794x
        5000        0.6            593.6     961x

The fast API provides orders of magnitude speedup for model construction. This matters when:

  • Building models inside optimization loops (e.g., parametric sweeps)

  • Generating models from large datasets

  • Production systems where latency is critical

For small models (< 100 constraints), the difference is negligible — use whichever API is more readable.

Part 5: Decision Guide#

Use the Expressive API when:#

  • Nonlinear functions are needed (exp, log, sin, bilinear products)

  • Readability matters (publications, teaching, stakeholder review)

  • Prototyping a new formulation

  • The model is small to medium (< 1000 constraints)

  • You need GDP features (if_then, either_or, sos1, sos2)

Use the Fast API when:#

  • The model is large-scale LP/QP/MILP/MIQP (1000+ constraints)

  • Data comes from numerical sources (numpy arrays, scipy sparse matrices)

  • Construction time is a bottleneck (parametric sweeps, real-time systems)

  • The formulation is naturally matrix-oriented (network flows, regression)

Use Hybrid mode when:#

  • Most constraints are linear (fast API) but a few are nonlinear (expressions)

  • You want the best of both: construction speed + modeling flexibility

Quick reference#

Task

Expressive

Fast

min x^2 + y^2

m.minimize(x**2 + y**2)

m.add_quadratic_objective(I, zeros, v)

Ax <= b (1000 rows)

Loop + m.subject_to()

m.add_linear_constraints(A, x, '<=', b)

exp(x) + log(y) <= 10

m.subject_to(dm.exp(x) + dm.log(y) <= 10)

Not supported

x * y >= 1 (bilinear)

m.subject_to(x * y >= 1)

Not supported

GDP / indicator

m.if_then(y, [...])

Not supported

API Reference Summary#

Model creation#

import discopt.modeling as dm
m = dm.Model("name")

Variables#

x = m.continuous("x", shape=(n,), lb=0, ub=100)
y = m.binary("y", shape=(k,))
z = m.integer("z", lb=0, ub=10)

Expressive objective and constraints#

m.minimize(expr)           # or m.maximize(expr)
m.subject_to(expr <= rhs, name="label")

Mathematical functions#

dm.exp(x), dm.log(x), dm.sqrt(x)
dm.sin(x), dm.cos(x), dm.tan(x)
dm.abs(x), dm.log2(x), dm.log10(x)
dm.sum(x), dm.sum(lambda i: ..., over=range(n))

Fast construction#

m.add_linear_constraints(A, x, "<=", b, name="prefix")
m.add_linear_objective(c, x, constant=0.0, sense="minimize")
m.add_quadratic_objective(Q, c, x, constant=0.0, sense="minimize")

Solving#

result = m.solve(
    time_limit=3600,
    max_nodes=100000,
    gap_tolerance=1e-4,
    nlp_solver="ipm",        # "ipm", "ipopt", or "ripopt"
)
print(result.status, result.objective, result.x)