Merge pull request Pyomo#3293 from jsiirola/constraint-store-expr-only

Constraint: only store the original expression (not lower/body/upper)

doc/OnlineDocs/contributed_packages/gdpopt.rst

@@ -93,10 +93,10 @@ An example that includes the modeling approach may be found below.
     Variables:
       x : Size=1, Index=None
           Key  : Lower : Value : Upper : Fixed : Stale : Domain
-          None :  -1.2 :   0.0 :     2 : False : False :  Reals
+          None :  -1.2 :     0 :     2 : False : False :  Reals
       y : Size=1, Index=None
           Key  : Lower : Value : Upper : Fixed : Stale : Domain
-          None :   -10 :   1.0 :    10 : False : False :  Reals
+          None :   -10 :     1 :    10 : False : False :  Reals
   <BLANKLINE>
     Objectives:
       objective : Size=1, Index=None, Active=True
@@ -106,7 +106,7 @@ An example that includes the modeling approach may be found below.
     Constraints:
       c : Size=1
           Key  : Lower : Body : Upper
-          None :   1.0 :  1.0 :   1.0
+          None :   1.0 :    1 :   1.0
 
 .. note:: 
 

pyomo/contrib/appsi/base.py

@@ -1007,7 +1007,7 @@ def add_constraints(self, cons: List[ConstraintData]):
                 raise ValueError(
                     'constraint {name} has already been added'.format(name=con.name)
                 )
-            self._active_constraints[con] = (con.lower, con.body, con.upper)
+            self._active_constraints[con] = con.expr
             if self.use_extensions and cmodel_available:
                 tmp = cmodel.prep_for_repn(con.body, self._expr_types)
             else:
@@ -1363,40 +1363,13 @@ def update(self, timer: HierarchicalTimer = None):
         cons_to_remove_and_add = dict()
         need_to_set_objective = False
         if config.update_constraints:
-            cons_to_update = list()
-            sos_to_update = list()
             for c in current_cons_dict.keys():
-                if c not in new_cons_set:
-                    cons_to_update.append(c)
+                if c not in new_cons_set and c.expr is not self._active_constraints[c]:
+                    cons_to_remove_and_add[c] = None
+            sos_to_update = []
             for c in current_sos_dict.keys():
                 if c not in new_sos_set:
                     sos_to_update.append(c)
-            for c in cons_to_update:
-                lower, body, upper = self._active_constraints[c]
-                new_lower, new_body, new_upper = c.lower, c.body, c.upper
-                if new_body is not body:
-                    cons_to_remove_and_add[c] = None
-                    continue
-                if new_lower is not lower:
-                    if (
-                        type(new_lower) is NumericConstant
-                        and type(lower) is NumericConstant
-                        and new_lower.value == lower.value
-                    ):
-                        pass
-                    else:
-                        cons_to_remove_and_add[c] = None
-                        continue
-                if new_upper is not upper:
-                    if (
-                        type(new_upper) is NumericConstant
-                        and type(upper) is NumericConstant
-                        and new_upper.value == upper.value
-                    ):
-                        pass
-                    else:
-                        cons_to_remove_and_add[c] = None
-                        continue
             self.remove_sos_constraints(sos_to_update)
             self.add_sos_constraints(sos_to_update)
         timer.stop('cons')

pyomo/contrib/appsi/cmodel/src/fbbt_model.cpp

@@ -205,7 +205,7 @@ void process_fbbt_constraints(FBBTModel *model, PyomoExprTypes &expr_types,
   py::handle con_body;
 
   for (py::handle c : cons) {
-    lower_body_upper = active_constraints[c];
+    lower_body_upper = c.attr("to_bounded_expression")();
     con_lb = lower_body_upper[0];
     con_body = lower_body_upper[1];
     con_ub = lower_body_upper[2];

pyomo/contrib/appsi/cmodel/src/lp_writer.cpp

@@ -289,7 +289,7 @@ void process_lp_constraints(py::list cons, py::object writer) {
   py::object nonlinear_expr;
   PyomoExprTypes expr_types = PyomoExprTypes();
   for (py::handle c : cons) {
-    lower_body_upper = active_constraints[c];
+    lower_body_upper = c.attr("to_bounded_expression")();
     cname = getSymbol(c, labeler);
     repn = generate_standard_repn(
         lower_body_upper[1], "compute_values"_a = false, "quadratic"_a = true);

pyomo/contrib/appsi/cmodel/src/nl_writer.cpp

@@ -527,7 +527,7 @@ void process_nl_constraints(NLWriter *nl_writer, PyomoExprTypes &expr_types,
   py::handle repn_nonlinear_expr;
 
   for (py::handle c : cons) {
-    lower_body_upper = active_constraints[c];
+    lower_body_upper = c.attr("to_bounded_expression")();
     repn = generate_standard_repn(
         lower_body_upper[1], "compute_values"_a = false, "quadratic"_a = false);
     _const = appsi_expr_from_pyomo_expr(repn.attr("constant"), var_map,

pyomo/contrib/community_detection/community_graph.py

@@ -123,7 +123,7 @@ def generate_model_graph(
         # Create a list of the variable numbers that occur in the given constraint equation
         numbered_variables_in_constraint_equation = [
             component_number_map[constraint_variable]
-            for constraint_variable in identify_variables(model_constraint.body)
+            for constraint_variable in identify_variables(model_constraint.expr)
         ]
 
         # Update constraint_variable_map

pyomo/contrib/fbbt/expression_bounds_walker.py

@@ -232,15 +232,15 @@ def _handle_unknowable_bounds(visitor, node, arg):
 
 
 def _handle_equality(visitor, node, arg1, arg2):
-    return eq(*arg1, *arg2)
+    return eq(*arg1, *arg2, feasibility_tol=visitor.feasibility_tol)
 
 
 def _handle_inequality(visitor, node, arg1, arg2):
-    return ineq(*arg1, *arg2)
+    return ineq(*arg1, *arg2, feasibility_tol=visitor.feasibility_tol)
 
 
 def _handle_ranged(visitor, node, arg1, arg2, arg3):
-    return ranged(*arg1, *arg2, *arg3)
+    return ranged(*arg1, *arg2, *arg3, feasibility_tol=visitor.feasibility_tol)
 
 
 def _handle_expr_if(visitor, node, arg1, arg2, arg3):

pyomo/contrib/fbbt/fbbt.py

@@ -12,6 +12,7 @@
 from collections import defaultdict
 from pyomo.common.collections import ComponentMap, ComponentSet
 from pyomo.contrib.fbbt.expression_bounds_walker import ExpressionBoundsVisitor
+import pyomo.core.expr.relational_expr as relational_expr
 import pyomo.core.expr.numeric_expr as numeric_expr
 from pyomo.core.expr.visitor import (
     ExpressionValueVisitor,
@@ -80,6 +81,27 @@ class FBBTException(PyomoException):
     pass
 
 
+def _prop_bnds_leaf_to_root_equality(visitor, node, arg1, arg2):
+    bnds_dict = visitor.bnds_dict
+    bnds_dict[node] = interval.eq(
+        *bnds_dict[arg1], *bnds_dict[arg2], visitor.feasibility_tol
+    )
+
+
+def _prop_bnds_leaf_to_root_inequality(visitor, node, arg1, arg2):
+    bnds_dict = visitor.bnds_dict
+    bnds_dict[node] = interval.ineq(
+        *bnds_dict[arg1], *bnds_dict[arg2], visitor.feasibility_tol
+    )
+
+
+def _prop_bnds_leaf_to_root_ranged(visitor, node, arg1, arg2, arg3):
+    bnds_dict = visitor.bnds_dict
+    bnds_dict[node] = interval.ranged(
+        *bnds_dict[arg1], *bnds_dict[arg2], *bnds_dict[arg3], visitor.feasibility_tol
+    )
+
+
 def _prop_bnds_leaf_to_root_ProductExpression(visitor, node, arg1, arg2):
     """
 
@@ -367,6 +389,9 @@ def _prop_bnds_leaf_to_root_NamedExpression(visitor, node, expr):
         numeric_expr.UnaryFunctionExpression: _prop_bnds_leaf_to_root_UnaryFunctionExpression,
         numeric_expr.LinearExpression: _prop_bnds_leaf_to_root_SumExpression,
         numeric_expr.AbsExpression: _prop_bnds_leaf_to_root_abs,
+        relational_expr.EqualityExpression: _prop_bnds_leaf_to_root_equality,
+        relational_expr.InequalityExpression: _prop_bnds_leaf_to_root_inequality,
+        relational_expr.RangedExpression: _prop_bnds_leaf_to_root_ranged,
         ExpressionData: _prop_bnds_leaf_to_root_NamedExpression,
         ScalarExpression: _prop_bnds_leaf_to_root_NamedExpression,
         ObjectiveData: _prop_bnds_leaf_to_root_NamedExpression,
@@ -375,6 +400,43 @@ def _prop_bnds_leaf_to_root_NamedExpression(visitor, node, expr):
 )
 
 
+def _prop_bnds_root_to_leaf_equality(node, bnds_dict, feasibility_tol):
+    assert bnds_dict[node][1]  # This expression is feasible
+    arg1, arg2 = node.args
+    lb1, ub1 = bnds_dict[arg1]
+    lb2, ub2 = bnds_dict[arg2]
+    bnds_dict[arg1] = bnds_dict[arg2] = max(lb1, lb2), min(ub1, ub2)
+
+
+def _prop_bnds_root_to_leaf_inequality(node, bnds_dict, feasibility_tol):
+    assert bnds_dict[node][1]  # This expression is feasible
+    arg1, arg2 = node.args
+    lb1, ub1 = bnds_dict[arg1]
+    lb2, ub2 = bnds_dict[arg2]
+    if lb1 > lb2:
+        bnds_dict[arg2] = lb1, ub2
+    if ub1 > ub2:
+        bnds_dict[arg1] = lb1, ub2
+
+
+def _prop_bnds_root_to_leaf_ranged(node, bnds_dict, feasibility_tol):
+    assert bnds_dict[node][1]  # This expression is feasible
+    arg1, arg2, arg3 = node.args
+    lb1, ub1 = bnds_dict[arg1]
+    lb2, ub2 = bnds_dict[arg2]
+    lb3, ub3 = bnds_dict[arg3]
+    if lb1 > lb2:
+        bnds_dict[arg2] = lb1, ub2
+        lb2 = lb1
+    if lb2 > lb3:
+        bnds_dict[arg3] = lb2, ub3
+    if ub2 > ub3:
+        bnds_dict[arg2] = lb2, ub3
+        ub2 = ub3
+    if ub1 > ub2:
+        bnds_dict[arg1] = lb1, ub2
+
+
 def _prop_bnds_root_to_leaf_ProductExpression(node, bnds_dict, feasibility_tol):
     """
 
@@ -953,6 +1015,16 @@ def _prop_bnds_root_to_leaf_NamedExpression(node, bnds_dict, feasibility_tol):
 _prop_bnds_root_to_leaf_map[ObjectiveData] = _prop_bnds_root_to_leaf_NamedExpression
 _prop_bnds_root_to_leaf_map[ScalarObjective] = _prop_bnds_root_to_leaf_NamedExpression
 
+_prop_bnds_root_to_leaf_map[relational_expr.EqualityExpression] = (
+    _prop_bnds_root_to_leaf_equality
+)
+_prop_bnds_root_to_leaf_map[relational_expr.InequalityExpression] = (
+    _prop_bnds_root_to_leaf_inequality
+)
+_prop_bnds_root_to_leaf_map[relational_expr.RangedExpression] = (
+    _prop_bnds_root_to_leaf_ranged
+)
+
 
 def _check_and_reset_bounds(var, lb, ub):
     """
@@ -1250,44 +1322,27 @@ def _fbbt_con(con, config):
 
     # a walker to propagate bounds from the variables to the root
     visitorA = _FBBTVisitorLeafToRoot(bnds_dict, feasibility_tol=config.feasibility_tol)
-    visitorA.walk_expression(con.body)
+    visitorA.walk_expression(con.expr)
 
-    # Now we need to replace the bounds in bnds_dict for the root
-    # node with the bounds on the constraint (if those bounds are
-    # better).
-    _lb = value(con.lower)
-    _ub = value(con.upper)
-    if _lb is None:
-        _lb = -interval.inf
-    if _ub is None:
-        _ub = interval.inf
-
-    lb, ub = bnds_dict[con.body]
+    always_feasible, possibly_feasible = bnds_dict[con.expr]
 
     # check if the constraint is infeasible
-    if lb > _ub + config.feasibility_tol or ub < _lb - config.feasibility_tol:
+    if not possibly_feasible:
         raise InfeasibleConstraintException(
             'Detected an infeasible constraint during FBBT: {0}'.format(str(con))
         )
 
     # check if the constraint is always satisfied
-    if config.deactivate_satisfied_constraints:
-        if lb >= _lb - config.feasibility_tol and ub <= _ub + config.feasibility_tol:
-            con.deactivate()
-
-    if _lb > lb:
-        lb = _lb
-    if _ub < ub:
-        ub = _ub
-    bnds_dict[con.body] = (lb, ub)
+    if config.deactivate_satisfied_constraints and always_feasible:
+        con.deactivate()
 
     # Now, propagate bounds back from the root to the variables
     visitorB = _FBBTVisitorRootToLeaf(
         bnds_dict,
         integer_tol=config.integer_tol,
         feasibility_tol=config.feasibility_tol,
     )
-    visitorB.dfs_postorder_stack(con.body)
+    visitorB.dfs_postorder_stack(con.expr)
 
     new_var_bounds = ComponentMap()
     for _node, _bnds in bnds_dict.items():
@@ -1334,7 +1389,7 @@ def _fbbt_block(m, config):
     for c in m.component_data_objects(
         ctype=Constraint, active=True, descend_into=config.descend_into, sort=True
     ):
-        for v in identify_variables(c.body):
+        for v in identify_variables(c.expr):
             if v not in var_to_con_map:
                 var_to_con_map[v] = list()
             if v.lb is None:
@@ -1521,14 +1576,14 @@ def __init__(self, comp):
         if comp.ctype == Constraint:
             if comp.is_indexed():
                 for c in comp.values():
-                    self._vars.update(identify_variables(c.body))
+                    self._vars.update(identify_variables(c.expr))
             else:
-                self._vars.update(identify_variables(comp.body))
+                self._vars.update(identify_variables(comp.expr))
         else:
             for c in comp.component_data_objects(
                 Constraint, descend_into=True, active=True, sort=True
             ):
-                self._vars.update(identify_variables(c.body))
+                self._vars.update(identify_variables(c.expr))
 
     def save_bounds(self):
         bnds = ComponentMap()

pyomo/contrib/fbbt/interval.py

@@ -57,7 +57,7 @@ def BoolFlag(val):
     return _true if val else _false
 
 
-def ineq(xl, xu, yl, yu):
+def ineq(xl, xu, yl, yu, feasibility_tol):
     """Compute the "bounds" on an InequalityExpression
 
     Note this is *not* performing interval arithmetic: we are
@@ -67,17 +67,17 @@ def ineq(xl, xu, yl, yu):
 
     """
     ans = []
-    if yl < xu:
+    if yl < xu - feasibility_tol:
         ans.append(_false)
-    if xl <= yu:
+    if xl <= yu + feasibility_tol:
         ans.append(_true)
     assert ans
     if len(ans) == 1:
         ans.append(ans[0])
     return tuple(ans)
 
 
-def eq(xl, xu, yl, yu):
+def eq(xl, xu, yl, yu, feasibility_tol):
     """Compute the "bounds" on an EqualityExpression
 
     Note this is *not* performing interval arithmetic: we are
@@ -87,17 +87,21 @@ def eq(xl, xu, yl, yu):
 
     """
     ans = []
-    if xl != xu or yl != yu or xl != yl:
+    if (
+        abs(xl - xu) > feasibility_tol
+        or abs(yl - yu) > feasibility_tol
+        or abs(xl - yl) > feasibility_tol
+    ):
         ans.append(_false)
-    if xl <= yu and yl <= xu:
+    if xl <= yu + feasibility_tol and yl <= xu + feasibility_tol:
         ans.append(_true)
     assert ans
     if len(ans) == 1:
         ans.append(ans[0])
     return tuple(ans)
 
 
-def ranged(xl, xu, yl, yu, zl, zu):
+def ranged(xl, xu, yl, yu, zl, zu, feasibility_tol):
     """Compute the "bounds" on a RangedExpression
 
     Note this is *not* performing interval arithmetic: we are
@@ -106,8 +110,8 @@ def ranged(xl, xu, yl, yu, zl, zu):
     `z` and `z`, `y` can be outside the range `x` and `z`, or both.
 
     """
-    lb = ineq(xl, xu, yl, yu)
-    ub = ineq(yl, yu, zl, zu)
+    lb = ineq(xl, xu, yl, yu, feasibility_tol)
+    ub = ineq(yl, yu, zl, zu, feasibility_tol)
     ans = []
     if not lb[0] or not ub[0]:
         ans.append(_false)