diff --git a/pyomo/contrib/fbbt/interval.py b/pyomo/contrib/fbbt/interval.py index 472f2cbe03f..3359dd9c029 100644 --- a/pyomo/contrib/fbbt/interval.py +++ b/pyomo/contrib/fbbt/interval.py @@ -205,6 +205,19 @@ def div(xl, xu, yl, yu, feasibility_tol): return mul(xl, xu, *inv(yl, yu, feasibility_tol)) +def pow_or_inf(x, y): + try: + z = x**y + except OverflowError: + if x > 0: + z = inf + elif x < 0: + z = -inf + else: + raise ValueError(f"Unexpected overflow error: {x}**{y}") + return z + + def power(xl, xu, yl, yu, feasibility_tol): """ Compute bounds on x**y. @@ -213,18 +226,18 @@ def power(xl, xu, yl, yu, feasibility_tol): # If x is always positive, things are simple. We only need to # worry about the sign of y. if yl < 0 < yu: - lb = min(xu**yl, xl**yu) - ub = max(xl**yl, xu**yu) + lb = min(pow_or_inf(xu, yl), pow_or_inf(xl, yu)) + ub = max(pow_or_inf(xl, yl), pow_or_inf(xu, yu)) elif yl >= 0: - lb = min(xl**yl, xl**yu) - ub = max(xu**yl, xu**yu) + lb = min(pow_or_inf(xl, yl), pow_or_inf(xl, yu)) + ub = max(pow_or_inf(xu, yl), pow_or_inf(xu, yu)) else: # yu <= 0: - lb = min(xu**yl, xu**yu) - ub = max(xl**yl, xl**yu) + lb = min(pow_or_inf(xu, yl), pow_or_inf(xu, yu)) + ub = max(pow_or_inf(xl, yl), pow_or_inf(xl, yu)) elif xl == 0: if yl >= 0: - lb = min(xl**yl, xl**yu) - ub = max(xu**yl, xu**yu) + lb = min(pow_or_inf(xl, yl), pow_or_inf(xl, yu)) + ub = max(pow_or_inf(xu, yl), pow_or_inf(xu, yu)) elif yu <= 0: lb, ub = inv( *power(xl, xu, *sub(0, 0, yl, yu), feasibility_tol), feasibility_tol diff --git a/pyomo/contrib/piecewise/transform/factorable.py b/pyomo/contrib/piecewise/transform/factorable.py index c6b7e328838..530cdcf602a 100644 --- a/pyomo/contrib/piecewise/transform/factorable.py +++ b/pyomo/contrib/piecewise/transform/factorable.py @@ -27,6 +27,7 @@ NPV_UnaryFunctionExpression, ) from pyomo.core.base.units_container import _PyomoUnit +import pyomo.environ as pyo from pyomo.repn.util import ExitNodeDispatcher from pyomo.core.base import ( @@ -272,9 +273,12 @@ def _handle_pow(node, data, visitor): def _handle_named_expression(node, data, visitor): assert len(data) == 1 - node.expr = data[0] - visitor.substitution_map[node] = node - return node + res = data[0] + # node.expr = data[0] + visitor.substitution_map[node] = res + # visitor.node_to_var_map[res] = visitor.node_to_var_map[data[0]] + # visitor.degree_map[res] = visitor.degree_map[data[0]] + return res def _handle_negation(node, data, visitor): @@ -386,8 +390,14 @@ def create_aux_var(self, expr): return expr else: x = self.block.x.add() + # initialize from the current expression value, if possible + # try: + # x.set_value(pyo.value(expr, exception=True)) + # except: + # x.set_value(None) self.substitution_map[expr] = x c = self.block.c.add(x == expr) + # c.pprint() # we need to compute bounds on x now because some of the # handlers depend on variable bounds (e.g., division) xl, xu = self._interval_visitor.walk_expression(expr) diff --git a/pyomo/contrib/piecewise/transform/nonlinear_to_pwl.py b/pyomo/contrib/piecewise/transform/nonlinear_to_pwl.py index 6912811eb00..466f577902b 100644 --- a/pyomo/contrib/piecewise/transform/nonlinear_to_pwl.py +++ b/pyomo/contrib/piecewise/transform/nonlinear_to_pwl.py @@ -353,6 +353,7 @@ def __init__(self, expr, expr_vars): def __call__(self, *args): for i, v in enumerate(self.expr_vars): v.value = args[i] + return value(self.expr) diff --git a/pyomo/devel/initialization/bounds/bound_variables.py b/pyomo/devel/initialization/bounds/bound_variables.py index ca6eabbc7f8..abef099abf3 100644 --- a/pyomo/devel/initialization/bounds/bound_variables.py +++ b/pyomo/devel/initialization/bounds/bound_variables.py @@ -23,14 +23,41 @@ def bound_all_nonlinear_variables(m: BlockData, default_bound: float = 1.0e8): """ fbbt(m) for v in get_vars(m): - if v.lb is None or v.lb < -default_bound: - logger.debug( - f'Could not obtain a lower bound for {str(v)} better than {-default_bound}; setting the lower bound to {-default_bound}' - ) - v.setlb(-default_bound) - if v.ub is None or v.ub > default_bound: + # If bounds are equal, treat first + if v.lb == v.ub: + if v.lb == None: + v.setlb(-default_bound) + v.setub(default_bound) + if abs(v.ub) < abs(default_bound): + continue + else: + v.setlb(-default_bound) + v.setub(default_bound) + # If bound is none or outside default bound, set to default + else: + if v.lb is None or v.lb < -default_bound: + logger.debug( + f'Could not obtain a lower bound for {str(v)} better than {-default_bound}; setting the lower bound to {-default_bound}' + ) + v.setlb(-default_bound) + + if v.ub is None or v.ub > default_bound: + logger.debug( + f'Could not obtain an upper bound for {str(v)} better than {default_bound}; setting the upper bound to {default_bound}' + ) + v.setub(default_bound) + # If changing v.lb makes it larger than v.lb, set them equal. + if v.lb > v.ub: logger.debug( - f'Could not obtain an upper bound for {str(v)} better than {default_bound}; setting the upper bound to {default_bound}' + f'Lower bound was set higher than upper bound, which is not allowed.' + 'Setting upper bound equal to lower bound.' ) - v.setub(default_bound) + v.setub(v.lb) fbbt(m) + # Slightly shift the bounds to prevent math domain error in log or exp functions + for v in get_vars(m): + d = v.ub - v.lb + d *= 1e-6 + d = min(d, 1e-6) + v.setlb(v.lb + d) + v.setub(v.ub - d) diff --git a/pyomo/devel/initialization/initialize.py b/pyomo/devel/initialization/initialize.py index 674e6547f38..d29afcaa467 100644 --- a/pyomo/devel/initialization/initialize.py +++ b/pyomo/devel/initialization/initialize.py @@ -40,9 +40,7 @@ def _get_solver(sname, reason): def _setup(nlp): # get all variable bounds, domains, etc. to restore them later orig_vars = get_vars(nlp) - orig_var_data = [ - (v, (v.lower, v.upper, v.domain, v.fixed, v.value)) for v in orig_vars - ] + orig_var_data = [(v, (v.lb, v.ub, v.domain, v.fixed, v.value)) for v in orig_vars] for v, vdata in orig_var_data: if vdata[2].isdiscrete(): raise RuntimeError(