From 40222b2f10a9a8f9cc5caa19ed226192e1969828 Mon Sep 17 00:00:00 2001
From: Michael Bynum <mlbynum@sandia.gov>
Date: Thu, 9 Jul 2026 11:08:53 -0600
Subject: [PATCH 01/12] better handling of named expressions in the factorable
 programming visitor

---
 pyomo/contrib/piecewise/transform/factorable.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/pyomo/contrib/piecewise/transform/factorable.py b/pyomo/contrib/piecewise/transform/factorable.py
index c6b7e328838..a978c716e59 100644
--- a/pyomo/contrib/piecewise/transform/factorable.py
+++ b/pyomo/contrib/piecewise/transform/factorable.py
@@ -274,6 +274,8 @@ def _handle_named_expression(node, data, visitor):
     assert len(data) == 1
     node.expr = data[0]
     visitor.substitution_map[node] = node
+    visitor.node_to_var_map[node] = visitor.node_to_var_map[data[0]]
+    visitor.degree_map[node] = visitor.degree_map[data[0]]
     return node
 
 

From 725c756ca05983ed81a41dd4d12ac47a78c0e3e6 Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Thu, 23 Jul 2026 12:12:13 -0600
Subject: [PATCH 02/12] Changes during Michael visit Thursday morning

---
 .../contrib/piecewise/transform/factorable.py  | 18 +++++++++++++-----
 1 file changed, 13 insertions(+), 5 deletions(-)

diff --git a/pyomo/contrib/piecewise/transform/factorable.py b/pyomo/contrib/piecewise/transform/factorable.py
index a978c716e59..7ae8462c71c 100644
--- a/pyomo/contrib/piecewise/transform/factorable.py
+++ b/pyomo/contrib/piecewise/transform/factorable.py
@@ -105,6 +105,11 @@ def _handle_product(node, data, visitor):
     arg1_degree = visitor.degree_map[arg1]
     arg2_degree = visitor.degree_map[arg2]
 
+    print(arg1, arg2)
+    print(arg1_degree, arg2_degree)
+    print(arg1_nvars, arg2_nvars)
+
+
     if arg1_degree == 0:
         res = arg1 * arg2
         visitor.node_to_var_map[res] = arg2_vars
@@ -125,7 +130,9 @@ def _handle_product(node, data, visitor):
         arg1_nvars = 1
         arg1_degree = 1
     if arg2_nvars > 1 or visitor.aggressive_substitution:
+        print('creating aux var for arg2: ', arg2)
         arg2 = visitor.create_aux_var(arg2)
+        print('new arg2: ', arg2)
         arg2_vars = (arg2,)
         arg2_nvars = 1
         arg2_degree = 1
@@ -272,11 +279,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
-    visitor.node_to_var_map[node] = visitor.node_to_var_map[data[0]]
-    visitor.degree_map[node] = visitor.degree_map[data[0]]
-    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):

From 918197ebc84cfb492182e52739208d16b8c62e4f Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Wed, 29 Jul 2026 09:38:07 -0600
Subject: [PATCH 03/12] Attempting to fix math domain issue, in progress

---
 .../contrib/piecewise/transform/factorable.py | 18 +++++++++++++-----
 .../piecewise/transform/nonlinear_to_pwl.py   | 19 +++++++++++++++----
 pyomo/devel/initialization/pwl_init.py        |  6 ++++++
 3 files changed, 34 insertions(+), 9 deletions(-)

diff --git a/pyomo/contrib/piecewise/transform/factorable.py b/pyomo/contrib/piecewise/transform/factorable.py
index 7ae8462c71c..51814cbb66f 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 (
@@ -105,9 +106,9 @@ def _handle_product(node, data, visitor):
     arg1_degree = visitor.degree_map[arg1]
     arg2_degree = visitor.degree_map[arg2]
 
-    print(arg1, arg2)
-    print(arg1_degree, arg2_degree)
-    print(arg1_nvars, arg2_nvars)
+    # print(arg1, arg2)
+    # print(arg1_degree, arg2_degree)
+    # print(arg1_nvars, arg2_nvars)
 
 
     if arg1_degree == 0:
@@ -130,9 +131,9 @@ def _handle_product(node, data, visitor):
         arg1_nvars = 1
         arg1_degree = 1
     if arg2_nvars > 1 or visitor.aggressive_substitution:
-        print('creating aux var for arg2: ', arg2)
+        # print('creating aux var for arg2: ', arg2)
         arg2 = visitor.create_aux_var(arg2)
-        print('new arg2: ', arg2)
+        # print('new arg2: ', arg2)
         arg2_vars = (arg2,)
         arg2_nvars = 1
         arg2_degree = 1
@@ -396,8 +397,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)
@@ -493,6 +500,7 @@ def _apply_to(self, model, **kwds):
         )
 
         for con in constraints:
+            # con.pprint()
             lower, body, upper = con.to_bounded_expression(evaluate_bounds=True)
             new_body = visitor.walk_expression(body)
             if lower is not None and lower == upper:
diff --git a/pyomo/contrib/piecewise/transform/nonlinear_to_pwl.py b/pyomo/contrib/piecewise/transform/nonlinear_to_pwl.py
index 6912811eb00..d11d91f82d4 100644
--- a/pyomo/contrib/piecewise/transform/nonlinear_to_pwl.py
+++ b/pyomo/contrib/piecewise/transform/nonlinear_to_pwl.py
@@ -352,7 +352,13 @@ def __init__(self, expr, expr_vars):
 
     def __call__(self, *args):
         for i, v in enumerate(self.expr_vars):
-            v.value = args[i]
+            # try:
+                v.value = args[i]
+            # except:
+            #     print(self.expr)
+            #     import pdb 
+            #     pdb.set_trace() 
+                
         return value(self.expr)
 
 
@@ -598,9 +604,14 @@ def _transform_constraint(self, cons, config):
         src_data_dict = cons.parent_block().private_data()
         constraints = cons.values() if cons.is_indexed() else (cons,)
         for c in constraints:
-            pw_approx, expr_type = self._approximate_expression(
-                c.body, c, trans_block, config, config.approximate_quadratic_constraints
-            )
+            try:
+                pw_approx, expr_type = self._approximate_expression(
+                    c.body, c, trans_block, config, config.approximate_quadratic_constraints
+                )
+            except:
+                c.pprint()
+                import pdb
+                pdb.set_trace()
 
             if pw_approx is None:
                 # Didn't need approximated, nothing to do
diff --git a/pyomo/devel/initialization/pwl_init.py b/pyomo/devel/initialization/pwl_init.py
index 94f2a19d5b5..5492afcbe96 100644
--- a/pyomo/devel/initialization/pwl_init.py
+++ b/pyomo/devel/initialization/pwl_init.py
@@ -283,11 +283,15 @@ def _initialize_with_piecewise_linear_approximation(
     bound_all_nonlinear_variables(pwl, default_bound=default_bound)
     logger.info('bounded nonlinear variables')
 
+    # pwl.pprint()
+
     # Now, we need to fix variables with equal (or nearly equal) bounds.
     # Otherwise, the PWL transformation complains
     fix_vars_with_equal_bounds(pwl)
     logger.info('fixed variables with equal bounds')
 
+
+
     # now we modify the model by introducing slacks to make sure the PWL
     # approximation is feasible
     # all of the slacks appear linearly, so we don't need to worry about
@@ -295,6 +299,8 @@ def _initialize_with_piecewise_linear_approximation(
     _minimize_infeasibility(pwl)
     logger.info('reformulated model to minimize infeasibility')
 
+    pwl.pprint()
+
     # build the PWL approximation
     trans = pyo.TransformationFactory('contrib.piecewise.nonlinear_to_pwl')
     trans.apply_to(pwl, num_points=2, additively_decompose=False)

From ad58745237a1fc504371932b1f2e59cd34e58a57 Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Thu, 30 Jul 2026 14:26:48 -0600
Subject: [PATCH 04/12] Added correction to shift bounds from 0.

---
 .../piecewise/transform/nonlinear_to_pwl.py   | 21 +++++++++++--------
 .../initialization/bounds/bound_variables.py  |  8 ++++++-
 2 files changed, 19 insertions(+), 10 deletions(-)

diff --git a/pyomo/contrib/piecewise/transform/nonlinear_to_pwl.py b/pyomo/contrib/piecewise/transform/nonlinear_to_pwl.py
index d11d91f82d4..08d63630416 100644
--- a/pyomo/contrib/piecewise/transform/nonlinear_to_pwl.py
+++ b/pyomo/contrib/piecewise/transform/nonlinear_to_pwl.py
@@ -186,10 +186,13 @@ def _get_pwl_function_approximation(func, config, bounds):
     bounds: list of tuples giving upper and lower bounds and a boolean indicating
        if the variable's domain is discrete or not, for each of func's arguments
     """
+    print(func.expr)
+    print(bounds)
+  
     method = config.domain_partitioning_method
     n = config.num_points
     points = _partition_method_dispatcher[method](bounds, n, func, config)
-
+    
     # Don't confuse PiecewiseLinearFunction constructor...
     dim = len(points[0])
     if dim == 1:
@@ -604,14 +607,14 @@ def _transform_constraint(self, cons, config):
         src_data_dict = cons.parent_block().private_data()
         constraints = cons.values() if cons.is_indexed() else (cons,)
         for c in constraints:
-            try:
-                pw_approx, expr_type = self._approximate_expression(
-                    c.body, c, trans_block, config, config.approximate_quadratic_constraints
-                )
-            except:
-                c.pprint()
-                import pdb
-                pdb.set_trace()
+            # try:
+            pw_approx, expr_type = self._approximate_expression(
+                c.body, c, trans_block, config, config.approximate_quadratic_constraints
+            )
+            # except:
+                # c.pprint()
+                # import pdb
+                # pdb.set_trace()
 
             if pw_approx is None:
                 # Didn't need approximated, nothing to do
diff --git a/pyomo/devel/initialization/bounds/bound_variables.py b/pyomo/devel/initialization/bounds/bound_variables.py
index ca6eabbc7f8..d7b798a8b21 100644
--- a/pyomo/devel/initialization/bounds/bound_variables.py
+++ b/pyomo/devel/initialization/bounds/bound_variables.py
@@ -32,5 +32,11 @@ def bound_all_nonlinear_variables(m: BlockData, default_bound: float = 1.0e8):
             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)
+            v.setub(default_bound)    
     fbbt(m)
+    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)

From b09ae644509675939865e29bb946296835e4ded5 Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Wed, 5 Aug 2026 09:08:56 -0600
Subject: [PATCH 05/12] Fix issue with overflow error for power functions in
 fbbt

---
 pyomo/contrib/fbbt/interval.py | 27 +++++++++++++++++++--------
 1 file changed, 19 insertions(+), 8 deletions(-)

diff --git a/pyomo/contrib/fbbt/interval.py b/pyomo/contrib/fbbt/interval.py
index 472f2cbe03f..bd92d334f94 100644
--- a/pyomo/contrib/fbbt/interval.py
+++ b/pyomo/contrib/fbbt/interval.py
@@ -204,6 +204,17 @@ def inv(xl, xu, feasibility_tol):
 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):
     """
@@ -213,18 +224,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

From 9362fb2265538b535454c656dee4c3b1e0516b07 Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Wed, 5 Aug 2026 09:16:49 -0600
Subject: [PATCH 06/12] Changed to fix lb > ub issue and added comments

---
 .../initialization/bounds/bound_variables.py  | 39 +++++++++++++------
 1 file changed, 28 insertions(+), 11 deletions(-)

diff --git a/pyomo/devel/initialization/bounds/bound_variables.py b/pyomo/devel/initialization/bounds/bound_variables.py
index d7b798a8b21..9a323ad6dbf 100644
--- a/pyomo/devel/initialization/bounds/bound_variables.py
+++ b/pyomo/devel/initialization/bounds/bound_variables.py
@@ -23,20 +23,37 @@ 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:
-            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 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)
+
+            elif 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:
+                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)
+        v.setub(v.ub-d)
\ No newline at end of file

From 9c7626feaf7a7e513600613f2e4da6e14c69354e Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Wed, 5 Aug 2026 09:19:01 -0600
Subject: [PATCH 07/12] Ran black

---
 pyomo/contrib/fbbt/interval.py                   | 10 ++++++----
 .../initialization/bounds/bound_variables.py     | 16 ++++++++--------
 pyomo/devel/initialization/pwl_init.py           |  2 --
 3 files changed, 14 insertions(+), 14 deletions(-)

diff --git a/pyomo/contrib/fbbt/interval.py b/pyomo/contrib/fbbt/interval.py
index bd92d334f94..3359dd9c029 100644
--- a/pyomo/contrib/fbbt/interval.py
+++ b/pyomo/contrib/fbbt/interval.py
@@ -204,6 +204,7 @@ def inv(xl, xu, feasibility_tol):
 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
@@ -213,9 +214,10 @@ def pow_or_inf(x, y):
         elif x < 0:
             z = -inf
         else:
-            raise ValueError(f"Unexpected overflow error: {x}**{y}")    
+            raise ValueError(f"Unexpected overflow error: {x}**{y}")
     return z
 
+
 def power(xl, xu, yl, yu, feasibility_tol):
     """
     Compute bounds on x**y.
@@ -225,16 +227,16 @@ def power(xl, xu, yl, yu, feasibility_tol):
         # worry about the sign of y.
         if yl < 0 < 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))
+            ub = max(pow_or_inf(xl, yl), pow_or_inf(xu, yu))
         elif yl >= 0:
-            lb = min(pow_or_inf(xl, yl), pow_or_inf(xl,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(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(pow_or_inf(xl, yl), pow_or_inf(xl,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(
diff --git a/pyomo/devel/initialization/bounds/bound_variables.py b/pyomo/devel/initialization/bounds/bound_variables.py
index 9a323ad6dbf..ad0674b8000 100644
--- a/pyomo/devel/initialization/bounds/bound_variables.py
+++ b/pyomo/devel/initialization/bounds/bound_variables.py
@@ -34,7 +34,7 @@ def bound_all_nonlinear_variables(m: BlockData, default_bound: float = 1.0e8):
                 v.setlb(-default_bound)
                 v.setub(default_bound)
         # If bound is none or outside default bound, set to default
-        else:                              
+        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}'
@@ -42,10 +42,10 @@ def bound_all_nonlinear_variables(m: BlockData, default_bound: float = 1.0e8):
                 v.setlb(-default_bound)
 
             elif 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)
+                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:
                 v.setub(v.lb)
@@ -53,7 +53,7 @@ def bound_all_nonlinear_variables(m: BlockData, default_bound: float = 1.0e8):
     # 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 *= 1e-6
         d = min(d, 1e-6)
-        v.setlb(v.lb+d)
-        v.setub(v.ub-d)
\ No newline at end of file
+        v.setlb(v.lb + d)
+        v.setub(v.ub - d)
diff --git a/pyomo/devel/initialization/pwl_init.py b/pyomo/devel/initialization/pwl_init.py
index 5492afcbe96..32d8e9cdb95 100644
--- a/pyomo/devel/initialization/pwl_init.py
+++ b/pyomo/devel/initialization/pwl_init.py
@@ -290,8 +290,6 @@ def _initialize_with_piecewise_linear_approximation(
     fix_vars_with_equal_bounds(pwl)
     logger.info('fixed variables with equal bounds')
 
-
-
     # now we modify the model by introducing slacks to make sure the PWL
     # approximation is feasible
     # all of the slacks appear linearly, so we don't need to worry about

From 759d7ddd72623078571841a06d0efca8a0f3498e Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Wed, 5 Aug 2026 09:34:55 -0600
Subject: [PATCH 08/12] Removed extra prints, ran black

---
 .../contrib/piecewise/transform/factorable.py |  8 --------
 .../piecewise/transform/nonlinear_to_pwl.py   | 20 ++++++-------------
 pyomo/devel/initialization/pwl_init.py        |  4 ----
 3 files changed, 6 insertions(+), 26 deletions(-)

diff --git a/pyomo/contrib/piecewise/transform/factorable.py b/pyomo/contrib/piecewise/transform/factorable.py
index 51814cbb66f..530cdcf602a 100644
--- a/pyomo/contrib/piecewise/transform/factorable.py
+++ b/pyomo/contrib/piecewise/transform/factorable.py
@@ -106,11 +106,6 @@ def _handle_product(node, data, visitor):
     arg1_degree = visitor.degree_map[arg1]
     arg2_degree = visitor.degree_map[arg2]
 
-    # print(arg1, arg2)
-    # print(arg1_degree, arg2_degree)
-    # print(arg1_nvars, arg2_nvars)
-
-
     if arg1_degree == 0:
         res = arg1 * arg2
         visitor.node_to_var_map[res] = arg2_vars
@@ -131,9 +126,7 @@ def _handle_product(node, data, visitor):
         arg1_nvars = 1
         arg1_degree = 1
     if arg2_nvars > 1 or visitor.aggressive_substitution:
-        # print('creating aux var for arg2: ', arg2)
         arg2 = visitor.create_aux_var(arg2)
-        # print('new arg2: ', arg2)
         arg2_vars = (arg2,)
         arg2_nvars = 1
         arg2_degree = 1
@@ -500,7 +493,6 @@ def _apply_to(self, model, **kwds):
         )
 
         for con in constraints:
-            # con.pprint()
             lower, body, upper = con.to_bounded_expression(evaluate_bounds=True)
             new_body = visitor.walk_expression(body)
             if lower is not None and lower == upper:
diff --git a/pyomo/contrib/piecewise/transform/nonlinear_to_pwl.py b/pyomo/contrib/piecewise/transform/nonlinear_to_pwl.py
index 08d63630416..1d77e4dbf96 100644
--- a/pyomo/contrib/piecewise/transform/nonlinear_to_pwl.py
+++ b/pyomo/contrib/piecewise/transform/nonlinear_to_pwl.py
@@ -186,13 +186,10 @@ def _get_pwl_function_approximation(func, config, bounds):
     bounds: list of tuples giving upper and lower bounds and a boolean indicating
        if the variable's domain is discrete or not, for each of func's arguments
     """
-    print(func.expr)
-    print(bounds)
-  
     method = config.domain_partitioning_method
     n = config.num_points
     points = _partition_method_dispatcher[method](bounds, n, func, config)
-    
+
     # Don't confuse PiecewiseLinearFunction constructor...
     dim = len(points[0])
     if dim == 1:
@@ -355,13 +352,8 @@ def __init__(self, expr, expr_vars):
 
     def __call__(self, *args):
         for i, v in enumerate(self.expr_vars):
-            # try:
-                v.value = args[i]
-            # except:
-            #     print(self.expr)
-            #     import pdb 
-            #     pdb.set_trace() 
-                
+            v.value = args[i]
+
         return value(self.expr)
 
 
@@ -612,9 +604,9 @@ def _transform_constraint(self, cons, config):
                 c.body, c, trans_block, config, config.approximate_quadratic_constraints
             )
             # except:
-                # c.pprint()
-                # import pdb
-                # pdb.set_trace()
+            # c.pprint()
+            # import pdb
+            # pdb.set_trace()
 
             if pw_approx is None:
                 # Didn't need approximated, nothing to do
diff --git a/pyomo/devel/initialization/pwl_init.py b/pyomo/devel/initialization/pwl_init.py
index 32d8e9cdb95..94f2a19d5b5 100644
--- a/pyomo/devel/initialization/pwl_init.py
+++ b/pyomo/devel/initialization/pwl_init.py
@@ -283,8 +283,6 @@ def _initialize_with_piecewise_linear_approximation(
     bound_all_nonlinear_variables(pwl, default_bound=default_bound)
     logger.info('bounded nonlinear variables')
 
-    # pwl.pprint()
-
     # Now, we need to fix variables with equal (or nearly equal) bounds.
     # Otherwise, the PWL transformation complains
     fix_vars_with_equal_bounds(pwl)
@@ -297,8 +295,6 @@ def _initialize_with_piecewise_linear_approximation(
     _minimize_infeasibility(pwl)
     logger.info('reformulated model to minimize infeasibility')
 
-    pwl.pprint()
-
     # build the PWL approximation
     trans = pyo.TransformationFactory('contrib.piecewise.nonlinear_to_pwl')
     trans.apply_to(pwl, num_points=2, additively_decompose=False)

From 68f5ae49053dbb0405db6280e2c5da0f40cf2811 Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Wed, 5 Aug 2026 09:36:05 -0600
Subject: [PATCH 09/12] Removed pdb

---
 pyomo/contrib/piecewise/transform/nonlinear_to_pwl.py | 5 -----
 1 file changed, 5 deletions(-)

diff --git a/pyomo/contrib/piecewise/transform/nonlinear_to_pwl.py b/pyomo/contrib/piecewise/transform/nonlinear_to_pwl.py
index 1d77e4dbf96..466f577902b 100644
--- a/pyomo/contrib/piecewise/transform/nonlinear_to_pwl.py
+++ b/pyomo/contrib/piecewise/transform/nonlinear_to_pwl.py
@@ -599,14 +599,9 @@ def _transform_constraint(self, cons, config):
         src_data_dict = cons.parent_block().private_data()
         constraints = cons.values() if cons.is_indexed() else (cons,)
         for c in constraints:
-            # try:
             pw_approx, expr_type = self._approximate_expression(
                 c.body, c, trans_block, config, config.approximate_quadratic_constraints
             )
-            # except:
-            # c.pprint()
-            # import pdb
-            # pdb.set_trace()
 
             if pw_approx is None:
                 # Didn't need approximated, nothing to do

From d41c7af7f653dc54c9f95774338eb427500ecd68 Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Thu, 6 Aug 2026 13:15:28 -0600
Subject: [PATCH 10/12] Adjust indent and add debug statement

---
 pyomo/devel/initialization/bounds/bound_variables.py | 12 ++++++++----
 1 file changed, 8 insertions(+), 4 deletions(-)

diff --git a/pyomo/devel/initialization/bounds/bound_variables.py b/pyomo/devel/initialization/bounds/bound_variables.py
index ad0674b8000..5216cf7dd38 100644
--- a/pyomo/devel/initialization/bounds/bound_variables.py
+++ b/pyomo/devel/initialization/bounds/bound_variables.py
@@ -41,14 +41,18 @@ def bound_all_nonlinear_variables(m: BlockData, default_bound: float = 1.0e8):
                 )
                 v.setlb(-default_bound)
 
-            elif v.ub is None or v.ub > 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:
-                v.setub(v.lb)
+        # If changing v.lb makes it larger than v.lb, set them equal.
+        if v.lb > v.ub:
+            logger.debug(
+                    f'Lower bound was set higher than upper bound, which is not allowed.' \
+                    'Setting upper bound equal to lower 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):

From a3f68e42ccf8709f0554dbcefb0f2552d2ee27ed Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Fri, 7 Aug 2026 10:54:23 -0600
Subject: [PATCH 11/12] Lower instead of lb

---
 pyomo/devel/initialization/initialize.py | 8 +++++++-
 1 file changed, 7 insertions(+), 1 deletion(-)

diff --git a/pyomo/devel/initialization/initialize.py b/pyomo/devel/initialization/initialize.py
index 7b1117cf5b9..51eac00054e 100644
--- a/pyomo/devel/initialization/initialize.py
+++ b/pyomo/devel/initialization/initialize.py
@@ -41,9 +41,13 @@ 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
+        (v, (v.lb, v.ub, v.domain, v.fixed, v.value)) for v in orig_vars
     ]
+    print(f"Len of orig_vars: {len(orig_vars)}")
     for v, vdata in orig_var_data:
+        print(v.name)
+        print(v.lower)
+        print(v.lb)
         if vdata[2].isdiscrete():
             raise RuntimeError(
                 'Initialization module currently only supports continuous models.'
@@ -54,6 +58,8 @@ def _setup(nlp):
 def _cleanup(orig_var_data):
     # restore variable bounds, domain, etc.
     for v, (lb, ub, domain, fixed, value) in orig_var_data:
+        print(v.name)
+        print(lb)
         v.setlb(lb)
         v.setub(ub)
         v.domain = domain

From e444ee97d9b4e5484507dd6245e70d160aa9f2cb Mon Sep 17 00:00:00 2001
From: Stephen Cini <sscini@sandia.gov>
Date: Fri, 7 Aug 2026 10:57:40 -0600
Subject: [PATCH 12/12] Removed extra prints, ran black

---
 pyomo/devel/initialization/bounds/bound_variables.py |  6 +++---
 pyomo/devel/initialization/initialize.py             | 10 +---------
 2 files changed, 4 insertions(+), 12 deletions(-)

diff --git a/pyomo/devel/initialization/bounds/bound_variables.py b/pyomo/devel/initialization/bounds/bound_variables.py
index 5216cf7dd38..abef099abf3 100644
--- a/pyomo/devel/initialization/bounds/bound_variables.py
+++ b/pyomo/devel/initialization/bounds/bound_variables.py
@@ -49,9 +49,9 @@ def bound_all_nonlinear_variables(m: BlockData, default_bound: float = 1.0e8):
         # If changing v.lb makes it larger than v.lb, set them equal.
         if v.lb > v.ub:
             logger.debug(
-                    f'Lower bound was set higher than upper bound, which is not allowed.' \
-                    'Setting upper bound equal to lower bound.'
-                            )
+                f'Lower bound was set higher than upper bound, which is not allowed.'
+                'Setting upper bound equal to lower bound.'
+            )
             v.setub(v.lb)
     fbbt(m)
     # Slightly shift the bounds to prevent math domain error in log or exp functions
diff --git a/pyomo/devel/initialization/initialize.py b/pyomo/devel/initialization/initialize.py
index 51eac00054e..76f6d556e4c 100644
--- a/pyomo/devel/initialization/initialize.py
+++ b/pyomo/devel/initialization/initialize.py
@@ -40,14 +40,8 @@ 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.lb, v.ub, v.domain, v.fixed, v.value)) for v in orig_vars
-    ]
-    print(f"Len of orig_vars: {len(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:
-        print(v.name)
-        print(v.lower)
-        print(v.lb)
         if vdata[2].isdiscrete():
             raise RuntimeError(
                 'Initialization module currently only supports continuous models.'
@@ -58,8 +52,6 @@ def _setup(nlp):
 def _cleanup(orig_var_data):
     # restore variable bounds, domain, etc.
     for v, (lb, ub, domain, fixed, value) in orig_var_data:
-        print(v.name)
-        print(lb)
         v.setlb(lb)
         v.setub(ub)
         v.domain = domain