reoptimize with DualReductions=0 if Gurobi INF_OR_UNBD

Following guidance at https://www.gurobi.com/documentation/current/refman/optimization_status_codes.html

polytope/solvers.py

@@ -175,8 +175,19 @@ def _solve_lp_using_gurobi(c, G, h):
         return result
     elif m.Status == gurobi.GRB.INFEASIBLE:
         result['status'] = 2
-    elif m.Status == gurobi.GRB.INF_OR_UNBD or m.Status == gurobi.GRB.UNBOUNDED:
+    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}')
 

tests/polytope_test.py

@@ -1,6 +1,7 @@
 #!/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, assert_array_equal, assert_equal
@@ -625,13 +626,30 @@ def test_reduce():
     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