diff --git a/pyomo/contrib/fbbt/expression_bounds_walker.py b/pyomo/contrib/fbbt/expression_bounds_walker.py new file mode 100644 index 00000000000..426d30f0ee6 --- /dev/null +++ b/pyomo/contrib/fbbt/expression_bounds_walker.py @@ -0,0 +1,283 @@ +# ___________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2022 +# National Technology and Engineering Solutions of Sandia, LLC +# Under the terms of Contract DE-NA0003525 with National Technology and +# Engineering Solutions of Sandia, LLC, the U.S. Government retains certain +# rights in this software. +# This software is distributed under the 3-clause BSD License. +# ___________________________________________________________________________ + +from math import pi +from pyomo.common.collections import ComponentMap +from pyomo.contrib.fbbt.interval import ( + add, + acos, + asin, + atan, + cos, + div, + exp, + interval_abs, + log, + log10, + mul, + power, + sin, + sub, + tan, +) +from pyomo.core.base.expression import Expression +from pyomo.core.expr.numeric_expr import ( + NegationExpression, + ProductExpression, + DivisionExpression, + PowExpression, + AbsExpression, + UnaryFunctionExpression, + MonomialTermExpression, + LinearExpression, + SumExpression, + ExternalFunctionExpression, +) +from pyomo.core.expr.numvalue import native_numeric_types, native_types, value +from pyomo.core.expr.visitor import StreamBasedExpressionVisitor +from pyomo.repn.util import BeforeChildDispatcher, ExitNodeDispatcher + +inf = float('inf') + + +class ExpressionBoundsBeforeChildDispatcher(BeforeChildDispatcher): + __slots__ = () + + def __init__(self): + self[ExternalFunctionExpression] = self._before_external_function + + @staticmethod + def _before_external_function(visitor, child): + # [ESJ 10/6/23]: If external functions ever implement callbacks to help with + # this then this should use them + return False, (-inf, inf) + + @staticmethod + def _before_var(visitor, child): + leaf_bounds = visitor.leaf_bounds + if child in leaf_bounds: + pass + elif child.is_fixed() and visitor.use_fixed_var_values_as_bounds: + val = child.value + if val is None: + raise ValueError( + "Var '%s' is fixed to None. This value cannot be used to " + "calculate bounds." % child.name + ) + leaf_bounds[child] = (child.value, child.value) + else: + lb = child.lb + ub = child.ub + if lb is None: + lb = -inf + if ub is None: + ub = inf + leaf_bounds[child] = (lb, ub) + return False, leaf_bounds[child] + + @staticmethod + def _before_named_expression(visitor, child): + leaf_bounds = visitor.leaf_bounds + if child in leaf_bounds: + return False, leaf_bounds[child] + else: + return True, None + + @staticmethod + def _before_param(visitor, child): + return False, (child.value, child.value) + + @staticmethod + def _before_native(visitor, child): + return False, (child, child) + + @staticmethod + def _before_string(visitor, child): + raise ValueError( + f"{child!r} ({type(child)}) is not a valid numeric type. " + f"Cannot compute bounds on expression." + ) + + @staticmethod + def _before_invalid(visitor, child): + raise ValueError( + f"{child!r} ({type(child)}) is not a valid numeric type. " + f"Cannot compute bounds on expression." + ) + + @staticmethod + def _before_complex(visitor, child): + raise ValueError( + f"Cannot compute bounds on expressions containing " + f"complex numbers. Encountered when processing {child}" + ) + + @staticmethod + def _before_npv(visitor, child): + val = value(child) + return False, (val, val) + + +_before_child_handlers = ExpressionBoundsBeforeChildDispatcher() + + +def _handle_ProductExpression(visitor, node, arg1, arg2): + if arg1 is arg2: + return power(*arg1, 2, 2, feasibility_tol=visitor.feasibility_tol) + return mul(*arg1, *arg2) + + +def _handle_SumExpression(visitor, node, *args): + bnds = (0, 0) + for arg in args: + bnds = add(*bnds, *arg) + return bnds + + +def _handle_DivisionExpression(visitor, node, arg1, arg2): + return div(*arg1, *arg2, feasibility_tol=visitor.feasibility_tol) + + +def _handle_PowExpression(visitor, node, arg1, arg2): + return power(*arg1, *arg2, feasibility_tol=visitor.feasibility_tol) + + +def _handle_NegationExpression(visitor, node, arg): + return sub(0, 0, *arg) + + +def _handle_exp(visitor, node, arg): + return exp(*arg) + + +def _handle_log(visitor, node, arg): + return log(*arg) + + +def _handle_log10(visitor, node, arg): + return log10(*arg) + + +def _handle_sin(visitor, node, arg): + return sin(*arg) + + +def _handle_cos(visitor, node, arg): + return cos(*arg) + + +def _handle_tan(visitor, node, arg): + return tan(*arg) + + +def _handle_asin(visitor, node, arg): + return asin(*arg, -pi / 2, pi / 2, visitor.feasibility_tol) + + +def _handle_acos(visitor, node, arg): + return acos(*arg, 0, pi, visitor.feasibility_tol) + + +def _handle_atan(visitor, node, arg): + return atan(*arg, -pi / 2, pi / 2) + + +def _handle_sqrt(visitor, node, arg): + return power(*arg, 0.5, 0.5, feasibility_tol=visitor.feasibility_tol) + + +def _handle_AbsExpression(visitor, node, arg): + return interval_abs(*arg) + + +def _handle_UnaryFunctionExpression(visitor, node, arg): + return _unary_function_dispatcher[node.getname()](visitor, node, arg) + + +def _handle_named_expression(visitor, node, arg): + visitor.leaf_bounds[node] = arg + return arg + + +_unary_function_dispatcher = { + 'exp': _handle_exp, + 'log': _handle_log, + 'log10': _handle_log10, + 'sin': _handle_sin, + 'cos': _handle_cos, + 'tan': _handle_tan, + 'asin': _handle_asin, + 'acos': _handle_acos, + 'atan': _handle_atan, + 'sqrt': _handle_sqrt, +} + + +_operator_dispatcher = ExitNodeDispatcher( + { + ProductExpression: _handle_ProductExpression, + DivisionExpression: _handle_DivisionExpression, + PowExpression: _handle_PowExpression, + AbsExpression: _handle_AbsExpression, + SumExpression: _handle_SumExpression, + MonomialTermExpression: _handle_ProductExpression, + NegationExpression: _handle_NegationExpression, + UnaryFunctionExpression: _handle_UnaryFunctionExpression, + LinearExpression: _handle_SumExpression, + Expression: _handle_named_expression, + } +) + + +class ExpressionBoundsVisitor(StreamBasedExpressionVisitor): + """ + Walker to calculate bounds on an expression, from leaf to root, with + caching of terminal node bounds (Vars and Expressions) + + NOTE: If anything changes on the model (e.g., Var bounds, fixing, mutable + Param values, etc), then you need to either create a new instance of this + walker, or clear self.leaf_bounds! + + Parameters + ---------- + leaf_bounds: ComponentMap in which to cache bounds at leaves of the expression + tree + feasibility_tol: float, feasibility tolerance for interval arithmetic + calculations + use_fixed_var_values_as_bounds: bool, whether or not to use the values of + fixed Vars as the upper and lower bounds for those Vars or to instead + ignore fixed status and use the bounds. Set to 'True' if you do not + anticipate the fixed status of Variables to change for the duration that + the computed bounds should be valid. + """ + + def __init__( + self, + leaf_bounds=None, + feasibility_tol=1e-8, + use_fixed_var_values_as_bounds=False, + ): + super().__init__() + self.leaf_bounds = leaf_bounds if leaf_bounds is not None else ComponentMap() + self.feasibility_tol = feasibility_tol + self.use_fixed_var_values_as_bounds = use_fixed_var_values_as_bounds + + def initializeWalker(self, expr): + walk, result = self.beforeChild(None, expr, 0) + if not walk: + return False, result + return True, expr + + def beforeChild(self, node, child, child_idx): + return _before_child_handlers[child.__class__](self, child) + + def exitNode(self, node, data): + return _operator_dispatcher[node.__class__](self, node, *data) diff --git a/pyomo/contrib/fbbt/fbbt.py b/pyomo/contrib/fbbt/fbbt.py index 5c486488540..db33c27dd96 100644 --- a/pyomo/contrib/fbbt/fbbt.py +++ b/pyomo/contrib/fbbt/fbbt.py @@ -9,9 +9,15 @@ # This software is distributed under the 3-clause BSD License. # ___________________________________________________________________________ +from collections import defaultdict from pyomo.common.collections import ComponentMap, ComponentSet +from pyomo.contrib.fbbt.expression_bounds_walker import ExpressionBoundsVisitor import pyomo.core.expr.numeric_expr as numeric_expr -from pyomo.core.expr.visitor import ExpressionValueVisitor, identify_variables +from pyomo.core.expr.visitor import ( + ExpressionValueVisitor, + identify_variables, + StreamBasedExpressionVisitor, +) from pyomo.core.expr.numvalue import nonpyomo_leaf_types, value from pyomo.core.expr.numvalue import is_fixed import pyomo.contrib.fbbt.interval as interval @@ -24,7 +30,7 @@ import logging from pyomo.common.errors import InfeasibleConstraintException, PyomoException from pyomo.common.config import ( - ConfigBlock, + ConfigDict, ConfigValue, In, NonNegativeFloat, @@ -73,377 +79,274 @@ class FBBTException(PyomoException): pass -def _prop_bnds_leaf_to_root_ProductExpression(node, bnds_dict, feasibility_tol): +def _prop_bnds_leaf_to_root_ProductExpression(visitor, node, arg1, arg2): """ Parameters ---------- + visitor: _FBBTVisitorLeafToRoot node: pyomo.core.expr.numeric_expr.ProductExpression - bnds_dict: ComponentMap - feasibility_tol: float - If the bounds computed on the body of a constraint violate the bounds of the constraint by more than - feasibility_tol, then the constraint is considered infeasible and an exception is raised. This tolerance - is also used when performing certain interval arithmetic operations to ensure that none of the feasible - region is removed due to floating point arithmetic and to prevent math domain errors (a larger value - is more conservative). + arg1: First arg in product expression + arg2: Second arg in product expression """ - assert len(node.args) == 2 - arg1, arg2 = node.args - lb1, ub1 = bnds_dict[arg1] - lb2, ub2 = bnds_dict[arg2] + bnds_dict = visitor.bnds_dict if arg1 is arg2: - bnds_dict[node] = interval.power(lb1, ub1, 2, 2, feasibility_tol) + bnds_dict[node] = interval.power( + *bnds_dict[arg1], 2, 2, visitor.feasibility_tol + ) else: - bnds_dict[node] = interval.mul(lb1, ub1, lb2, ub2) + bnds_dict[node] = interval.mul(*bnds_dict[arg1], *bnds_dict[arg2]) -def _prop_bnds_leaf_to_root_SumExpression(node, bnds_dict, feasibility_tol): +def _prop_bnds_leaf_to_root_SumExpression(visitor, node, *args): """ Parameters ---------- + visitor: _FBBTVisitorLeafToRoot node: pyomo.core.expr.numeric_expr.SumExpression - bnds_dict: ComponentMap - feasibility_tol: float - If the bounds computed on the body of a constraint violate the bounds of the constraint by more than - feasibility_tol, then the constraint is considered infeasible and an exception is raised. This tolerance - is also used when performing certain interval arithmetic operations to ensure that none of the feasible - region is removed due to floating point arithmetic and to prevent math domain errors (a larger value - is more conservative). + args: summands in SumExpression """ + bnds_dict = visitor.bnds_dict bnds = (0, 0) - for arg in node.args: + for arg in args: bnds = interval.add(*bnds, *bnds_dict[arg]) bnds_dict[node] = bnds -def _prop_bnds_leaf_to_root_DivisionExpression(node, bnds_dict, feasibility_tol): +def _prop_bnds_leaf_to_root_DivisionExpression(visitor, node, arg1, arg2): """ Parameters ---------- + visitor: _FBBTVisitorLeafToRoot node: pyomo.core.expr.numeric_expr.DivisionExpression - bnds_dict: ComponentMap - feasibility_tol: float - If the bounds computed on the body of a constraint violate the bounds of the constraint by more than - feasibility_tol, then the constraint is considered infeasible and an exception is raised. This tolerance - is also used when performing certain interval arithmetic operations to ensure that none of the feasible - region is removed due to floating point arithmetic and to prevent math domain errors (a larger value - is more conservative). + arg1: dividend + arg2: divisor """ - assert len(node.args) == 2 - arg1, arg2 = node.args - lb1, ub1 = bnds_dict[arg1] - lb2, ub2 = bnds_dict[arg2] - bnds_dict[node] = interval.div(lb1, ub1, lb2, ub2, feasibility_tol=feasibility_tol) + bnds_dict = visitor.bnds_dict + bnds_dict[node] = interval.div( + *bnds_dict[arg1], *bnds_dict[arg2], feasibility_tol=visitor.feasibility_tol + ) -def _prop_bnds_leaf_to_root_PowExpression(node, bnds_dict, feasibility_tol): +def _prop_bnds_leaf_to_root_PowExpression(visitor, node, arg1, arg2): """ Parameters ---------- + visitor: _FBBTVisitorLeafToRoot node: pyomo.core.expr.numeric_expr.PowExpression - bnds_dict: ComponentMap - feasibility_tol: float - If the bounds computed on the body of a constraint violate the bounds of the constraint by more than - feasibility_tol, then the constraint is considered infeasible and an exception is raised. This tolerance - is also used when performing certain interval arithmetic operations to ensure that none of the feasible - region is removed due to floating point arithmetic and to prevent math domain errors (a larger value - is more conservative). + arg1: base + arg2: exponent """ - assert len(node.args) == 2 - arg1, arg2 = node.args - lb1, ub1 = bnds_dict[arg1] - lb2, ub2 = bnds_dict[arg2] + bnds_dict = visitor.bnds_dict bnds_dict[node] = interval.power( - lb1, ub1, lb2, ub2, feasibility_tol=feasibility_tol + *bnds_dict[arg1], *bnds_dict[arg2], feasibility_tol=visitor.feasibility_tol ) -def _prop_bnds_leaf_to_root_NegationExpression(node, bnds_dict, feasibility_tol): +def _prop_bnds_leaf_to_root_NegationExpression(visitor, node, arg): """ Parameters ---------- - node: pyomo.core.expr.numeric_expr.UnaryFunctionExpression - bnds_dict: ComponentMap - feasibility_tol: float - If the bounds computed on the body of a constraint violate the bounds of the constraint by more than - feasibility_tol, then the constraint is considered infeasible and an exception is raised. This tolerance - is also used when performing certain interval arithmetic operations to ensure that none of the feasible - region is removed due to floating point arithmetic and to prevent math domain errors (a larger value - is more conservative). + visitor: _FBBTVisitorLeafToRoot + node: pyomo.core.expr.numeric_expr.NegationExpression + arg: NegationExpression arg """ - assert len(node.args) == 1 - arg = node.args[0] - lb1, ub1 = bnds_dict[arg] - bnds_dict[node] = interval.sub(0, 0, lb1, ub1) + bnds_dict = visitor.bnds_dict + bnds_dict[node] = interval.sub(0, 0, *bnds_dict[arg]) -def _prop_bnds_leaf_to_root_exp(node, bnds_dict, feasibility_tol): +def _prop_bnds_leaf_to_root_exp(visitor, node, arg): """ Parameters ---------- + visitor: _FBBTVisitorLeafToRoot node: pyomo.core.expr.numeric_expr.UnaryFunctionExpression - bnds_dict: ComponentMap - feasibility_tol: float - If the bounds computed on the body of a constraint violate the bounds of the constraint by more than - feasibility_tol, then the constraint is considered infeasible and an exception is raised. This tolerance - is also used when performing certain interval arithmetic operations to ensure that none of the feasible - region is removed due to floating point arithmetic and to prevent math domain errors (a larger value - is more conservative). + arg: UnaryFunctionExpression arg """ - assert len(node.args) == 1 - arg = node.args[0] - lb1, ub1 = bnds_dict[arg] - bnds_dict[node] = interval.exp(lb1, ub1) + bnds_dict = visitor.bnds_dict + bnds_dict[node] = interval.exp(*bnds_dict[arg]) -def _prop_bnds_leaf_to_root_log(node, bnds_dict, feasibility_tol): +def _prop_bnds_leaf_to_root_log(visitor, node, arg): """ Parameters ---------- + visitor: _FBBTVisitorLeafToRoot node: pyomo.core.expr.numeric_expr.UnaryFunctionExpression - bnds_dict: ComponentMap - feasibility_tol: float - If the bounds computed on the body of a constraint violate the bounds of the constraint by more than - feasibility_tol, then the constraint is considered infeasible and an exception is raised. This tolerance - is also used when performing certain interval arithmetic operations to ensure that none of the feasible - region is removed due to floating point arithmetic and to prevent math domain errors (a larger value - is more conservative). + arg: UnaryFunctionExpression arg """ - assert len(node.args) == 1 - arg = node.args[0] - lb1, ub1 = bnds_dict[arg] - bnds_dict[node] = interval.log(lb1, ub1) + bnds_dict = visitor.bnds_dict + bnds_dict[node] = interval.log(*bnds_dict[arg]) -def _prop_bnds_leaf_to_root_log10(node, bnds_dict, feasibility_tol): +def _prop_bnds_leaf_to_root_log10(visitor, node, arg): """ Parameters ---------- + visitor: _FBBTVisitorLeafToRoot node: pyomo.core.expr.numeric_expr.UnaryFunctionExpression - bnds_dict: ComponentMap - feasibility_tol: float - If the bounds computed on the body of a constraint violate the bounds of the constraint by more than - feasibility_tol, then the constraint is considered infeasible and an exception is raised. This tolerance - is also used when performing certain interval arithmetic operations to ensure that none of the feasible - region is removed due to floating point arithmetic and to prevent math domain errors (a larger value - is more conservative). + arg: UnaryFunctionExpression arg """ - assert len(node.args) == 1 - arg = node.args[0] - lb1, ub1 = bnds_dict[arg] - bnds_dict[node] = interval.log10(lb1, ub1) + bnds_dict = visitor.bnds_dict + bnds_dict[node] = interval.log10(*bnds_dict[arg]) -def _prop_bnds_leaf_to_root_sin(node, bnds_dict, feasibility_tol): +def _prop_bnds_leaf_to_root_sin(visitor, node, arg): """ Parameters ---------- + visitor: _FBBTVisitorLeafToRoot node: pyomo.core.expr.numeric_expr.UnaryFunctionExpression - bnds_dict: ComponentMap - feasibility_tol: float - If the bounds computed on the body of a constraint violate the bounds of the constraint by more than - feasibility_tol, then the constraint is considered infeasible and an exception is raised. This tolerance - is also used when performing certain interval arithmetic operations to ensure that none of the feasible - region is removed due to floating point arithmetic and to prevent math domain errors (a larger value - is more conservative). + arg: UnaryFunctionExpression arg """ - assert len(node.args) == 1 - arg = node.args[0] - lb1, ub1 = bnds_dict[arg] - bnds_dict[node] = interval.sin(lb1, ub1) + bnds_dict = visitor.bnds_dict + bnds_dict[node] = interval.sin(*bnds_dict[arg]) -def _prop_bnds_leaf_to_root_cos(node, bnds_dict, feasibility_tol): +def _prop_bnds_leaf_to_root_cos(visitor, node, arg): """ Parameters ---------- + visitor: _FBBTVisitorLeafToRoot node: pyomo.core.expr.numeric_expr.UnaryFunctionExpression - bnds_dict: ComponentMap - feasibility_tol: float - If the bounds computed on the body of a constraint violate the bounds of the constraint by more than - feasibility_tol, then the constraint is considered infeasible and an exception is raised. This tolerance - is also used when performing certain interval arithmetic operations to ensure that none of the feasible - region is removed due to floating point arithmetic and to prevent math domain errors (a larger value - is more conservative). + arg: UnaryFunctionExpression arg """ - assert len(node.args) == 1 - arg = node.args[0] - lb1, ub1 = bnds_dict[arg] - bnds_dict[node] = interval.cos(lb1, ub1) + bnds_dict = visitor.bnds_dict + bnds_dict[node] = interval.cos(*bnds_dict[arg]) -def _prop_bnds_leaf_to_root_tan(node, bnds_dict, feasibility_tol): +def _prop_bnds_leaf_to_root_tan(visitor, node, arg): """ Parameters ---------- + visitor: _FBBTVisitorLeafToRoot node: pyomo.core.expr.numeric_expr.UnaryFunctionExpression - bnds_dict: ComponentMap - feasibility_tol: float - If the bounds computed on the body of a constraint violate the bounds of the constraint by more than - feasibility_tol, then the constraint is considered infeasible and an exception is raised. This tolerance - is also used when performing certain interval arithmetic operations to ensure that none of the feasible - region is removed due to floating point arithmetic and to prevent math domain errors (a larger value - is more conservative). + arg: UnaryFunctionExpression arg """ - assert len(node.args) == 1 - arg = node.args[0] - lb1, ub1 = bnds_dict[arg] - bnds_dict[node] = interval.tan(lb1, ub1) + bnds_dict = visitor.bnds_dict + bnds_dict[node] = interval.tan(*bnds_dict[arg]) -def _prop_bnds_leaf_to_root_asin(node, bnds_dict, feasibility_tol): +def _prop_bnds_leaf_to_root_asin(visitor, node, arg): """ Parameters ---------- + visitor: _FBBTVisitorLeafToRoot node: pyomo.core.expr.numeric_expr.UnaryFunctionExpression - bnds_dict: ComponentMap - feasibility_tol: float - If the bounds computed on the body of a constraint violate the bounds of the constraint by more than - feasibility_tol, then the constraint is considered infeasible and an exception is raised. This tolerance - is also used when performing certain interval arithmetic operations to ensure that none of the feasible - region is removed due to floating point arithmetic and to prevent math domain errors (a larger value - is more conservative). + arg: UnaryFunctionExpression arg """ - assert len(node.args) == 1 - arg = node.args[0] - lb1, ub1 = bnds_dict[arg] + bnds_dict = visitor.bnds_dict bnds_dict[node] = interval.asin( - lb1, ub1, -interval.inf, interval.inf, feasibility_tol + *bnds_dict[arg], -interval.inf, interval.inf, visitor.feasibility_tol ) -def _prop_bnds_leaf_to_root_acos(node, bnds_dict, feasibility_tol): +def _prop_bnds_leaf_to_root_acos(visitor, node, arg): """ Parameters ---------- + visitor: _FBBTVisitorLeafToRoot node: pyomo.core.expr.numeric_expr.UnaryFunctionExpression - bnds_dict: ComponentMap - feasibility_tol: float - If the bounds computed on the body of a constraint violate the bounds of the constraint by more than - feasibility_tol, then the constraint is considered infeasible and an exception is raised. This tolerance - is also used when performing certain interval arithmetic operations to ensure that none of the feasible - region is removed due to floating point arithmetic and to prevent math domain errors (a larger value - is more conservative). + arg: UnaryFunctionExpression arg """ - assert len(node.args) == 1 - arg = node.args[0] - lb1, ub1 = bnds_dict[arg] + bnds_dict = visitor.bnds_dict bnds_dict[node] = interval.acos( - lb1, ub1, -interval.inf, interval.inf, feasibility_tol + *bnds_dict[arg], -interval.inf, interval.inf, visitor.feasibility_tol ) -def _prop_bnds_leaf_to_root_atan(node, bnds_dict, feasibility_tol): +def _prop_bnds_leaf_to_root_atan(visitor, node, arg): """ Parameters ---------- + visitor: _FBBTVisitorLeafToRoot node: pyomo.core.expr.numeric_expr.UnaryFunctionExpression - bnds_dict: ComponentMap - feasibility_tol: float - If the bounds computed on the body of a constraint violate the bounds of the constraint by more than - feasibility_tol, then the constraint is considered infeasible and an exception is raised. This tolerance - is also used when performing certain interval arithmetic operations to ensure that none of the feasible - region is removed due to floating point arithmetic and to prevent math domain errors (a larger value - is more conservative). """ - assert len(node.args) == 1 - arg = node.args[0] - lb1, ub1 = bnds_dict[arg] - bnds_dict[node] = interval.atan(lb1, ub1, -interval.inf, interval.inf) + bnds_dict = visitor.bnds_dict + bnds_dict[node] = interval.atan(*bnds_dict[arg], -interval.inf, interval.inf) -def _prop_bnds_leaf_to_root_sqrt(node, bnds_dict, feasibility_tol): +def _prop_bnds_leaf_to_root_sqrt(visitor, node, arg): """ Parameters ---------- + visitor: _FBBTVisitorLeafToRoot node: pyomo.core.expr.numeric_expr.UnaryFunctionExpression - bnds_dict: ComponentMap - feasibility_tol: float - If the bounds computed on the body of a constraint violate the bounds of the constraint by more than - feasibility_tol, then the constraint is considered infeasible and an exception is raised. This tolerance - is also used when performing certain interval arithmetic operations to ensure that none of the feasible - region is removed due to floating point arithmetic and to prevent math domain errors (a larger value - is more conservative). + arg: UnaryFunctionExpression arg """ - assert len(node.args) == 1 - arg = node.args[0] - lb1, ub1 = bnds_dict[arg] + bnds_dict = visitor.bnds_dict bnds_dict[node] = interval.power( - lb1, ub1, 0.5, 0.5, feasibility_tol=feasibility_tol + *bnds_dict[arg], 0.5, 0.5, feasibility_tol=visitor.feasibility_tol ) -def _prop_bnds_leaf_to_root_abs(node, bnds_dict, feasibility_tol): - assert len(node.args) == 1 - arg = node.args[0] - lb1, ub1 = bnds_dict[arg] - bnds_dict[node] = interval.interval_abs(lb1, ub1) +def _prop_bnds_leaf_to_root_abs(visitor, node, arg): + bnds_dict = visitor.bnds_dict + bnds_dict[node] = interval.interval_abs(*bnds_dict[arg]) -_unary_leaf_to_root_map = dict() -_unary_leaf_to_root_map['exp'] = _prop_bnds_leaf_to_root_exp -_unary_leaf_to_root_map['log'] = _prop_bnds_leaf_to_root_log -_unary_leaf_to_root_map['log10'] = _prop_bnds_leaf_to_root_log10 -_unary_leaf_to_root_map['sin'] = _prop_bnds_leaf_to_root_sin -_unary_leaf_to_root_map['cos'] = _prop_bnds_leaf_to_root_cos -_unary_leaf_to_root_map['tan'] = _prop_bnds_leaf_to_root_tan -_unary_leaf_to_root_map['asin'] = _prop_bnds_leaf_to_root_asin -_unary_leaf_to_root_map['acos'] = _prop_bnds_leaf_to_root_acos -_unary_leaf_to_root_map['atan'] = _prop_bnds_leaf_to_root_atan -_unary_leaf_to_root_map['sqrt'] = _prop_bnds_leaf_to_root_sqrt -_unary_leaf_to_root_map['abs'] = _prop_bnds_leaf_to_root_abs +def _prop_no_bounds(visitor, node, *args): + visitor.bnds_dict[node] = (-interval.inf, interval.inf) -def _prop_bnds_leaf_to_root_UnaryFunctionExpression(node, bnds_dict, feasibility_tol): +_unary_leaf_to_root_map = defaultdict( + lambda: _prop_no_bounds, + { + 'exp': _prop_bnds_leaf_to_root_exp, + 'log': _prop_bnds_leaf_to_root_log, + 'log10': _prop_bnds_leaf_to_root_log10, + 'sin': _prop_bnds_leaf_to_root_sin, + 'cos': _prop_bnds_leaf_to_root_cos, + 'tan': _prop_bnds_leaf_to_root_tan, + 'asin': _prop_bnds_leaf_to_root_asin, + 'acos': _prop_bnds_leaf_to_root_acos, + 'atan': _prop_bnds_leaf_to_root_atan, + 'sqrt': _prop_bnds_leaf_to_root_sqrt, + 'abs': _prop_bnds_leaf_to_root_abs, + }, +) + + +def _prop_bnds_leaf_to_root_UnaryFunctionExpression(visitor, node, arg): """ Parameters ---------- + visitor: _FBBTVisitorLeafToRoot node: pyomo.core.expr.numeric_expr.UnaryFunctionExpression - bnds_dict: ComponentMap - feasibility_tol: float - If the bounds computed on the body of a constraint violate the bounds of the constraint by more than - feasibility_tol, then the constraint is considered infeasible and an exception is raised. This tolerance - is also used when performing certain interval arithmetic operations to ensure that none of the feasible - region is removed due to floating point arithmetic and to prevent math domain errors (a larger value - is more conservative). + arg: UnaryFunctionExpression arg """ - if node.getname() in _unary_leaf_to_root_map: - _unary_leaf_to_root_map[node.getname()](node, bnds_dict, feasibility_tol) - else: - bnds_dict[node] = (-interval.inf, interval.inf) + _unary_leaf_to_root_map[node.getname()](visitor, node, arg) -def _prop_bnds_leaf_to_root_GeneralExpression(node, bnds_dict, feasibility_tol): +def _prop_bnds_leaf_to_root_GeneralExpression(visitor, node, expr): """ Propagate bounds from children to parent Parameters ---------- + visitor: _FBBTVisitorLeafToRoot node: pyomo.core.base.expression._GeneralExpressionData - bnds_dict: ComponentMap - feasibility_tol: float - If the bounds computed on the body of a constraint violate the bounds of the constraint by more than - feasibility_tol, then the constraint is considered infeasible and an exception is raised. This tolerance - is also used when performing certain interval arithmetic operations to ensure that none of the feasible - region is removed due to floating point arithmetic and to prevent math domain errors (a larger value - is more conservative). + expr: GeneralExpression arg """ - (expr,) = node.args + bnds_dict = visitor.bnds_dict + if node in bnds_dict: + return + if expr.__class__ in native_types: expr_lb = expr_ub = expr else: @@ -451,39 +354,22 @@ def _prop_bnds_leaf_to_root_GeneralExpression(node, bnds_dict, feasibility_tol): bnds_dict[node] = (expr_lb, expr_ub) -_prop_bnds_leaf_to_root_map = dict() -_prop_bnds_leaf_to_root_map[ - numeric_expr.ProductExpression -] = _prop_bnds_leaf_to_root_ProductExpression -_prop_bnds_leaf_to_root_map[ - numeric_expr.DivisionExpression -] = _prop_bnds_leaf_to_root_DivisionExpression -_prop_bnds_leaf_to_root_map[ - numeric_expr.PowExpression -] = _prop_bnds_leaf_to_root_PowExpression -_prop_bnds_leaf_to_root_map[ - numeric_expr.SumExpression -] = _prop_bnds_leaf_to_root_SumExpression -_prop_bnds_leaf_to_root_map[ - numeric_expr.MonomialTermExpression -] = _prop_bnds_leaf_to_root_ProductExpression -_prop_bnds_leaf_to_root_map[ - numeric_expr.NegationExpression -] = _prop_bnds_leaf_to_root_NegationExpression -_prop_bnds_leaf_to_root_map[ - numeric_expr.UnaryFunctionExpression -] = _prop_bnds_leaf_to_root_UnaryFunctionExpression -_prop_bnds_leaf_to_root_map[ - numeric_expr.LinearExpression -] = _prop_bnds_leaf_to_root_SumExpression -_prop_bnds_leaf_to_root_map[numeric_expr.AbsExpression] = _prop_bnds_leaf_to_root_abs - -_prop_bnds_leaf_to_root_map[ - _GeneralExpressionData -] = _prop_bnds_leaf_to_root_GeneralExpression -_prop_bnds_leaf_to_root_map[ - ScalarExpression -] = _prop_bnds_leaf_to_root_GeneralExpression +_prop_bnds_leaf_to_root_map = defaultdict( + lambda: _prop_no_bounds, + { + numeric_expr.ProductExpression: _prop_bnds_leaf_to_root_ProductExpression, + numeric_expr.DivisionExpression: _prop_bnds_leaf_to_root_DivisionExpression, + numeric_expr.PowExpression: _prop_bnds_leaf_to_root_PowExpression, + numeric_expr.SumExpression: _prop_bnds_leaf_to_root_SumExpression, + numeric_expr.MonomialTermExpression: _prop_bnds_leaf_to_root_ProductExpression, + numeric_expr.NegationExpression: _prop_bnds_leaf_to_root_NegationExpression, + 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, + _GeneralExpressionData: _prop_bnds_leaf_to_root_GeneralExpression, + ScalarExpression: _prop_bnds_leaf_to_root_GeneralExpression, + }, +) def _prop_bnds_root_to_leaf_ProductExpression(node, bnds_dict, feasibility_tol): @@ -1071,8 +957,8 @@ def _check_and_reset_bounds(var, lb, ub): """ This function ensures that lb is not less than var.lb and that ub is not greater than var.ub. """ - orig_lb = value(var.lb) - orig_ub = value(var.ub) + orig_lb = var.lb + orig_ub = var.ub if orig_lb is None: orig_lb = -interval.inf if orig_ub is None: @@ -1084,7 +970,77 @@ def _check_and_reset_bounds(var, lb, ub): return lb, ub -class _FBBTVisitorLeafToRoot(ExpressionValueVisitor): +def _before_constant(visitor, child): + if child in visitor.bnds_dict: + pass + else: + visitor.bnds_dict[child] = (child, child) + return False, None + + +def _before_var(visitor, child): + if child in visitor.bnds_dict: + return False, None + elif child.is_fixed() and not visitor.ignore_fixed: + lb = value(child.value) + ub = lb + else: + lb = child.lb + ub = child.ub + if lb is None: + lb = -interval.inf + if ub is None: + ub = interval.inf + if lb - visitor.feasibility_tol > ub: + raise InfeasibleConstraintException( + 'Variable has a lower bound that is larger than its ' + 'upper bound: {0}'.format(str(child)) + ) + visitor.bnds_dict[child] = (lb, ub) + return False, None + + +def _before_NPV(visitor, child): + if child in visitor.bnds_dict: + return False, None + val = value(child) + visitor.bnds_dict[child] = (val, val) + return False, None + + +def _before_other(visitor, child): + return True, None + + +def _before_external_function(visitor, child): + # TODO: provide some mechanism for users to provide interval + # arithmetic callback functions for general external + # functions + visitor.bnds_dict[child] = (-interval.inf, interval.inf) + return False, None + + +def _register_new_before_child_handler(visitor, child): + handlers = _before_child_handlers + child_type = child.__class__ + if child.is_variable_type(): + handlers[child_type] = _before_var + elif not child.is_potentially_variable(): + handlers[child_type] = _before_NPV + else: + handlers[child_type] = _before_other + return handlers[child_type](visitor, child) + + +_before_child_handlers = defaultdict(lambda: _register_new_before_child_handler) +_before_child_handlers[ + numeric_expr.ExternalFunctionExpression +] = _before_external_function +for _type in nonpyomo_leaf_types: + _before_child_handlers[_type] = _before_constant + + +class _FBBTVisitorLeafToRoot(StreamBasedExpressionVisitor): """ This walker propagates bounds from the variables to each node in the expression tree (all the way to the root node). @@ -1099,68 +1055,30 @@ def __init__( bnds_dict: ComponentMap integer_tol: float feasibility_tol: float - If the bounds computed on the body of a constraint violate the bounds of the constraint by more than - feasibility_tol, then the constraint is considered infeasible and an exception is raised. This tolerance - is also used when performing certain interval arithmetic operations to ensure that none of the feasible - region is removed due to floating point arithmetic and to prevent math domain errors (a larger value - is more conservative). + If the bounds computed on the body of a constraint violate the bounds of + the constraint by more than feasibility_tol, then the constraint is + considered infeasible and an exception is raised. This tolerance is also + used when performing certain interval arithmetic operations to ensure that + none of the feasible region is removed due to floating point arithmetic and + to prevent math domain errors (a larger value is more conservative). """ + super().__init__() self.bnds_dict = bnds_dict self.integer_tol = integer_tol self.feasibility_tol = feasibility_tol self.ignore_fixed = ignore_fixed - def visit(self, node, values): - if node.__class__ in _prop_bnds_leaf_to_root_map: - _prop_bnds_leaf_to_root_map[node.__class__]( - node, self.bnds_dict, self.feasibility_tol - ) - else: - self.bnds_dict[node] = (-interval.inf, interval.inf) - return None + def initializeWalker(self, expr): + walk, result = self.beforeChild(None, expr, 0) + if not walk: + return False, result + return True, expr - def visiting_potential_leaf(self, node): - if node.__class__ in nonpyomo_leaf_types: - self.bnds_dict[node] = (node, node) - return True, None + def beforeChild(self, node, child, child_idx): + return _before_child_handlers[child.__class__](self, child) - if node.is_variable_type(): - if node in self.bnds_dict: - return True, None - if node.is_fixed() and not self.ignore_fixed: - lb = value(node.value) - ub = lb - else: - lb = value(node.lb) - ub = value(node.ub) - if lb is None: - lb = -interval.inf - if ub is None: - ub = interval.inf - if lb - self.feasibility_tol > ub: - raise InfeasibleConstraintException( - 'Variable has a lower bound that is larger than its upper bound: {0}'.format( - str(node) - ) - ) - self.bnds_dict[node] = (lb, ub) - return True, None - - if not node.is_potentially_variable(): - # NPV nodes are effectively constant leaves. Evaluate it - # and return the value. - val = value(node) - self.bnds_dict[node] = (val, val) - return True, None - - if node.__class__ is numeric_expr.ExternalFunctionExpression: - # TODO: provide some mechanism for users to provide interval - # arithmetic callback functions for general external - # functions - self.bnds_dict[node] = (-interval.inf, interval.inf) - return True, None - - return False, None + def exitNode(self, node, data): + _prop_bnds_leaf_to_root_map[node.__class__](self, node, *node.args) class _FBBTVisitorRootToLeaf(ExpressionValueVisitor): @@ -1313,7 +1231,7 @@ def _fbbt_con(con, config): ---------- con: pyomo.core.base.constraint.Constraint constraint on which to perform fbbt - config: ConfigBlock + config: ConfigDict see documentation for fbbt Returns @@ -1331,7 +1249,7 @@ 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.dfs_postorder_stack(con.body) + visitorA.walk_expression(con.body) # Now we need to replace the bounds in bnds_dict for the root # node with the bounds on the constraint (if those bounds are @@ -1398,7 +1316,7 @@ def _fbbt_block(m, config): Parameters ---------- m: pyomo.core.base.block.Block or pyomo.core.base.PyomoModel.ConcreteModel - config: ConfigBlock + config: ConfigDict See the docs for fbbt Returns @@ -1421,11 +1339,11 @@ def _fbbt_block(m, config): if v.lb is None: var_lbs[v] = -interval.inf else: - var_lbs[v] = value(v.lb) + var_lbs[v] = v.lb if v.ub is None: var_ubs[v] = interval.inf else: - var_ubs[v] = value(v.ub) + var_ubs[v] = v.ub var_to_con_map[v].append(c) n_cons += 1 @@ -1529,7 +1447,7 @@ def fbbt( A ComponentMap mapping from variables a tuple containing the lower and upper bounds, respectively, computed from FBBT. """ - config = ConfigBlock() + config = ConfigDict() dsc_config = ConfigValue( default=deactivate_satisfied_constraints, domain=In({True, False}) ) @@ -1569,21 +1487,23 @@ def fbbt( def compute_bounds_on_expr(expr, ignore_fixed=False): """ - Compute bounds on an expression based on the bounds on the variables in the expression. + Compute bounds on an expression based on the bounds on the variables in + the expression. Parameters ---------- expr: pyomo.core.expr.numeric_expr.NumericExpression + ignore_fixed: bool, treats fixed Vars as constants if False, else treats + them as Vars Returns ------- lb: float ub: float """ - bnds_dict = ComponentMap() - visitor = _FBBTVisitorLeafToRoot(bnds_dict, ignore_fixed=ignore_fixed) - visitor.dfs_postorder_stack(expr) - lb, ub = bnds_dict[expr] + lb, ub = ExpressionBoundsVisitor( + use_fixed_var_values_as_bounds=not ignore_fixed + ).walk_expression(expr) if lb == -interval.inf: lb = None if ub == interval.inf: diff --git a/pyomo/contrib/fbbt/interval.py b/pyomo/contrib/fbbt/interval.py index aca6531c8df..fd86af4c106 100644 --- a/pyomo/contrib/fbbt/interval.py +++ b/pyomo/contrib/fbbt/interval.py @@ -26,13 +26,15 @@ def sub(xl, xu, yl, yu): def mul(xl, xu, yl, yu): - options = [xl * yl, xl * yu, xu * yl, xu * yu] - if any(math.isnan(i) for i in options): - lb = -inf - ub = inf - else: - lb = min(options) - ub = max(options) + lb = inf + ub = -inf + for i in (xl * yl, xu * yu, xu * yl, xl * yu): + if i < lb: + lb = i + if i > ub: + ub = i + if i != i: # math.isnan(i) + return (-inf, inf) return lb, ub @@ -79,8 +81,7 @@ def inv(xl, xu, feasibility_tol): def div(xl, xu, yl, yu, feasibility_tol): - lb, ub = mul(xl, xu, *inv(yl, yu, feasibility_tol)) - return lb, ub + return mul(xl, xu, *inv(yl, yu, feasibility_tol)) def power(xl, xu, yl, yu, feasibility_tol): @@ -146,7 +147,7 @@ def power(xl, xu, yl, yu, feasibility_tol): else: lb = xl**y ub = xu**y - else: + else: # xu is positive if y < 0: if y % 2 == 0: lb = min(xl**y, xu**y) @@ -154,8 +155,9 @@ def power(xl, xu, yl, yu, feasibility_tol): else: lb = -inf ub = inf - else: + else: # exponent is nonnegative if y % 2 == 0: + # xl is negative and xu is positive, so lb is 0 lb = 0 ub = max(xl**y, xu**y) else: @@ -320,7 +322,7 @@ def _inverse_power2(zl, zu, xl, xu, feasiblity_tol): def interval_abs(xl, xu): abs_xl = abs(xl) abs_xu = abs(xu) - if xl <= 0 <= xu: + if xl <= 0 and 0 <= xu: res_lb = 0 res_ub = max(abs_xl, abs_xu) else: diff --git a/pyomo/contrib/fbbt/tests/test_expression_bounds_walker.py b/pyomo/contrib/fbbt/tests/test_expression_bounds_walker.py new file mode 100644 index 00000000000..c51230155a7 --- /dev/null +++ b/pyomo/contrib/fbbt/tests/test_expression_bounds_walker.py @@ -0,0 +1,305 @@ +# ___________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2022 +# National Technology and Engineering Solutions of Sandia, LLC +# Under the terms of Contract DE-NA0003525 with National Technology and +# Engineering Solutions of Sandia, LLC, the U.S. Government retains certain +# rights in this software. +# This software is distributed under the 3-clause BSD License. +# ___________________________________________________________________________ + +import math +from pyomo.environ import exp, log, log10, sin, cos, tan, asin, acos, atan, sqrt +import pyomo.common.unittest as unittest +from pyomo.contrib.fbbt.expression_bounds_walker import ExpressionBoundsVisitor +from pyomo.core import Any, ConcreteModel, Expression, Param, Var + + +class TestExpressionBoundsWalker(unittest.TestCase): + def make_model(self): + m = ConcreteModel() + m.x = Var(bounds=(-2, 4)) + m.y = Var(bounds=(3, 5)) + m.z = Var(bounds=(0.5, 0.75)) + return m + + def test_sum_bounds(self): + m = self.make_model() + visitor = ExpressionBoundsVisitor() + lb, ub = visitor.walk_expression(m.x + m.y) + self.assertEqual(lb, 1) + self.assertEqual(ub, 9) + + self.assertEqual(len(visitor.leaf_bounds), 2) + self.assertIn(m.x, visitor.leaf_bounds) + self.assertIn(m.y, visitor.leaf_bounds) + self.assertEqual(visitor.leaf_bounds[m.x], (-2, 4)) + self.assertEqual(visitor.leaf_bounds[m.y], (3, 5)) + + def test_fixed_var(self): + m = self.make_model() + m.x.fix(3) + + visitor = ExpressionBoundsVisitor() + lb, ub = visitor.walk_expression(m.x + m.y) + self.assertEqual(lb, 1) + self.assertEqual(ub, 9) + + self.assertEqual(len(visitor.leaf_bounds), 2) + self.assertIn(m.x, visitor.leaf_bounds) + self.assertIn(m.y, visitor.leaf_bounds) + self.assertEqual(visitor.leaf_bounds[m.x], (-2, 4)) + self.assertEqual(visitor.leaf_bounds[m.y], (3, 5)) + + def test_fixed_var_value_used_for_bounds(self): + m = self.make_model() + m.x.fix(3) + + visitor = ExpressionBoundsVisitor(use_fixed_var_values_as_bounds=True) + lb, ub = visitor.walk_expression(m.x + m.y) + self.assertEqual(lb, 6) + self.assertEqual(ub, 8) + + self.assertEqual(len(visitor.leaf_bounds), 2) + self.assertIn(m.x, visitor.leaf_bounds) + self.assertIn(m.y, visitor.leaf_bounds) + self.assertEqual(visitor.leaf_bounds[m.x], (3, 3)) + self.assertEqual(visitor.leaf_bounds[m.y], (3, 5)) + + def test_product_bounds(self): + m = self.make_model() + visitor = ExpressionBoundsVisitor() + lb, ub = visitor.walk_expression(m.x * m.y) + self.assertEqual(lb, -10) + self.assertEqual(ub, 20) + + def test_division_bounds(self): + m = self.make_model() + visitor = ExpressionBoundsVisitor() + lb, ub = visitor.walk_expression(m.x / m.y) + self.assertAlmostEqual(lb, -2 / 3) + self.assertAlmostEqual(ub, 4 / 3) + + def test_power_bounds(self): + m = self.make_model() + visitor = ExpressionBoundsVisitor() + lb, ub = visitor.walk_expression(m.y**m.x) + self.assertEqual(lb, 5 ** (-2)) + self.assertEqual(ub, 5**4) + + def test_sums_of_squares_bounds(self): + m = ConcreteModel() + m.x = Var([1, 2], bounds=(-2, 6)) + visitor = ExpressionBoundsVisitor() + lb, ub = visitor.walk_expression(m.x[1] * m.x[1] + m.x[2] * m.x[2]) + self.assertEqual(lb, 0) + self.assertEqual(ub, 72) + + def test_negation_bounds(self): + m = self.make_model() + visitor = ExpressionBoundsVisitor() + lb, ub = visitor.walk_expression(-(m.y + 3 * m.x)) + self.assertEqual(lb, -17) + self.assertEqual(ub, 3) + + def test_exp_bounds(self): + m = self.make_model() + visitor = ExpressionBoundsVisitor() + lb, ub = visitor.walk_expression(exp(m.y)) + self.assertAlmostEqual(lb, math.e**3) + self.assertAlmostEqual(ub, math.e**5) + + def test_log_bounds(self): + m = self.make_model() + visitor = ExpressionBoundsVisitor() + lb, ub = visitor.walk_expression(log(m.y)) + self.assertAlmostEqual(lb, log(3)) + self.assertAlmostEqual(ub, log(5)) + + def test_log10_bounds(self): + m = self.make_model() + visitor = ExpressionBoundsVisitor() + lb, ub = visitor.walk_expression(log10(m.y)) + self.assertAlmostEqual(lb, log10(3)) + self.assertAlmostEqual(ub, log10(5)) + + def test_sin_bounds(self): + m = self.make_model() + visitor = ExpressionBoundsVisitor() + lb, ub = visitor.walk_expression(sin(m.y)) + self.assertAlmostEqual(lb, -1) # reaches -1 at 3*pi/2 \approx 4.712 + self.assertAlmostEqual(ub, sin(3)) # it's positive here + + def test_cos_bounds(self): + m = self.make_model() + visitor = ExpressionBoundsVisitor() + lb, ub = visitor.walk_expression(cos(m.y)) + self.assertAlmostEqual(lb, -1) # reaches -1 at pi + self.assertAlmostEqual(ub, cos(5)) # it's positive here + + def test_tan_bounds(self): + m = self.make_model() + visitor = ExpressionBoundsVisitor() + lb, ub = visitor.walk_expression(tan(m.y)) + self.assertEqual(lb, -float('inf')) + self.assertEqual(ub, float('inf')) + + def test_asin_bounds(self): + m = self.make_model() + visitor = ExpressionBoundsVisitor() + lb, ub = visitor.walk_expression(asin(m.z)) + self.assertAlmostEqual(lb, asin(0.5)) + self.assertAlmostEqual(ub, asin(0.75)) + + def test_acos_bounds(self): + m = self.make_model() + visitor = ExpressionBoundsVisitor() + lb, ub = visitor.walk_expression(acos(m.z)) + self.assertAlmostEqual(lb, acos(0.75)) + self.assertAlmostEqual(ub, acos(0.5)) + + def test_atan_bounds(self): + m = self.make_model() + visitor = ExpressionBoundsVisitor() + lb, ub = visitor.walk_expression(atan(m.z)) + self.assertAlmostEqual(lb, atan(0.5)) + self.assertAlmostEqual(ub, atan(0.75)) + + def test_sqrt_bounds(self): + m = self.make_model() + visitor = ExpressionBoundsVisitor() + lb, ub = visitor.walk_expression(sqrt(m.y)) + self.assertAlmostEqual(lb, sqrt(3)) + self.assertAlmostEqual(ub, sqrt(5)) + + def test_abs_bounds(self): + m = self.make_model() + visitor = ExpressionBoundsVisitor() + lb, ub = visitor.walk_expression(abs(m.x)) + self.assertEqual(lb, 0) + self.assertEqual(ub, 4) + + def test_leaf_bounds_cached(self): + m = self.make_model() + visitor = ExpressionBoundsVisitor() + lb, ub = visitor.walk_expression(m.x - m.y) + self.assertEqual(lb, -7) + self.assertEqual(ub, 1) + + self.assertIn(m.x, visitor.leaf_bounds) + self.assertEqual(visitor.leaf_bounds[m.x], m.x.bounds) + self.assertIn(m.y, visitor.leaf_bounds) + self.assertEqual(visitor.leaf_bounds[m.y], m.y.bounds) + + # This should exercise the code that uses the cache. + lb, ub = visitor.walk_expression(m.x**2 + 3) + self.assertEqual(lb, 3) + self.assertEqual(ub, 19) + + def test_var_fixed_to_None(self): + m = self.make_model() + m.x.fix(None) + + visitor = ExpressionBoundsVisitor(use_fixed_var_values_as_bounds=True) + with self.assertRaisesRegex( + ValueError, + "Var 'x' is fixed to None. This value cannot be " + "used to calculate bounds.", + ): + lb, ub = visitor.walk_expression(m.x - m.y) + + def test_var_with_no_lb(self): + m = self.make_model() + m.x.setlb(None) + + visitor = ExpressionBoundsVisitor() + lb, ub = visitor.walk_expression(m.x - m.y) + self.assertEqual(lb, -float('inf')) + self.assertEqual(ub, 1) + + def test_var_with_no_ub(self): + m = self.make_model() + m.y.setub(None) + + visitor = ExpressionBoundsVisitor() + lb, ub = visitor.walk_expression(m.x - m.y) + self.assertEqual(lb, -float('inf')) + self.assertEqual(ub, 1) + + def test_param(self): + m = self.make_model() + m.p = Param(initialize=6) + + visitor = ExpressionBoundsVisitor() + lb, ub = visitor.walk_expression(m.p**m.y) + self.assertEqual(lb, 6**3) + self.assertEqual(ub, 6**5) + + def test_mutable_param(self): + m = self.make_model() + m.p = Param(initialize=6, mutable=True) + + visitor = ExpressionBoundsVisitor() + lb, ub = visitor.walk_expression(m.p**m.y) + self.assertEqual(lb, 6**3) + self.assertEqual(ub, 6**5) + + def test_named_expression(self): + m = self.make_model() + m.e = Expression(expr=sqrt(m.x**2 + m.y**2)) + visitor = ExpressionBoundsVisitor() + + lb, ub = visitor.walk_expression(m.e + 4) + self.assertEqual(lb, 7) + self.assertAlmostEqual(ub, sqrt(41) + 4) + + self.assertIn(m.e, visitor.leaf_bounds) + self.assertEqual(visitor.leaf_bounds[m.e][0], 3) + self.assertAlmostEqual(visitor.leaf_bounds[m.e][1], sqrt(41)) + + # exercise the using of the cached bounds + lb, ub = visitor.walk_expression(m.e) + self.assertEqual(lb, 3) + self.assertAlmostEqual(ub, sqrt(41)) + + def test_npv_expression(self): + m = self.make_model() + m.p = Param(initialize=4, mutable=True) + visitor = ExpressionBoundsVisitor() + lb, ub = visitor.walk_expression(1 / m.p) + self.assertEqual(lb, 0.25) + self.assertEqual(ub, 0.25) + + def test_invalid_numeric_type(self): + m = self.make_model() + m.p = Param(initialize=True, domain=Any) + visitor = ExpressionBoundsVisitor() + with self.assertRaisesRegex( + ValueError, + r"True \(\) is not a valid numeric type. " + r"Cannot compute bounds on expression.", + ): + lb, ub = visitor.walk_expression(m.p + m.y) + + def test_invalid_string(self): + m = self.make_model() + m.p = Param(initialize='True', domain=Any) + visitor = ExpressionBoundsVisitor() + with self.assertRaisesRegex( + ValueError, + r"'True' \(\) is not a valid numeric type. " + r"Cannot compute bounds on expression.", + ): + lb, ub = visitor.walk_expression(m.p + m.y) + + def test_invalid_complex(self): + m = self.make_model() + m.p = Param(initialize=complex(4, 5), domain=Any) + visitor = ExpressionBoundsVisitor() + with self.assertRaisesRegex( + ValueError, + r"Cannot compute bounds on expressions containing " + r"complex numbers. Encountered when processing \(4\+5j\)", + ): + lb, ub = visitor.walk_expression(m.p + m.y) diff --git a/pyomo/gdp/plugins/bigm.py b/pyomo/gdp/plugins/bigm.py index c3f65a7ad70..b960b5087ea 100644 --- a/pyomo/gdp/plugins/bigm.py +++ b/pyomo/gdp/plugins/bigm.py @@ -15,6 +15,7 @@ from pyomo.common.collections import ComponentMap from pyomo.common.config import ConfigDict, ConfigValue +from pyomo.common.gc_manager import PauseGC from pyomo.common.modeling import unique_component_name from pyomo.common.deprecation import deprecated, deprecation_warning from pyomo.contrib.cp.transform.logical_to_disjunctive_program import ( @@ -161,6 +162,7 @@ class BigM_Transformation(GDP_to_MIP_Transformation, _BigM_MixIn): def __init__(self): super().__init__(logger) + self._set_up_expr_bound_visitor() def _apply_to(self, instance, **kwds): self.used_args = ComponentMap() # If everything was sure to go well, @@ -169,14 +171,19 @@ def _apply_to(self, instance, **kwds): # as a key in bigMargs, I need the error # not to be when I try to put it into # this map! - try: - self._apply_to_impl(instance, **kwds) - finally: - self._restore_state() - self.used_args.clear() + with PauseGC(): + try: + self._apply_to_impl(instance, **kwds) + finally: + self._restore_state() + self.used_args.clear() + self._expr_bound_visitor.leaf_bounds.clear() + self._expr_bound_visitor.use_fixed_var_values_as_bounds = False def _apply_to_impl(self, instance, **kwds): self._process_arguments(instance, **kwds) + if self._config.assume_fixed_vars_permanent: + self._expr_bound_visitor.use_fixed_var_values_as_bounds = True # filter out inactive targets and handle case where targets aren't # specified. diff --git a/pyomo/gdp/plugins/bigm_mixin.py b/pyomo/gdp/plugins/bigm_mixin.py index ba25dfeffd0..a4df641c8c6 100644 --- a/pyomo/gdp/plugins/bigm_mixin.py +++ b/pyomo/gdp/plugins/bigm_mixin.py @@ -11,7 +11,8 @@ from pyomo.gdp import GDP_Error from pyomo.common.collections import ComponentSet -from pyomo.contrib.fbbt.fbbt import compute_bounds_on_expr +from pyomo.contrib.fbbt.expression_bounds_walker import ExpressionBoundsVisitor +import pyomo.contrib.fbbt.interval as interval from pyomo.core import Suffix @@ -103,6 +104,13 @@ def _get_bigM_arg_list(self, bigm_args, block): block = block.parent_block() return arg_list + def _set_up_expr_bound_visitor(self): + # we assume the default config arg for 'assume_fixed_vars_permanent,` + # and we will change it during apply_to if we need to + self._expr_bound_visitor = ExpressionBoundsVisitor( + use_fixed_var_values_as_bounds=False + ) + def _process_M_value( self, m, @@ -210,10 +218,8 @@ def _get_M_from_args(self, constraint, bigMargs, arg_list, lower, upper): return lower, upper def _estimate_M(self, expr, constraint): - expr_lb, expr_ub = compute_bounds_on_expr( - expr, ignore_fixed=not self._config.assume_fixed_vars_permanent - ) - if expr_lb is None or expr_ub is None: + expr_lb, expr_ub = self._expr_bound_visitor.walk_expression(expr) + if expr_lb == -interval.inf or expr_ub == interval.inf: raise GDP_Error( "Cannot estimate M for unbounded " "expressions.\n\t(found while processing " diff --git a/pyomo/gdp/plugins/multiple_bigm.py b/pyomo/gdp/plugins/multiple_bigm.py index ca0762987a8..18f159c7ca2 100644 --- a/pyomo/gdp/plugins/multiple_bigm.py +++ b/pyomo/gdp/plugins/multiple_bigm.py @@ -14,6 +14,7 @@ from pyomo.common.collections import ComponentMap from pyomo.common.config import ConfigDict, ConfigValue +from pyomo.common.gc_manager import PauseGC from pyomo.common.modeling import unique_component_name from pyomo.core import ( @@ -202,18 +203,24 @@ def __init__(self): super().__init__(logger) self.handlers[Suffix] = self._warn_for_active_suffix self._arg_list = {} + self._set_up_expr_bound_visitor() def _apply_to(self, instance, **kwds): self.used_args = ComponentMap() - try: - self._apply_to_impl(instance, **kwds) - finally: - self._restore_state() - self.used_args.clear() - self._arg_list.clear() + with PauseGC(): + try: + self._apply_to_impl(instance, **kwds) + finally: + self._restore_state() + self.used_args.clear() + self._arg_list.clear() + self._expr_bound_visitor.leaf_bounds.clear() + self._expr_bound_visitor.use_fixed_var_values_as_bounds = False def _apply_to_impl(self, instance, **kwds): self._process_arguments(instance, **kwds) + if self._config.assume_fixed_vars_permanent: + self._bound_visitor.use_fixed_var_values_as_bounds = True if ( self._config.only_mbigm_bound_constraints