approximation of linear, quadratic disjunction constraints

Author: ctdunc

src/pyscipopt/scip.pxd

@@ -1493,6 +1493,23 @@ cdef extern from "scip/cons_sos2.h":
                                    SCIP_CONS* cons,
                                    SCIP_VAR* var)
 
+cdef extern from "scip/cons_disjunction.h":
+    SCIP_RETCODE SCIPcreateConsDisjunction(SCIP *scip,
+                                            SCIP_CONS **cons,
+                                            const char *name,
+                                            int nconss,
+                                            SCIP_CONS **conss,
+                                            SCIP_CONS *relaxcons,
+                                            SCIP_Bool initial,
+                                            SCIP_Bool enforce,
+                                            SCIP_Bool check,
+                                            SCIP_Bool local,
+                                            SCIP_Bool modifiable,
+                                            SCIP_Bool dynamic)
+    SCIP_RETCODE SCIPaddConsElemDisjunction(SCIP *scip,
+                                            SCIP_CONS *cons,
+                                            SCIP_CONS *addcons)
+
 cdef extern from "scip/cons_and.h":
     SCIP_RETCODE SCIPcreateConsAnd(SCIP* scip,
                                          SCIP_CONS** cons,

src/pyscipopt/scip.pxi

@@ -2122,7 +2122,7 @@ cdef class Model:
             return self._addGenNonlinearCons(cons, **kwargs)
         else:
             return self._addNonlinearCons(cons, **kwargs)
-
+    
     def addConss(self, conss, name='', initial=True, separate=True,
                  enforce=True, check=True, propagate=True, local=False,
                  modifiable=False, dynamic=False, removable=False,
@@ -2188,6 +2188,106 @@ cdef class Model:
 
         return constraints
 
+    def addConsDisjunction(self, conss, name = '', initial = True, 
+        relaxcons = None, enforce=True, check =True, 
+        local=False, modifiable = False, dynamic = False):
+        
+        def ensure_iterable(elem, length):
+            if isinstance(elem, Iterable):
+                return elem
+            else:
+                return list(repeat(elem, length))
+        assert isinstance(conss, Iterable), "Given constraint list is not iterable"
+
+        conss = list(conss)
+        n_conss = len(conss)
+        assert n_conss >= 2, "Given constraint list contains fewer than 2 items!"
+
+        cdef SCIP_CONS* disj_cons
+
+        cdef SCIP_CONS* scip_cons
+        
+        cdef SCIP_EXPR* scip_expr
+
+        PY_SCIP_CALL(SCIPcreateConsDisjunction(
+            self._scip, &disj_cons, str_conversion(name), 0, &scip_cons, NULL, 
+            initial, enforce, check, local, modifiable, dynamic
+        ))
+
+        #PY_SCIP_CALL(SCIPreleaseCons(self._scip, &scip_cons))
+
+        for i, cons in enumerate(conss):
+            deg = cons.expr.degree()
+            assert isinstance(cons, ExprCons), "given constraint is not ExprCons but %s" % cons.__class__.__name__     
+
+            kwargs = dict(name='c'+str(SCIPgetNConss(self._scip)+1),initial=True,separate=True,
+                enforce=True, check=True, propagate=True, local=False,
+                modifiable=False, dynamic=False, removable=False,
+                stickingatnode=False)
+
+            kwargs['lhs'] = -SCIPinfinity(self._scip) if cons._lhs is None else cons._lhs
+            kwargs['rhs'] =  SCIPinfinity(self._scip) if cons._rhs is None else cons._rhs
+            terms = cons.expr.terms
+
+            if deg <= 1:
+                nvars = len(terms.items())
+                vars_array = <SCIP_VAR**> malloc(nvars * sizeof(SCIP_VAR*))
+                coeffs_array = <SCIP_Real*> malloc(nvars * sizeof(SCIP_Real))
+
+                for i, (key, coeff) in enumerate(terms.items()):
+                    vars_array[i] = <SCIP_VAR*>(<Variable>key[0]).scip_var
+                    coeffs_array[i] = <SCIP_Real>coeff
+
+                PY_SCIP_CALL(SCIPcreateConsLinear(
+                    self._scip, &scip_cons, str_conversion(kwargs['name']), nvars, vars_array, coeffs_array,
+                    kwargs['lhs'], kwargs['rhs'], kwargs['initial'],
+                    kwargs['separate'], kwargs['enforce'], kwargs['check'],
+                    kwargs['propagate'], kwargs['local'], kwargs['modifiable'],
+                    kwargs['dynamic'], kwargs['removable'], kwargs['stickingatnode']))
+                PY_SCIP_CALL(SCIPaddConsElemDisjunction(self._scip, disj_cons, scip_cons)) 
+                PY_SCIP_CALL(SCIPreleaseCons(self._scip, &scip_cons))
+                free(vars_array)
+                free(coeffs_array)
+            elif deg <=2:
+                PY_SCIP_CALL(SCIPcreateConsQuadraticNonlinear(
+                    self._scip, &scip_cons, str_conversion(kwargs['name']),
+                    0, NULL, NULL,        # linear
+                    0, NULL, NULL, NULL,  # quadratc
+                    kwargs['lhs'], kwargs['rhs'],
+                    kwargs['initial'], kwargs['separate'], kwargs['enforce'],
+                    kwargs['check'], kwargs['propagate'], kwargs['local'],
+                    kwargs['modifiable'], kwargs['dynamic'], kwargs['removable']))
+                for v, c in terms.items():
+                    if len(v) == 1: # linear
+                        var = <Variable>v[0]
+                        PY_SCIP_CALL(SCIPaddLinearVarNonlinear(self._scip, scip_cons, var.scip_var, c))
+                    else: # quadratic
+                        assert len(v) == 2, 'term length must be 1 or 2 but it is %s' % len(v)
+
+                        varexprs = <SCIP_EXPR**> malloc(2 * sizeof(SCIP_EXPR*))
+                        var1, var2 = <Variable>v[0], <Variable>v[1]
+                        PY_SCIP_CALL( SCIPcreateExprVar(self._scip, &varexprs[0], var1.scip_var, NULL, NULL) )
+                        PY_SCIP_CALL( SCIPcreateExprVar(self._scip, &varexprs[1], var2.scip_var, NULL, NULL) )
+                        PY_SCIP_CALL( SCIPcreateExprProduct(self._scip, &scip_expr, 2, varexprs, 1.0, NULL, NULL) )
+
+                        PY_SCIP_CALL( SCIPaddExprNonlinear(self._scip, scip_cons, scip_expr, c) )
+
+                        PY_SCIP_CALL(SCIPaddConsElemDisjunction(self._scip, disj_cons, scip_cons)) 
+                        PY_SCIP_CALL( SCIPreleaseExpr(self._scip, &scip_expr) )
+                        PY_SCIP_CALL(SCIPreleaseCons(self._scip, &scip_cons))
+                        PY_SCIP_CALL( SCIPreleaseExpr(self._scip, &varexprs[1]) )
+                        PY_SCIP_CALL( SCIPreleaseExpr(self._scip, &varexprs[0]) )
+                        free(varexprs)
+            else:
+                raise NotImplementedError("Only Linear Expressions are currently supported")
+
+
+        PY_SCIP_CALL(SCIPaddCons(self._scip, disj_cons))
+        PyCons = Constraint.create(disj_cons)
+        PY_SCIP_CALL(SCIPreleaseCons(self._scip, &disj_cons))
+        return PyCons
+
+
     def getConsNVars(self, Constraint constraint):
         """
         Gets number of variables in a constraint.

tests/test_cons.py

@@ -131,6 +131,23 @@ def test_cons_indicator_fail():
     m.setSolVal(sol, binvar, 0)
     assert m.checkSol(sol)  # solution should be feasible
 
+def test_addConsDisjunction():
+    m = Model()
+
+    x1 = m.addVar(vtype="C", lb=-1, ub=1)
+    x2 = m.addVar(vtype="C", lb=-3, ub=3)
+    o = m.addVar(vtype="C")
+
+    c = m.addConsDisjunction([x1 <= 0, x2 <= 0])
+    m.addCons(o <= x1 + x2)
+
+    m.setObjective(o, "maximize")
+    m.optimize()
+
+    assert m.isEQ(m.getVal(x1), 0.0)
+    assert m.isEQ(m.getVal(x2), 3.0)
+    assert m.isEQ(m.getVal(o), 3.0)
+    #assert c.getConshdlrName() == "disjunction"
 
 def test_addConsCardinality():
     m = Model()

tests/test_disjunction.py

@@ -0,0 +1,2 @@
+
+