From 4863d1fbf196d607c71ac1114ba75f28a7431536 Mon Sep 17 00:00:00 2001 From: Michael Bynum Date: Fri, 16 Jan 2026 20:50:58 -0700 Subject: [PATCH 1/5] initial draft of factorable programming transformation for pwl approximations --- .../contrib/piecewise/transform/factorable.py | 448 ++++++++++++++++++ 1 file changed, 448 insertions(+) create mode 100644 pyomo/contrib/piecewise/transform/factorable.py diff --git a/pyomo/contrib/piecewise/transform/factorable.py b/pyomo/contrib/piecewise/transform/factorable.py new file mode 100644 index 00000000000..8527702baae --- /dev/null +++ b/pyomo/contrib/piecewise/transform/factorable.py @@ -0,0 +1,448 @@ +from pyomo.core.expr.visitor import StreamBasedExpressionVisitor +from pyomo.common.collections import ComponentMap, ComponentSet +from pyomo.common.numeric_types import native_numeric_types +from pyomo.core.expr.numeric_expr import ( + NegationExpression, + PowExpression, + ProductExpression, + MonomialTermExpression, + DivisionExpression, + SumExpression, + LinearExpression, + UnaryFunctionExpression, +) + +from pyomo.repn.util import ExitNodeDispatcher +from pyomo.core.base import ( + VarData, + ParamData, + ExpressionData, + VarList, + ConstraintList, + Block, + Constraint, + Objective, +) +from pyomo.core.base.var import ScalarVar +from pyomo.contrib.fbbt.fbbt import compute_bounds_on_expr, fbbt +from pyomo.core.base.expression import ScalarExpression +from pyomo.core.base.transformation import Transformation, TransformationFactory +from pyomo.common.modeling import unique_component_name +from pyomo.core.base.component import ActiveComponent +from pyomo.core.base.suffix import Suffix +from pyomo.gdp import Disjunct +from pyomo.repn.linear import LinearRepn, LinearRepnVisitor + + +""" +The purpose of this module/transformation is to convert any nonlinear model +to the following form: + +min/max f(x_i) +s.t. + g_j(x_i)*h_j(x_k) + a_j^T*x >/// 1: + arg1 = visitor.create_aux_var(arg1) + arg1_vars = (arg1,) + arg1_nvars = 1 + arg1_degree = 1 + if arg2_nvars > 1: + arg2 = visitor.create_aux_var(arg2) + arg2_vars = (arg2,) + arg2_nvars = 1 + arg2_degree = 1 + res = arg1 * arg2 + # at this point arg1 should have at most 1 variable + # and arg2 should have at most 1 variable + if arg1_nvars == 0: + visitor.node_to_var_map[res] = arg2_vars + elif arg2_nvars == 0: + visitor.node_to_var_map[res] = arg1_vars + else: + x = arg1_vars[0] + y = arg2_vars[0] + if x is y: + visitor.node_to_var_map[res] = (x,) + else: + visitor.node_to_var_map[res] = (x, y) + if arg1_degree == 0: + visitor.degree_map[res] = arg2_degree + elif arg2_degree == 0: + visitor.degree_map[res] = arg1_degree + else: + visitor.degree_map[res] = -1 + visitor.substitution_map[node] = res + return res + + +def _handle_sum(node, data, visitor): + # remember to separate the linear parts of the root expression + # from the nonlinear parts prior to using this visitor + # otherwise, we will get auxilliary variables that we don't + # strictly need. If this is not done, we don't do anything + # incorrect, we just get extra variables and constraints + # I'm not sure we should even worry about it, because pretty much + # any presolve should remove them + arg_list = [] + new_degree = 0 + vset = ComponentSet() + for arg in data: + arg_vars = visitor.node_to_var_map[arg] + arg_degree = visitor.degree_map[arg] + if arg_degree == -1: + arg = visitor.create_aux_var(arg) + arg_vars = (arg,) + arg_degree = 1 + arg_list.append(arg) + if arg_degree != 0: + new_degree = 1 + vset.update(arg_vars) + res = sum(arg_list) + visitor.node_to_var_map[res] = tuple(vset) + visitor.degree_map[res] = new_degree + visitor.substitution_map[node] = res + return res + + +def _handle_division(node, data, visitor): + """ + This one is a bit tricky. If we encounter both x/z and y/z + at different places in the model, we only want one auxilliary + variable for 1/z. + """ + arg1, arg2 = data + arg1_vars = visitor.node_to_var_map[arg1] + arg2_vars = visitor.node_to_var_map[arg2] + arg1_nvars = len(arg1_vars) + arg2_nvars = len(arg2_vars) + arg1_degree = visitor.degree_map[arg1] + arg2_degree = visitor.degree_map[arg2] + + if arg2_degree == 0: + res = arg1 / arg2 + visitor.node_to_var_map[res] = arg1_vars + visitor.degree_map[res] = arg1_degree + visitor.substitution_map[node] = res + return res + + if arg1_nvars > 1: + arg1 = visitor.create_aux_var(arg1) + arg1_vars = (arg1,) + arg1_nvars = 1 + arg1_degree = 1 + + if arg2_nvars > 1: + arg2 = visitor.create_aux_var(arg2) + arg2_vars = (arg2,) + arg2_nvars = 1 + arg2_degree = 1 + + if "div" not in visitor.substitution_map: + visitor.substitution_map["div"] = ComponentMap() + + # now we need to figure out if we have seen 1/arg2 before + if arg2 in visitor.substitution_map["div"]: + aux = visitor.substitution_map["div"][arg2] + else: + aux = visitor.block.x.add() + visitor.substitution_map["div"][arg2] = aux + """ + we can only create a piecewise linear function of 1/arg2 if arg2 is either + strictly greater than 0 or strictly less than 0 + + otherwise, we do + aux = 1 / arg2 + aux * arg2 = 1 + """ + arg2_lb, arg2_ub = compute_bounds_on_expr(arg2) + if (arg2_lb is not None and arg2_lb > 0) or (arg2_ub is not None and arg2_ub < 0): + c = visitor.block.c.add(aux == 1/arg2) # keep it univariate if we can + else: + c = visitor.block.c.add(aux * arg2 == 1) + fbbt(c) + + arg2 = aux + arg2_vars = (arg2,) + arg2_nvars = 1 + arg2_degree = 1 + res = arg1 * arg2 + # at this point arg1 should have at most 1 variable + # and arg2 should have exactly 1 variable + if arg1_nvars == 0: + visitor.node_to_var_map[res] = arg2_vars + else: + x = arg1_vars[0] + y = arg2_vars[0] + if x is y: + visitor.node_to_var_map[res] = (x,) + else: + visitor.node_to_var_map[res] = (x, y) + if arg1_degree == 0: + visitor.degree_map[res] = arg2_degree + else: + visitor.degree_map[res] = -1 + visitor.substitution_map[node] = res + return res + + +def _handle_pow(node, data, visitor): + # arg1 ** arg2 + # exp(arg2 * log(arg1)) + arg1, arg2 = data + arg1_vars = visitor.node_to_var_map[arg1] + arg2_vars = visitor.node_to_var_map[arg2] + arg1_nvars = len(arg1_vars) + arg2_nvars = len(arg2_vars) + arg1_degree = visitor.degree_map[arg1] + arg2_degree = visitor.degree_map[arg2] + + if arg1_nvars > 1: + arg1 = visitor.create_aux_var(arg1) + arg1_vars = (arg1,) + arg1_nvars = 1 + arg1_degree = 1 + + if arg2_nvars > 1: + arg2 = visitor.create_aux_var(arg2) + arg2_vars = (arg2,) + arg2_nvars = 1 + arg2_degree = 1 + + res = arg1**arg2 + # at this point arg1 should have at most 1 variable + # and arg2 should have at most 1 variable + if arg1_nvars == 0: + visitor.node_to_var_map[res] = arg2_vars + elif arg2_nvars == 0: + visitor.node_to_var_map[res] = arg1_vars + else: + x = arg1_vars[0] + y = arg2_vars[0] + if x is y: + visitor.node_to_var_map[res] = (x,) + else: + visitor.node_to_var_map[res] = (x, y) + if arg1_degree == 0 and arg2_degree == 0: + visitor.degree_map[res] = 0 + else: + visitor.degree_map[res] = -1 + visitor.substitution_map[node] = res + return res + + +def _handle_named_expression(node, data, visitor): + assert len(data) == 1 + res = data[0] + visitor.substitution_map[node] = res + return res + + +def _handle_negation(node, data, visitor): + arg = data[0] + res = -arg + visitor.node_to_var_map[res] = visitor.node_to_var_map[arg] + visitor.degree_map[res] = visitor.degree_map[arg] + visitor.substitution_map[node] = res + return res + + +def _handle_unary(node, data, visitor): + arg = data[0] + arg_vars = visitor.node_to_var_map[arg] + arg_nvars = len(arg_vars) + arg_degree = visitor.degree_map[arg] + + if arg_nvars > 1: + arg = visitor.create_aux_var(arg) + arg_vars = (arg,) + arg_nvars = 1 + arg_degree = 1 + res = node.create_node_with_local_data((arg,)) + visitor.node_to_var_map[res] = arg_vars + if arg_degree == 0: + visitor.degree_map[res] = 0 + else: + visitor.degree_map[res] = -1 + visitor.substitution_map[node] = res + return res + + +handlers = ExitNodeDispatcher() +handlers[VarData] = _handle_var +handlers[ScalarVar] = _handle_var +handlers[ProductExpression] = _handle_product +handlers[SumExpression] = _handle_sum +handlers[DivisionExpression] = _handle_division +handlers[PowExpression] = _handle_pow +handlers[MonomialTermExpression] = _handle_product +handlers[LinearExpression] = _handle_sum +handlers[ParamData] = _handle_param +handlers[ExpressionData] = _handle_named_expression +handlers[ScalarExpression] = _handle_named_expression +handlers[NegationExpression] = _handle_negation +handlers[UnaryFunctionExpression] = _handle_unary +handlers[int] = _handle_float +handlers[float] = _handle_float + + +class _UnivariateNonlinearDecompositionVisitor(StreamBasedExpressionVisitor): + def __init__(self, **kwds): + self.block = kwds.pop('aux_block') + super().__init__(**kwds) + self.node_to_var_map = ComponentMap() + self.degree_map = ComponentMap() # values will be 0 (constant), 1 (linear), or -1 (nonlinear) + + self.substitution_map = ComponentMap() + + self.block.x = VarList() + self.block.c = ConstraintList() + + def initializeWalker(self, expr): + if expr in self.substitution_map: + return False, self.substitution_map[expr] + return True, None + + def beforeChild(self, node, child, child_idx): + if child in self.substitution_map: + return False, self.substitution_map[child] + return True, None + + def exitNode(self, node, data): + nt = type(node) + if nt in handlers: + return handlers[type(node)](node, data, self) + elif nt in native_numeric_types: + handlers[nt] = _handle_float + return _handle_float(node, data, self) + else: + raise NotImplementedError(f'unrecognized expression type: {nt}') + + def create_aux_var(self, expr): + if expr in self.substitution_map: + x = self.substitution_map[expr] + else: + x = self.block.x.add() + self.substitution_map[expr] = x + c = self.block.c.add(x == expr) + # we need to compute bounds on x now because some of the + # handlers depend on variable bounds (e.g., division) + fbbt(c) + return x + + +class UnivariateNonlinearDecompositionTransformation(Transformation): + def __init__(self): + super().__init__() + + def _check_for_unknown_active_components(self, model): + known_ctypes = {Constraint, Objective, Block} + for ctype in model.collect_ctypes(active=True, descend_into=True): + if not issubclass(ctype, ActiveComponent): + continue + if ctype in known_ctypes: + continue + if ctype is Suffix: + continue + raise NotImplementedError( + f'UnivariateNonlinearDecompositionTransformation does not know how to ' + f'handle components with ctype {ctype}' + ) + + def _apply_to(self, model, **kwds): + if kwds: + raise ValueError('UnivariateNonlinearDecompositionTransformation does not take any keyword arguments') + + self._check_for_unknown_active_components(model) + + objectives = list(model.component_data_objects(Objective, active=True, descend_into=True)) + constraints = list(model.component_data_objects(Constraint, active=True, descend_into=True)) + + bname = unique_component_name(model, 'auxiliary') + setattr(model, bname, Block()) + block = getattr(model, bname) + visitor = _UnivariateNonlinearDecompositionVisitor(aux_block=block) + linear_repn_visitor = LinearRepnVisitor(subexpression_cache={}) + + for con in constraints: + lower, body, upper = con.to_bounded_expression() + repn = linear_repn_visitor.walk_expression(body) + if repn.nonlinear is None: + continue + nonlinear_part = repn.nonlinear + repn.nonlinear = None + linear_part = repn.to_expression(linear_repn_visitor) + nonlinear_part = visitor.walk_expression(nonlinear_part) + new_body = linear_part + nonlinear_part + con.set_value((lower, new_body, upper)) + + for obj in objectives: + repn = linear_repn_visitor.walk_expression(obj.expr) + if repn.nonlinear is None: + continue + nonlinear_part = repn.nonlinear + repn.nonlinear = None + linear_part = repn.to_expression(linear_repn_visitor) + nonlinear_part = visitor.walk_expression(nonlinear_part) + obj.expr = linear_part + nonlinear_part From c8223818033ed13cf87fc1f6fea8fd71c00e2e5e Mon Sep 17 00:00:00 2001 From: Michael Bynum Date: Tue, 20 Jan 2026 09:59:21 -0700 Subject: [PATCH 2/5] tests for piecewise univariate nonlinear decomposition --- pyomo/contrib/piecewise/__init__.py | 3 + ...test_univariate_nonlinear_decomposition.py | 57 +++++++++++++++++++ .../contrib/piecewise/transform/factorable.py | 42 ++++++-------- 3 files changed, 77 insertions(+), 25 deletions(-) create mode 100644 pyomo/contrib/piecewise/tests/test_univariate_nonlinear_decomposition.py diff --git a/pyomo/contrib/piecewise/__init__.py b/pyomo/contrib/piecewise/__init__.py index e402628b423..3362fb0122d 100644 --- a/pyomo/contrib/piecewise/__init__.py +++ b/pyomo/contrib/piecewise/__init__.py @@ -37,6 +37,9 @@ DomainPartitioningMethod, NonlinearToPWL, ) +from pyomo.contrib.piecewise.transform.factorable import ( + UnivariateNonlinearDecompositionTransformation, +) from pyomo.contrib.piecewise.transform.nested_inner_repn import ( NestedInnerRepresentationGDPTransformation, ) diff --git a/pyomo/contrib/piecewise/tests/test_univariate_nonlinear_decomposition.py b/pyomo/contrib/piecewise/tests/test_univariate_nonlinear_decomposition.py new file mode 100644 index 00000000000..da247d56335 --- /dev/null +++ b/pyomo/contrib/piecewise/tests/test_univariate_nonlinear_decomposition.py @@ -0,0 +1,57 @@ +from pyomo.common.unittest import TestCase +import pyomo.environ as pyo +from pyomo.contrib import piecewise +from pyomo.core.expr.compare import assertExpressionsEqual + + +pe = pyo + + +def _get_trans(): + return pyo.TransformationFactory('contrib.piecewise.univariate_nonlinear_decomposition') + + +class TestUnivariateNonlinearDecomposition(TestCase): + def test_multiterm(self): + m = pe.ConcreteModel() + m.x = pe.Var() + m.y = pe.Var() + m.z = pe.Var() + m.c = pe.Constraint(expr=m.x + pe.log(m.y + m.z) + 1/pe.exp(m.x**0.5) <= 0) + + trans = _get_trans() + trans.apply_to(m) + aux = m.auxiliary + + assertExpressionsEqual( + self, + m.c.body, + aux.x[3] + aux.x[2] + m.x, + ) + assertExpressionsEqual( + self, + aux.c[1].expr, + aux.x[1] == (m.y + m.z), + ) + assertExpressionsEqual( + self, + aux.c[2].expr, + aux.x[2] == 1/pyo.exp(m.x**0.5), + ) + assertExpressionsEqual( + self, + aux.c[3].expr, + aux.x[3] == pyo.log(aux.x[1]), + ) + self.assertEqual(m.x.lb, 0) + self.assertIsNone(m.x.ub) + self.assertIsNone(m.y.lb) + self.assertIsNone(m.y.ub) + self.assertIsNone(m.z.lb) + self.assertIsNone(m.z.ub) + self.assertEqual(aux.x[1].lb, 0) + self.assertIsNone(aux.x[1].ub) + self.assertEqual(aux.x[2].lb, 0) + self.assertEqual(aux.x[2].ub, 1) + self.assertIsNone(aux.x[3].lb) + self.assertIsNone(aux.x[3].ub) diff --git a/pyomo/contrib/piecewise/transform/factorable.py b/pyomo/contrib/piecewise/transform/factorable.py index 8527702baae..a9c3e337379 100644 --- a/pyomo/contrib/piecewise/transform/factorable.py +++ b/pyomo/contrib/piecewise/transform/factorable.py @@ -34,31 +34,6 @@ from pyomo.repn.linear import LinearRepn, LinearRepnVisitor -""" -The purpose of this module/transformation is to convert any nonlinear model -to the following form: - -min/max f(x_i) -s.t. - g_j(x_i)*h_j(x_k) + a_j^T*x >////// Date: Tue, 20 Jan 2026 12:39:43 -0700 Subject: [PATCH 3/5] testing for univariate nonlinear decomposition --- ...test_univariate_nonlinear_decomposition.py | 213 +++++++++++++++++- .../contrib/piecewise/transform/factorable.py | 51 ++--- 2 files changed, 227 insertions(+), 37 deletions(-) diff --git a/pyomo/contrib/piecewise/tests/test_univariate_nonlinear_decomposition.py b/pyomo/contrib/piecewise/tests/test_univariate_nonlinear_decomposition.py index da247d56335..aa027570313 100644 --- a/pyomo/contrib/piecewise/tests/test_univariate_nonlinear_decomposition.py +++ b/pyomo/contrib/piecewise/tests/test_univariate_nonlinear_decomposition.py @@ -1,7 +1,9 @@ -from pyomo.common.unittest import TestCase +from pyomo.common.unittest import TestCase, skipUnless import pyomo.environ as pyo from pyomo.contrib import piecewise from pyomo.core.expr.compare import assertExpressionsEqual +from pyomo.common.dependencies import numpy_available, numpy +from pyomo.core.expr.numeric_expr import ProductExpression pe = pyo @@ -26,7 +28,7 @@ def test_multiterm(self): assertExpressionsEqual( self, m.c.body, - aux.x[3] + aux.x[2] + m.x, + m.x + aux.x[3] + aux.x[2], ) assertExpressionsEqual( self, @@ -55,3 +57,210 @@ def test_multiterm(self): self.assertEqual(aux.x[2].ub, 1) self.assertIsNone(aux.x[3].lb) self.assertIsNone(aux.x[3].ub) + + def test_common_subexpressions(self): + m = pe.ConcreteModel() + m.x = pe.Var() + m.y = pe.Var() + m.z1 = pe.Var() + m.z2 = pe.Var() + e = -pe.log(m.x + m.y) + m.c1 = pe.Constraint(expr=m.z1 + e == 0) + m.c2 = pe.Constraint(expr=m.z2 + e == 0) + + trans = _get_trans() + trans.apply_to(m) + aux = m.auxiliary + + assertExpressionsEqual( + self, + m.c1.expr, + m.z1 + aux.x[2] == 0, + ) + assertExpressionsEqual( + self, + m.c2.expr, + m.z2 + aux.x[2] == 0, + ) + assertExpressionsEqual( + self, + aux.c[1].expr, + aux.x[1] == m.x + m.y, + ) + assertExpressionsEqual( + self, + aux.c[2].expr, + aux.x[2] == -pe.log(aux.x[1]), + ) + + def test_product_fixed_variable(self): + m = pe.ConcreteModel() + m.x = pe.Var() + m.y = pe.Var() + m.z = pe.Var() + m.c = pe.Constraint(expr=2*pe.log(m.x + m.y) <= 0) + + trans = _get_trans() + trans.apply_to(m) + aux = m.auxiliary + + assertExpressionsEqual( + self, + m.c.expr, + 2*pe.log(aux.x[1]) <= 0, + ) + assertExpressionsEqual( + self, + aux.c[1].expr, + aux.x[1] == m.x + m.y, + ) + + def test_product_variable_fixed(self): + m = pe.ConcreteModel() + m.x = pe.Var() + m.y = pe.Var() + m.z = pe.Var() + m.c = pe.Constraint(expr=pe.log(m.x + m.y)*2 <= 0) + + trans = _get_trans() + trans.apply_to(m) + aux = m.auxiliary + + assertExpressionsEqual( + self, + m.c.expr, + pe.log(aux.x[1])*2 <= 0, + ) + assertExpressionsEqual( + self, + aux.c[1].expr, + aux.x[1] == m.x + m.y, + ) + + def test_prod_sum_sum(self): + m = pe.ConcreteModel() + m.x1 = pe.Var() + m.x2 = pe.Var() + m.x3 = pe.Var() + m.x4 = pe.Var() + m.c = pe.Constraint(expr=(m.x1 + m.x2) * (m.x3 + m.x4) <= 1) + + trans = _get_trans() + trans.apply_to(m) + aux = m.auxiliary + + assertExpressionsEqual( + self, + m.c.expr, + aux.x[1] * aux.x[2] <= 1, + ) + assertExpressionsEqual( + self, + aux.c[1].expr, + aux.x[1] == m.x1 + m.x2, + ) + assertExpressionsEqual( + self, + aux.c[2].expr, + aux.x[2] == m.x3 + m.x4, + ) + + def test_pow_sum_sum(self): + m = pe.ConcreteModel() + m.x1 = pe.Var() + m.x2 = pe.Var() + m.x3 = pe.Var() + m.x4 = pe.Var() + m.c = pe.Constraint(expr=(m.x1 + m.x2) ** (m.x3 + m.x4) <= 1) + + trans = _get_trans() + trans.apply_to(m) + aux = m.auxiliary + + assertExpressionsEqual( + self, + m.c.expr, + aux.x[1] ** aux.x[2] <= 1, + ) + assertExpressionsEqual( + self, + aux.c[1].expr, + aux.x[1] == m.x1 + m.x2, + ) + assertExpressionsEqual( + self, + aux.c[2].expr, + aux.x[2] == m.x3 + m.x4, + ) + + def test_division_var_const(self): + m = pe.ConcreteModel() + m.x = pe.Var() + m.y = pe.Var() + m.c = pe.Constraint(expr=(m.x + m.y) / 2 <= 0) + + trans = _get_trans() + trans.apply_to(m) + aux = m.auxiliary + + assertExpressionsEqual( + self, + m.c.expr, + (m.x + m.y) / 2 <= 0, + ) + + def test_division_sum_sum(self): + m = pe.ConcreteModel() + m.x1 = pe.Var() + m.x2 = pe.Var() + m.x3 = pe.Var() + m.x4 = pe.Var() + m.c = pe.Constraint(expr=(m.x1 + m.x2) / (m.x3 + m.x4) <= 1) + + trans = _get_trans() + trans.apply_to(m) + aux = m.auxiliary + + assertExpressionsEqual( + self, + m.c.expr, + aux.x[1] * aux.x[3] <= 1, + ) + assertExpressionsEqual( + self, + aux.c[1].expr, + aux.x[1] == m.x1 + m.x2, + ) + assertExpressionsEqual( + self, + aux.c[2].expr, + aux.x[2] == m.x3 + m.x4, + ) + assertExpressionsEqual( + self, + aux.c[3].expr, + aux.x[3] * aux.x[2] == 1, + ) + + @skipUnless(numpy_available, "Numpy is not available") + def test_numpy_float(self): + m = pe.ConcreteModel() + m.x = pe.Var() + m.y = pe.Var() + m.z = pe.Var() + m.c = pe.Constraint(expr=ProductExpression((numpy.float64(2.5), pe.log(m.x + m.y))) <= 0) + + trans = _get_trans() + trans.apply_to(m) + aux = m.auxiliary + + assertExpressionsEqual( + self, + m.c.expr, + 2.5*pe.log(aux.x[1]) <= 0, + ) + assertExpressionsEqual( + self, + aux.c[1].expr, + aux.x[1] == m.x + m.y, + ) diff --git a/pyomo/contrib/piecewise/transform/factorable.py b/pyomo/contrib/piecewise/transform/factorable.py index a9c3e337379..3e6af881735 100644 --- a/pyomo/contrib/piecewise/transform/factorable.py +++ b/pyomo/contrib/piecewise/transform/factorable.py @@ -24,31 +24,29 @@ Objective, ) from pyomo.core.base.var import ScalarVar +from pyomo.core.base.param import ScalarParam from pyomo.contrib.fbbt.fbbt import compute_bounds_on_expr, fbbt from pyomo.core.base.expression import ScalarExpression from pyomo.core.base.transformation import Transformation, TransformationFactory from pyomo.common.modeling import unique_component_name from pyomo.core.base.component import ActiveComponent from pyomo.core.base.suffix import Suffix -from pyomo.gdp import Disjunct -from pyomo.repn.linear import LinearRepn, LinearRepnVisitor def _handle_var(node, data, visitor): + if node.fixed: + return _handle_float(node.value, data, visitor) visitor.node_to_var_map[node] = (node,) visitor.degree_map[node] = 1 visitor.substitution_map[node] = node return node -def _handle_float(node, data, visitor): - visitor.node_to_var_map[node] = tuple() - visitor.degree_map[node] = 0 - visitor.substitution_map[node] = node - return node +def _handle_param(node, data, visitor): + return _handle_float(node.value, data, visitor) -def _handle_param(node, data, visitor): +def _handle_float(node, data, visitor): visitor.node_to_var_map[node] = tuple() visitor.degree_map[node] = 0 visitor.substitution_map[node] = node @@ -113,13 +111,6 @@ def _handle_product(node, data, visitor): def _handle_sum(node, data, visitor): - # remember to separate the linear parts of the root expression - # from the nonlinear parts prior to using this visitor - # otherwise, we will get auxilliary variables that we don't - # strictly need. If this is not done, we don't do anything - # incorrect, we just get extra variables and constraints - # I'm not sure we should even worry about it, because pretty much - # any presolve should remove them arg_list = [] new_degree = 0 vset = ComponentSet() @@ -307,13 +298,14 @@ def _handle_unary(node, data, visitor): handlers = ExitNodeDispatcher() handlers[VarData] = _handle_var handlers[ScalarVar] = _handle_var +handlers[ParamData] = _handle_param +handlers[ScalarParam] = _handle_param handlers[ProductExpression] = _handle_product handlers[SumExpression] = _handle_sum handlers[DivisionExpression] = _handle_division handlers[PowExpression] = _handle_pow handlers[MonomialTermExpression] = _handle_product handlers[LinearExpression] = _handle_sum -handlers[ParamData] = _handle_param handlers[ExpressionData] = _handle_named_expression handlers[ScalarExpression] = _handle_named_expression handlers[NegationExpression] = _handle_negation @@ -415,26 +407,15 @@ def _apply_to(self, model, **kwds): setattr(model, bname, Block()) block = getattr(model, bname) visitor = _UnivariateNonlinearDecompositionVisitor(aux_block=block) - linear_repn_visitor = LinearRepnVisitor(subexpression_cache={}) for con in constraints: - lower, body, upper = con.to_bounded_expression() - repn = linear_repn_visitor.walk_expression(body) - if repn.nonlinear is None: - continue - nonlinear_part = repn.nonlinear - repn.nonlinear = None - linear_part = repn.to_expression(linear_repn_visitor) - nonlinear_part = visitor.walk_expression(nonlinear_part) - new_body = linear_part + nonlinear_part - con.set_value((lower, new_body, upper)) + lower, body, upper = con.to_bounded_expression(evaluate_bounds=True) + new_body = visitor.walk_expression(body) + if lower == upper: + con.set_value(new_body == lower) + else: + con.set_value((lower, new_body, upper)) for obj in objectives: - repn = linear_repn_visitor.walk_expression(obj.expr) - if repn.nonlinear is None: - continue - nonlinear_part = repn.nonlinear - repn.nonlinear = None - linear_part = repn.to_expression(linear_repn_visitor) - nonlinear_part = visitor.walk_expression(nonlinear_part) - obj.expr = linear_part + nonlinear_part + new_expr = visitor.walk_expression(obj.expr) + obj.expr = new_expr From ee96e0754662e90f66072285f0d216e25d39c59a Mon Sep 17 00:00:00 2001 From: John Siirola Date: Thu, 5 Mar 2026 21:25:03 -0700 Subject: [PATCH 4/5] NFC: apply black --- ...test_univariate_nonlinear_decomposition.py | 165 ++++-------------- .../contrib/piecewise/transform/factorable.py | 56 +++--- 2 files changed, 67 insertions(+), 154 deletions(-) diff --git a/pyomo/contrib/piecewise/tests/test_univariate_nonlinear_decomposition.py b/pyomo/contrib/piecewise/tests/test_univariate_nonlinear_decomposition.py index aa027570313..fdeb5ae6f05 100644 --- a/pyomo/contrib/piecewise/tests/test_univariate_nonlinear_decomposition.py +++ b/pyomo/contrib/piecewise/tests/test_univariate_nonlinear_decomposition.py @@ -5,12 +5,13 @@ from pyomo.common.dependencies import numpy_available, numpy from pyomo.core.expr.numeric_expr import ProductExpression - pe = pyo def _get_trans(): - return pyo.TransformationFactory('contrib.piecewise.univariate_nonlinear_decomposition') + return pyo.TransformationFactory( + 'contrib.piecewise.univariate_nonlinear_decomposition' + ) class TestUnivariateNonlinearDecomposition(TestCase): @@ -19,32 +20,16 @@ def test_multiterm(self): m.x = pe.Var() m.y = pe.Var() m.z = pe.Var() - m.c = pe.Constraint(expr=m.x + pe.log(m.y + m.z) + 1/pe.exp(m.x**0.5) <= 0) + m.c = pe.Constraint(expr=m.x + pe.log(m.y + m.z) + 1 / pe.exp(m.x**0.5) <= 0) trans = _get_trans() trans.apply_to(m) aux = m.auxiliary - assertExpressionsEqual( - self, - m.c.body, - m.x + aux.x[3] + aux.x[2], - ) - assertExpressionsEqual( - self, - aux.c[1].expr, - aux.x[1] == (m.y + m.z), - ) - assertExpressionsEqual( - self, - aux.c[2].expr, - aux.x[2] == 1/pyo.exp(m.x**0.5), - ) - assertExpressionsEqual( - self, - aux.c[3].expr, - aux.x[3] == pyo.log(aux.x[1]), - ) + assertExpressionsEqual(self, m.c.body, m.x + aux.x[3] + aux.x[2]) + assertExpressionsEqual(self, aux.c[1].expr, aux.x[1] == (m.y + m.z)) + assertExpressionsEqual(self, aux.c[2].expr, aux.x[2] == 1 / pyo.exp(m.x**0.5)) + assertExpressionsEqual(self, aux.c[3].expr, aux.x[3] == pyo.log(aux.x[1])) self.assertEqual(m.x.lb, 0) self.assertIsNone(m.x.ub) self.assertIsNone(m.y.lb) @@ -72,70 +57,38 @@ def test_common_subexpressions(self): trans.apply_to(m) aux = m.auxiliary - assertExpressionsEqual( - self, - m.c1.expr, - m.z1 + aux.x[2] == 0, - ) - assertExpressionsEqual( - self, - m.c2.expr, - m.z2 + aux.x[2] == 0, - ) - assertExpressionsEqual( - self, - aux.c[1].expr, - aux.x[1] == m.x + m.y, - ) - assertExpressionsEqual( - self, - aux.c[2].expr, - aux.x[2] == -pe.log(aux.x[1]), - ) + assertExpressionsEqual(self, m.c1.expr, m.z1 + aux.x[2] == 0) + assertExpressionsEqual(self, m.c2.expr, m.z2 + aux.x[2] == 0) + assertExpressionsEqual(self, aux.c[1].expr, aux.x[1] == m.x + m.y) + assertExpressionsEqual(self, aux.c[2].expr, aux.x[2] == -pe.log(aux.x[1])) def test_product_fixed_variable(self): m = pe.ConcreteModel() m.x = pe.Var() m.y = pe.Var() m.z = pe.Var() - m.c = pe.Constraint(expr=2*pe.log(m.x + m.y) <= 0) + m.c = pe.Constraint(expr=2 * pe.log(m.x + m.y) <= 0) trans = _get_trans() trans.apply_to(m) aux = m.auxiliary - assertExpressionsEqual( - self, - m.c.expr, - 2*pe.log(aux.x[1]) <= 0, - ) - assertExpressionsEqual( - self, - aux.c[1].expr, - aux.x[1] == m.x + m.y, - ) + assertExpressionsEqual(self, m.c.expr, 2 * pe.log(aux.x[1]) <= 0) + assertExpressionsEqual(self, aux.c[1].expr, aux.x[1] == m.x + m.y) def test_product_variable_fixed(self): m = pe.ConcreteModel() m.x = pe.Var() m.y = pe.Var() m.z = pe.Var() - m.c = pe.Constraint(expr=pe.log(m.x + m.y)*2 <= 0) + m.c = pe.Constraint(expr=pe.log(m.x + m.y) * 2 <= 0) trans = _get_trans() trans.apply_to(m) aux = m.auxiliary - assertExpressionsEqual( - self, - m.c.expr, - pe.log(aux.x[1])*2 <= 0, - ) - assertExpressionsEqual( - self, - aux.c[1].expr, - aux.x[1] == m.x + m.y, - ) + assertExpressionsEqual(self, m.c.expr, pe.log(aux.x[1]) * 2 <= 0) + assertExpressionsEqual(self, aux.c[1].expr, aux.x[1] == m.x + m.y) def test_prod_sum_sum(self): m = pe.ConcreteModel() @@ -149,21 +102,9 @@ def test_prod_sum_sum(self): trans.apply_to(m) aux = m.auxiliary - assertExpressionsEqual( - self, - m.c.expr, - aux.x[1] * aux.x[2] <= 1, - ) - assertExpressionsEqual( - self, - aux.c[1].expr, - aux.x[1] == m.x1 + m.x2, - ) - assertExpressionsEqual( - self, - aux.c[2].expr, - aux.x[2] == m.x3 + m.x4, - ) + assertExpressionsEqual(self, m.c.expr, aux.x[1] * aux.x[2] <= 1) + assertExpressionsEqual(self, aux.c[1].expr, aux.x[1] == m.x1 + m.x2) + assertExpressionsEqual(self, aux.c[2].expr, aux.x[2] == m.x3 + m.x4) def test_pow_sum_sum(self): m = pe.ConcreteModel() @@ -177,21 +118,9 @@ def test_pow_sum_sum(self): trans.apply_to(m) aux = m.auxiliary - assertExpressionsEqual( - self, - m.c.expr, - aux.x[1] ** aux.x[2] <= 1, - ) - assertExpressionsEqual( - self, - aux.c[1].expr, - aux.x[1] == m.x1 + m.x2, - ) - assertExpressionsEqual( - self, - aux.c[2].expr, - aux.x[2] == m.x3 + m.x4, - ) + assertExpressionsEqual(self, m.c.expr, aux.x[1] ** aux.x[2] <= 1) + assertExpressionsEqual(self, aux.c[1].expr, aux.x[1] == m.x1 + m.x2) + assertExpressionsEqual(self, aux.c[2].expr, aux.x[2] == m.x3 + m.x4) def test_division_var_const(self): m = pe.ConcreteModel() @@ -203,11 +132,7 @@ def test_division_var_const(self): trans.apply_to(m) aux = m.auxiliary - assertExpressionsEqual( - self, - m.c.expr, - (m.x + m.y) / 2 <= 0, - ) + assertExpressionsEqual(self, m.c.expr, (m.x + m.y) / 2 <= 0) def test_division_sum_sum(self): m = pe.ConcreteModel() @@ -221,26 +146,10 @@ def test_division_sum_sum(self): trans.apply_to(m) aux = m.auxiliary - assertExpressionsEqual( - self, - m.c.expr, - aux.x[1] * aux.x[3] <= 1, - ) - assertExpressionsEqual( - self, - aux.c[1].expr, - aux.x[1] == m.x1 + m.x2, - ) - assertExpressionsEqual( - self, - aux.c[2].expr, - aux.x[2] == m.x3 + m.x4, - ) - assertExpressionsEqual( - self, - aux.c[3].expr, - aux.x[3] * aux.x[2] == 1, - ) + assertExpressionsEqual(self, m.c.expr, aux.x[1] * aux.x[3] <= 1) + assertExpressionsEqual(self, aux.c[1].expr, aux.x[1] == m.x1 + m.x2) + assertExpressionsEqual(self, aux.c[2].expr, aux.x[2] == m.x3 + m.x4) + assertExpressionsEqual(self, aux.c[3].expr, aux.x[3] * aux.x[2] == 1) @skipUnless(numpy_available, "Numpy is not available") def test_numpy_float(self): @@ -248,19 +157,13 @@ def test_numpy_float(self): m.x = pe.Var() m.y = pe.Var() m.z = pe.Var() - m.c = pe.Constraint(expr=ProductExpression((numpy.float64(2.5), pe.log(m.x + m.y))) <= 0) + m.c = pe.Constraint( + expr=ProductExpression((numpy.float64(2.5), pe.log(m.x + m.y))) <= 0 + ) trans = _get_trans() trans.apply_to(m) aux = m.auxiliary - assertExpressionsEqual( - self, - m.c.expr, - 2.5*pe.log(aux.x[1]) <= 0, - ) - assertExpressionsEqual( - self, - aux.c[1].expr, - aux.x[1] == m.x + m.y, - ) + assertExpressionsEqual(self, m.c.expr, 2.5 * pe.log(aux.x[1]) <= 0) + assertExpressionsEqual(self, aux.c[1].expr, aux.x[1] == m.x + m.y) diff --git a/pyomo/contrib/piecewise/transform/factorable.py b/pyomo/contrib/piecewise/transform/factorable.py index 3e6af881735..5704b062a32 100644 --- a/pyomo/contrib/piecewise/transform/factorable.py +++ b/pyomo/contrib/piecewise/transform/factorable.py @@ -14,11 +14,11 @@ from pyomo.repn.util import ExitNodeDispatcher from pyomo.core.base import ( - VarData, - ParamData, - ExpressionData, - VarList, - ConstraintList, + VarData, + ParamData, + ExpressionData, + VarList, + ConstraintList, Block, Constraint, Objective, @@ -68,7 +68,7 @@ def _handle_product(node, data, visitor): visitor.degree_map[res] = arg2_degree visitor.substitution_map[node] = res return res - + if arg2_degree == 0: res = arg1 * arg2 visitor.node_to_var_map[res] = arg1_vars @@ -87,7 +87,7 @@ def _handle_product(node, data, visitor): arg2_nvars = 1 arg2_degree = 1 res = arg1 * arg2 - # at this point arg1 should have at most 1 variable + # at this point arg1 should have at most 1 variable # and arg2 should have at most 1 variable if arg1_nvars == 0: visitor.node_to_var_map[res] = arg2_vars @@ -134,8 +134,8 @@ def _handle_sum(node, data, visitor): def _handle_division(node, data, visitor): """ - This one is a bit tricky. If we encounter both x/z and y/z - at different places in the model, we only want one auxilliary + This one is a bit tricky. If we encounter both x/z and y/z + at different places in the model, we only want one auxilliary variable for 1/z. """ arg1, arg2 = data @@ -183,8 +183,10 @@ def _handle_division(node, data, visitor): aux * arg2 = 1 """ arg2_lb, arg2_ub = compute_bounds_on_expr(arg2) - if (arg2_lb is not None and arg2_lb > 0) or (arg2_ub is not None and arg2_ub < 0): - c = visitor.block.c.add(aux == 1/arg2) # keep it univariate if we can + if (arg2_lb is not None and arg2_lb > 0) or ( + arg2_ub is not None and arg2_ub < 0 + ): + c = visitor.block.c.add(aux == 1 / arg2) # keep it univariate if we can else: c = visitor.block.c.add(aux * arg2 == 1) fbbt(c) @@ -194,7 +196,7 @@ def _handle_division(node, data, visitor): arg2_nvars = 1 arg2_degree = 1 res = arg1 * arg2 - # at this point arg1 should have at most 1 variable + # at this point arg1 should have at most 1 variable # and arg2 should have exactly 1 variable if arg1_nvars == 0: visitor.node_to_var_map[res] = arg2_vars @@ -237,7 +239,7 @@ def _handle_pow(node, data, visitor): arg2_degree = 1 res = arg1**arg2 - # at this point arg1 should have at most 1 variable + # at this point arg1 should have at most 1 variable # and arg2 should have at most 1 variable if arg1_nvars == 0: visitor.node_to_var_map[res] = arg2_vars @@ -319,7 +321,9 @@ def __init__(self, **kwds): self.block = kwds.pop('aux_block') super().__init__(**kwds) self.node_to_var_map = ComponentMap() - self.degree_map = ComponentMap() # values will be 0 (constant), 1 (linear), or -1 (nonlinear) + self.degree_map = ( + ComponentMap() + ) # values will be 0 (constant), 1 (linear), or -1 (nonlinear) self.substitution_map = ComponentMap() @@ -330,12 +334,12 @@ def initializeWalker(self, expr): if expr in self.substitution_map: return False, self.substitution_map[expr] return True, None - + def beforeChild(self, node, child, child_idx): if child in self.substitution_map: return False, self.substitution_map[child] return True, None - + def exitNode(self, node, data): nt = type(node) if nt in handlers: @@ -353,7 +357,7 @@ def create_aux_var(self, expr): x = self.block.x.add() self.substitution_map[expr] = x c = self.block.c.add(x == expr) - # we need to compute bounds on x now because some of the + # we need to compute bounds on x now because some of the # handlers depend on variable bounds (e.g., division) fbbt(c) return x @@ -374,7 +378,7 @@ def create_aux_var(self, expr): By doing so, each nonlinear function is only a function of one or two variables. If this transformation is used prior to the nonlinear_to_pwl transformation, it can significantly reduce the complexity of the PWL approximation. - """ + """, ) class UnivariateNonlinearDecompositionTransformation(Transformation): def __init__(self): @@ -396,13 +400,19 @@ def _check_for_unknown_active_components(self, model): def _apply_to(self, model, **kwds): if kwds: - raise ValueError('UnivariateNonlinearDecompositionTransformation does not take any keyword arguments') - + raise ValueError( + 'UnivariateNonlinearDecompositionTransformation does not take any keyword arguments' + ) + self._check_for_unknown_active_components(model) - objectives = list(model.component_data_objects(Objective, active=True, descend_into=True)) - constraints = list(model.component_data_objects(Constraint, active=True, descend_into=True)) - + objectives = list( + model.component_data_objects(Objective, active=True, descend_into=True) + ) + constraints = list( + model.component_data_objects(Constraint, active=True, descend_into=True) + ) + bname = unique_component_name(model, 'auxiliary') setattr(model, bname, Block()) block = getattr(model, bname) From ce55a10c23979074db67bae060015b9bdc97b708 Mon Sep 17 00:00:00 2001 From: John Siirola Date: Thu, 5 Mar 2026 22:48:49 -0700 Subject: [PATCH 5/5] NFC: fix typo --- pyomo/contrib/piecewise/transform/factorable.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyomo/contrib/piecewise/transform/factorable.py b/pyomo/contrib/piecewise/transform/factorable.py index 5704b062a32..8e2661f42a0 100644 --- a/pyomo/contrib/piecewise/transform/factorable.py +++ b/pyomo/contrib/piecewise/transform/factorable.py @@ -135,7 +135,7 @@ def _handle_sum(node, data, visitor): def _handle_division(node, data, visitor): """ This one is a bit tricky. If we encounter both x/z and y/z - at different places in the model, we only want one auxilliary + at different places in the model, we only want one auxiliary variable for 1/z. """ arg1, arg2 = data