Logical Modeling and Generalized Disjunctive Programs - Optimization for Decision Science

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

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.

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 = "appsi_highs"

Motivating Example: Separation¶

The following excerpts are from Section 15.7 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.

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 conjunctive normal form ( $Q_1 \land Q_2 \ \land ... \land \ Q_s$ ).

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$

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$

Logical Statements and Conjunctive 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 conjunctive normal form:

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

P_{1} \Rightarrow P_{2} \Leftrightarrow \neg P_{1} \lor P_{2} \tag{15.18}

(1)

Let’s enumerate 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:

\neg (P_{1} \land P_{2}) \Leftrightarrow \neg P_{1} \lor \neg P_{2}\tag{15.19}

(2)

Let’s enumerate 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.

\neg (P_{1} \lor P_{2}) \Leftrightarrow \neg P_{1} \land \neg P_{2}\tag{15.20}

(3)

Let’s enumerate 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:

(P_{1} \land P_{2}) \lor P_{3} \Leftrightarrow (P_{1} \lor P_{3}) \land (P_{2} \lor P_{3})\tag{15.21}

(4)

Let’s enumerate 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	true	true
true	false	true	true	true
true	false	false	false	false
false	true	true	true	true
false	true	false	false	false
false	false	true	true	true
false	false	false	false	false

Again, the example holds.

Once in conjunctive 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$

Example: Separation Sequence¶

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

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

(5)

Step 1:

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

(6)

Step 2:

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

(7)

Step 3:

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

(8)

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$

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)

(9)

Step 1:

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

(10)

which simplifies to:

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

(11)

Step 2:

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

(12)

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)

(13)

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

(14)

Simplifying:

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

(15)

Rearranging:

y_1 + y_2 - y_4 - y_5 \leq 1

(16)

On the right of the $\land$ is equivalent to:

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

(17)

Simplifying:

y_3 - y_4 - y_5 \leq 0

(18)

Modeling Disjunctions¶

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

General Notation¶

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

(19)

where:

Notation	Definition
$\lor$	The OR operator that connects a finite collection of disjunctive clauses
$D_{R}$	The set of disjunctive terms
$Y_{ik}$	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

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:

\lor_{i \in D} \begin{bmatrix} A_{i}x \leq b_{i} \end{bmatrix}

(20)

Applied to the Reactor Problem:

\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}

(21)

where $y_{1}$ represents reactor 1 and $y_{2}$ represents reactor 2.

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 (
    check_optimal_termination,
    ConcreteModel,
    Param,
    Set,
    SolverFactory,
    TransformationFactory,
    Var,
)

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(b):
        return b.P

    ## Disjunction
    # ONE disjunction over the reactors. Each DISJUNCT is itself a list -- the
    # two pressure bounds that hold together when that reactor is selected.
    # 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(
        doc="Select exactly one reactor; each brings its own pressure window"
    )
    def pressure_bounds(b):
        return [
            [b.P <= b.max_pressure[r], b.P >= b.min_pressure[r]] for r in b.reactors
        ]

    return model

Warning

The nesting matters. @model.Disjunction returns a list of disjuncts, and each disjunct is itself a list of constraints that hold together. Here that is a list of two lists: one per reactor, each carrying both of that reactor’s pressure bounds.

Flatten it, or index the decorator over model.reactors, and Pyomo builds a different model. Indexing over the reactors gives two disjunctions --- “for reactor 1, either $P \leq 10$ or $P \geq 5$ ” and “for reactor 2, either $P \leq 30$ or $P \geq 20$ ” --- neither of which restricts $P$ at all, since $P \in [5, 30]$ already satisfies one branch of each. That model admits $P = 15$ , which the disjunction forbids.

It is worth knowing that this error is invisible in the answer below: the objective is $\min P$ , and both the right model and the wrong one return $P = 5$ . Check a model by asking what it excludes, not only what it returns.

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.

Big-M Relaxation Approach¶

General Notation:

A_{i}x \leq b_{i} + M_{i}(1-y_{i}) \ , \ \forall i \in D

(22)

\sum_{i \in D} y_{i} = 1

(23)

y_{i} \in \{0,1\} \ , \ \forall i \in D

(24)

Applied to the Reactor Problem:

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

(25)

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.

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=2, 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=2, Index={1, 2}, Active=True
                    Key : Lower : Body : Upper : Active
                      1 :  -Inf :    P :  10.0 :   True
                      2 :   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[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=2, Index={1, 2}, Active=True
                    Key : Lower : Body : Upper : Active
                      1 :  -Inf :    P :  30.0 :   True
                      2 :  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 : Select exactly one reactor; each brings its own pressure window
        Size=1, Index=None, Active=True
        Key  : Disjuncts                                                        : Active : XOR
        None : ['pressure_bounds_disjuncts[0]', 'pressure_bounds_disjuncts[1]'] :   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()

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

# Solve and print the solution
results = SolverFactory(milp_solver).solve(model, tee=True)
assert check_optimal_termination(results), (
    f"Solve failed: status={results.solver.status}, "
    f"termination={results.solver.termination_condition}"
)

model.P.display()

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

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.

Convex Hull Relaxation Approach¶

General Notation:

x = \sum_{i \in D} z_{i}

(26)

A_{i} z_{i} \leq b_{i}y_{i} \ , \ \forall i \in D

(27)

\sum_{i \in D} y_{i} = 1

(28)

0 \leq z_{i} \leq Uy_{i} \ , \ \forall i \in D

(29)

y_{i} \in \{0,1\} \ , \ \forall i \in D

(30)

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

Applied to the Reactor Problem:

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

(31)

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.

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
results = SolverFactory(milp_solver).solve(model, tee=True)
assert check_optimal_termination(results), (
    f"Solve failed: status={results.solver.status}, "
    f"termination={results.solver.termination_condition}"
)

model.P.display()

Discussion Questions¶

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