Logical Modeling and Generalized Disjunctive Programs

2.1. Logical Modeling and Generalized Disjunctive Programs#

Prepared by: Prof. Alexander Dowling (adowling@nd.edu), Hailey Lynch (hlynch@nd.edu, 2023)

2.1.1. Introduction and Learning Objectives#

This notebook introduces generalized disjunctive programs through an example in Pyomo.GDP. Students will learn concepts related to Logical Modeling and Modeling Disjunctions in this notebook. These techniques will be applied to the Reactor Problem and then implemented into Pyomo. Critical thinking discussion questions will be included to connect concepts from CBE 60499.

2.1.2. Import Modules#

# Imports
import sys
if "google.colab" in sys.modules:
    !wget "https://raw.githubusercontent.com/ndcbe/optimization/main/notebooks/helper.py"
    import helper
    helper.easy_install()
else:
    sys.path.insert(0, '../')
    import helper
helper.set_plotting_style()

milp_solver = 'cbc'

2.1.3. Motivating Example: Separation#

The following excerpts are from Section 15.8 in Biegler, Grossmann, and Westerberg (1997).

From Integer Programs, we saw that modeling “choose only 1 item” is straightforward:

\( \sum_{i \in R} y_i = 1 \)

Let’s consider something more complex:

“If the absorber to recover the product is selected or the membrane is selected, then do not use cryogenic distillation.”

Yes/No binary decisions:

\(y_A\): absorber
\(y_M\): membrane
\(y_{CS}\): cryogenic separation

How to translate this logical statement into a linear constraint?

We seek a formal system to go from logical statements to linear constraints.

2.1.4. Logical Modeling#

First we will look at important logic notation that is commonly used in logical modeling. This will enable us to convert logical expressions such as disjunctive clauses (\(Q_i = P_1 \lor P_2 \ \lor ... \lor \ P_r\)) into conjuctive normal form (\(Q_1 \land Q_2 \ \land ... \land \ Q_s\)).

2.1.4.1. Symbolic Logic Notation#

Logical Operation	Logical Symbol
\(\text{OR}\)	\(\lor\)
\(\text{AND}\)	\(\land\)
\(\text{IMPLICATION}\)	\(\Rightarrow\)
\(\text{NEGATION}\)	\(\neg\)
\(\text{EQUIVALENCE}\)	\(\Leftrightarrow\)
\(\text{EXCLUSIVE OR}\)	\(\veebar\)

2.1.4.2. Vocabulary#

Literal \(P_i\) is a selection or action and \(y_i\) is the associated binary (true/false variable):

\( y_i = \begin{cases} 1 & \text{if~} P_i \text{~is~true} \\ 0 & \text{otherwise} \end{cases} \)

Negation or complement \(\neg P_i\) implies \(1 - y_i\).

Conjunctive normal form: sequence of clauses connected by AND operators

Atom or atomic formula: no deeper structure (e.g., no connectives or subformulas)

Literal: atom or its negation

Conjunctive clause: finite collection of literals connected with \( \wedge \). (The clause is true when all literals are true.)

Disjunctive clause: finite collection connected with \( \vee \). (The clause is true when at least one literal is true.)

Example: \( Q_i \) is a disjunctive clause, \( P_1 \vee P_2 \vee \dots \vee P_r \), such that \( y_1 + \dots + y_r \geq 1 \).

Need a system to convert logical statements into conjunctive normal form, e.g., \(Q_1 \wedge Q_2 \wedge ... \wedge Q_n\)

2.1.4.3. Logical Statements and Conjuctive Normal Form#

The following excerpts are from Section 15.7 in Biegler, Grossmann, and Westerberg (1997).

The three step procedure below shows how to convert logical statements into conjuctive normal form:

Step 1: Replace the implication by its equivalent disjunction. Example:

(2.1)#\[\begin{equation} P_{1} \Rightarrow P_{2} \Leftrightarrow \neg P_{1} \lor P_{2} \tag{15.18} \end{equation}\]

Click to further explore this example

Let’s enumrate to better understand this example. We will look at all possible outcomes for \(P_1\) and \(P_2\) and then assess if the left and right statements are true or false.

\(P_1\)	\(P_2\)	\(P_{1} \Rightarrow P_{2}\)	\(\neg P_{1} \lor P_{2}\)
true	true	true	true
true	false	false	false
false	true	true	true
false	false	true	true

Step 2: Distribute the negation by applying DeMorgan’s Theorem. Examples:

(2.2)#\[\begin{equation} \neg (P_{1} \land P_{2}) \Leftrightarrow \neg P_{1} \lor \neg P_{2}\tag{15.19} \end{equation}\]

Click to further explore this example

Let’s enumrate to better understand this example

\(P_1\)	\(P_2\)	\(\neg (P_{1} \land P_{2})\)	\(\neg P_{1} \lor \neg P_{2}\)
true	true	false	false
true	false	true	true
false	true	true	true
false	false	true	true

The example holds.

(2.3)#\[\begin{equation} \neg (P_{1} \lor P_{2}) \Leftrightarrow \neg P_{1} \land \neg P_{2}\tag{15.20} \end{equation}\]

Click to further explore this example

Let’s enumrate to better understand this example

\(P_1\)	\(P_2\)	\(\neg (P_{1} \lor P_{2})\)	\(\neg P_{1} \land \neg P_{2}\)
true	true	false	false
true	false	false	false
false	true	false	false
false	false	true	true

The example holds.

Step 3: Recursively distribute the \(\textbf{OR}\) over the \(\textbf{AND}\) by using the following equivalence. Example:

(2.4)#\[\begin{equation} (P_{1} \land P_{2}) \lor P_{3} \Leftrightarrow (P_{1} \lor P_{3}) \land (P_{2} \lor P_{3})\tag{15.21} \end{equation}\]

Click to further explore this example

Let’s enumrate to better understand this example

\(P_1\)	\(P_2\)	\(P_3\)	\((P_{1} \land P_{2}) \lor P_{3}\)	\((P_{1} \lor P_{3}) \land (P_{2} \lor P_{3})\)
true	true	true	true	true
true	true	false	false	false
true	false	true	true	true
true	false	false	false	false
false	true	true	true	true
false	true	false	false	false
false	false	true	false	false
false	false	false	false	false

Again, the example holds.

Once in conjuctive normal form, we can apply the following rules:

Logical Relation	Comments	Boolean Expression	Representation as Linear Inequalities
Logical OR		\(P_{1}\) \(\lor\) \(P_{2}\) \(\lor \ .. \lor\) \(P_{r}\)	\(y_{1} + y_{2}\) \( + \ ..\ +\) \(y_{r} \geq 1\)
Logical AND		\(P_{1}\) \(\land\) \(P_{2}\) \(\land \ .. \land\) \(P_{r}\)	\(y_{1} \geq 1\\ y_{2} \geq 1\\ .. \\ y_{r} \geq 1\)
Implication	\(P_{1} \Rightarrow P_{2}\)	\(\neg P_{1} \lor P_{2}\)	\(1-y_{1} + y_{2} \geq 1\)
Equivalence	\(P_{1}\) iff \(P_{2} (P_{1} \Rightarrow P_{2}) \land (P_{2} \Rightarrow P_{1})\)	\((\neg P_{1} \lor P_{2})\) \(\land\) \((\neg P_{2} \lor P_{1})\)	\(y_{1} = y_{2}\)
Exclusive OR	Exactly one of the variables is true	\(P_{1} \veebar P_{2} \) \(\veebar \ .. \ \veebar\) \(P_{r}\)	\(y_{1} + y_{2}\) \(+ \ ..\ +\) \(y_{r} = 1\)

2.1.4.4. Example: Separation Sequence#

Let’s revist the example from the top of the notebook. Reformulate:

\[ P_A \lor P_M \implies \neg P_{CS} \]

Step 1:

\[ \neg (P_A \lor P_M) \lor \neg P_{CS} \]

Step 2:

\[ (\neg P_A \land \neg P_M) \lor \neg P_{CS} \]

Step 3:

\[ (\neg P_A \lor \neg P_{CS}) \land (\neg P_M \lor \neg P_{CS}) \]

Now substituting:

\( \neg P_A \rightarrow 1 - y_A, \quad \neg P_{CS} \rightarrow 1 - y_{CS} \)

\( \neg P_M \rightarrow 1 - y_M, \quad \neg P_{CS} \rightarrow 1 - y_{CS} \)

We get:

\( (1 - y_A) + (1 - y_{CS}) \geq 1 \)

\( (1 - y_M) + (1 - y_{CS}) \geq 1 \)

Rearrange:

\( y_A + y_{CS} \leq 1 \)

\( y_M + y_{CS} \leq 1 \)

2.1.4.5. Example: Assembling Components#

If you use (parts 1 and 2) or part 3, then you must also use parts 4 or 5.

\(P_i\) true means “use part \(i\)” (corresponding to \(y_i\)).

\[ (P_1 \land P_2) \lor P_3 \implies (P_4 \lor P_5) \]

Step 1:

\[ \neg ((P_1 \land P_2) \lor P_3) \lor (P_4 \lor P_5) \]

which simplifies to:

\[ \neg (P_1 \land P_2) \land \neg P_3 \lor (P_4 \lor P_5) \]

Step 2:

\[ (\neg P_1 \lor \neg P_2) \land \neg P_3 \lor (P_4 \lor P_5) \]

Step 3:

\[ (\neg P_1 \lor \neg P_2 \lor P_4 \lor P_5) \land (\neg P_3 \lor P_4 \lor P_5) \]

We now have our statement in conjunctive normal form.

On the left of the \(\land\) is equivalent to:

\[ (1 - y_1) + (1 - y_2) + y_4 + y_5 \geq 1 \]

Simplifying:

\[ -y_1 - y_2 + y_4 + y_5 \geq -1 \]

Rearranging:

\[ y_1 + y_2 - y_4 - y_5 \leq 1 \]

On the right of the \(\land\) is equivalent to:

\[ (1 - y_3) + y_4 + y_5 \geq 1 \]

Simplifying:

\[ y_3 - y_4 - y_5 \leq 0 \]

2.1.5. Modeling Disjunctions#

When modeling disjunctions, we will have to represent logical constraints that involve continuous variables.

2.1.5.1. General Notation#

(2.5)#\[\begin{equation} \lor_{i \in D_{R}} \begin{bmatrix} Y_{ik}\\ n_{ik}(x)\leq0\\ c_{k} = \mu_{ik} \end{bmatrix} \ \ \ , \ \ \ \Omega(Y)= \text{True} \end{equation}\]

where:

Notation	Definition
\(\lor\)	The OR operator that connects a finite collection of disjunctive clauses
\(D_{R}\)	The set of disjunctive terms
\(Y_{ij}\)	Boolean “indicator variable”
\(n_{ik}(x)\leq0\)	Constraint enforced when \(Y_{ik}\) is true
\(\mu_{ik}\)	Parameter values when indicator is true
\(\Omega(Y)\)	Additional logical constraints

2.1.6. Example: The Reactor Problem#

The following excerpts are from Section 15.8 in Biegler, Grossmann, and Westerberg (1997).

This modeling disjunctions example involves selecting between two reactors:

If reactor 1 is selected, then pressure \(P\) must be between \(5\) and \(10\) atmospheres.
If reactor 2 is selected, then pressure \(P\) must be between \(20\) and \(30\) atmospheres.

Linear Disjunction Form:

(2.6)#\[\begin{equation} \lor_{i \in D} \begin{bmatrix} A_{i}x \leq b_{i}(x) \end{bmatrix} \end{equation}\]

Applied to the Reactor Problem:

(2.7)#\[\begin{equation} \begin{bmatrix} y_1\\ P \leq 10\\ -P \leq -5 \end{bmatrix} \lor \begin{bmatrix} y_2\\ P \leq 30\\ -P \leq -20 \end{bmatrix} \end{equation}\]

where \(y_{1}\) represents reactor 1 and \(y_{2}\) represents reactor 2.

2.1.6.1. Define Model in Pyomo with GDP#

First we will define the model (including disjunctions) for the Reactor Problem in Pyomo.

'''
Instead of using
# import pyomo.environ as pyo
We can import specific functions/objects
'''
from pyomo.environ import (ConcreteModel, NonNegativeReals, Objective, Param,
                           Set, SolverFactory, TransformationFactory, Var, value)

def create_model():
    '''
    Build the reactor problem model.
    
    Return:
    model: Pyomo model
    
    '''
    ## Model
    model = ConcreteModel(name="Selecting a reactor")

    ## Sets
    # Initialized for reactor 1 (1) and reactor 2 (2)
    model.reactors = Set(initialize=[1,2])

    ## Parameters
    # Initialized with a dictionary where the keys are 1 and 2 (the reactors)
    # for the minimum and maximum pressure values (atm)
    model.min_pressure = Param(model.reactors,initialize={1:5,2:20}) 
    model.max_pressure = Param(model.reactors,initialize={1:10,2:30})
    
    ## Variables
    # Reactor pressure bounded between the lower bound (5 atm) and upper bound (30 atm)
    model.P = Var(bounds=(5,30),doc='Reactor pressure (atm)')

    ## Adding an objective for the example
    @model.Objective()
    def objective(model):
      return model.P 

    ## Disjunctions
    # Objective is bounded by the maximum and minimum pressure values
    # Note: Pyomo.GDP by default treats the disjunction as a xor (choose only one)
    # https://pyomo.readthedocs.io/en/latest/modeling_extensions/gdp/modeling.html
    @model.Disjunction(
        model.reactors,
        doc="Pressure bounds for different reactor selections")
    def pressure_bounds(m, r):
        return [
            m.P <= m.max_pressure[r],
            m.P >= m.min_pressure[r],
        ]

    return model

2.1.6.2. Transform and Solve with Big-M Relaxation#

The following excerpts are from Section 15.8 in Biegler, Grossmann, and Westerberg (1997).

Use “Big-M” constraints to convert linear disjunctions into mixed-integer constraints to represent logic with continuous variables.

2.1.6.3. Big-M Relaxation Approach#

General Notation:

(2.8)#\[\begin{equation} A_{i}(x) \leq b_{i} + M_{i}(1-y_{i}) \ , \ \forall i \in D \end{equation}\]

(2.9)#\[\begin{equation} \sum_{i \in D} y_{i} = 1 \end{equation}\]

(2.10)#\[\begin{equation} y_{i} \in \{0,1\} \ , \ \forall i \in D \end{equation}\]

Applied to the Reactor Problem:

(2.11)#\[\begin{equation} P \leq 10 + M_{1}(1-y_{1})\\ -P \leq -5 + M_{1}(1-y_{1})\\ P \leq 30 + M_{2}(1-y_{2})\\ -P \leq -20 + M_{2}(1-y_{2})\\ y_{1} + y_{2} = 1 \end{equation}\]

When the \(y\)’s are considered continuous variables, weak bounds for the objective function are formed for large values such as:

\(M_{1} = 100 \ \ \text{and} \ \ M_{2} = 100\)

Main Idea: Considering the special case where \(h_{i}(x) = A_{i}x - b_{i} \leq 0,\)
\(M_{i}\) is sufficiently large to relax \(h_{i}(x) \leq 0\) when \(y_{i}=0\)

Key Takeaways:

If \(M\) is too large, we can get a “weak relaxation” because integer programming algorithms need more iterations.
If \(M\) is too small, we can get unintended bounds.

Big-M is the best to use if the problem is small.

2.1.6.4. Big-M Implementation in Pyomo#

First we will create and print the model.

# Creating the model
model = create_model()

# Printing the model
model.pprint()

1 Set Declarations
    reactors : Size=1, Index=None, Ordered=Insertion
        Key  : Dimen : Domain : Size : Members
        None :     1 :    Any :    2 : {1, 2}

2 Param Declarations
    max_pressure : Size=2, Index=reactors, Domain=Any, Default=None, Mutable=False
        Key : Value
          1 :    10
          2 :    30
    min_pressure : Size=2, Index=reactors, Domain=Any, Default=None, Mutable=False
        Key : Value
          1 :     5
          2 :    20

1 Var Declarations
    P : Reactor pressure (atm)
        Size=1, Index=None
        Key  : Lower : Value : Upper : Fixed : Stale : Domain
        None :     5 :  None :    30 : False :  True :  Reals

1 Objective Declarations
    objective : Size=1, Index=None, Active=True
        Key  : Active : Sense    : Expression
        None :   True : minimize :          P

1 Disjunct Declarations
    pressure_bounds_disjuncts : Size=4, Index=Any, Active=True
        pressure_bounds_disjuncts[0] : Active=True
            1 Var Declarations
                binary_indicator_var : Size=1, Index=None
                    Key  : Lower : Value : Upper : Fixed : Stale : Domain
                    None :     0 :  None :     1 : False :  True : Binary

            1 Constraint Declarations
                constraint : Size=1, Index={1}, Active=True
                    Key : Lower : Body : Upper : Active
                      1 :  -Inf :    P :  10.0 :   True

            1 BooleanVar Declarations
                indicator_var : Size=1, Index=None
                    Key  : Value : Fixed : Stale
                    None :  None : False :  True

            1 LogicalConstraint Declarations
                propositions : Size=0, Index={}, Active=True
                    Key : Body : Active

            4 Declarations: indicator_var binary_indicator_var constraint propositions
        pressure_bounds_disjuncts[1] : Active=True
            1 Var Declarations
                binary_indicator_var : Size=1, Index=None
                    Key  : Lower : Value : Upper : Fixed : Stale : Domain
                    None :     0 :  None :     1 : False :  True : Binary

            1 Constraint Declarations
                constraint : Size=1, Index={1}, Active=True
                    Key : Lower : Body : Upper : Active
                      1 :   5.0 :    P :  +Inf :   True

            1 BooleanVar Declarations
                indicator_var : Size=1, Index=None
                    Key  : Value : Fixed : Stale
                    None :  None : False :  True

            1 LogicalConstraint Declarations
                propositions : Size=0, Index={}, Active=True
                    Key : Body : Active

            4 Declarations: indicator_var binary_indicator_var constraint propositions
        pressure_bounds_disjuncts[2] : Active=True
            1 Var Declarations
                binary_indicator_var : Size=1, Index=None
                    Key  : Lower : Value : Upper : Fixed : Stale : Domain
                    None :     0 :  None :     1 : False :  True : Binary

            1 Constraint Declarations
                constraint : Size=1, Index={1}, Active=True
                    Key : Lower : Body : Upper : Active
                      1 :  -Inf :    P :  30.0 :   True

            1 BooleanVar Declarations
                indicator_var : Size=1, Index=None
                    Key  : Value : Fixed : Stale
                    None :  None : False :  True

            1 LogicalConstraint Declarations
                propositions : Size=0, Index={}, Active=True
                    Key : Body : Active

            4 Declarations: indicator_var binary_indicator_var constraint propositions
        pressure_bounds_disjuncts[3] : Active=True
            1 Var Declarations
                binary_indicator_var : Size=1, Index=None
                    Key  : Lower : Value : Upper : Fixed : Stale : Domain
                    None :     0 :  None :     1 : False :  True : Binary

            1 Constraint Declarations
                constraint : Size=1, Index={1}, Active=True
                    Key : Lower : Body : Upper : Active
                      1 :  20.0 :    P :  +Inf :   True

            1 BooleanVar Declarations
                indicator_var : Size=1, Index=None
                    Key  : Value : Fixed : Stale
                    None :  None : False :  True

            1 LogicalConstraint Declarations
                propositions : Size=0, Index={}, Active=True
                    Key : Body : Active

            4 Declarations: indicator_var binary_indicator_var constraint propositions

1 Disjunction Declarations
    pressure_bounds : Pressure bounds for different reactor selections
        Size=2, Index=reactors, Active=True
        Key : Disjuncts                                                        : Active : XOR
          1 : ['pressure_bounds_disjuncts[0]', 'pressure_bounds_disjuncts[1]'] :   True : True
          2 : ['pressure_bounds_disjuncts[2]', 'pressure_bounds_disjuncts[3]'] :   True : True

7 Declarations: reactors min_pressure max_pressure P objective pressure_bounds pressure_bounds_disjuncts

Next, let’s transform using Big-M and print the model again.

# Applying Big-M relaxation to the model
# Add your solution here

# Printing
model.pprint()

1 Set Declarations
    reactors : Size=1, Index=None, Ordered=Insertion
        Key  : Dimen : Domain : Size : Members
        None :     1 :    Any :    2 : {1, 2}

2 Param Declarations
    max_pressure : Size=2, Index=reactors, Domain=Any, Default=None, Mutable=False
        Key : Value
          1 :    10
          2 :    30
    min_pressure : Size=2, Index=reactors, Domain=Any, Default=None, Mutable=False
        Key : Value
          1 :     5
          2 :    20

1 Var Declarations
    P : Reactor pressure (atm)
        Size=1, Index=None
        Key  : Lower : Value : Upper : Fixed : Stale : Domain
        None :     5 :  None :    30 : False :  True :  Reals

1 Objective Declarations
    objective : Size=1, Index=None, Active=True
        Key  : Active : Sense    : Expression
        None :   True : minimize :          P

1 Block Declarations
    _pyomo_gdp_bigm_reformulation : Size=1, Index=None, Active=True
        1 Constraint Declarations
            pressure_bounds_xor : Size=2, Index=Any, Active=True
                Key : Lower : Body                                                                                                  : Upper : Active
                  1 :   1.0 : pressure_bounds_disjuncts[0].binary_indicator_var + pressure_bounds_disjuncts[1].binary_indicator_var :   1.0 :   True
                  2 :   1.0 : pressure_bounds_disjuncts[2].binary_indicator_var + pressure_bounds_disjuncts[3].binary_indicator_var :   1.0 :   True

        1 Block Declarations
            relaxedDisjuncts : Size=4, Index=NonNegativeIntegers, Active=True
                _pyomo_gdp_bigm_reformulation.relaxedDisjuncts[0] : Active=True
                    1 Constraint Declarations
                        transformedConstraints : Size=1, Index=Any, Active=True
                            Key                          : Lower : Body                                                             : Upper : Active
                            ('constraint[1]_0', 1, 'ub') :  -Inf : P - 20.0*(1 - pressure_bounds_disjuncts[0].binary_indicator_var) :  10.0 :   True

                    1 Block Declarations
                        localVarReferences : Size=1, Index=None, Active=True
                            1 Var Declarations
                                binary_indicator_var : Size=1, Index=UnindexedComponent_ReferenceSet, ReferenceTo=pressure_bounds_disjuncts[0].binary_indicator_var
                                    Key  : Lower : Value : Upper : Fixed : Stale : Domain
                                    None :     0 :  None :     1 : False :  True : Binary

                            1 Declarations: binary_indicator_var

                    2 Declarations: transformedConstraints localVarReferences
                _pyomo_gdp_bigm_reformulation.relaxedDisjuncts[1] : Active=True
                    1 Constraint Declarations
                        transformedConstraints : Size=1, Index=Any, Active=True
                            Key                          : Lower : Body                                                            : Upper : Active
                            ('constraint[1]_0', 1, 'lb') :   5.0 : P - 0.0*(1 - pressure_bounds_disjuncts[1].binary_indicator_var) :  +Inf :   True

                    1 Block Declarations
                        localVarReferences : Size=1, Index=None, Active=True
                            1 Var Declarations
                                binary_indicator_var : Size=1, Index=UnindexedComponent_ReferenceSet, ReferenceTo=pressure_bounds_disjuncts[1].binary_indicator_var
                                    Key  : Lower : Value : Upper : Fixed : Stale : Domain
                                    None :     0 :  None :     1 : False :  True : Binary

                            1 Declarations: binary_indicator_var

                    2 Declarations: transformedConstraints localVarReferences
                _pyomo_gdp_bigm_reformulation.relaxedDisjuncts[2] : Active=True
                    1 Constraint Declarations
                        transformedConstraints : Size=1, Index=Any, Active=True
                            Key                          : Lower : Body                                                            : Upper : Active
                            ('constraint[1]_0', 1, 'ub') :  -Inf : P - 0.0*(1 - pressure_bounds_disjuncts[2].binary_indicator_var) :  30.0 :   True

                    1 Block Declarations
                        localVarReferences : Size=1, Index=None, Active=True
                            1 Var Declarations
                                binary_indicator_var : Size=1, Index=UnindexedComponent_ReferenceSet, ReferenceTo=pressure_bounds_disjuncts[2].binary_indicator_var
                                    Key  : Lower : Value : Upper : Fixed : Stale : Domain
                                    None :     0 :  None :     1 : False :  True : Binary

                            1 Declarations: binary_indicator_var

                    2 Declarations: transformedConstraints localVarReferences
                _pyomo_gdp_bigm_reformulation.relaxedDisjuncts[3] : Active=True
                    1 Constraint Declarations
                        transformedConstraints : Size=1, Index=Any, Active=True
                            Key                          : Lower : Body                                                             : Upper : Active
                            ('constraint[1]_0', 1, 'lb') :  20.0 : P + 15.0*(1 - pressure_bounds_disjuncts[3].binary_indicator_var) :  +Inf :   True

                    1 Block Declarations
                        localVarReferences : Size=1, Index=None, Active=True
                            1 Var Declarations
                                binary_indicator_var : Size=1, Index=UnindexedComponent_ReferenceSet, ReferenceTo=pressure_bounds_disjuncts[3].binary_indicator_var
                                    Key  : Lower : Value : Upper : Fixed : Stale : Domain
                                    None :     0 :  None :     1 : False :  True : Binary

                            1 Declarations: binary_indicator_var

                    2 Declarations: transformedConstraints localVarReferences

        2 Declarations: relaxedDisjuncts pressure_bounds_xor

1 Disjunct Declarations
    pressure_bounds_disjuncts : Size=4, Index=Any, Active=False
        pressure_bounds_disjuncts[0] : Active=False
            1 Var Declarations
                binary_indicator_var : Size=1, Index=None
                    Key  : Lower : Value : Upper : Fixed : Stale : Domain
                    None :     0 :  None :     1 : False :  True : Binary

            1 Constraint Declarations
                constraint : Size=1, Index={1}, Active=True
                    Key : Lower : Body : Upper : Active
                      1 :  -Inf :    P :  10.0 :  False

            1 BooleanVar Declarations
                indicator_var : Size=1, Index=None
                    Key  : Value : Fixed : Stale
                    None :  None : False :  True

            1 LogicalConstraint Declarations
                propositions : Size=0, Index={}, Active=False
                    Key : Body : Active

            4 Declarations: indicator_var binary_indicator_var constraint propositions
        pressure_bounds_disjuncts[1] : Active=False
            1 Var Declarations
                binary_indicator_var : Size=1, Index=None
                    Key  : Lower : Value : Upper : Fixed : Stale : Domain
                    None :     0 :  None :     1 : False :  True : Binary

            1 Constraint Declarations
                constraint : Size=1, Index={1}, Active=True
                    Key : Lower : Body : Upper : Active
                      1 :   5.0 :    P :  +Inf :  False

            1 BooleanVar Declarations
                indicator_var : Size=1, Index=None
                    Key  : Value : Fixed : Stale
                    None :  None : False :  True

            1 LogicalConstraint Declarations
                propositions : Size=0, Index={}, Active=False
                    Key : Body : Active

            4 Declarations: indicator_var binary_indicator_var constraint propositions
        pressure_bounds_disjuncts[2] : Active=False
            1 Var Declarations
                binary_indicator_var : Size=1, Index=None
                    Key  : Lower : Value : Upper : Fixed : Stale : Domain
                    None :     0 :  None :     1 : False :  True : Binary

            1 Constraint Declarations
                constraint : Size=1, Index={1}, Active=True
                    Key : Lower : Body : Upper : Active
                      1 :  -Inf :    P :  30.0 :  False

            1 BooleanVar Declarations
                indicator_var : Size=1, Index=None
                    Key  : Value : Fixed : Stale
                    None :  None : False :  True

            1 LogicalConstraint Declarations
                propositions : Size=0, Index={}, Active=False
                    Key : Body : Active

            4 Declarations: indicator_var binary_indicator_var constraint propositions
        pressure_bounds_disjuncts[3] : Active=False
            1 Var Declarations
                binary_indicator_var : Size=1, Index=None
                    Key  : Lower : Value : Upper : Fixed : Stale : Domain
                    None :     0 :  None :     1 : False :  True : Binary

            1 Constraint Declarations
                constraint : Size=1, Index={1}, Active=True
                    Key : Lower : Body : Upper : Active
                      1 :  20.0 :    P :  +Inf :  False

            1 BooleanVar Declarations
                indicator_var : Size=1, Index=None
                    Key  : Value : Fixed : Stale
                    None :  None : False :  True

            1 LogicalConstraint Declarations
                propositions : Size=0, Index={}, Active=False
                    Key : Body : Active

            4 Declarations: indicator_var binary_indicator_var constraint propositions

1 Disjunction Declarations
    pressure_bounds : Pressure bounds for different reactor selections
        Size=2, Index=reactors, Active=False
        Key : Disjuncts                                                        : Active : XOR
          1 : ['pressure_bounds_disjuncts[0]', 'pressure_bounds_disjuncts[1]'] :  False : True
          2 : ['pressure_bounds_disjuncts[2]', 'pressure_bounds_disjuncts[3]'] :  False : True

8 Declarations: reactors min_pressure max_pressure P objective pressure_bounds pressure_bounds_disjuncts _pyomo_gdp_bigm_reformulation

Finally, we’ll solve the model and examine the solution.

# Solve and print the solution
SolverFactory(milp_solver).solve(model, tee=True)

model.P.display()

Welcome to the CBC MILP Solver 
Version: 2.10.10 
Build Date: Jun  7 2023 

command line - /Users/adowling/.idaes/bin/cbc -printingOptions all -import /var/folders/3w/vr4xmyqs451dg23xk88pqcg00000gq/T/tmpx99o357i.pyomo.lp -stat=1 -solve -solu /var/folders/3w/vr4xmyqs451dg23xk88pqcg00000gq/T/tmpx99o357i.pyomo.soln (default strategy 1)
Option for printingOptions changed from normal to all
Presolve 0 (-6) rows, 0 (-5) columns and 0 (-10) elements
Statistics for presolved model
Original problem has 4 integers (4 of which binary)


Problem has 0 rows, 0 columns (0 with objective) and 0 elements
Column breakdown:
0 of type 0.0->inf, 0 of type 0.0->up, 0 of type lo->inf, 
0 of type lo->up, 0 of type free, 0 of type fixed, 
0 of type -inf->0.0, 0 of type -inf->up, 0 of type 0.0->1.0 
Row breakdown:
0 of type E 0.0, 0 of type E 1.0, 0 of type E -1.0, 
0 of type E other, 0 of type G 0.0, 0 of type G 1.0, 
0 of type G other, 0 of type L 0.0, 0 of type L 1.0, 
0 of type L other, 0 of type Range 0.0->1.0, 0 of type Range other, 
0 of type Free 
Continuous objective value is 5 - 0.00 seconds
Cgl0004I processed model has 0 rows, 0 columns (0 integer (0 of which binary)) and 0 elements
Cbc3007W No integer variables - nothing to do
Cuts at root node changed objective from 5 to -1.79769e+308
Probing was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Gomory was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Knapsack was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Clique was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
MixedIntegerRounding2 was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
FlowCover was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
TwoMirCuts was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
ZeroHalf was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)

Result - Optimal solution found

Objective value:                5.00000000
Enumerated nodes:               0
Total iterations:               0
Time (CPU seconds):             0.00
Time (Wallclock seconds):       0.00

Total time (CPU seconds):       0.00   (Wallclock seconds):       0.01

P : Reactor pressure (atm)
    Size=1, Index=None
    Key  : Lower : Value : Upper : Fixed : Stale : Domain
    None :     5 :   5.0 :    30 : False : False :  Reals

2.1.6.5. Transform and Solve with Convex Hull Relaxation#

The following excerpts are from Section 15.8 in Biegler, Grossmann, and Westerberg (1997).

Convex hull can be used if we don’t want to implement Big-M parameters. This approach requires separating the continuous variables into its components.

2.1.6.6. Convex Hull Relaxation Approach#

General Notation:

\[ \begin{equation} x = \sum_{i \in D} z_{i} \end{equation} \]

\[ \begin{equation} A_{i} z_{i} \leq b_{i}y_{i} \ , \ \forall i \in D \end{equation} \]

\[ \begin{equation} \sum_{i \in D} y_{i} = 1 \end{equation} \]

\[ \begin{equation} 0 \leq z_{i} \leq Uy_{i} \ , \ \forall i \in D \end{equation} \]

\[ \begin{equation} y_{i} = \{0,1\} \ , \ \forall i \in D \end{equation} \]

\(z_{i}\): continuous variables separated into as many new variables as there are terms for the disjunctions.

Applied to the Reactor Problem:

(2.12)#\[\begin{equation} P = P_{1} + P_{2} \\ P_{1} \leq 10y_{1}\\ P_{2} \leq 30y_{2} \\ -P_{1} \leq -5y_{1}\\ -P_{2} \leq -20y_{2}\\ y_{1} + y_{2} = 1 \end{equation}\]

Key Takeaways:

(+) Constraints do not require Big-M parameters which produce a tight linear programming relaxation.
(–) A larger number of variables and constraints is required.

Convex hull is better to use over Big-M if the problem is large.

2.1.6.7. Convex Hull Implementation in Pyomo#

We will repeat the procedure above but using Convex Hull now.

# Creating the model
model = create_model()

# Applying convex hull relaxation to the model
# Add your solution here

# Solve and print the solution
SolverFactory(milp_solver).solve(model, tee=True)

model.P.display()

Welcome to the CBC MILP Solver 
Version: 2.10.10 
Build Date: Jun  7 2023 

command line - /Users/adowling/.idaes/bin/cbc -printingOptions all -import /var/folders/3w/vr4xmyqs451dg23xk88pqcg00000gq/T/tmppgpvdtjh.pyomo.lp -stat=1 -solve -solu /var/folders/3w/vr4xmyqs451dg23xk88pqcg00000gq/T/tmppgpvdtjh.pyomo.soln (default strategy 1)
Option for printingOptions changed from normal to all
Presolve 10 (-6) rows, 5 (-4) columns and 24 (-10) elements
Statistics for presolved model
Original problem has 4 integers (4 of which binary)
Presolved problem has 2 integers (2 of which binary)
==== 4 zero objective 2 different
4 variables have objective of 0
1 variables have objective of 1
==== absolute objective values 2 different
4 variables have objective of 0
1 variables have objective of 1
==== for integers 2 zero objective 1 different
2 variables have objective of 0
==== for integers absolute objective values 1 different
2 variables have objective of 0
===== end objective counts


Problem has 10 rows, 5 columns (1 with objective) and 24 elements
Column breakdown:
0 of type 0.0->inf, 2 of type 0.0->up, 0 of type lo->inf, 
1 of type lo->up, 0 of type free, 0 of type fixed, 
0 of type -inf->0.0, 0 of type -inf->up, 2 of type 0.0->1.0 
Row breakdown:
0 of type E 0.0, 0 of type E 1.0, 0 of type E -1.0, 
0 of type E other, 0 of type G 0.0, 0 of type G 1.0, 
0 of type G other, 5 of type L 0.0, 0 of type L 1.0, 
5 of type L other, 0 of type Range 0.0->1.0, 0 of type Range other, 
0 of type Free 
Continuous objective value is 5 - 0.00 seconds
Cgl0003I 0 fixed, 0 tightened bounds, 1 strengthened rows, 0 substitutions
Cgl0003I 0 fixed, 0 tightened bounds, 1 strengthened rows, 0 substitutions
Cgl0004I processed model has 8 rows, 5 columns (2 integer (2 of which binary)) and 21 elements
Cbc0038I Initial state - 0 integers unsatisfied sum - 0
Cbc0038I Solution found of 5
Cbc0038I Relaxing continuous gives 5
Cbc0038I Before mini branch and bound, 2 integers at bound fixed and 3 continuous
Cbc0038I Mini branch and bound did not improve solution (0.00 seconds)
Cbc0038I After 0.00 seconds - Feasibility pump exiting with objective of 5 - took 0.00 seconds
Cbc0012I Integer solution of 5 found by feasibility pump after 0 iterations and 0 nodes (0.00 seconds)
Cbc0001I Search completed - best objective 5, took 0 iterations and 0 nodes (0.00 seconds)
Cbc0035I Maximum depth 0, 0 variables fixed on reduced cost
Cuts at root node changed objective from 5 to 5
Probing was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Gomory was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Knapsack was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Clique was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
MixedIntegerRounding2 was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
FlowCover was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
TwoMirCuts was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
ZeroHalf was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)

Result - Optimal solution found

Objective value:                5.00000000
Enumerated nodes:               0
Total iterations:               0
Time (CPU seconds):             0.00
Time (Wallclock seconds):       0.00

Total time (CPU seconds):       0.00   (Wallclock seconds):       0.00

P : Reactor pressure (atm)
    Size=1, Index=None
    Key  : Lower : Value : Upper : Fixed : Stale : Domain
    None :     5 :   5.0 :    30 : False : False :  Reals

2.1.7. Discussion Questions#

How do we create a system to go from logical expressions to linear constraints?
Are conjuctive or disjunctive clauses more common? Why might this be the case?
If \(y_{i}=0\) and \(y_{i}=1\), what happens when a Big-M parameter is introduced in the general notation?
When will the convex hull formulation simplify?

Logical Modeling and Generalized Disjunctive Programs

Contents

2.1. Logical Modeling and Generalized Disjunctive Programs#

2.1.1. Introduction and Learning Objectives#

2.1.2. Import Modules#

2.1.3. Motivating Example: Separation#

2.1.4. Logical Modeling#

2.1.4.1. Symbolic Logic Notation#

2.1.4.2. Vocabulary#

2.1.4.3. Logical Statements and Conjuctive Normal Form#

2.1.4.4. Example: Separation Sequence#

2.1.4.5. Example: Assembling Components#

2.1.5. Modeling Disjunctions#

2.1.5.1. General Notation#

2.1.6. Example: The Reactor Problem#

2.1.6.1. Define Model in Pyomo with GDP#

2.1.6.2. Transform and Solve with Big-M Relaxation#

2.1.6.3. Big-M Relaxation Approach#

2.1.6.4. Big-M Implementation in Pyomo#

2.1.6.5. Transform and Solve with Convex Hull Relaxation#

2.1.6.6. Convex Hull Relaxation Approach#

2.1.6.7. Convex Hull Implementation in Pyomo#

2.1.7. Discussion Questions#