ConvertTight: Detect and fix sub-optimality

miplearn/components/steps/convert_tight.py

@@ -52,10 +52,11 @@ class ConvertTightIneqsIntoEqsStep(Component):
         )
         y = self.predict(x)
 
-        self.total_converted = 0
-        self.total_restored = 0
-        self.total_kept = 0
-        self.total_iterations = 0
+        self.n_converted = 0
+        self.n_restored = 0
+        self.n_kept = 0
+        self.n_infeasible_iterations = 0
+        self.n_suboptimal_iterations = 0
         for category in y.keys():
             for i in range(len(y[category])):
                 if y[category][i][0] == 1:
@@ -64,17 +65,18 @@ class ConvertTightIneqsIntoEqsStep(Component):
                     self.original_sense[cid] = s
                     solver.internal_solver.set_constraint_sense(cid, "=")
                     self.converted += [cid]
-                    self.total_converted += 1
+                    self.n_converted += 1
                 else:
-                    self.total_kept += 1
-        logger.info(f"Converted {self.total_converted} inequalities")
+                    self.n_kept += 1
+        logger.info(f"Converted {self.n_converted} inequalities")
 
     def after_solve(self, solver, instance, model, results):
         instance.slacks = solver.internal_solver.get_inequality_slacks()
-        results["ConvertTight: Kept"] = self.total_kept
-        results["ConvertTight: Converted"] = self.total_converted
-        results["ConvertTight: Restored"] = self.total_restored
-        results["ConvertTight: Iterations"] = self.total_iterations
+        results["ConvertTight: Kept"] = self.n_kept
+        results["ConvertTight: Converted"] = self.n_converted
+        results["ConvertTight: Restored"] = self.n_restored
+        results["ConvertTight: Inf iterations"] = self.n_infeasible_iterations
+        results["ConvertTight: Subopt iterations"] = self.n_suboptimal_iterations
 
     def fit(self, training_instances):
         logger.debug("Extracting x and y...")
@@ -173,21 +175,56 @@ class ConvertTightIneqsIntoEqsStep(Component):
     def iteration_cb(self, solver, instance, model):
         if not self.check_converted:
             return False
+
         logger.debug("Checking converted inequalities...")
+        is_infeasible, is_suboptimal = False, False
         restored = []
+
+        def check_pi(msense, csense, pi):
+            if csense == "=":
+                return True
+            if msense == "max":
+                if csense == "<":
+                    return pi >= 0
+                else:
+                    return pi <= 0
+            else:
+                if csense == ">":
+                    return pi >= 0
+                else:
+                    return pi <= 0
+
+        def restore(cid):
+            nonlocal restored
+            csense = self.original_sense[cid]
+            solver.internal_solver.set_constraint_sense(cid, csense)
+            restored += [cid]
+
         if solver.internal_solver.is_infeasible():
             for cid in self.converted:
-                f = solver.internal_solver.get_farkas_dual(cid)
-                if abs(f) > 0:
-                    s = self.original_sense[cid]
-                    solver.internal_solver.set_constraint_sense(cid, s)
-                    restored += [cid]
-            for cid in restored:
-                self.converted.remove(cid)
+                pi = solver.internal_solver.get_dual(cid)
+                if abs(pi) > 0:
+                    is_infeasible = True
+                    restore(cid)
+        else:
+            for cid in self.converted:
+                pi = solver.internal_solver.get_dual(cid)
+                csense = self.original_sense[cid]
+                msense = solver.internal_solver.get_sense()
+                if not check_pi(msense, csense, pi):
+                    is_suboptimal = True
+                    restore(cid)
+
+        for cid in restored:
+            self.converted.remove(cid)
+
         if len(restored) > 0:
-            self.total_restored += len(restored)
+            self.n_restored += len(restored)
+            if is_infeasible:
+                self.n_infeasible_iterations += 1
+            if is_suboptimal:
+                self.n_suboptimal_iterations += 1
             logger.info(f"Restored {len(restored)} inequalities")
-            self.total_iterations += 1
             return True
         else:
             return False

miplearn/components/steps/tests/convert_tight_test.py

@@ -72,3 +72,26 @@ def test_convert_tight_infeasibility():
     instance = TestInstance()
     solver.solve(instance)
     assert instance.lower_bound == 5.0
+
+
+def test_convert_tight_suboptimality():
+    comp = ConvertTightIneqsIntoEqsStep(
+        check_converted=True,
+    )
+    comp.classifiers = {
+        "c1": Mock(spec=Classifier),
+        "c2": Mock(spec=Classifier),
+        "c3": Mock(spec=Classifier),
+    }
+    comp.classifiers["c1"].predict_proba = Mock(return_value=[[0, 1]])
+    comp.classifiers["c2"].predict_proba = Mock(return_value=[[1, 0]])
+    comp.classifiers["c3"].predict_proba = Mock(return_value=[[0, 1]])
+
+    solver = LearningSolver(
+        solver=GurobiSolver(params={}),
+        components=[comp],
+        solve_lp_first=False,
+    )
+    instance = TestInstance()
+    solver.solve(instance)
+    assert instance.lower_bound == 5.0

miplearn/solvers/gurobi.py

@@ -162,6 +162,12 @@ class GurobiSolver(InternalSolver):
             "Warm start value": self._extract_warm_start_value(log),
         }
 
+    def get_sense(self):
+        if self.model.modelSense == 1:
+            return "min"
+        else:
+            return "max"
+
     def get_solution(self):
         self._raise_if_callback()
 
@@ -179,9 +185,12 @@ class GurobiSolver(InternalSolver):
     def is_infeasible(self):
         return self.model.status in [self.GRB.INFEASIBLE, self.GRB.INF_OR_UNBD]
 
-    def get_farkas_dual(self, cid):
+    def get_dual(self, cid):
         c = self.model.getConstrByName(cid)
-        return c.farkasDual
+        if self.is_infeasible():
+            return c.farkasDual
+        else:
+            return c.pi
 
     def _get_value(self, var):
         if self.cb_where == self.GRB.Callback.MIPSOL:

miplearn/solvers/internal.py

@@ -200,11 +200,20 @@ class InternalSolver(ABC):
         pass
 
     @abstractmethod
-    def get_farkas_dual(self, cid):
+    def get_dual(self, cid):
         """
-        If the model is infeasible, returns a portion of the infeasibility certificate
-        corresponding to the given constraint. If the model is feasible, calling this
-        function raises an error.
+        If the model is feasible and has been solved to optimality, returns the optimal
+        value of the dual variable associated with this constraint. If the model is infeasible,
+        returns a portion of the infeasibility certificate corresponding to the given constraint.
+
+        Solve must be called prior to this method.
+        """
+        pass
+
+    @abstractmethod
+    def get_sense(self):
+        """
+        Returns the sense of the problem (either "min" or "max").
         """
         pass
 

miplearn/solvers/pyomo/base.py

@@ -267,5 +267,8 @@ class BasePyomoSolver(InternalSolver):
     def is_infeasible(self):
         raise Exception("Not implemented")
 
-    def get_farkas_dual(self, cid):
+    def get_dual(self, cid):
+        raise Exception("Not implemented")
+
+    def get_sense(self):
         raise Exception("Not implemented")