Merge pull request #2627 from emma58/logical-to-linear-rewrite

Adding logical_to_disjunctive, moving GDP transformations from logical-to-linear to logical-to-disjunctive

doc/OnlineDocs/modeling_extensions/gdp/concepts.rst

@@ -115,7 +115,7 @@ These logical propositions can include:
 .. |equiv| replace:: :math:`Y_1 \Leftrightarrow Y_2`
 .. |land| replace:: :math:`Y_1 \land Y_2`
 .. |lor| replace:: :math:`Y_1 \lor Y_2`
-.. |xor| replace:: :math:`Y_1 \underline{\lor} Y_2`
+.. |xor| replace:: :math:`Y_1 \veebar Y_2`
 .. |impl| replace:: :math:`Y_1 \Rightarrow Y_2`
 
 +-----------------+---------+-------------+-------------+-------------+

doc/OnlineDocs/modeling_extensions/gdp/index.rst

@@ -54,7 +54,7 @@ These can be expressed as a disjunction as follows:
     \text{constraints} \\
     \text{for }\textit{else}
     \end{gathered}\right] \\
-    Y_1 \underline{\vee} Y_2
+    Y_1 \veebar Y_2
     \end{gather*}
 
 Here, if the Boolean :math:`Y_1` is ``True``, then the constraints in the first disjunct are enforced; otherwise, the constraints in the second disjunct are enforced.

doc/OnlineDocs/modeling_extensions/gdp/modeling.rst

@@ -86,7 +86,7 @@ When the ``Disjunction`` object constructor is passed a list of lists, the outer
 
     By default, Pyomo.GDP ``Disjunction`` objects enforce an implicit "exactly one" relationship among the selection of the disjuncts (generalization of exclusive-OR).
     That is, exactly one of the ``Disjunct`` indicator variables should take a ``True`` value.
-    This can be seen as an implicit logical proposition, in our example, :math:`Y_1 \underline{\lor} Y_2`.
+    This can be seen as an implicit logical proposition, in our example, :math:`Y_1 \veebar Y_2`.
 
 Logical Propositions
 ====================
@@ -139,6 +139,11 @@ Pyomo.GDP logical expression system supported operators and their usage are list
 | Equivalence  |                        | :code:`Y[1].equivalent_to(Y[2])`  | :code:`equivalent(Y[1], Y[2])` |
 +--------------+------------------------+-----------------------------------+--------------------------------+
 
+.. note::
+
+    We omit support for most infix operators, e.g. :code:`Y[1] >> Y[2]`, due to concerns about non-intuitive Python operator precedence.
+    That is :code:`Y[1] | Y[2] >> Y[3]` would translate to :math:`Y_1 \lor (Y_2 \Rightarrow Y_3)` rather than :math:`(Y_1 \lor Y_2) \Rightarrow Y_3`
+
 In addition, the following constraint-programming-inspired operators are provided: ``exactly``, ``atmost``, and ``atleast``.
 These predicates enforce, respectively, that exactly, at most, or at least N of their ``BooleanVar`` arguments are ``True``.
 
@@ -148,26 +153,23 @@ Usage:
 - :code:`atmost(3, Y)`
 - :code:`exactly(3, Y)`
 
-.. note::
-
-    We omit support for most infix operators, e.g. :code:`Y[1] >> Y[2]`, due to concerns about non-intuitive Python operator precedence.
-    That is :code:`Y[1] | Y[2] >> Y[3]` would translate to :math:`Y_1 \lor (Y_2 \Rightarrow Y_3)` rather than :math:`(Y_1 \lor Y_2) \Rightarrow Y_3`
-
 .. doctest::
 
     >>> m = ConcreteModel()
     >>> m.my_set = RangeSet(4)
     >>> m.Y = BooleanVar(m.my_set)
     >>> m.p = LogicalConstraint(expr=atleast(3, m.Y))
-    >>> TransformationFactory('core.logical_to_linear').apply_to(m)
-    >>> m.logic_to_linear.transformed_constraints.pprint()  # constraint auto-generated by transformation
-    transformed_constraints : Size=1, Index=logic_to_linear.transformed_constraints_index, Active=True
-        Key : Lower : Body                                                          : Upper : Active
-          1 :   3.0 : Y_asbinary[1] + Y_asbinary[2] + Y_asbinary[3] + Y_asbinary[4] :  +Inf :   True
     >>> m.p.pprint()
     p : Size=1, Index=None, Active=False
         Key  : Body                                 : Active
         None : atleast(3: [Y[1], Y[2], Y[3], Y[4]]) :  False
+    >>> TransformationFactory('core.logical_to_linear').apply_to(m)
+        ...
+    >>> # constraint auto-generated by transformation
+    >>> m.logic_to_linear.transformed_constraints.pprint()
+    transformed_constraints : Size=1, Index=logic_to_linear.transformed_constraints_index, Active=True
+        Key : Lower : Body                                                          : Upper : Active
+          1 :   3.0 : Y_asbinary[1] + Y_asbinary[2] + Y_asbinary[3] + Y_asbinary[4] :  +Inf :   True
 
 We elaborate on the ``logical_to_linear`` transformation :ref:`on the next page <gdp-reformulations>`.
 
@@ -223,8 +225,8 @@ Here, we demonstrate this capability with a toy example:
     \min~&x\\
     \text{s.t.}~&\left[\begin{gathered}Y_1\\x \geq 2\end{gathered}\right] \vee \left[\begin{gathered}Y_2\\x \geq 3\end{gathered}\right]\\
     &\left[\begin{gathered}Y_3\\x \leq 8\end{gathered}\right] \vee \left[\begin{gathered}Y_4\\x = 2.5\end{gathered}\right] \\
-    &Y_1 \underline{\vee} Y_2\\
-    &Y_3 \underline{\vee} Y_4\\
+    &Y_1 \veebar Y_2\\
+    &Y_3 \veebar Y_4\\
     &Y_1 \Rightarrow Y_4
 
 .. doctest::
@@ -359,7 +361,7 @@ In the ``logical_to_linear`` transformation, we automatically convert these spec
 Additional Examples
 ===================
 
-The following models all work and are equivalent for :math:`\left[x = 0\right] \underline{\lor} \left[y = 0\right]`:
+The following models all work and are equivalent for :math:`\left[x = 0\right] \veebar \left[y = 0\right]`:
 
 .. doctest::
 

doc/OnlineDocs/modeling_extensions/gdp/solving.rst

@@ -29,15 +29,25 @@ Logical constraints
 
 .. note::
 
-    Historically it was required to convert logical propositions to
-    algebraic form prior to use of the MI(N)LP reformulations and the
-    GDPopt solver. However, this is mathematically incorrect since these
-    reformulations convert logical formulations to algebraic formulations.
-    It is therefore recommended to use both the MI(N)LP reformulations
-    and GDPopt directly to transform or solve GDPs that include logical
-    propositions.
+    Historically users needed to explicitly convert logical propositions
+    to algebraic form prior to invoking the GDP MI(N)LP reformulations
+    or the GDPopt solver. However, this is mathematically incorrect
+    since the GDP MI(N)LP reformulations themselves convert logical
+    formulations to algebraic formulations.  The current recommended
+    practice is to pass the entire (mixed logical / algebraic) model to
+    the MI(N)LP reformulations or GDPopt directly.
+
+There are several approaches to convert logical constraints into
+algebraic form.
+
+Conjunctive Normal Form
+^^^^^^^^^^^^^^^^^^^^^^^
 
-The following transforms logical propositions on the model to algebraic form:
+The first transformation (`core.logical_to_linear`) leverages the
+`sympy` package to generate the conjunctive normal form of the logical
+constraints and then adds the equivalent as a list algebraic
+constraints.  The following transforms logical propositions on the model
+to algebraic form:
 
 .. code::
 
@@ -61,6 +71,29 @@ Following solution of the GDP model, values of the Boolean variables may be upda
 
 .. autofunction:: pyomo.core.plugins.transform.logical_to_linear.update_boolean_vars_from_binary
 
+Factorable Programming
+^^^^^^^^^^^^^^^^^^^^^^
+
+The second transformation (`contrib.logical_to_disjunctive`) leverages
+ideas from factorable programming to first generate an equivalent set of
+"factored" logical constraints form by traversing each logical
+proposition and replacing each logical operator with an additional
+Boolean variable and then adding the "simple" logical constraint that
+equates the new Boolean variable with the single logical operator.
+
+The resulting "simple" logical constraints are converted to either MIP
+or GDP form: if the constraint contains only Boolean variables, then
+then MIP representation is emitted.  Logical constraints with mixed
+integer-Boolean arguments (e.g., `atmost`, `atleast`, `exactly`, etc.)
+are converted to a disjunctive representation.
+
+As this transformation both avoids the conversion into `sympy` and only
+requires a single traversal of each logical constraint,
+`contrib.logical_to_disjunctive` is significantly faster than
+`core.logical_to_linear` at the cost of a larger model.  In practice,
+the cost of the larger model is negated by the effectiveness of the MIP
+presolve in most solvers.
+
 Reformulation to MI(N)LP
 ------------------------
 

pyomo/common/errors.py

@@ -70,3 +70,23 @@ class TempfileContextError(PyomoException, IndexError):
 
     """
     pass
+
+
+class MouseTrap(PyomoException, NotImplementedError):
+    """
+    Exception class used to throw errors for not-implemented functionality
+    that might be rational to support (i.e., we already gave you a cookie)
+    but risks taking Pyomo's flexibility a step beyond what is sane,
+    or solvable, or communicable to a solver, etc. (i.e., Really? Now you
+    want a glass of milk too?)
+    """
+    def __init__(self, val):
+        self.parameter = val
+
+    def __str__(self):
+        return ("Sorry, mouse, no cookies here!\n\t%s\n"
+                "\tThis is functionality we think may be rational to "
+                "support, but is not yet implemented since it might go "
+                "beyond what can practically be solved. However, please "
+                "feed the mice: pull requests are always welcome!"
+                % (repr(self.parameter), ))

pyomo/contrib/cp/__init__.py

@@ -5,3 +5,5 @@
     DocplexWriter, CPOptimizerSolver)
 from pyomo.contrib.cp.scheduling_expr.step_function_expressions import (
     AlwaysIn, Step, Pulse)
+# register logical_to_disjunctive transformation
+import pyomo.contrib.cp.transform.logical_to_disjunctive_program

pyomo/contrib/cp/plugins.py

@@ -0,0 +1,15 @@
+#  ___________________________________________________________________________
+#
+#  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.
+#  ___________________________________________________________________________
+
+def load():
+    from . import interval_var
+    from .repn import docplex_writer
+    from .transform import logical_to_disjunctive_program

pyomo/contrib/cp/repn/docplex_writer.py

@@ -10,8 +10,6 @@
 #  ___________________________________________________________________________
 
 from pyomo.common.dependencies import attempt_import
-cp, docplex_available = attempt_import('docplex.cp.model')
-cp_solver, docplex_available = attempt_import('docplex.cp.solver')
 
 import itertools
 import logging
@@ -68,6 +66,33 @@
 from pyomo.network import Port
 ###
 
+def _finalize_docplex(module, available):
+    if not available:
+        return
+    _deferred_element_getattr_dispatcher['start_time'] = module.start_of
+    _deferred_element_getattr_dispatcher['end_time'] = module.end_of
+    _deferred_element_getattr_dispatcher['length'] = module.length_of
+    _deferred_element_getattr_dispatcher['is_present'] = module.presence_of
+
+    # Scheduling dispatchers
+    _before_dispatchers[_START_TIME, _START_TIME] = module.start_before_start
+    _before_dispatchers[_START_TIME, _END_TIME] = module.start_before_end
+    _before_dispatchers[_END_TIME, _START_TIME] = module.end_before_start
+    _before_dispatchers[_END_TIME, _END_TIME] = module.end_before_end
+
+    _at_dispatchers[_START_TIME, _START_TIME] = module.start_at_start
+    _at_dispatchers[_START_TIME, _END_TIME] = module.start_at_end
+    _at_dispatchers[_END_TIME, _START_TIME] = module.end_at_start
+    _at_dispatchers[_END_TIME, _END_TIME] = module.end_at_end
+
+    _time_point_dispatchers[_START_TIME] = module.start_of
+    _time_point_dispatchers[_END_TIME] = module.end_of
+
+
+cp, docplex_available = attempt_import(
+    'docplex.cp.model', callback=_finalize_docplex)
+cp_solver, docplex_available = attempt_import('docplex.cp.solver')
+
 logger = logging.getLogger('pyomo.contrib.cp')
 
 # These are things that don't need special handling:
@@ -217,13 +242,8 @@ def _handle_getitem(visitor, node, *data):
     'xor': _XOR,
     'equivalent_to': _EQUIVALENT_TO,
 }
-if docplex_available:
-    _deferred_element_getattr_dispatcher = {
-        'start_time': cp.start_of,
-        'end_time': cp.end_of,
-        'length': cp.length_of,
-        'is_present': cp.presence_of,
-    }
+# This will get populated when cp is finally imported
+_deferred_element_getattr_dispatcher = {}
 
 def _handle_getattr(visitor, node, obj, attr):
     # We either end up here because we do not yet know the list of variables to
@@ -676,25 +696,13 @@ def _handle_call(visitor, node, *args):
 # Scheduling
 ##
 
-if docplex_available:
-    _before_dispatchers = {
-        (_START_TIME, _START_TIME): cp.start_before_start,
-        (_START_TIME, _END_TIME): cp.start_before_end,
-        (_END_TIME, _START_TIME): cp.end_before_start,
-        (_END_TIME, _END_TIME): cp.end_before_end,
-    }
-    _at_dispatchers = {
-        (_START_TIME, _START_TIME): cp.start_at_start,
-        (_START_TIME, _END_TIME): cp.start_at_end,
-        (_END_TIME, _START_TIME): cp.end_at_start,
-        (_END_TIME, _END_TIME): cp.end_at_end
-    }
-    _time_point_dispatchers = {
-        _START_TIME: cp.start_of,
-        _END_TIME: cp.end_of,
-        _GENERAL: lambda x: x,
-        _ELEMENT_CONSTRAINT: lambda x: x,
-    }
+# This will get populated when cp is finally imported
+_before_dispatchers = {}
+_at_dispatchers = {}
+_time_point_dispatchers = {
+    _GENERAL: lambda x: x,
+    _ELEMENT_CONSTRAINT: lambda x: x,
+}
 
 _non_precedence_types = {_GENERAL, _ELEMENT_CONSTRAINT}