From 19af1b443b41045d1cccfda6d807e79f219f85ec Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Thu, 9 Mar 2023 16:18:34 -0700 Subject: [PATCH 01/11] Putting everything except the '_transform_pw_linear_expr' method in a base class that all the PW to GDP transformations can use. --- .../transform/inner_representation_gdp.py | 174 +--------------- .../piecewise_to_gdp_transformation.py | 188 ++++++++++++++++++ 2 files changed, 194 insertions(+), 168 deletions(-) create mode 100644 pyomo/contrib/piecewise/transform/piecewise_to_gdp_transformation.py diff --git a/pyomo/contrib/piecewise/transform/inner_representation_gdp.py b/pyomo/contrib/piecewise/transform/inner_representation_gdp.py index b76ba070200..30d39bfe61d 100644 --- a/pyomo/contrib/piecewise/transform/inner_representation_gdp.py +++ b/pyomo/contrib/piecewise/transform/inner_representation_gdp.py @@ -9,30 +9,21 @@ # This software is distributed under the 3-clause BSD License. # ___________________________________________________________________________ -from pyomo.common.collections import ComponentMap from pyomo.common.config import ConfigDict, ConfigValue -from pyomo.common.modeling import unique_component_name from pyomo.contrib.fbbt.fbbt import compute_bounds_on_expr -from pyomo.contrib.piecewise import PiecewiseLinearFunction -from pyomo.contrib.piecewise.transform.piecewise_to_mip_visitor import ( - PiecewiseLinearToMIP) -from pyomo.core import ( - Constraint, Objective, Var, BooleanVar, Expression, Suffix, Param, Set, - SetOf, RangeSet, ExternalFunction, Connector, SortComponents, Any, - NonNegativeIntegers, NonNegativeReals) -from pyomo.core.base import Transformation, TransformationFactory -from pyomo.core.base.block import _BlockData, Block +from pyomo.contrib.piecewise.transform.piecewise_to_gdp_transformation import ( + PiecewiseLinearToGDP) +from pyomo.core import Constraint, NonNegativeIntegers, Suffix, Var +from pyomo.core.base import TransformationFactory from pyomo.core.util import target_list from pyomo.gdp import Disjunct, Disjunction -from pyomo.gdp.util import is_child_of -from pyomo.network import Port @TransformationFactory.register('contrib.inner_repn_gdp', doc="Convert piecewise-linear model to a GDP " "using an inner representation of the " "simplices that are the domains of the linear " "functions.") -class InnerRepresentationGDPTransformation(Transformation): +class InnerRepresentationGDPTransformation(PiecewiseLinearToGDP): """ Convert a model involving piecewise linear expressions into a GDP by representing the piecewise linear functions as Disjunctions where the @@ -85,115 +76,7 @@ class InnerRepresentationGDPTransformation(Transformation): transformed. """ )) - def __init__(self): - super().__init__() - self.handlers = { - Constraint: self._transform_constraint, - Objective: self._transform_objective, - Var: False, - BooleanVar: False, - Connector: False, - Expression: False, - Suffix: False, - Param: False, - Set: False, - SetOf: False, - RangeSet: False, - Disjunction: False, - Disjunct: False, - Block: self._transform_block, - ExternalFunction: False, - Port: False, - PiecewiseLinearFunction: self._transform_piecewise_linear_function, - } - self._transformation_blocks = {} - - def _apply_to(self, instance, **kwds): - try: - self._apply_to_impl(instance, **kwds) - finally: - self._transformation_blocks.clear() - - def _apply_to_impl(self, instance, **kwds): - config = self.CONFIG(kwds.pop('options', {})) - config.set_value(kwds) - - targets = config.targets - if targets is None: - targets = (instance, ) - - knownBlocks = {} - not_walking_exprs_msg = ( - "When not descending into expressions, Constraints " - "and Objectives are not valid targets. Please specify " - "PiecewiseLinearFunction component and the Blocks " - "containing them, or (at the cost of some performance " - "in this transformation), set the 'descend_into_expressions' " - "option to 'True'.") - for t in targets: - if not is_child_of(parent=instance, child=t, - knownBlocks=knownBlocks): - raise ValueError("Target '%s' is not a component on instance " - "'%s'!" % (t.name, instance.name)) - if t.ctype is PiecewiseLinearFunction: - if config.descend_into_expressions: - raise ValueError( - "When descending into expressions, the transformation " - "cannot take PiecewiseLinearFunction components as " - "targets. Please instead specify the Blocks, " - "Constraints, and Objectives where your " - "PiecewiseLinearFunctions have been used in " - "expressions.") - self._transform_piecewise_linear_function( - t, config.descend_into_expressions) - elif t.ctype is Block or isinstance(t, _BlockData): - self._transform_block(t, config.descend_into_expressions) - elif t.ctype is Constraint: - if not config.descend_into_expressions: - raise ValueError( - "Encountered Constraint target '%s':\n%s" - % (t.name, not_walking_exprs_msg)) - self._transform_constraint(t, config.descend_into_expressions) - elif t.ctype is Objective: - if not config.descend_into_expressions: - raise ValueError( - "Encountered Objective target '%s':\n%s" - % (t.name, not_walking_exprs_msg)) - self._transform_objective(t, config.descend_into_expressions) - else: - raise ValueError( - "Target '%s' is not a PiecewiseLinearFunction, Block or " - "Constraint. It was of type '%s' and can't be transformed." - % (t.name, type(t))) - - def _get_transformation_block(self, parent): - if parent in self._transformation_blocks: - return self._transformation_blocks[parent] - - nm = unique_component_name( - parent, - '_pyomo_contrib_pw_linear_inner_repn') - self._transformation_blocks[parent] = transBlock = Block() - parent.add_component(nm, transBlock) - - transBlock.transformed_functions = Block(Any) - return transBlock - - def _transform_block(self, block, descend_into_expressions): - blocks = block.values() if block.is_indexed() else (block,) - for b in blocks: - for obj in b.component_objects( - active=True, - descend_into=(Block, Disjunct), - sort=SortComponents.deterministic): - handler = self.handlers.get(obj.ctype, None) - if not handler: - if handler is None: - raise RuntimeError( - "No transformation handler registered for modeling " - "components of type '%s'." % obj.ctype) - continue - handler(obj, descend_into_expressions) + _transformation_name = 'pw_linear_inner_repn' def _transform_pw_linear_expr(self, pw_expr, pw_linear_func, transformation_block): @@ -245,48 +128,3 @@ def linear_combo(disj, i): expr=[d for d in transBlock.disjuncts.values()]) return transBlock.substitute_var - - def _transform_piecewise_linear_function(self, pw_linear_func, - descend_into_expressions): - if descend_into_expressions: - return - - transBlock = self._get_transformation_block( - pw_linear_func.parent_block()) - _functions = pw_linear_func.values() if pw_linear_func.is_indexed() \ - else (pw_linear_func,) - for pw_func in _functions: - for pw_expr in pw_func._expressions.values(): - substitute_var = self._transform_pw_linear_expr(pw_expr.expr, - pw_func, - transBlock) - # We change the named expression to point to the variable that - # will take the appropriate value of the piecewise linear - # function. - pw_expr.expr = substitute_var - - def _transform_constraint(self, constraint, descend_into_expressions): - if not descend_into_expressions: - return - - transBlock = self._get_transformation_block(constraint.parent_block()) - visitor = PiecewiseLinearToMIP(self._transform_pw_linear_expr, - transBlock) - - _constraints = constraint.values() if constraint.is_indexed() else \ - (constraint,) - for c in _constraints: - visitor.walk_expression((c.expr, c, 0)) - - def _transform_objective(self, objective, descend_into_expressions): - if not descend_into_expressions: - return - - transBlock = self._get_transformation_block(objective.parent_block()) - visitor = PiecewiseLinearToMIP(self._transform_pw_linear_expr, - transBlock) - - _objectives = objective.values() if objective.is_indexed() else \ - (objective,) - for o in _objectives: - visitor.walk_expression((o.expr, o, 0)) diff --git a/pyomo/contrib/piecewise/transform/piecewise_to_gdp_transformation.py b/pyomo/contrib/piecewise/transform/piecewise_to_gdp_transformation.py new file mode 100644 index 00000000000..4b09398f524 --- /dev/null +++ b/pyomo/contrib/piecewise/transform/piecewise_to_gdp_transformation.py @@ -0,0 +1,188 @@ +# ___________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2022 +# 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.errors import DeveloperError +from pyomo.common.modeling import unique_component_name +from pyomo.contrib.piecewise import PiecewiseLinearFunction +from pyomo.contrib.piecewise.transform.piecewise_to_mip_visitor import ( + PiecewiseLinearToMIP) +from pyomo.core import ( + Constraint, Objective, Var, BooleanVar, Expression, Suffix, Param, Set, + SetOf, RangeSet, ExternalFunction, Connector, SortComponents, Any) +from pyomo.core.base import Transformation +from pyomo.core.base.block import _BlockData, Block +from pyomo.gdp import Disjunct, Disjunction +from pyomo.gdp.util import is_child_of +from pyomo.network import Port + +class PiecewiseLinearToGDP(Transformation): + """ + Base class for transformations of piecewise-linear models to GDPs + """ + def __init__(self): + super().__init__() + self.handlers = { + Constraint: self._transform_constraint, + Objective: self._transform_objective, + Var: False, + BooleanVar: False, + Connector: False, + Expression: False, + Suffix: False, + Param: False, + Set: False, + SetOf: False, + RangeSet: False, + Disjunction: False, + Disjunct: False, + Block: self._transform_block, + ExternalFunction: False, + Port: False, + PiecewiseLinearFunction: self._transform_piecewise_linear_function, + } + self._transformation_blocks = {} + + def _apply_to(self, instance, **kwds): + try: + self._apply_to_impl(instance, **kwds) + finally: + self._transformation_blocks.clear() + + def _apply_to_impl(self, instance, **kwds): + config = self.CONFIG(kwds.pop('options', {})) + config.set_value(kwds) + + targets = config.targets + if targets is None: + targets = (instance, ) + + knownBlocks = {} + not_walking_exprs_msg = ( + "When not descending into expressions, Constraints " + "and Objectives are not valid targets. Please specify " + "PiecewiseLinearFunction component and the Blocks " + "containing them, or (at the cost of some performance " + "in this transformation), set the 'descend_into_expressions' " + "option to 'True'.") + for t in targets: + if not is_child_of(parent=instance, child=t, + knownBlocks=knownBlocks): + raise ValueError("Target '%s' is not a component on instance " + "'%s'!" % (t.name, instance.name)) + if t.ctype is PiecewiseLinearFunction: + if config.descend_into_expressions: + raise ValueError( + "When descending into expressions, the transformation " + "cannot take PiecewiseLinearFunction components as " + "targets. Please instead specify the Blocks, " + "Constraints, and Objectives where your " + "PiecewiseLinearFunctions have been used in " + "expressions.") + self._transform_piecewise_linear_function( + t, config.descend_into_expressions) + elif t.ctype is Block or isinstance(t, _BlockData): + self._transform_block(t, config.descend_into_expressions) + elif t.ctype is Constraint: + if not config.descend_into_expressions: + raise ValueError( + "Encountered Constraint target '%s':\n%s" + % (t.name, not_walking_exprs_msg)) + self._transform_constraint(t, config.descend_into_expressions) + elif t.ctype is Objective: + if not config.descend_into_expressions: + raise ValueError( + "Encountered Objective target '%s':\n%s" + % (t.name, not_walking_exprs_msg)) + self._transform_objective(t, config.descend_into_expressions) + else: + raise ValueError( + "Target '%s' is not a PiecewiseLinearFunction, Block or " + "Constraint. It was of type '%s' and can't be transformed." + % (t.name, type(t))) + + def _get_transformation_block(self, parent): + if parent in self._transformation_blocks: + return self._transformation_blocks[parent] + + nm = unique_component_name( + parent, + '_pyomo_contrib_%s' % self._transformation_name) + self._transformation_blocks[parent] = transBlock = Block() + parent.add_component(nm, transBlock) + + transBlock.transformed_functions = Block(Any) + return transBlock + + def _transform_block(self, block, descend_into_expressions): + blocks = block.values() if block.is_indexed() else (block,) + for b in blocks: + for obj in b.component_objects( + active=True, + descend_into=(Block, Disjunct), + sort=SortComponents.deterministic): + handler = self.handlers.get(obj.ctype, None) + if not handler: + if handler is None: + raise RuntimeError( + "No transformation handler registered for modeling " + "components of type '%s'." % obj.ctype) + continue + handler(obj, descend_into_expressions) + + def _transform_piecewise_linear_function(self, pw_linear_func, + descend_into_expressions): + if descend_into_expressions: + return + + transBlock = self._get_transformation_block( + pw_linear_func.parent_block()) + _functions = pw_linear_func.values() if pw_linear_func.is_indexed() \ + else (pw_linear_func,) + for pw_func in _functions: + for pw_expr in pw_func._expressions.values(): + substitute_var = self._transform_pw_linear_expr(pw_expr.expr, + pw_func, + transBlock) + # We change the named expression to point to the variable that + # will take the appropriate value of the piecewise linear + # function. + pw_expr.expr = substitute_var + + def _transform_constraint(self, constraint, descend_into_expressions): + if not descend_into_expressions: + return + + transBlock = self._get_transformation_block(constraint.parent_block()) + visitor = PiecewiseLinearToMIP(self._transform_pw_linear_expr, + transBlock) + + _constraints = constraint.values() if constraint.is_indexed() else \ + (constraint,) + for c in _constraints: + visitor.walk_expression((c.expr, c, 0)) + + def _transform_objective(self, objective, descend_into_expressions): + if not descend_into_expressions: + return + + transBlock = self._get_transformation_block(objective.parent_block()) + visitor = PiecewiseLinearToMIP(self._transform_pw_linear_expr, + transBlock) + + _objectives = objective.values() if objective.is_indexed() else \ + (objective,) + for o in _objectives: + visitor.walk_expression((o.expr, o, 0)) + + def _transform_pw_linear_expr(self, pw_expr, pw_linear_func, + transformation_block): + raise DeveloperError("Derived class failed to implement " + "'_transform_pw_linear_expr'") From 82bff935cb8f94e586f8bb1613107196951873fe Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Fri, 10 Mar 2023 21:03:45 -0700 Subject: [PATCH 02/11] Creating common place for reused tests and testing models --- pyomo/contrib/piecewise/tests/common_tests.py | 23 +++++ pyomo/contrib/piecewise/tests/models.py | 63 ++++++++++++++ .../piecewise/tests/test_inner_repn_gdp.py | 83 +++---------------- 3 files changed, 97 insertions(+), 72 deletions(-) create mode 100644 pyomo/contrib/piecewise/tests/common_tests.py create mode 100644 pyomo/contrib/piecewise/tests/models.py diff --git a/pyomo/contrib/piecewise/tests/common_tests.py b/pyomo/contrib/piecewise/tests/common_tests.py new file mode 100644 index 00000000000..7649f0e094a --- /dev/null +++ b/pyomo/contrib/piecewise/tests/common_tests.py @@ -0,0 +1,23 @@ +# ___________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2022 +# 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.core import Var +from pyomo.gdp import Disjunct, Disjunction + +def check_trans_block_structure(test_case, block): + # One (indexed) disjunct + test_case.assertEqual(len(block.component_map(Disjunct)), 1) + # One disjunction + test_case.assertEqual(len(block.component_map(Disjunction)), 1) + # The 'z' var (that we will substitute in for the function being + # approximated) is here: + test_case.assertEqual(len(block.component_map(Var)), 1) + test_case.assertIsInstance(block.substitute_var, Var) diff --git a/pyomo/contrib/piecewise/tests/models.py b/pyomo/contrib/piecewise/tests/models.py new file mode 100644 index 00000000000..9639a95fbb8 --- /dev/null +++ b/pyomo/contrib/piecewise/tests/models.py @@ -0,0 +1,63 @@ +# ___________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2022 +# 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.contrib.piecewise import PiecewiseLinearFunction +from pyomo.environ import ConcreteModel, Constraint, log, Objective, Var + +def make_log_x_model(): + m = ConcreteModel() + m.x = Var(bounds=(1, 10)) + def log_function(x): + return log(x) + m.log_function = log_function + m.pw_log = PiecewiseLinearFunction(points=[1, 3, 6, 10], + function=m.log_function) + + # Here are the linear functions, for safe keeping. + def f1(x): + return (log(3)/2)*x - log(3)/2 + m.f1 = f1 + def f2(x): + return (log(2)/3)*x + log(3/2) + m.f2 = f2 + def f3(x): + return (log(5/3)/4)*x + log(6/((5/3)**(3/2))) + m.f3 = f3 + + m.log_expr = m.pw_log(m.x) + m.obj = Objective(expr=m.log_expr) + + m.x1 = Var(bounds=(0, 3)) + m.x2 = Var(bounds=(1, 7)) + + ## apprximates paraboloid x1**2 + x2**2 + def g1(x1, x2): + return 3*x1 + 5*x2 - 4 + m.g1 = g1 + def g2(x1, x2): + return 3*x1 + 11*x2 - 28 + m.g2 = g2 + simplices = [[(0, 1), (0, 4), (3, 4)], + [(0, 1), (3, 4), (3, 1)], + [(3, 4), (3, 7), (0, 7)], + [(0, 7), (0, 4), (3, 4)]] + m.pw_paraboloid = PiecewiseLinearFunction(simplices=simplices, + linear_functions=[g1, g1, g2, + g2]) + m.paraboloid_expr = m.pw_paraboloid(m.x1, m.x2) + def c_rule(m, i): + if i == 0: + return m.x >= m.paraboloid_expr + else: + return (1, m.x1, 2) + m.indexed_c = Constraint([0, 1], rule=c_rule) + + return m diff --git a/pyomo/contrib/piecewise/tests/test_inner_repn_gdp.py b/pyomo/contrib/piecewise/tests/test_inner_repn_gdp.py index 78ced23b5eb..a3bea4a44b0 100644 --- a/pyomo/contrib/piecewise/tests/test_inner_repn_gdp.py +++ b/pyomo/contrib/piecewise/tests/test_inner_repn_gdp.py @@ -10,76 +10,15 @@ # ___________________________________________________________________________ import pyomo.common.unittest as unittest -from pyomo.contrib.piecewise import ( - PiecewiseLinearFunction, PiecewiseLinearExpression) +from pyomo.contrib.piecewise.tests import models +import pyomo.contrib.piecewise.tests.common_tests as ct from pyomo.core.base import TransformationFactory from pyomo.core.expr.compare import ( assertExpressionsEqual, assertExpressionsStructurallyEqual) from pyomo.gdp import Disjunct, Disjunction -from pyomo.environ import ( - ConcreteModel, Constraint, log, SolverFactory, Objective, value, Var) +from pyomo.environ import Constraint, SolverFactory, value, Var class TestTransformPiecewiseModelToInnerRepnGDP(unittest.TestCase): - def make_model(self): - m = ConcreteModel() - m.x = Var(bounds=(1, 10)) - def log_function(x): - return log(x) - m.log_function = log_function - m.pw_log = PiecewiseLinearFunction(points=[1, 3, 6, 10], - function=m.log_function) - - # Here are the linear functions, for safe keeping. - def f1(x): - return (log(3)/2)*x - log(3)/2 - m.f1 = f1 - def f2(x): - return (log(2)/3)*x + log(3/2) - m.f2 = f2 - def f3(x): - return (log(5/3)/4)*x + log(6/((5/3)**(3/2))) - m.f3 = f3 - - m.log_expr = m.pw_log(m.x) - m.obj = Objective(expr=m.log_expr) - - m.x1 = Var(bounds=(0, 3)) - m.x2 = Var(bounds=(1, 7)) - - ## apprximates paraboloid x1**2 + x2**2 - def g1(x1, x2): - return 3*x1 + 5*x2 - 4 - m.g1 = g1 - def g2(x1, x2): - return 3*x1 + 11*x2 - 28 - m.g2 = g2 - simplices = [[(0, 1), (0, 4), (3, 4)], - [(0, 1), (3, 4), (3, 1)], - [(3, 4), (3, 7), (0, 7)], - [(0, 7), (0, 4), (3, 4)]] - m.pw_paraboloid = PiecewiseLinearFunction(simplices=simplices, - linear_functions=[g1, g1, g2, - g2]) - m.paraboloid_expr = m.pw_paraboloid(m.x1, m.x2) - def c_rule(m, i): - if i == 0: - return m.x >= m.paraboloid_expr - else: - return (1, m.x1, 2) - m.indexed_c = Constraint([0, 1], rule=c_rule) - - return m - - def check_trans_block_structure(self, block): - # One (indexed) disjunct - self.assertEqual(len(block.component_map(Disjunct)), 1) - # One disjunction - self.assertEqual(len(block.component_map(Disjunction)), 1) - # The 'z' var (that we will substitute in for the function being - # approximated) is here: - self.assertEqual(len(block.component_map(Var)), 1) - self.assertIsInstance(block.substitute_var, Var) - def check_log_disjunct(self, d, pts, f, substitute_var, x): self.assertEqual(len(d.component_map(Constraint)), 3) # lambdas and indicator_var @@ -135,7 +74,7 @@ def check_pw_log(self, m): self.assertIsInstance(z, Var) # Now we can use those Vars to check on what the transformation created log_block = z.parent_block() - self.check_trans_block_structure(log_block) + ct.check_trans_block_structure(self, log_block) # Check that all of the Disjuncts have what they should self.assertEqual(len(log_block.disjuncts), 3) @@ -164,7 +103,7 @@ def check_pw_paraboloid(self, m): z = m.pw_paraboloid.get_transformation_var(m.paraboloid_expr) self.assertIsInstance(z, Var) paraboloid_block = z.parent_block() - self.check_trans_block_structure(paraboloid_block) + ct.check_trans_block_structure(self, paraboloid_block) self.assertEqual(len(paraboloid_block.disjuncts), 4) disjuncts_dict = { @@ -190,7 +129,7 @@ def check_pw_paraboloid(self, m): paraboloid_block.substitute_var) def test_transformation_do_not_descend(self): - m = self.make_model() + m = models.make_log_x_model() inner_repn = TransformationFactory('contrib.inner_repn_gdp') inner_repn.apply_to(m) @@ -198,7 +137,7 @@ def test_transformation_do_not_descend(self): self.check_pw_paraboloid(m) def test_transformation_PiecewiseLinearFunction_targets(self): - m = self.make_model() + m = models.make_log_x_model() inner_repn = TransformationFactory('contrib.inner_repn_gdp') inner_repn.apply_to(m, targets=[m.pw_log]) @@ -209,7 +148,7 @@ def test_transformation_PiecewiseLinearFunction_targets(self): m.pw_paraboloid.get_transformation_var(m.paraboloid_expr)) def test_descend_into_expressions(self): - m = self.make_model() + m = models.make_log_x_model() inner_repn = TransformationFactory('contrib.inner_repn_gdp') inner_repn.apply_to(m, descend_into_expressions=True) @@ -218,7 +157,7 @@ def test_descend_into_expressions(self): self.check_pw_paraboloid(m) def test_descend_into_expressions_constraint_target(self): - m = self.make_model() + m = models.make_log_x_model() inner_repn = TransformationFactory('contrib.inner_repn_gdp') inner_repn.apply_to(m, descend_into_expressions=True, targets=[m.indexed_c]) @@ -228,7 +167,7 @@ def test_descend_into_expressions_constraint_target(self): self.assertIsNone(m.pw_log.get_transformation_var(m.log_expr)) def test_descend_into_expressions_objective_target(self): - m = self.make_model() + m = models.make_log_x_model() inner_repn = TransformationFactory('contrib.inner_repn_gdp') inner_repn.apply_to(m, descend_into_expressions=True, targets=[m.obj]) @@ -241,7 +180,7 @@ def test_descend_into_expressions_objective_target(self): @unittest.skipUnless(SolverFactory('gurobi').available(), 'Gurobi is not available') def test_solve_disaggregated_convex_combo_model(self): - m = self.make_model() + m = models.make_log_x_model() TransformationFactory( 'contrib.disaggregated_convex_combination').apply_to(m) From b221ba6beca7698126937547e600eae22691ba17 Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Fri, 10 Mar 2023 21:04:15 -0700 Subject: [PATCH 03/11] Adding outer rep'n transformation --- .../piecewise/tests/test_outer_repn_gdp.py | 189 ++++++++++++++++++ .../transform/outer_representation_gdp.py | 146 ++++++++++++++ 2 files changed, 335 insertions(+) create mode 100644 pyomo/contrib/piecewise/tests/test_outer_repn_gdp.py create mode 100644 pyomo/contrib/piecewise/transform/outer_representation_gdp.py diff --git a/pyomo/contrib/piecewise/tests/test_outer_repn_gdp.py b/pyomo/contrib/piecewise/tests/test_outer_repn_gdp.py new file mode 100644 index 00000000000..8118e33e8db --- /dev/null +++ b/pyomo/contrib/piecewise/tests/test_outer_repn_gdp.py @@ -0,0 +1,189 @@ +# ___________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2022 +# 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 math import sqrt +import pyomo.common.unittest as unittest +from pyomo.contrib.piecewise.tests import models +import pyomo.contrib.piecewise.tests.common_tests as ct +from pyomo.core.base import TransformationFactory +from pyomo.core.expr.compare import ( + assertExpressionsEqual, assertExpressionsStructurallyEqual) +from pyomo.gdp import Disjunct, Disjunction +from pyomo.environ import Constraint, SolverFactory, value, Var + +class TestTransformPiecewiseModelToOuterRepnGDP(unittest.TestCase): + def check_log_disjunct(self, d, pts, f, substitute_var, x): + # We can fit both bounds constraints in one constraint, then we have the + # linear function + self.assertEqual(len(d.component_map(Constraint)), 2) + # indicator_var + self.assertEqual(len(d.component_map(Var)), 1) + self.assertIsInstance(d.simplex_halfspaces, Constraint) + self.assertEqual(d.simplex_halfspaces.lower, pts[0]) + self.assertEqual(d.simplex_halfspaces.upper, pts[1]) + self.assertIs(d.simplex_halfspaces.body, x) + + self.assertIsInstance(d.set_substitute, Constraint) + assertExpressionsEqual(self, d.set_substitute.expr, + substitute_var == f(x), places=7) + + def check_paraboloid_disjunct(self, d, constraint_coefs, f, + substitute_var, x1, x2): + self.assertEqual(len(d.component_map(Constraint)), 2) + # just indicator_var + self.assertEqual(len(d.component_map(Var)), 1) + for i, cons in d.simplex_halfspaces.items(): + coefs = constraint_coefs[i] + assertExpressionsEqual(self, cons.expr, coefs[0]*x1 + coefs[1]*x2 + + coefs[2] <= 0, places=6) + + self.assertIsInstance(d.set_substitute, Constraint) + assertExpressionsEqual(self, d.set_substitute.expr, + substitute_var == f(x1, x2), places=7) + + def check_pw_log(self, m): + ## + # Check the transformation of the approximation of log(x) + ## + z = m.pw_log.get_transformation_var(m.log_expr) + self.assertIsInstance(z, Var) + # Now we can use those Vars to check on what the transformation created + log_block = z.parent_block() + ct.check_trans_block_structure(self, log_block) + + # Check that all of the Disjuncts have what they should + self.assertEqual(len(log_block.disjuncts), 3) + disjuncts_dict = { + log_block.disjuncts[0]: ((1, 3), m.f1), + log_block.disjuncts[1]: ((3, 6), m.f2), + log_block.disjuncts[2]: ((6, 10), m.f3), + } + for d, (pts, f) in disjuncts_dict.items(): + self.check_log_disjunct(d, pts, f, log_block.substitute_var, m.x) + + # Check the Disjunction + self.assertIsInstance(log_block.pick_a_piece, Disjunction) + self.assertEqual(len(log_block.pick_a_piece.disjuncts), 3) + for i in range(2): + self.assertIs(log_block.pick_a_piece.disjuncts[i], + log_block.disjuncts[i]) + + # And check the substitute Var is in the objective now. + self.assertIs(m.obj.expr.expr, log_block.substitute_var) + + def check_pw_paraboloid(self, m): + ## + # Check the approximation of the transformation of the paraboloid + ## + z = m.pw_paraboloid.get_transformation_var(m.paraboloid_expr) + self.assertIsInstance(z, Var) + paraboloid_block = z.parent_block() + ct.check_trans_block_structure(self, paraboloid_block) + + self.assertEqual(len(paraboloid_block.disjuncts), 4) + disjuncts_dict = { + # the normal vectors of the faces are normalized when we get + # them from scipy: + paraboloid_block.disjuncts[0]: ([[sqrt(2)/2, -sqrt(2)/2, sqrt(2)/2], + [-1.0, 0.0, 0.0], + [0.0, 1.0, -4.0]], m.g1), + paraboloid_block.disjuncts[1]: ([[-sqrt(2)/2, sqrt(2)/2, + -sqrt(2)/2], + [0.0, -1.0, 1.0], + [1.0, 0.0, -3.0], + ], m.g1), + paraboloid_block.disjuncts[2]: ([[-sqrt(2)/2, -sqrt(2)/2, + 7*sqrt(2)/2], + [0.0, 1.0, -7.0], + [1.0, 0.0, -3.0], + ], m.g2), + paraboloid_block.disjuncts[3]: ([[sqrt(2)/2, sqrt(2)/2, + -7*sqrt(2)/2], + [-1.0, 0.0, 0.0], + [0.0, -1.0, 4.0]], m.g2), + } + for d, (constraint_coefs, f) in disjuncts_dict.items(): + self.check_paraboloid_disjunct(d, constraint_coefs, f, + paraboloid_block.substitute_var, + m.x1, m.x2) + + # Check the Disjunction + self.assertIsInstance(paraboloid_block.pick_a_piece, Disjunction) + self.assertEqual(len(paraboloid_block.pick_a_piece.disjuncts), 4) + for i in range(3): + self.assertIs(paraboloid_block.pick_a_piece.disjuncts[i], + paraboloid_block.disjuncts[i]) + + # And check the substitute Var is in the objective now. + self.assertIs(m.indexed_c[0].body.args[0].expr, + paraboloid_block.substitute_var) + + def test_transformation_do_not_descend(self): + m = models.make_log_x_model() + outer_repn = TransformationFactory('contrib.outer_repn_gdp') + outer_repn.apply_to(m) + + self.check_pw_log(m) + self.check_pw_paraboloid(m) + + def test_transformation_PiecewiseLinearFunction_targets(self): + m = models.make_log_x_model() + outer_repn = TransformationFactory('contrib.outer_repn_gdp') + outer_repn.apply_to(m, targets=[m.pw_log]) + + self.check_pw_log(m) + + # And check that the paraboloid was *not* transformed. + self.assertIsNone( + m.pw_paraboloid.get_transformation_var(m.paraboloid_expr)) + + def test_descend_into_expressions(self): + m = models.make_log_x_model() + outer_repn = TransformationFactory('contrib.outer_repn_gdp') + outer_repn.apply_to(m, descend_into_expressions=True) + + # Everything should be transformed + self.check_pw_log(m) + self.check_pw_paraboloid(m) + + def test_descend_into_expressions_constraint_target(self): + m = models.make_log_x_model() + outer_repn = TransformationFactory('contrib.outer_repn_gdp') + outer_repn.apply_to(m, descend_into_expressions=True, + targets=[m.indexed_c]) + + self.check_pw_paraboloid(m) + # And check that the log was *not* transformed. + self.assertIsNone(m.pw_log.get_transformation_var(m.log_expr)) + + def test_descend_into_expressions_objective_target(self): + m = models.make_log_x_model() + outer_repn = TransformationFactory('contrib.outer_repn_gdp') + outer_repn.apply_to(m, descend_into_expressions=True, + targets=[m.obj]) + + self.check_pw_log(m) + # And check that the paraboloid was *not* transformed. + self.assertIsNone( + m.pw_paraboloid.get_transformation_var(m.paraboloid_expr)) + + @unittest.skipUnless(SolverFactory('gurobi').available(), + 'Gurobi is not available') + def test_solve_multiple_choice_model(self): + m = models.make_log_x_model() + TransformationFactory('contrib.piecewise.multiple_choice').apply_to(m) + + SolverFactory('gurobi').solve(m, tee=True) + + self.assertAlmostEqual(value(m.x), 4) + self.assertAlmostEqual(value(m.x1), 1) + self.assertAlmostEqual(value(m.x2), 1) + self.assertAlmostEqual(value(m.obj), m.f2(4)) diff --git a/pyomo/contrib/piecewise/transform/outer_representation_gdp.py b/pyomo/contrib/piecewise/transform/outer_representation_gdp.py new file mode 100644 index 00000000000..f32304f23d4 --- /dev/null +++ b/pyomo/contrib/piecewise/transform/outer_representation_gdp.py @@ -0,0 +1,146 @@ +# ___________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2022 +# 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 +import pyomo.common.dependencies.numpy as np +from pyomo.common.dependencies.scipy import spatial +from pyomo.contrib.fbbt.fbbt import compute_bounds_on_expr +from pyomo.contrib.piecewise.transform.piecewise_to_gdp_transformation import ( + PiecewiseLinearToGDP) +from pyomo.core import Constraint, NonNegativeIntegers, Suffix, Var +from pyomo.core.base import TransformationFactory +from pyomo.core.util import target_list +from pyomo.gdp import Disjunct, Disjunction + +## DEBUG +from pytest import set_trace + +@TransformationFactory.register('contrib.outer_repn_gdp', + doc="Convert piecewise-linear model to a GDP " + "using an outer (Ax <= b) representation of " + "the simplices that are the domains of the " + "linear functions.") +class OuterRepresentationGDPTransformation(PiecewiseLinearToGDP): + """ + Convert a model involving piecewise linear expressions into a GDP by + representing the piecewise linear functions as Disjunctions where the + simplices over which the linear functions are defined are represented + in an "outer" representation--in sets of constraints of the form Ax <= b. + + This transformation can be called in one of two ways: + 1) The default, where 'descend_into_expressions' is False. This is + more computationally efficient, but relies on the + PiecewiseLinearFunctions being declared on the same Block in which + they are used in Expressions (if you are hoping to maintain the + original hierarchical structure of the model). In this mode, + targets must be Blocks and/or PiecewiseLinearFunctions. + 2) With 'descend_into_expressions' True. This is less computationally + efficient, but will respect hierarchical structure by finding + uses of PiecewiseLinearFunctions in Constraint and Obective + expressions and putting their transformed counterparts on the same + parent Block as the component owning their parent expression. In + this mode, targets must be Blocks, Constraints, and/or Objectives. + """ + CONFIG = ConfigDict('piecewise.outer_repn_gdp') + CONFIG.declare('targets', ConfigValue( + default=None, + domain=target_list, + description="target or list of targets that will be transformed", + doc=""" + This specifies the list of components to transform. If None (default), + the entire model is transformed. Note that if the transformation is + done out of place, the list of targets should be attached to the model + before it is cloned, and the list will specify the targets on the cloned + instance.""" + )) + CONFIG.declare('descend_into_expressions', ConfigValue( + default=False, + domain=bool, + description="Whether to look for uses of PiecewiseLinearFunctions in " + "the Constraint and Objective expressions, rather than assuming " + "all PiecewiseLinearFunctions are on the active tree(s) of 'instance' " + "and 'targets.'", + doc=""" + It is *strongly* recommended that, in hierarchical models, the + PiecewiseLinearFunction components are on the same Block as where + they are used in expressions. If you follow this recommendation, + this option can remain False, which will make this transformation + more efficient. However, if you do not follow the recommendation, + unless you know what you are doing, turn this option to 'True' to + ensure that all of the uses of PiecewiseLinearFunctions are + transformed. + """ + )) + _transformation_name = 'pw_linear_outer_repn' + + 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() + transBlock.disjuncts = Disjunct(NonNegativeIntegers) + substitute_var = transBlock.substitute_var = Var() + pw_linear_func.map_transformation_var(pw_expr, + substitute_var) + substitute_var_lb = float('inf') + substitute_var_ub = -float('inf') + if dimension > 1: + A = np.ones((dimension + 1, dimension + 1)) + b = np.zeros(dimension + 1) + b[-1] = 1 + + for simplex, linear_func in zip(pw_linear_func._simplices, + pw_linear_func._linear_functions): + disj = transBlock.disjuncts[len(transBlock.disjuncts)] + + if dimension == 1: + # We don't need scipy, and the polytopes are 1-dimensional + # simplices, so they are defined by two bounds constraints: + disj.simplex_halfspaces = Constraint( + expr=(pw_linear_func._points[simplex[0]][0], + pw_expr.args[0], + pw_linear_func._points[simplex[1]][0])) + else: + disj.simplex_halfspaces = Constraint(range(dimension + 1)) + # we will use scipy to get the convex hull of the extreme pts of + # the simplex + extreme_pts = [] + for idx in simplex: + extreme_pts.append(pw_linear_func._points[idx]) + chull = spatial.ConvexHull(extreme_pts) + vars = pw_expr.args + for i, eqn in enumerate(chull.equations): + # The equations are given as normal vectors (A) followed by + # offsets (b) such that Ax + b <= 0 gives the halfspaces + # defining the simplex. (See Qhull documentation) + disj.simplex_halfspaces[i] = sum(eqn[j]*v for j, v in + enumerate(vars)) + \ + float(eqn[dimension]) <= 0 + + linear_func_expr = linear_func(*pw_expr.args) + disj.set_substitute = Constraint(expr=substitute_var == + linear_func_expr) + (lb, ub) = compute_bounds_on_expr(linear_func_expr) + if lb is not None and lb < substitute_var_lb: + substitute_var_lb = lb + if ub is not None and ub > substitute_var_ub: + substitute_var_ub = ub + + if substitute_var_lb < float('inf'): + transBlock.substitute_var.setlb(substitute_var_lb) + if substitute_var_ub > -float('inf'): + transBlock.substitute_var.setub(substitute_var_ub) + transBlock.pick_a_piece = Disjunction( + expr=[d for d in transBlock.disjuncts.values()]) + + return transBlock.substitute_var From 61828d2ec5506dbc179e0fc25447f6fd7d1be15f Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Fri, 10 Mar 2023 21:04:36 -0700 Subject: [PATCH 04/11] Adding multiple choice transformation --- pyomo/contrib/piecewise/__init__.py | 5 +++ .../piecewise/transform/multiple_choice.py | 35 +++++++++++++++++++ 2 files changed, 40 insertions(+) create mode 100644 pyomo/contrib/piecewise/transform/multiple_choice.py diff --git a/pyomo/contrib/piecewise/__init__.py b/pyomo/contrib/piecewise/__init__.py index c0ed7306abb..3fc578e0870 100644 --- a/pyomo/contrib/piecewise/__init__.py +++ b/pyomo/contrib/piecewise/__init__.py @@ -7,3 +7,8 @@ InnerRepresentationGDPTransformation) from pyomo.contrib.piecewise.transform.disaggregated_convex_combination import ( DisaggregatedConvexCombinationTransformation) +from pyomo.contrib.piecewise.transform.outer_representation_gdp import ( + OuterRepresentationGDPTransformation) +from pyomo.contrib.piecewise.transform.multiple_choice import ( + MultipleChoiceTransformation) + diff --git a/pyomo/contrib/piecewise/transform/multiple_choice.py b/pyomo/contrib/piecewise/transform/multiple_choice.py new file mode 100644 index 00000000000..9f546311a87 --- /dev/null +++ b/pyomo/contrib/piecewise/transform/multiple_choice.py @@ -0,0 +1,35 @@ +# ___________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2022 +# 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.core.base import Transformation, TransformationFactory +import pyomo.gdp.plugins.hull + +@TransformationFactory.register('contrib.piecewise.multiple_choice', + doc="Convert piecewise-linear model to a GDP " + "to 'Multiple Choice' MIP " + "formulation.") +class MultipleChoiceTransformation(Transformation): + """ + Converts a model containing PiecewiseLinearFunctions to a an equivalent + MIP via the Multiple Choice method from [1]. Note that, + while this model probably resolves to the model described in [1] after + presolve, the Pyomo version is not as simplified. + + References + ---------- + [1] J.P. Vielma, S. Ahmed, and G. Nemhauser, "Mixed-integer models + for nonseparable piecewise-linear optimization: unifying framework + and extensions," Operations Research, vol. 58, no. 2, pp. 305-315, + 2010. + """ + def _apply_to(self, instance, **kwds): + TransformationFactory('contrib.outer_repn_gdp').apply_to(instance) + TransformationFactory('gdp.hull').apply_to(instance) From 4cad405ce2eafd474001308cd4346072c05900e5 Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Fri, 10 Mar 2023 21:13:22 -0700 Subject: [PATCH 05/11] Skipping tests involving the paraboloid if scipy is not available --- pyomo/contrib/piecewise/tests/test_outer_repn_gdp.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/pyomo/contrib/piecewise/tests/test_outer_repn_gdp.py b/pyomo/contrib/piecewise/tests/test_outer_repn_gdp.py index 8118e33e8db..3eb44029b32 100644 --- a/pyomo/contrib/piecewise/tests/test_outer_repn_gdp.py +++ b/pyomo/contrib/piecewise/tests/test_outer_repn_gdp.py @@ -10,6 +10,7 @@ # ___________________________________________________________________________ from math import sqrt +from pyomo.common.dependencies import scipy_available import pyomo.common.unittest as unittest from pyomo.contrib.piecewise.tests import models import pyomo.contrib.piecewise.tests.common_tests as ct @@ -126,6 +127,7 @@ def check_pw_paraboloid(self, m): self.assertIs(m.indexed_c[0].body.args[0].expr, paraboloid_block.substitute_var) + @unittest.skipUnless(scipy_available, "Scipy is not available") def test_transformation_do_not_descend(self): m = models.make_log_x_model() outer_repn = TransformationFactory('contrib.outer_repn_gdp') @@ -145,6 +147,7 @@ def test_transformation_PiecewiseLinearFunction_targets(self): self.assertIsNone( m.pw_paraboloid.get_transformation_var(m.paraboloid_expr)) + @unittest.skipUnless(scipy_available, "Scipy is not available") def test_descend_into_expressions(self): m = models.make_log_x_model() outer_repn = TransformationFactory('contrib.outer_repn_gdp') @@ -154,6 +157,7 @@ def test_descend_into_expressions(self): self.check_pw_log(m) self.check_pw_paraboloid(m) + @unittest.skipUnless(scipy_available, "Scipy is not available") def test_descend_into_expressions_constraint_target(self): m = models.make_log_x_model() outer_repn = TransformationFactory('contrib.outer_repn_gdp') @@ -175,8 +179,9 @@ def test_descend_into_expressions_objective_target(self): self.assertIsNone( m.pw_paraboloid.get_transformation_var(m.paraboloid_expr)) - @unittest.skipUnless(SolverFactory('gurobi').available(), - 'Gurobi is not available') + @unittest.skipUnless(SolverFactory('gurobi').available() and + scipy_available, + 'Gurobi and/or scipy is not available') def test_solve_multiple_choice_model(self): m = models.make_log_x_model() TransformationFactory('contrib.piecewise.multiple_choice').apply_to(m) From a0ddb42de2a86d9abb1f06fcb46e1aeb654ed513 Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Fri, 10 Mar 2023 21:34:11 -0700 Subject: [PATCH 06/11] Moving the common tests to the common tests file --- pyomo/contrib/piecewise/tests/common_tests.py | 66 +++++++++++++++++-- .../piecewise/tests/test_inner_repn_gdp.py | 62 +++++------------ .../piecewise/tests/test_outer_repn_gdp.py | 57 ++++------------ 3 files changed, 88 insertions(+), 97 deletions(-) diff --git a/pyomo/contrib/piecewise/tests/common_tests.py b/pyomo/contrib/piecewise/tests/common_tests.py index 7649f0e094a..e6c681d47b3 100644 --- a/pyomo/contrib/piecewise/tests/common_tests.py +++ b/pyomo/contrib/piecewise/tests/common_tests.py @@ -9,15 +9,71 @@ # This software is distributed under the 3-clause BSD License. # ___________________________________________________________________________ +import pyomo.contrib.piecewise.tests.models as models from pyomo.core import Var +from pyomo.core.base import TransformationFactory +from pyomo.environ import value from pyomo.gdp import Disjunct, Disjunction -def check_trans_block_structure(test_case, block): +def check_trans_block_structure(test, block): # One (indexed) disjunct - test_case.assertEqual(len(block.component_map(Disjunct)), 1) + test.assertEqual(len(block.component_map(Disjunct)), 1) # One disjunction - test_case.assertEqual(len(block.component_map(Disjunction)), 1) + test.assertEqual(len(block.component_map(Disjunction)), 1) # The 'z' var (that we will substitute in for the function being # approximated) is here: - test_case.assertEqual(len(block.component_map(Var)), 1) - test_case.assertIsInstance(block.substitute_var, Var) + test.assertEqual(len(block.component_map(Var)), 1) + test.assertIsInstance(block.substitute_var, Var) + +def check_log_x_model_soln(test, m): + test.assertAlmostEqual(value(m.x), 4) + test.assertAlmostEqual(value(m.x1), 1) + test.assertAlmostEqual(value(m.x2), 1) + test.assertAlmostEqual(value(m.obj), m.f2(4)) + +def check_transformation_do_not_descend(test, transformation): + m = models.make_log_x_model() + transform = TransformationFactory(transformation) + transform.apply_to(m) + + test.check_pw_log(m) + test.check_pw_paraboloid(m) + +def check_transformation_PiecewiseLinearFunction_targets(test, transformation): + m = models.make_log_x_model() + transform = TransformationFactory(transformation) + transform.apply_to(m, targets=[m.pw_log]) + + test.check_pw_log(m) + + # And check that the paraboloid was *not* transformed. + test.assertIsNone( + m.pw_paraboloid.get_transformation_var(m.paraboloid_expr)) + +def check_descend_into_expressions(test, transformation): + m = models.make_log_x_model() + transform = TransformationFactory(transformation) + transform.apply_to(m, descend_into_expressions=True) + + # Everything should be transformed + test.check_pw_log(m) + test.check_pw_paraboloid(m) + +def check_descend_into_expressions_constraint_target(test, transformation): + m = models.make_log_x_model() + transform = TransformationFactory(transformation) + transform.apply_to(m, descend_into_expressions=True, targets=[m.indexed_c]) + + test.check_pw_paraboloid(m) + # And check that the log was *not* transformed. + test.assertIsNone(m.pw_log.get_transformation_var(m.log_expr)) + +def check_descend_into_expressions_objective_target(test, transformation): + m = models.make_log_x_model() + transform = TransformationFactory(transformation) + transform.apply_to(m, descend_into_expressions=True, targets=[m.obj]) + + test.check_pw_log(m) + # And check that the paraboloid was *not* transformed. + test.assertIsNone( + m.pw_paraboloid.get_transformation_var(m.paraboloid_expr)) diff --git a/pyomo/contrib/piecewise/tests/test_inner_repn_gdp.py b/pyomo/contrib/piecewise/tests/test_inner_repn_gdp.py index 70e1d30f736..5b701ae4ac1 100644 --- a/pyomo/contrib/piecewise/tests/test_inner_repn_gdp.py +++ b/pyomo/contrib/piecewise/tests/test_inner_repn_gdp.py @@ -16,7 +16,7 @@ from pyomo.core.expr.compare import ( assertExpressionsEqual, assertExpressionsStructurallyEqual) from pyomo.gdp import Disjunct, Disjunction -from pyomo.environ import Constraint, SolverFactory, value, Var +from pyomo.environ import Constraint, SolverFactory, Var class TestTransformPiecewiseModelToInnerRepnGDP(unittest.TestCase): def check_log_disjunct(self, d, pts, f, substitute_var, x): @@ -129,53 +129,27 @@ def check_pw_paraboloid(self, m): paraboloid_block.substitute_var) def test_transformation_do_not_descend(self): - m = models.make_log_x_model() - inner_repn = TransformationFactory('contrib.piecewise.inner_repn_gdp') - inner_repn.apply_to(m) - - self.check_pw_log(m) - self.check_pw_paraboloid(m) + ct.check_transformation_do_not_descend(self, + 'contrib.piecewise.inner_repn_gdp') def test_transformation_PiecewiseLinearFunction_targets(self): - m = models.make_log_x_model() - inner_repn = TransformationFactory('contrib.piecewise.inner_repn_gdp') - inner_repn.apply_to(m, targets=[m.pw_log]) - - self.check_pw_log(m) - - # And check that the paraboloid was *not* transformed. - self.assertIsNone( - m.pw_paraboloid.get_transformation_var(m.paraboloid_expr)) + ct.check_transformation_PiecewiseLinearFunction_targets( + self, + 'contrib.piecewise.inner_repn_gdp') def test_descend_into_expressions(self): - m = models.make_log_x_model() - inner_repn = TransformationFactory('contrib.piecewise.inner_repn_gdp') - inner_repn.apply_to(m, descend_into_expressions=True) - - # Everything should be transformed - self.check_pw_log(m) - self.check_pw_paraboloid(m) + ct.check_descend_into_expressions(self, + 'contrib.piecewise.inner_repn_gdp') def test_descend_into_expressions_constraint_target(self): - m = models.make_log_x_model() - inner_repn = TransformationFactory('contrib.piecewise.inner_repn_gdp') - inner_repn.apply_to(m, descend_into_expressions=True, - targets=[m.indexed_c]) - - self.check_pw_paraboloid(m) - # And check that the log was *not* transformed. - self.assertIsNone(m.pw_log.get_transformation_var(m.log_expr)) + ct.check_descend_into_expressions_constraint_target( + self, + 'contrib.piecewise.inner_repn_gdp') def test_descend_into_expressions_objective_target(self): - m = models.make_log_x_model() - inner_repn = TransformationFactory('contrib.piecewise.inner_repn_gdp') - inner_repn.apply_to(m, descend_into_expressions=True, - targets=[m.obj]) - - self.check_pw_log(m) - # And check that the paraboloid was *not* transformed. - self.assertIsNone( - m.pw_paraboloid.get_transformation_var(m.paraboloid_expr)) + ct.check_descend_into_expressions_objective_target( + self, + 'contrib.piecewise.inner_repn_gdp') @unittest.skipUnless(SolverFactory('gurobi').available(), 'Gurobi is not available') @@ -183,10 +157,6 @@ def test_solve_disaggregated_convex_combo_model(self): m = models.make_log_x_model() TransformationFactory( 'contrib.piecewise.disaggregated_convex_combination').apply_to(m) + SolverFactory('gurobi').solve(m) - SolverFactory('gurobi').solve(m, tee=True) - - self.assertAlmostEqual(value(m.x), 4) - self.assertAlmostEqual(value(m.x1), 1) - self.assertAlmostEqual(value(m.x2), 1) - self.assertAlmostEqual(value(m.obj), m.f2(4)) + ct.check_log_x_model_soln(self, m) diff --git a/pyomo/contrib/piecewise/tests/test_outer_repn_gdp.py b/pyomo/contrib/piecewise/tests/test_outer_repn_gdp.py index 3eb44029b32..04a994edd5c 100644 --- a/pyomo/contrib/piecewise/tests/test_outer_repn_gdp.py +++ b/pyomo/contrib/piecewise/tests/test_outer_repn_gdp.py @@ -18,7 +18,7 @@ from pyomo.core.expr.compare import ( assertExpressionsEqual, assertExpressionsStructurallyEqual) from pyomo.gdp import Disjunct, Disjunction -from pyomo.environ import Constraint, SolverFactory, value, Var +from pyomo.environ import Constraint, SolverFactory, Var class TestTransformPiecewiseModelToOuterRepnGDP(unittest.TestCase): def check_log_disjunct(self, d, pts, f, substitute_var, x): @@ -129,55 +129,24 @@ def check_pw_paraboloid(self, m): @unittest.skipUnless(scipy_available, "Scipy is not available") def test_transformation_do_not_descend(self): - m = models.make_log_x_model() - outer_repn = TransformationFactory('contrib.outer_repn_gdp') - outer_repn.apply_to(m) - - self.check_pw_log(m) - self.check_pw_paraboloid(m) + ct.check_transformation_do_not_descend(self, 'contrib.outer_repn_gdp') def test_transformation_PiecewiseLinearFunction_targets(self): - m = models.make_log_x_model() - outer_repn = TransformationFactory('contrib.outer_repn_gdp') - outer_repn.apply_to(m, targets=[m.pw_log]) - - self.check_pw_log(m) - - # And check that the paraboloid was *not* transformed. - self.assertIsNone( - m.pw_paraboloid.get_transformation_var(m.paraboloid_expr)) + ct.check_transformation_PiecewiseLinearFunction_targets( + self, 'contrib.outer_repn_gdp') @unittest.skipUnless(scipy_available, "Scipy is not available") def test_descend_into_expressions(self): - m = models.make_log_x_model() - outer_repn = TransformationFactory('contrib.outer_repn_gdp') - outer_repn.apply_to(m, descend_into_expressions=True) - - # Everything should be transformed - self.check_pw_log(m) - self.check_pw_paraboloid(m) + ct.check_descend_into_expressions(self, 'contrib.outer_repn_gdp') @unittest.skipUnless(scipy_available, "Scipy is not available") def test_descend_into_expressions_constraint_target(self): - m = models.make_log_x_model() - outer_repn = TransformationFactory('contrib.outer_repn_gdp') - outer_repn.apply_to(m, descend_into_expressions=True, - targets=[m.indexed_c]) - - self.check_pw_paraboloid(m) - # And check that the log was *not* transformed. - self.assertIsNone(m.pw_log.get_transformation_var(m.log_expr)) + ct.check_descend_into_expressions_constraint_target( + self, 'contrib.outer_repn_gdp') def test_descend_into_expressions_objective_target(self): - m = models.make_log_x_model() - outer_repn = TransformationFactory('contrib.outer_repn_gdp') - outer_repn.apply_to(m, descend_into_expressions=True, - targets=[m.obj]) - - self.check_pw_log(m) - # And check that the paraboloid was *not* transformed. - self.assertIsNone( - m.pw_paraboloid.get_transformation_var(m.paraboloid_expr)) + ct.check_descend_into_expressions_objective_target( + self, 'contrib.outer_repn_gdp') @unittest.skipUnless(SolverFactory('gurobi').available() and scipy_available, @@ -185,10 +154,6 @@ def test_descend_into_expressions_objective_target(self): def test_solve_multiple_choice_model(self): m = models.make_log_x_model() TransformationFactory('contrib.piecewise.multiple_choice').apply_to(m) + SolverFactory('gurobi').solve(m) - SolverFactory('gurobi').solve(m, tee=True) - - self.assertAlmostEqual(value(m.x), 4) - self.assertAlmostEqual(value(m.x1), 1) - self.assertAlmostEqual(value(m.x2), 1) - self.assertAlmostEqual(value(m.obj), m.f2(4)) + ct.check_log_x_model_soln(self, m) From 9eba234a8ee413151308bb063e7312cfdbb5cbcb Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Fri, 10 Mar 2023 21:37:36 -0700 Subject: [PATCH 07/11] Regestering the outer repn transformation under contrib.piecewise also --- .../contrib/piecewise/tests/test_outer_repn_gdp.py | 13 ++++++++----- .../contrib/piecewise/transform/multiple_choice.py | 3 ++- .../piecewise/transform/outer_representation_gdp.py | 5 +---- 3 files changed, 11 insertions(+), 10 deletions(-) diff --git a/pyomo/contrib/piecewise/tests/test_outer_repn_gdp.py b/pyomo/contrib/piecewise/tests/test_outer_repn_gdp.py index 04a994edd5c..0a22d2c528e 100644 --- a/pyomo/contrib/piecewise/tests/test_outer_repn_gdp.py +++ b/pyomo/contrib/piecewise/tests/test_outer_repn_gdp.py @@ -129,24 +129,27 @@ def check_pw_paraboloid(self, m): @unittest.skipUnless(scipy_available, "Scipy is not available") def test_transformation_do_not_descend(self): - ct.check_transformation_do_not_descend(self, 'contrib.outer_repn_gdp') + ct.check_transformation_do_not_descend( + self, + 'contrib.piecewise.outer_repn_gdp') def test_transformation_PiecewiseLinearFunction_targets(self): ct.check_transformation_PiecewiseLinearFunction_targets( - self, 'contrib.outer_repn_gdp') + self, 'contrib.piecewise.outer_repn_gdp') @unittest.skipUnless(scipy_available, "Scipy is not available") def test_descend_into_expressions(self): - ct.check_descend_into_expressions(self, 'contrib.outer_repn_gdp') + ct.check_descend_into_expressions(self, + 'contrib.piecewise.outer_repn_gdp') @unittest.skipUnless(scipy_available, "Scipy is not available") def test_descend_into_expressions_constraint_target(self): ct.check_descend_into_expressions_constraint_target( - self, 'contrib.outer_repn_gdp') + self, 'contrib.piecewise.outer_repn_gdp') def test_descend_into_expressions_objective_target(self): ct.check_descend_into_expressions_objective_target( - self, 'contrib.outer_repn_gdp') + self, 'contrib.piecewise.outer_repn_gdp') @unittest.skipUnless(SolverFactory('gurobi').available() and scipy_available, diff --git a/pyomo/contrib/piecewise/transform/multiple_choice.py b/pyomo/contrib/piecewise/transform/multiple_choice.py index 9f546311a87..c3da8ceea00 100644 --- a/pyomo/contrib/piecewise/transform/multiple_choice.py +++ b/pyomo/contrib/piecewise/transform/multiple_choice.py @@ -31,5 +31,6 @@ class MultipleChoiceTransformation(Transformation): 2010. """ def _apply_to(self, instance, **kwds): - TransformationFactory('contrib.outer_repn_gdp').apply_to(instance) + TransformationFactory('contrib.piecewise.outer_repn_gdp').apply_to( + instance) TransformationFactory('gdp.hull').apply_to(instance) diff --git a/pyomo/contrib/piecewise/transform/outer_representation_gdp.py b/pyomo/contrib/piecewise/transform/outer_representation_gdp.py index f32304f23d4..5a07f63c2a7 100644 --- a/pyomo/contrib/piecewise/transform/outer_representation_gdp.py +++ b/pyomo/contrib/piecewise/transform/outer_representation_gdp.py @@ -20,10 +20,7 @@ from pyomo.core.util import target_list from pyomo.gdp import Disjunct, Disjunction -## DEBUG -from pytest import set_trace - -@TransformationFactory.register('contrib.outer_repn_gdp', +@TransformationFactory.register('contrib.piecewise.outer_repn_gdp', doc="Convert piecewise-linear model to a GDP " "using an outer (Ax <= b) representation of " "the simplices that are the domains of the " From af67991668538a45944925fbb7d6fc9ec3a8bbb0 Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Sat, 11 Mar 2023 08:09:21 -0700 Subject: [PATCH 08/11] Moving ConfigDict to base class --- .../transform/inner_representation_gdp.py | 33 +------------------ .../transform/outer_representation_gdp.py | 33 +------------------ .../piecewise_to_gdp_transformation.py | 33 +++++++++++++++++++ 3 files changed, 35 insertions(+), 64 deletions(-) diff --git a/pyomo/contrib/piecewise/transform/inner_representation_gdp.py b/pyomo/contrib/piecewise/transform/inner_representation_gdp.py index 9b05def19ba..beaf88ed9aa 100644 --- a/pyomo/contrib/piecewise/transform/inner_representation_gdp.py +++ b/pyomo/contrib/piecewise/transform/inner_representation_gdp.py @@ -9,13 +9,11 @@ # This software is distributed under the 3-clause BSD License. # ___________________________________________________________________________ -from pyomo.common.config import ConfigDict, ConfigValue from pyomo.contrib.fbbt.fbbt import compute_bounds_on_expr from pyomo.contrib.piecewise.transform.piecewise_to_gdp_transformation import ( PiecewiseLinearToGDP) from pyomo.core import Constraint, NonNegativeIntegers, Suffix, Var from pyomo.core.base import TransformationFactory -from pyomo.core.util import target_list from pyomo.gdp import Disjunct, Disjunction @TransformationFactory.register('contrib.piecewise.inner_repn_gdp', @@ -46,36 +44,7 @@ class InnerRepresentationGDPTransformation(PiecewiseLinearToGDP): parent Block as the component owning their parent expression. In this mode, targets must be Blocks, Constraints, and/or Objectives. """ - CONFIG = ConfigDict('piecewise.inner_repn_gdp') - CONFIG.declare('targets', ConfigValue( - default=None, - domain=target_list, - description="target or list of targets that will be transformed", - doc=""" - This specifies the list of components to transform. If None (default), - the entire model is transformed. Note that if the transformation is - done out of place, the list of targets should be attached to the model - before it is cloned, and the list will specify the targets on the cloned - instance.""" - )) - CONFIG.declare('descend_into_expressions', ConfigValue( - default=False, - domain=bool, - description="Whether to look for uses of PiecewiseLinearFunctions in " - "the Constraint and Objective expressions, rather than assuming " - "all PiecewiseLinearFunctions are on the active tree(s) of 'instance' " - "and 'targets.'", - doc=""" - It is *strongly* recommended that, in hierarchical models, the - PiecewiseLinearFunction components are on the same Block as where - they are used in expressions. If you follow this recommendation, - this option can remain False, which will make this transformation - more efficient. However, if you do not follow the recommendation, - unless you know what you are doing, turn this option to 'True' to - ensure that all of the uses of PiecewiseLinearFunctions are - transformed. - """ - )) + CONFIG = PiecewiseLinearToGDP.CONFIG() _transformation_name = 'pw_linear_inner_repn' def _transform_pw_linear_expr(self, pw_expr, pw_linear_func, diff --git a/pyomo/contrib/piecewise/transform/outer_representation_gdp.py b/pyomo/contrib/piecewise/transform/outer_representation_gdp.py index 5a07f63c2a7..8aaf616ea8b 100644 --- a/pyomo/contrib/piecewise/transform/outer_representation_gdp.py +++ b/pyomo/contrib/piecewise/transform/outer_representation_gdp.py @@ -9,7 +9,6 @@ # This software is distributed under the 3-clause BSD License. # ___________________________________________________________________________ -from pyomo.common.config import ConfigDict, ConfigValue import pyomo.common.dependencies.numpy as np from pyomo.common.dependencies.scipy import spatial from pyomo.contrib.fbbt.fbbt import compute_bounds_on_expr @@ -17,7 +16,6 @@ PiecewiseLinearToGDP) from pyomo.core import Constraint, NonNegativeIntegers, Suffix, Var from pyomo.core.base import TransformationFactory -from pyomo.core.util import target_list from pyomo.gdp import Disjunct, Disjunction @TransformationFactory.register('contrib.piecewise.outer_repn_gdp', @@ -46,36 +44,7 @@ class OuterRepresentationGDPTransformation(PiecewiseLinearToGDP): parent Block as the component owning their parent expression. In this mode, targets must be Blocks, Constraints, and/or Objectives. """ - CONFIG = ConfigDict('piecewise.outer_repn_gdp') - CONFIG.declare('targets', ConfigValue( - default=None, - domain=target_list, - description="target or list of targets that will be transformed", - doc=""" - This specifies the list of components to transform. If None (default), - the entire model is transformed. Note that if the transformation is - done out of place, the list of targets should be attached to the model - before it is cloned, and the list will specify the targets on the cloned - instance.""" - )) - CONFIG.declare('descend_into_expressions', ConfigValue( - default=False, - domain=bool, - description="Whether to look for uses of PiecewiseLinearFunctions in " - "the Constraint and Objective expressions, rather than assuming " - "all PiecewiseLinearFunctions are on the active tree(s) of 'instance' " - "and 'targets.'", - doc=""" - It is *strongly* recommended that, in hierarchical models, the - PiecewiseLinearFunction components are on the same Block as where - they are used in expressions. If you follow this recommendation, - this option can remain False, which will make this transformation - more efficient. However, if you do not follow the recommendation, - unless you know what you are doing, turn this option to 'True' to - ensure that all of the uses of PiecewiseLinearFunctions are - transformed. - """ - )) + CONFIG = PiecewiseLinearToGDP.CONFIG() _transformation_name = 'pw_linear_outer_repn' def _transform_pw_linear_expr(self, pw_expr, pw_linear_func, diff --git a/pyomo/contrib/piecewise/transform/piecewise_to_gdp_transformation.py b/pyomo/contrib/piecewise/transform/piecewise_to_gdp_transformation.py index 4b09398f524..4cbf3990eb3 100644 --- a/pyomo/contrib/piecewise/transform/piecewise_to_gdp_transformation.py +++ b/pyomo/contrib/piecewise/transform/piecewise_to_gdp_transformation.py @@ -9,6 +9,7 @@ # This software is distributed under the 3-clause BSD License. # ___________________________________________________________________________ +from pyomo.common.config import ConfigDict, ConfigValue from pyomo.common.errors import DeveloperError from pyomo.common.modeling import unique_component_name from pyomo.contrib.piecewise import PiecewiseLinearFunction @@ -19,6 +20,7 @@ SetOf, RangeSet, ExternalFunction, Connector, SortComponents, Any) from pyomo.core.base import Transformation from pyomo.core.base.block import _BlockData, Block +from pyomo.core.util import target_list from pyomo.gdp import Disjunct, Disjunction from pyomo.gdp.util import is_child_of from pyomo.network import Port @@ -27,6 +29,37 @@ class PiecewiseLinearToGDP(Transformation): """ Base class for transformations of piecewise-linear models to GDPs """ + CONFIG = ConfigDict('contrib.piecewise_to_gdp') + CONFIG.declare('targets', ConfigValue( + default=None, + domain=target_list, + description="target or list of targets that will be transformed", + doc=""" + This specifies the list of components to transform. If None (default), + the entire model is transformed. Note that if the transformation is + done out of place, the list of targets should be attached to the model + before it is cloned, and the list will specify the targets on the cloned + instance.""" + )) + CONFIG.declare('descend_into_expressions', ConfigValue( + default=False, + domain=bool, + description="Whether to look for uses of PiecewiseLinearFunctions in " + "the Constraint and Objective expressions, rather than assuming " + "all PiecewiseLinearFunctions are on the active tree(s) of 'instance' " + "and 'targets.'", + doc=""" + It is *strongly* recommended that, in hierarchical models, the + PiecewiseLinearFunction components are on the same Block as where + they are used in expressions. If you follow this recommendation, + this option can remain False, which will make this transformation + more efficient. However, if you do not follow the recommendation, + unless you know what you are doing, turn this option to 'True' to + ensure that all of the uses of PiecewiseLinearFunctions are + transformed. + """ + )) + def __init__(self): super().__init__() self.handlers = { From d1b235b9c877888c03d321ad35448f8a729d1d9b Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Mon, 13 Mar 2023 11:08:58 -0600 Subject: [PATCH 09/11] Adding note about the simplification to get Vielma et al.'s multiple choice model from our multiple choice transformation --- pyomo/contrib/piecewise/transform/multiple_choice.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/pyomo/contrib/piecewise/transform/multiple_choice.py b/pyomo/contrib/piecewise/transform/multiple_choice.py index c3da8ceea00..fcc727c5510 100644 --- a/pyomo/contrib/piecewise/transform/multiple_choice.py +++ b/pyomo/contrib/piecewise/transform/multiple_choice.py @@ -21,7 +21,11 @@ class MultipleChoiceTransformation(Transformation): Converts a model containing PiecewiseLinearFunctions to a an equivalent MIP via the Multiple Choice method from [1]. Note that, while this model probably resolves to the model described in [1] after - presolve, the Pyomo version is not as simplified. + presolve, the Pyomo version is not as simplified. Specifically, in [1], the + the 'z' variables (representing the value of the piecewise-linear function + in each Disjunct) are not disaggregated. In this transformation's output + they will be, but a linear combination of inequalities yields a model + equivalent to the Multiple Choice model in [1]. References ---------- From 4f053dbb24061fb2659d18f2446a6ba96a85fefb Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Mon, 13 Mar 2023 11:17:27 -0600 Subject: [PATCH 10/11] Fixing a comment --- pyomo/contrib/piecewise/transform/outer_representation_gdp.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pyomo/contrib/piecewise/transform/outer_representation_gdp.py b/pyomo/contrib/piecewise/transform/outer_representation_gdp.py index 8aaf616ea8b..d1872c379fc 100644 --- a/pyomo/contrib/piecewise/transform/outer_representation_gdp.py +++ b/pyomo/contrib/piecewise/transform/outer_representation_gdp.py @@ -78,8 +78,8 @@ def _transform_pw_linear_expr(self, pw_expr, pw_linear_func, pw_linear_func._points[simplex[1]][0])) else: disj.simplex_halfspaces = Constraint(range(dimension + 1)) - # we will use scipy to get the convex hull of the extreme pts of - # the simplex + # we will use scipy to get the convex hull of the extreme + # points of the simplex extreme_pts = [] for idx in simplex: extreme_pts.append(pw_linear_func._points[idx]) From 8164d7af7d8a281fa16b6de42f369a87c416a286 Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Mon, 13 Mar 2023 11:22:12 -0600 Subject: [PATCH 11/11] Just passing the log function directly in the test model --- pyomo/contrib/piecewise/tests/models.py | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/pyomo/contrib/piecewise/tests/models.py b/pyomo/contrib/piecewise/tests/models.py index 9639a95fbb8..1f326eb1793 100644 --- a/pyomo/contrib/piecewise/tests/models.py +++ b/pyomo/contrib/piecewise/tests/models.py @@ -15,11 +15,7 @@ def make_log_x_model(): m = ConcreteModel() m.x = Var(bounds=(1, 10)) - def log_function(x): - return log(x) - m.log_function = log_function - m.pw_log = PiecewiseLinearFunction(points=[1, 3, 6, 10], - function=m.log_function) + m.pw_log = PiecewiseLinearFunction(points=[1, 3, 6, 10], function=log) # Here are the linear functions, for safe keeping. def f1(x):