Merge pull request #3145 from emma58/fix-nested-everywhere

pyomo.gdp: handle nested GDPs correctly in all the transformations

pyomo/gdp/plugins/bigm.py

@@ -224,11 +224,10 @@ def _transform_disjunctionData(
             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:
-            xorConstraint[index] = or_expr == rhs
+            xorConstraint[index] = or_expr == 1
         else:
-            xorConstraint[index] = or_expr >= rhs
+            xorConstraint[index] = or_expr >= 1
         # Mark the DisjunctionData as transformed by mapping it to its XOR
         # constraint.
         obj._algebraic_constraint = weakref_ref(xorConstraint[index])

pyomo/gdp/plugins/binary_multiplication.py

@@ -92,11 +92,10 @@ def _transform_disjunctionData(
             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
+            xorConstraint[index] = or_expr == 1
         else:
-            xorConstraint[index] = or_expr >= rhs
+            xorConstraint[index] = or_expr >= 1
         # Mark the DisjunctionData as transformed by mapping it to its XOR
         # constraint.
         obj._algebraic_constraint = weakref_ref(xorConstraint[index])

pyomo/gdp/plugins/gdp_to_mip_transformation.py

@@ -213,21 +213,26 @@ def _setup_transform_disjunctionData(self, obj, root_disjunct):
                 "likely indicative of a modeling error." % obj.name
             )
 
-        # Create or fetch the transformation block
+        # We always need to create or fetch a transformation block on the parent block.
+        trans_block, new_block = self._add_transformation_block(obj.parent_block())
+        # This is where we put exactly_one/or constraint
+        algebraic_constraint = self._add_xor_constraint(
+            obj.parent_component(), trans_block
+        )
+
+        # If requested, create or fetch the transformation block above the
+        # nested hierarchy
         if root_disjunct is not None:
-            # We want to put all the transformed things on the root
-            # Disjunct's parent's block so that they do not get
-            # re-transformed
-            transBlock, new_block = self._add_transformation_block(
+            # We want to put some transformed things on the root Disjunct's
+            # parent's block so that they do not get re-transformed. (Note this
+            # is never true for hull, but it calls this method with
+            # root_disjunct=None. BigM can't put the exactly-one constraint up
+            # here, but it can put everything else.)
+            trans_block, new_block = self._add_transformation_block(
                 root_disjunct.parent_block()
             )
-        else:
-            # This isn't nested--just put it on the parent block.
-            transBlock, new_block = self._add_transformation_block(obj.parent_block())
-
-        xorConstraint = self._add_xor_constraint(obj.parent_component(), transBlock)
 
-        return transBlock, xorConstraint
+        return trans_block, algebraic_constraint
 
     def _get_disjunct_transformation_block(self, disjunct, transBlock):
         if disjunct.transformation_block is not None:

pyomo/gdp/plugins/multiple_bigm.py

@@ -336,8 +336,7 @@ def _transform_disjunctionData(self, obj, index, parent_disjunct, root_disjunct)
         for disjunct in active_disjuncts:
             or_expr += disjunct.indicator_var.get_associated_binary()
             self._transform_disjunct(disjunct, transBlock, active_disjuncts, Ms)
-        rhs = 1 if parent_disjunct is None else parent_disjunct.binary_indicator_var
-        algebraic_constraint.add(index, (or_expr, rhs))
+        algebraic_constraint.add(index, or_expr == 1)
         # map the DisjunctionData to its XOR constraint to mark it as
         # transformed
         obj._algebraic_constraint = weakref_ref(algebraic_constraint[index])

pyomo/gdp/tests/common_tests.py

@@ -1944,3 +1944,17 @@ def check_nested_disjuncts_in_flat_gdp(self, transformation):
     for t in m.T:
         self.assertTrue(value(m.disj1[t].indicator_var))
         self.assertTrue(value(m.disj1[t].sub1.indicator_var))
+
+
+def check_do_not_assume_nested_indicators_local(self, transformation):
+    m = models.why_indicator_vars_are_not_always_local()
+    TransformationFactory(transformation).apply_to(m)
+
+    results = SolverFactory('gurobi').solve(m)
+    self.assertEqual(results.solver.termination_condition, TerminationCondition.optimal)
+    self.assertAlmostEqual(value(m.obj), 9)
+    self.assertAlmostEqual(value(m.x), 9)
+    self.assertTrue(value(m.Y2.indicator_var))
+    self.assertFalse(value(m.Y1.indicator_var))
+    self.assertTrue(value(m.Z1.indicator_var))
+    self.assertTrue(value(m.Z1.indicator_var))

pyomo/gdp/tests/models.py

@@ -563,6 +563,44 @@ def makeNestedDisjunctions_NestedDisjuncts():
     return m
 
 
+def why_indicator_vars_are_not_always_local():
+    m = ConcreteModel()
+    m.x = Var(bounds=(1, 10))
+
+    @m.Disjunct()
+    def Z1(d):
+        m = d.model()
+        d.c = Constraint(expr=m.x >= 1.1)
+
+    @m.Disjunct()
+    def Z2(d):
+        m = d.model()
+        d.c = Constraint(expr=m.x >= 1.2)
+
+    @m.Disjunct()
+    def Y1(d):
+        m = d.model()
+        d.c = Constraint(expr=(1.15, m.x, 8))
+        d.disjunction = Disjunction(expr=[m.Z1, m.Z2])
+
+    @m.Disjunct()
+    def Y2(d):
+        m = d.model()
+        d.c = Constraint(expr=m.x == 9)
+
+    m.disjunction = Disjunction(expr=[m.Y1, m.Y2])
+
+    m.logical_cons = LogicalConstraint(
+        expr=m.Y2.indicator_var.implies(m.Z1.indicator_var.land(m.Z2.indicator_var))
+    )
+
+    # optimal value is 9, but it will be 8 if we wrongly assume that the nested
+    # indicator_vars are local.
+    m.obj = Objective(expr=m.x, sense=maximize)
+
+    return m
+
+
 def makeTwoSimpleDisjunctions():
     """Two SimpleDisjunctions on the same model."""
     m = ConcreteModel()

pyomo/gdp/tests/test_bigm.py

@@ -33,6 +33,7 @@
     assertExpressionsStructurallyEqual,
 )
 from pyomo.repn import generate_standard_repn
+from pyomo.repn.linear import LinearRepnVisitor
 from pyomo.common.log import LoggingIntercept
 import logging
 
@@ -1764,22 +1765,19 @@ def test_transformation_block_structure(self):
         # we have the XOR constraints for both the outer and inner disjunctions
         self.assertIsInstance(transBlock.component("disjunction_xor"), Constraint)
 
-    def test_transformation_block_on_inner_disjunct_empty(self):
-        m = models.makeNestedDisjunctions()
-        TransformationFactory('gdp.bigm').apply_to(m)
-        self.assertIsNone(m.disjunct[1].component("_pyomo_gdp_bigm_reformulation"))
-
     def test_mappings_between_disjunctions_and_xors(self):
         m = models.makeNestedDisjunctions()
         transform = TransformationFactory('gdp.bigm')
         transform.apply_to(m)
 
         transBlock1 = m.component("_pyomo_gdp_bigm_reformulation")
+        transBlock2 = m.disjunct[1].component("_pyomo_gdp_bigm_reformulation")
+        transBlock3 = m.simpledisjunct.component("_pyomo_gdp_bigm_reformulation")
 
         disjunctionPairs = [
             (m.disjunction, transBlock1.disjunction_xor),
-            (m.disjunct[1].innerdisjunction[0], transBlock1.innerdisjunction_xor_4[0]),
-            (m.simpledisjunct.innerdisjunction, transBlock1.innerdisjunction_xor),
+            (m.disjunct[1].innerdisjunction[0], transBlock2.innerdisjunction_xor[0]),
+            (m.simpledisjunct.innerdisjunction, transBlock3.innerdisjunction_xor),
         ]
 
         # check disjunction mappings
@@ -1899,19 +1897,32 @@ def check_bigM_constraint(self, cons, variable, M, indicator_var):
         ct.check_linear_coef(self, repn, variable, 1)
         ct.check_linear_coef(self, repn, indicator_var, M)
 
-    def check_inner_xor_constraint(
-        self, inner_disjunction, outer_disjunct, inner_disjuncts
-    ):
-        self.assertIsNotNone(inner_disjunction.algebraic_constraint)
-        cons = inner_disjunction.algebraic_constraint
-        self.assertEqual(cons.lower, 0)
-        self.assertEqual(cons.upper, 0)
-        repn = generate_standard_repn(cons.body)
-        self.assertTrue(repn.is_linear())
-        self.assertEqual(repn.constant, 0)
-        for disj in inner_disjuncts:
-            ct.check_linear_coef(self, repn, disj.binary_indicator_var, 1)
-        ct.check_linear_coef(self, repn, outer_disjunct.binary_indicator_var, -1)
+    def check_inner_xor_constraint(self, inner_disjunction, outer_disjunct, bigm):
+        inner_xor = inner_disjunction.algebraic_constraint
+        sum_indicators = sum(
+            d.binary_indicator_var for d in inner_disjunction.disjuncts
+        )
+        assertExpressionsEqual(self, inner_xor.expr, sum_indicators == 1)
+        # this guy has been transformed
+        self.assertFalse(inner_xor.active)
+        cons = bigm.get_transformed_constraints(inner_xor)
+        self.assertEqual(len(cons), 2)
+        lb = cons[0]
+        ct.check_obj_in_active_tree(self, lb)
+        lb_expr = self.simplify_cons(lb, leq=False)
+        assertExpressionsEqual(
+            self,
+            lb_expr,
+            1.0 <= sum_indicators - outer_disjunct.binary_indicator_var + 1.0,
+        )
+        ub = cons[1]
+        ct.check_obj_in_active_tree(self, ub)
+        ub_expr = self.simplify_cons(ub, leq=True)
+        assertExpressionsEqual(
+            self,
+            ub_expr,
+            sum_indicators + outer_disjunct.binary_indicator_var - 1 <= 1.0,
+        )
 
     def test_transformed_constraints(self):
         # We'll check all the transformed constraints to make sure
@@ -1993,26 +2004,9 @@ def test_transformed_constraints(self):
 
         # Here we check that the xor constraint from
         # simpledisjunct.innerdisjunction is transformed.
-        cons5 = m.simpledisjunct.innerdisjunction.algebraic_constraint
-        self.assertIsNotNone(cons5)
         self.check_inner_xor_constraint(
-            m.simpledisjunct.innerdisjunction,
-            m.simpledisjunct,
-            [m.simpledisjunct.innerdisjunct0, m.simpledisjunct.innerdisjunct1],
-        )
-        self.assertIsInstance(cons5, Constraint)
-        self.assertEqual(cons5.lower, 0)
-        self.assertEqual(cons5.upper, 0)
-        repn = generate_standard_repn(cons5.body)
-        self.assertTrue(repn.is_linear())
-        self.assertEqual(repn.constant, 0)
-        ct.check_linear_coef(
-            self, repn, m.simpledisjunct.innerdisjunct0.binary_indicator_var, 1
-        )
-        ct.check_linear_coef(
-            self, repn, m.simpledisjunct.innerdisjunct1.binary_indicator_var, 1
+            m.simpledisjunct.innerdisjunction, m.simpledisjunct, bigm
         )
-        ct.check_linear_coef(self, repn, m.simpledisjunct.binary_indicator_var, -1)
 
         cons6 = bigm.get_transformed_constraints(m.disjunct[0].c)
         self.assertEqual(len(cons6), 2)
@@ -2028,9 +2022,7 @@ def test_transformed_constraints(self):
         # now we check that the xor constraint from disjunct[1].innerdisjunction
         # is correct.
         self.check_inner_xor_constraint(
-            m.disjunct[1].innerdisjunction[0],
-            m.disjunct[1],
-            [m.disjunct[1].innerdisjunct[0], m.disjunct[1].innerdisjunct[1]],
+            m.disjunct[1].innerdisjunction[0], m.disjunct[1], bigm
         )
 
         cons8 = bigm.get_transformed_constraints(m.disjunct[1].c)
@@ -2136,33 +2128,27 @@ def check_second_disjunct_constraint(self, disj2c, x, ind_var):
             ct.check_squared_term_coef(self, repn, x[i], 1)
             ct.check_linear_coef(self, repn, x[i], -6)
 
+    def simplify_cons(self, cons, leq):
+        visitor = LinearRepnVisitor({}, {}, {}, None)
+        repn = visitor.walk_expression(cons.body)
+        self.assertIsNone(repn.nonlinear)
+        if leq:
+            self.assertIsNone(cons.lower)
+            ub = cons.upper
+            return ub >= repn.to_expression(visitor)
+        else:
+            self.assertIsNone(cons.upper)
+            lb = cons.lower
+            return lb <= repn.to_expression(visitor)
+
     def check_hierarchical_nested_model(self, m, bigm):
         outer_xor = m.disjunction_block.disjunction.algebraic_constraint
         ct.check_two_term_disjunction_xor(
             self, outer_xor, m.disj1, m.disjunct_block.disj2
         )
 
-        inner_xor = m.disjunct_block.disj2.disjunction.algebraic_constraint
-        self.assertEqual(inner_xor.lower, 0)
-        self.assertEqual(inner_xor.upper, 0)
-        repn = generate_standard_repn(inner_xor.body)
-        self.assertTrue(repn.is_linear())
-        self.assertEqual(len(repn.linear_vars), 3)
-        self.assertEqual(repn.constant, 0)
-        ct.check_linear_coef(
-            self,
-            repn,
-            m.disjunct_block.disj2.disjunction_disjuncts[0].binary_indicator_var,
-            1,
-        )
-        ct.check_linear_coef(
-            self,
-            repn,
-            m.disjunct_block.disj2.disjunction_disjuncts[1].binary_indicator_var,
-            1,
-        )
-        ct.check_linear_coef(
-            self, repn, m.disjunct_block.disj2.binary_indicator_var, -1
+        self.check_inner_xor_constraint(
+            m.disjunct_block.disj2.disjunction, m.disjunct_block.disj2, bigm
         )
 
         # outer disjunction constraints
@@ -2214,6 +2200,10 @@ def test_decl_order_opposite_instantiation_order(self):
         # the same check to make sure everything is transformed correctly.
         self.check_hierarchical_nested_model(m, bigm)
 
+    @unittest.skipUnless(gurobi_available, "Gurobi is not available")
+    def test_do_not_assume_nested_indicators_local(self):
+        ct.check_do_not_assume_nested_indicators_local(self, 'gdp.bigm')
+
 
 class IndexedDisjunction(unittest.TestCase):
     # this tests that if the targets are a subset of the

pyomo/gdp/tests/test_binary_multiplication.py

@@ -18,6 +18,7 @@
     ConcreteModel,
     Var,
     Any,
+    SolverFactory,
 )
 from pyomo.gdp import Disjunct, Disjunction
 from pyomo.core.expr.compare import assertExpressionsEqual
@@ -30,6 +31,11 @@
 
 import random
 
+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):
@@ -297,5 +303,13 @@ def test_local_var(self):
         self.assertEqual(eq.ub, 0)
 
 
+class TestNestedGDP(unittest.TestCase):
+    @unittest.skipUnless(gurobi_available, "Gurobi is not available")
+    def test_do_not_assume_nested_indicators_local(self):
+        ct.check_do_not_assume_nested_indicators_local(
+            self, 'gdp.binary_multiplication'
+        )
+
+
 if __name__ == '__main__':
     unittest.main()

pyomo/gdp/tests/test_hull.py

@@ -2030,6 +2030,10 @@ def test_nested_with_var_that_skips_a_level(self):
         cons_expr = self.simplify_cons(cons)
         assertExpressionsEqual(self, cons_expr, m.y - y_y2 - y_y1 == 0.0)
 
+    @unittest.skipUnless(gurobi_available, "Gurobi is not available")
+    def test_do_not_assume_nested_indicators_local(self):
+        ct.check_do_not_assume_nested_indicators_local(self, 'gdp.hull')
+
 
 class TestSpecialCases(unittest.TestCase):
     def test_local_vars(self):