diff --git a/Cargo.lock b/Cargo.lock
index 193706e..9e80e43 100644
--- a/Cargo.lock
+++ b/Cargo.lock
@@ -432,7 +432,6 @@ dependencies = [
"oximo-core",
"oximo-expr",
"oximo-solver",
- "rayon",
"rustc-hash",
"thiserror",
]
diff --git a/README.md b/README.md
index 173bd12..67b623d 100644
--- a/README.md
+++ b/README.md
@@ -11,9 +11,7 @@
-oximo is a Rust algebraic modeling library for mathematical optimization. Build LP and MILP models with a clean builder API, then solve them with bundled or commercial solvers.
-
-> Support for nonlinear programming (NLP) and mixed-integer nonlinear programming (MINLP) is planned.
+oximo is a Rust algebraic modeling library for mathematical optimization. Build LP, MILP, QP/MIQP, NLP, and MINLP models with a clean builder API, then solve them with bundled or commercial solvers.
```rust,no_run
use oximo::prelude::*;
@@ -37,12 +35,12 @@ println!("y = {:?}", result.value_of(y)); // 4.0
## Features
-| Feature | What it adds | Default |
-|----------|---------------------------------------------------|---------|
-| `highs` | HiGHS LP/MILP solver (bundled, no install) | yes |
-| `io` | MPS and LP file writers | yes |
-| `gurobi` | Gurobi LP/MILP solver (requires licensed install) | no |
-| `gams` | GAMS solver bridge (requires GAMS on PATH) | no |
+| Feature | What it adds | Default |
+|----------|-----------------------------------------------------------------------|---------|
+| `highs` | HiGHS - LP/MILP solver (bundled, no install) | yes |
+| `io` | MPS and LP file writers | yes |
+| `gurobi` | Gurobi - LP/MILP/QP/MIQP/NLP/MINLP solver (requires licensed install) | no |
+| `gams` | GAMS bridge - LP/MILP/QP/MIQP/NLP/MINLP depending on solver | no |
```toml
[dependencies]
@@ -170,6 +168,23 @@ m.add_constraints_over("supply", &plants, |p: String| {
m.add_constraints_over("c", &set, |k: IndexKey| x[&k].le(1.0));
```
+### Nonlinear expressions
+
+`Pow`, `Sin`, `Cos`, `Exp`, `Log`, and bilinear products are first-class. The
+model's kind (`LP`/`MILP`/`QP`/`MIQP`/`NLP`/`MINLP`) is inferred from the
+expressions.
+
+```rust,ignore
+// Rosenbrock NLP
+m.minimize((1.0 - x).powi(2) + 100.0 * (y - x.powi(2)).powi(2));
+
+// Quadratic constraint
+m.constraint("disk", (x.powi(2) + y.powi(2)).le(1.0));
+
+// Transcendental utility (MINLP when any variable is integer/binary)
+m.maximize(sum_over(&items, |i: usize| u[i] * (1.0 + w[i] * x[i]).log()));
+```
+
## Solving
All backends implement the `Solver` trait:
diff --git a/crates/oximo-core/README.md b/crates/oximo-core/README.md
index b772d7c..a436aef 100644
--- a/crates/oximo-core/README.md
+++ b/crates/oximo-core/README.md
@@ -181,8 +181,6 @@ Inferred automatically from variables and expressions, cached and invalidated on
| `NLP` | All continuous, `Pow`/`Sin`/`Cos`/`Exp`/`Log` |
| `MINLP` | Any integer/binary + nonlinear |
-For now, we only support linear constraints, so `QP` and `NLP` are not possible. But the API is designed to allow nonlinear constraints in the future without breaking changes.
-
## License
MIT OR Apache-2.0
diff --git a/crates/oximo-core/src/model.rs b/crates/oximo-core/src/model.rs
index 7fabf2a..a59ab42 100644
--- a/crates/oximo-core/src/model.rs
+++ b/crates/oximo-core/src/model.rs
@@ -1,6 +1,6 @@
use std::cell::{Ref, RefCell};
-use oximo_expr::{Expr, ExprArena, VarId};
+use oximo_expr::{Expr, ExprArena, ExprClass, VarId, classify};
use rustc_hash::FxHashMap;
use smol_str::SmolStr;
@@ -269,38 +269,30 @@ impl Model {
}
let arena = self.arena.borrow();
let has_int = self.variables.borrow().iter().any(|v| v.domain.is_integer());
- let nonlinear = self.constraints.borrow().iter().any(|c| has_nonlinear(&arena, c.lhs))
- || self.objective.borrow().as_ref().is_some_and(|o| has_nonlinear(&arena, o.expr));
- let k = match (has_int, nonlinear) {
- (false, false) => ModelKind::LP,
- (true, false) => ModelKind::MILP,
- (false, true) => ModelKind::NLP,
- (true, true) => ModelKind::MINLP,
- };
- *self.cached_kind.borrow_mut() = Some(k);
- k
- }
-}
-fn has_nonlinear(arena: &ExprArena, id: oximo_expr::ExprId) -> bool {
- use oximo_expr::ExprNode as N;
- match arena.get(id) {
- N::Const(_) | N::Var(_) | N::Param(_) | N::Linear { .. } => false,
- N::Neg(inner) => has_nonlinear(arena, *inner),
- N::Add(children) => children.iter().any(|c| has_nonlinear(arena, *c)),
- N::Mul(children) => {
- let mut nonconst = 0;
- for c in children {
- if !matches!(arena.get(*c), N::Const(_)) {
- nonconst += 1;
- }
- if has_nonlinear(arena, *c) {
- return true;
- }
+ // Highest expression class across the objective and every constraint
+ // determines the model class
+ let mut class = ExprClass::Linear;
+ if let Some(o) = self.objective.borrow().as_ref() {
+ class = class.max(classify(&arena, o.expr));
+ }
+ for c in self.constraints.borrow().iter() {
+ if class == ExprClass::Nonlinear {
+ break;
}
- nonconst >= 2
+ class = class.max(classify(&arena, c.lhs));
}
- N::Pow(_, _) | N::Sin(_) | N::Cos(_) | N::Exp(_) | N::Log(_) => true,
+
+ let k = match (has_int, class) {
+ (false, ExprClass::Linear) => ModelKind::LP,
+ (true, ExprClass::Linear) => ModelKind::MILP,
+ (false, ExprClass::Quadratic) => ModelKind::QP,
+ (true, ExprClass::Quadratic) => ModelKind::MIQP,
+ (false, ExprClass::Nonlinear) => ModelKind::NLP,
+ (true, ExprClass::Nonlinear) => ModelKind::MINLP,
+ };
+ *self.cached_kind.borrow_mut() = Some(k);
+ k
}
}
diff --git a/crates/oximo-core/tests/model.rs b/crates/oximo-core/tests/model.rs
index dab88db..a8bd569 100644
--- a/crates/oximo-core/tests/model.rs
+++ b/crates/oximo-core/tests/model.rs
@@ -20,14 +20,53 @@ fn classifies_milp() {
assert_eq!(m.kind(), ModelKind::MILP);
}
+#[test]
+fn classifies_qp() {
+ let m = Model::new("qp");
+ let x = m.var("x").lb(0.0).build();
+ m.minimize(x.powi(2));
+ assert_eq!(m.kind(), ModelKind::QP);
+}
+
+#[test]
+fn classifies_miqp() {
+ let m = Model::new("miqp");
+ let x = m.var("x").lb(0.0).build();
+ let y = m.var("y").lb(0.0).ub(1.0).integer().build();
+ // Bilinear term keeps it quadratic, the integer var makes QP -> MIQP.
+ m.minimize(x * y);
+ assert_eq!(m.kind(), ModelKind::MIQP);
+}
+
+#[test]
+fn quadratic_constraint_classifies_qp() {
+ let m = Model::new("qp_con");
+ let x = m.var("x").lb(0.0).build();
+ // Linear objective but a quadratic constraint still makes the model a QP.
+ m.constraint("c", x.powi(2).le(4.0));
+ m.minimize(x);
+ assert_eq!(m.kind(), ModelKind::QP);
+}
+
#[test]
fn classifies_nlp() {
let m = Model::new("nlp");
let x = m.var("x").lb(0.0).build();
- m.minimize(x.powi(2));
+ // Degree-3, so it falls through to the nonlinear path.
+ m.minimize(x.powi(3));
assert_eq!(m.kind(), ModelKind::NLP);
}
+#[test]
+fn classifies_minlp_with_division() {
+ let m = Model::new("minlp_div");
+ let x = m.var("x").lb(1.0).build();
+ let y = m.var("y").lb(0.0).ub(1.0).integer().build();
+ // x / y is nonlinear (non-constant denominator)
+ m.minimize(x / y);
+ assert_eq!(m.kind(), ModelKind::MINLP);
+}
+
#[test]
fn variable_count_matches_register() {
let m = Model::new("vars");
diff --git a/crates/oximo-expr/README.md b/crates/oximo-expr/README.md
index 7421780..5dea422 100644
--- a/crates/oximo-expr/README.md
+++ b/crates/oximo-expr/README.md
@@ -28,6 +28,7 @@ Add(Children) // generic n-ary add
Mul(Children) // generic n-ary mul
Neg(ExprId)
Pow(ExprId, ExprId)
+Div(ExprId, ExprId) // numerator / denominator
Sin(ExprId) / Cos(ExprId) / Exp(ExprId) / Log(ExprId)
Linear { coeffs: Vec<(VarId, f64)>, constant: f64 } // LP fast-path
```
@@ -36,13 +37,14 @@ Linear { coeffs: Vec<(VarId, f64)>, constant: f64 } // LP fast-path
## Operator overloads
-`Expr` implements `Add`, `Sub`, `Mul`, `Neg` against other `Expr` values and against `f64`. All operations that stay linear produce a `Linear` node. For example:
+`Expr` implements `Add`, `Sub`, `Mul`, `Div`, `Neg` against other `Expr` values and against `f64`. All operations that stay linear produce a `Linear` node. For example:
```rust,ignore
// All of these produce a single Linear node, not an Add/Mul tree:
let e = 2.0 * x + 3.0 * y - 1.0;
let e = x + y;
let e = -x;
+let e = x / 2.0; // constant denominator: stays linear (x*0.5)
```
## Nonlinear methods on `Expr`
@@ -55,6 +57,7 @@ expr.sin()
expr.cos()
expr.exp()
expr.log()
+expr/expr
```
## Utilities
diff --git a/crates/oximo-expr/src/arena.rs b/crates/oximo-expr/src/arena.rs
index b56b74e..730999d 100644
--- a/crates/oximo-expr/src/arena.rs
+++ b/crates/oximo-expr/src/arena.rs
@@ -45,6 +45,7 @@ pub enum ExprNode {
Mul(Children),
Neg(ExprId),
Pow(ExprId, ExprId),
+ Div(ExprId, ExprId),
Sin(ExprId),
Cos(ExprId),
Exp(ExprId),
diff --git a/crates/oximo-expr/src/classify.rs b/crates/oximo-expr/src/classify.rs
new file mode 100644
index 0000000..ae72d8c
--- /dev/null
+++ b/crates/oximo-expr/src/classify.rs
@@ -0,0 +1,199 @@
+use crate::arena::{ExprArena, ExprId, ExprNode};
+
+/// Highest-degree polynomial class an expression belongs to, ignoring constant
+/// folding. Used by backends to pick between linear, quadratic, and general
+/// nonlinear translation paths.
+///
+/// Variants are ordered by increasing degree, so `max` of two classes yields the
+/// dominating one (e.g. a model with a quadratic objective and a nonlinear
+/// constraint is `Nonlinear`).
+#[derive(Copy, Clone, Debug, PartialEq, Eq, PartialOrd, Ord)]
+pub enum ExprClass {
+ Linear,
+ Quadratic,
+ Nonlinear,
+}
+
+/// Polynomial-degree bucket. `Higher` is a saturating sentinel for "anything
+/// above quadratic". Both polynomial degree > 2 and transcendentals collapse
+/// into it, since neither fits a QP solver's quadratic API.
+#[derive(Copy, Clone, Debug, PartialEq, Eq, PartialOrd, Ord)]
+enum Degree {
+ Zero,
+ One,
+ Two,
+ Higher,
+}
+
+impl Degree {
+ /// `+` on a sum: take the maximum, saturating at `Higher`.
+ fn add(self, other: Degree) -> Degree {
+ self.max(other)
+ }
+
+ /// `*` on a product: add ordinal degrees, saturating at `Higher`.
+ fn mul(self, other: Degree) -> Degree {
+ match (self, other) {
+ (Degree::Higher, _) | (_, Degree::Higher) => Degree::Higher,
+ (Degree::Zero, x) | (x, Degree::Zero) => x,
+ (Degree::One, Degree::One) => Degree::Two,
+ _ => Degree::Higher,
+ }
+ }
+
+ /// `^n` on a power: multiply by `n`, saturating at `Higher`.
+ fn pow(self, n: u32) -> Degree {
+ match (self, n) {
+ (_, 0) | (Degree::Zero, _) => Degree::Zero,
+ (d, 1) => d,
+ (Degree::One, 2) => Degree::Two,
+ _ => Degree::Higher,
+ }
+ }
+}
+
+fn degree(arena: &ExprArena, id: ExprId) -> Degree {
+ match arena.get(id) {
+ ExprNode::Const(_) => Degree::Zero,
+ ExprNode::Var(_) | ExprNode::Param(_) | ExprNode::Linear { .. } => Degree::One,
+ ExprNode::Neg(inner) => degree(arena, *inner),
+ ExprNode::Add(children) => {
+ let mut d = Degree::Zero;
+ for c in children {
+ d = d.add(degree(arena, *c));
+ if d == Degree::Higher {
+ return d;
+ }
+ }
+ d
+ }
+ ExprNode::Mul(children) => {
+ let mut d = Degree::Zero;
+ for c in children {
+ d = d.mul(degree(arena, *c));
+ if d == Degree::Higher {
+ return d;
+ }
+ }
+ d
+ }
+ ExprNode::Pow(base, exp) => {
+ let ExprNode::Const(e) = arena.get(*exp) else { return Degree::Higher };
+ if (*e - e.round()).abs() >= f64::EPSILON || *e < 0.0 {
+ return Degree::Higher;
+ }
+ // Bucket the exponent into the only values `Degree::pow` treats
+ // distinctly.
+ let n = match e.round() {
+ v if v < 0.5 => 0,
+ v if v < 1.5 => 1,
+ v if v < 2.5 => 2,
+ _ => 3,
+ };
+ degree(arena, *base).pow(n)
+ }
+ // Transcendentals are always > quadratic. Division is too: `div_into`
+ // folds the only degree-preserving case (constant denominator) before a
+ // `Div` node is created, so any other `Div` has a non-constant
+ // denominator.
+ ExprNode::Div(_, _)
+ | ExprNode::Sin(_)
+ | ExprNode::Cos(_)
+ | ExprNode::Exp(_)
+ | ExprNode::Log(_) => Degree::Higher,
+ }
+}
+
+/// Classify an expression as Linear, Quadratic (polynomial degree <= 2 with at
+/// least one degree-2 term), or Nonlinear (transcendentals, non-integer powers,
+/// or polynomial degree > 2).
+pub fn classify(arena: &ExprArena, id: ExprId) -> ExprClass {
+ match degree(arena, id) {
+ Degree::Zero | Degree::One => ExprClass::Linear,
+ Degree::Two => ExprClass::Quadratic,
+ Degree::Higher => ExprClass::Nonlinear,
+ }
+}
+
+#[cfg(test)]
+mod tests {
+ use super::*;
+ use crate::arena::{ExprArena, ExprNode, VarId};
+ use smallvec::smallvec;
+
+ fn var(arena: &mut ExprArena, i: u32) -> ExprId {
+ arena.push(ExprNode::Var(VarId(i)))
+ }
+
+ #[test]
+ fn linear_var_sum() {
+ let mut a = ExprArena::new();
+ let x = var(&mut a, 0);
+ let y = var(&mut a, 1);
+ let sum = a.push(ExprNode::Add(smallvec![x, y]));
+ assert_eq!(classify(&a, sum), ExprClass::Linear);
+ }
+
+ #[test]
+ fn quadratic_mul_two_vars() {
+ let mut a = ExprArena::new();
+ let x = var(&mut a, 0);
+ let y = var(&mut a, 1);
+ let xy = a.push(ExprNode::Mul(smallvec![x, y]));
+ assert_eq!(classify(&a, xy), ExprClass::Quadratic);
+ }
+
+ #[test]
+ fn quadratic_pow_two() {
+ let mut a = ExprArena::new();
+ let x = var(&mut a, 0);
+ let two = a.push(ExprNode::Const(2.0));
+ let sq = a.push(ExprNode::Pow(x, two));
+ assert_eq!(classify(&a, sq), ExprClass::Quadratic);
+ }
+
+ #[test]
+ fn nonlinear_pow_three() {
+ let mut a = ExprArena::new();
+ let x = var(&mut a, 0);
+ let three = a.push(ExprNode::Const(3.0));
+ let cube = a.push(ExprNode::Pow(x, three));
+ assert_eq!(classify(&a, cube), ExprClass::Nonlinear);
+ }
+
+ #[test]
+ fn nonlinear_div() {
+ let mut a = ExprArena::new();
+ let x = var(&mut a, 0);
+ let y = var(&mut a, 1);
+ let q = a.push(ExprNode::Div(x, y));
+ assert_eq!(classify(&a, q), ExprClass::Nonlinear);
+ }
+
+ #[test]
+ fn nonlinear_sin() {
+ let mut a = ExprArena::new();
+ let x = var(&mut a, 0);
+ let s = a.push(ExprNode::Sin(x));
+ assert_eq!(classify(&a, s), ExprClass::Nonlinear);
+ }
+
+ #[test]
+ fn nonlinear_triple_mul() {
+ let mut arena = ExprArena::new();
+ let x = var(&mut arena, 0);
+ let y = var(&mut arena, 1);
+ let z = var(&mut arena, 2);
+ let prod = arena.push(ExprNode::Mul(smallvec![x, y, z]));
+ assert_eq!(classify(&arena, prod), ExprClass::Nonlinear);
+ }
+
+ #[test]
+ fn linear_promoted_by_const_mul() {
+ let mut a = ExprArena::new();
+ let x = var(&mut a, 0);
+ let c = a.push(ExprNode::Const(3.0));
+ let m = a.push(ExprNode::Mul(smallvec![c, x]));
+ assert_eq!(classify(&a, m), ExprClass::Linear);
+ }
+}
diff --git a/crates/oximo-expr/src/eval.rs b/crates/oximo-expr/src/eval.rs
index 900608d..ec26676 100644
--- a/crates/oximo-expr/src/eval.rs
+++ b/crates/oximo-expr/src/eval.rs
@@ -43,6 +43,7 @@ pub fn evaluate(arena: &ExprArena, id: ExprId, ctx: &C) -> Resul
.try_fold(1.0, |acc, c| Ok::<_, EvalError>(acc * evaluate(arena, *c, ctx)?))?,
ExprNode::Neg(inner) => -evaluate(arena, *inner, ctx)?,
ExprNode::Pow(base, exp) => evaluate(arena, *base, ctx)?.powf(evaluate(arena, *exp, ctx)?),
+ ExprNode::Div(num, den) => evaluate(arena, *num, ctx)? / evaluate(arena, *den, ctx)?,
ExprNode::Sin(inner) => evaluate(arena, *inner, ctx)?.sin(),
ExprNode::Cos(inner) => evaluate(arena, *inner, ctx)?.cos(),
ExprNode::Exp(inner) => evaluate(arena, *inner, ctx)?.exp(),
diff --git a/crates/oximo-expr/src/lib.rs b/crates/oximo-expr/src/lib.rs
index 9731104..dd34f52 100644
--- a/crates/oximo-expr/src/lib.rs
+++ b/crates/oximo-expr/src/lib.rs
@@ -2,6 +2,7 @@
#![forbid(unsafe_code)]
mod arena;
+mod classify;
mod eval;
mod handle;
mod linear;
@@ -10,6 +11,7 @@ mod simplify;
mod visit;
pub use arena::{ExprArena, ExprId, ExprNode, ParamId, VarId};
+pub use classify::{ExprClass, classify};
pub use eval::{EvalContext, EvalError, evaluate};
pub use handle::Expr;
pub use linear::{LinearTerms, extract_linear};
diff --git a/crates/oximo-expr/src/linear.rs b/crates/oximo-expr/src/linear.rs
index 5c48ca8..853859b 100644
--- a/crates/oximo-expr/src/linear.rs
+++ b/crates/oximo-expr/src/linear.rs
@@ -118,6 +118,25 @@ pub(crate) fn mul_into(arena: &mut ExprArena, lhs: ExprId, rhs: ExprId) -> ExprI
arena.push(ExprNode::Mul(smallvec![lhs, rhs]))
}
+/// Build `num / den`. If `den` is a nonzero constant `c`, fold to `num * (1/c)`
+/// so a constant-denominator division stays on the linear fast-path. Otherwise
+/// produce a `Div` node (always nonlinear, even when the numerator is linear).
+pub(crate) fn div_into(arena: &mut ExprArena, num: ExprId, den: ExprId) -> ExprId {
+ if let ExprNode::Const(c) = *arena.get(den) {
+ if c != 0.0 {
+ if let Some(mut t) = as_linear(arena, num) {
+ let inv = 1.0 / c;
+ t.coeffs.iter_mut().for_each(|(_, co)| *co *= inv);
+ t.constant *= inv;
+ return push_linear(arena, t);
+ }
+ let inv = arena.push(ExprNode::Const(1.0 / c));
+ return mul_into(arena, num, inv);
+ }
+ }
+ arena.push(ExprNode::Div(num, den))
+}
+
/// Build `-rhs`, preserving linearity.
pub(crate) fn neg_into(arena: &mut ExprArena, rhs: ExprId) -> ExprId {
if let Some(mut t) = as_linear(arena, rhs) {
diff --git a/crates/oximo-expr/src/ops.rs b/crates/oximo-expr/src/ops.rs
index b5d85ba..880d931 100644
--- a/crates/oximo-expr/src/ops.rs
+++ b/crates/oximo-expr/src/ops.rs
@@ -1,7 +1,7 @@
-use std::ops::{Add, Mul, Neg, Sub};
+use std::ops::{Add, Div, Mul, Neg, Sub};
use crate::handle::Expr;
-use crate::linear::{add_into, mul_into, neg_into, sub_into};
+use crate::linear::{add_into, div_into, mul_into, neg_into, sub_into};
// -----------------------------------------------------------------------------
// Expr Expr
@@ -31,6 +31,14 @@ impl<'a> Mul for Expr<'a> {
}
}
+impl<'a> Div for Expr<'a> {
+ type Output = Self;
+ fn div(self, rhs: Self) -> Self {
+ let id = div_into(&mut self.arena.borrow_mut(), self.id, rhs.id);
+ Self::new(id, self.arena)
+ }
+}
+
impl<'a> Neg for Expr<'a> {
type Output = Self;
fn neg(self) -> Self {
@@ -45,14 +53,13 @@ impl<'a> Neg for Expr<'a> {
// -----------------------------------------------------------------------------
macro_rules! impl_scalar_ops {
- ($scalar:ty) => {
+ ($scalar:ty, $to_f64:expr) => {
impl<'a> Add<$scalar> for Expr<'a> {
type Output = Self;
fn add(self, rhs: $scalar) -> Self {
- #[allow(clippy::cast_lossless)]
let id = {
let mut a = self.arena.borrow_mut();
- let rhs_id = a.constant(rhs as f64);
+ let rhs_id = a.constant($to_f64(rhs));
add_into(&mut a, self.id, rhs_id)
};
Self::new(id, self.arena)
@@ -69,10 +76,9 @@ macro_rules! impl_scalar_ops {
impl<'a> Sub<$scalar> for Expr<'a> {
type Output = Self;
fn sub(self, rhs: $scalar) -> Self {
- #[allow(clippy::cast_lossless)]
let id = {
let mut a = self.arena.borrow_mut();
- let rhs_id = a.constant(rhs as f64);
+ let rhs_id = a.constant($to_f64(rhs));
sub_into(&mut a, self.id, rhs_id)
};
Self::new(id, self.arena)
@@ -82,10 +88,9 @@ macro_rules! impl_scalar_ops {
impl<'a> Sub> for $scalar {
type Output = Expr<'a>;
fn sub(self, rhs: Expr<'a>) -> Expr<'a> {
- #[allow(clippy::cast_lossless)]
let id = {
let mut a = rhs.arena.borrow_mut();
- let lhs_id = a.constant(self as f64);
+ let lhs_id = a.constant($to_f64(self));
sub_into(&mut a, lhs_id, rhs.id)
};
Expr::new(id, rhs.arena)
@@ -95,10 +100,9 @@ macro_rules! impl_scalar_ops {
impl<'a> Mul<$scalar> for Expr<'a> {
type Output = Self;
fn mul(self, rhs: $scalar) -> Self {
- #[allow(clippy::cast_lossless)]
let id = {
let mut a = self.arena.borrow_mut();
- let rhs_id = a.constant(rhs as f64);
+ let rhs_id = a.constant($to_f64(rhs));
mul_into(&mut a, self.id, rhs_id)
};
Self::new(id, self.arena)
@@ -111,11 +115,35 @@ macro_rules! impl_scalar_ops {
rhs * self
}
}
+
+ impl<'a> Div<$scalar> for Expr<'a> {
+ type Output = Self;
+ fn div(self, rhs: $scalar) -> Self {
+ let id = {
+ let mut a = self.arena.borrow_mut();
+ let rhs_id = a.constant($to_f64(rhs));
+ div_into(&mut a, self.id, rhs_id)
+ };
+ Self::new(id, self.arena)
+ }
+ }
+
+ impl<'a> Div> for $scalar {
+ type Output = Expr<'a>;
+ fn div(self, rhs: Expr<'a>) -> Expr<'a> {
+ let id = {
+ let mut a = rhs.arena.borrow_mut();
+ let lhs_id = a.constant($to_f64(self));
+ div_into(&mut a, lhs_id, rhs.id)
+ };
+ Expr::new(id, rhs.arena)
+ }
+ }
};
}
-impl_scalar_ops!(f64);
-impl_scalar_ops!(i32);
+impl_scalar_ops!(f64, core::convert::identity);
+impl_scalar_ops!(i32, f64::from);
// -----------------------------------------------------------------------------
// std::iter::Sum: the first element of the iterator carries the arena handle,
diff --git a/crates/oximo-expr/src/simplify.rs b/crates/oximo-expr/src/simplify.rs
index d39cb34..19742bf 100644
--- a/crates/oximo-expr/src/simplify.rs
+++ b/crates/oximo-expr/src/simplify.rs
@@ -18,6 +18,10 @@ pub fn simplify(arena: &mut ExprArena, id: ExprId) -> ExprId {
(ExprNode::Const(b), ExprNode::Const(e)) => Some(ExprNode::Const(b.powf(*e))),
_ => None,
},
+ ExprNode::Div(num, den) => match (arena.get(num), arena.get(den)) {
+ (ExprNode::Const(n), ExprNode::Const(d)) => Some(ExprNode::Const(n / d)),
+ _ => None,
+ },
ExprNode::Sin(inner)
| ExprNode::Cos(inner)
| ExprNode::Exp(inner)
diff --git a/crates/oximo-expr/src/visit.rs b/crates/oximo-expr/src/visit.rs
index 690e294..aae04d2 100644
--- a/crates/oximo-expr/src/visit.rs
+++ b/crates/oximo-expr/src/visit.rs
@@ -31,6 +31,12 @@ pub fn walk(arena: &ExprArena, id: ExprId, visitor: &mut V) {
walk(arena, base, visitor);
walk(arena, exp, visitor);
}
+ ExprNode::Div(num, den) => {
+ let num = *num;
+ let den = *den;
+ walk(arena, num, visitor);
+ walk(arena, den, visitor);
+ }
ExprNode::Const(_) | ExprNode::Var(_) | ExprNode::Param(_) | ExprNode::Linear { .. } => {}
}
}
diff --git a/crates/oximo-expr/tests/arithmetic.rs b/crates/oximo-expr/tests/arithmetic.rs
index 5278a46..e4a7409 100644
--- a/crates/oximo-expr/tests/arithmetic.rs
+++ b/crates/oximo-expr/tests/arithmetic.rs
@@ -2,7 +2,9 @@
use std::cell::RefCell;
-use oximo_expr::{Expr, ExprArena, ExprNode, VarId, dot, evaluate, extract_linear};
+use oximo_expr::{
+ Expr, ExprArena, ExprClass, ExprNode, VarId, classify, dot, evaluate, extract_linear,
+};
fn make_var(arena: &RefCell, idx: u32) -> Expr<'_> {
Expr::from_var(arena, VarId(idx))
@@ -105,6 +107,66 @@ fn dot_panics_on_empty() {
let _ = dot(&xs, &coeffs);
}
+#[test]
+fn div_by_constant_stays_linear() {
+ let arena = RefCell::new(ExprArena::new());
+ let x = make_var(&arena, 0);
+ let combo = x / 2.0;
+
+ let snapshot = arena.borrow().get(combo.id).clone();
+ match snapshot {
+ ExprNode::Linear { coeffs, constant } => {
+ assert_eq!(constant, 0.0);
+ assert_eq!(coeffs, vec![(VarId(0), 0.5)]);
+ }
+ n => panic!("expected Linear node, got {n:?}"),
+ }
+ assert_eq!(classify(&arena.borrow(), combo.id), ExprClass::Linear);
+}
+
+#[test]
+fn div_two_vars_is_nonlinear() {
+ let arena = RefCell::new(ExprArena::new());
+ let a = make_var(&arena, 0);
+ let b = make_var(&arena, 1);
+ let q = a / b;
+
+ assert!(matches!(arena.borrow().get(q.id), ExprNode::Div(_, _)));
+ assert_eq!(classify(&arena.borrow(), q.id), ExprClass::Nonlinear);
+ assert!(extract_linear(&arena.borrow(), q.id).is_none());
+}
+
+#[test]
+fn scalar_over_var_is_nonlinear() {
+ let arena = RefCell::new(ExprArena::new());
+ let x = make_var(&arena, 0);
+ let recip = 1.0 / x;
+ assert!(matches!(arena.borrow().get(recip.id), ExprNode::Div(_, _)));
+ assert_eq!(classify(&arena.borrow(), recip.id), ExprClass::Nonlinear);
+}
+
+#[test]
+fn evaluate_division() {
+ let arena = RefCell::new(ExprArena::new());
+ let a = make_var(&arena, 0);
+ let b = make_var(&arena, 1);
+ let q = a / b;
+ let values: &[f64] = &[12.0, 4.0];
+ let arena_ref = arena.borrow();
+ assert_eq!(evaluate(&arena_ref, q.id, &values).unwrap(), 3.0);
+}
+
+#[test]
+fn evaluate_division_by_zero_is_infinite() {
+ let arena = RefCell::new(ExprArena::new());
+ let a = make_var(&arena, 0);
+ let b = make_var(&arena, 1);
+ let q = a / b;
+ let values: &[f64] = &[1.0, 0.0];
+ let arena_ref = arena.borrow();
+ assert!(evaluate(&arena_ref, q.id, &values).unwrap().is_infinite());
+}
+
#[test]
fn large_sum_extracts_correctly() {
let arena = RefCell::new(ExprArena::new());
diff --git a/crates/oximo-gams/README.md b/crates/oximo-gams/README.md
index 8b25798..3c3dc49 100644
--- a/crates/oximo-gams/README.md
+++ b/crates/oximo-gams/README.md
@@ -1,8 +1,8 @@
# oximo-gams
-GAMS LP/MILP backend for [oximo](https://github.com/germanheim/oximo).
+GAMS backend for [oximo](https://github.com/germanheim/oximo).
-Writes an oximo`Model`] to a temporary `.gms` file, invokes the GAMS executable via `std::process::Command`, and parses the solution from a PUT-generated text file. Supports `LP` and `MILP` model kinds. **QP/NLP/MINLP return `SolverError::UnsupportedKind` for now**.
+Writes an oximo `Model` to a temporary `.gms` file, invokes the GAMS executable via `std::process::Command`, and parses the solution from a PUT-generated text file. Supports `LP`, `MILP`, `QP`, `MIQP`, `NLP`, and `MINLP` model kinds.
The sub-solver is determined by the GAMS installation (default) or set explicitly via `GamsOptions::solver`. Any solver available in your GAMS distribution can be targeted, see [Sub-solver selection](#sub-solver-selection) below.
@@ -100,7 +100,8 @@ let opts = GamsOptions::default()
## Sub-solver selection
Pass a `GamsSolverConfig` to `.solver(...)` to select a GAMS sub-solver. This emits
-`option {LP|MIP} = ;` in the generated `.gms` file.
+`option {LP|MIP|NLP|MINLP|QCP|MIQCP} = ;` in the generated `.gms` file, scoped to
+the solve type resolved from `Model::kind()` (`QP` -> `QCP`, `MIQP` -> `MIQCP`).
```rust
use oximo_gams::{GamsOptions, GamsSolver, GamsSolverConfig};
diff --git a/crates/oximo-gams/src/lib.rs b/crates/oximo-gams/src/lib.rs
index d39737e..7026ccb 100644
--- a/crates/oximo-gams/src/lib.rs
+++ b/crates/oximo-gams/src/lib.rs
@@ -50,7 +50,15 @@ impl Solver for Gams {
}
fn supports(&self, kind: ModelKind) -> bool {
- matches!(kind, ModelKind::LP | ModelKind::MILP)
+ matches!(
+ kind,
+ ModelKind::LP
+ | ModelKind::MILP
+ | ModelKind::QP
+ | ModelKind::MIQP
+ | ModelKind::NLP
+ | ModelKind::MINLP
+ )
}
fn solve(&mut self, model: &Model, opts: &GamsOptions) -> Result {
diff --git a/crates/oximo-gams/src/options.rs b/crates/oximo-gams/src/options.rs
index 5cc093e..fe8a25f 100644
--- a/crates/oximo-gams/src/options.rs
+++ b/crates/oximo-gams/src/options.rs
@@ -25,15 +25,15 @@ pub struct GamsOptions {
/// Reference:
#[derive(Clone, Debug, Eq, PartialEq)]
pub enum GamsSolver {
- /// ALPHAECP: MINLP
+ /// ALPHAECP: MINLP, MIQCP
AlphaEcp,
- /// ANTIGONE: LP, MIP, NLP, MINLP, QCP, MIQCP, Global
+ /// ANTIGONE: NLP, CNS, DNLP, MINLP, QCP, MIQCP, Global
Antigone,
- /// BARON: LP, MIP, NLP, MCP, MPEC, CNS, MINLP, QCP, MIQCP, Global
+ /// BARON: LP, MIP, NLP, CNS, DNLP, MINLP, QCP, MIQCP, Global
Baron,
/// CBC: LP, MIP
Cbc,
- /// CONOPT: NLP, DNLP, CNS, MPEC
+ /// CONOPT: LP, NLP, CNS, DNLP, QCP
Conopt,
/// COPT: LP, MIP, QCP, MIQCP
Copt,
@@ -41,55 +41,55 @@ pub enum GamsSolver {
Cplex,
/// DECIS: LP, Stochastic
Decis,
- /// DICOPT: MINLP
+ /// DICOPT: MINLP, MIQCP
Dicopt,
/// GLPK: LP, MIP (not in GAMS docs but recognized)
Glpk,
- /// GUROBI: LP, MIP, NLP, MINLP, QCP, MIQCP, Global
+ /// GUROBI: LP, MIP, NLP, DNLP, MINLP, QCP, MIQCP, Global
Gurobi,
- /// GUSS: LP, MIP, NLP, MCP, MPEC, CNS, DNLP, MINLP, QCP, MIQCP
+ /// GUSS: LP, MIP, NLP, MCP, CNS, DNLP, MINLP, QCP, MIQCP
Guss,
/// HiGHS: LP, MIP
Highs,
- /// IPOPT: NLP, DNLP, CNS, MPEC
+ /// IPOPT: LP, NLP, CNS, DNLP, QCP
Ipopt,
- /// JAMS: MPEC
+ /// JAMS: EMP
Jams,
/// KESTREL: all model types (remote solver submission)
Kestrel,
- /// KNITRO: LP, MIP, NLP, MCP, MPEC, CNS, DNLP, MINLP, QCP, MIQCP
+ /// KNITRO: LP, NLP, MCP, MPEC, CNS, DNLP, MINLP, QCP, MIQCP
Knitro,
- /// LINDO: LP, MIP, NLP, MCP, MPEC, CNS, DNLP, MINLP, QCP, Stochastic, Global
+ /// LINDO: LP, MIP, NLP, DNLP, MINLP, QCP, MIQCP, Stochastic, Global
Lindo,
- /// LINDOGLOBAL: LP, MIP, NLP, MINLP, QCP, MIQCP, Global
+ /// LINDOGLOBAL: LP, MIP, NLP, DNLP, MINLP, QCP, MIQCP, Global
LindoGlobal,
/// MILES: MCP
Miles,
- /// MINOS: NLP, DNLP, CNS, MPEC
+ /// MINOS: LP, NLP, CNS, DNLP, QCP
Minos,
- /// MOSEK: LP, MIP, NLP, QCP, MIQCP
+ /// MOSEK: LP, MIP, NLP, DNLP, MINLP, QCP, MIQCP
Mosek,
- /// NLPEC: NLP, MPEC
+ /// NLPEC: MCP, MPEC
Nlpec,
- /// ODHCPLEX: LP, MIP
+ /// ODHCPLEX: MIP, MIQCP
OdhCplex,
- /// PATH: MCP, MPEC
+ /// PATH: MCP, CNS
Path,
- /// QUADMINOS: NLP
+ /// QUADMINOS: LP
QuadMinos,
- /// RESHOP: NLP
+ /// RESHOP: EMP
Reshop,
- /// SBB: NLP, MINLP
+ /// SBB: MINLP, MIQCP
Sbb,
- /// SCIP: LP, MIP, NLP, MINLP, QCP, MIQCP, Global
+ /// SCIP: MIP, NLP, CNS, DNLP, MINLP, QCP, MIQCP, Global
Scip,
- /// SHOT: MINLP
+ /// SHOT: MINLP, MIQCP
Shot,
- /// SNOPT: NLP, DNLP, CNS, MPEC
+ /// SNOPT: LP, NLP, CNS, DNLP, QCP
Snopt,
/// SOPLEX: LP
Soplex,
- /// XPRESS: LP, MIP, NLP, MINLP, QCP, MIQCP, Global
+ /// XPRESS: LP, MIP, NLP, CNS, DNLP, MINLP, QCP, MIQCP, Global
Xpress,
/// Any other GAMS-recognized solver name, emitted verbatim.
Custom(String),
@@ -170,7 +170,8 @@ impl HasUniversal for GamsOptions {
/// Emit GAMS option statements into `gms` before the `Solve` statement.
///
-/// `solve_type` is `"LP"` or `"MIP"`, used to scope the `solver` option.
+/// `solve_type` is the GAMS model type (`"LP"` / `"MIP"` / `"NLP"` / `"MINLP"`
+/// / `"QCP"` / `"MIQCP"`), used to scope the `solver` option.
pub fn write_options(gms: &mut String, o: &GamsOptions, solve_type: &str) {
if let Some(d) = o.universal.time_limit {
writeln!(gms, "option ResLim = {};", d.as_secs_f64()).unwrap();
diff --git a/crates/oximo-gams/src/solver_options.rs b/crates/oximo-gams/src/solver_options.rs
index e924925..5a9239e 100644
--- a/crates/oximo-gams/src/solver_options.rs
+++ b/crates/oximo-gams/src/solver_options.rs
@@ -3,7 +3,10 @@
use std::fmt::Write as FmtWrite;
+use oximo_core::ModelKind;
+
use crate::options::GamsSolver;
+use crate::translate::gams_solve_type;
// - Config enum
@@ -73,6 +76,53 @@ impl GamsSolverConfig {
Self::Named(_) => false,
}
}
+
+ /// Whether this solver can handle `kind` under oximo's GAMS translation,
+ /// which emits `QP` as a `QCP` solve and `MIQP` as a `MIQCP` solve.
+ ///
+ /// [`GamsSolver::Custom`] and any unrecognized name return `true`: their
+ /// capabilities are unknown, so they are left for GAMS to accept or reject.
+ #[must_use]
+ pub fn supports(&self, kind: ModelKind) -> bool {
+ solver_supports_type(self.gams_name(), gams_solve_type(kind))
+ }
+}
+
+/// GAMS solve types a named solver supports, restricted to the six oximo can
+/// emit: `LP` / `MIP` / `NLP` / `MINLP` / `QCP` / `MIQCP`. `None` means the
+/// name is unrecognized and cannot be validated.
+///
+/// Transcribed from the GAMS solver/model-type matrix (other model types —
+/// `MCP`, `MPEC`, `CNS`, `DNLP`, `EMP`, stochastic — are omitted because oximo
+/// never emits them):
+/// - "GAMS Solver Manuals," GAMS Development Corporation.
+/// (accessed May 14, 2026).
+fn supported_solve_types(gams_name: &str) -> Option<&'static [&'static str]> {
+ Some(match gams_name {
+ "ALPHAECP" | "DICOPT" | "SBB" | "SHOT" => &["MINLP", "MIQCP"],
+ "CONOPT" | "CONOPT3" | "CONOPT4" | "IPOPT" | "MINOS" | "SNOPT" => &["LP", "NLP", "QCP"],
+ "DECIS" | "SOPLEX" | "QUADMINOS" => &["LP"],
+ "CBC" | "GLPK" | "HIGHS" => &["LP", "MIP"],
+ "ODHCPLEX" => &["MIP", "MIQCP"],
+ "COPT" | "CPLEX" => &["LP", "MIP", "QCP", "MIQCP"],
+ "ANTIGONE" => &["NLP", "MINLP", "QCP", "MIQCP"],
+ "KNITRO" => &["LP", "NLP", "MINLP", "QCP", "MIQCP"],
+ "SCIP" => &["MIP", "NLP", "MINLP", "QCP", "MIQCP"],
+ "BARON" | "GUROBI" | "GUSS" | "KESTREL" | "LINDO" | "LINDOGLOBAL" | "MOSEK" | "XPRESS" => {
+ &["LP", "MIP", "NLP", "MINLP", "QCP", "MIQCP"]
+ }
+ // JAMS (EMP), MILES (MCP), NLPEC (MCP/MPEC), PATH (MCP/MPEC/CNS),
+ // RESHOP (EMP) support none of the model types oximo emits.
+ "JAMS" | "MILES" | "NLPEC" | "PATH" | "RESHOP" => &[],
+ _ => return None,
+ })
+}
+
+/// Whether the GAMS solver named `gams_name` supports `gams_type`
+/// (`"LP"` / `"MIP"` / `"NLP"` / `"MINLP"` / `"QCP"` / `"MIQCP"`). Unrecognized
+/// names return `true`.
+pub(crate) fn solver_supports_type(gams_name: &str, gams_type: &str) -> bool {
+ supported_solve_types(gams_name).is_none_or(|types| types.contains(&gams_type))
}
impl From for GamsSolverConfig {
@@ -1017,4 +1067,51 @@ mod tests {
assert!(matches!(cfg, GamsSolverConfig::Named(GamsSolver::Gurobi)));
assert_eq!(cfg.gams_name(), "GUROBI");
}
+
+ #[test]
+ fn supports_matches_solver_capabilities() {
+ // CPLEX: LP, MIP, QCP, MIQCP. QP/MIQP route through QCP/MIQCP, so the
+ // quadratic kinds pass while the general (MI)NLP kinds fail.
+ let cplex = GamsSolverConfig::Cplex(GamsCplexOptions::default());
+ assert!(cplex.supports(ModelKind::LP));
+ assert!(cplex.supports(ModelKind::MILP));
+ assert!(cplex.supports(ModelKind::QP), "QP routes through QCP");
+ assert!(cplex.supports(ModelKind::MIQP), "MIQP routes through MIQCP");
+ assert!(!cplex.supports(ModelKind::NLP));
+ assert!(!cplex.supports(ModelKind::MINLP));
+
+ // IPOPT: LP, NLP, QCP. LP, NLP, and QP pass, the integer kinds fail.
+ let ipopt = GamsSolverConfig::Ipopt(GamsIpoptOptions::default());
+ assert!(ipopt.supports(ModelKind::LP));
+ assert!(ipopt.supports(ModelKind::NLP));
+ assert!(ipopt.supports(ModelKind::QP), "QP routes through QCP");
+ assert!(!ipopt.supports(ModelKind::MIQP), "MIQP routes through MIQCP");
+ assert!(!ipopt.supports(ModelKind::MINLP));
+
+ // BARON handles all six oximo solve types.
+ let baron = GamsSolverConfig::Named(GamsSolver::Baron);
+ for k in [
+ ModelKind::LP,
+ ModelKind::MILP,
+ ModelKind::QP,
+ ModelKind::MIQP,
+ ModelKind::NLP,
+ ModelKind::MINLP,
+ ] {
+ assert!(baron.supports(k), "BARON should support {k:?}");
+ }
+
+ // HiGHS: LP/MIP only, no quadratic or nonlinear support THROUGH GAMS.
+ let highs = GamsSolverConfig::Highs(GamsHighsOptions::default());
+ assert!(highs.supports(ModelKind::LP));
+ assert!(!highs.supports(ModelKind::QP));
+ assert!(!highs.supports(ModelKind::NLP));
+ }
+
+ #[test]
+ fn supports_is_permissive_for_unknown_names() {
+ let custom = GamsSolverConfig::Named(GamsSolver::Custom("MYSOLVER".into()));
+ assert!(custom.supports(ModelKind::MINLP));
+ assert!(custom.supports(ModelKind::LP));
+ }
}
diff --git a/crates/oximo-gams/src/translate.rs b/crates/oximo-gams/src/translate.rs
index 85bc7f6..a874d97 100644
--- a/crates/oximo-gams/src/translate.rs
+++ b/crates/oximo-gams/src/translate.rs
@@ -6,10 +6,12 @@ use std::{fs, io};
static SOLVE_ID: AtomicU64 = AtomicU64::new(0);
-use oximo_core::{ConstraintId, Domain, Model, ModelKind, ObjectiveSense, Sense, VarId};
-use oximo_expr::{ExprArena, LinearTerms, extract_linear};
+use oximo_core::{
+ Constraint, ConstraintId, Domain, Model, ModelKind, Objective, ObjectiveSense, Sense, VarId,
+ Variable,
+};
+use oximo_expr::{ExprArena, ExprId, ExprNode, LinearTerms, extract_linear};
use oximo_solver::{SolverError, SolverResult, SolverStatus};
-use rayon::prelude::*;
use rustc_hash::FxHashMap;
use crate::GamsOptions;
@@ -36,167 +38,28 @@ pub fn solve(
exec: Option<&str>,
) -> Result {
let kind = model.kind();
- if !matches!(kind, ModelKind::LP | ModelKind::MILP) {
- return Err(SolverError::UnsupportedKind(kind));
- }
-
+ validate_solver(opts, kind)?;
let arena = model.arena();
let vars = model.variables();
let constraints = model.constraints();
let objective = model.try_objective().map_err(SolverError::Core)?;
- let obj_terms = extract_linear(&arena, objective.expr).ok_or(SolverError::Nonlinear)?;
-
- let arena_ref: &ExprArena = &arena;
- let con_terms: Vec = constraints
- .par_iter()
- .map(|c| extract_linear(arena_ref, c.lhs).ok_or(SolverError::Nonlinear))
- .collect::, _>>()?;
-
- let solve_type = if matches!(kind, ModelKind::MILP) { "MIP" } else { "LP" };
let sense_kw = match objective.sense {
ObjectiveSense::Minimize => "minimizing",
ObjectiveSense::Maximize => "maximizing",
};
- // Pre-compute solver opt file content so we know whether to inject optfile=1.
- let solver_opt: Option<(String, String)> = opts.solver.as_ref().and_then(|cfg| {
- let mut buf = String::new();
- if cfg.write_opt_file(&mut buf) {
- let fname = format!("{}.opt", cfg.gams_name().to_ascii_lowercase());
- Some((fname, buf))
- } else {
- None
- }
- });
-
- // - Build the .gms file
let mut gms = String::with_capacity(4096);
-
- writeln!(gms, "$title oximo_model").unwrap();
- writeln!(gms, "$offSymList").unwrap();
- writeln!(gms, "$offSymXRef").unwrap();
- writeln!(gms, "option solprint = off;").unwrap();
- writeln!(gms, "option limrow = 0;").unwrap();
- writeln!(gms, "option limcol = 0;").unwrap();
- writeln!(gms).unwrap();
-
- // Variable declarations, split by domain
- let (mut cont_vars, mut bin_vars, mut int_vars) = (Vec::new(), Vec::new(), Vec::new());
- for v in vars.iter() {
- match v.domain {
- Domain::Binary => bin_vars.push(v),
- Domain::Integer | Domain::SemiInteger { .. } => int_vars.push(v),
- _ => cont_vars.push(v),
- }
- }
-
- // Continuous + objective variable
- write!(gms, "Variables\n v_obj").unwrap();
- for v in &cont_vars {
- write!(gms, ", v{}", v.id.index()).unwrap();
- }
- writeln!(gms, ";").unwrap();
-
- if !bin_vars.is_empty() {
- write!(gms, "Binary Variables").unwrap();
- for (k, v) in bin_vars.iter().enumerate() {
- if k == 0 {
- write!(gms, "\n v{}", v.id.index()).unwrap();
- } else {
- write!(gms, ", v{}", v.id.index()).unwrap();
- }
- }
- writeln!(gms, ";").unwrap();
- }
-
- if !int_vars.is_empty() {
- write!(gms, "Integer Variables").unwrap();
- for (k, v) in int_vars.iter().enumerate() {
- if k == 0 {
- write!(gms, "\n v{}", v.id.index()).unwrap();
- } else {
- write!(gms, ", v{}", v.id.index()).unwrap();
- }
- }
- writeln!(gms, ";").unwrap();
- }
- writeln!(gms).unwrap();
-
- // Bounds
- for v in vars.iter() {
- let i = v.id.index();
- if matches!(v.domain, Domain::Binary) {
- // Default binary bounds are [0, 1], only emit when overridden.
- if (v.lb - v.ub).abs() < f64::EPSILON {
- writeln!(gms, "v{i}.fx = {};", fmt(v.lb)).unwrap();
- } else {
- if v.lb.abs() > f64::EPSILON {
- writeln!(gms, "v{i}.lo = {};", fmt(v.lb)).unwrap();
- }
- if (v.ub - 1.0).abs() > f64::EPSILON {
- writeln!(gms, "v{i}.up = {};", fmt(v.ub)).unwrap();
- }
- }
- continue;
- }
- // Lower bound
- if v.lb == f64::NEG_INFINITY {
- writeln!(gms, "v{i}.lo = -Inf;").unwrap();
- } else if v.lb.is_finite() {
- writeln!(gms, "v{i}.lo = {};", fmt(v.lb)).unwrap();
- }
- // Upper bound (+Inf is the GAMS default, only write when finite)
- if v.ub.is_finite() {
- writeln!(gms, "v{i}.up = {};", fmt(v.ub)).unwrap();
- }
- }
-
- // Initial levels (warm start)
- for v in vars.iter() {
- if let Some(val) = v.initial {
- writeln!(gms, "v{}.l = {};", v.id.index(), fmt(val)).unwrap();
- }
- }
- writeln!(gms).unwrap();
-
- // Equations declaration
- write!(gms, "Equations\n eq_obj").unwrap();
- for i in 0..constraints.len() {
- write!(gms, ", eq_c{i}").unwrap();
- }
- writeln!(gms, ";").unwrap();
- writeln!(gms).unwrap();
-
- // Objective equation: v_obj =e=
- write!(gms, "eq_obj.. v_obj =e=").unwrap();
- write_expr(&mut gms, &obj_terms, true);
- writeln!(gms, ";").unwrap();
-
- // Constraint equations: variable terms only, constant folded into RHS
- for (ci, (c, t)) in constraints.iter().zip(con_terms.iter()).enumerate() {
- let adjusted_rhs = c.rhs - t.constant;
- let sense_str = match c.sense {
- Sense::Le => "=l=",
- Sense::Ge => "=g=",
- Sense::Eq => "=e=",
- };
- write!(gms, "eq_c{ci}..").unwrap();
- write_expr(&mut gms, t, false);
- writeln!(gms, " {sense_str} {};", fmt(adjusted_rhs)).unwrap();
- }
- writeln!(gms).unwrap();
-
- // Options (time limit, MIP gap, sub-solver, etc.)
- write_options(&mut gms, opts, solve_type);
-
- // Model and solve statements
- writeln!(gms, "Model oximo_m / all /;").unwrap();
- if solver_opt.is_some() {
- writeln!(gms, "oximo_m.optfile = 1;").unwrap();
- }
- writeln!(gms, "Solve oximo_m using {solve_type} {sense_kw} v_obj;").unwrap();
- writeln!(gms).unwrap();
+ let (solve_type, solver_opt) = build_model_section(
+ &mut gms,
+ kind,
+ &arena,
+ &vars,
+ &constraints,
+ &objective,
+ sense_kw,
+ opts,
+ );
// - Temp directory
// Combine timestamp with a per-process atomic counter so concurrent
@@ -432,11 +295,232 @@ fn map_status(modelstat: i32, solvestat: i32) -> SolverStatus {
// - Helpers
+/// Write the formulation portion of the `.gms` file: title, variables, bounds,
+/// equations, options, model, and solve statement. Returns the solve type
+/// (`"LP"` / `"MIP"` / `"NLP"` / `"MINLP"` / `"QCP"` / `"MIQCP"`) and any
+/// solver-options file pair `(filename, content)` the caller should also
+/// persist alongside the `.gms`.
+#[allow(clippy::too_many_arguments)]
+fn build_model_section(
+ gms: &mut String,
+ kind: ModelKind,
+ arena: &ExprArena,
+ vars: &[Variable],
+ constraints: &[Constraint],
+ objective: &Objective,
+ sense_kw: &str,
+ opts: &GamsOptions,
+) -> (&'static str, Option<(String, String)>) {
+ let solve_type = gams_solve_type(kind);
+ let solver_opt = build_solver_opt(opts);
+
+ write_preamble(gms);
+ write_var_declarations(gms, vars);
+ write_bounds_and_initials(gms, vars);
+ write_equations(gms, arena, constraints, objective);
+ write_options(gms, opts, solve_type);
+ write_model_and_solve(gms, solve_type, sense_kw, solver_opt.is_some());
+
+ (solve_type, solver_opt)
+}
+
+pub(crate) fn gams_solve_type(kind: ModelKind) -> &'static str {
+ match kind {
+ ModelKind::LP => "LP",
+ ModelKind::MILP => "MIP",
+ ModelKind::QP => "QCP",
+ ModelKind::MIQP => "MIQCP",
+ ModelKind::NLP => "NLP",
+ ModelKind::MINLP => "MINLP",
+ }
+}
+
+/// Reject an explicitly selected sub-solver that cannot handle `kind` before
+/// invoking GAMS, so the caller gets a clear error naming the solver and model
+/// type instead of a downstream GAMS compilation failure.
+fn validate_solver(opts: &GamsOptions, kind: ModelKind) -> Result<(), SolverError> {
+ if let Some(cfg) = &opts.solver {
+ if !cfg.supports(kind) {
+ let solve_type = gams_solve_type(kind);
+ return Err(SolverError::Backend(format!(
+ "GAMS solver {} does not support {solve_type} models (model kind {kind:?}); \
+ select a solver that supports {solve_type}",
+ cfg.gams_name()
+ )));
+ }
+ }
+ Ok(())
+}
+
+fn build_solver_opt(opts: &GamsOptions) -> Option<(String, String)> {
+ opts.solver.as_ref().and_then(|cfg| {
+ let mut buf = String::new();
+ cfg.write_opt_file(&mut buf)
+ .then(|| (format!("{}.opt", cfg.gams_name().to_ascii_lowercase()), buf))
+ })
+}
+
+fn write_preamble(gms: &mut String) {
+ writeln!(gms, "$title oximo_model").unwrap();
+ writeln!(gms, "$offSymList").unwrap();
+ writeln!(gms, "$offSymXRef").unwrap();
+ writeln!(gms, "option solprint = off;").unwrap();
+ writeln!(gms, "option limrow = 0;").unwrap();
+ writeln!(gms, "option limcol = 0;").unwrap();
+ writeln!(gms).unwrap();
+}
+
+/// Emit `Variables`, `Binary Variables`, `Integer Variables` sections.
+fn write_var_declarations(gms: &mut String, vars: &[Variable]) {
+ let (mut cont, mut bin, mut int) = (Vec::new(), Vec::new(), Vec::new());
+ for v in vars {
+ match v.domain {
+ Domain::Binary => bin.push(v),
+ Domain::Integer | Domain::SemiInteger { .. } => int.push(v),
+ _ => cont.push(v),
+ }
+ }
+
+ write!(gms, "Variables\n v_obj").unwrap();
+ for v in &cont {
+ write!(gms, ", v{}", v.id.index()).unwrap();
+ }
+ writeln!(gms, ";").unwrap();
+
+ write_typed_var_section(gms, "Binary Variables", &bin);
+ write_typed_var_section(gms, "Integer Variables", &int);
+ writeln!(gms).unwrap();
+}
+
+fn write_typed_var_section(gms: &mut String, header: &str, vars: &[&Variable]) {
+ if vars.is_empty() {
+ return;
+ }
+ write!(gms, "{header}\n ").unwrap();
+ for (k, v) in vars.iter().enumerate() {
+ if k > 0 {
+ write!(gms, ", ").unwrap();
+ }
+ write!(gms, "v{}", v.id.index()).unwrap();
+ }
+ writeln!(gms, ";").unwrap();
+}
+
+fn write_bounds_and_initials(gms: &mut String, vars: &[Variable]) {
+ for v in vars {
+ write_var_bounds(gms, v);
+ }
+ for v in vars {
+ if let Some(val) = v.initial {
+ writeln!(gms, "v{}.l = {};", v.id.index(), fmt(val)).unwrap();
+ }
+ }
+ writeln!(gms).unwrap();
+}
+
+fn write_var_bounds(gms: &mut String, v: &Variable) {
+ let i = v.id.index();
+ if matches!(v.domain, Domain::Binary) {
+ // Default binary bounds are [0, 1], only emit when overridden or fixed.
+ if (v.lb - v.ub).abs() < f64::EPSILON {
+ writeln!(gms, "v{i}.fx = {};", fmt(v.lb)).unwrap();
+ return;
+ }
+ if v.lb.abs() > f64::EPSILON {
+ writeln!(gms, "v{i}.lo = {};", fmt(v.lb)).unwrap();
+ }
+ if (v.ub - 1.0).abs() > f64::EPSILON {
+ writeln!(gms, "v{i}.up = {};", fmt(v.ub)).unwrap();
+ }
+ return;
+ }
+ if v.lb == f64::NEG_INFINITY {
+ writeln!(gms, "v{i}.lo = -Inf;").unwrap();
+ } else if v.lb.is_finite() {
+ writeln!(gms, "v{i}.lo = {};", fmt(v.lb)).unwrap();
+ }
+ if v.ub.is_finite() {
+ writeln!(gms, "v{i}.up = {};", fmt(v.ub)).unwrap();
+ }
+}
+
+fn write_equations(
+ gms: &mut String,
+ arena: &ExprArena,
+ constraints: &[Constraint],
+ objective: &Objective,
+) {
+ write!(gms, "Equations\n eq_obj").unwrap();
+ for i in 0..constraints.len() {
+ write!(gms, ", eq_c{i}").unwrap();
+ }
+ writeln!(gms, ";").unwrap();
+ writeln!(gms).unwrap();
+
+ let obj_form = ExprForm::from(arena, objective.expr);
+ write!(gms, "eq_obj.. v_obj =e=").unwrap();
+ write_form(gms, arena, &obj_form, true);
+ writeln!(gms, ";").unwrap();
+
+ for (ci, c) in constraints.iter().enumerate() {
+ let sense_str = match c.sense {
+ Sense::Le => "=l=",
+ Sense::Ge => "=g=",
+ Sense::Eq => "=e=",
+ };
+ write!(gms, "eq_c{ci}..").unwrap();
+ match ExprForm::from(arena, c.lhs) {
+ ExprForm::Linear(t) => {
+ let adjusted_rhs = c.rhs - t.constant;
+ write_linear(gms, &t, false);
+ writeln!(gms, " {sense_str} {};", fmt(adjusted_rhs)).unwrap();
+ }
+ ExprForm::Nonlinear(id) => {
+ write_gams_expr(gms, arena, id, true);
+ writeln!(gms, " {sense_str} {};", fmt(c.rhs)).unwrap();
+ }
+ }
+ }
+ writeln!(gms).unwrap();
+}
+
+fn write_model_and_solve(gms: &mut String, solve_type: &str, sense_kw: &str, has_opt: bool) {
+ writeln!(gms, "Model oximo_m / all /;").unwrap();
+ if has_opt {
+ writeln!(gms, "oximo_m.optfile = 1;").unwrap();
+ }
+ writeln!(gms, "Solve oximo_m using {solve_type} {sense_kw} v_obj;").unwrap();
+ writeln!(gms).unwrap();
+}
+
+/// Captured form of an expression for GAMS emission.
+enum ExprForm {
+ Linear(LinearTerms),
+ Nonlinear(ExprId),
+}
+
+impl ExprForm {
+ fn from(arena: &ExprArena, id: ExprId) -> Self {
+ match extract_linear(arena, id) {
+ Some(t) => ExprForm::Linear(t),
+ None => ExprForm::Nonlinear(id),
+ }
+ }
+}
+
+/// Append a captured expression form to `gms`.
+fn write_form(gms: &mut String, arena: &ExprArena, form: &ExprForm, include_constant: bool) {
+ match form {
+ ExprForm::Linear(t) => write_linear(gms, t, include_constant),
+ ExprForm::Nonlinear(id) => write_gams_expr(gms, arena, *id, true),
+ }
+}
+
/// Append the linear expression `t` to `gms`.
/// When `include_constant` is true, the constant term is included; otherwise
/// only variable terms are emitted (used for constraints where the constant is
/// folded into the RHS).
-fn write_expr(gms: &mut String, t: &LinearTerms, include_constant: bool) {
+fn write_linear(gms: &mut String, t: &LinearTerms, include_constant: bool) {
let mut first = true;
for (v, coef) in &t.coeffs {
if *coef == 0.0 {
@@ -467,6 +551,101 @@ fn write_expr(gms: &mut String, t: &LinearTerms, include_constant: bool) {
}
}
+/// Recursive infix printer for a GAMS-compatible expression.
+fn write_gams_expr(gms: &mut String, arena: &ExprArena, id: ExprId, leading_space: bool) {
+ if leading_space {
+ write!(gms, " ").unwrap();
+ }
+ match arena.get(id) {
+ ExprNode::Const(c) => write!(gms, "{}", fmt(*c)).unwrap(),
+ ExprNode::Var(v) => write!(gms, "v{}", v.index()).unwrap(),
+ ExprNode::Param(_) => {
+ // Since params not yet passed into GAMS emission, we emit a placeholder
+ // so downstream errors are clear.
+ write!(gms, "0 /* unsupported: param */").unwrap();
+ }
+ ExprNode::Linear { coeffs, constant } => {
+ let t = LinearTerms { coeffs: coeffs.clone(), constant: *constant };
+ write!(gms, "(").unwrap();
+ write_linear(gms, &t, true);
+ write!(gms, " )").unwrap();
+ }
+ ExprNode::Neg(inner) => {
+ write!(gms, "(-").unwrap();
+ write_gams_expr(gms, arena, *inner, true);
+ write!(gms, ")").unwrap();
+ }
+ ExprNode::Add(children) => {
+ write!(gms, "(").unwrap();
+ for (i, c) in children.iter().enumerate() {
+ if i > 0 {
+ write!(gms, " +").unwrap();
+ }
+ write_gams_expr(gms, arena, *c, true);
+ }
+ write!(gms, ")").unwrap();
+ }
+ ExprNode::Mul(children) => {
+ write!(gms, "(").unwrap();
+ for (i, c) in children.iter().enumerate() {
+ if i > 0 {
+ write!(gms, " *").unwrap();
+ }
+ write_gams_expr(gms, arena, *c, true);
+ }
+ write!(gms, ")").unwrap();
+ }
+ ExprNode::Pow(base, exp) => {
+ // GAMS's `**` lowers to `rPower(x, r)`, which rejects negative
+ // bases. For small integer constant exponents emit `power(x, n)`
+ // (accepts any real base), otherwise fall back to `**`.
+ //
+ // The 1e9 cap keeps the cast safe and rejects nonsense huge exponents
+ // that would still satisfy the integer check after f64 rounding.
+ if let ExprNode::Const(c) = arena.get(*exp) {
+ if (c - c.round()).abs() < f64::EPSILON && c.abs() <= 1e9 {
+ write!(gms, "power(").unwrap();
+ write_gams_expr(gms, arena, *base, false);
+ write!(gms, ", {:.0})", c.round()).unwrap();
+ return;
+ }
+ }
+ write!(gms, "(").unwrap();
+ write_gams_expr(gms, arena, *base, false);
+ write!(gms, " **").unwrap();
+ write_gams_expr(gms, arena, *exp, true);
+ write!(gms, ")").unwrap();
+ }
+ ExprNode::Div(num, den) => {
+ write!(gms, "(").unwrap();
+ write_gams_expr(gms, arena, *num, false);
+ write!(gms, " /").unwrap();
+ write_gams_expr(gms, arena, *den, true);
+ write!(gms, ")").unwrap();
+ }
+ ExprNode::Sin(a) => {
+ write!(gms, "sin(").unwrap();
+ write_gams_expr(gms, arena, *a, false);
+ write!(gms, ")").unwrap();
+ }
+ ExprNode::Cos(a) => {
+ write!(gms, "cos(").unwrap();
+ write_gams_expr(gms, arena, *a, false);
+ write!(gms, ")").unwrap();
+ }
+ ExprNode::Exp(a) => {
+ write!(gms, "exp(").unwrap();
+ write_gams_expr(gms, arena, *a, false);
+ write!(gms, ")").unwrap();
+ }
+ ExprNode::Log(a) => {
+ write!(gms, "log(").unwrap();
+ write_gams_expr(gms, arena, *a, false);
+ write!(gms, ")").unwrap();
+ }
+ }
+}
+
/// Format an `f64` for use in a GAMS file.
fn fmt(v: f64) -> String {
if v == f64::INFINITY {
@@ -481,12 +660,10 @@ fn fmt(v: f64) -> String {
/// Parse a GAMS-formatted integer (may be written as `"1"` or `"1.000"`).
fn parse_gams_int(s: &str) -> Option {
let trimmed = s.trim();
- if let Ok(n) = trimmed.parse::() {
- return Some(n);
- }
- // Fall back through f64 for GAMS formats like "1.000".
- #[allow(clippy::cast_possible_truncation)]
- trimmed.parse::().ok().map(|f| f.round() as i32)
+ // GAMS writes modelstat/solvestat with the `:0:0` PUT format, so we
+ // normally see a bare integer.
+ let head = trimmed.split_once('.').map_or(trimmed, |(int, _)| int);
+ head.parse::().ok()
}
/// Parse a GAMS-formatted float, tolerating GAMS special tokens.
@@ -498,3 +675,101 @@ fn parse_gams_float(s: &str) -> Option {
other => other.parse().ok(),
}
}
+
+#[cfg(test)]
+mod tests {
+ use super::*;
+ use oximo_core::prelude::*;
+
+ fn render(model: &Model, opts: &GamsOptions) -> String {
+ let arena = model.arena();
+ let vars = model.variables();
+ let constraints = model.constraints();
+ let objective = model.try_objective().expect("objective set");
+ let sense_kw = match objective.sense {
+ ObjectiveSense::Minimize => "minimizing",
+ ObjectiveSense::Maximize => "maximizing",
+ };
+ let mut gms = String::new();
+ build_model_section(
+ &mut gms,
+ model.kind(),
+ &arena,
+ &vars,
+ &constraints,
+ &objective,
+ sense_kw,
+ opts,
+ );
+ gms
+ }
+
+ #[test]
+ fn linear_objective_uses_lp_solve_type() {
+ let m = Model::new("lp");
+ let x = m.var("x").lb(0.0).ub(10.0).build();
+ let y = m.var("y").lb(0.0).ub(10.0).build();
+ m.constraint("c", (x + y).le(5.0));
+ m.minimize(x + 2.0 * y);
+ let gms = render(&m, &GamsOptions::default());
+ assert!(gms.contains("Solve oximo_m using LP minimizing v_obj;"), "got:\n{gms}");
+ }
+
+ #[test]
+ fn nlp_uses_transcendental_and_picks_nlp_solve_type() {
+ let m = Model::new("nlp");
+ let x = m.var("x").lb(-std::f64::consts::PI).ub(std::f64::consts::PI).build();
+ m.minimize(x.sin() + x.exp());
+ let gms = render(&m, &GamsOptions::default());
+ assert!(gms.contains("Solve oximo_m using NLP minimizing v_obj;"), "got:\n{gms}");
+ assert!(gms.contains("sin("), "expected sin(...) in objective:\n{gms}");
+ assert!(gms.contains("exp("), "expected exp(...) in objective:\n{gms}");
+ }
+
+ #[test]
+ fn minlp_nonlinear_knapsack_routes_to_minlp_solve_type() {
+ let m = Model::new("minlp");
+ let x = m.var("x").binary().build();
+ let y = m.var("y").lb(0.0).ub(10.0).build();
+ m.constraint("budget", (x + y).le(8.0));
+ let one = Expr::constant(x.arena, 1.0);
+ m.maximize((one + y).log() + 2.0 * x);
+ let gms = render(&m, &GamsOptions::default());
+ assert!(gms.contains("Solve oximo_m using MINLP maximizing v_obj;"), "got:\n{gms}");
+ assert!(gms.contains("log("), "expected log(...) in objective:\n{gms}");
+ }
+
+ #[test]
+ fn quadratic_constraint_emits_full_expression_against_rhs() {
+ let m = Model::new("qcp");
+ let x = m.var("x").lb(0.0).ub(5.0).build();
+ let y = m.var("y").lb(0.0).ub(5.0).build();
+ m.constraint("xy", (x * y).le(4.0));
+ m.minimize(x + y);
+ let gms = render(&m, &GamsOptions::default());
+ assert!(gms.contains("Solve oximo_m using QCP minimizing v_obj;"), "got:\n{gms}");
+ // The product term must appear on the LHS, the user RHS untouched.
+ assert!(gms.contains("v0") && gms.contains("v1"), "vars missing:\n{gms}");
+ assert!(gms.contains("=l= 4"), "expected =l= 4 on the right:\n{gms}");
+ }
+
+ #[test]
+ fn integer_power_uses_power_func() {
+ let m = Model::new("pow");
+ let x = m.var("x").lb(-10.0).ub(10.0).build();
+ m.minimize(x.powi(3));
+ let gms = render(&m, &GamsOptions::default());
+ assert!(gms.contains("power("), "expected power(...) for int Pow:\n{gms}");
+ assert!(gms.contains(", 3)"), "expected exponent 3:\n{gms}");
+ assert!(gms.contains("Solve oximo_m using NLP minimizing v_obj;"), "got:\n{gms}");
+ }
+
+ #[test]
+ fn real_power_falls_back_to_double_star() {
+ let m = Model::new("rpow");
+ let x = m.var("x").lb(0.1).ub(10.0).build();
+ m.minimize(x.powf(0.5));
+ let gms = render(&m, &GamsOptions::default());
+ assert!(gms.contains(" **"), "expected ** for real Pow:\n{gms}");
+ }
+}
diff --git a/crates/oximo-gurobi/Cargo.toml b/crates/oximo-gurobi/Cargo.toml
index d0b524d..260feca 100644
--- a/crates/oximo-gurobi/Cargo.toml
+++ b/crates/oximo-gurobi/Cargo.toml
@@ -14,7 +14,6 @@ oximo-expr.workspace = true
oximo-core.workspace = true
oximo-solver.workspace = true
grb = { workspace = true, features = ["gurobi12"] }
-rayon.workspace = true
rustc-hash.workspace = true
thiserror.workspace = true
diff --git a/crates/oximo-gurobi/README.md b/crates/oximo-gurobi/README.md
index 8a4eb7e..6f6b3e1 100644
--- a/crates/oximo-gurobi/README.md
+++ b/crates/oximo-gurobi/README.md
@@ -1,8 +1,8 @@
# oximo-gurobi
-Gurobi LP/MILP backend for [oximo](https://github.com/germanheim/oximo).
+Gurobi backend for [oximo](https://github.com/germanheim/oximo).
-Wraps the [`grb`](https://crates.io/crates/grb) crate (`gurobi12` feature). Supports `LP` and `MILP` model kinds. **QP/NLP return `SolverError::UnsupportedKind` for now**.
+Wraps the [`grb`](https://crates.io/crates/grb) crate (`gurobi12` feature). Supports `LP`, `MILP`, `QP`, `MIQP`, `NLP`, and `MINLP` model kinds. Nonlinear expressions (`Pow`, `Sin`, `Cos`, `Exp`, `Log`, bilinear `Mul`) are lowered to auxiliary variables wired together with Gurobi's native `add_genconstr_*` and `add_qconstr` APIs. `NonConvex = 2` is enabled automatically when needed.
Gurobi v13.0+ works with the `gurobi12` feature, which is the default for this crate. For more information see [grb/issue#31](https://github.com/ykrist/rust-grb/issues/31). To use older versions, enable the `gurobi11` feature instead (see [grb docs](https://docs.rs/grb/latest/grb/#features) for details).
diff --git a/crates/oximo-gurobi/src/lib.rs b/crates/oximo-gurobi/src/lib.rs
index c0b3ea2..f25e6ae 100644
--- a/crates/oximo-gurobi/src/lib.rs
+++ b/crates/oximo-gurobi/src/lib.rs
@@ -1,6 +1,7 @@
#![doc = include_str!("../README.md")]
#![forbid(unsafe_code)]
+mod nonlinear;
mod options;
mod translate;
@@ -21,7 +22,15 @@ impl Solver for Gurobi {
}
fn supports(&self, kind: ModelKind) -> bool {
- matches!(kind, ModelKind::LP | ModelKind::MILP)
+ matches!(
+ kind,
+ ModelKind::LP
+ | ModelKind::MILP
+ | ModelKind::QP
+ | ModelKind::MIQP
+ | ModelKind::NLP
+ | ModelKind::MINLP
+ )
}
fn solve(&mut self, model: &Model, opts: &GurobiOptions) -> Result {
diff --git a/crates/oximo-gurobi/src/nonlinear.rs b/crates/oximo-gurobi/src/nonlinear.rs
new file mode 100644
index 0000000..31de1ed
--- /dev/null
+++ b/crates/oximo-gurobi/src/nonlinear.rs
@@ -0,0 +1,440 @@
+//! Lower an oximo expression tree onto Gurobi's native types.
+//!
+//! Gurobi 12 supports linear, quadratic, and a fixed set of general nonlinear
+//! function constraints (exp/log/sin/cos/pow/...) where each one takes the form
+//! `y = f(x)` between two variables. Anything more complex must be flattened
+//! into a DAG of auxiliary variables connected by those primitive equalities.
+//! This module is that flattening pass.
+
+// TODO: This will change once we can get support for V13 in the grb crate
+// since we will be able to use generic expresion trees
+
+use grb::expr::{LinExpr, QuadExpr};
+use grb::prelude::*;
+use oximo_expr::{ExprArena, ExprId, ExprNode, VarId};
+use oximo_solver::SolverError;
+
+pub(crate) enum LoweredExpr {
+ Linear(LinExpr),
+ Quadratic(QuadExpr),
+ Var(Var),
+}
+
+pub(crate) struct LoweringCtx<'a> {
+ pub model: &'a mut Model,
+ pub gurobi_vars: &'a [Var],
+ pub aux_counter: u32,
+}
+
+impl LoweringCtx<'_> {
+ fn next_name(&mut self, tag: &str) -> String {
+ let n = self.aux_counter;
+ self.aux_counter += 1;
+ format!("aux_{tag}_{n}")
+ }
+
+ fn new_aux(&mut self, tag: &str, lb: f64, ub: f64) -> Result {
+ let name = self.next_name(tag);
+ let m = &mut *self.model;
+ #[allow(clippy::unnecessary_cast)]
+ let v = add_ctsvar!(m, name: &name, bounds: lb..ub)?;
+ Ok(v)
+ }
+}
+
+fn map_grb(e: grb::Error) -> SolverError {
+ SolverError::Backend(e.to_string())
+}
+
+fn linear_from_var(v: Var) -> LinExpr {
+ let mut e = LinExpr::new();
+ e.add_term(1.0, v);
+ e
+}
+
+fn linear_constant(c: f64) -> LinExpr {
+ let mut e = LinExpr::new();
+ e.add_constant(c);
+ e
+}
+
+fn quad_from_linear(e: LinExpr) -> QuadExpr {
+ let mut q = QuadExpr::new();
+ q.add_constant(e.get_offset());
+ for (v, c) in e.into_parts().0 {
+ q.add_term(c, v);
+ }
+ q
+}
+
+fn add_linears(mut a: LinExpr, b: LinExpr) -> LinExpr {
+ let (coeffs, offset) = b.into_parts();
+ a.add_constant(offset);
+ for (v, c) in coeffs {
+ a.add_term(c, v);
+ }
+ a
+}
+
+fn add_quads(mut a: QuadExpr, b: QuadExpr) -> QuadExpr {
+ let (qcoeffs, linexpr) = b.into_parts();
+ a.add_constant(linexpr.get_offset());
+ for (v, c) in linexpr.into_parts().0 {
+ a.add_term(c, v);
+ }
+ for ((x, y), c) in qcoeffs {
+ a.add_qterm(c, x, y);
+ }
+ a
+}
+
+fn scale_linear(mut e: LinExpr, k: f64) -> LinExpr {
+ e.mul_scalar(k);
+ e
+}
+
+fn scale_quad(mut e: QuadExpr, k: f64) -> QuadExpr {
+ e.mul_scalar(k);
+ e
+}
+
+/// Convert any lowered form to a `LinExpr`.
+///
+/// Panics if quadratic, callers must only invoke when the value
+/// is known linear (i.e. degree <= 1).
+fn into_linexpr(l: LoweredExpr) -> LinExpr {
+ match l {
+ LoweredExpr::Linear(e) => e,
+ LoweredExpr::Var(v) => linear_from_var(v),
+ LoweredExpr::Quadratic(_) => {
+ panic!("internal: into_linexpr called on Quadratic LoweredExpr")
+ }
+ }
+}
+
+/// Combine two lowered values additively, promoting to the highest order.
+fn lowered_add(a: LoweredExpr, b: LoweredExpr) -> LoweredExpr {
+ use LoweredExpr::{Linear, Quadratic, Var};
+ match (a, b) {
+ (Linear(x), Linear(y)) => Linear(add_linears(x, y)),
+ (Var(v), Linear(y)) | (Linear(y), Var(v)) => Linear(add_linears(y, linear_from_var(v))),
+ (Var(v), Var(w)) => Linear(add_linears(linear_from_var(v), linear_from_var(w))),
+ (Quadratic(x), Quadratic(y)) => Quadratic(add_quads(x, y)),
+ (Quadratic(x), other) | (other, Quadratic(x)) => {
+ let y = into_linexpr(other);
+ Quadratic(add_quads(x, quad_from_linear(y)))
+ }
+ }
+}
+
+fn lowered_scale(l: LoweredExpr, k: f64) -> LoweredExpr {
+ use LoweredExpr::{Linear, Quadratic, Var};
+ if k == 0.0 {
+ return Linear(linear_constant(0.0));
+ }
+ match l {
+ Linear(e) => Linear(scale_linear(e, k)),
+ Quadratic(e) => Quadratic(scale_quad(e, k)),
+ Var(v) => Linear(scale_linear(linear_from_var(v), k)),
+ }
+}
+
+fn lowered_neg(l: LoweredExpr) -> LoweredExpr {
+ lowered_scale(l, -1.0)
+}
+
+/// Materialize `lowered` into a single Gurobi variable, returning the existing
+/// one if it is already opaque. Bounds default to +-inf, callers should tighten
+/// when they know more about the value's range.
+fn as_aux_var(lowered: LoweredExpr, ctx: &mut LoweringCtx<'_>) -> Result {
+ match lowered {
+ LoweredExpr::Var(v) => Ok(v),
+ LoweredExpr::Linear(e) => {
+ let aux = ctx.new_aux("v", f64::NEG_INFINITY, f64::INFINITY).map_err(map_grb)?;
+ let name = ctx.next_name("eq");
+ ctx.model.add_constr(&name, c!(aux == e)).map_err(map_grb)?;
+ Ok(aux)
+ }
+ LoweredExpr::Quadratic(e) => {
+ let aux = ctx.new_aux("v", f64::NEG_INFINITY, f64::INFINITY).map_err(map_grb)?;
+ let name = ctx.next_name("qeq");
+ ctx.model.add_qconstr(&name, c!(aux == e)).map_err(map_grb)?;
+ Ok(aux)
+ }
+ }
+}
+
+/// True if `id` is a literal constant, returns the value if so.
+fn as_const(arena: &ExprArena, id: ExprId) -> Option {
+ match arena.get(id) {
+ ExprNode::Const(c) => Some(*c),
+ _ => None,
+ }
+}
+
+pub(crate) fn lower(
+ arena: &ExprArena,
+ id: ExprId,
+ ctx: &mut LoweringCtx<'_>,
+) -> Result {
+ match arena.get(id) {
+ ExprNode::Const(c) => Ok(LoweredExpr::Linear(linear_constant(*c))),
+ ExprNode::Var(v) => Ok(LoweredExpr::Linear(linear_from_var(grb_var(ctx, *v)))),
+ ExprNode::Param(_) => Err(SolverError::Backend(
+ "Param nodes are not supported by oximo-gurobi nonlinear lowering".into(),
+ )),
+ ExprNode::Linear { coeffs, constant } => {
+ let mut e = linear_constant(*constant);
+ for (v, c) in coeffs {
+ e.add_term(*c, grb_var(ctx, *v));
+ }
+ Ok(LoweredExpr::Linear(e))
+ }
+ ExprNode::Neg(inner) => {
+ let l = lower(arena, *inner, ctx)?;
+ Ok(lowered_neg(l))
+ }
+ ExprNode::Add(children) => {
+ let mut acc = LoweredExpr::Linear(linear_constant(0.0));
+ for c in children {
+ let l = lower(arena, *c, ctx)?;
+ acc = lowered_add(acc, l);
+ }
+ Ok(acc)
+ }
+ ExprNode::Mul(children) => lower_mul(arena, children, ctx),
+ ExprNode::Pow(base, exp) => lower_pow(arena, *base, *exp, ctx),
+ ExprNode::Div(num, den) => lower_div(arena, *num, *den, ctx),
+ ExprNode::Sin(a) => lower_unary(arena, *a, ctx, UnaryFn::Sin),
+ ExprNode::Cos(a) => lower_unary(arena, *a, ctx, UnaryFn::Cos),
+ ExprNode::Exp(a) => lower_unary(arena, *a, ctx, UnaryFn::Exp),
+ ExprNode::Log(a) => lower_unary(arena, *a, ctx, UnaryFn::Log),
+ }
+}
+
+fn grb_var(ctx: &LoweringCtx<'_>, v: VarId) -> Var {
+ ctx.gurobi_vars[v.index()]
+}
+
+fn lower_mul(
+ arena: &ExprArena,
+ children: &[ExprId],
+ ctx: &mut LoweringCtx<'_>,
+) -> Result {
+ let mut scalar = 1.0_f64;
+ let mut non_consts: Vec = Vec::new();
+ for c in children {
+ if let Some(k) = as_const(arena, *c) {
+ scalar *= k;
+ } else {
+ non_consts.push(*c);
+ }
+ }
+ if non_consts.is_empty() {
+ return Ok(LoweredExpr::Linear(linear_constant(scalar)));
+ }
+ if non_consts.len() == 1 {
+ let l = lower(arena, non_consts[0], ctx)?;
+ return Ok(lowered_scale(l, scalar));
+ }
+ if non_consts.len() == 2 {
+ let a = lower(arena, non_consts[0], ctx)?;
+ let b = lower(arena, non_consts[1], ctx)?;
+ // Linear * Linear -> quadratic, if either side has variable terms.
+ if let (LoweredExpr::Linear(la), LoweredExpr::Linear(lb)) = (&a, &b) {
+ let q = multiply_linears(la, lb, scalar);
+ return Ok(LoweredExpr::Quadratic(q));
+ }
+ // Mixed with Quadratic or Var -> materialize aux vars and recompose.
+ let va = as_aux_var(a, ctx)?;
+ let vb = as_aux_var(b, ctx)?;
+ let mut q = QuadExpr::new();
+ q.add_qterm(scalar, va, vb);
+ return Ok(LoweredExpr::Quadratic(q));
+ }
+ // 3+ non-constant factors: degree > 2. Fold left, materializing aux vars.
+ let mut acc_var = {
+ let a = lower(arena, non_consts[0], ctx)?;
+ let b = lower(arena, non_consts[1], ctx)?;
+ let va = as_aux_var(a, ctx)?;
+ let vb = as_aux_var(b, ctx)?;
+ let mut q = QuadExpr::new();
+ q.add_qterm(1.0, va, vb);
+ as_aux_var(LoweredExpr::Quadratic(q), ctx)?
+ };
+ for c in &non_consts[2..] {
+ let next = lower(arena, *c, ctx)?;
+ let vn = as_aux_var(next, ctx)?;
+ let mut q = QuadExpr::new();
+ q.add_qterm(1.0, acc_var, vn);
+ acc_var = as_aux_var(LoweredExpr::Quadratic(q), ctx)?;
+ }
+ Ok(lowered_scale(LoweredExpr::Var(acc_var), scalar))
+}
+
+fn multiply_linears(a: &LinExpr, b: &LinExpr, scalar: f64) -> QuadExpr {
+ let a_off = a.get_offset();
+ let b_off = b.get_offset();
+ let mut q = QuadExpr::new();
+ q.add_constant(scalar * a_off * b_off);
+ for (va, ca) in a.iter_terms() {
+ q.add_term(scalar * ca * b_off, *va);
+ }
+ for (vb, cb) in b.iter_terms() {
+ q.add_term(scalar * a_off * cb, *vb);
+ }
+ for (va, ca) in a.iter_terms() {
+ for (vb, cb) in b.iter_terms() {
+ q.add_qterm(scalar * ca * cb, *va, *vb);
+ }
+ }
+ q
+}
+
+fn lower_pow(
+ arena: &ExprArena,
+ base: ExprId,
+ exp: ExprId,
+ ctx: &mut LoweringCtx<'_>,
+) -> Result {
+ // Constant exponent -> either expand to a Mul chain (small ints) or use
+ // Gurobi's add_genconstr_pow. Non-constant exponent -> exp(b*log(a)).
+ if let Some(alpha) = as_const(arena, exp) {
+ if alpha == 0.0 {
+ return Ok(LoweredExpr::Linear(linear_constant(1.0)));
+ }
+ if (alpha - 1.0).abs() < f64::EPSILON {
+ return lower(arena, base, ctx);
+ }
+ if (alpha - alpha.round()).abs() < f64::EPSILON && alpha > 0.0 && alpha <= 4.0 {
+ // Pre-checks guarantee alpha in {1.0, 2.0, 3.0, 4.0}, alpha == 1 was
+ // already handled above, so bucket the remaining three values.
+ let n: u32 = match alpha.round() {
+ v if v < 2.5 => 2,
+ v if v < 3.5 => 3,
+ _ => 4,
+ };
+ let base_lowered = lower(arena, base, ctx)?;
+ let v = as_aux_var(base_lowered, ctx)?;
+ let mut q = QuadExpr::new();
+ q.add_qterm(1.0, v, v);
+ let mut acc = LoweredExpr::Quadratic(q);
+ for _ in 2..n {
+ let aux = as_aux_var(acc, ctx)?;
+ let mut next = QuadExpr::new();
+ next.add_qterm(1.0, aux, v);
+ acc = LoweredExpr::Quadratic(next);
+ }
+ return Ok(acc);
+ }
+ // General constant exponent via Gurobi's genconstr_pow
+ let base_lowered = lower(arena, base, ctx)?;
+ let x = as_aux_var(base_lowered, ctx)?;
+ let y = ctx.new_aux("pow", f64::NEG_INFINITY, f64::INFINITY).map_err(map_grb)?;
+ let name = ctx.next_name("pow");
+ ctx.model.add_genconstr_pow(&name, x, y, alpha, "").map_err(map_grb)?;
+ return Ok(LoweredExpr::Var(y));
+ }
+ // exp(b*log(a))
+ let base_lowered = lower(arena, base, ctx)?;
+ let x = as_aux_var(base_lowered, ctx)?;
+ let log_y = ctx.new_aux("log", f64::NEG_INFINITY, f64::INFINITY).map_err(map_grb)?;
+ let log_name = ctx.next_name("log");
+ ctx.model.add_genconstr_natural_log(&log_name, x, log_y, "").map_err(map_grb)?;
+ let b_lowered = lower(arena, exp, ctx)?;
+ let b_var = as_aux_var(b_lowered, ctx)?;
+ let prod = ctx.new_aux("mul", f64::NEG_INFINITY, f64::INFINITY).map_err(map_grb)?;
+ let mut q = QuadExpr::new();
+ q.add_qterm(1.0, b_var, log_y);
+ let prod_eq = ctx.next_name("mul_eq");
+ ctx.model.add_qconstr(&prod_eq, c!(prod == q)).map_err(map_grb)?;
+ let result = ctx.new_aux("exp", 0.0, f64::INFINITY).map_err(map_grb)?;
+ let exp_name = ctx.next_name("exp");
+ ctx.model.add_genconstr_natural_exp(&exp_name, prod, result, "").map_err(map_grb)?;
+ Ok(LoweredExpr::Var(result))
+}
+
+/// Lower `num / den` as `num * recip`, introducing a reciprocal variable
+/// `recip` pinned by the bilinear equality `den * recip == 1`.
+///
+/// This avoids a `pow(den, -1)` lowering. Gurobi's `genconstr_pow` has a pole
+/// at 0 and so requires its base to stay in a strictly positive domain,
+/// which rules out negative or zero-straddling denominators. The bilinear pin
+/// carries no such restriction: it is a plain non-convex quadratic
+/// (handled via `NonConvex = 2`) valid for `den` of either sign,
+/// and is infeasible only at `den == 0`.
+fn lower_div(
+ arena: &ExprArena,
+ num: ExprId,
+ den: ExprId,
+ ctx: &mut LoweringCtx<'_>,
+) -> Result {
+ // `div_into` folds every nonzero constant denominator into the linear path
+ // at construction, so the only constant `den` that reaches here is a literal
+ // zero.
+ if as_const(arena, den) == Some(0.0) {
+ return Err(SolverError::Backend("division by zero: constant denominator is 0".into()));
+ }
+
+ // recip = 1 / den, pinned by `den * recip == 1`.
+ let den_lowered = lower(arena, den, ctx)?;
+ let x = as_aux_var(den_lowered, ctx)?;
+ let recip = ctx.new_aux("recip", f64::NEG_INFINITY, f64::INFINITY).map_err(map_grb)?;
+ let mut pin = QuadExpr::new();
+ pin.add_qterm(1.0, x, recip);
+ let name = ctx.next_name("recip_eq");
+ ctx.model.add_qconstr(&name, c!(pin == 1.0)).map_err(map_grb)?;
+
+ // Constant numerator stays linear in `recip`, otherwise form `num * recip`.
+ if let Some(k) = as_const(arena, num) {
+ return Ok(lowered_scale(LoweredExpr::Var(recip), k));
+ }
+ let num_lowered = lower(arena, num, ctx)?;
+ let vn = as_aux_var(num_lowered, ctx)?;
+ let mut q = QuadExpr::new();
+ q.add_qterm(1.0, vn, recip);
+ Ok(LoweredExpr::Quadratic(q))
+}
+
+enum UnaryFn {
+ Sin,
+ Cos,
+ Exp,
+ Log,
+}
+
+fn lower_unary(
+ arena: &ExprArena,
+ inner: ExprId,
+ ctx: &mut LoweringCtx<'_>,
+ f: UnaryFn,
+) -> Result {
+ let lowered = lower(arena, inner, ctx)?;
+ let x = as_aux_var(lowered, ctx)?;
+ let (lb, ub, tag) = match f {
+ UnaryFn::Sin | UnaryFn::Cos => (-1.0, 1.0, "trig"),
+ UnaryFn::Exp => (0.0, f64::INFINITY, "exp"),
+ UnaryFn::Log => (f64::NEG_INFINITY, f64::INFINITY, "log"),
+ };
+ let y = ctx.new_aux(tag, lb, ub).map_err(map_grb)?;
+ let name = ctx.next_name(tag);
+ match f {
+ UnaryFn::Sin => ctx.model.add_genconstr_sin(&name, x, y, "").map_err(map_grb)?,
+ UnaryFn::Cos => ctx.model.add_genconstr_cos(&name, x, y, "").map_err(map_grb)?,
+ UnaryFn::Exp => ctx.model.add_genconstr_natural_exp(&name, x, y, "").map_err(map_grb)?,
+ UnaryFn::Log => ctx.model.add_genconstr_natural_log(&name, x, y, "").map_err(map_grb)?,
+ };
+ Ok(LoweredExpr::Var(y))
+}
+
+/// Helpers used by translate.rs to materialize a lowered constraint or
+/// objective expression.
+impl LoweredExpr {
+ pub(crate) fn into_expr_for_objective(self) -> grb::expr::Expr {
+ match self {
+ LoweredExpr::Linear(e) => grb::expr::Expr::from(e),
+ LoweredExpr::Quadratic(e) => grb::expr::Expr::from(e),
+ LoweredExpr::Var(v) => grb::expr::Expr::from(v),
+ }
+ }
+}
diff --git a/crates/oximo-gurobi/src/translate.rs b/crates/oximo-gurobi/src/translate.rs
index f84cf4e..3de03d9 100644
--- a/crates/oximo-gurobi/src/translate.rs
+++ b/crates/oximo-gurobi/src/translate.rs
@@ -3,12 +3,12 @@ use std::time::Instant;
use grb::expr::LinExpr;
use grb::prelude::*;
use oximo_core::{ConstraintId, Domain, Model, ModelKind, ObjectiveSense, Sense, VarId};
-use oximo_expr::{ExprArena, LinearTerms, extract_linear};
+use oximo_expr::{ExprArena, ExprId, LinearTerms, extract_linear};
use oximo_solver::{SolverError, SolverResult, SolverStatus};
-use rayon::prelude::*;
use rustc_hash::FxHashMap;
use crate::GurobiOptions;
+use crate::nonlinear::{LoweredExpr, LoweringCtx, lower};
use crate::options::apply as apply_options;
fn map_grb_err(e: grb::Error) -> SolverError {
@@ -21,34 +21,22 @@ fn map_grb_err(e: grb::Error) -> SolverError {
/// # Errors
///
/// Returns a [`SolverError`] if the model is unsupported, contains nonlinear
-/// expressions, or if Gurobi reports an error during setup or optimization.
+/// expressions Gurobi cannot represent, or if Gurobi reports an error during
+/// setup or optimization.
///
/// # Panics
///
/// Panics if model variable or constraint indices overflow `u32`.
-#[allow(clippy::unnecessary_cast, clippy::cast_possible_truncation, clippy::cast_sign_loss)]
pub fn solve(model: &Model, opts: &GurobiOptions) -> Result {
let kind = model.kind();
- if !matches!(kind, ModelKind::LP | ModelKind::MILP) {
- return Err(SolverError::UnsupportedKind(kind));
- }
+ let nonlinear_kind =
+ matches!(kind, ModelKind::QP | ModelKind::MIQP | ModelKind::NLP | ModelKind::MINLP);
let arena = model.arena();
let vars = model.variables();
let constraints = model.constraints();
let objective = model.try_objective().map_err(SolverError::Core)?;
- let mut obj_coeffs = vec![0.0; vars.len()];
- let obj_constant = match extract_linear(&arena, objective.expr) {
- Some(t) => {
- for (v, c) in t.coeffs {
- obj_coeffs[v.index()] = c;
- }
- t.constant
- }
- None => return Err(SolverError::Nonlinear),
- };
-
let env = Env::new("").map_err(|e| SolverError::Backend(format!("Gurobi env: {e}")))?;
let mut grb_model = grb::Model::with_env("oximo", &env).map_err(map_grb_err)?;
@@ -61,7 +49,9 @@ pub fn solve(model: &Model, opts: &GurobiOptions) -> Result VarType::SemiCont,
Domain::SemiInteger { .. } => VarType::SemiInt,
};
- let gvar = add_var!(grb_model, vtype, obj: obj_coeffs[i], bounds: v.lb..v.ub, name: &format!("x{i}"))
+ // `add_var!` expands the f64 bounds with an `as f64` cast.
+ #[allow(clippy::unnecessary_cast)]
+ let gvar = add_var!(grb_model, vtype, bounds: v.lb..v.ub, name: &format!("x{i}"))
.map_err(map_grb_err)?;
gurobi_vars.push(gvar);
if let Some(val) = v.initial {
@@ -70,46 +60,69 @@ pub fn solve(model: &Model, opts: &GurobiOptions) -> Result = constraints
- .par_iter()
- .map(|c| extract_linear(arena_ref, c.lhs).ok_or(SolverError::Nonlinear))
- .collect::, _>>()?;
-
- let mut gurobi_constrs = Vec::with_capacity(constraints.len());
- for (c_id, (c, t)) in constraints.iter().zip(con_terms).enumerate() {
- let adjusted_rhs = c.rhs - t.constant;
-
- let mut expr = LinExpr::new();
- for (v, co) in t.coeffs {
- expr.add_term(co, gurobi_vars[v.index()]);
- }
- let name = format!("c{c_id}");
- let constr = match c.sense {
- Sense::Le => {
- grb_model.add_constr(&name, c!(expr <= adjusted_rhs)).map_err(map_grb_err)?
- }
- Sense::Ge => {
- grb_model.add_constr(&name, c!(expr >= adjusted_rhs)).map_err(map_grb_err)?
- }
- Sense::Eq => {
- grb_model.add_constr(&name, c!(expr == adjusted_rhs)).map_err(map_grb_err)?
+ // Linear fast path per constraint. We fall back to the general lowering only
+ // for those that need it. Keep linear constraints tracked so duals/Pi can
+ // still be reported in the LP/MILP case.
+ let con_lin_terms: Vec> =
+ constraints.iter().map(|c| extract_linear(arena_ref, c.lhs)).collect();
+
+ let mut gurobi_constrs: Vec> = Vec::with_capacity(constraints.len());
+ let mut aux_counter = 0_u32;
+
+ for (c_id, (c, lin)) in constraints.iter().zip(con_lin_terms).enumerate() {
+ if let Some(t) = lin {
+ let adjusted_rhs = c.rhs - t.constant;
+ let mut expr = LinExpr::new();
+ for (v, co) in t.coeffs {
+ expr.add_term(co, gurobi_vars[v.index()]);
}
- };
- gurobi_constrs.push(constr);
+ let name = format!("c{c_id}");
+ let constr = match c.sense {
+ Sense::Le => {
+ grb_model.add_constr(&name, c!(expr <= adjusted_rhs)).map_err(map_grb_err)?
+ }
+ Sense::Ge => {
+ grb_model.add_constr(&name, c!(expr >= adjusted_rhs)).map_err(map_grb_err)?
+ }
+ Sense::Eq => {
+ grb_model.add_constr(&name, c!(expr == adjusted_rhs)).map_err(map_grb_err)?
+ }
+ };
+ gurobi_constrs.push(Some(constr));
+ } else {
+ add_nonlinear_constraint(
+ &arena,
+ c.lhs,
+ c.sense,
+ c.rhs,
+ c_id,
+ &mut grb_model,
+ &gurobi_vars,
+ &mut aux_counter,
+ )?;
+ gurobi_constrs.push(None);
+ }
}
- grb_model
- .set_attr(
- attr::ModelSense,
- match objective.sense {
- ObjectiveSense::Minimize => 1,
- ObjectiveSense::Maximize => -1,
- },
- )
- .map_err(map_grb_err)?;
+ let obj_constant = set_objective(
+ &arena,
+ objective.expr,
+ objective.sense,
+ &mut grb_model,
+ &gurobi_vars,
+ &mut aux_counter,
+ )?;
apply_options(&mut grb_model, opts).map_err(map_grb_err)?;
+ if nonlinear_kind {
+ // Gurobi requires NonConvex=2 for general nonlinear constraints and
+ // bilinear non-convex objectives. Skip if the user already set it.
+ let current = grb_model.get_param(grb::param::NonConvex).map_err(map_grb_err)?;
+ if current < 2 {
+ grb_model.set_param(grb::param::NonConvex, 2).map_err(map_grb_err)?;
+ }
+ }
let started = Instant::now();
grb_model.optimize().map_err(map_grb_err)?;
@@ -117,9 +130,10 @@ pub fn solve(model: &Model, opts: &GurobiOptions) -> Result Result Result<(), SolverError> {
+ let mut ctx = LoweringCtx { model: grb_model, gurobi_vars, aux_counter: *aux_counter };
+ let lowered = lower(arena, lhs, &mut ctx)?;
+ *aux_counter = ctx.aux_counter;
+ let name = format!("c{c_id}");
+ match lowered {
+ LoweredExpr::Linear(e) => {
+ match sense {
+ Sense::Le => grb_model.add_constr(&name, c!(e <= rhs)),
+ Sense::Ge => grb_model.add_constr(&name, c!(e >= rhs)),
+ Sense::Eq => grb_model.add_constr(&name, c!(e == rhs)),
+ }
+ .map_err(map_grb_err)?;
+ }
+ LoweredExpr::Quadratic(e) => {
+ match sense {
+ Sense::Le => grb_model.add_qconstr(&name, c!(e <= rhs)),
+ Sense::Ge => grb_model.add_qconstr(&name, c!(e >= rhs)),
+ Sense::Eq => grb_model.add_qconstr(&name, c!(e == rhs)),
+ }
+ .map_err(map_grb_err)?;
+ }
+ LoweredExpr::Var(v) => {
+ match sense {
+ Sense::Le => grb_model.add_constr(&name, c!(v <= rhs)),
+ Sense::Ge => grb_model.add_constr(&name, c!(v >= rhs)),
+ Sense::Eq => grb_model.add_constr(&name, c!(v == rhs)),
+ }
+ .map_err(map_grb_err)?;
+ }
+ }
+ Ok(())
+}
+
+fn set_objective(
+ arena: &ExprArena,
+ obj_expr: ExprId,
+ sense: ObjectiveSense,
+ grb_model: &mut grb::Model,
+ gurobi_vars: &[grb::Var],
+ aux_counter: &mut u32,
+) -> Result {
+ let grb_sense = match sense {
+ ObjectiveSense::Minimize => ModelSense::Minimize,
+ ObjectiveSense::Maximize => ModelSense::Maximize,
+ };
+ if let Some(t) = extract_linear(arena, obj_expr) {
+ let mut e = LinExpr::new();
+ for (v, c) in t.coeffs {
+ e.add_term(c, gurobi_vars[v.index()]);
+ }
+ // Gurobi's set_objective absorbs LinExpr offsets into ObjCon, so we do
+ // not need to track the constant separately.
+ e.add_constant(t.constant);
+ grb_model.set_objective(e, grb_sense).map_err(map_grb_err)?;
+ return Ok(0.0);
+ }
+ let mut ctx = LoweringCtx { model: grb_model, gurobi_vars, aux_counter: *aux_counter };
+ let lowered = lower(arena, obj_expr, &mut ctx)?;
+ *aux_counter = ctx.aux_counter;
+ grb_model.set_objective(lowered.into_expr_for_objective(), grb_sense).map_err(map_grb_err)?;
+ Ok(0.0)
+}
+
fn collect_solution(
status: &SolverStatus,
+ kind: ModelKind,
model: &grb::Model,
vars: &[grb::Var],
- constrs: &[grb::Constr],
+ constrs: &[Option],
) -> (FxHashMap, FxHashMap, FxHashMap) {
// `has_solution` only flags Optimal/Feasible, but Gurobi often holds an
// incumbent on TimeLimit/IterationLimit/NodeLimit too.
@@ -148,29 +237,41 @@ fn collect_solution(
}
let primal_vals = model.get_obj_attr_batch(attr::X, vars.iter().copied()).ok();
- let rc_vals = model.get_obj_attr_batch(attr::RC, vars.iter().copied()).ok();
- let pi_vals = model.get_obj_attr_batch(attr::Pi, constrs.iter().copied()).ok();
-
- let to_var_map = |vals: Option>| -> FxHashMap {
- vals.map(|v| {
+ let primal = primal_vals
+ .map(|v| {
v.into_iter()
.enumerate()
.map(|(i, val)| (VarId(u32::try_from(i).unwrap()), val))
.collect()
})
- .unwrap_or_default()
- };
- let to_constr_map = |vals: Option>| -> FxHashMap {
- vals.map(|v| {
+ .unwrap_or_default();
+
+ // Skip retrieval of duals and reduced costs for any model
+ // class where Gurobi will either return zeros or refuse the attribute.
+ if !matches!(kind, ModelKind::LP) {
+ return (primal, FxHashMap::default(), FxHashMap::default());
+ }
+
+ let rc_vals = model.get_obj_attr_batch(attr::RC, vars.iter().copied()).ok();
+ let reduced_costs = rc_vals
+ .map(|v| {
v.into_iter()
.enumerate()
- .map(|(i, val)| (ConstraintId(u32::try_from(i).unwrap()), val))
+ .map(|(i, val)| (VarId(u32::try_from(i).unwrap()), val))
.collect()
})
- .unwrap_or_default()
- };
+ .unwrap_or_default();
+
+ let mut dual = FxHashMap::default();
+ for (i, c) in constrs.iter().enumerate() {
+ if let Some(c) = c {
+ if let Ok(pi) = model.get_obj_attr(attr::Pi, c) {
+ dual.insert(ConstraintId(u32::try_from(i).unwrap()), pi);
+ }
+ }
+ }
- (to_var_map(primal_vals), to_var_map(rc_vals), to_constr_map(pi_vals))
+ (primal, reduced_costs, dual)
}
fn map_status(model: &grb::Model) -> Result {
diff --git a/crates/oximo-gurobi/tests/nonlinear.rs b/crates/oximo-gurobi/tests/nonlinear.rs
new file mode 100644
index 0000000..94beeac
--- /dev/null
+++ b/crates/oximo-gurobi/tests/nonlinear.rs
@@ -0,0 +1,144 @@
+//! Live Gurobi tests for QP, NLP, and MINLP models.
+
+use oximo_core::prelude::*;
+use oximo_gurobi::{Gurobi, GurobiOptions};
+use oximo_solver::{Solver, SolverStatus};
+
+fn close(a: f64, b: f64, tol: f64) -> bool {
+ (a - b).abs() < tol
+}
+
+fn assert_solved(r: &oximo_solver::SolverResult) {
+ assert!(
+ matches!(r.status, SolverStatus::Optimal | SolverStatus::Feasible),
+ "status = {:?}",
+ r.status
+ );
+}
+
+#[test]
+fn qp_min_sum_of_squares() {
+ // min x^2 + y^2 s.t. x + y >= 1.
+ // Optimum at x = y = 0.5, objective = 0.5.
+ let m = Model::new("qp");
+ let x = m.var("x").lb(-10.0).ub(10.0).build();
+ let y = m.var("y").lb(-10.0).ub(10.0).build();
+ m.constraint("c", (x + y).ge(1.0));
+ m.minimize(x.powi(2) + y.powi(2));
+
+ let r = Gurobi.solve(&m, &GurobiOptions::default()).expect("solve");
+ assert!(matches!(r.status, SolverStatus::Optimal | SolverStatus::Feasible));
+ let obj = r.objective.expect("obj");
+ assert!(close(obj, 0.5, 1e-4), "obj = {obj}");
+}
+
+#[test]
+fn nlp_with_sin_objective() {
+ // min (x - 1)^2 + 0.1 * sin(x)^2 over x in [-3, 3].
+ // Local minimum near x = 1, objective near 0.
+ let m = Model::new("nlp_sin");
+ let x = m.var("x").lb(-3.0).ub(3.0).initial(0.5).build();
+ let one = Expr::constant(x.arena, 1.0);
+ m.minimize((x - one).powi(2) + Expr::constant(x.arena, 0.1) * x.sin().powi(2));
+
+ let r = Gurobi.solve(&m, &GurobiOptions::default()).expect("solve");
+ assert!(matches!(r.status, SolverStatus::Optimal | SolverStatus::Feasible));
+ let primal_x = r.primal.get(&VarId(0)).copied().expect("primal");
+ assert!(close(primal_x, 1.0, 0.1), "x = {primal_x}");
+}
+
+#[test]
+fn minlp_binary_with_log() {
+ // Binary b, continuous x in [0.1, 10]. Min (x - 1)^2 + b * log(1 + x).
+ // Optimal: b = 0, x = 1, objective = 0.
+ let m = Model::new("minlp_log");
+ let b = m.var("b").binary().build();
+ let x = m.var("x").lb(0.1).ub(10.0).initial(0.5).build();
+ let one = Expr::constant(x.arena, 1.0);
+ m.minimize((x - one).powi(2) + b * (one + x).log());
+
+ let r = Gurobi.solve(&m, &GurobiOptions::default()).expect("solve");
+ assert!(matches!(r.status, SolverStatus::Optimal | SolverStatus::Feasible));
+ let obj = r.objective.expect("obj");
+ assert!(close(obj, 0.0, 1e-3), "obj = {obj}");
+}
+
+// Division lowering
+
+#[test]
+fn div_by_linear_denominator() {
+ // x / (y + z) == 3, with x = 12 and z = 1 fixed -> y + 1 = 4 -> y = 3.
+ let m = Model::new("div_linear");
+ let x = m.var("x").fix(12.0).build();
+ let y = m.var("y").lb(0.1).ub(100.0).build();
+ let z = m.var("z").fix(1.0).build();
+ m.constraint("c", (x / (y + z)).eq(3.0));
+ m.minimize(y);
+
+ let sol = Gurobi.solve(&m, &GurobiOptions::default()).expect("solve");
+ assert_solved(&sol);
+ let yv = sol.primal.get(&VarId(1)).copied().expect("primal y");
+ assert!(close(yv, 3.0, 1e-4), "y = {yv}");
+}
+
+#[test]
+fn div_by_negative_denominator() {
+ // x / d == -3, with x = 12 fixed and d in [-100, -0.1] -> d = -4.
+ // A `pow(den, -1)` lowering could not represent a negative denominator,
+ // the bilinear `d * recip == 1` pin can.
+ let m = Model::new("div_negative");
+ let x = m.var("x").fix(12.0).build();
+ let d = m.var("d").lb(-100.0).ub(-0.1).build();
+ m.constraint("c", (x / d).eq(-3.0));
+ m.minimize(d);
+
+ let r = Gurobi.solve(&m, &GurobiOptions::default()).expect("solve");
+ assert_solved(&r);
+ let dv = r.primal.get(&VarId(1)).copied().expect("primal d");
+ assert!(close(dv, -4.0, 1e-4), "d = {dv}");
+}
+
+#[test]
+fn div_constant_numerator() {
+ // 6 / d == 2 -> d = 3. Exercises the constant-numerator fold (`6 * recip`,
+ // which stays linear in `recip` rather than materializing a product).
+ let m = Model::new("div_const_num");
+ let d = m.var("d").lb(0.1).ub(100.0).build();
+ m.constraint("c", (6.0 / d).eq(2.0));
+ m.minimize(d);
+
+ let r = Gurobi.solve(&m, &GurobiOptions::default()).expect("solve");
+ assert_solved(&r);
+ let dv = r.primal.get(&VarId(0)).copied().expect("primal d");
+ assert!(close(dv, 3.0, 1e-4), "d = {dv}");
+}
+
+#[test]
+fn div_by_quadratic_denominator() {
+ // x / (x * y) == 0.5 reduces to 1 / y == 0.5 -> y = 2 for any nonzero x.
+ // The quadratic denominator is first materialized into an aux variable so
+ // the reciprocal pin stays bilinear rather than cubic.
+ let m = Model::new("div_quadratic");
+ let x = m.var("x").lb(1.0).ub(10.0).build();
+ let y = m.var("y").lb(0.1).ub(10.0).build();
+ m.constraint("c", (x / (x * y)).eq(0.5));
+ m.minimize(x);
+
+ let r = Gurobi.solve(&m, &GurobiOptions::default()).expect("solve");
+ assert_solved(&r);
+ let yv = r.primal.get(&VarId(1)).copied().expect("primal y");
+ assert!(close(yv, 2.0, 1e-3), "y = {yv}");
+}
+
+#[test]
+fn div_by_zero_constant_errors() {
+ // A literal zero denominator survives construction as a `Div` node (only
+ // nonzero constants are folded into the linear path), so lowering must
+ // reject it rather than emit an infeasible `0 * recip == 1`.
+ let m = Model::new("div_zero");
+ let x = m.var("x").lb(0.0).ub(10.0).build();
+ m.minimize(x / 0.0);
+
+ let err = Gurobi.solve(&m, &GurobiOptions::default()).expect_err("expected error");
+ assert!(err.to_string().contains("division by zero"), "err = {err}");
+}
diff --git a/crates/oximo/README.md b/crates/oximo/README.md
index 173bd12..67b623d 100644
--- a/crates/oximo/README.md
+++ b/crates/oximo/README.md
@@ -11,9 +11,7 @@
-oximo is a Rust algebraic modeling library for mathematical optimization. Build LP and MILP models with a clean builder API, then solve them with bundled or commercial solvers.
-
-> Support for nonlinear programming (NLP) and mixed-integer nonlinear programming (MINLP) is planned.
+oximo is a Rust algebraic modeling library for mathematical optimization. Build LP, MILP, QP/MIQP, NLP, and MINLP models with a clean builder API, then solve them with bundled or commercial solvers.
```rust,no_run
use oximo::prelude::*;
@@ -37,12 +35,12 @@ println!("y = {:?}", result.value_of(y)); // 4.0
## Features
-| Feature | What it adds | Default |
-|----------|---------------------------------------------------|---------|
-| `highs` | HiGHS LP/MILP solver (bundled, no install) | yes |
-| `io` | MPS and LP file writers | yes |
-| `gurobi` | Gurobi LP/MILP solver (requires licensed install) | no |
-| `gams` | GAMS solver bridge (requires GAMS on PATH) | no |
+| Feature | What it adds | Default |
+|----------|-----------------------------------------------------------------------|---------|
+| `highs` | HiGHS - LP/MILP solver (bundled, no install) | yes |
+| `io` | MPS and LP file writers | yes |
+| `gurobi` | Gurobi - LP/MILP/QP/MIQP/NLP/MINLP solver (requires licensed install) | no |
+| `gams` | GAMS bridge - LP/MILP/QP/MIQP/NLP/MINLP depending on solver | no |
```toml
[dependencies]
@@ -170,6 +168,23 @@ m.add_constraints_over("supply", &plants, |p: String| {
m.add_constraints_over("c", &set, |k: IndexKey| x[&k].le(1.0));
```
+### Nonlinear expressions
+
+`Pow`, `Sin`, `Cos`, `Exp`, `Log`, and bilinear products are first-class. The
+model's kind (`LP`/`MILP`/`QP`/`MIQP`/`NLP`/`MINLP`) is inferred from the
+expressions.
+
+```rust,ignore
+// Rosenbrock NLP
+m.minimize((1.0 - x).powi(2) + 100.0 * (y - x.powi(2)).powi(2));
+
+// Quadratic constraint
+m.constraint("disk", (x.powi(2) + y.powi(2)).le(1.0));
+
+// Transcendental utility (MINLP when any variable is integer/binary)
+m.maximize(sum_over(&items, |i: usize| u[i] * (1.0 + w[i] * x[i]).log()));
+```
+
## Solving
All backends implement the `Solver` trait:
diff --git a/crates/oximo/examples/process_selection.rs b/crates/oximo/examples/process_selection.rs
new file mode 100644
index 0000000..409adc9
--- /dev/null
+++ b/crates/oximo/examples/process_selection.rs
@@ -0,0 +1,149 @@
+//! Structural optimization of a process flowsheet
+//! This is the GAMS model library problem `PROCSEL`
+//! (SEQ=116), from:
+//!
+//! Kocis & Grossmann (1987), "Relaxation Strategy for the Structural
+//! Optimization of Process Flow Sheets", Ind. Eng. Chem. Res. 26(9),
+//! 1869-1880; also Morari & Grossmann (eds.), "Chemical Engineering
+//! Optimization Models with GAMS" (1991).
+//!
+//! Chemical C is produced from B in unit 1. B is either purchased on the
+//! external market (`bp`) or produced from raw material A through one of two
+//! competing units (2 or 3). Binaries `y1,y2,y3` switch the three units on/off;
+//! the goal is to maximise annual profit.
+//!
+//! A2 +-----+ B2 BP
+//! +----->| 2 |----->+ |
+//! A | +-----+ | | B1 +-----+ C1
+//! ---->| +----+------->| 1 |-------->
+//! | +-----+ | +-----+
+//! +----->| 3 |----->+
+//! A3 +-----+ B3
+//!
+//! The input-output relations of units 2 and 3 are the (convexified) nonlinear
+//! laws `exp(b2) - 1 = a2` and `exp(b3/1.2) - 1 = a3`, which make this an MINLP.
+//!
+//! Run with one nonlinear backend (enable exactly one, never both at once):
+//! cargo run --example process_selection --features gams
+//! cargo run --example process_selection --features gurobi
+
+#[cfg(any(feature = "gams", feature = "gurobi"))]
+use oximo::prelude::*;
+
+/// Build the flowsheet model, solve it with `solver`, and print the solution.
+/// Generic over the [`Solver`] backend and its options type.
+#[cfg(any(feature = "gams", feature = "gurobi"))]
+fn solve_and_report(
+ label: &str,
+ mut solver: S,
+ opts: &S::Options,
+) -> Result<(), Box> {
+ let m = Model::new("procsel");
+
+ // Positive variables (consumptions, capacities, purchases).
+ let a2 = m.var("a2").lb(0.0).build();
+ let a3 = m.var("a3").lb(0.0).build();
+ let b2 = m.var("b2").lb(0.0).build();
+ let b3 = m.var("b3").lb(0.0).build();
+ let bp = m.var("bp").lb(0.0).build();
+ let b1 = m.var("b1").lb(0.0).build();
+ let c1 = m.var("c1").lb(0.0).ub(1.0).build();
+
+ // Binaries: existence of each process unit.
+ let y1 = m.var("y1").binary().build();
+ let y2 = m.var("y2").binary().build();
+ let y3 = m.var("y3").binary().build();
+
+ m.constraint("inout1", c1.eq(0.9 * b1));
+ m.constraint("inout2", (b2.exp() - 1.0).eq(a2));
+ m.constraint("inout3", ((b3 / 1.2).exp() - 1.0).eq(a3));
+ m.constraint("mbalb", b1.eq(b2 + b3 + bp));
+ m.constraint("log1", c1.le(2.0 * y1));
+ m.constraint("log2", b2.le(4.0 * y2));
+ m.constraint("log3", b3.le(5.0 * y3));
+
+ // profit = sales - fixed investment - operating cost - purchases
+ m.maximize(11.0 * c1 - 3.5 * y1 - y2 - 1.5 * y3 - b2 - 1.2 * b3 - 1.8 * (a2 + a3) - 7.0 * bp);
+
+ let result = solver.solve(&m, opts)?;
+
+ if let Some(obj) = result.objective {
+ println!("--- {label} ---");
+ println!("status = {:?}", result.status);
+ println!("profit = {obj:.4} M$/yr");
+
+ println!("units selected:");
+ for (name, y) in [("process 1", y1), ("process 2", y2), ("process 3", y3)] {
+ println!(
+ " {name}: {}",
+ if (result.value_of(y).unwrap_or(0.0) - 1.0).abs() < f64::EPSILON {
+ "ON"
+ } else {
+ "OFF"
+ }
+ );
+ }
+
+ println!("flows:");
+ println!(" c1 (C produced) = {:.4}", result.value_of(c1).unwrap_or(0.0));
+ println!(" b1 (B into 1) = {:.4}", result.value_of(b1).unwrap_or(0.0));
+ println!(
+ " b2, b3 (B from A) = {:.4}, {:.4}",
+ result.value_of(b2).unwrap_or(0.0),
+ result.value_of(b3).unwrap_or(0.0)
+ );
+ println!(" bp (B purchased) = {:.4}", result.value_of(bp).unwrap_or(0.0));
+ println!(
+ " a2, a3 (A used) = {:.4}, {:.4}",
+ result.value_of(a2).unwrap_or(0.0),
+ result.value_of(a3).unwrap_or(0.0)
+ );
+ } else {
+ println!("--- {label} ---");
+ println!("status = {:?}", result.status);
+ println!("no objective value");
+ }
+
+ Ok(())
+}
+
+#[cfg(feature = "gams")]
+fn run_gams() -> Result<(), Box> {
+ use oximo::gams::{GamsBaronOptions, GamsSolverConfig};
+ use oximo::solvers::Gams;
+ use std::time::Duration;
+
+ let opts = GamsOptions::default()
+ .time_limit(Duration::from_secs(120))
+ .solver(GamsSolverConfig::Baron(GamsBaronOptions {
+ eps_r: Some(1e-4),
+ ..Default::default()
+ }))
+ .verbose(true);
+ solve_and_report("GAMS + BARON", Gams::new(), &opts)
+}
+
+#[cfg(feature = "gurobi")]
+fn run_gurobi() -> Result<(), Box> {
+ use oximo::solvers::Gurobi;
+ use std::time::Duration;
+
+ let opts = GurobiOptions::default().time_limit(Duration::from_secs(120));
+ solve_and_report("Gurobi", Gurobi, &opts)
+}
+
+#[cfg(any(feature = "gams", feature = "gurobi"))]
+fn main() -> Result<(), Box> {
+ #[cfg(feature = "gams")]
+ run_gams()?;
+ #[cfg(feature = "gurobi")]
+ run_gurobi()?;
+ Ok(())
+}
+
+#[cfg(not(any(feature = "gams", feature = "gurobi")))]
+fn main() {
+ println!("Enable a nonlinear-capable backend feature:");
+ println!(" cargo run --example process_selection --features gams");
+ println!(" cargo run --example process_selection --features gurobi");
+}