diff --git a/pyomo/contrib/piecewise/piecewise_linear_function.py b/pyomo/contrib/piecewise/piecewise_linear_function.py index df0bbd3bf93..80fd658989a 100644 --- a/pyomo/contrib/piecewise/piecewise_linear_function.py +++ b/pyomo/contrib/piecewise/piecewise_linear_function.py @@ -256,6 +256,10 @@ class PiecewiseLinearFunction(Block): as specified, trusting the user in the case of AssumeValid. When no argument or None is passed, the default is Triangulation.Unknown + convex (optional): If True, serves as a user-provided guarantee that the + piecewise-linear function is convex. If False, serves as a user-provided + guarantee that the piecewise-linear function is concave. If None (the + default), the function may be concave or convex. """ _ComponentDataClass = PiecewiseLinearFunctionData @@ -283,6 +287,7 @@ def __init__(self, *args, **kwargs): _tabular_data_arg = kwargs.pop('tabular_data', None) _tabular_data_rule_arg = kwargs.pop('tabular_data_rule', None) _triangulation_rule_arg = kwargs.pop('triangulation', None) + _convex_rule_arg = kwargs.pop('convex', None) kwargs.setdefault('ctype', PiecewiseLinearFunction) Block.__init__(self, *args, **kwargs) @@ -304,6 +309,9 @@ def __init__(self, *args, **kwargs): self._triangulation_rule = Initializer( _triangulation_rule_arg, treat_sequences_as_mappings=False ) + self._convex_rule = Initializer( + _convex_rule_arg, treat_sequences_as_mappings=False + ) def _get_dimension_from_points(self, points): if len(points) < 1: @@ -606,6 +614,11 @@ def _getitem_when_not_present(self, index): ) obj = handler(self, obj, parent, nonlinear_function) + # Update convexity info + obj.convex = None + if self._convex_rule is not None: + obj.convex = self._convex_rule(parent, index) + return obj diff --git a/pyomo/contrib/piecewise/tests/test_epigraph_hypograph.py b/pyomo/contrib/piecewise/tests/test_epigraph_hypograph.py new file mode 100644 index 00000000000..d37a056897b --- /dev/null +++ b/pyomo/contrib/piecewise/tests/test_epigraph_hypograph.py @@ -0,0 +1,335 @@ +# ____________________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2026 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 pyomo.common.unittest as unittest +from pyomo.contrib.piecewise import PiecewiseLinearFunction +from pyomo.contrib.piecewise.transform.epigraph_hypograph import ( + PWLToEpigraphOrHypograph, +) +from pyomo.core.expr.compare import assertExpressionsEqual +from pyomo.environ import ( + ConcreteModel, + Constraint, + Set, + Var, + Objective, + TransformationFactory, + SolverFactory, + maximize, + minimize, + value, +) +from pyomo.contrib.piecewise.transform.piecewise_linear_transformation_base import ( + PiecewiseLinearTransformationBase, +) + +gurobi_available = ( + SolverFactory('gurobi').available(exception_flag=False) + and SolverFactory('gurobi').license_is_valid() +) + + +def x_squared(x): + return x**2 + + +def neg_x_squared(x): + return -(x**2) + + +class TestEpigraphHypographTransformation(unittest.TestCase): + """Test the epigraph/hypograph transformation for piecewise linear functions.""" + + def pwl_x_squared(self): + # Create a model with a convex piecewise linear function + m = ConcreteModel() + m.x = Var(bounds=(0, 4)) + + # Define breakpoints for x in [0, 4] + breakpoints = [0, 1, 2, 3, 4] + + # Create piecewise linear approximation of x² + m.pw = PiecewiseLinearFunction( + points=breakpoints, + function=x_squared, + convex=True, # Specify that this is convex + ) + + m.c = Constraint(expr=m.x >= 1.5) + + # Use the piecewise function in an objective - minimize to test the piecewise part + m.obj = Objective(expr=m.pw(m.x), sense=minimize) + + return m + + def test_convex_pwl_minimization(self): + """Test transformation of a convex piecewise linear function (x² approximation).""" + m = self.pwl_x_squared() + # Apply the epigraph transformation + transformation = TransformationFactory('contrib.piecewise.epigraph_hypograph') + transformation.apply_to(m) + + # Verify the transformation created the expected components + # The transformation should have created: + # 1. A substitute variable + # 2. Epigraphical constraints (one for each linear segment) + + # Get the substitute var and the transformed function block + sub_var = m.pw.get_transformation_var(m.obj.expr) + self.assertIsInstance(sub_var, Var) + pw_block = sub_var.parent_block() + + # Check that epigraphical constraints exist + epigraphical_constraints = pw_block.component('epigraphical_constraints') + self.assertIsInstance(epigraphical_constraints, Constraint) + + # For 5 breakpoints, we have 4 segments, so we should have 4 constraints + self.assertEqual(len(epigraphical_constraints), 4) + + # Verify the objective now uses the substitute variable + # The original pw(m.x) should be replaced with substitute_var + obj_expr = m.obj.expr.expr + self.assertIs(obj_expr, sub_var) + + exprs = [ + m.x - sub_var, + 3.0 * m.x - 2.0 - sub_var, + 5.0 * m.x - 6.0 - sub_var, + 7.0 * m.x - 12.0 - sub_var, + ] + + # Verify each constraint is of the form: linear_func_expr <= substitute_var + # (since this is a convex function, we use epigraph) + for idx, cons in epigraphical_constraints.items(): + constraint = pw_block.epigraphical_constraints[idx] + # For epigraph of convex function, we expect <= constraints + # The constraint is stored as: linear_func_expr <= substitute_var + # which Pyomo represents as: linear_func_expr - substitute_var <= 0 + # So the upper bound should be 0 + self.assertEqual(constraint.upper, 0) + self.assertIsNone(constraint.lower) + print(constraint.body) + assertExpressionsEqual(self, constraint.body, exprs[idx]) + + @unittest.skipUnless(gurobi_available, "Gurobi is not available") + def test_solve_pwl_minimization(self): + m = self.pwl_x_squared() + # Apply the epigraph transformation + transformation = TransformationFactory('contrib.piecewise.epigraph_hypograph') + transformation.apply_to(m) + + # Now solve the model to verify it produces the correct solution + # The minimum of x² for x >= 1.5 is at x = 1.5, with value 2.25 + solver = SolverFactory('gurobi') + results = solver.solve(m) + + # Check that the solve was successful + self.assertTrue( + results.solver.termination_condition == 'optimal' + or results.solver.termination_condition == 'feasible' + ) + + # Check the solution + # The optimal x should be 1.5 (the lower bound) + self.assertAlmostEqual(m.x.value, 1.5, places=4) + + # The objective value should be approximately 2.25 (1.5²) + expected_obj = 3 * 1.5 - 2 + self.assertAlmostEqual(value(m.obj.expr), expected_obj, places=4) + + def pwl_neg_x_squared(self): + # Create a model with a concave piecewise linear function + m = ConcreteModel() + m.x = Var(bounds=(0, 4)) + + # Define breakpoints for x in [0, 4] + breakpoints = [0, 1, 2, 3, 4] + + # Create piecewise linear approximation of -x² + m.pw = PiecewiseLinearFunction( + points=breakpoints, + function=neg_x_squared, + convex=False, # Specify that this is concave + ) + + m.c = Constraint(expr=m.x <= 2.5) + + # Use the piecewise function in an objective - maximize to test the + # piecewise part + m.obj = Objective(expr=m.pw(m.x), sense=maximize) + + return m + + def test_concave_pwl_maximization(self): + """Test transformation of a concave piecewise linear function (-x² + approximation).""" + m = self.pwl_neg_x_squared() + # Apply the hypograph transformation + transformation = TransformationFactory('contrib.piecewise.epigraph_hypograph') + transformation.apply_to(m) + + # Verify the transformation created the expected components + # The transformation should have created: + # 1. A substitute variable + # 2. Epigraphical constraints (one for each linear segment) + + # Get the substitute var and the transformed function block + sub_var = m.pw.get_transformation_var(m.obj.expr) + self.assertIsInstance(sub_var, Var) + pw_block = sub_var.parent_block() + + # Check that epigraphical constraints exist + epigraphical_constraints = pw_block.component('epigraphical_constraints') + self.assertIsInstance(epigraphical_constraints, Constraint) + + # For 5 breakpoints, we have 4 segments, so we should have 4 constraints + self.assertEqual(len(epigraphical_constraints), 4) + + # Verify the objective now uses the substitute variable + # The original pw(m.x) should be replaced with substitute_var + obj_expr = m.obj.expr.expr + self.assertIs(obj_expr, sub_var) + + exprs = [ + sub_var - (-1.0 * m.x + 0.0), + sub_var - (-3.0 * m.x + 2.0), + sub_var - (-5.0 * m.x + 6.0), + sub_var - (-7.0 * m.x + 12.0), + ] + + # Verify each constraint is of the form: linear_func_expr >= substitute_var + # (since this is a concave function, we use hypograph) + for idx, cons in epigraphical_constraints.items(): + constraint = pw_block.epigraphical_constraints[idx] + # For hypograph of concave function, we expect >= constraints + # The constraint is stored as: linear_func_expr >= substitute_var + # which Pyomo represents as: substitute_var - linear_func_expr <= 0 + # So the upper bound should be 0 + self.assertEqual(constraint.upper, 0) + self.assertIsNone(constraint.lower) + assertExpressionsEqual(self, constraint.body, exprs[idx]) + + def pwl_indexed_convex_concave(self): + # Create a model with an IndexedPiecewiseLinearFunction with two + # indices: one convex (x²) and one concave (-x²). This time the + # piecewise linear functions are used in Constraints rather than + # an Objective. + m = ConcreteModel() + m.idx = Set(initialize=[1, 2]) + m.x = Var(m.idx, bounds=(0, 4)) + + breakpoints = [0, 1, 2, 3, 4] + + def func_rule(m, i): + return x_squared if i == 1 else neg_x_squared + + def convex_rule(m, i): + return True if i == 1 else False + + m.pw = PiecewiseLinearFunction( + m.idx, points=breakpoints, function_rule=func_rule, convex=convex_rule + ) + + m.c1 = Constraint(expr=m.pw[1](m.x[1]) <= 10) + m.c2 = Constraint(expr=m.pw[2](m.x[2]) >= -10) + + return m + + def test_indexed_convex_concave_constraints(self): + """Test transformation of an IndexedPiecewiseLinearFunction with one + convex and one concave index, used in Constraints.""" + m = self.pwl_indexed_convex_concave() + + transformation = TransformationFactory('contrib.piecewise.epigraph_hypograph') + transformation.apply_to(m) + + # pw[1] is convex, so it gets an epigraph. pw[2] is concave, so it + # gets a hypograph. + sub_var1 = m.pw[1].get_transformation_var(m.c1.body) + sub_var2 = m.pw[2].get_transformation_var(m.c2.body) + self.assertIsInstance(sub_var1, Var) + self.assertIsInstance(sub_var2, Var) + + # Verify the constraints now use the substitute variables + self.assertIs(m.c1.body.expr, sub_var1) + self.assertIsNone(m.c1.lower) + self.assertEqual(m.c1.upper, 10) + self.assertIs(m.c2.body.expr, sub_var2) + self.assertEqual(m.c2.lower, -10) + self.assertIsNone(m.c2.upper) + + pw_block1 = sub_var1.parent_block() + pw_block2 = sub_var2.parent_block() + + epigraphical_constraints1 = pw_block1.component('epigraphical_constraints') + epigraphical_constraints2 = pw_block2.component('epigraphical_constraints') + self.assertIsInstance(epigraphical_constraints1, Constraint) + self.assertIsInstance(epigraphical_constraints2, Constraint) + self.assertEqual(len(epigraphical_constraints1), 4) + self.assertEqual(len(epigraphical_constraints2), 4) + + exprs1 = [ + m.x[1] - sub_var1, + 3.0 * m.x[1] - 2.0 - sub_var1, + 5.0 * m.x[1] - 6.0 - sub_var1, + 7.0 * m.x[1] - 12.0 - sub_var1, + ] + for idx, cons in epigraphical_constraints1.items(): + self.assertEqual(cons.upper, 0) + self.assertIsNone(cons.lower) + assertExpressionsEqual(self, cons.body, exprs1[idx]) + + exprs2 = [ + sub_var2 - (-1.0 * m.x[2] + 0.0), + sub_var2 - (-3.0 * m.x[2] + 2.0), + sub_var2 - (-5.0 * m.x[2] + 6.0), + sub_var2 - (-7.0 * m.x[2] + 12.0), + ] + for idx, cons in epigraphical_constraints2.items(): + self.assertEqual(cons.upper, 0) + self.assertIsNone(cons.lower) + assertExpressionsEqual(self, cons.body, exprs2[idx]) + + def test_indexed_convex_concave_targets(self): + """Test that the 'targets' argument can be used to transform only + one of the indices of an IndexedPiecewiseLinearFunction.""" + m = self.pwl_indexed_convex_concave() + + transformation = TransformationFactory('contrib.piecewise.epigraph_hypograph') + transformation.apply_to(m, targets=[m.pw[1]]) + + # pw[1] should be transformed... + self.assertFalse(m.pw[1].active) + sub_var1 = m.pw[1].get_transformation_var(m.c1.body) + self.assertIsInstance(sub_var1, Var) + self.assertIs(m.c1.body.expr, sub_var1) + + pw_block1 = sub_var1.parent_block() + epigraphical_constraints1 = pw_block1.component('epigraphical_constraints') + self.assertIsInstance(epigraphical_constraints1, Constraint) + self.assertEqual(len(epigraphical_constraints1), 4) + + exprs1 = [ + m.x[1] - sub_var1, + 3.0 * m.x[1] - 2.0 - sub_var1, + 5.0 * m.x[1] - 6.0 - sub_var1, + 7.0 * m.x[1] - 12.0 - sub_var1, + ] + for idx, cons in epigraphical_constraints1.items(): + self.assertEqual(cons.upper, 0) + self.assertIsNone(cons.lower) + assertExpressionsEqual(self, cons.body, exprs1[idx]) + + # ...but pw[2] should *not* be transformed. + self.assertTrue(m.pw[2].active) + self.assertIsNone(m.pw[2].get_transformation_var(m.c2.body)) + + +if __name__ == '__main__': + unittest.main() diff --git a/pyomo/contrib/piecewise/tests/test_piecewise_linear_function.py b/pyomo/contrib/piecewise/tests/test_piecewise_linear_function.py index 5bfc564fbad..d636699091b 100644 --- a/pyomo/contrib/piecewise/tests/test_piecewise_linear_function.py +++ b/pyomo/contrib/piecewise/tests/test_piecewise_linear_function.py @@ -416,6 +416,72 @@ def g2(x1, x2): self.assertAlmostEqual(m.pw(0.2, 4.3), g2(0.2, 4.3)) +class TestConvexOption(unittest.TestCase): + def test_scalar_pw_convex_option(self): + """Test that convex option is correctly set for ScalarPiecewiseLinearFunction""" + m = ConcreteModel() + m.x = Var(bounds=(1, 10)) + + # Test convex=True + m.pw_convex = PiecewiseLinearFunction( + points=[1, 3, 6, 10], function=log, convex=True + ) + self.assertTrue(m.pw_convex.convex) + + # Test convex=False + m.pw_concave = PiecewiseLinearFunction( + points=[1, 3, 6, 10], function=log, convex=False + ) + self.assertFalse(m.pw_concave.convex) + + # Test convex=None (default) + m.pw_none = PiecewiseLinearFunction(points=[1, 3, 6, 10], function=log) + self.assertIsNone(m.pw_none.convex) + + def test_indexed_pw_convex_option_points(self): + """Test that convex option is correctly set for + IndexedPiecewiseLinearFunction""" + m = ConcreteModel() + m.x = Var([1, 2], bounds=(1, 10)) + + def g1(x): + return x**2 + + def g2(x): + return log(x) + + m.funcs = {1: g1, 2: g2} + + # True and False + m.pw_convex = PiecewiseLinearFunction( + [1, 2], + points=[1, 3, 6, 10], + function_rule=m.funcs, + convex={1: True, 2: False}, + ) + self.assertTrue(m.pw_convex[1].convex) + self.assertFalse(m.pw_convex[2].convex) + + def convex_rule(m, i): + if i == 1: + return True + return None + + # Test convex=False for both indices + m.pw_concave = PiecewiseLinearFunction( + [1, 2], points=[1, 3, 6, 10], function_rule=m.funcs, convex=convex_rule + ) + self.assertTrue(m.pw_concave[1].convex) + self.assertIsNone(m.pw_concave[2].convex) + + # Test convex=None (default) for both indices + m.pw_none = PiecewiseLinearFunction( + [1, 2], points=[1, 3, 6, 10], function_rule=lambda m, i: m.funcs[i] + ) + self.assertIsNone(m.pw_none[1].convex) + self.assertIsNone(m.pw_none[2].convex) + + class TestTriangulationProducesDegenerateSimplices(unittest.TestCase): cube_extreme_pt_indices = [ {10, 11, 13, 14, 19, 20, 22, 23}, # right bottom back diff --git a/pyomo/contrib/piecewise/transform/epigraph_hypograph.py b/pyomo/contrib/piecewise/transform/epigraph_hypograph.py new file mode 100644 index 00000000000..5d1845d56c8 --- /dev/null +++ b/pyomo/contrib/piecewise/transform/epigraph_hypograph.py @@ -0,0 +1,70 @@ +# ____________________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2026 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 pyomo.common.config import ConfigDict, ConfigValue +from pyomo.core.base import TransformationFactory +from pyomo.contrib.piecewise.transform.piecewise_linear_transformation_base import ( + PiecewiseLinearTransformationBase, +) +from pyomo.contrib.piecewise.transform.piecewise_to_mip_visitor import ( + PiecewiseLinearToMIP, +) +from pyomo.environ import Constraint, NonNegativeIntegers, Var + + +@TransformationFactory.register( + 'contrib.piecewise.epigraph_hypograph', + doc="Transforms convex/concave piecewise linear functions to their epigraphical/" + "hypographical LP formulations", +) +class PWLToEpigraphOrHypograph(PiecewiseLinearTransformationBase): + CONFIG = PiecewiseLinearTransformationBase.CONFIG() + _transformation_name = 'pw_linear_epigraph' + + def _transform_pw_linear_expr(self, pw_expr, pw_linear_func, transformation_block): + transBlock = transformation_block.transformed_functions[ + len(transformation_block.transformed_functions) + ] + # get the PiecewiseLinearFunctionExpression + dimension = pw_expr.nargs() + + # Create variable to substitute + substitute_var = transBlock.substitute_var = Var() + pw_linear_func.map_transformation_var(pw_expr, substitute_var) + + transBlock.epigraphical_constraints = Constraint(NonNegativeIntegers) + + epigraph = False + if pw_linear_func.convex is None: + # we should autodetect if the dimension isn't insane, yell otherwise + raise NotImplementedError( + "It is (quadratically, in the number of pieces) possible to " + "auto-detect if a piecewise-linear " + "function is convex or concave, but we don't currently have an " + "implementation. PRs are welcome, or manually specify in the " + "PiecewiseLinearFunction constructor the convexity/concavity " + "using the 'convex' argument." + ) + if pw_linear_func.convex: + # This will be epigraph + epigraph = True + # Else this is the hypograph (the function is concave) + + for idx, linear_func in enumerate(pw_linear_func._linear_functions): + linear_func_expr = linear_func(*pw_expr.args) + if epigraph: + transBlock.epigraphical_constraints[idx] = ( + linear_func_expr <= substitute_var + ) + else: + transBlock.epigraphical_constraints[idx] = ( + linear_func_expr >= substitute_var + ) + + return substitute_var