Merge pull request Pyomo#3009 from emma58/fix-2702

Various bug fixes in GDP transformations

pyomo/gdp/plugins/bigm.py

@@ -195,27 +195,27 @@ def _apply_to_impl(self, instance, **kwds):
                 self._transform_disjunctionData(
                     t,
                     t.index(),
+                    bigM,
                     parent_disjunct=gdp_tree.parent(t),
                     root_disjunct=gdp_tree.root_disjunct(t),
                 )
-            else:  # We know t is a Disjunct after preprocessing
-                self._transform_disjunct(
-                    t, bigM, root_disjunct=gdp_tree.root_disjunct(t)
-                )
 
         # issue warnings about anything that was in the bigM args dict that we
         # didn't use
         _warn_for_unused_bigM_args(bigM, self.used_args, logger)
 
     def _transform_disjunctionData(
-        self, obj, index, parent_disjunct=None, root_disjunct=None
+        self, obj, index, bigM, parent_disjunct=None, root_disjunct=None
     ):
         (transBlock, xorConstraint) = self._setup_transform_disjunctionData(
             obj, root_disjunct
         )
 
         # add or (or xor) constraint
-        or_expr = sum(disjunct.binary_indicator_var for disjunct in obj.disjuncts)
+        or_expr = 0
+        for disjunct in obj.disjuncts:
+            or_expr += disjunct.binary_indicator_var
+            self._transform_disjunct(disjunct, bigM, transBlock)
 
         rhs = 1 if parent_disjunct is None else parent_disjunct.binary_indicator_var
         if obj.xor:
@@ -229,13 +229,13 @@ def _transform_disjunctionData(
         # and deactivate for the writers
         obj.deactivate()
 
-    def _transform_disjunct(self, obj, bigM, root_disjunct):
-        root = (
-            root_disjunct.parent_block()
-            if root_disjunct is not None
-            else obj.parent_block()
-        )
-        transBlock = self._add_transformation_block(root)[0]
+    def _transform_disjunct(self, obj, bigM, 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
+
         suffix_list = _get_bigM_suffix_list(obj)
         arg_list = self._get_bigM_arg_list(bigM, obj)
 

pyomo/gdp/plugins/multiple_bigm.py

@@ -477,10 +477,9 @@ def _transform_bound_constraints(self, active_disjuncts, transBlock, Ms):
         # Now we actually construct the constraints. We do this separately so
         # that we can make sure that we have a term for every active disjunct in
         # the disjunction (falling back on the variable bounds if they are there
-        transformed = transBlock.transformed_bound_constraints = Constraint(
-            NonNegativeIntegers, ['lb', 'ub']
-        )
-        for idx, (v, (lower_dict, upper_dict)) in enumerate(bounds_cons.items()):
+        transformed = transBlock.transformed_bound_constraints
+        offset = len(transformed)
+        for i, (v, (lower_dict, upper_dict)) in enumerate(bounds_cons.items()):
             lower_rhs = 0
             upper_rhs = 0
             for disj in active_disjuncts:
@@ -515,6 +514,7 @@ def _transform_bound_constraints(self, active_disjuncts, transBlock, Ms):
                             "one of these." % (v.name, disj.name)
                         )
                     upper_rhs += M * disj.indicator_var.get_associated_binary()
+            idx = i + offset
             if len(lower_dict) > 0:
                 transformed.add((idx, 'lb'), v >= lower_rhs)
                 relaxationBlock._constraintMap['srcConstraints'][
@@ -559,6 +559,9 @@ def _add_transformation_block(self, to_block):
         if new_block:
             # Will store M values as we transform
             transBlock._mbm_values = {}
+            transBlock.transformed_bound_constraints = Constraint(
+                NonNegativeIntegers, ['lb', 'ub']
+            )
         return transBlock, new_block
 
     def _get_all_var_objects(self, active_disjuncts):

pyomo/gdp/tests/common_tests.py

@@ -1865,3 +1865,15 @@ def check_transformed_model_pickles_with_dill(self, transformation):
     unpickle = dill.loads(dill.dumps(m))
 
     check_pprint_equal(self, m, unpickle)
+
+
+def check_nested_disjuncts_in_flat_gdp(self, transformation):
+    m = models.make_non_nested_model_declaring_Disjuncts_on_each_other()
+    TransformationFactory('gdp.%s' % transformation).apply_to(m)
+    SolverFactory('gurobi').solve(m)
+    self.assertAlmostEqual(value(m.obj), 1020)
+
+    # check the Boolean solution
+    for t in m.T:
+        self.assertTrue(value(m.disj1[t].indicator_var))
+        self.assertTrue(value(m.disj1[t].sub1.indicator_var))

pyomo/gdp/tests/models.py

@@ -11,6 +11,8 @@
     Any,
     Expression,
     maximize,
+    minimize,
+    NonNegativeReals,
     TransformationFactory,
     BooleanVar,
     LogicalConstraint,
@@ -1137,3 +1139,69 @@ def makeBooleanVarsOnDisjuncts():
     m.bwahaha = LogicalConstraint(expr=m.d[1].Y[1].xor(m.d[1].Y[2]))
 
     return m
+
+
+def make_non_nested_model_declaring_Disjuncts_on_each_other():
+    """
+    T = {1, 2, ..., 10}
+
+    min  sum(x_t + y_t for t in T)
+
+    s.t. 1 <= x_t <= 10, for all t in T
+         1 <= y_t <= 100, for all t in T
+
+         [y_t = 100] v [y_t = 1000], for all t in T
+         [x_t = 2] v [y_t = 10], for all t in T.
+
+
+    We can't choose y_t = 10 because then the first Disjunction is infeasible.
+    so in the optimal solution we choose x_t = 2 and y_t = 100 for all t in T.
+    That gives us an optimal value of (100 + 2)*10 = 1020.
+    """
+    model = ConcreteModel()
+    model.T = RangeSet(10)
+    model.x = Var(model.T, bounds=(1, 10))
+    model.y = Var(model.T, bounds=(1, 100))
+
+    def _op_mode_sub(m, t):
+        m.disj1[t].c1 = Constraint(expr=m.x[t] == 2)
+        m.disj1[t].sub1 = Disjunct()
+        m.disj1[t].sub1.c1 = Constraint(expr=m.y[t] == 100)
+        m.disj1[t].sub2 = Disjunct()
+        m.disj1[t].sub2.c1 = Constraint(expr=m.y[t] == 1000)
+        return [m.disj1[t].sub1, m.disj1[t].sub2]
+
+    def _op_mode(m, t):
+        m.disj2[t].c1 = Constraint(expr=m.y[t] == 10)
+        return [m.disj1[t], m.disj2[t]]
+
+    model.disj1 = Disjunct(model.T)
+    model.disj2 = Disjunct(model.T)
+    model.disjunction1sub = Disjunction(model.T, rule=_op_mode_sub)
+    model.disjunction1 = Disjunction(model.T, rule=_op_mode)
+
+    def obj_rule(m, t):
+        return sum(m.x[t] + m.y[t] for t in m.T)
+
+    model.obj = Objective(rule=obj_rule)
+
+    return model
+
+
+def make_indexed_equality_model():
+    """
+    min  x_1 + x_2
+    s.t. [x_1 = 1] v [x_1 = 2]
+         [x_2 = 1] v [x_2 = 2]
+    """
+
+    def disj_rule(m, t):
+        return [[m.x[t] == 1], [m.x[t] == 2]]
+
+    m = ConcreteModel()
+    m.T = RangeSet(2)
+    m.x = Var(m.T, within=NonNegativeReals, bounds=(0, 5))
+    m.d = Disjunction(m.T, rule=disj_rule)
+    m.obj = Objective(expr=m.x[1] + m.x[2], sense=minimize)
+
+    return m

pyomo/gdp/tests/test_bigm.py

@@ -19,6 +19,7 @@
     Set,
     Constraint,
     ComponentMap,
+    SolverFactory,
     Suffix,
     ConcreteModel,
     Var,
@@ -45,6 +46,11 @@
 
 from io import StringIO
 
+gurobi_available = (
+    SolverFactory('gurobi').available(exception_flag=False)
+    and SolverFactory('gurobi').license_is_valid()
+)
+
 
 class CommonTests:
     def diff_apply_to_and_create_using(self, model):
@@ -1729,8 +1735,7 @@ def check_disjunction_transformation_block_structure(self, transBlock, pairs):
         self.assertEqual(len(disjBlock), len(pairs))
 
         # This test will also rely on the disjunctions being relaxed in the same
-        # order every time (and moved up to the new transformation block in the
-        # same order)
+        # order every time
         bigm = TransformationFactory('gdp.bigm')
         for i, j in pairs:
             for comp in j:
@@ -1750,9 +1755,9 @@ def test_transformation_block_structure(self):
             (1, [m.simpledisjunct.innerdisjunct1.c]),
             (2, []),  # No constraints, just a reference to simpledisjunct's
             # indicator_var
-            (3, [m.disjunct[0].c]),
-            (4, [m.disjunct[1].innerdisjunct[0].c]),
-            (5, [m.disjunct[1].innerdisjunct[1].c]),
+            (5, [m.disjunct[0].c]),
+            (2, [m.disjunct[1].innerdisjunct[0].c]),
+            (3, [m.disjunct[1].innerdisjunct[1].c]),
             (6, []),  # Again no constraints, just indicator var ref
         ]
         self.check_disjunction_transformation_block_structure(transBlock, pairs)
@@ -1792,14 +1797,14 @@ def test_disjunct_mappings(self):
         # I want to check that I correctly updated the pointers to the
         # transformation blocks on the inner Disjuncts.
         self.assertIs(
-            m.disjunct[1].innerdisjunct[0].transformation_block, disjunctBlocks[4]
+            m.disjunct[1].innerdisjunct[0].transformation_block, disjunctBlocks[2]
         )
-        self.assertIs(disjunctBlocks[4]._src_disjunct(), m.disjunct[1].innerdisjunct[0])
+        self.assertIs(disjunctBlocks[2]._src_disjunct(), m.disjunct[1].innerdisjunct[0])
 
         self.assertIs(
-            m.disjunct[1].innerdisjunct[1].transformation_block, disjunctBlocks[5]
+            m.disjunct[1].innerdisjunct[1].transformation_block, disjunctBlocks[3]
         )
-        self.assertIs(disjunctBlocks[5]._src_disjunct(), m.disjunct[1].innerdisjunct[1])
+        self.assertIs(disjunctBlocks[3]._src_disjunct(), m.disjunct[1].innerdisjunct[1])
 
         self.assertIs(
             m.simpledisjunct.innerdisjunct0.transformation_block, disjunctBlocks[0]
@@ -2844,5 +2849,15 @@ def test_dill_pickle(self):
         ct.check_transformed_model_pickles_with_dill(self, 'bigm')
 
 
+@unittest.skipUnless(gurobi_available, "Gurobi is not available")
+class NestedDisjunctsInFlatGDP(unittest.TestCase):
+    """
+    This class tests the fix for #2702
+    """
+
+    def test_declare_disjuncts_in_disjunction_rule(self):
+        ct.check_nested_disjuncts_in_flat_gdp(self, 'bigm')
+
+
 if __name__ == '__main__':
     unittest.main()

pyomo/gdp/tests/test_hull.py

@@ -57,6 +57,11 @@
 EPS = TransformationFactory('gdp.hull').CONFIG.EPS
 linear_solvers = ct.linear_solvers
 
+gurobi_available = (
+    SolverFactory('gurobi').available(exception_flag=False)
+    and SolverFactory('gurobi').license_is_valid()
+)
+
 
 class CommonTests:
     def setUp(self):
@@ -2695,3 +2700,13 @@ def test_pickle(self):
     @unittest.skipIf(not dill_available, "Dill is not available")
     def test_dill_pickle(self):
         ct.check_transformed_model_pickles_with_dill(self, 'hull')
+
+
+@unittest.skipUnless(gurobi_available, "Gurobi is not available")
+class NestedDisjunctsInFlatGDP(unittest.TestCase):
+    """
+    This class tests the fix for #2702
+    """
+
+    def test_declare_disjuncts_in_disjunction_rule(self):
+        ct.check_nested_disjuncts_in_flat_gdp(self, 'hull')

pyomo/gdp/tests/test_mbigm.py

@@ -36,9 +36,11 @@
 from pyomo.gdp import Disjunct, Disjunction, GDP_Error
 from pyomo.gdp.tests.common_tests import (
     check_linear_coef,
+    check_nested_disjuncts_in_flat_gdp,
     check_obj_in_active_tree,
     check_pprint_equal,
 )
+from pyomo.gdp.tests.models import make_indexed_equality_model
 from pyomo.repn import generate_standard_repn
 
 gurobi_available = (
@@ -863,3 +865,68 @@ def test_only_multiple_bigm_bound_constraints_arg_Ms(self):
         self.check_traditionally_bigmed_constraints(
             m, mbm, {m.d1: (-1050, 1050), m.d2: (-2000, 1200), m.d3: (-4000, 4000)}
         )
+
+
+@unittest.skipUnless(gurobi_available, "Gurobi is not available")
+class NestedDisjunctsInFlatGDP(unittest.TestCase):
+    """
+    This class tests the fix for #2702
+    """
+
+    def test_declare_disjuncts_in_disjunction_rule(self):
+        check_nested_disjuncts_in_flat_gdp(self, 'bigm')
+
+
+@unittest.skipUnless(gurobi_available, "Gurobi is not available")
+class IndexedDisjunction(unittest.TestCase):
+    def test_two_term_indexed_disjunction(self):
+        """
+        This test checks that we don't do anything silly with transformation Blocks in
+        the case that the Disjunction is indexed.
+        """
+        m = make_indexed_equality_model()
+
+        mbm = TransformationFactory('gdp.mbigm')
+        mbm.apply_to(m)
+
+        cons = mbm.get_transformed_constraints(m.d[1].disjuncts[0].constraint[1])
+        self.assertEqual(len(cons), 2)
+        assertExpressionsEqual(
+            self,
+            cons[0].expr,
+            m.x[1]
+            >= m.d[1].disjuncts[0].binary_indicator_var
+            + 2.0 * m.d[1].disjuncts[1].binary_indicator_var,
+        )
+        assertExpressionsEqual(
+            self,
+            cons[1].expr,
+            m.x[1]
+            <= m.d[1].disjuncts[0].binary_indicator_var
+            + 2.0 * m.d[1].disjuncts[1].binary_indicator_var,
+        )
+        cons_again = mbm.get_transformed_constraints(m.d[1].disjuncts[1].constraint[1])
+        self.assertEqual(len(cons_again), 2)
+        self.assertIs(cons_again[0], cons[0])
+        self.assertIs(cons_again[1], cons[1])
+
+        cons = mbm.get_transformed_constraints(m.d[2].disjuncts[0].constraint[1])
+        self.assertEqual(len(cons), 2)
+        assertExpressionsEqual(
+            self,
+            cons[0].expr,
+            m.x[2]
+            >= m.d[2].disjuncts[0].binary_indicator_var
+            + 2.0 * m.d[2].disjuncts[1].binary_indicator_var,
+        )
+        assertExpressionsEqual(
+            self,
+            cons[1].expr,
+            m.x[2]
+            <= m.d[2].disjuncts[0].binary_indicator_var
+            + 2.0 * m.d[2].disjuncts[1].binary_indicator_var,
+        )
+        cons_again = mbm.get_transformed_constraints(m.d[2].disjuncts[1].constraint[1])
+        self.assertEqual(len(cons_again), 2)
+        self.assertIs(cons_again[0], cons[0])
+        self.assertIs(cons_again[1], cons[1])