From 161ce1daf364084298a77c671fbf727aecc1c8d9 Mon Sep 17 00:00:00 2001
From: Albert Lee <lee4382@purdue.edu>
Date: Mon, 13 Nov 2023 10:37:59 -0500
Subject: [PATCH 1/9] Added the small batch on the gdplib

---
 gdplib/small_batch/gdp_small_batch.py | 460 ++++++++++++++++++++++++++
 1 file changed, 460 insertions(+)
 create mode 100644 gdplib/small_batch/gdp_small_batch.py

diff --git a/gdplib/small_batch/gdp_small_batch.py b/gdplib/small_batch/gdp_small_batch.py
new file mode 100644
index 0000000..b9cdce5
--- /dev/null
+++ b/gdplib/small_batch/gdp_small_batch.py
@@ -0,0 +1,460 @@
+"""
+gdp_small_batch.py
+The gdp_small_batch.py module contains the GDP model for the small batch problem based on the Kocis and Grossmann (1988) paper.
+The problem is based on the Example 4 of the paper.
+
+References:
+    - Kocis, G. R.; Grossmann, I. E. Global Optimization of Nonconvex Mixed-Integer Nonlinear Programming (MINLP) Problems in Process Synthesis. Ind. Eng. Chem. Res. 1988, 27 (8), 1407–1421. 
+"""
+import os
+
+import pyomo.environ as pe
+from pyomo.core.base.misc import display
+from pyomo.core.plugins.transform.logical_to_linear import (
+    update_boolean_vars_from_binary,
+)
+from pyomo.gdp import Disjunct, Disjunction
+from pyomo.opt.base.solvers import SolverFactory
+
+
+def build_small_batch():
+    """
+    The function build the GDP model for the small batch problem.
+
+    References:
+    - Kocis, G. R.; Grossmann, I. E. Global Optimization of Nonconvex Mixed-Integer Nonlinear Programming (MINLP) Problems in Process Synthesis. Ind. Eng. Chem. Res. 1988, 27 (8), 1407–1421.
+
+    Args:
+        None
+
+    Returns:
+        m (pyomo.ConcreteModel): The GDP model for the small batch problem is created.
+    """
+    NK = 3
+
+    # Model
+    m = pe.ConcreteModel()
+
+    # Sets
+    m.i = pe.Set(
+        initialize=['a', 'b'], doc='Set of products'
+    )  # Set of products, i = a, b
+    m.j = pe.Set(
+        initialize=['mixer', 'reactor', 'centrifuge']
+    )  # Set of stages, j = mixer, reactor, centrifuge
+    m.k = pe.RangeSet(NK)  # Set of potential number of parallel units, k = 1, 2, 3
+
+    # Parameters and Scalars
+
+    m.h = pe.Param(
+        initialize=6000, doc='Horizon time [hr]'
+    )  # Horizon time  (available time) [hr]
+    m.vlow = pe.Param(
+        initialize=250, doc='Lower bound for size of batch unit [L]'
+    )  # Lower bound for size of batch unit [L]
+    m.vupp = pe.Param(
+        initialize=2500, doc='Upper bound for size of batch unit [L]'
+    )  # Upper bound for size of batch unit [L]
+
+    # Demand of product i
+    m.q = pe.Param(
+        m.i,
+        initialize={'a': 200000, 'b': 150000},
+        doc='Production rate of the product [kg]',
+    )
+    # Cost coefficient for batch units
+    m.alpha = pe.Param(
+        m.j,
+        initialize={'mixer': 250, 'reactor': 500, 'centrifuge': 340},
+        doc='Cost coefficient for batch units [$/L^beta*No. of units]]',
+    )
+    # Cost exponent for batch units
+    m.beta = pe.Param(
+        m.j,
+        initialize={'mixer': 0.6, 'reactor': 0.6, 'centrifuge': 0.6},
+        doc='Cost exponent for batch units',
+    )
+
+    def coeff_init(m, k):
+        """
+        Coefficient for number of parallel units.
+
+        Args:
+            m (pyomo.ConcreteModel): small batch GDP model
+            k (int): number of parallel units
+
+        Returns:
+            Coefficient for number of parallel units.
+        """
+        return pe.log(k)
+
+    # Represent number of parallel units
+    m.coeff = pe.Param(
+        m.k, initialize=coeff_init, doc='Coefficient for number of parallel units'
+    )
+
+    s_init = {
+        ('a', 'mixer'): 2,
+        ('a', 'reactor'): 3,
+        ('a', 'centrifuge'): 4,
+        ('b', 'mixer'): 4,
+        ('b', 'reactor'): 6,
+        ('b', 'centrifuge'): 3,
+    }
+
+    # Size factor for product i in stage j [kg/L]
+    m.s = pe.Param(
+        m.i, m.j, initialize=s_init, doc='Size factor for product i in stage j [kg/L]'
+    )
+
+    t_init = {
+        ('a', 'mixer'): 8,
+        ('a', 'reactor'): 20,
+        ('a', 'centrifuge'): 4,
+        ('b', 'mixer'): 10,
+        ('b', 'reactor'): 12,
+        ('b', 'centrifuge'): 3,
+    }
+
+    # Processing time of product i in batch j [hr]
+    m.t = pe.Param(
+        m.i, m.j, initialize=t_init, doc='Processing time of product i in batch j [hr]'
+    )
+
+    # Variables
+    m.Y = pe.BooleanVar(m.k, m.j, doc='Stage existence')  # Stage existence
+    m.coeffval = pe.Var(
+        m.k,
+        m.j,
+        within=pe.NonNegativeReals,
+        bounds=(0, pe.log(NK)),
+        doc='Activation of Coefficient',
+    )  # Activation of coeff
+    m.v = pe.Var(
+        m.j,
+        within=pe.NonNegativeReals,
+        bounds=(pe.log(m.vlow), pe.log(m.vupp)),
+        doc='Colume of stage j [L]',
+    )  # Volume of stage j [L]
+    m.b = pe.Var(
+        m.i, within=pe.NonNegativeReals, doc='Batch size of product i [L]'
+    )  # Batch size of product i [L]
+    m.tl = pe.Var(
+        m.i, within=pe.NonNegativeReals, doc='Cycle time of product i [hr]'
+    )  # Cycle time of product i [hr]
+    # Number of units in parallel stage j
+    m.n = pe.Var(
+        m.j, within=pe.NonNegativeReals, doc='Number of units in parallel stage j'
+    )
+
+    # Constraints
+
+    # Volume requirement in stage j
+    @m.Constraint(m.i, m.j)
+    def vol(m, i, j):
+        """
+        Volume Requirement for Stage j.
+        Equation:
+            v_j \geq log(s_ij) + b_i for i = a, b and j = mixer, reactor, centrifuge
+
+        Args:
+            m (pyomo.ConcreteModel): small batch GDP model
+            i (str): product
+            j (str): stage
+
+        Returns:
+            Algebraic Constraint
+        """
+        return m.v[j] >= pe.log(m.s[i, j]) + m.b[i]
+
+    # Cycle time for each product i
+    @m.Constraint(m.i, m.j)
+    def cycle(m, i, j):
+        """
+        Cycle time for each product i.
+        Equation:
+            n_j + tl_i \geq log(t_ij) for i = a, b and j = mixer, reactor, centrifuge
+
+        Args:
+            m (pyomo.ConcreteModel): small batch GDP model
+            i (str): product
+            j (str): stage
+
+        Returns:
+            Algebraic Constraint
+        """
+        return m.n[j] + m.tl[i] >= pe.log(m.t[i, j])
+
+    # Constraint for production time
+    @m.Constraint()
+    def time(m):
+        """
+        Production time constraint.
+        Equation:
+            \sum_{i \in I} q_i * \exp(tl_i - b_i) \leq h
+
+        Args:
+            m (pyomo.ConcreteModel): small batch GDP model
+
+        Returns:
+            Algebraic Constraint
+        """
+        return sum(m.q[i] * pe.exp(m.tl[i] - m.b[i]) for i in m.i) <= m.h
+
+    # Relating number of units to 0-1 variables
+    @m.Constraint(m.j)
+    def units(m, j):
+        """
+        Relating number of units to 0-1 variables.
+        Equation:
+            n_j = \sum_{k \in K} coeffval_{k,j} for j = mixer, reactor, centrifuge
+
+        Args:
+            m (pyomo.ConcreteModel): small batch GDP model
+            j (str): stage
+            k (int): number of parallel units
+
+        Returns:
+            Algebraic Constraint
+        """
+        return m.n[j] == sum(m.coeffval[k, j] for k in m.k)
+
+    # Only one choice for parallel units is feasible
+    @m.LogicalConstraint(m.j)
+    def lim(m, j):
+        """
+        Only one choice for parallel units is feasible.
+        Equation:
+            \sum_{k \in K} Y_{k,j} = 1 for j = mixer, reactor, centrifuge
+
+        Args:
+            m (pyomo.ConcreteModel): small batch GDP model
+            j (str): stage
+
+        Returns:
+            Logical Constraint
+        """
+        return pe.exactly(1, m.Y[1, j], m.Y[2, j], m.Y[3, j])
+
+    # _______ Disjunction_________
+
+    def build_existence_equations(disjunct, k, j):
+        """
+        Build the Logic Proposition (euqations) for the existence of the stage.
+
+        Args:
+            disjunct (pyomo.gdp.Disjunct): Disjunct block
+            k (int): number of parallel units
+            j (str): stage
+
+        Returns:
+            None, the proposition is built inside the function
+        """
+        m = disjunct.model()
+
+        # Coeffval activation
+        @disjunct.Constraint()
+        def coeffval_act(disjunct):
+            """
+            Coeffval activation.
+            m.coeffval[k,j] = m.coeff[k] = log(k)
+
+            Args:
+                disjunct (pyomo.gdp.Disjunct): Disjunct block
+
+            Returns:
+               Logical Constraint
+            """
+            return m.coeffval[k, j] == m.coeff[k]
+
+    def build_not_existence_equations(disjunct, k, j):
+        """
+        Build the Logic Proposition (euqations) for the unexistence of the stage.
+
+        Args:
+            disjunct (pyomo.gdp.Disjunct): Disjunct block
+            k (int): number of parallel units
+            j (str): stage
+
+        Returns:
+            None, the proposition is built inside the function.
+        """
+        m = disjunct.model()
+
+        # Coeffval deactivation
+        @disjunct.Constraint()
+        def coeffval_deact(disjunct):
+            """
+            Coeffval deactivation.
+            m.coeffval[k,j] = 0
+
+            Args:
+                disjunct (pyomo.gdp.Disjunct): Disjunct block
+
+            Returns:
+                Logical Constraint
+            """
+            return m.coeffval[k, j] == 0
+
+    # Create disjunction block
+    m.Y_exists = Disjunct(m.k, m.j, rule=build_existence_equations)
+    m.Y_not_exists = Disjunct(m.k, m.j, rule=build_not_existence_equations)
+
+    # Create disjunction
+
+    @m.Disjunction(m.k, m.j)
+    def Y_exists_or_not(m, k, j):
+        """
+        Build the Logical Disjunctions of the GDP model for the small batch problem.
+
+        Args:
+            m (pyomo.ConcreteModel): small batch GDP model
+            k (int): number of parallel units
+            j (str): stage
+
+        Returns:
+            Y_exists_or_not (list): List of disjuncts
+        """
+        return [m.Y_exists[k, j], m.Y_not_exists[k, j]]
+
+    # Associate Boolean variables with with disjunction
+    for k in m.k:
+        for j in m.j:
+            m.Y[k, j].associate_binary_var(m.Y_exists[k, j].indicator_var)
+
+    # ____________________________
+
+    # Objective
+    def obj_rule(m):
+        """
+        Objective: mininimize the investment cost [$].
+        Equation:
+            min z = sum(alpha[j] * exp(n[j] + beta[j]*v[j])) for j = mixer, reactor, centrifuge
+
+        Args:
+            m (pyomo.ConcreteModel): small batch GDP model
+
+        Returns:
+            Objective function (pyomo.Objective): Objective function to minimize the investment cost [$].
+        """
+        return sum(m.alpha[j] * (pe.exp(m.n[j] + m.beta[j] * m.v[j])) for j in m.j)
+
+    m.obj = pe.Objective(rule=obj_rule, sense=pe.minimize)
+
+    return m
+
+
+def external_ref(m, x, logic_expr=None):
+    """
+    Add the external variables to the GDP optimization problem.
+
+    Args:
+        m (pyomo.ConcreteModel): GDP optimization model
+        x (list): External variables
+        logic_expr (list, optional): Logic expressions to be used in the disjunctive constraints
+
+    Returns:
+        m (pyomo.ConcreteModel): GDP optimization model with the external variables
+    """
+    ext_var = {}
+    p = 0
+    for j in m.j:
+        ext_var[j] = x[p]
+        p = p + 1
+
+    for k in m.k:
+        for j in m.j:
+            if k == ext_var[j]:
+                m.Y[k, j].fix(True)
+                # m.Y_exists[k, j].indicator_var.fix(
+                #     True
+                # )  # Is this necessary?: m.Y_exists[k, j].indicator_var.fix(True).
+                # m.Y_not_exists[k, j].indicator_var.fix(
+                #     False
+                # )  # Is this necessary?: m.Y_not_exists[k, j].indicator_var.fix(True),
+            else:
+                m.Y[k, j].fix(False)
+                # m.Y_exists[k, j].indicator_var.fix(
+                #     False
+                # )  # Is this necessary?: m.Y_exists[k, j].indicator_var.fix(True),
+                # m.Y_not_exists[k, j].indicator_var.fix(
+                #     True
+                # )  # Is this necessary?: m.Y_not_exists[k, j].indicator_var.fix(True),
+
+    pe.TransformationFactory('core.logical_to_linear').apply_to(m)
+    pe.TransformationFactory('gdp.fix_disjuncts').apply_to(m)
+    pe.TransformationFactory('contrib.deactivate_trivial_constraints').apply_to(
+        m, tmp=False, ignore_infeasible=True
+    )
+
+    return m
+
+
+def solve_with_minlp(m, transformation='bigm', minlp='baron', timelimit=10):
+    """
+    Solve the GDP optimization problem with a MINLP solver.
+    The function applies the big-M Reformulation on the GDP and solve the MINLP problem with BARON.
+
+    Args:
+        m (pyomo.ConcreteModel): GDP optimization model
+        transformation (str, optional): Reformulation applied to the GDP.
+        minlp (str, optional): MINLP solver.
+        timelimit (float, optional): Time limit for the MINLP solver.
+
+    Returns:
+        m (pyomo.ConcreteModel): GDP optimization model with the solution.
+    """
+    # Transformation step
+    pe.TransformationFactory('core.logical_to_linear').apply_to(m)
+    transformation_string = 'gdp.' + transformation
+    pe.TransformationFactory(transformation_string).apply_to(m)
+
+    # Solution step
+    dir_path = os.path.dirname(os.path.abspath(__file__))
+    gams_path = os.path.join(dir_path, "gamsfiles/")
+    if not (os.path.exists(gams_path)):
+        print(
+            'Directory for automatically generated files '
+            + gams_path
+            + ' does not exist. We will create it'
+        )
+        os.makedirs(gams_path)
+
+    solvername = 'gams'
+    opt = SolverFactory(solvername, solver=minlp)
+    m.results = opt.solve(
+        m,
+        tee=True,
+        # Uncomment the following lines if you want to save GAMS models
+        # keepfiles=True,
+        # tmpdir=gams_path,
+        # symbolic_solver_labels=True,
+        add_options=[
+            'option reslim = ' + str(timelimit) + ';'
+            'option optcr = 0.0;'
+            # Uncomment the following lines to setup IIS computation of BARON through option file
+            # 'GAMS_MODEL.optfile = 1;'
+            # '\n'
+            # '$onecho > baron.opt \n'
+            # 'CompIIS 1 \n'
+            # '$offecho'
+            # 'display(execError);'
+        ],
+    )
+    update_boolean_vars_from_binary(m)
+    return m
+
+
+if __name__ == "__main__":
+    m = build_small_batch()
+    m_solved = solve_with_minlp(m, transformation='bigm', minlp='baron', timelimit=120)
+
+    # EXTERNAL REF TEST (this thest can be deleted)
+    newmodel = external_ref(m, [1, 2, 3], logic_expr=None)
+    # print('External Ref Test')
+    # print('Y[1, mixer] = ', newmodel.Y[1, 'mixer'].value)
+    # print('Y_exists[1, mixer] = ', newmodel.Y_exists[1, 'mixer'].indicator_var.value)
+    # print('Y_not_exists[1, mixer] = ', newmodel.Y_not_exists[1, 'mixer'].indicator_var.value)
+    # print('Y[2, mixer] = ', newmodel.Y[2, 'mixer'].value)
+    # print('Y_exists[2, mixer] = ', newmodel.Y_exists[2, 'mixer'].indicator_var.value)
+    # print('Y_not_exists[2, mixer] = ', newmodel.Y_not_exists[2, 'mixer'].indicator_var.value)

From 0c500440cac7e011fc303714152cf9ed8323180a Mon Sep 17 00:00:00 2001
From: Albert Lee <lee4382@purdue.edu>
Date: Tue, 21 May 2024 14:17:53 -0400
Subject: [PATCH 2/9] Managed the comments.

---
 gdplib/small_batch/gdp_small_batch.py | 85 +++------------------------
 1 file changed, 7 insertions(+), 78 deletions(-)

diff --git a/gdplib/small_batch/gdp_small_batch.py b/gdplib/small_batch/gdp_small_batch.py
index b9cdce5..314cdff 100644
--- a/gdplib/small_batch/gdp_small_batch.py
+++ b/gdplib/small_batch/gdp_small_batch.py
@@ -2,9 +2,11 @@
 gdp_small_batch.py
 The gdp_small_batch.py module contains the GDP model for the small batch problem based on the Kocis and Grossmann (1988) paper.
 The problem is based on the Example 4 of the paper.
+The objective is to minimize the investment cost of the batch units.
 
 References:
-    - Kocis, G. R.; Grossmann, I. E. Global Optimization of Nonconvex Mixed-Integer Nonlinear Programming (MINLP) Problems in Process Synthesis. Ind. Eng. Chem. Res. 1988, 27 (8), 1407–1421. 
+    [1] Kocis, G. R.; Grossmann, I. E. Global Optimization of Nonconvex Mixed-Integer Nonlinear Programming (MINLP) Problems in Process Synthesis. Ind. Eng. Chem. Res. 1988, 27 (8), 1407-1421. https://doi.org/10.1021/ie00080a013
+    [2] Ovalle, D., Liñán, D. A., Lee, A., Gómez, J. M., Ricardez-Sandoval, L., Grossmann, I. E., & Neira, D. E. B. (2024). Logic-Based Discrete-Steepest Descent: A Solution Method for Process Synthesis Generalized Disjunctive Programs. arXiv preprint arXiv:2405.05358. https://doi.org/10.48550/arXiv.2405.05358
 """
 import os
 
@@ -22,7 +24,8 @@ def build_small_batch():
     The function build the GDP model for the small batch problem.
 
     References:
-    - Kocis, G. R.; Grossmann, I. E. Global Optimization of Nonconvex Mixed-Integer Nonlinear Programming (MINLP) Problems in Process Synthesis. Ind. Eng. Chem. Res. 1988, 27 (8), 1407–1421.
+    [1] Kocis, G. R.; Grossmann, I. E. Global Optimization of Nonconvex Mixed-Integer Nonlinear Programming (MINLP) Problems in Process Synthesis. Ind. Eng. Chem. Res. 1988, 27 (8), 1407-1421. https://doi.org/10.1021/ie00080a013
+    [2] Ovalle, D., Liñán, D. A., Lee, A., Gómez, J. M., Ricardez-Sandoval, L., Grossmann, I. E., & Neira, D. E. B. (2024). Logic-Based Discrete-Steepest Descent: A Solution Method for Process Synthesis Generalized Disjunctive Programs. arXiv preprint arXiv:2405.05358. https://doi.org/10.48550/arXiv.2405.05358
 
     Args:
         None
@@ -366,20 +369,8 @@ def external_ref(m, x, logic_expr=None):
         for j in m.j:
             if k == ext_var[j]:
                 m.Y[k, j].fix(True)
-                # m.Y_exists[k, j].indicator_var.fix(
-                #     True
-                # )  # Is this necessary?: m.Y_exists[k, j].indicator_var.fix(True).
-                # m.Y_not_exists[k, j].indicator_var.fix(
-                #     False
-                # )  # Is this necessary?: m.Y_not_exists[k, j].indicator_var.fix(True),
             else:
                 m.Y[k, j].fix(False)
-                # m.Y_exists[k, j].indicator_var.fix(
-                #     False
-                # )  # Is this necessary?: m.Y_exists[k, j].indicator_var.fix(True),
-                # m.Y_not_exists[k, j].indicator_var.fix(
-                #     True
-                # )  # Is this necessary?: m.Y_not_exists[k, j].indicator_var.fix(True),
 
     pe.TransformationFactory('core.logical_to_linear').apply_to(m)
     pe.TransformationFactory('gdp.fix_disjuncts').apply_to(m)
@@ -390,71 +381,9 @@ def external_ref(m, x, logic_expr=None):
     return m
 
 
-def solve_with_minlp(m, transformation='bigm', minlp='baron', timelimit=10):
-    """
-    Solve the GDP optimization problem with a MINLP solver.
-    The function applies the big-M Reformulation on the GDP and solve the MINLP problem with BARON.
-
-    Args:
-        m (pyomo.ConcreteModel): GDP optimization model
-        transformation (str, optional): Reformulation applied to the GDP.
-        minlp (str, optional): MINLP solver.
-        timelimit (float, optional): Time limit for the MINLP solver.
-
-    Returns:
-        m (pyomo.ConcreteModel): GDP optimization model with the solution.
-    """
-    # Transformation step
-    pe.TransformationFactory('core.logical_to_linear').apply_to(m)
-    transformation_string = 'gdp.' + transformation
-    pe.TransformationFactory(transformation_string).apply_to(m)
-
-    # Solution step
-    dir_path = os.path.dirname(os.path.abspath(__file__))
-    gams_path = os.path.join(dir_path, "gamsfiles/")
-    if not (os.path.exists(gams_path)):
-        print(
-            'Directory for automatically generated files '
-            + gams_path
-            + ' does not exist. We will create it'
-        )
-        os.makedirs(gams_path)
-
-    solvername = 'gams'
-    opt = SolverFactory(solvername, solver=minlp)
-    m.results = opt.solve(
-        m,
-        tee=True,
-        # Uncomment the following lines if you want to save GAMS models
-        # keepfiles=True,
-        # tmpdir=gams_path,
-        # symbolic_solver_labels=True,
-        add_options=[
-            'option reslim = ' + str(timelimit) + ';'
-            'option optcr = 0.0;'
-            # Uncomment the following lines to setup IIS computation of BARON through option file
-            # 'GAMS_MODEL.optfile = 1;'
-            # '\n'
-            # '$onecho > baron.opt \n'
-            # 'CompIIS 1 \n'
-            # '$offecho'
-            # 'display(execError);'
-        ],
-    )
-    update_boolean_vars_from_binary(m)
-    return m
 
 
 if __name__ == "__main__":
     m = build_small_batch()
-    m_solved = solve_with_minlp(m, transformation='bigm', minlp='baron', timelimit=120)
-
-    # EXTERNAL REF TEST (this thest can be deleted)
-    newmodel = external_ref(m, [1, 2, 3], logic_expr=None)
-    # print('External Ref Test')
-    # print('Y[1, mixer] = ', newmodel.Y[1, 'mixer'].value)
-    # print('Y_exists[1, mixer] = ', newmodel.Y_exists[1, 'mixer'].indicator_var.value)
-    # print('Y_not_exists[1, mixer] = ', newmodel.Y_not_exists[1, 'mixer'].indicator_var.value)
-    # print('Y[2, mixer] = ', newmodel.Y[2, 'mixer'].value)
-    # print('Y_exists[2, mixer] = ', newmodel.Y_exists[2, 'mixer'].indicator_var.value)
-    # print('Y_not_exists[2, mixer] = ', newmodel.Y_not_exists[2, 'mixer'].indicator_var.value)
+
+

From 80903ee7ffab0c49fe6b968e103d5786c22120ed Mon Sep 17 00:00:00 2001
From: Albert Lee <lee4382@purdue.edu>
Date: Tue, 21 May 2024 15:57:33 -0400
Subject: [PATCH 3/9] Numpy Documentation Style

---
 gdplib/small_batch/gdp_small_batch.py | 324 ++++++++++++++------------
 1 file changed, 175 insertions(+), 149 deletions(-)

diff --git a/gdplib/small_batch/gdp_small_batch.py b/gdplib/small_batch/gdp_small_batch.py
index 314cdff..3961bdd 100644
--- a/gdplib/small_batch/gdp_small_batch.py
+++ b/gdplib/small_batch/gdp_small_batch.py
@@ -10,7 +10,7 @@
 """
 import os
 
-import pyomo.environ as pe
+import pyomo.environ as pyo
 from pyomo.core.base.misc import display
 from pyomo.core.plugins.transform.logical_to_linear import (
     update_boolean_vars_from_binary,
@@ -21,58 +21,58 @@
 
 def build_small_batch():
     """
-    The function build the GDP model for the small batch problem.
+    Build the GDP model for the small batch problem.
 
-    References:
-    [1] Kocis, G. R.; Grossmann, I. E. Global Optimization of Nonconvex Mixed-Integer Nonlinear Programming (MINLP) Problems in Process Synthesis. Ind. Eng. Chem. Res. 1988, 27 (8), 1407-1421. https://doi.org/10.1021/ie00080a013
-    [2] Ovalle, D., Liñán, D. A., Lee, A., Gómez, J. M., Ricardez-Sandoval, L., Grossmann, I. E., & Neira, D. E. B. (2024). Logic-Based Discrete-Steepest Descent: A Solution Method for Process Synthesis Generalized Disjunctive Programs. arXiv preprint arXiv:2405.05358. https://doi.org/10.48550/arXiv.2405.05358
+    Returns
+    -------
+    m : Pyomo.ConcreteModel
+        The GDP model for the small batch problem is created.
 
-    Args:
-        None
-
-    Returns:
-        m (pyomo.ConcreteModel): The GDP model for the small batch problem is created.
+    References
+    ----------
+    [1] Kocis, G. R.; Grossmann, I. E. (1988). Global Optimization of Nonconvex Mixed-Integer Nonlinear Programming (MINLP) Problems in Process Synthesis. Ind. Eng. Chem. Res., 27(8), 1407-1421. https://doi.org/10.1021/ie00080a013
+    [2] Ovalle, D., Liñán, D. A., Lee, A., Gómez, J. M., Ricardez-Sandoval, L., Grossmann, I. E., & Neira, D. E. B. (2024). Logic-Based Discrete-Steepest Descent: A Solution Method for Process Synthesis Generalized Disjunctive Programs. arXiv preprint arXiv:2405.05358. https://doi.org/10.48550/arXiv.2405.05358
     """
     NK = 3
 
     # Model
-    m = pe.ConcreteModel()
+    m = pyo.ConcreteModel()
 
     # Sets
-    m.i = pe.Set(
+    m.i = pyo.Set(
         initialize=['a', 'b'], doc='Set of products'
     )  # Set of products, i = a, b
-    m.j = pe.Set(
+    m.j = pyo.Set(
         initialize=['mixer', 'reactor', 'centrifuge']
     )  # Set of stages, j = mixer, reactor, centrifuge
-    m.k = pe.RangeSet(NK)  # Set of potential number of parallel units, k = 1, 2, 3
+    m.k = pyo.RangeSet(NK, doc="Set of potential number of parallel units")  # Set of potential number of parallel units, k = 1, 2, 3
 
     # Parameters and Scalars
 
-    m.h = pe.Param(
+    m.h = pyo.Param(
         initialize=6000, doc='Horizon time [hr]'
     )  # Horizon time  (available time) [hr]
-    m.vlow = pe.Param(
+    m.vlow = pyo.Param(
         initialize=250, doc='Lower bound for size of batch unit [L]'
     )  # Lower bound for size of batch unit [L]
-    m.vupp = pe.Param(
+    m.vupp = pyo.Param(
         initialize=2500, doc='Upper bound for size of batch unit [L]'
     )  # Upper bound for size of batch unit [L]
 
     # Demand of product i
-    m.q = pe.Param(
+    m.q = pyo.Param(
         m.i,
         initialize={'a': 200000, 'b': 150000},
         doc='Production rate of the product [kg]',
     )
     # Cost coefficient for batch units
-    m.alpha = pe.Param(
+    m.alpha = pyo.Param(
         m.j,
         initialize={'mixer': 250, 'reactor': 500, 'centrifuge': 340},
         doc='Cost coefficient for batch units [$/L^beta*No. of units]]',
     )
     # Cost exponent for batch units
-    m.beta = pe.Param(
+    m.beta = pyo.Param(
         m.j,
         initialize={'mixer': 0.6, 'reactor': 0.6, 'centrifuge': 0.6},
         doc='Cost exponent for batch units',
@@ -82,17 +82,22 @@ def coeff_init(m, k):
         """
         Coefficient for number of parallel units.
 
-        Args:
-            m (pyomo.ConcreteModel): small batch GDP model
-            k (int): number of parallel units
-
-        Returns:
-            Coefficient for number of parallel units.
+        Parameters
+        ----------
+        m : Pyomo.ConcreteModel
+            The small batch GDP model.
+        k : int
+            The number of parallel units.
+
+        Returns
+        -------
+        pyomo.Param
+            Coefficient for number of parallel units. logarithm of k is applied to convexify the model.
         """
-        return pe.log(k)
+        return pyo.log(k)
 
     # Represent number of parallel units
-    m.coeff = pe.Param(
+    m.coeff = pyo.Param(
         m.k, initialize=coeff_init, doc='Coefficient for number of parallel units'
     )
 
@@ -106,7 +111,7 @@ def coeff_init(m, k):
     }
 
     # Size factor for product i in stage j [kg/L]
-    m.s = pe.Param(
+    m.s = pyo.Param(
         m.i, m.j, initialize=s_init, doc='Size factor for product i in stage j [kg/L]'
     )
 
@@ -120,34 +125,34 @@ def coeff_init(m, k):
     }
 
     # Processing time of product i in batch j [hr]
-    m.t = pe.Param(
+    m.t = pyo.Param(
         m.i, m.j, initialize=t_init, doc='Processing time of product i in batch j [hr]'
     )
 
     # Variables
-    m.Y = pe.BooleanVar(m.k, m.j, doc='Stage existence')  # Stage existence
-    m.coeffval = pe.Var(
+    m.Y = pyo.BooleanVar(m.k, m.j, doc='Stage existence')  # Stage existence
+    m.coeffval = pyo.Var(
         m.k,
         m.j,
-        within=pe.NonNegativeReals,
-        bounds=(0, pe.log(NK)),
+        within=pyo.NonNegativeReals,
+        bounds=(0, pyo.log(NK)),
         doc='Activation of Coefficient',
     )  # Activation of coeff
-    m.v = pe.Var(
+    m.v = pyo.Var(
         m.j,
-        within=pe.NonNegativeReals,
-        bounds=(pe.log(m.vlow), pe.log(m.vupp)),
+        within=pyo.NonNegativeReals,
+        bounds=(pyo.log(m.vlow), pyo.log(m.vupp)),
         doc='Colume of stage j [L]',
     )  # Volume of stage j [L]
-    m.b = pe.Var(
-        m.i, within=pe.NonNegativeReals, doc='Batch size of product i [L]'
+    m.b = pyo.Var(
+        m.i, within=pyo.NonNegativeReals, doc='Batch size of product i [L]'
     )  # Batch size of product i [L]
-    m.tl = pe.Var(
-        m.i, within=pe.NonNegativeReals, doc='Cycle time of product i [hr]'
+    m.tl = pyo.Var(
+        m.i, within=pyo.NonNegativeReals, doc='Cycle time of product i [hr]'
     )  # Cycle time of product i [hr]
     # Number of units in parallel stage j
-    m.n = pe.Var(
-        m.j, within=pe.NonNegativeReals, doc='Number of units in parallel stage j'
+    m.n = pyo.Var(
+        m.j, within=pyo.NonNegativeReals, doc='Number of units in parallel stage j'
     )
 
     # Constraints
@@ -160,15 +165,21 @@ def vol(m, i, j):
         Equation:
             v_j \geq log(s_ij) + b_i for i = a, b and j = mixer, reactor, centrifuge
 
-        Args:
-            m (pyomo.ConcreteModel): small batch GDP model
-            i (str): product
-            j (str): stage
-
-        Returns:
+        Parameters
+        ----------
+        m : pyomo.ConcreteModel
+            The small batch GDP model.
+        i : str
+            Index of Product.
+        j : str
+            Stage.
+
+        Returns
+        -------
+        Pyomo.Constraint
             Algebraic Constraint
         """
-        return m.v[j] >= pe.log(m.s[i, j]) + m.b[i]
+        return m.v[j] >= pyo.log(m.s[i, j]) + m.b[i]
 
     # Cycle time for each product i
     @m.Constraint(m.i, m.j)
@@ -178,15 +189,21 @@ def cycle(m, i, j):
         Equation:
             n_j + tl_i \geq log(t_ij) for i = a, b and j = mixer, reactor, centrifuge
 
-        Args:
-            m (pyomo.ConcreteModel): small batch GDP model
-            i (str): product
-            j (str): stage
-
-        Returns:
+        Parameters
+        ----------
+        m : pyomo.ConcreteModel
+            The small batch GDP model.
+        i : str
+            Index of Product.
+        j : str
+            Index of Stage.
+
+        Returns
+        -------
+        Pyomo.Constraint
             Algebraic Constraint
         """
-        return m.n[j] + m.tl[i] >= pe.log(m.t[i, j])
+        return m.n[j] + m.tl[i] >= pyo.log(m.t[i, j])
 
     # Constraint for production time
     @m.Constraint()
@@ -196,13 +213,17 @@ def time(m):
         Equation:
             \sum_{i \in I} q_i * \exp(tl_i - b_i) \leq h
 
-        Args:
-            m (pyomo.ConcreteModel): small batch GDP model
+        Parameters
+        ----------
+        m : pyomo.ConcreteModel
+            The small batch GDP model.
 
-        Returns:
+        Returns
+        -------
+        Pyomo.Constraint
             Algebraic Constraint
         """
-        return sum(m.q[i] * pe.exp(m.tl[i] - m.b[i]) for i in m.i) <= m.h
+        return sum(m.q[i] * pyo.exp(m.tl[i] - m.b[i]) for i in m.i) <= m.h
 
     # Relating number of units to 0-1 variables
     @m.Constraint(m.j)
@@ -212,12 +233,16 @@ def units(m, j):
         Equation:
             n_j = \sum_{k \in K} coeffval_{k,j} for j = mixer, reactor, centrifuge
 
-        Args:
-            m (pyomo.ConcreteModel): small batch GDP model
-            j (str): stage
-            k (int): number of parallel units
+        Parameters
+        ----------
+        m : pyomo.ConcreteModel
+            The small batch GDP model.
+        j : str
+            Index of Stage.
 
-        Returns:
+        Returns
+        -------
+        Pyomo.Constraint
             Algebraic Constraint
         """
         return m.n[j] == sum(m.coeffval[k, j] for k in m.k)
@@ -230,28 +255,39 @@ def lim(m, j):
         Equation:
             \sum_{k \in K} Y_{k,j} = 1 for j = mixer, reactor, centrifuge
 
-        Args:
-            m (pyomo.ConcreteModel): small batch GDP model
-            j (str): stage
+        Parameters
+        ----------
+        m : pyomo.ConcreteModel
+            The small batch GDP model.
+        j : str
+            Index of Stage.
 
-        Returns:
-            Logical Constraint
+        Returns
+        -------
+        Pyomo.LogicalConstraint
+            Algebraic Constraint
         """
-        return pe.exactly(1, m.Y[1, j], m.Y[2, j], m.Y[3, j])
+        return pyo.exactly(1, m.Y[1, j], m.Y[2, j], m.Y[3, j])
 
     # _______ Disjunction_________
 
     def build_existence_equations(disjunct, k, j):
         """
-        Build the Logic Proposition (euqations) for the existence of the stage.
-
-        Args:
-            disjunct (pyomo.gdp.Disjunct): Disjunct block
-            k (int): number of parallel units
-            j (str): stage
-
-        Returns:
-            None, the proposition is built inside the function
+        Build the Logic  Disjunct Constraints (equation) for the existence of the stage.
+
+        Parameters
+        ----------
+        disjunct : Pyomo.Disjunct
+            Disjunct block for the existence of the stage.
+        k : int
+            Number of parallel units.
+        j : str
+            Index of Stage.
+
+        Returns
+        -------
+        None
+            None, the constraints are built inside the disjunct.
         """
         m = disjunct.model()
 
@@ -262,25 +298,35 @@ def coeffval_act(disjunct):
             Coeffval activation.
             m.coeffval[k,j] = m.coeff[k] = log(k)
 
-            Args:
-                disjunct (pyomo.gdp.Disjunct): Disjunct block
+            Parameters
+            ----------
+            disjunct : Pyomo.Disjunct
+                Disjunct block for the existence of the stage.
 
-            Returns:
-               Logical Constraint
+            Returns
+            -------
+            Pyomo.Constraint
+                A algebraic constraint
             """
             return m.coeffval[k, j] == m.coeff[k]
 
     def build_not_existence_equations(disjunct, k, j):
         """
-        Build the Logic Proposition (euqations) for the unexistence of the stage.
-
-        Args:
-            disjunct (pyomo.gdp.Disjunct): Disjunct block
-            k (int): number of parallel units
-            j (str): stage
-
-        Returns:
-            None, the proposition is built inside the function.
+        Build the Logic Disjunct Constraints (equations) for the absence of the stage.
+
+        Parameters
+        ----------
+        disjunct : Pyomo.Disjunct
+            Disjunct block for the absence of the stage.
+        k : int
+            Number of parallel units.
+        j : str
+            Index of Stage.
+
+        Returns
+        -------
+        None
+            None, the constraints are built inside the disjunct..
         """
         m = disjunct.model()
 
@@ -291,8 +337,10 @@ def coeffval_deact(disjunct):
             Coeffval deactivation.
             m.coeffval[k,j] = 0
 
-            Args:
-                disjunct (pyomo.gdp.Disjunct): Disjunct block
+            Parameters
+            ----------
+            disjunct : Pyomo.Disjunct
+                Disjunct block for the absence of the stage.
 
             Returns:
                 Logical Constraint
@@ -300,8 +348,8 @@ def coeffval_deact(disjunct):
             return m.coeffval[k, j] == 0
 
     # Create disjunction block
-    m.Y_exists = Disjunct(m.k, m.j, rule=build_existence_equations)
-    m.Y_not_exists = Disjunct(m.k, m.j, rule=build_not_existence_equations)
+    m.Y_exists = Disjunct(m.k, m.j, rule=build_existence_equations, doc="Existence of the stage")
+    m.Y_not_exists = Disjunct(m.k, m.j, rule=build_not_existence_equations, doc="Absence of the stage")
 
     # Create disjunction
 
@@ -310,13 +358,19 @@ def Y_exists_or_not(m, k, j):
         """
         Build the Logical Disjunctions of the GDP model for the small batch problem.
 
-        Args:
-            m (pyomo.ConcreteModel): small batch GDP model
-            k (int): number of parallel units
-            j (str): stage
-
-        Returns:
-            Y_exists_or_not (list): List of disjuncts
+        Parameters
+        ----------
+        m : pyomo.ConcreteModel
+            The small batch GDP model.
+        k : int
+            Number of parallel units.
+        j : str
+            Index of Stage.
+
+        Returns
+        -------
+        list
+            List of disjuncts. The disjunction is between the existence and absence of the stage.
         """
         return [m.Y_exists[k, j], m.Y_not_exists[k, j]]
 
@@ -330,60 +384,32 @@ def Y_exists_or_not(m, k, j):
     # Objective
     def obj_rule(m):
         """
-        Objective: mininimize the investment cost [$].
+        Objective: minimize the investment cost [$].
+
         Equation:
             min z = sum(alpha[j] * exp(n[j] + beta[j]*v[j])) for j = mixer, reactor, centrifuge
 
-        Args:
-            m (pyomo.ConcreteModel): small batch GDP model
+        Parameters
+        ----------
+        m : pyomo.ConcreteModel
+            The small batch GDP model.
 
-        Returns:
-            Objective function (pyomo.Objective): Objective function to minimize the investment cost [$].
+        Returns
+        -------
+        Pyomo.Objective
+            Objective function to minimize the investment cost [$].
         """
-        return sum(m.alpha[j] * (pe.exp(m.n[j] + m.beta[j] * m.v[j])) for j in m.j)
-
-    m.obj = pe.Objective(rule=obj_rule, sense=pe.minimize)
-
-    return m
-
-
-def external_ref(m, x, logic_expr=None):
-    """
-    Add the external variables to the GDP optimization problem.
-
-    Args:
-        m (pyomo.ConcreteModel): GDP optimization model
-        x (list): External variables
-        logic_expr (list, optional): Logic expressions to be used in the disjunctive constraints
+        return sum(m.alpha[j] * (pyo.exp(m.n[j] + m.beta[j] * m.v[j])) for j in m.j)
 
-    Returns:
-        m (pyomo.ConcreteModel): GDP optimization model with the external variables
-    """
-    ext_var = {}
-    p = 0
-    for j in m.j:
-        ext_var[j] = x[p]
-        p = p + 1
-
-    for k in m.k:
-        for j in m.j:
-            if k == ext_var[j]:
-                m.Y[k, j].fix(True)
-            else:
-                m.Y[k, j].fix(False)
-
-    pe.TransformationFactory('core.logical_to_linear').apply_to(m)
-    pe.TransformationFactory('gdp.fix_disjuncts').apply_to(m)
-    pe.TransformationFactory('contrib.deactivate_trivial_constraints').apply_to(
-        m, tmp=False, ignore_infeasible=True
-    )
+    m.obj = pyo.Objective(rule=obj_rule, sense=pyo.minimize)
 
     return m
 
 
-
-
 if __name__ == "__main__":
     m = build_small_batch()
-
+    pyo.TransformationFactory('core.logical_to_linear').apply_to(m)
+    pyo.TransformationFactory('gdp.bigm').apply_to(m)
+    pyo.SolverFactory('gams').solve(m, solver='baron', tee=True, add_options=['option optcr=1e-6;'])
+    display(m)
 

From 0b75697193cca20d961322358f52cecd87f97bdf Mon Sep 17 00:00:00 2001
From: Albert Lee <lee4382@purdue.edu>
Date: Tue, 21 May 2024 16:19:55 -0400
Subject: [PATCH 4/9] renew the documentation

---
 gdplib/small_batch/gdp_small_batch.py | 54 ++++++++++++++++-----------
 1 file changed, 33 insertions(+), 21 deletions(-)

diff --git a/gdplib/small_batch/gdp_small_batch.py b/gdplib/small_batch/gdp_small_batch.py
index 3961bdd..853bbcc 100644
--- a/gdplib/small_batch/gdp_small_batch.py
+++ b/gdplib/small_batch/gdp_small_batch.py
@@ -162,8 +162,9 @@ def coeff_init(m, k):
     def vol(m, i, j):
         """
         Volume Requirement for Stage j.
-        Equation:
-            v_j \geq log(s_ij) + b_i for i = a, b and j = mixer, reactor, centrifuge
+        Equation
+        --------
+        v_j \geq log(s_ij) + b_i for i = a, b and j = mixer, reactor, centrifuge
 
         Parameters
         ----------
@@ -177,7 +178,7 @@ def vol(m, i, j):
         Returns
         -------
         Pyomo.Constraint
-            Algebraic Constraint
+            A Pyomo constraint object representing the volume requirement for a given stage.
         """
         return m.v[j] >= pyo.log(m.s[i, j]) + m.b[i]
 
@@ -186,8 +187,10 @@ def vol(m, i, j):
     def cycle(m, i, j):
         """
         Cycle time for each product i.
-        Equation:
-            n_j + tl_i \geq log(t_ij) for i = a, b and j = mixer, reactor, centrifuge
+
+        Equation
+        --------
+        n_j + tl_i \geq log(t_ij) for i = a, b and j = mixer, reactor, centrifuge
 
         Parameters
         ----------
@@ -201,7 +204,7 @@ def cycle(m, i, j):
         Returns
         -------
         Pyomo.Constraint
-            Algebraic Constraint
+            A Pyomo constraint object representing the cycle time requirement for each product in each stage.
         """
         return m.n[j] + m.tl[i] >= pyo.log(m.t[i, j])
 
@@ -211,7 +214,7 @@ def time(m):
         """
         Production time constraint.
         Equation:
-            \sum_{i \in I} q_i * \exp(tl_i - b_i) \leq h
+            sum_{i \in I} q_i * \exp(tl_i - b_i) \leq h
 
         Parameters
         ----------
@@ -221,7 +224,7 @@ def time(m):
         Returns
         -------
         Pyomo.Constraint
-            Algebraic Constraint
+            A Pyomo constraint object representing the production time constraint.
         """
         return sum(m.q[i] * pyo.exp(m.tl[i] - m.b[i]) for i in m.i) <= m.h
 
@@ -231,7 +234,7 @@ def units(m, j):
         """
         Relating number of units to 0-1 variables.
         Equation:
-            n_j = \sum_{k \in K} coeffval_{k,j} for j = mixer, reactor, centrifuge
+            n_j = sum_{k \in K} coeffval_{k,j} for j = mixer, reactor, centrifuge
 
         Parameters
         ----------
@@ -243,7 +246,7 @@ def units(m, j):
         Returns
         -------
         Pyomo.Constraint
-            Algebraic Constraint
+            A Pyomo constraint object representing the relationship between the number of units and the binary variables.
         """
         return m.n[j] == sum(m.coeffval[k, j] for k in m.k)
 
@@ -253,7 +256,7 @@ def lim(m, j):
         """
         Only one choice for parallel units is feasible.
         Equation:
-            \sum_{k \in K} Y_{k,j} = 1 for j = mixer, reactor, centrifuge
+            sum_{k \in K} Y_{k,j} = 1 for j = mixer, reactor, centrifuge
 
         Parameters
         ----------
@@ -265,7 +268,7 @@ def lim(m, j):
         Returns
         -------
         Pyomo.LogicalConstraint
-            Algebraic Constraint
+            A Pyomo logical constraint ensuring only one choice for parallel units is feasible.
         """
         return pyo.exactly(1, m.Y[1, j], m.Y[2, j], m.Y[3, j])
 
@@ -291,11 +294,14 @@ def build_existence_equations(disjunct, k, j):
         """
         m = disjunct.model()
 
-        # Coeffval activation
+        # Coefficient value activation
         @disjunct.Constraint()
         def coeffval_act(disjunct):
             """
-            Coeffval activation.
+            Coefficien value activation.
+
+            Equation
+            --------
             m.coeffval[k,j] = m.coeff[k] = log(k)
 
             Parameters
@@ -306,7 +312,7 @@ def coeffval_act(disjunct):
             Returns
             -------
             Pyomo.Constraint
-                A algebraic constraint
+                A Pyomo constraint object representing the activation of the coefficient value.
             """
             return m.coeffval[k, j] == m.coeff[k]
 
@@ -330,11 +336,14 @@ def build_not_existence_equations(disjunct, k, j):
         """
         m = disjunct.model()
 
-        # Coeffval deactivation
+        # Coefficient value deactivation
         @disjunct.Constraint()
         def coeffval_deact(disjunct):
             """
-            Coeffval deactivation.
+            Coefficient value deactivation.
+
+            Equation
+            --------
             m.coeffval[k,j] = 0
 
             Parameters
@@ -342,8 +351,10 @@ def coeffval_deact(disjunct):
             disjunct : Pyomo.Disjunct
                 Disjunct block for the absence of the stage.
 
-            Returns:
-                Logical Constraint
+            Returns
+            -------
+            Pyomo.Constraint
+                A Pyomo constraint object representing the deactivation of the coefficient value.
             """
             return m.coeffval[k, j] == 0
 
@@ -386,8 +397,9 @@ def obj_rule(m):
         """
         Objective: minimize the investment cost [$].
 
-        Equation:
-            min z = sum(alpha[j] * exp(n[j] + beta[j]*v[j])) for j = mixer, reactor, centrifuge
+        Equation
+        --------
+        min z = sum(alpha[j] * exp(n[j] + beta[j]*v[j])) for j = mixer, reactor, centrifuge
 
         Parameters
         ----------

From 87186cbd702692ac7c4996570ab68f120af28946 Mon Sep 17 00:00:00 2001
From: Albert Lee <lee4382@purdue.edu>
Date: Tue, 21 May 2024 16:22:44 -0400
Subject: [PATCH 5/9] black format.

---
 gdplib/small_batch/gdp_small_batch.py | 85 +++++++++++++++------------
 1 file changed, 46 insertions(+), 39 deletions(-)

diff --git a/gdplib/small_batch/gdp_small_batch.py b/gdplib/small_batch/gdp_small_batch.py
index 853bbcc..ff910d8 100644
--- a/gdplib/small_batch/gdp_small_batch.py
+++ b/gdplib/small_batch/gdp_small_batch.py
@@ -40,42 +40,44 @@ def build_small_batch():
 
     # Sets
     m.i = pyo.Set(
-        initialize=['a', 'b'], doc='Set of products'
+        initialize=["a", "b"], doc="Set of products"
     )  # Set of products, i = a, b
     m.j = pyo.Set(
-        initialize=['mixer', 'reactor', 'centrifuge']
+        initialize=["mixer", "reactor", "centrifuge"]
     )  # Set of stages, j = mixer, reactor, centrifuge
-    m.k = pyo.RangeSet(NK, doc="Set of potential number of parallel units")  # Set of potential number of parallel units, k = 1, 2, 3
+    m.k = pyo.RangeSet(
+        NK, doc="Set of potential number of parallel units"
+    )  # Set of potential number of parallel units, k = 1, 2, 3
 
     # Parameters and Scalars
 
     m.h = pyo.Param(
-        initialize=6000, doc='Horizon time [hr]'
+        initialize=6000, doc="Horizon time [hr]"
     )  # Horizon time  (available time) [hr]
     m.vlow = pyo.Param(
-        initialize=250, doc='Lower bound for size of batch unit [L]'
+        initialize=250, doc="Lower bound for size of batch unit [L]"
     )  # Lower bound for size of batch unit [L]
     m.vupp = pyo.Param(
-        initialize=2500, doc='Upper bound for size of batch unit [L]'
+        initialize=2500, doc="Upper bound for size of batch unit [L]"
     )  # Upper bound for size of batch unit [L]
 
     # Demand of product i
     m.q = pyo.Param(
         m.i,
-        initialize={'a': 200000, 'b': 150000},
-        doc='Production rate of the product [kg]',
+        initialize={"a": 200000, "b": 150000},
+        doc="Production rate of the product [kg]",
     )
     # Cost coefficient for batch units
     m.alpha = pyo.Param(
         m.j,
-        initialize={'mixer': 250, 'reactor': 500, 'centrifuge': 340},
-        doc='Cost coefficient for batch units [$/L^beta*No. of units]]',
+        initialize={"mixer": 250, "reactor": 500, "centrifuge": 340},
+        doc="Cost coefficient for batch units [$/L^beta*No. of units]]",
     )
     # Cost exponent for batch units
     m.beta = pyo.Param(
         m.j,
-        initialize={'mixer': 0.6, 'reactor': 0.6, 'centrifuge': 0.6},
-        doc='Cost exponent for batch units',
+        initialize={"mixer": 0.6, "reactor": 0.6, "centrifuge": 0.6},
+        doc="Cost exponent for batch units",
     )
 
     def coeff_init(m, k):
@@ -98,61 +100,61 @@ def coeff_init(m, k):
 
     # Represent number of parallel units
     m.coeff = pyo.Param(
-        m.k, initialize=coeff_init, doc='Coefficient for number of parallel units'
+        m.k, initialize=coeff_init, doc="Coefficient for number of parallel units"
     )
 
     s_init = {
-        ('a', 'mixer'): 2,
-        ('a', 'reactor'): 3,
-        ('a', 'centrifuge'): 4,
-        ('b', 'mixer'): 4,
-        ('b', 'reactor'): 6,
-        ('b', 'centrifuge'): 3,
+        ("a", "mixer"): 2,
+        ("a", "reactor"): 3,
+        ("a", "centrifuge"): 4,
+        ("b", "mixer"): 4,
+        ("b", "reactor"): 6,
+        ("b", "centrifuge"): 3,
     }
 
     # Size factor for product i in stage j [kg/L]
     m.s = pyo.Param(
-        m.i, m.j, initialize=s_init, doc='Size factor for product i in stage j [kg/L]'
+        m.i, m.j, initialize=s_init, doc="Size factor for product i in stage j [kg/L]"
     )
 
     t_init = {
-        ('a', 'mixer'): 8,
-        ('a', 'reactor'): 20,
-        ('a', 'centrifuge'): 4,
-        ('b', 'mixer'): 10,
-        ('b', 'reactor'): 12,
-        ('b', 'centrifuge'): 3,
+        ("a", "mixer"): 8,
+        ("a", "reactor"): 20,
+        ("a", "centrifuge"): 4,
+        ("b", "mixer"): 10,
+        ("b", "reactor"): 12,
+        ("b", "centrifuge"): 3,
     }
 
     # Processing time of product i in batch j [hr]
     m.t = pyo.Param(
-        m.i, m.j, initialize=t_init, doc='Processing time of product i in batch j [hr]'
+        m.i, m.j, initialize=t_init, doc="Processing time of product i in batch j [hr]"
     )
 
     # Variables
-    m.Y = pyo.BooleanVar(m.k, m.j, doc='Stage existence')  # Stage existence
+    m.Y = pyo.BooleanVar(m.k, m.j, doc="Stage existence")  # Stage existence
     m.coeffval = pyo.Var(
         m.k,
         m.j,
         within=pyo.NonNegativeReals,
         bounds=(0, pyo.log(NK)),
-        doc='Activation of Coefficient',
+        doc="Activation of Coefficient",
     )  # Activation of coeff
     m.v = pyo.Var(
         m.j,
         within=pyo.NonNegativeReals,
         bounds=(pyo.log(m.vlow), pyo.log(m.vupp)),
-        doc='Colume of stage j [L]',
+        doc="Colume of stage j [L]",
     )  # Volume of stage j [L]
     m.b = pyo.Var(
-        m.i, within=pyo.NonNegativeReals, doc='Batch size of product i [L]'
+        m.i, within=pyo.NonNegativeReals, doc="Batch size of product i [L]"
     )  # Batch size of product i [L]
     m.tl = pyo.Var(
-        m.i, within=pyo.NonNegativeReals, doc='Cycle time of product i [hr]'
+        m.i, within=pyo.NonNegativeReals, doc="Cycle time of product i [hr]"
     )  # Cycle time of product i [hr]
     # Number of units in parallel stage j
     m.n = pyo.Var(
-        m.j, within=pyo.NonNegativeReals, doc='Number of units in parallel stage j'
+        m.j, within=pyo.NonNegativeReals, doc="Number of units in parallel stage j"
     )
 
     # Constraints
@@ -359,8 +361,12 @@ def coeffval_deact(disjunct):
             return m.coeffval[k, j] == 0
 
     # Create disjunction block
-    m.Y_exists = Disjunct(m.k, m.j, rule=build_existence_equations, doc="Existence of the stage")
-    m.Y_not_exists = Disjunct(m.k, m.j, rule=build_not_existence_equations, doc="Absence of the stage")
+    m.Y_exists = Disjunct(
+        m.k, m.j, rule=build_existence_equations, doc="Existence of the stage"
+    )
+    m.Y_not_exists = Disjunct(
+        m.k, m.j, rule=build_not_existence_equations, doc="Absence of the stage"
+    )
 
     # Create disjunction
 
@@ -420,8 +426,9 @@ def obj_rule(m):
 
 if __name__ == "__main__":
     m = build_small_batch()
-    pyo.TransformationFactory('core.logical_to_linear').apply_to(m)
-    pyo.TransformationFactory('gdp.bigm').apply_to(m)
-    pyo.SolverFactory('gams').solve(m, solver='baron', tee=True, add_options=['option optcr=1e-6;'])
+    pyo.TransformationFactory("core.logical_to_linear").apply_to(m)
+    pyo.TransformationFactory("gdp.bigm").apply_to(m)
+    pyo.SolverFactory("gams").solve(
+        m, solver="baron", tee=True, add_options=["option optcr=1e-6;"]
+    )
     display(m)
-

From 27063a3f9eb3b0c375a50c11133522ed3fea57ea Mon Sep 17 00:00:00 2001
From: Albert Lee <lee4382@purdue.edu>
Date: Tue, 21 May 2024 16:26:17 -0400
Subject: [PATCH 6/9] fixed documentation

---
 gdplib/small_batch/gdp_small_batch.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/gdplib/small_batch/gdp_small_batch.py b/gdplib/small_batch/gdp_small_batch.py
index ff910d8..1922aed 100644
--- a/gdplib/small_batch/gdp_small_batch.py
+++ b/gdplib/small_batch/gdp_small_batch.py
@@ -4,12 +4,12 @@
 The problem is based on the Example 4 of the paper.
 The objective is to minimize the investment cost of the batch units.
 
-References:
-    [1] Kocis, G. R.; Grossmann, I. E. Global Optimization of Nonconvex Mixed-Integer Nonlinear Programming (MINLP) Problems in Process Synthesis. Ind. Eng. Chem. Res. 1988, 27 (8), 1407-1421. https://doi.org/10.1021/ie00080a013
-    [2] Ovalle, D., Liñán, D. A., Lee, A., Gómez, J. M., Ricardez-Sandoval, L., Grossmann, I. E., & Neira, D. E. B. (2024). Logic-Based Discrete-Steepest Descent: A Solution Method for Process Synthesis Generalized Disjunctive Programs. arXiv preprint arXiv:2405.05358. https://doi.org/10.48550/arXiv.2405.05358
+References
+----------
+[1] Kocis, G. R.; Grossmann, I. E. Global Optimization of Nonconvex Mixed-Integer Nonlinear Programming (MINLP) Problems in Process Synthesis. Ind. Eng. Chem. Res. 1988, 27 (8), 1407-1421. https://doi.org/10.1021/ie00080a013
+[2] Ovalle, D., Liñán, D. A., Lee, A., Gómez, J. M., Ricardez-Sandoval, L., Grossmann, I. E., & Neira, D. E. B. (2024). Logic-Based Discrete-Steepest Descent: A Solution Method for Process Synthesis Generalized Disjunctive Programs. arXiv preprint arXiv:2405.05358. https://doi.org/10.48550/arXiv.2405.05358
 """
 import os
-
 import pyomo.environ as pyo
 from pyomo.core.base.misc import display
 from pyomo.core.plugins.transform.logical_to_linear import (

From 171fc1c2eddcc4210ffbcc67d7ef980954e07cf1 Mon Sep 17 00:00:00 2001
From: ZedongPeng <peng_zedong@126.com>
Date: Wed, 12 Jun 2024 16:24:32 +0800
Subject: [PATCH 7/9] black format

---
 gdplib/small_batch/gdp_small_batch.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/gdplib/small_batch/gdp_small_batch.py b/gdplib/small_batch/gdp_small_batch.py
index 1922aed..67f7467 100644
--- a/gdplib/small_batch/gdp_small_batch.py
+++ b/gdplib/small_batch/gdp_small_batch.py
@@ -9,6 +9,7 @@
 [1] Kocis, G. R.; Grossmann, I. E. Global Optimization of Nonconvex Mixed-Integer Nonlinear Programming (MINLP) Problems in Process Synthesis. Ind. Eng. Chem. Res. 1988, 27 (8), 1407-1421. https://doi.org/10.1021/ie00080a013
 [2] Ovalle, D., Liñán, D. A., Lee, A., Gómez, J. M., Ricardez-Sandoval, L., Grossmann, I. E., & Neira, D. E. B. (2024). Logic-Based Discrete-Steepest Descent: A Solution Method for Process Synthesis Generalized Disjunctive Programs. arXiv preprint arXiv:2405.05358. https://doi.org/10.48550/arXiv.2405.05358
 """
+
 import os
 import pyomo.environ as pyo
 from pyomo.core.base.misc import display

From 7a05fb033de74f0e962286231b44207960d0d369 Mon Sep 17 00:00:00 2001
From: Albert Lee <lee4382@purdue.edu>
Date: Fri, 16 Aug 2024 15:51:56 -0400
Subject: [PATCH 8/9] Added README.md for gdp_small_batch.py

---
 gdplib/small_batch/README.md | 15 +++++++++++++++
 1 file changed, 15 insertions(+)
 create mode 100644 gdplib/small_batch/README.md

diff --git a/gdplib/small_batch/README.md b/gdplib/small_batch/README.md
new file mode 100644
index 0000000..a583c06
--- /dev/null
+++ b/gdplib/small_batch/README.md
@@ -0,0 +1,15 @@
+## gdp_small_batch.py
+
+The gdp_small_batch.py module contains the GDP model for the small batch problem based on the Kocis and Grossmann (1988) paper.
+
+The problem is based on the Example 4 of the paper.
+
+The objective is to minimize the investment cost of the batch units.
+
+The solution is 167427.65711.
+
+### References
+
+[1] Kocis, G. R.; Grossmann, I. E. Global Optimization of Nonconvex Mixed-Integer Nonlinear Programming (MINLP) Problems in Process Synthesis. Ind. Eng. Chem. Res. 1988, 27 (8), 1407-1421. https://doi.org/10.1021/ie00080a013
+
+[2] Ovalle, D., Liñán, D. A., Lee, A., Gómez, J. M., Ricardez-Sandoval, L., Grossmann, I. E., & Neira, D. E. B. (2024). Logic-Based Discrete-Steepest Descent: A Solution Method for Process Synthesis Generalized Disjunctive Programs. arXiv preprint arXiv:2405.05358. https://doi.org/10.48550/arXiv.2405.05358

From 493dc9b7e8f08d2f4527b1e005df0b4cad6f24ea Mon Sep 17 00:00:00 2001
From: Albert Lee <lee4382@purdue.edu>
Date: Mon, 19 Aug 2024 12:18:21 -0400
Subject: [PATCH 9/9] Fixed the References.

---
 gdplib/small_batch/README.md          | 2 +-
 gdplib/small_batch/gdp_small_batch.py | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/gdplib/small_batch/README.md b/gdplib/small_batch/README.md
index a583c06..837d76e 100644
--- a/gdplib/small_batch/README.md
+++ b/gdplib/small_batch/README.md
@@ -12,4 +12,4 @@ The solution is 167427.65711.
 
 [1] Kocis, G. R.; Grossmann, I. E. Global Optimization of Nonconvex Mixed-Integer Nonlinear Programming (MINLP) Problems in Process Synthesis. Ind. Eng. Chem. Res. 1988, 27 (8), 1407-1421. https://doi.org/10.1021/ie00080a013
 
-[2] Ovalle, D., Liñán, D. A., Lee, A., Gómez, J. M., Ricardez-Sandoval, L., Grossmann, I. E., & Neira, D. E. B. (2024). Logic-Based Discrete-Steepest Descent: A Solution Method for Process Synthesis Generalized Disjunctive Programs. arXiv preprint arXiv:2405.05358. https://doi.org/10.48550/arXiv.2405.05358
+[2] Ovalle, D., Liñán, D. A., Lee, A., Gómez, J. M., Ricardez-Sandoval, L., Grossmann, I. E., & Bernal Neira, D. E. (2024). Logic-Based Discrete-Steepest Descent: A Solution Method for Process Synthesis Generalized Disjunctive Programs. arXiv preprint arXiv:2405.05358. https://doi.org/10.48550/arXiv.2405.05358
diff --git a/gdplib/small_batch/gdp_small_batch.py b/gdplib/small_batch/gdp_small_batch.py
index 67f7467..b7d2b48 100644
--- a/gdplib/small_batch/gdp_small_batch.py
+++ b/gdplib/small_batch/gdp_small_batch.py
@@ -7,7 +7,7 @@
 References
 ----------
 [1] Kocis, G. R.; Grossmann, I. E. Global Optimization of Nonconvex Mixed-Integer Nonlinear Programming (MINLP) Problems in Process Synthesis. Ind. Eng. Chem. Res. 1988, 27 (8), 1407-1421. https://doi.org/10.1021/ie00080a013
-[2] Ovalle, D., Liñán, D. A., Lee, A., Gómez, J. M., Ricardez-Sandoval, L., Grossmann, I. E., & Neira, D. E. B. (2024). Logic-Based Discrete-Steepest Descent: A Solution Method for Process Synthesis Generalized Disjunctive Programs. arXiv preprint arXiv:2405.05358. https://doi.org/10.48550/arXiv.2405.05358
+[2] Ovalle, D., Liñán, D. A., Lee, A., Gómez, J. M., Ricardez-Sandoval, L., Grossmann, I. E., & Bernal Neira, D. E. (2024). Logic-Based Discrete-Steepest Descent: A Solution Method for Process Synthesis Generalized Disjunctive Programs. arXiv preprint arXiv:2405.05358. https://doi.org/10.48550/arXiv.2405.05358
 """
 
 import os