Merge pull request #857 from Joao-Dionisio/get_infeasibilities

Add recipe for getting infeasible constraints

Author: Opt-Mucca

CHANGELOG.md

@@ -2,6 +2,7 @@
 
 ## Unreleased
 ### Added
+- Added recipe with reformulation for detecting infeasible constraints
 - Wrapped SCIPcreateOrigSol and added tests 
 - Added verbose option for writeProblem and writeParams
 - Expanded locale test

src/pyscipopt/recipes/infeasibilities.py

@@ -0,0 +1,50 @@
+from pyscipopt import Model, quicksum
+
+
+def get_infeasible_constraints(orig_model: Model, verbose=False):
+    """
+    Given a model, adds slack variables to all the constraints and minimizes a binary variable that indicates if they're positive.
+    Positive slack variables correspond to infeasible constraints.
+    """
+    
+    model = Model(sourceModel=orig_model, origcopy=True) # to preserve the model
+
+    slack      = {}
+    aux        = {}
+    binary     = {}
+    aux_binary = {}
+
+    for c in model.getConss():
+
+        slack[c.name] = model.addVar(lb=-float("inf"), name="s_"+c.name) 
+        model.addConsCoeff(c, slack[c.name], 1)
+        binary[c.name] = model.addVar(vtype="B") # Binary variable to get minimum infeasible constraints. See PR #857.
+
+        # getting the absolute value because of <= and >= constraints 
+        aux[c.name] = model.addVar()
+        model.addCons(aux[c.name] >= slack[c.name])
+        model.addCons(aux[c.name] >= -slack[c.name])
+        
+        # modeling aux > 0 => binary = 1 constraint. See https://or.stackexchange.com/q/12142/5352 for an explanation
+        aux_binary[c.name] = model.addVar(vtype="B")
+        model.addCons(binary[c.name]+aux_binary[c.name] == 1)
+        model.addConsSOS1([aux[c.name], aux_binary[c.name]])
+
+    model.setObjective(quicksum(binary[c.name] for c in orig_model.getConss()))
+    model.hideOutput()
+    model.optimize()
+
+    n_infeasibilities_detected = 0
+    for c in binary:
+        if model.isGT(model.getVal(binary[c]), 0):
+            n_infeasibilities_detected += 1
+            print("Constraint %s is causing an infeasibility." % c)
+
+    if verbose:
+        if n_infeasibilities_detected > 0:
+            print("If the constraint names are unhelpful, consider giving them\
+                a suitable name when creating the model with model.addCons(..., name=\"the_name_you_want\")")
+        else:
+            print("Model is feasible.")
+    
+    return n_infeasibilities_detected, aux

src/pyscipopt/scip.pxi

@@ -1607,7 +1607,6 @@ cdef class Model:
         PY_SCIP_CALL(SCIPtightenVarLb(self._scip, var.scip_var, lb, force, &infeasible, &tightened))
         return infeasible, tightened
 
-
     def tightenVarUb(self, Variable var, ub, force=False):
         """Tighten the upper bound in preprocessing or current node, if the bound is tighter.
 
@@ -1624,7 +1623,6 @@ cdef class Model:
         PY_SCIP_CALL(SCIPtightenVarUb(self._scip, var.scip_var, ub, force, &infeasible, &tightened))
         return infeasible, tightened
 
-
     def tightenVarUbGlobal(self, Variable var, ub, force=False):
         """Tighten the global upper bound, if the bound is tighter.
 
@@ -2556,7 +2554,6 @@ cdef class Model:
         PY_SCIP_CALL(SCIPreleaseCons(self._scip, &(<Constraint>cons).scip_cons))
         return disj_cons
 
-
     def getConsNVars(self, Constraint constraint):
         """
         Gets number of variables in a constraint.
@@ -2699,6 +2696,7 @@ cdef class Model:
 
         PY_SCIP_CALL(SCIPaddCons(self._scip, scip_cons))
         return Constraint.create(scip_cons)
+
     def addConsSOS2(self, vars, weights=None, name="SOS2cons",
                 initial=True, separate=True, enforce=True, check=True,
                 propagate=True, local=False, dynamic=False,

tests/test_recipe_infeasibilities.py

@@ -0,0 +1,29 @@
+from pyscipopt import Model
+from pyscipopt.recipes.infeasibilities import get_infeasible_constraints
+
+
+def test_get_infeasible_constraints():
+    m = Model()
+
+    x = m.addVar(lb=0)
+    m.setObjective(2*x)
+
+    m.addCons(x <= 4)
+    
+    n_infeasibilities_detected = get_infeasible_constraints(m)[0]
+    assert n_infeasibilities_detected == 0
+
+    m.addCons(x <= -1)
+
+    n_infeasibilities_detected = get_infeasible_constraints(m)[0]
+    assert n_infeasibilities_detected == 1
+
+    m.addCons(x == 2)
+
+    n_infeasibilities_detected = get_infeasible_constraints(m)[0]
+    assert n_infeasibilities_detected == 1
+
+    m.addCons(x == -4)
+
+    n_infeasibilities_detected = get_infeasible_constraints(m)[0]
+    assert n_infeasibilities_detected == 2