Fixing bug in piecewise/factorable with math domain and range error, plus named expressions

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

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)

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)
 
 

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)

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(