diff --git a/pyomo/gdp/plugins/__init__.py b/pyomo/gdp/plugins/__init__.py index 1222ce500f1..2edb99bbe1b 100644 --- a/pyomo/gdp/plugins/__init__.py +++ b/pyomo/gdp/plugins/__init__.py @@ -22,3 +22,4 @@ def load(): import pyomo.gdp.plugins.multiple_bigm import pyomo.gdp.plugins.transform_current_disjunctive_state import pyomo.gdp.plugins.bound_pretransformation + import pyomo.gdp.plugins.binary_multiplication diff --git a/pyomo/gdp/plugins/binary_multiplication.py b/pyomo/gdp/plugins/binary_multiplication.py new file mode 100644 index 00000000000..5afb661aaa8 --- /dev/null +++ b/pyomo/gdp/plugins/binary_multiplication.py @@ -0,0 +1,177 @@ +# ___________________________________________________________________________ +# +# 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 .gdp_to_mip_transformation import GDP_to_MIP_Transformation +from pyomo.common.config import ConfigDict, ConfigValue +from pyomo.core.base import TransformationFactory +from pyomo.core.util import target_list +from pyomo.gdp import Disjunction +from weakref import ref as weakref_ref +import logging + + +logger = logging.getLogger(__name__) + + +@TransformationFactory.register( + 'gdp.binary_multiplication', + doc="Reformulate the GDP as an MINLP by multiplying f(x) <= 0 by y to get " + "f(x) * y <= 0 where y is the binary corresponding to the Boolean indicator " + "var of the Disjunct containing f(x) <= 0.", +) +class GDPBinaryMultiplicationTransformation(GDP_to_MIP_Transformation): + CONFIG = ConfigDict("gdp.binary_multiplication") + 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.""", + ), + ) + + transformation_name = 'binary_multiplication' + + def __init__(self): + super().__init__(logger) + + def _apply_to(self, instance, **kwds): + try: + self._apply_to_impl(instance, **kwds) + finally: + self._restore_state() + + def _apply_to_impl(self, instance, **kwds): + self._process_arguments(instance, **kwds) + + # filter out inactive targets and handle case where targets aren't + # specified. + targets = self._filter_targets(instance) + # transform logical constraints based on targets + self._transform_logical_constraints(instance, targets) + # we need to preprocess targets to make sure that if there are any + # disjunctions in targets that their disjuncts appear before them in + # the list. + gdp_tree = self._get_gdp_tree_from_targets(instance, targets) + preprocessed_targets = gdp_tree.reverse_topological_sort() + + for t in preprocessed_targets: + if t.ctype is Disjunction: + self._transform_disjunctionData( + t, + t.index(), + parent_disjunct=gdp_tree.parent(t), + root_disjunct=gdp_tree.root_disjunct(t), + ) + + def _transform_disjunctionData( + self, obj, index, parent_disjunct=None, root_disjunct=None + ): + (transBlock, xorConstraint) = self._setup_transform_disjunctionData( + obj, root_disjunct + ) + + # add or (or xor) constraint + or_expr = 0 + for disjunct in obj.disjuncts: + or_expr += disjunct.binary_indicator_var + self._transform_disjunct(disjunct, transBlock) + + rhs = 1 if parent_disjunct is None else parent_disjunct.binary_indicator_var + if obj.xor: + xorConstraint[index] = or_expr == rhs + else: + xorConstraint[index] = or_expr >= rhs + # Mark the DisjunctionData as transformed by mapping it to its XOR + # constraint. + obj._algebraic_constraint = weakref_ref(xorConstraint[index]) + + # and deactivate for the writers + obj.deactivate() + + def _transform_disjunct(self, obj, transBlock): + # We're not using the preprocessed list here, so this could be + # inactive. We've already done the error checking in preprocessing, so + # we just skip it here. + if not obj.active: + return + + relaxationBlock = self._get_disjunct_transformation_block(obj, transBlock) + + # Transform each component within this disjunct + self._transform_block_components(obj, obj) + + # deactivate disjunct to keep the writers happy + obj._deactivate_without_fixing_indicator() + + def _transform_constraint(self, obj, disjunct): + # add constraint to the transformation block, we'll transform it there. + transBlock = disjunct._transformation_block() + constraintMap = transBlock._constraintMap + + disjunctionRelaxationBlock = transBlock.parent_block() + + # We will make indexes from ({obj.local_name} x obj.index_set() x ['lb', + # 'ub']), but don't bother construct that set here, as taking Cartesian + # products is kind of expensive (and redundant since we have the + # original model) + newConstraint = transBlock.transformedConstraints + + for i in sorted(obj.keys()): + c = obj[i] + if not c.active: + continue + + self._add_constraint_expressions( + c, i, disjunct.binary_indicator_var, newConstraint, constraintMap + ) + + # deactivate because we relaxed + c.deactivate() + + def _add_constraint_expressions( + self, c, i, indicator_var, newConstraint, constraintMap + ): + # Since we are both combining components from multiple blocks and using + # local names, we need to make sure that the first index for + # transformedConstraints is guaranteed to be unique. We just grab the + # current length of the list here since that will be monotonically + # increasing and hence unique. We'll append it to the + # slightly-more-human-readable constraint name for something familiar + # but unique. (Note that we really could do this outside of the loop + # over the constraint indices, but I don't think it matters a lot.) + unique = len(newConstraint) + name = c.local_name + "_%s" % unique + transformed = constraintMap['transformedConstraints'][c] = [] + + lb, ub = c.lower, c.upper + if (c.equality or lb is ub) and lb is not None: + # equality + newConstraint.add((name, i, 'eq'), (c.body - lb) * indicator_var == 0) + transformed.append(newConstraint[name, i, 'eq']) + constraintMap['srcConstraints'][newConstraint[name, i, 'eq']] = c + else: + # inequality + if lb is not None: + newConstraint.add((name, i, 'lb'), 0 <= (c.body - lb) * indicator_var) + transformed.append(newConstraint[name, i, 'lb']) + constraintMap['srcConstraints'][newConstraint[name, i, 'lb']] = c + if ub is not None: + newConstraint.add((name, i, 'ub'), (c.body - ub) * indicator_var <= 0) + transformed.append(newConstraint[name, i, 'ub']) + constraintMap['srcConstraints'][newConstraint[name, i, 'ub']] = c diff --git a/pyomo/gdp/tests/common_tests.py b/pyomo/gdp/tests/common_tests.py index b475334981b..4a772a7ae56 100644 --- a/pyomo/gdp/tests/common_tests.py +++ b/pyomo/gdp/tests/common_tests.py @@ -58,6 +58,24 @@ def check_linear_coef(self, repn, var, coef): self.assertAlmostEqual(repn.linear_coefs[var_id], coef) +def check_quadratic_coef(self, repn, v1, v2, coef): + if isinstance(v1, BooleanVar): + v1 = v1.get_associated_binary() + if isinstance(v2, BooleanVar): + v2 = v2.get_associated_binary() + + v1id = id(v1) + v2id = id(v2) + + qcoef_map = dict() + for (_v1, _v2), _coef in zip(repn.quadratic_vars, repn.quadratic_coefs): + qcoef_map[id(_v1), id(_v2)] = _coef + qcoef_map[id(_v2), id(_v1)] = _coef + + self.assertIn((v1id, v2id), qcoef_map) + self.assertAlmostEqual(qcoef_map[v1id, v2id], coef) + + def check_squared_term_coef(self, repn, var, coef): var_id = None for i, (v1, v2) in enumerate(repn.quadratic_vars): diff --git a/pyomo/gdp/tests/test_binary_multiplication.py b/pyomo/gdp/tests/test_binary_multiplication.py new file mode 100644 index 00000000000..6b3ba87fa21 --- /dev/null +++ b/pyomo/gdp/tests/test_binary_multiplication.py @@ -0,0 +1,301 @@ +# ___________________________________________________________________________ +# +# 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. +# ___________________________________________________________________________ + +import pyomo.common.unittest as unittest + +from pyomo.environ import ( + TransformationFactory, + Block, + Constraint, + ConcreteModel, + Var, + Any, +) +from pyomo.gdp import Disjunct, Disjunction +from pyomo.core.expr.compare import assertExpressionsEqual +from pyomo.repn import generate_standard_repn +from pyomo.core.expr.compare import assertExpressionsEqual + +import pyomo.core.expr as EXPR +import pyomo.gdp.tests.models as models +import pyomo.gdp.tests.common_tests as ct + +import random + + +class CommonTests: + def diff_apply_to_and_create_using(self, model): + ct.diff_apply_to_and_create_using(self, model, 'gdp.binary_multiplication') + + +class TwoTermDisj(unittest.TestCase, CommonTests): + def setUp(self): + # set seed so we can test name collisions predictably + random.seed(666) + + def test_new_block_created(self): + m = models.makeTwoTermDisj() + TransformationFactory('gdp.binary_multiplication').apply_to(m) + + # we have a transformation block + transBlock = m.component("_pyomo_gdp_binary_multiplication_reformulation") + self.assertIsInstance(transBlock, Block) + + disjBlock = transBlock.component("relaxedDisjuncts") + self.assertIsInstance(disjBlock, Block) + self.assertEqual(len(disjBlock), 2) + # it has the disjuncts on it + self.assertIs(m.d[0].transformation_block, disjBlock[0]) + self.assertIs(m.d[1].transformation_block, disjBlock[1]) + + def test_disjunction_deactivated(self): + ct.check_disjunction_deactivated(self, 'binary_multiplication') + + def test_disjunctDatas_deactivated(self): + ct.check_disjunctDatas_deactivated(self, 'binary_multiplication') + + def test_do_not_transform_twice_if_disjunction_reactivated(self): + ct.check_do_not_transform_twice_if_disjunction_reactivated( + self, 'binary_multiplication' + ) + + def test_xor_constraint_mapping(self): + ct.check_xor_constraint_mapping(self, 'binary_multiplication') + + def test_xor_constraint_mapping_two_disjunctions(self): + ct.check_xor_constraint_mapping_two_disjunctions(self, 'binary_multiplication') + + def test_disjunct_mapping(self): + ct.check_disjunct_mapping(self, 'binary_multiplication') + + def test_disjunct_and_constraint_maps(self): + """Tests the actual data structures used to store the maps.""" + m = models.makeTwoTermDisj() + binary_multiplication = TransformationFactory('gdp.binary_multiplication') + binary_multiplication.apply_to(m) + disjBlock = m._pyomo_gdp_binary_multiplication_reformulation.relaxedDisjuncts + oldblock = m.component("d") + + # we are counting on the fact that the disjuncts get relaxed in the + # same order every time. + for i in [0, 1]: + self.assertIs(oldblock[i].transformation_block, disjBlock[i]) + self.assertIs( + binary_multiplication.get_src_disjunct(disjBlock[i]), oldblock[i] + ) + + # check constraint dict has right mapping + c1_list = binary_multiplication.get_transformed_constraints(oldblock[1].c1) + # this is an equality + self.assertEqual(len(c1_list), 1) + self.assertIs(c1_list[0].parent_block(), disjBlock[1]) + self.assertIs( + binary_multiplication.get_src_constraint(c1_list[0]), oldblock[1].c1 + ) + + c2_list = binary_multiplication.get_transformed_constraints(oldblock[1].c2) + # just ub + self.assertEqual(len(c2_list), 1) + self.assertIs(c2_list[0].parent_block(), disjBlock[1]) + self.assertIs( + binary_multiplication.get_src_constraint(c2_list[0]), oldblock[1].c2 + ) + + c_list = binary_multiplication.get_transformed_constraints(oldblock[0].c) + # just lb + self.assertEqual(len(c_list), 1) + self.assertIs(c_list[0].parent_block(), disjBlock[0]) + self.assertIs( + binary_multiplication.get_src_constraint(c_list[0]), oldblock[0].c + ) + + def test_new_block_nameCollision(self): + ct.check_transformation_block_name_collision(self, 'binary_multiplication') + + def test_indicator_vars(self): + ct.check_indicator_vars(self, 'binary_multiplication') + + def test_xor_constraints(self): + ct.check_xor_constraint(self, 'binary_multiplication') + + def test_or_constraints(self): + m = models.makeTwoTermDisj() + m.disjunction.xor = False + TransformationFactory('gdp.binary_multiplication').apply_to(m) + + # check or constraint is an or (upper bound is None) + orcons = m._pyomo_gdp_binary_multiplication_reformulation.component( + "disjunction_xor" + ) + self.assertIsInstance(orcons, Constraint) + assertExpressionsEqual( + self, + orcons.body, + EXPR.LinearExpression( + [ + EXPR.MonomialTermExpression((1, m.d[0].binary_indicator_var)), + EXPR.MonomialTermExpression((1, m.d[1].binary_indicator_var)), + ] + ), + ) + self.assertEqual(orcons.lower, 1) + self.assertIsNone(orcons.upper) + + def test_deactivated_constraints(self): + ct.check_deactivated_constraints(self, 'binary_multiplication') + + def test_transformed_constraints(self): + m = models.makeTwoTermDisj() + binary_multiplication = TransformationFactory('gdp.binary_multiplication') + binary_multiplication.apply_to(m) + self.check_transformed_constraints(m, binary_multiplication, -3, 2, 7, 2) + + def test_do_not_transform_userDeactivated_disjuncts(self): + ct.check_user_deactivated_disjuncts(self, 'binary_multiplication') + + def test_improperly_deactivated_disjuncts(self): + ct.check_improperly_deactivated_disjuncts(self, 'binary_multiplication') + + def test_do_not_transform_userDeactivated_IndexedDisjunction(self): + ct.check_do_not_transform_userDeactivated_indexedDisjunction( + self, 'binary_multiplication' + ) + + def check_transformed_constraints( + self, model, binary_multiplication, cons1lb, cons2lb, cons2ub, cons3ub + ): + disjBlock = ( + model._pyomo_gdp_binary_multiplication_reformulation.relaxedDisjuncts + ) + + # first constraint + c = binary_multiplication.get_transformed_constraints(model.d[0].c) + self.assertEqual(len(c), 1) + c_lb = c[0] + self.assertTrue(c[0].active) + ind_var = model.d[0].indicator_var + assertExpressionsEqual( + self, c[0].body, (model.a - model.d[0].c.lower) * ind_var + ) + self.assertEqual(c[0].lower, 0) + self.assertIsNone(c[0].upper) + + # second constraint + c = binary_multiplication.get_transformed_constraints(model.d[1].c1) + self.assertEqual(len(c), 1) + c_eq = c[0] + self.assertTrue(c[0].active) + ind_var = model.d[1].indicator_var + assertExpressionsEqual(self, c[0].body, model.a * ind_var) + self.assertEqual(c[0].lower, 0) + self.assertEqual(c[0].upper, 0) + + # third constraint + c = binary_multiplication.get_transformed_constraints(model.d[1].c2) + self.assertEqual(len(c), 1) + c_ub = c[0] + self.assertTrue(c_ub.active) + assertExpressionsEqual( + self, c_ub.body, (model.x - model.d[1].c2.upper) * ind_var + ) + self.assertIsNone(c_ub.lower) + self.assertEqual(c_ub.upper, 0) + + def test_create_using(self): + m = models.makeTwoTermDisj() + self.diff_apply_to_and_create_using(m) + + def test_indexed_constraints_in_disjunct(self): + m = ConcreteModel() + m.I = [1, 2, 3] + m.x = Var(m.I, bounds=(0, 10)) + + def c_rule(b, i): + m = b.model() + return m.x[i] >= i + + def d_rule(d, j): + m = d.model() + d.c = Constraint(m.I[:j], rule=c_rule) + + m.d = Disjunct(m.I, rule=d_rule) + m.disjunction = Disjunction(expr=[m.d[i] for i in m.I]) + + TransformationFactory('gdp.binary_multiplication').apply_to(m) + transBlock = m._pyomo_gdp_binary_multiplication_reformulation + + # 2 blocks: the original Disjunct and the transformation block + self.assertEqual(len(list(m.component_objects(Block, descend_into=False))), 1) + self.assertEqual(len(list(m.component_objects(Disjunct))), 1) + + # Each relaxed disjunct should have 1 var (the reference to the + # indicator var), and i "d[i].c" Constraints + for i in [1, 2, 3]: + relaxed = transBlock.relaxedDisjuncts[i - 1] + self.assertEqual(len(list(relaxed.component_objects(Var))), 1) + self.assertEqual(len(list(relaxed.component_data_objects(Var))), 1) + self.assertEqual(len(list(relaxed.component_objects(Constraint))), 1) + self.assertEqual(len(list(relaxed.component_data_objects(Constraint))), i) + + def test_virtual_indexed_constraints_in_disjunct(self): + m = ConcreteModel() + m.I = [1, 2, 3] + m.x = Var(m.I, bounds=(0, 10)) + + def d_rule(d, j): + m = d.model() + d.c = Constraint(Any) + for k in range(j): + d.c[k + 1] = m.x[k + 1] >= k + 1 + + m.d = Disjunct(m.I, rule=d_rule) + m.disjunction = Disjunction(expr=[m.d[i] for i in m.I]) + + TransformationFactory('gdp.binary_multiplication').apply_to(m) + transBlock = m._pyomo_gdp_binary_multiplication_reformulation + + # 2 blocks: the original Disjunct and the transformation block + self.assertEqual(len(list(m.component_objects(Block, descend_into=False))), 1) + self.assertEqual(len(list(m.component_objects(Disjunct))), 1) + + # Each relaxed disjunct should have 1 var (the reference to the + # indicator var), and i "d[i].c" Constraints + for i in [1, 2, 3]: + relaxed = transBlock.relaxedDisjuncts[i - 1] + self.assertEqual(len(list(relaxed.component_objects(Var))), 1) + self.assertEqual(len(list(relaxed.component_data_objects(Var))), 1) + self.assertEqual(len(list(relaxed.component_objects(Constraint))), 1) + self.assertEqual(len(list(relaxed.component_data_objects(Constraint))), i) + + def test_local_var(self): + m = models.localVar() + binary_multiplication = TransformationFactory('gdp.binary_multiplication') + binary_multiplication.apply_to(m) + + # we just need to make sure that constraint was transformed correctly, + # which just means that the M values were correct. + transformedC = binary_multiplication.get_transformed_constraints(m.disj2.cons) + self.assertEqual(len(transformedC), 1) + eq = transformedC[0] + repn = generate_standard_repn(eq.body) + self.assertIsNone(repn.nonlinear_expr) + self.assertEqual(len(repn.linear_coefs), 1) + self.assertEqual(len(repn.quadratic_coefs), 2) + ct.check_linear_coef(self, repn, m.disj2.indicator_var, -3) + ct.check_quadratic_coef(self, repn, m.x, m.disj2.indicator_var, 1) + ct.check_quadratic_coef(self, repn, m.disj2.y, m.disj2.indicator_var, 1) + self.assertEqual(repn.constant, 0) + self.assertEqual(eq.lb, 0) + self.assertEqual(eq.ub, 0) + + +if __name__ == '__main__': + unittest.main()