discopt MINLP Examples#
This notebook demonstrates solving Mixed-Integer Nonlinear Programs (MINLPs) using the discopt modeling API. discopt solves MINLPs via spatial Branch and Bound [Belotti et al., 2013, Land and Doig, 1960] with McCormick relaxations. We show:
Pure NLP problems (continuous variables only) solved via the default IPM backend
MINLP problems (with binary/integer variables) solved via spatial Branch & Bound
MINLPLib instances loaded from
.nlfiles
For each example, we verify the solution against known optimal values.
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. Unconstrained Quadratic (NLP)#
\[\min_{x,y} \; (x - 2)^2 + (y + 1)^2, \quad x, y \in [-5, 5]\]
Known optimal: \(x^* = 2, \; y^* = -1, \; f^* = 0\)
m = dm.Model("unconstrained_quadratic")
x = m.continuous("x", lb=-5, ub=5)
y = m.continuous("y", lb=-5, ub=5)
m.minimize((x - 2) ** 2 + (y + 1) ** 2)
result = m.solve()
print(f"Status: {result.status}")
print(f"Objective: {result.objective:.8f} (expected: 0.0)")
print(f"x = {result.value(x):.6f} (expected: 2.0)")
print(f"y = {result.value(y):.6f} (expected: -1.0)")
print(f"Time: {result.wall_time:.3f}s")
Status: optimal
Objective: 0.00000000 (expected: 0.0)
x = 2.000000 (expected: 2.0)
y = -1.000000 (expected: -1.0)
Time: 0.453s
2. Rosenbrock Function (NLP)#
\[\min_{x,y} \; (1 - x)^2 + 100(y - x^2)^2, \quad x, y \in [-5, 5]\]
The classic banana-shaped valley. Known optimal: \(x^* = 1, \; y^* = 1, \; f^* = 0\)
m = dm.Model("rosenbrock")
x = m.continuous("x", lb=-5, ub=5)
y = m.continuous("y", lb=-5, ub=5)
m.minimize((1 - x) ** 2 + 100 * (y - x**2) ** 2)
result = m.solve()
print(f"Status: {result.status}")
print(f"Objective: {result.objective:.8f} (expected: 0.0)")
print(f"x = {result.value(x):.6f} (expected: 1.0)")
print(f"y = {result.value(y):.6f} (expected: 1.0)")
print(f"Time: {result.wall_time:.3f}s")
******************************************************************************
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: 0.00000000 (expected: 0.0)
x = 1.000000 (expected: 1.0)
y = 1.000000 (expected: 1.0)
Time: 0.081s
3. Constrained Quadratic with Inequality (NLP)#
\[\min_{x,y} \; x^2 + y^2 \quad \text{s.t.} \quad x + y \geq 1, \quad x, y \in [-5, 5]\]
Known optimal: \(x^* = y^* = 0.5, \; f^* = 0.5\)
m = dm.Model("constrained_quadratic")
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, name="sum_lb")
result = m.solve()
print(f"Status: {result.status}")
print(f"Objective: {result.objective:.8f} (expected: 0.5)")
print(f"x = {result.value(x):.6f} (expected: 0.5)")
print(f"y = {result.value(y):.6f} (expected: 0.5)")
print(f"Time: {result.wall_time:.3f}s")
Status: optimal
Objective: 0.50000000 (expected: 0.5)
x = 0.500000 (expected: 0.5)
y = 0.500000 (expected: 0.5)
Time: 0.026s
4. Nonlinear Equality Constraint (NLP)#
\[\min_{x,y} \; x^2 + y^2 \quad \text{s.t.} \quad xy = 1, \quad x, y \in [0.1, 5]\]
Known optimal: \(x^* = y^* = 1, \; f^* = 2\)
m = dm.Model("nonlinear_equality")
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, name="product_eq")
result = m.solve()
print(f"Status: {result.status}")
print(f"Objective: {result.objective:.8f} (expected: 2.0)")
print(f"x = {result.value(x):.6f} (expected: 1.0)")
print(f"y = {result.value(y):.6f} (expected: 1.0)")
print(f"Time: {result.wall_time:.3f}s")
Status: optimal
Objective: 2.00000000 (expected: 2.0)
x = 1.000000 (expected: 1.0)
y = 1.000000 (expected: 1.0)
Time: 0.090s
5. Exponential NLP#
\[\min_{x,y} \; e^x + y^2 \quad \text{s.t.} \quad x + y \geq 1, \quad x, y \in [-2, 2]\]
From KKT conditions: \(e^x = 2y\) and \(x + y = 1\), giving \(f^* \approx 1.8395\).
m = dm.Model("exp_nlp")
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, name="sum_lb")
result = m.solve()
print(f"Status: {result.status}")
print(f"Objective: {result.objective:.8f} (expected: ~1.8395)")
print(f"x = {result.value(x):.6f}")
print(f"y = {result.value(y):.6f}")
print(f"Time: {result.wall_time:.3f}s")
Status: optimal
Objective: 1.83948429 (expected: ~1.8395)
x = 0.314923
y = 0.685077
Time: 0.085s
6. Simple MINLP (Branch & Bound)#
\[\min_{x_1, x_2, x_3} \; x_1^2 + x_2^2 + x_3 \quad \text{s.t.} \quad x_1 + x_2 \geq 1, \; x_1^2 + x_2 \leq 3, \; x_3 \in \{0, 1\}, \; x_1, x_2 \in [0, 5]\]
The binary variable \(x_3\) is minimized, so \(x_3^* = 0\). The continuous part reduces to a constrained quadratic.
Known optimal: \(x_1^* = x_2^* = 0.5, \; x_3^* = 0, \; f^* = 0.5\)
m = dm.Model("simple_minlp")
x1 = m.continuous("x1", lb=0, ub=5)
x2 = m.continuous("x2", lb=0, ub=5)
x3 = m.binary("x3")
m.minimize(x1**2 + x2**2 + x3)
m.subject_to(x1 + x2 >= 1, name="sum_lb")
m.subject_to(x1**2 + x2 <= 3, name="quad_ub")
result = m.solve()
print(f"Status: {result.status}")
print(f"Objective: {result.objective:.8f} (expected: 0.5)")
print(f"x1 = {result.value(x1):.6f} (expected: 0.5)")
print(f"x2 = {result.value(x2):.6f} (expected: 0.5)")
print(f"x3 = {result.value(x3):.6f} (expected: 0.0)")
print(f"Time: {result.wall_time:.3f}s | Nodes: {result.node_count}")
Status: optimal
Objective: 0.49999999 (expected: 0.5)
x1 = 0.500000 (expected: 0.5)
x2 = 0.500000 (expected: 0.5)
x3 = -0.000000 (expected: 0.0)
Time: 1.177s | Nodes: 1
7. Binary Knapsack-like MINLP#
\[\min_{x, y_1, y_2} \; (x - 3)^2 + 2 y_1 + 3 y_2 \quad \text{s.t.} \quad x \leq 5 y_1, \; x \leq 4 y_2, \; x \in [0, 5], \; y_1, y_2 \in \{0, 1\}\]
Known optimal: \(x^* = 3, \; y_1^* = 1, \; y_2^* = 1, \; f^* = 5.0\)
m = dm.Model("binary_knapsack_minlp")
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, name="link_y1")
m.subject_to(x <= 4 * y2, name="link_y2")
result = m.solve()
print(f"Status: {result.status}")
print(f"Objective: {result.objective:.8f} (expected: 5.0)")
print(f"x = {result.value(x):.6f} (expected: 3.0)")
print(f"y1 = {result.value(y1):.6f} (expected: 1.0)")
print(f"y2 = {result.value(y2):.6f} (expected: 1.0)")
print(f"Time: {result.wall_time:.3f}s | Nodes: {result.node_count}")
Status: optimal
Objective: 5.00000000 (expected: 5.0)
x = 3.000000 (expected: 3.0)
y1 = 1.000000 (expected: 1.0)
y2 = 1.000000 (expected: 1.0)
Time: 0.987s | Nodes: 5
8. Integer Quadratic#
\[\min_{n, x} \; (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\}\]
Known optimal: \(x^* = 2.7, \; n^* = 3, \; f^* = 0.0\)
m = dm.Model("integer_quadratic")
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, name="sum_ub")
result = m.solve()
print(f"Status: {result.status}")
print(f"Objective: {result.objective:.8f} (expected: 0.0)")
print(f"x = {result.value(x):.6f} (expected: 2.7)")
print(f"n = {result.value(n):.6f} (expected: 3.0)")
print(f"Time: {result.wall_time:.3f}s | Nodes: {result.node_count}")
Status: optimal
Objective: -0.00000000 (expected: 0.0)
x = 2.700000 (expected: 2.7)
n = 3.000000 (expected: 3.0)
Time: 0.227s | Nodes: 1
9. Logarithmic NLP#
\[\min_{x,y} \; -\log(x) - \log(y) \quad \text{s.t.} \quad x + y \leq 1, \; 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_nlp")
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, name="simplex")
expected = 2 * math.log(2)
result = m.solve()
print(f"Status: {result.status}")
print(f"Objective: {result.objective:.8f} (expected: {expected:.8f})")
print(f"x = {result.value(x):.6f} (expected: 0.5)")
print(f"y = {result.value(y):.6f} (expected: 0.5)")
print(f"Time: {result.wall_time:.3f}s")
Status: optimal
Objective: 1.38629434 (expected: 1.38629436)
x = 0.500000 (expected: 0.5)
y = 0.500000 (expected: 0.5)
Time: 0.082s
10. MINLPLib Instances (from .nl files)#
discopt can also load standard .nl format files from the MINLPLib [Bussieck et al., 2003] benchmark library. These are parsed by the Rust backend and solved through the same B&B engine.
Below we solve several MINLPLib instances and compare against published optimal values.
from pathlib import Path
# Known optima from MINLPLib (verified by BARON/ANTIGONE/SCIP)
minlplib_instances = {
"ex1221": {"opt": 7.66718007, "type": "MINLP"},
"ex1226": {"opt": -17.0, "type": "MINLP"},
"st_e13": {"opt": 2.0, "type": "NLP"},
"st_e15": {"opt": 7.66718007, "type": "NLP"},
"st_e27": {"opt": 2.0, "type": "NLP"},
"nvs03": {"opt": 16.0, "type": "MINLP"},
"nvs04": {"opt": 0.72, "type": "MINLP"},
"nvs07": {"opt": 4.0, "type": "MINLP"},
"nvs15": {"opt": 1.0, "type": "MINLP"},
"nvs16": {"opt": 0.70312500, "type": "MINLP"},
"prob03": {"opt": 10.0, "type": "MINLP"},
"gear": {"opt": 0.0, "type": "MINLP"},
}
nl_dir = Path("../python/tests/data/minlplib")
header = (
f"{'Instance':<12s} {'Type':<8s} {'Status':<10s} "
f"{'Objective':>14s} {'Expected':>14s} "
f"{'Error':>10s} {'Time':>8s}"
)
print(header)
print("-" * len(header))
n_correct = 0
n_total = 0
for name, info in minlplib_instances.items():
nl_path = nl_dir / f"{name}.nl"
if not nl_path.exists():
print(f"{name:<12s} SKIP .nl file not found")
continue
model = dm.from_nl(str(nl_path))
result = model.solve(time_limit=60, max_nodes=10000)
n_total += 1
expected = info["opt"]
if result.objective is not None:
error = abs(result.objective - expected)
tol = 1e-4 + 1e-3 * abs(expected)
correct = error <= tol
if correct:
n_correct += 1
error_str = f"{error:.2e}"
else:
error_str = "N/A"
correct = False
marker = " *" if not correct else ""
obj_str = f"{result.objective:.6f}" if result.objective is not None else "None"
print(
f"{name:<12s} {info['type']:<8s} "
f"{result.status:<10s} {obj_str:>14s} "
f"{expected:>14.6f} {error_str:>10s} "
f"{result.wall_time:>7.2f}s{marker}"
)
print(f"\nCorrect: {n_correct}/{n_total}")
Instance Type Status Objective Expected Error Time
----------------------------------------------------------------------------------
ex1221 SKIP .nl file not found
ex1226 SKIP .nl file not found
st_e13 SKIP .nl file not found
st_e15 SKIP .nl file not found
st_e27 SKIP .nl file not found
nvs03 SKIP .nl file not found
nvs04 SKIP .nl file not found
nvs07 SKIP .nl file not found
nvs15 SKIP .nl file not found
nvs16 SKIP .nl file not found
prob03 SKIP .nl file not found
gear SKIP .nl file not found
Correct: 0/0
11. Direct Comparison: discopt vs. cyipopt#
For continuous NLP problems, we can compare discopt’s solution against a raw cyipopt (Ipopt [Wächter and Biegler, 2006]) solve to verify correctness.
import cyipopt
# ===== Solve with discopt =====
m = dm.Model("comparison")
x = m.continuous("x", lb=-5, ub=5)
y = m.continuous("y", lb=-5, ub=5)
m.minimize((x - 2) ** 2 + (y + 1) ** 2 + x * y)
m.subject_to(x + y >= 0.5, name="sum_lb")
m.subject_to(x - y <= 3, name="diff_ub")
discopt_result = m.solve()
print("=== discopt ===")
print(f"Status: {discopt_result.status}")
print(f"Objective: {discopt_result.objective:.10f}")
print(f"x = {discopt_result.value(x):.8f}")
print(f"y = {discopt_result.value(y):.8f}")
print(f"Time: {discopt_result.wall_time:.4f}s")
print()
# ===== Solve with raw cyipopt =====
class NLPProblem:
def objective(self, x):
return (x[0] - 2) ** 2 + (x[1] + 1) ** 2 + x[0] * x[1]
def gradient(self, x):
return np.array(
[
2 * (x[0] - 2) + x[1],
2 * (x[1] + 1) + x[0],
]
)
def constraints(self, x):
return np.array(
[
x[0] + x[1], # >= 0.5
x[0] - x[1], # <= 3
]
)
def jacobian(self, x):
return np.array(
[
1.0,
1.0, # d(c1)/d(x0), d(c1)/d(x1)
1.0,
-1.0, # d(c2)/d(x0), d(c2)/d(x1)
]
)
nlp = cyipopt.Problem(
n=2,
m=2,
problem_obj=NLPProblem(),
lb=np.array([-5.0, -5.0]),
ub=np.array([5.0, 5.0]),
cl=np.array([0.5, -1e20]),
cu=np.array([1e20, 3.0]),
)
nlp.add_option("print_level", 0)
x_sol, info = nlp.solve(np.array([0.0, 0.0]))
print("=== raw cyipopt ===")
print(f"Status: {info['status_msg']}")
print(f"Objective: {info['obj_val']:.10f}")
print(f"x = {x_sol[0]:.8f}")
print(f"y = {x_sol[1]:.8f}")
print()
# ===== Compare =====
obj_diff = abs(discopt_result.objective - info["obj_val"])
x_diff = abs(float(discopt_result.value(x)) - x_sol[0])
y_diff = abs(float(discopt_result.value(y)) - x_sol[1])
print("=== Comparison ===")
print(f"Objective difference: {obj_diff:.2e}")
print(f"x difference: {x_diff:.2e}")
print(f"y difference: {y_diff:.2e}")
print(f"Match: {'YES' if obj_diff < 1e-6 else 'NO'}")
=== discopt ===
Status: optimal
Objective: -2.0833333333
x = 1.83333332
y = -1.16666668
Time: 0.1386s
=== raw cyipopt ===
Status: b'Algorithm terminated successfully at a locally optimal point, satisfying the convergence tolerances (can be specified by options).'
Objective: -2.0833333783
x = 1.83333336
y = -1.16666667
=== Comparison ===
Objective difference: 4.50e-08
x difference: 3.58e-08
y difference: 5.81e-09
Match: YES
12. Summary of the discopt API#
The discopt.modeling API provides a clean, Pythonic interface for MINLP formulation:
import discopt.modeling as dm
# Create model
m = dm.Model("my_problem")
# Variables
x = m.continuous("x", shape=(3,), lb=0, ub=10)
y = m.binary("y", shape=(2,))
n = m.integer("n", lb=0, ub=100)
# Objective
m.minimize(dm.sum(x) + dm.exp(x[0]) + 10*dm.sum(y))
# Constraints (named for clarity)
m.subject_to(x[0] + x[1] >= 1, name="lb")
m.subject_to(x[0] * x[1] <= 5 * y[0], name="bilinear") # bilinear
m.subject_to(dm.log(x[2] + 1) <= 2, name="log_ub") # nonlinear
# Solve
result = m.solve(time_limit=60)
print(result.status) # 'optimal', 'feasible', 'infeasible'
print(result.objective) # optimal objective value
print(result.value(x)) # solution value for variable x
print(result.value(y)) # solution value for variable y
print(result.wall_time) # solve time in seconds
print(result.node_count) # B&B nodes explored (0 for pure NLP)
Available nonlinear functions: dm.exp, dm.log, dm.sqrt, dm.sin, dm.cos, dm.tan, dm.abs
Load from files: dm.from_nl("problem.nl")