Support nonlinear constraints with Gurobi 12 (#95)

Closes #93

* Implement support for NLExpr objects in add_constrs
* Add a page to the users guide covering how to add nonlinear constraints
to a model using add_constrs.
* Show how to create nonlinear constraints not of the form y = f(x)
by introducing temporary variables.

Author: simonbowly

docs/source/index.rst

@@ -78,6 +78,7 @@ If you encounter issues with Gurobi or ``gurobipy`` please contact
 
    performance
    naming
+   nonlinear
    advanced
    typing
 

docs/source/nonlinear.rst

@@ -0,0 +1,119 @@
+Adding Nonlinear Constraints
+============================
+
+Gurobi 12 supports adding nonlinear constraints, using the ``NLExpr`` object to
+capture nonlinear expressions. ``gurobipy-pandas`` supports adding a ``Series``
+of nonlinear constraints to a model via ``add_constrs``. Note that ``gurobipy``
+only supports nonlinear constraints of the form :math:`y = f(\bar{x})` and
+``gurobipy-pandas`` applies the same restriction.
+
+.. note::
+
+   To add nonlinear constraints, you must have at least ``gurobipy`` version
+   12.0.0 installed.
+
+Nonlinear Equality Constraints
+------------------------------
+
+This example builds the constraint set :math:`y_i = \log(\frac{1}{x_i})`, for
+each :math:`i` in the index.
+
+.. doctest:: [nonlinear]
+
+    >>> import pandas as pd
+    >>> import gurobipy as gp
+    >>> from gurobipy import GRB
+    >>> from gurobipy import nlfunc
+    >>> import gurobipy_pandas as gppd
+    >>> gppd.set_interactive()
+
+    >>> model = gp.Model()
+    >>> index = pd.RangeIndex(5)
+    >>> x = gppd.add_vars(model, index, lb=1.0, name="x")
+    >>> y = gppd.add_vars(model, index, lb=-GRB.INFINITY, name="y")
+
+You can create nonlinear expressions using standard Python operators. A
+nonlinear expression is any expression which is not a polynomial of degree at
+most 2. For example, dividing a constant by a series of ``Var`` objects produces
+a series of nonlinear expressions.
+
+.. doctest:: [nonlinear]
+
+   >>> 1 / x
+   0    1.0 / x[0]
+   1    1.0 / x[1]
+   2    1.0 / x[2]
+   3    1.0 / x[3]
+   4    1.0 / x[4]
+   Name: x, dtype: object
+
+The ``nlfunc`` module of gurobipy is used to create nonlinear expressions
+involving mathematical functions. You can use ``apply`` to construct a series
+representing the logarithm of the above result.
+
+.. doctest:: [nonlinear]
+
+   >>> (1 / x).apply(nlfunc.log)
+   0    log(1.0 / x[0])
+   1    log(1.0 / x[1])
+   2    log(1.0 / x[2])
+   3    log(1.0 / x[3])
+   4    log(1.0 / x[4])
+   Name: x, dtype: object
+
+This series of expressions can be added as constraints using ``add_constrs``.
+Note that you can only provide nonlinear constraints in the form :math:`y =
+f(\bar{x})` where :math:`f` may be a univariate or multivariate compound
+function. Therefore the ``lhs`` argument must be a series of ``Var`` objects,
+while the ``sense`` argument must be ``GRB.EQUAL``.
+
+.. doctest:: [nonlinear]
+
+   >>> gppd.add_constrs(model, y, GRB.EQUAL, (1 / x).apply(nlfunc.log), name="log_x")
+   0    <gurobi.GenConstr 0>
+   1    <gurobi.GenConstr 1>
+   2    <gurobi.GenConstr 2>
+   3    <gurobi.GenConstr 3>
+   4    <gurobi.GenConstr 4>
+   Name: log_x, dtype: object
+
+Nonlinear Inequality Constraints
+--------------------------------
+
+As noted above, nonlinear constraints can only be provided in the form :math:`y=
+f(\bar{x})`. To formulate inequality constraints, you must introduce bounded
+intermediate variables. The following example formulates the set of constraints
+:math:`\log(x_i^2 + 1) \le 2` by introducing intermediate variables :math:`z_i`
+with no lower bound and an upper bound of :math:`2.0`.
+
+.. doctest:: [nonlinear]
+
+   >>> z = gppd.add_vars(model, index, lb=-GRB.INFINITY, ub=2.0, name="z")
+   >>> constrs = gppd.add_constrs(model, z, GRB.EQUAL, (x**2 + 1).apply(nlfunc.log))
+
+To see the effect of this constraint, you can set a maximization objective
+:math:`\sum_{i=0}^{4} x_i` and solve the resulting model. You will find that the
+original variables :math:`x_i` are maximized to :math:`\sqrt{e^2 - 1}` in
+the optimal solution. The intermediate variables :math:`z_i` are at their upper
+bounds, indicating that the constraint is satisfied with equality.
+
+.. doctest:: [nonlinear]
+
+   >>> model.setObjective(x.sum(), sense=GRB.MAXIMIZE)
+   >>> model.optimize()  # doctest: +ELLIPSIS
+   Gurobi Optimizer ...
+   Optimal solution found ...
+   >>> x.gppd.X.round(3)
+   0    2.528
+   1    2.528
+   2    2.528
+   3    2.528
+   4    2.528
+   Name: x, dtype: float64
+   >>> z.gppd.X.round(3)
+   0    2.0
+   1    2.0
+   2    2.0
+   3    2.0
+   4    2.0
+   Name: z, dtype: float64

src/gurobipy_pandas/constraints.py

@@ -13,6 +13,7 @@
 from gurobipy_pandas.util import create_names, gppd_global_options
 
 CONSTRAINT_SENSES = frozenset([GRB.LESS_EQUAL, GRB.EQUAL, GRB.GREATER_EQUAL])
+GUROBIPY_MAJOR_VERSION, *_ = gp.gurobi.version()
 
 
 def add_constrs_from_dataframe(
@@ -138,6 +139,15 @@ def _add_constr(model, lhs, sense, rhs, name):
         name = ""
     if not isinstance(sense, str) or sense[0] not in CONSTRAINT_SENSES:
         raise ValueError(f"'{sense}' is not a valid constraint sense")
+    if GUROBIPY_MAJOR_VERSION >= 12:
+        # Check for nonlinear constraints; accept only y = f(x)
+        if isinstance(lhs, gp.NLExpr):
+            raise ValueError("Nonlinear constraints must be in the form y = f(x)")
+        if isinstance(rhs, gp.NLExpr):
+            if isinstance(lhs, gp.Var) and sense == GRB.EQUAL:
+                return model.addGenConstrNL(lhs, rhs, name=name)
+            else:
+                raise ValueError("Nonlinear constraints must be in the form y = f(x)")
     if isinstance(lhs, gp.QuadExpr) or isinstance(rhs, gp.QuadExpr):
         return model.addQConstr(lhs, sense, rhs, name=name)
     return model.addLConstr(lhs, sense, rhs, name=name)

tests/test_api.py

@@ -3,13 +3,19 @@
 tests of data types, errors, etc, are done on the lower-level functions.
 """
 
+import math
+import unittest
+
+import gurobipy as gp
 import pandas as pd
 from gurobipy import GRB
 from pandas.testing import assert_index_equal, assert_series_equal
 
 import gurobipy_pandas as gppd
 from tests.utils import GurobiModelTestCase
 
+GUROBIPY_MAJOR_VERSION, *_ = gp.gurobi.version()
+
 
 class TestAddVars(GurobiModelTestCase):
     def test_from_dataframe(self):
@@ -269,6 +275,136 @@ def test_sense_series(self):
             self.assertEqual(constr.RHS, -1.0)
 
 
+@unittest.skipIf(
+    GUROBIPY_MAJOR_VERSION < 12,
+    "Nonlinear constraints are only supported for Gurobi 12 and later",
+)
+class TestNonlinear(GurobiModelTestCase):
+    def assert_approx_equal(self, value, expected, tolerance=1e-6):
+        difference = abs(value - expected)
+        self.assertLessEqual(difference, tolerance)
+
+    def test_log(self):
+        # max sum y_i
+        # s.t. y_i = log(x_i)
+        #      1.0 <= x <= 2.0
+        from gurobipy import nlfunc
+
+        index = pd.RangeIndex(3)
+        x = gppd.add_vars(self.model, index, lb=1.0, ub=2.0, name="x")
+        y = gppd.add_vars(self.model, index, obj=1.0, name="y")
+        self.model.ModelSense = GRB.MAXIMIZE
+
+        gppd.add_constrs(self.model, y, GRB.EQUAL, x.apply(nlfunc.log), name="log_x")
+
+        self.model.optimize()
+        self.assert_approx_equal(self.model.ObjVal, 3 * math.log(2.0))
+
+    def test_inequality(self):
+        # max  sum x_i
+        # s.t. log2(x_i^2 + 1) <= 2.0
+        #      0.0 <= x <= 1.0
+        #
+        # Formulated as
+        #
+        # max  sum x_i
+        # s.t. log2(x_i^2 + 1) == z_i
+        #      0.0 <= x <= 1.0
+        #      -GRB.INFINITY <= z_i <= 2
+        from gurobipy import nlfunc
+
+        index = pd.RangeIndex(3)
+        x = gppd.add_vars(self.model, index, name="x")
+        z = gppd.add_vars(self.model, index, lb=-GRB.INFINITY, ub=2.0, name="z")
+        self.model.setObjective(x.sum(), sense=GRB.MAXIMIZE)
+
+        gppd.add_constrs(self.model, z, GRB.EQUAL, (x**2 + 1).apply(nlfunc.log))
+
+        self.model.optimize()
+        self.model.write("model.lp")
+        self.model.write("model.sol")
+
+        x_sol = x.gppd.X
+        z_sol = z.gppd.X
+        for i in range(3):
+            self.assert_approx_equal(x_sol[i], math.sqrt(math.exp(2.0) - 1))
+            self.assert_approx_equal(z_sol[i], 2.0)
+
+    def test_wrong_usage(self):
+        index = pd.RangeIndex(3)
+        x = gppd.add_vars(self.model, index, name="x")
+        y = gppd.add_vars(self.model, index, name="y")
+
+        with self.assertRaisesRegex(
+            gp.GurobiError, "Objective must be linear or quadratic"
+        ):
+            self.model.setObjective((x / y).sum())
+
+        with self.assertRaisesRegex(
+            ValueError, "Nonlinear constraints must be in the form"
+        ):
+            gppd.add_constrs(self.model, y, GRB.LESS_EQUAL, x**4)
+
+        with self.assertRaisesRegex(
+            ValueError, "Nonlinear constraints must be in the form"
+        ):
+            gppd.add_constrs(self.model, y + x**4, GRB.EQUAL, 1)
+
+        with self.assertRaisesRegex(
+            ValueError, "Nonlinear constraints must be in the form"
+        ):
+            gppd.add_constrs(self.model, y**4, GRB.EQUAL, x)
+
+        with self.assertRaisesRegex(
+            ValueError, "Nonlinear constraints must be in the form"
+        ):
+            x.to_frame().gppd.add_constrs(self.model, "x**3 == 1", name="bad")
+
+    def test_eval(self):
+        index = pd.RangeIndex(3)
+        df = (
+            gppd.add_vars(self.model, index, name="x")
+            .to_frame()
+            .gppd.add_vars(self.model, name="y")
+            .gppd.add_constrs(self.model, "y == x**3", name="nlconstr")
+        )
+
+        self.model.update()
+        for row in df.itertuples(index=False):
+            self.assert_nlconstr_equal(
+                row.nlconstr,
+                row.y,
+                [
+                    (GRB.OPCODE_POW, -1.0, -1),
+                    (GRB.OPCODE_VARIABLE, row.x, 0),
+                    (GRB.OPCODE_CONSTANT, 3.0, 0),
+                ],
+            )
+
+    def test_frame(self):
+        from gurobipy import nlfunc
+
+        index = pd.RangeIndex(3)
+        df = (
+            gppd.add_vars(self.model, index, name="x")
+            .to_frame()
+            .gppd.add_vars(self.model, name="y")
+            .assign(exp_x=lambda df: df["x"].apply(nlfunc.exp))
+            .gppd.add_constrs(self.model, "y", GRB.EQUAL, "exp_x", name="nlconstr")
+        )
+
+        self.model.update()
+        for row in df.itertuples(index=False):
+            self.assert_nlconstr_equal(
+                row.nlconstr,
+                row.y,
+                [
+                    (GRB.OPCODE_EXP, -1.0, -1),
+                    (GRB.OPCODE_VARIABLE, row.x, 0),
+                ],
+            )
+
+
 class TestDataValidation(GurobiModelTestCase):
     # Test that we throw some more informative errors, instead of obscure
     # ones from the underlying gurobipy library