Add Gurobi as optional solver

For discussion that concluded with this commit,
#98
#60

Gurobi will not automatically become the default solver due to its
restrictive license. However, users can manually select it by

    import polytope as pc
    pc.solvers.default_solver = 'gurobi'

Then, `polytope` can be used without modification.

If Gurobi returns INF_OR_UNBD, the optimization is run again with
DualReductions=0 as recommended at
https://www.gurobi.com/documentation/current/refman/optimization_status_codes.html
and a notification about this is included under the `"message"` key in
the result returned by `polytope.solvers.lpsolve`.

Tests that require Gurobi are skipped and a warning is shown if
gurobipy is not installed.


Co-authored-by: Richard Oberdieck <[email protected]>

.github/workflows/main.yml

@@ -81,6 +81,7 @@ jobs:
                 set -o posix
                 echo 'Exported environment variables:'
                 export -p
+                cd tests/
                 pytest \
                     -v \
                     --continue-on-collection-errors \
@@ -90,6 +91,9 @@ jobs:
                 export CVXOPT_BUILD_GLPK=1
                 pip install cvxopt
                 python -c 'import cvxopt.glpk'
+            - name: Install optional solvers with restrictive licenses
+              run: |
+                pip install gurobipy
             - name: Run all tests, using `cvxopt.glpk`
               run: |
                 set -o posix

polytope/solvers.py

@@ -10,6 +10,12 @@
 You can change this default at runtime by setting the variable
 `default_solver` in the module `solvers`.
 
+Nonfree (i.e., having restrictive licenses) solvers are also
+supported but require extra packages:
+
+* Gurobi, https://pypi.org/project/gurobipy/
+* MOSEK, https://pypi.org/project/Mosek/
+
 For example:
 
 ```python
@@ -49,6 +55,12 @@
 except ImportError:
     logger.info('MOSEK solver not found.')
 
+try:
+    import gurobipy as gurobi
+    installed_solvers.add('gurobi')
+except ImportError:
+    logger.info('GUROBI solver not found')
+
 
 # choose default from installed choices
 if 'glpk' in installed_solvers:
@@ -61,7 +73,6 @@
         "`installed_solvers` wasn't empty above?")
 
 
-
 def lpsolve(c, G, h, solver=None):
     """Try to solve linear program with given or default solver.
 
@@ -87,6 +98,8 @@ def lpsolve(c, G, h, solver=None):
         result = _solve_lp_using_cvxopt(c, G, h, solver=solver)
     elif solver == 'scipy':
         result = _solve_lp_using_scipy(c, G, h)
+    elif solver == 'gurobi':
+        result = _solve_lp_using_gurobi(c, G, h)
     else:
         raise Exception(
             'unknown LP solver "{s}".'.format(s=solver))
@@ -145,6 +158,45 @@ def _solve_lp_using_scipy(c, G, h):
         fun=sol.fun)
 
 
+def _solve_lp_using_gurobi(c, G, h):
+    """Attempt linear optimization using gurobipy."""
+    _assert_have_solver('gurobi')
+    m = gurobi.Model()
+    x = m.addMVar(G.shape[1], lb=-gurobi.GRB.INFINITY)
+    m.addConstr(G@x <= h)
+    m.setObjective(c@x)
+    m.optimize()
+
+    result = dict()
+    if m.Status == gurobi.GRB.OPTIMAL:
+        result['status'] = 0
+        result['x'] = x.x
+        result['fun'] = m.ObjVal
+        return result
+    elif m.Status == gurobi.GRB.INFEASIBLE:
+        result['status'] = 2
+    elif m.Status == gurobi.GRB.UNBOUNDED:
+        result['status'] = 3
+    elif m.Status == gurobi.GRB.INF_OR_UNBD:
+        m.reset(0)
+        m.Params.DualReductions = 0
+        m.optimize()
+        result['message'] = 'Gurobi optimization status was INF_OR_UNBD, so reoptimized with DualReductions=0'
+        if m.Status == gurobi.GRB.INFEASIBLE:
+            result['status'] = 2
+        elif m.Status == gurobi.GRB.UNBOUNDED:
+            result['status'] = 3
+        else:
+            raise ValueError(f'`gurobipy` returned unexpected status value after reoptimizing with DualReductions=0: {m.Status}')
+    else:
+        raise ValueError(f'`gurobipy` returned unexpected status value: {m.Status}')
+
+    result['x'] = None
+    result['fun'] = None
+    return result
+
+
+
 def _assert_have_solver(solver):
     """Raise `RuntimeError` if `solver` is absent."""
     if solver in installed_solvers:

tests/plot_test.py

@@ -1,6 +1,6 @@
 #!/usr/bin/env python
 """Tests for plotting."""
-import unittest
+import pytest
 
 try:
     import matplotlib as mpl
@@ -19,9 +19,9 @@ def add_patch(self, x):
         pass
 
 
-@unittest.skipIf(
+@pytest.mark.skipif(
     mpl is None,
-    '`matplotlib` is not installed')
+    reason='`matplotlib` is not installed')
 def test_plot_transition_arrow():
     p0 = pc.box2poly([[0.0, 1.0], [0.0, 2.0]])
     p1 = pc.box2poly([[0.1, 2.0], [0.0, 2.0]])

tests/polytope_test.py

@@ -1,13 +1,19 @@
 #!/usr/bin/env python
 """Tests for the polytope package."""
 import logging
+from warnings import warn
 
 import numpy as np
-from numpy.testing import assert_allclose
-from numpy.testing import assert_array_equal
+from numpy.testing import assert_allclose, assert_array_equal, assert_equal
+
 import pytest
 import scipy.optimize
 
+try:
+    import gurobipy
+except ImportError:
+    gurobipy = None
+
 import polytope as pc
 import polytope.polytope as alg
 from polytope import solvers
@@ -616,5 +622,34 @@ def test_reduce():
     assert_allclose(u, np.array([[50.], [1.]]), rtol=1e-07, atol=1e-07)
 
 
+@pytest.mark.skipif(
+    gurobipy is None,
+    reason='`gurobipy` is not installed')
+def test_gurobipy_return_same_result_as_scipy():
+    # Try 1-D problem with solution
+    c, A, b = example_1d()
+    result_gurobi = solvers.lpsolve(c, A, b, solver='gurobi')
+    result_scipy = solvers.lpsolve(c, A, b, solver='scipy')
+    assert_equal(result_gurobi['status'], result_scipy['status'])
+    assert_allclose(result_gurobi['x'][0], result_scipy['x'][0])
+    assert_allclose(result_gurobi['fun'], result_scipy['fun'])
+
+    # Try 1-D unbounded problem that may trigger INF_OR_UNBD Gurobi status
+    c = np.array([1])
+    A = np.array([[1]])
+    b = np.array([1])
+    result = solvers.lpsolve(c, A, b, solver='gurobi')
+    result_gurobi = solvers.lpsolve(c, A, b, solver='gurobi')
+    result_scipy = solvers.lpsolve(c, A, b, solver='scipy')
+    assert_equal(result_gurobi['status'], result_scipy['status'])
+    assert_equal(
+        result_gurobi['status'],
+        3,
+        'Optimization status expected to be unbounded, but is not'
+    )
+    if 'message' not in result_gurobi or 'INF_OR_UNBD' not in result_gurobi['message']:
+        warn('Test with Gurobi often results in INF_OR_UNBD but did not')
+
+
 if __name__ == '__main__':
     pass

tests/pytest.ini

@@ -0,0 +1,2 @@
+[pytest]
+addopts = --strict-markers -rs