.. image:: /../logos/gdp/Pyomo-GDP-150.png
    :scale: 20%
    :class: no-scaled-link
    :align: right

*********************
Modeling in Pyomo.GDP
*********************

.. testsetup::

    from pyomo.environ import (
        ConcreteModel, RangeSet, BooleanVar, LogicalConstraint,
        TransformationFactory, atleast, SolverFactory, Objective,
        Constraint, Var, land, Reference
    )
    from pyomo.gdp import Disjunct, Disjunction
    from pyomo.core.plugins.transform.logical_to_linear import update_boolean_vars_from_binary


Disjunctions
============

To demonstrate modeling with disjunctions in Pyomo.GDP, we revisit the small example from :ref:`the previous page <gdp-disjunctions-concept>`.

.. math::

    \left[\begin{gathered}
    Y_1 \\
    \exp(x_2) - 1 = x_1 \\
    x_3 = x_4 = 0
    \end{gathered}
    \right] \bigvee \left[\begin{gathered}
    Y_2 \\
    \exp\left(\frac{x_4}{1.2}\right) - 1 = x_3 \\
    x_1 = x_2 = 0
    \end{gathered}
    \right]

Explicit syntax: more descriptive
---------------------------------

Pyomo.GDP explicit syntax (see below) provides more clarity in the declaration of each modeling object,  and gives the user explicit control over the ``Disjunct`` names.
Assuming the ``ConcreteModel`` object :code:`m` and variables have been defined, lines 1 and 5 declare the ``Disjunct`` objects corresponding to selection of unit 1 and 2, respectively.
Lines 2 and 6 define the input-output relations for each unit, and lines 3-4 and 7-8 enforce zero flow through the unit that is not selected.
Finally, line 9 declares the logical disjunction between the two disjunctive terms.

.. code-block:: python
    :linenos:

    m.unit1 = Disjunct()
    m.unit1.inout = Constraint(expr=exp(m.x[2]) - 1 == m.x[1])
    m.unit1.no_unit2_flow1 = Constraint(expr=m.x[3] == 0)
    m.unit1.no_unit2_flow2 = Constraint(expr=m.x[4] == 0)
    m.unit2 = Disjunct()
    m.unit2.inout = Constraint(expr=exp(m.x[4] / 1.2) - 1 == m.x[3])
    m.unit2.no_unit1_flow1 = Constraint(expr=m.x[1] == 0)
    m.unit2.no_unit1_flow2 = Constraint(expr=m.x[2] == 0)
    m.use_unit1or2 = Disjunction(expr=[m.unit1, m.unit2])

The indicator variables for each disjunct :math:`Y_1` and :math:`Y_2` are automatically generated by Pyomo.GDP, accessible via :code:`m.unit1.indicator_var` and :code:`m.unit2.indicator_var`.

Compact syntax: more concise
----------------------------

For more advanced users, a compact syntax is also available below, taking advantage of the ability to declare disjuncts and constraints implicitly.
When the ``Disjunction`` object constructor is passed a list of lists, the outer list defines the disjuncts and the inner list defines the constraint expressions associated with the respective disjunct.

.. code-block:: python
    :linenos:

    m.use1or2 = Disjunction(expr=[
        # First disjunct
        [exp(m.x[2])-1 == m.x[1],
         m.x[3] == 0, m.x[4] == 0],
        # Second disjunct
        [exp(m.x[4]/1.2)-1 == m.x[3],
         m.x[1] == 0, m.x[2] == 0]])

.. note::

    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 \veebar Y_2`.

Logical Propositions
====================

Pyomo.GDP also supports the use of logical propositions through the use of the ``BooleanVar`` and ``LogicalConstraint`` objects.
The ``BooleanVar`` object in Pyomo represents Boolean variables, analogous to ``Var`` for numeric variables.
``BooleanVar`` can be indexed over a Pyomo ``Set``, as below:

.. doctest::

    >>> m = ConcreteModel()
    >>> m.my_set = RangeSet(4)
    >>> m.Y = BooleanVar(m.my_set)
    >>> m.Y.display()
    Y : Size=4, Index=my_set
        Key : Value : Fixed : Stale
          1 :  None : False :  True
          2 :  None : False :  True
          3 :  None : False :  True
          4 :  None : False :  True

Using these Boolean variables, we can define ``LogicalConstraint`` objects, analogous to algebraic ``Constraint`` objects.

.. doctest::

    >>> m.p = LogicalConstraint(expr=m.Y[1].implies(m.Y[2] & m.Y[3]) | m.Y[4])
    >>> m.p.pprint()
    p : Size=1, Index=None, Active=True
        Key  : Body                          : Active
        None : (Y[1] --> Y[2] ∧ Y[3]) ∨ Y[4] :   True

Supported Logical Operators
---------------------------

Pyomo.GDP logical expression system supported operators and their usage are listed in the table below.

+--------------+------------------------+-----------------------------------+--------------------------------+
| Operator     | Operator               | Method                            | Function                       |
+==============+========================+===================================+================================+
| Negation     | :code:`~Y[1]`          |                                   | :code:`lnot(Y[1])`             |
+--------------+------------------------+-----------------------------------+--------------------------------+
| Conjunction  | :code:`Y[1] & Y[2]`    | :code:`Y[1].land(Y[2])`           | :code:`land(Y[1],Y[2])`        |
+--------------+------------------------+-----------------------------------+--------------------------------+
| Disjunction  | :code:`Y[1] | Y[2]`    | :code:`Y[1].lor(Y[2])`            | :code:`lor(Y[1],Y[2])`         |
+--------------+------------------------+-----------------------------------+--------------------------------+
| Exclusive OR | :code:`Y[1] ^ Y[2]`    | :code:`Y[1].xor(Y[2])`            | :code:`xor(Y[1], Y[2])`        |
+--------------+------------------------+-----------------------------------+--------------------------------+
| Implication  |                        | :code:`Y[1].implies(Y[2])`        | :code:`implies(Y[1], Y[2])`    |
+--------------+------------------------+-----------------------------------+--------------------------------+
| Equivalence  |                        | :code:`Y[1].equivalent_to(Y[2])`  | :code:`equivalent(Y[1], Y[2])` |
+--------------+------------------------+-----------------------------------+--------------------------------+

.. note::

    We omit support for some 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``.

Usage:

- :code:`atleast(3, Y[1], Y[2], Y[3])`
- :code:`atmost(3, Y)`
- :code:`exactly(3, Y)`

.. doctest::

    >>> m = ConcreteModel()
    >>> m.my_set = RangeSet(4)
    >>> m.Y = BooleanVar(m.my_set)
    >>> m.p = LogicalConstraint(expr=atleast(3, m.Y))
    >>> m.p.pprint()
    p : Size=1, Index=None, Active=True
        Key  : Body                                 : Active
        None : atleast(3: [Y[1], Y[2], Y[3], Y[4]]) :   True
    >>> 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={1}, 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>`.

Indexed logical constraints
---------------------------

Like ``Constraint`` objects for algebraic expressions, ``LogicalConstraint`` objects can be indexed.
An example of this usage may be found below for the expression:

.. math::

    Y_{i+1} \Rightarrow Y_{i}, \quad i \in \{1, 2, \dots, n-1\}

.. doctest::

    >>> m = ConcreteModel()
    >>> n = 5
    >>> m.I = RangeSet(n)
    >>> m.Y = BooleanVar(m.I)

    >>> @m.LogicalConstraint(m.I)
    ... def p(m, i):
    ...     return m.Y[i+1].implies(m.Y[i]) if i < n else Constraint.Skip

    >>> m.p.pprint()
    p : Size=4, Index=I, Active=True
        Key : Body          : Active
          1 : Y[2] --> Y[1] :   True
          2 : Y[3] --> Y[2] :   True
          3 : Y[4] --> Y[3] :   True
          4 : Y[5] --> Y[4] :   True

Integration with Disjunctions
-----------------------------

.. note::

    Historically, the ``indicator_var`` on ``Disjunct`` objects was
    implemented as a binary ``Var``.  Beginning in Pyomo 6.0, that has
    been changed to the more mathematically correct ``BooleanVar``, with
    the associated binary variable available as
    ``binary_indicator_var``.

The logical expression system is designed to augment the previously
introduced ``Disjunct`` and ``Disjunction`` components.  Mathematically,
the disjunct indicator variable is Boolean, and can be used directly in
logical propositions.

Here, we demonstrate this capability with a toy example:

.. math::
   :nowrap:

   \[\begin{array}{ll}
    \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 \veebar Y_2 \\
    & Y_3 \veebar Y_4 \\
    & Y_1 \Rightarrow Y_4
   \end{array}\]

.. doctest::
    :skipif: not glpk_available

    >>> m = ConcreteModel()
    >>> m.s = RangeSet(4)
    >>> m.ds = RangeSet(2)
    >>> m.d = Disjunct(m.s)
    >>> m.djn = Disjunction(m.ds)
    >>> m.djn[1] = [m.d[1], m.d[2]]
    >>> m.djn[2] = [m.d[3], m.d[4]]
    >>> m.x = Var(bounds=(-2, 10))
    >>> m.d[1].c = Constraint(expr=m.x >= 2)
    >>> m.d[2].c = Constraint(expr=m.x >= 3)
    >>> m.d[3].c = Constraint(expr=m.x <= 8)
    >>> m.d[4].c = Constraint(expr=m.x == 2.5)
    >>> m.o = Objective(expr=m.x)

    >>> # Add the logical proposition
    >>> m.p = LogicalConstraint(
    ...    expr=m.d[1].indicator_var.implies(m.d[4].indicator_var))
    >>> # Note: the implicit XOR enforced by m.djn[1] and m.djn[2] still apply

    >>> # Apply the Big-M reformulation: It will convert the logical
    >>> # propositions to algebraic expressions.
    >>> TransformationFactory('gdp.bigm').apply_to(m)

    >>> # Before solve, Boolean vars have no value
    >>> Reference(m.d[:].indicator_var).display()
    IndexedBooleanVar : Size=4, Index=s, ReferenceTo=d[:].indicator_var
        Key : Value : Fixed : Stale
          1 :  None : False :  True
          2 :  None : False :  True
          3 :  None : False :  True
          4 :  None : False :  True

    >>> # Solve the reformulated model
    >>> run_data = SolverFactory('glpk').solve(m)
    >>> Reference(m.d[:].indicator_var).display()
    IndexedBooleanVar : Size=4, Index=s, ReferenceTo=d[:].indicator_var
        Key : Value : Fixed : Stale
          1 :  True : False : False
          2 : False : False : False
          3 : False : False : False
          4 :  True : False : False

.. _gdp-advanced-examples:

Advanced LogicalConstraint Examples
===================================

Support for complex nested expressions is a key benefit of the logical expression system.
Below are examples of expressions that we support, and with some, an explanation of their implementation.

Composition of standard operators
---------------------------------

.. math::
    Y_1 \vee Y_2 \implies Y_3 \wedge \neg Y_4 \wedge (Y_5 \vee Y_6)

.. code::

    m.p = LogicalConstraint(expr=(m.Y[1] | m.Y[2]).implies(
        m.Y[3] & ~m.Y[4] & (m.Y[5] | m.Y[6]))
    )

Expressions within CP-type operators
------------------------------------

.. math::
    \text{atleast}(3, Y_1, Y_2 \vee Y_3, Y_4 \Rightarrow Y_5, Y_6)

Here, augmented variables may be automatically added to the model as follows:

.. math::
    \text{atleast}(3, Y_1, Y_A, Y_B, Y_6)\\
    Y_A \Leftrightarrow Y_2 \vee Y_3\\
    Y_B \Leftrightarrow (Y_4 \Rightarrow Y_5)

.. code::

    m.p = LogicalConstraint(
        expr=atleast(3, m.Y[1], Or(m.Y[2], m.Y[3]), m.Y[4].implies(m.Y[5]), m.Y[6]))

Nested CP-style operators
-------------------------

.. math::
    \text{atleast}(2, Y_1, \text{exactly}(2, Y_2, Y_3, Y_4), Y_5, Y_6)

Here, we again need to add augmented variables:

.. math::
    \text{atleast}(2, Y_1, Y_A, Y_5, Y_6)\\
    Y_A \Leftrightarrow \text{exactly}(2, Y_2, Y_3, Y_4)

However, we also need to further interpret the second statement as a disjunction:

.. math::
    :nowrap:

    \begin{gather*}
    \text{atleast}(2, Y_1, Y_A, Y_5, Y_6)\\
    \left[\begin{gathered}Y_A\\\text{exactly}(2, Y_2, Y_3, Y_4)\end{gathered}\right]
    \vee
    \left[\begin{gathered}\neg Y_A\\
    \left[\begin{gathered}Y_B\\\text{atleast}(3, Y_2, Y_3, Y_4)\end{gathered}\right] \vee \left[\begin{gathered}Y_C\\\text{atmost}(1, Y_2, Y_3, Y_4)\end{gathered}\right]
    \end{gathered}\right]
    \end{gather*}

or equivalently,

.. math::
    :nowrap:

    \begin{gather*}
    \text{atleast}(2, Y_1, Y_A, Y_5, Y_6)\\
    \text{exactly}(1, Y_A, Y_B, Y_C)\\
    \left[\begin{gathered}Y_A\\\text{exactly}(2, Y_2, Y_3, Y_4)\end{gathered}\right]
    \vee
    \left[\begin{gathered}Y_B\\\text{atleast}(3, Y_2, Y_3, Y_4)\end{gathered}\right] \vee \left[\begin{gathered}Y_C\\\text{atmost}(1, Y_2, Y_3, Y_4)\end{gathered}\right]
    \end{gather*}

.. code::

    m.p = LogicalConstraint(
        expr=atleast(2, m.Y[1], exactly(2, m.Y[2], m.Y[3], m.Y[4]), m.Y[5], m.Y[6]))

In the ``logical_to_linear`` transformation, we automatically convert these special disjunctions to linear form using a Big M reformulation.

Additional Examples
===================

The following models all work and are equivalent for :math:`\left[x = 0\right] \veebar \left[y = 0\right]`:

.. doctest::

   Option 1: Rule-based construction

   >>> import pyomo.environ as pyo
   >>> from pyomo.gdp import Disjunct, Disjunction
   >>> model = pyo.ConcreteModel()

   >>> model.x = pyo.Var()
   >>> model.y = pyo.Var()

   >>> # Two conditions
   >>> def _d(disjunct, flag):
   ...    model = disjunct.model()
   ...    if flag:
   ...       # x == 0
   ...       disjunct.c = pyo.Constraint(expr=model.x == 0)
   ...    else:
   ...       # y == 0
   ...       disjunct.c = pyo.Constraint(expr=model.y == 0)
   >>> model.d = Disjunct([0,1], rule=_d)

   >>> # Define the disjunction
   >>> def _c(model):
   ...    return [model.d[0], model.d[1]]
   >>> model.c = Disjunction(rule=_c)

   Option 2: Explicit disjuncts

   >>> import pyomo.environ as pyo
   >>> from pyomo.gdp import Disjunct, Disjunction
   >>> model = pyo.ConcreteModel()

   >>> model.x = pyo.Var()
   >>> model.y = pyo.Var()

   >>> model.fix_x = Disjunct()
   >>> model.fix_x.c = pyo.Constraint(expr=model.x == 0)

   >>> model.fix_y = Disjunct()
   >>> model.fix_y.c = pyo.Constraint(expr=model.y == 0)

   >>> model.c = Disjunction(expr=[model.fix_x, model.fix_y])

   Option 3: Implicit disjuncts (disjunction rule returns a list of
   expressions or a list of lists of expressions)

   >>> import pyomo.environ as pyo
   >>> from pyomo.gdp import Disjunction
   >>> model = pyo.ConcreteModel()

   >>> model.x = pyo.Var()
   >>> model.y = pyo.Var()

   >>> model.c = Disjunction(expr=[model.x == 0, model.y == 0])