InternalSolver: Better specify and test infeasibility

miplearn/instance.py

@@ -10,6 +10,7 @@ from typing import Any, List
 import numpy as np
 
 from miplearn.types import TrainingSample
+import pyomo.environ as pe
 
 
 class Instance(ABC):
@@ -30,7 +31,7 @@ class Instance(ABC):
     @abstractmethod
     def to_model(self) -> Any:
         """
-        Returns a concrete Pyomo model corresponding to this instance.
+        Returns the optimization model corresponding to this instance.
         """
         pass
 
@@ -163,3 +164,12 @@ class Instance(ABC):
         data = json.dumps(self.__dict__, indent=2).encode("utf-8")
         with gzip.GzipFile(filename, "w") as f:
             f.write(data)
+
+
+class PyomoInstance(Instance, ABC):
+    @abstractmethod
+    def to_model(self) -> pe.ConcreteModel:
+        """
+        Returns the concrete Pyomo model corresponding to this instance.
+        """
+        pass

miplearn/solvers/gurobi.py

@@ -128,8 +128,11 @@ class GurobiSolver(InternalSolver):
             for (idx, var) in vardict.items():
                 var.vtype = self.GRB.BINARY
         log = streams[0].getvalue()
+        opt_value = None
+        if not self.is_infeasible():
+            opt_value = self.model.objVal
         return {
-            "Optimal value": self.model.objVal,
+            "Optimal value": opt_value,
             "Log": log,
         }
 
@@ -173,14 +176,15 @@ class GurobiSolver(InternalSolver):
             if not should_repeat:
                 break
         log = streams[0].getvalue()
-        if self.model.modelSense == 1:
-            sense = "min"
-            lb = self.model.objBound
-            ub = self.model.objVal
-        else:
-            sense = "max"
-            lb = self.model.objVal
-            ub = self.model.objBound
+        ub, lb = None, None
+        sense = "min" if self.model.modelSense == 1 else "max"
+        if self.model.solCount > 0:
+            if self.model.modelSense == 1:
+                lb = self.model.objBound
+                ub = self.model.objVal
+            else:
+                lb = self.model.objVal
+                ub = self.model.objBound
         ws_value = self._extract_warm_start_value(log)
         stats: MIPSolveStats = {
             "Lower bound": lb,
@@ -194,8 +198,10 @@ class GurobiSolver(InternalSolver):
         }
         return stats
 
-    def get_solution(self) -> Dict:
+    def get_solution(self) -> Optional[Dict]:
         self._raise_if_callback()
+        if self.model.solCount == 0:
+            return None
         solution: Dict = {}
         for (varname, vardict) in self._all_vars.items():
             solution[varname] = {}
@@ -228,7 +234,7 @@ class GurobiSolver(InternalSolver):
         var = self._all_vars[var_name][index]
         return self._get_value(var)
 
-    def is_infeasible(self):
+    def is_infeasible(self) -> bool:
         return self.model.status in [self.GRB.INFEASIBLE, self.GRB.INF_OR_UNBD]
 
     def get_dual(self, cid):

miplearn/solvers/internal.py

@@ -4,7 +4,7 @@
 
 import logging
 from abc import ABC, abstractmethod
-from typing import Any, Dict, List
+from typing import Any, Dict, List, Optional
 
 from miplearn.instance import Instance
 from miplearn.types import (
@@ -39,15 +39,13 @@ class InternalSolver(ABC):
         Solves the LP relaxation of the currently loaded instance. After this
         method finishes, the solution can be retrieved by calling `get_solution`.
 
+        This method should not permanently modify the problem. That is, subsequent
+        calls to `solve` should solve the original MIP, not the LP relaxation.
+
         Parameters
         ----------
-        tee: bool
+        tee
             If true, prints the solver log to the screen.
-
-        Returns
-        -------
-        dict
-            A dictionary of solver statistics.
         """
         pass
 
@@ -64,34 +62,27 @@ class InternalSolver(ABC):
 
         Parameters
         ----------
-        iteration_cb: () -> Bool
+        iteration_cb:
             By default, InternalSolver makes a single call to the native `solve`
             method and returns the result. If an iteration callback is provided
             instead, InternalSolver enters a loop, where `solve` and `iteration_cb`
-            are called alternatively. To stop the loop, `iteration_cb` should
-            return False. Any other result causes the solver to loop again.
-        lazy_cb: (internal_solver, model) -> None
+            are called alternatively. To stop the loop, `iteration_cb` should return
+            False. Any other result causes the solver to loop again.
+        lazy_cb:
             This function is called whenever the solver finds a new candidate
-            solution and can be used to add lazy constraints to the model. Only
-            the following operations within the callback are allowed:
-                - Querying the value of a variable, through `get_value(var, idx)`
-                - Querying if a constraint is satisfied, through `is_constraint_satisfied(cobj)`
-                - Adding a new constraint to the problem, through `add_constraint`
+            solution and can be used to add lazy constraints to the model. Only the
+            following operations within the callback are allowed:
+                - Querying the value of a variable
+                - Querying if a constraint is satisfied
+                - Adding a new constraint to the problem
             Additional operations may be allowed by specific subclasses.
-        tee: Bool
+        tee
             If true, prints the solver log to the screen.
-
-        Returns
-        -------
-        dict
-            A dictionary of solver statistics containing the following keys:
-            "Lower bound", "Upper bound", "Wallclock time", "Nodes", "Sense",
-             "Log" and "Warm start value".
         """
         pass
 
     @abstractmethod
-    def get_solution(self) -> Dict:
+    def get_solution(self) -> Optional[Dict]:
         """
         Returns current solution found by the solver.
 
@@ -201,7 +192,7 @@ class InternalSolver(ABC):
         pass
 
     @abstractmethod
-    def set_constraint_rhs(self, cid: str, rhs: str) -> None:
+    def set_constraint_rhs(self, cid: str, rhs: float) -> None:
         pass
 
     @abstractmethod

miplearn/solvers/pyomo/base.py

@@ -11,6 +11,7 @@ from typing import Any, List, Dict, Optional
 import pyomo
 from pyomo import environ as pe
 from pyomo.core import Var, Constraint
+from pyomo.opt import TerminationCondition
 
 from miplearn.instance import Instance
 from miplearn.solvers import RedirectOutput
@@ -44,6 +45,7 @@ class BasePyomoSolver(InternalSolver):
         self._obj_sense = None
         self._varname_to_var = {}
         self._cname_to_constr = {}
+        self._termination_condition = None
         for (key, value) in params.items():
             self._pyomo_solver.options[key] = value
 
@@ -65,8 +67,11 @@ class BasePyomoSolver(InternalSolver):
         for var in self._bin_vars:
             var.domain = pyomo.core.base.set_types.Binary
             self._pyomo_solver.update_var(var)
+        opt_value = None
+        if not self.is_infeasible():
+            opt_value = results["Problem"][0]["Lower bound"]
         return {
-            "Optimal value": results["Problem"][0]["Lower bound"],
+            "Optimal value": opt_value,
             "Log": streams[0].getvalue(),
         }
 
@@ -100,9 +105,14 @@ class BasePyomoSolver(InternalSolver):
         log = streams[0].getvalue()
         node_count = self._extract_node_count(log)
         ws_value = self._extract_warm_start_value(log)
+        self._termination_condition = results["Solver"][0]["Termination condition"]
+        lb, ub = None, None
+        if not self.is_infeasible():
+            lb = results["Problem"][0]["Lower bound"]
+            ub = results["Problem"][0]["Upper bound"]
         stats: MIPSolveStats = {
-            "Lower bound": results["Problem"][0]["Lower bound"],
-            "Upper bound": results["Problem"][0]["Upper bound"],
+            "Lower bound": lb,
+            "Upper bound": ub,
             "Wallclock time": total_wallclock_time,
             "Sense": self._obj_sense,
             "Log": log,
@@ -112,7 +122,9 @@ class BasePyomoSolver(InternalSolver):
         }
         return stats
 
-    def get_solution(self) -> Dict:
+    def get_solution(self) -> Optional[Dict]:
+        if self.is_infeasible():
+            return None
         solution: Dict = {}
         for var in self.model.component_objects(Var):
             solution[str(var)] = {}
@@ -276,8 +288,8 @@ class BasePyomoSolver(InternalSolver):
     def set_constraint_rhs(self, cid, rhs):
         raise Exception("Not implemented")
 
-    def is_infeasible(self):
-        raise Exception("Not implemented")
+    def is_infeasible(self) -> bool:
+        return self._termination_condition == TerminationCondition.infeasible
 
     def get_dual(self, cid):
         raise Exception("Not implemented")

miplearn/solvers/tests/__init__.py

@@ -3,8 +3,11 @@
 #  Released under the modified BSD license. See COPYING.md for more details.
 
 from inspect import isclass
-from typing import List, Callable
+from typing import List, Callable, Any
 
+from pyomo import environ as pe
+
+from miplearn.instance import Instance, PyomoInstance
 from miplearn.problems.knapsack import KnapsackInstance, GurobiKnapsackInstance
 from miplearn.solvers.gurobi import GurobiSolver
 from miplearn.solvers.internal import InternalSolver
@@ -13,28 +16,55 @@ from miplearn.solvers.pyomo.gurobi import GurobiPyomoSolver
 from miplearn.solvers.pyomo.xpress import XpressPyomoSolver
 
 
-def _get_instance(solver):
-    def _is_subclass_or_instance(obj, parent_class):
-        return isinstance(obj, parent_class) or (
-            isclass(obj) and issubclass(obj, parent_class)
-        )
+class InfeasiblePyomoInstance(PyomoInstance):
+    def to_model(self) -> pe.ConcreteModel:
+        model = pe.ConcreteModel()
+        model.x = pe.Var(domain=pe.Binary)
+        model.OBJ = pe.Objective(expr=model.x, sense=pe.maximize)
+        model.eq = pe.Constraint(expr=model.x >= 2)
+        return model
+
+
+class InfeasibleGurobiInstance(Instance):
+    def to_model(self) -> Any:
+        import gurobipy as gp
+        from gurobipy import GRB
+
+        model = gp.Model()
+        x = model.addVars(1, vtype=GRB.BINARY, name="x")
+        model.addConstr(x[0] >= 2)
+        model.setObjective(x[0])
+        return model
+
+
+def _is_subclass_or_instance(obj, parent_class):
+    return isinstance(obj, parent_class) or (
+        isclass(obj) and issubclass(obj, parent_class)
+    )
 
+
+def _get_knapsack_instance(solver):
     if _is_subclass_or_instance(solver, BasePyomoSolver):
         return KnapsackInstance(
             weights=[23.0, 26.0, 20.0, 18.0],
             prices=[505.0, 352.0, 458.0, 220.0],
             capacity=67.0,
         )
-
     if _is_subclass_or_instance(solver, GurobiSolver):
         return GurobiKnapsackInstance(
             weights=[23.0, 26.0, 20.0, 18.0],
             prices=[505.0, 352.0, 458.0, 220.0],
             capacity=67.0,
         )
-
     assert False
 
 
+def _get_infeasible_instance(solver):
+    if _is_subclass_or_instance(solver, BasePyomoSolver):
+        return InfeasiblePyomoInstance()
+    if _is_subclass_or_instance(solver, GurobiSolver):
+        return InfeasibleGurobiInstance()
+
+
 def _get_internal_solvers() -> List[Callable[[], InternalSolver]]:
     return [GurobiPyomoSolver, GurobiSolver, XpressPyomoSolver]

miplearn/solvers/tests/test_internal_solver.py

@@ -11,7 +11,11 @@ import pyomo.environ as pe
 from miplearn.solvers import RedirectOutput
 from miplearn.solvers.gurobi import GurobiSolver
 from miplearn.solvers.pyomo.base import BasePyomoSolver
-from miplearn.solvers.tests import _get_instance, _get_internal_solvers
+from miplearn.solvers.tests import (
+    _get_knapsack_instance,
+    _get_internal_solvers,
+    _get_infeasible_instance,
+)
 
 logger = logging.getLogger(__name__)
 
@@ -30,7 +34,7 @@ def test_redirect_output():
 def test_internal_solver_warm_starts():
     for solver_class in _get_internal_solvers():
         logger.info("Solver: %s" % solver_class)
-        instance = _get_instance(solver_class)
+        instance = _get_knapsack_instance(solver_class)
         model = instance.to_model()
         solver = solver_class()
         solver.set_instance(instance, model)
@@ -82,7 +86,7 @@ def test_internal_solver():
     for solver_class in _get_internal_solvers():
         logger.info("Solver: %s" % solver_class)
 
-        instance = _get_instance(solver_class)
+        instance = _get_knapsack_instance(solver_class)
         model = instance.to_model()
         solver = solver_class()
         solver.set_instance(instance, model)
@@ -158,10 +162,26 @@ def test_internal_solver():
             assert round(stats["Lower bound"]) == 1179.0
 
 
+def test_infeasible_instance():
+    for solver_class in _get_internal_solvers():
+        instance = _get_infeasible_instance(solver_class)
+        solver = solver_class()
+        solver.set_instance(instance)
+        stats = solver.solve()
+
+        assert solver.get_solution() is None
+        assert stats["Upper bound"] is None
+        assert stats["Lower bound"] is None
+
+        stats = solver.solve_lp()
+        assert solver.get_solution() is None
+        assert stats["Optimal value"] is None
+
+
 def test_iteration_cb():
     for solver_class in _get_internal_solvers():
         logger.info("Solver: %s" % solver_class)
-        instance = _get_instance(solver_class)
+        instance = _get_knapsack_instance(solver_class)
         solver = solver_class()
         solver.set_instance(instance)
         count = 0

miplearn/solvers/tests/test_lazy_cb.py

@@ -5,14 +5,14 @@
 import logging
 
 from miplearn.solvers.gurobi import GurobiSolver
-from miplearn.solvers.tests import _get_instance
+from miplearn.solvers.tests import _get_knapsack_instance
 
 logger = logging.getLogger(__name__)
 
 
 def test_lazy_cb():
     solver = GurobiSolver()
-    instance = _get_instance(solver)
+    instance = _get_knapsack_instance(solver)
     model = instance.to_model()
 
     def lazy_cb(cb_solver, cb_model):

miplearn/solvers/tests/test_learning_solver.py

@@ -9,7 +9,7 @@ import os
 
 from miplearn.solvers.gurobi import GurobiSolver
 from miplearn.solvers.learning import LearningSolver
-from miplearn.solvers.tests import _get_instance, _get_internal_solvers
+from miplearn.solvers.tests import _get_knapsack_instance, _get_internal_solvers
 
 logger = logging.getLogger(__name__)
 
@@ -18,7 +18,7 @@ def test_learning_solver():
     for mode in ["exact", "heuristic"]:
         for internal_solver in _get_internal_solvers():
             logger.info("Solver: %s" % internal_solver)
-            instance = _get_instance(internal_solver)
+            instance = _get_knapsack_instance(internal_solver)
             solver = LearningSolver(
                 solver=internal_solver,
                 mode=mode,
@@ -50,7 +50,7 @@ def test_learning_solver():
 def test_solve_without_lp():
     for internal_solver in _get_internal_solvers():
         logger.info("Solver: %s" % internal_solver)
-        instance = _get_instance(internal_solver)
+        instance = _get_knapsack_instance(internal_solver)
         solver = LearningSolver(
             solver=internal_solver,
             solve_lp_first=False,
@@ -62,7 +62,7 @@ def test_solve_without_lp():
 
 def test_parallel_solve():
     for internal_solver in _get_internal_solvers():
-        instances = [_get_instance(internal_solver) for _ in range(10)]
+        instances = [_get_knapsack_instance(internal_solver) for _ in range(10)]
         solver = LearningSolver(solver=internal_solver)
         results = solver.parallel_solve(instances, n_jobs=3)
         assert len(results) == 10
@@ -76,7 +76,7 @@ def test_solve_fit_from_disk():
         # Create instances and pickle them
         filenames = []
         for k in range(3):
-            instance = _get_instance(internal_solver)
+            instance = _get_knapsack_instance(internal_solver)
             with tempfile.NamedTemporaryFile(suffix=".pkl", delete=False) as file:
                 filenames += [file.name]
                 pickle.dump(instance, file)
@@ -114,7 +114,7 @@ def test_solve_fit_from_disk():
 
 def test_simulate_perfect():
     internal_solver = GurobiSolver
-    instance = _get_instance(internal_solver)
+    instance = _get_knapsack_instance(internal_solver)
     with tempfile.NamedTemporaryFile(suffix=".pkl", delete=False) as tmp:
         pickle.dump(instance, tmp)
         tmp.flush()

miplearn/types.py

@@ -22,7 +22,7 @@ TrainingSample = TypedDict(
 LPSolveStats = TypedDict(
     "LPSolveStats",
     {
-        "Optimal value": float,
+        "Optimal value": Optional[float],
         "Log": str,
     },
 )