From b308fb05dffd30a601dda9c958662a75110a9c69 Mon Sep 17 00:00:00 2001 From: Albert Lee Date: Fri, 17 May 2024 17:33:53 -0400 Subject: [PATCH 01/15] Added the gdp_reactor file on the gdp lib. --- gdplib/cstr/gdp_reactor.py | 741 +++++++++++++++++++++++++++++++++++++ 1 file changed, 741 insertions(+) create mode 100644 gdplib/cstr/gdp_reactor.py diff --git a/gdplib/cstr/gdp_reactor.py b/gdplib/cstr/gdp_reactor.py new file mode 100644 index 0000000..c34d35e --- /dev/null +++ b/gdplib/cstr/gdp_reactor.py @@ -0,0 +1,741 @@ +import os +import sys +import pyomo.environ as pyo +from pyomo.core.base.misc import display +from pyomo.gdp import Disjunct, Disjunction +from pyomo.opt.base.solvers import SolverFactory + + +def build_cstrs(NT: int = 5) -> pyo.ConcreteModel(): + """ + Function that builds CSTR superstructure model of size NT. + The CSTRs have a single 1st order reaction A -> B and minimizes (TODO Check) + total reactor volume. The optimal solution should yield NT reactors with a recycle before reactor NT. + Reference: Paper Linhan 1. TODO Correct reference + + Args: + NT: int. Positive Integer defining the maximum number of CSTRs + Returns: + m = Pyomo GDP model + """ + + # PYOMO MODEL + m = pyo.ConcreteModel(name='gdp_reactors') + + # SETS + m.I = pyo.Set(initialize=['A', 'B'], doc='Set of components') + m.N = pyo.RangeSet(1, NT, doc='Set of units in the superstructure') + + # PARAMETERS + m.k = pyo.Param(initialize=2, ) # Kinetic constant [L/(mol*s)] + m.order1 = pyo.Param(initialize=1) # Partial order of reacton 1 + m.order2 = pyo.Param(initialize=1) # Partial order of reaction 2 + m.QF0 = pyo.Param(initialize=1) # Inlet volumetric flow [L/s] + C0_Def = {'A': 0.99, 'B': 0.01} + + # Initial concentration of reagents [mol/L] + m.C0 = pyo.Param(m.I, initialize=C0_Def) + + # Inlet molar flow [mol/s] + + def F0_Def(m, i): + """_summary_ + + Parameters + ---------- + m : _type_ + _description_ + i : _type_ + _description_ + + Returns + ------- + _type_ + _description_ + """ + return m.C0[i]*m.QF0 + m.F0 = pyo.Param(m.I, initialize=F0_Def) + + # BOOLEAN VARIABLES + + # Unreacted feed in reactor n + m.YF = pyo.BooleanVar(m.N) + + # Existence of recycle flow in unit n + m.YR = pyo.BooleanVar(m.N) + + # Unit operation in n (True if unit n is a CSTR, False if unit n is a bypass) + m.YP = pyo.BooleanVar(m.N) + + # REAL VARIABLES + + # Network Variables + # Outlet flow rate of the superstructure unit [L/s] + m.Q = pyo.Var(m.N, initialize=0, within=pyo.NonNegativeReals, bounds=(0, 10)) + + # Outlet flow rate recycle activation of the superstructure unit [L/s] + m.QFR = pyo.Var(m.N, initialize=0, + within=pyo.NonNegativeReals, bounds=(0, 10)) + + # Molar flow [mol/s] + m.F = pyo.Var(m.I, m.N, initialize=0, + within=pyo.NonNegativeReals, bounds=(0, 10)) + + # Molar flow recycle activation [mol/s] + m.FR = pyo.Var(m.I, m.N, initialize=0, + within=pyo.NonNegativeReals, bounds=(0, 10)) + + # Reaction rate [mol/(L*s)] + m.rate = pyo.Var(m.I, m.N, initialize=0, within=pyo.Reals, bounds=(-10, 10)) + + # Reactor volume [L] + m.V = pyo.Var(m.N, initialize=0, within=pyo.NonNegativeReals, bounds=(0, 10)) + + # Volume activation [L] + m.c = pyo.Var(m.N, initialize=0, within=pyo.NonNegativeReals, bounds=(0, 10)) + + # Splitter Variables + # Recycle flow rate [L/s] + m.QR = pyo.Var(initialize=0, within=pyo.NonNegativeReals, bounds=(0, 10)) + + # Product flow rate [L/s] + m.QP = pyo.Var(initialize=0, within=pyo.NonNegativeReals, bounds=(0, 10)) + + # Recycle molar flow [mol/s] + m.R = pyo.Var(m.I, initialize=0, within=pyo.NonNegativeReals, bounds=(0, 10)) + + # Product molar flow [mol/s] + m.P = pyo.Var(m.I, initialize=0, within=pyo.NonNegativeReals, bounds=(0, 10)) + + # CONSTRAINTS + + # Unreacted Feed Balances + # Unreacted feed unit mole balance + + def unreact_mole_rule(m, i, n): + """_summary_ + + Parameters + ---------- + m : _type_ + _description_ + i : _type_ + _description_ + n : _type_ + _description_ + + Returns + ------- + _type_ + _description_ + """ + if n == NT: + return m.F0[i] + m.FR[i, n] - m.F[i, n] + m.rate[i, n]*m.V[n] == 0 + else: + return pyo.Constraint.Skip + + m.unreact_mole = pyo.Constraint(m.I, m.N, rule=unreact_mole_rule) + + # Unreacted feed unit continuity + + def unreact_cont_rule(m, n): + """_summary_ + + Parameters + ---------- + m : _type_ + _description_ + n : _type_ + _description_ + + Returns + ------- + _type_ + _description_ + """ + if n == NT: + return m.QF0 + m.QFR[n] - m.Q[n] == 0 + else: + return pyo.Constraint.Skip + + m.unreact_cont = pyo.Constraint(m.N, rule=unreact_cont_rule) + + # Reactor Balances + # Reactor mole balance + + def react_mole_rule(m, i, n): + """_summary_ + + Parameters + ---------- + m : _type_ + _description_ + i : _type_ + _description_ + n : _type_ + _description_ + + Returns + ------- + _type_ + _description_ + """ + if n != NT: + return m.F[i, n+1] + m.FR[i, n] - m.F[i, n] + m.rate[i, n]*m.V[n] == 0 + else: + return pyo.Constraint.Skip + + m.react_mole = pyo.Constraint(m.I, m.N, rule=react_mole_rule) + + # Reactor continuity + + def react_cont_rule(m, n): + """_summary_ + + Parameters + ---------- + m : _type_ + _description_ + n : _type_ + _description_ + + Returns + ------- + _type_ + _description_ + """ + if n != NT: + return m.Q[n+1] + m.QFR[n] - m.Q[n] == 0 + else: + return pyo.Constraint.Skip + + m.react_cont = pyo.Constraint(m.N, rule=react_cont_rule) + + # Splitting Point Balances + # Splitting point mole balance + + def split_mole_rule(m, i): + """_summary_ + + Parameters + ---------- + m : _type_ + _description_ + i : _type_ + _description_ + + Returns + ------- + _type_ + _description_ + """ + return m.F[i, 1] - m.P[i] - m.R[i] == 0 + + m.split_mole = pyo.Constraint(m.I, rule=split_mole_rule) + + # Splitting point continuity + + def split_cont_rule(m): + """_summary_ + + Parameters + ---------- + m : _type_ + _description_ + + Returns + ------- + _type_ + _description_ + """ + return m.Q[1] - m.QP - m.QR == 0 + + m.split_cont = pyo.Constraint(rule=split_cont_rule) + + # Splitting point additional constraints + + def split_add_rule(m, i): + """_summary_ + + Parameters + ---------- + m : _type_ + _description_ + i : _type_ + _description_ + + Returns + ------- + _type_ + _description_ + """ + return m.P[i]*m.Q[1] - m.F[i, 1]*m.QP == 0 + + m.split_add = pyo.Constraint(m.I, rule=split_add_rule) + + # Product Specification + + def prod_spec_rule(m): + """_summary_ + + Parameters + ---------- + m : _type_ + _description_ + + Returns + ------- + _type_ + _description_ + """ + return m.QP*0.95 - m.P['B'] == 0 + + m.prod_spec = pyo.Constraint(rule=prod_spec_rule) + + # Volume Constraint + + def vol_cons_rule(m, n): + """_summary_ + + Parameters + ---------- + m : _type_ + _description_ + n : _type_ + _description_ + + Returns + ------- + _type_ + _description_ + """ + if n != 1: + return m.V[n] - m.V[n-1] == 0 + else: + return pyo.Constraint.Skip + + m.vol_cons = pyo.Constraint(m.N, rule=vol_cons_rule) + + # YD Disjunction block equation definition + + def build_cstr_equations(disjunct, n): + """_summary_ + + Parameters + ---------- + disjunct : _type_ + _description_ + n : _type_ + _description_ + + Returns + ------- + _type_ + _description_ + """ + m = disjunct.model() + + # Reaction rates calculation + @disjunct.Constraint() + def YPD_rate_calc(disjunct): + """_summary_ + + Parameters + ---------- + disjunct : _type_ + _description_ + + Returns + ------- + _type_ + _description_ + """ + return m.rate['A', n]*((m.Q[n])**m.order1)*((m.Q[n])**m.order2)+m.k*((m.F['A', n])**m.order1)*((m.F['B', n])**m.order2) == 0 + + # Reaction rate relation + @disjunct.Constraint() + def YPD_rate_rel(disjunct): + """_summary_ + + Parameters + ---------- + disjunct : _type_ + _description_ + + Returns + ------- + _type_ + _description_ + """ + return m.rate['B', n] + m.rate['A', n] == 0 + + # Volume activation + @disjunct.Constraint() + def YPD_vol_act(disjunct): + """_summary_ + + Parameters + ---------- + disjunct : _type_ + _description_ + + Returns + ------- + _type_ + _description_ + """ + return m.c[n] - m.V[n] == 0 + + def build_bypass_equations(disjunct, n): + """_summary_ + + Parameters + ---------- + disjunct : _type_ + _description_ + n : _type_ + _description_ + + Returns + ------- + _type_ + _description_ + """ + m = disjunct.model() + + # FR desactivation + @disjunct.Constraint(m.I) + def neg_YPD_FR_desact(disjunct, i): + """_summary_ + + Parameters + ---------- + disjunct : _type_ + _description_ + i : _type_ + _description_ + + Returns + ------- + _type_ + _description_ + """ + return m.FR[i, n] == 0 + + # Rate desactivation + @disjunct.Constraint(m.I) + def neg_YPD_rate_desact(disjunct, i): + """_summary_ + + Parameters + ---------- + disjunct : _type_ + _description_ + i : _type_ + _description_ + + Returns + ------- + _type_ + _description_ + """ + return m.rate[i, n] == 0 + + # QFR desactivation + @disjunct.Constraint() + def neg_YPD_QFR_desact(disjunct): + """_summary_ + + Parameters + ---------- + disjunct : _type_ + _description_ + + Returns + ------- + _type_ + _description_ + """ + return m.QFR[n] == 0 + + @disjunct.Constraint() + def neg_YPD_vol_desact(disjunct): + ''' + Volume desactivation function for defining pyomo model + args: + disjunct: pyomo block with disjunct to include the constraint + n: pyomo set with reactor index + return: + return constraint + ''' + return m.c[n] == 0 + + # YR Disjuction block equation definition + + def build_recycle_equations(disjunct, n): + """_summary_ + + Parameters + ---------- + disjunct : _type_ + _description_ + n : _type_ + _description_ + + Returns + ------- + _type_ + _description_ + """ + m = disjunct.model() + + # FR activation + @disjunct.Constraint(m.I) + def YRD_FR_act(disjunct, i): + """_summary_ + + Parameters + ---------- + disjunct : _type_ + _description_ + i : _type_ + _description_ + + Returns + ------- + _type_ + _description_ + """ + return m.FR[i, n] - m.R[i] == 0 + + # QFR activation + @disjunct.Constraint() + def YRD_QFR_act(disjunct): + """_summary_ + + Parameters + ---------- + disjunct : _type_ + _description_ + + Returns + ------- + _type_ + _description_ + """ + return m.QFR[n] - m.QR == 0 + + def build_no_recycle_equations(disjunct, n): + """_summary_ + + Parameters + ---------- + disjunct : _type_ + _description_ + n : _type_ + _description_ + + Returns + ------- + _type_ + _description_ + """ + m = disjunct.model() + + # FR desactivation + @disjunct.Constraint(m.I) + def neg_YRD_FR_desact(disjunct, i): + """_summary_ + + Parameters + ---------- + disjunct : _type_ + _description_ + i : _type_ + _description_ + + Returns + ------- + _type_ + _description_ + """ + return m.FR[i, n] == 0 + + # QFR desactivation + @disjunct.Constraint() + def neg_YRD_QFR_desact(disjunct): + """_summary_ + + Parameters + ---------- + disjunct : _type_ + _description_ + + Returns + ------- + _type_ + _description_ + """ + return m.QFR[n] == 0 + + # Create disjunction blocks + m.YR_is_recycle = Disjunct(m.N, rule=build_recycle_equations) + m.YR_is_not_recycle = Disjunct(m.N, rule=build_no_recycle_equations) + + m.YP_is_cstr = Disjunct(m.N, rule=build_cstr_equations) + m.YP_is_bypass = Disjunct(m.N, rule=build_bypass_equations) + + # Create disjunctions + + @m.Disjunction(m.N) + def YP_is_cstr_or_bypass(m, n): + """_summary_ + + Parameters + ---------- + m : _type_ + _description_ + n : _type_ + _description_ + + Returns + ------- + _type_ + _description_ + """ + return [m.YP_is_cstr[n], m.YP_is_bypass[n]] + + @m.Disjunction(m.N) + def YR_is_recycle_or_not(m, n): + """_summary_ + + Parameters + ---------- + m : _type_ + _description_ + n : _type_ + _description_ + + Returns + ------- + _type_ + _description_ + """ + return [m.YR_is_recycle[n], m.YR_is_not_recycle[n]] + + # Associate Boolean variables with with disjunctions + for n in m.N: + m.YP[n].associate_binary_var(m.YP_is_cstr[n].indicator_var) + m.YR[n].associate_binary_var(m.YR_is_recycle[n].indicator_var) + + # Logic Constraints + # Unit must be a CSTR to include a recycle + + def cstr_if_recycle_rule(m, n): + """_summary_ + + Parameters + ---------- + m : _type_ + _description_ + n : _type_ + _description_ + + Returns + ------- + _type_ + _description_ + """ + return m.YR[n].implies(m.YP[n]) + + m.cstr_if_recycle = pyo.LogicalConstraint(m.N, rule=cstr_if_recycle_rule) + + # There is only one unreacted feed + + def one_unreacted_feed_rule(m): + """_summary_ + + Parameters + ---------- + m : _type_ + _description_ + + Returns + ------- + _type_ + _description_ + """ + return pyo.exactly(1, m.YF) + + m.one_unreacted_feed = pyo.LogicalConstraint(rule=one_unreacted_feed_rule) + + # There is only one recycle stream + + def one_recycle_rule(m): + """_summary_ + + Parameters + ---------- + m : _type_ + _description_ + + Returns + ------- + _type_ + _description_ + """ + return pyo.exactly(1, m.YR) + + m.one_recycle = pyo.LogicalConstraint(rule=one_recycle_rule) + + # Unit operation in n constraint + + def unit_in_n_rule(m, n): + """_summary_ + + Parameters + ---------- + m : _type_ + _description_ + n : _type_ + _description_ + + Returns + ------- + _type_ + _description_ + """ + if n == 1: + return m.YP[n].equivalent_to(True) + else: + return m.YP[n].equivalent_to(pyo.lor(pyo.land(~m.YF[n2] for n2 in range(1, n)), m.YF[n])) + + m.unit_in_n = pyo.LogicalConstraint(m.N, rule=unit_in_n_rule) + + # OBJECTIVE + + def obj_rule(m): + """_summary_ + + Parameters + ---------- + m : _type_ + _description_ + + Returns + ------- + _type_ + _description_ + """ + return sum(m.c[n] for n in m.N) + + m.obj = pyo.Objective(rule=obj_rule, sense=pyo.minimize) + + return m + +if __name__ == "__main__": + m = build_cstrs() + 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) \ No newline at end of file From 0aa9e0c1c23db72b773831f8f8daab3fb45c38fe Mon Sep 17 00:00:00 2001 From: Albert Lee Date: Fri, 17 May 2024 17:34:13 -0400 Subject: [PATCH 02/15] Added the simple documentation of the cstr using readme. --- gdplib/cstr/README.md | 30 ++++++++++++++++++++++++++++++ 1 file changed, 30 insertions(+) create mode 100644 gdplib/cstr/README.md diff --git a/gdplib/cstr/README.md b/gdplib/cstr/README.md new file mode 100644 index 0000000..18882a4 --- /dev/null +++ b/gdplib/cstr/README.md @@ -0,0 +1,30 @@ +# GDP Reactor Series Design +Function that builds CSTR superstructure model of size NT (default = 5). +NT is the number of reactors in series. +The CSTRs have a single 1st order auto catalytic reaction A -> B and minimizes total reactors series volume. +The optimal solution should yield NT reactors with a recycle before reactor NT. + +Reference: +> Linan, D. A., Bernal, D. E., Gomez, J. M., & Ricardez-Sandoval, L. A. (2021). Optimal synthesis and design of catalytic distillation columns: A rate-based modeling approach. Chemical Engineering Science, 231, 116294. https://doi.org/10.1016/j.ces.2020.116294 + +### Solution + +Best known objective value: 3.06181298849707 + +### Size + +Number of reactors in series is 5. + +| Problem | vars | Bool | bin | int | cont | cons | nl | disj | disjtn | +|-----------|------|------|-----|-----|------|------|----|------|--------| +| gdp_reactors | 71 | 15 | 0 | 0 | 56 | 25 | 2 | 20 | 10 | + +- ``vars``: variables +- ``Bool``: Boolean variables +- ``bin``: binary variables +- ``int``: integer variables +- ``cont``: continuous variables +- ``cons``: constraints +- ``nl``: nonlinear constraints +- ``disj``: disjuncts +- ``disjtn``: disjunctions \ No newline at end of file From d91ca9d994823c2f5008a555206758d6eb32ac99 Mon Sep 17 00:00:00 2001 From: Albert Lee Date: Fri, 17 May 2024 18:01:34 -0400 Subject: [PATCH 03/15] gdp reactor series add docs. --- gdplib/cstr/gdp_reactor.py | 69 ++++++++++++++++++++++---------------- 1 file changed, 41 insertions(+), 28 deletions(-) diff --git a/gdplib/cstr/gdp_reactor.py b/gdplib/cstr/gdp_reactor.py index c34d35e..f55cd7d 100644 --- a/gdplib/cstr/gdp_reactor.py +++ b/gdplib/cstr/gdp_reactor.py @@ -7,16 +7,29 @@ def build_cstrs(NT: int = 5) -> pyo.ConcreteModel(): + # """ + # Function that builds CSTR superstructure model of size NT. + # The CSTRs have a single 1st order reaction A -> B and minimizes (TODO Check) + # total reactor volume. The optimal solution should yield NT reactors with a recycle before reactor NT. + # Reference: Paper Linhan 1. TODO Correct reference + + # Args: + # NT: int. Positive Integer defining the maximum number of CSTRs + # Returns: + # m = Pyomo GDP model + # """ """ - Function that builds CSTR superstructure model of size NT. - The CSTRs have a single 1st order reaction A -> B and minimizes (TODO Check) - total reactor volume. The optimal solution should yield NT reactors with a recycle before reactor NT. - Reference: Paper Linhan 1. TODO Correct reference - - Args: - NT: int. Positive Integer defining the maximum number of CSTRs - Returns: - m = Pyomo GDP model + Build the CSTR superstructure model of size NT. + + + Parameters + ---------- + NT : int + Number of possible reactors in the reactor series superstructure + Returns + ------- + m : Pyomo.ConcreteModel + _description_ """ # PYOMO MODEL @@ -27,8 +40,8 @@ def build_cstrs(NT: int = 5) -> pyo.ConcreteModel(): m.N = pyo.RangeSet(1, NT, doc='Set of units in the superstructure') # PARAMETERS - m.k = pyo.Param(initialize=2, ) # Kinetic constant [L/(mol*s)] - m.order1 = pyo.Param(initialize=1) # Partial order of reacton 1 + m.k = pyo.Param(initialize=2, doc="Kinetic constant [L/(mol*s)]") # Kinetic constant [L/(mol*s)] + m.order1 = pyo.Param(initialize=1, doc="Partial order of") # Partial order of reacton 1 m.order2 = pyo.Param(initialize=1) # Partial order of reaction 2 m.QF0 = pyo.Param(initialize=1) # Inlet volumetric flow [L/s] C0_Def = {'A': 0.99, 'B': 0.01} @@ -43,10 +56,10 @@ def F0_Def(m, i): Parameters ---------- - m : _type_ - _description_ - i : _type_ + m : Pyomo.ConcreteModel _description_ + i : float + Index of the component in the reactor series. Returns ------- @@ -59,53 +72,53 @@ def F0_Def(m, i): # BOOLEAN VARIABLES # Unreacted feed in reactor n - m.YF = pyo.BooleanVar(m.N) + m.YF = pyo.BooleanVar(m.N, doc="Unreacted feed in reactor n") # Existence of recycle flow in unit n - m.YR = pyo.BooleanVar(m.N) + m.YR = pyo.BooleanVar(m.N, doc="Existence of recycle flow in unit n") # Unit operation in n (True if unit n is a CSTR, False if unit n is a bypass) - m.YP = pyo.BooleanVar(m.N) + m.YP = pyo.BooleanVar(m.N, doc="Unit operation in n") # REAL VARIABLES # Network Variables # Outlet flow rate of the superstructure unit [L/s] - m.Q = pyo.Var(m.N, initialize=0, within=pyo.NonNegativeReals, bounds=(0, 10)) + m.Q = pyo.Var(m.N, initialize=0, within=pyo.NonNegativeReals, bounds=(0, 10), doc="Outlet flow rate of the superstructure unit [L/s]") # Outlet flow rate recycle activation of the superstructure unit [L/s] m.QFR = pyo.Var(m.N, initialize=0, - within=pyo.NonNegativeReals, bounds=(0, 10)) + within=pyo.NonNegativeReals, bounds=(0, 10), doc="Outlet flow rate recycle activation of the superstructure unit [L/s]") # Molar flow [mol/s] m.F = pyo.Var(m.I, m.N, initialize=0, - within=pyo.NonNegativeReals, bounds=(0, 10)) + within=pyo.NonNegativeReals, bounds=(0, 10), doc="Molar flow [mol/s]") # Molar flow recycle activation [mol/s] m.FR = pyo.Var(m.I, m.N, initialize=0, - within=pyo.NonNegativeReals, bounds=(0, 10)) + within=pyo.NonNegativeReals, bounds=(0, 10), doc="Molar flow recycle activation [mol/s]") # Reaction rate [mol/(L*s)] - m.rate = pyo.Var(m.I, m.N, initialize=0, within=pyo.Reals, bounds=(-10, 10)) + m.rate = pyo.Var(m.I, m.N, initialize=0, within=pyo.Reals, bounds=(-10, 10), doc="Reaction rate [mol/(L*s)]") # Reactor volume [L] - m.V = pyo.Var(m.N, initialize=0, within=pyo.NonNegativeReals, bounds=(0, 10)) + m.V = pyo.Var(m.N, initialize=0, within=pyo.NonNegativeReals, bounds=(0, 10), doc="Reactor volume [L]") # Volume activation [L] - m.c = pyo.Var(m.N, initialize=0, within=pyo.NonNegativeReals, bounds=(0, 10)) + m.c = pyo.Var(m.N, initialize=0, within=pyo.NonNegativeReals, bounds=(0, 10), doc="Volume activation [L]") # Splitter Variables # Recycle flow rate [L/s] - m.QR = pyo.Var(initialize=0, within=pyo.NonNegativeReals, bounds=(0, 10)) + m.QR = pyo.Var(initialize=0, within=pyo.NonNegativeReals, bounds=(0, 10), doc="Recycle flow rate [L/s]") # Product flow rate [L/s] - m.QP = pyo.Var(initialize=0, within=pyo.NonNegativeReals, bounds=(0, 10)) + m.QP = pyo.Var(initialize=0, within=pyo.NonNegativeReals, bounds=(0, 10), doc="Product flow rate [L/s]") # Recycle molar flow [mol/s] - m.R = pyo.Var(m.I, initialize=0, within=pyo.NonNegativeReals, bounds=(0, 10)) + m.R = pyo.Var(m.I, initialize=0, within=pyo.NonNegativeReals, bounds=(0, 10), doc="Recycle molar flow [mol/s]") # Product molar flow [mol/s] - m.P = pyo.Var(m.I, initialize=0, within=pyo.NonNegativeReals, bounds=(0, 10)) + m.P = pyo.Var(m.I, initialize=0, within=pyo.NonNegativeReals, bounds=(0, 10), doc="Product molar flow [mol/s]") # CONSTRAINTS From bdef4682c7e4198232ad3026032eb883c6c70c20 Mon Sep 17 00:00:00 2001 From: Albert Lee Date: Sat, 18 May 2024 17:55:47 -0400 Subject: [PATCH 04/15] add documentation on the inlet flow. --- gdplib/cstr/gdp_reactor.py | 210 ++++++++++++++++++------------------- 1 file changed, 103 insertions(+), 107 deletions(-) diff --git a/gdplib/cstr/gdp_reactor.py b/gdplib/cstr/gdp_reactor.py index f55cd7d..305424e 100644 --- a/gdplib/cstr/gdp_reactor.py +++ b/gdplib/cstr/gdp_reactor.py @@ -7,29 +7,21 @@ def build_cstrs(NT: int = 5) -> pyo.ConcreteModel(): - # """ - # Function that builds CSTR superstructure model of size NT. - # The CSTRs have a single 1st order reaction A -> B and minimizes (TODO Check) - # total reactor volume. The optimal solution should yield NT reactors with a recycle before reactor NT. - # Reference: Paper Linhan 1. TODO Correct reference - - # Args: - # NT: int. Positive Integer defining the maximum number of CSTRs - # Returns: - # m = Pyomo GDP model - # """ """ Build the CSTR superstructure model of size NT. - + NT is the number of reactors in series. + The CSTRs have a single 1st order auto catalytic reaction A -> B and minimizes total reactors series volume. + The optimal solution should yield NT reactors with a recycle before reactor NT. Parameters ---------- NT : int Number of possible reactors in the reactor series superstructure + Returns ------- m : Pyomo.ConcreteModel - _description_ + Pyomo GDP model which represents the superstructure model of size of NT reactors. """ # PYOMO MODEL @@ -41,30 +33,32 @@ def build_cstrs(NT: int = 5) -> pyo.ConcreteModel(): # PARAMETERS m.k = pyo.Param(initialize=2, doc="Kinetic constant [L/(mol*s)]") # Kinetic constant [L/(mol*s)] - m.order1 = pyo.Param(initialize=1, doc="Partial order of") # Partial order of reacton 1 - m.order2 = pyo.Param(initialize=1) # Partial order of reaction 2 - m.QF0 = pyo.Param(initialize=1) # Inlet volumetric flow [L/s] + m.order1 = pyo.Param(initialize=1, doc="Partial order of reaction 1") # Partial order of reacton 1 + m.order2 = pyo.Param(initialize=1, doc="Partial order of reaction 2") # Partial order of reaction 2 + m.QF0 = pyo.Param(initialize=1, doc="Inlet volumetric flow [L/s]") # Inlet volumetric flow [L/s] C0_Def = {'A': 0.99, 'B': 0.01} # Initial concentration of reagents [mol/L] - m.C0 = pyo.Param(m.I, initialize=C0_Def) + m.C0 = pyo.Param(m.I, initialize=C0_Def, doc="Initial concentration of reagents [mol/L]") # Inlet molar flow [mol/s] def F0_Def(m, i): - """_summary_ + """ + Inlet molar flow [mol/s] for component i. + The function multiplies the initial concentration of component i by the inlet volumetric flow. Parameters ---------- m : Pyomo.ConcreteModel - _description_ + Pyomo GDP model of the CSTR superstructure. i : float Index of the component in the reactor series. Returns ------- - _type_ - _description_ + Pyomo.Param + Inlet molar flow [mol/s] for component i """ return m.C0[i]*m.QF0 m.F0 = pyo.Param(m.I, initialize=F0_Def) @@ -126,20 +120,21 @@ def F0_Def(m, i): # Unreacted feed unit mole balance def unreact_mole_rule(m, i, n): - """_summary_ + """ + Unreacted feed unit mole balance. Parameters ---------- - m : _type_ - _description_ - i : _type_ - _description_ - n : _type_ - _description_ + m : Pyomo.ConcreteModel + Pyomo GDP model of the CSTR superstructure. + i : float + Index of the component in the reactor series. + n : int + Index of the reactor in the reactor series. The reactor index starts at 1. The number increases on the left direction Returns ------- - _type_ + Pyomo.Constraint or Pyomo.Constraint.Skip _description_ """ if n == NT: @@ -156,10 +151,10 @@ def unreact_cont_rule(m, n): Parameters ---------- - m : _type_ - _description_ - n : _type_ - _description_ + m : Pyomo.ConcreteModel + Pyomo GDP model of the CSTR superstructure. + n : int + Index of the reactor in the reactor series. Returns ------- @@ -181,12 +176,12 @@ def react_mole_rule(m, i, n): Parameters ---------- - m : _type_ - _description_ - i : _type_ - _description_ - n : _type_ - _description_ + m : Pyomo.ConcreteModel + Pyomo GDP model of the CSTR superstructure. + i : float + Index of the component in the reactor series. + n : int + Index of the reactor in the reactor series. Returns ------- @@ -207,10 +202,10 @@ def react_cont_rule(m, n): Parameters ---------- - m : _type_ - _description_ - n : _type_ - _description_ + m : Pyomo.ConcreteModel + Pyomo GDP model of the CSTR superstructure. + n : int + Index of the reactor in the reactor series. Returns ------- @@ -232,10 +227,10 @@ def split_mole_rule(m, i): Parameters ---------- - m : _type_ - _description_ - i : _type_ - _description_ + m : Pyomo.ConcreteModel + Pyomo GDP model of the CSTR superstructure. + i : float + Index of the component in the reactor series. Returns ------- @@ -253,8 +248,8 @@ def split_cont_rule(m): Parameters ---------- - m : _type_ - _description_ + m : Pyomo.ConcreteModel + Pyomo GDP model of the CSTR superstructure. Returns ------- @@ -272,10 +267,10 @@ def split_add_rule(m, i): Parameters ---------- - m : _type_ - _description_ - i : _type_ - _description_ + m : Pyomo.ConcreteModel + Pyomo GDP model of the CSTR superstructure. + i : float + Index of the component in the reactor series. Returns ------- @@ -293,8 +288,8 @@ def prod_spec_rule(m): Parameters ---------- - m : _type_ - _description_ + m : Pyomo.ConcreteModel + Pyomo GDP model of the CSTR superstructure. Returns ------- @@ -312,10 +307,10 @@ def vol_cons_rule(m, n): Parameters ---------- - m : _type_ - _description_ - n : _type_ - _description_ + m : Pyomo.ConcreteModel + Pyomo GDP model of the CSTR superstructure. + n : int + Index of the reactor in the reactor series. Returns ------- @@ -338,8 +333,8 @@ def build_cstr_equations(disjunct, n): ---------- disjunct : _type_ _description_ - n : _type_ - _description_ + n : int + Index of the reactor in the reactor series. Returns ------- @@ -406,8 +401,8 @@ def build_bypass_equations(disjunct, n): ---------- disjunct : _type_ _description_ - n : _type_ - _description_ + n : int + Index of the reactor in the reactor series. Returns ------- @@ -425,8 +420,8 @@ def neg_YPD_FR_desact(disjunct, i): ---------- disjunct : _type_ _description_ - i : _type_ - _description_ + i : float + Index of the component in the reactor series. Returns ------- @@ -444,8 +439,8 @@ def neg_YPD_rate_desact(disjunct, i): ---------- disjunct : _type_ _description_ - i : _type_ - _description_ + i : float + Index of the component in the reactor series. Returns ------- @@ -492,8 +487,8 @@ def build_recycle_equations(disjunct, n): ---------- disjunct : _type_ _description_ - n : _type_ - _description_ + n : int + Index of the reactor in the reactor series. Returns ------- @@ -511,8 +506,8 @@ def YRD_FR_act(disjunct, i): ---------- disjunct : _type_ _description_ - i : _type_ - _description_ + i : float + Index of the component in the reactor series. Returns ------- @@ -545,8 +540,8 @@ def build_no_recycle_equations(disjunct, n): ---------- disjunct : _type_ _description_ - n : _type_ - _description_ + n : int + Index of the reactor in the reactor series. Returns ------- @@ -564,9 +559,9 @@ def neg_YRD_FR_desact(disjunct, i): ---------- disjunct : _type_ _description_ - i : _type_ - _description_ - + i : float + Index of the component in the reactor series. + Returns ------- _type_ @@ -606,14 +601,14 @@ def YP_is_cstr_or_bypass(m, n): Parameters ---------- - m : _type_ - _description_ - n : _type_ - _description_ + m : Pyomo.ConcreteModel + Pyomo GDP model of the CSTR superstructure. + n : int + Index of the reactor in the reactor series. Returns ------- - _type_ + list _description_ """ return [m.YP_is_cstr[n], m.YP_is_bypass[n]] @@ -624,14 +619,14 @@ def YR_is_recycle_or_not(m, n): Parameters ---------- - m : _type_ - _description_ - n : _type_ - _description_ + m : Pyomo.ConcreteModel + Pyomo GDP model of the CSTR superstructure. + n : int + Index of the reactor in the reactor series. Returns ------- - _type_ + list _description_ """ return [m.YR_is_recycle[n], m.YR_is_not_recycle[n]] @@ -649,10 +644,10 @@ def cstr_if_recycle_rule(m, n): Parameters ---------- - m : _type_ - _description_ - n : _type_ - _description_ + m : Pyomo.ConcreteModel + Pyomo GDP model of the CSTR superstructure. + n : int + Index of the reactor in the reactor series. Returns ------- @@ -670,12 +665,12 @@ def one_unreacted_feed_rule(m): Parameters ---------- - m : _type_ - _description_ + m : Pyomo.ConcreteModel + Pyomo GDP model of the CSTR superstructure. Returns ------- - _type_ + Pyomo.LogicalConstraint _description_ """ return pyo.exactly(1, m.YF) @@ -689,12 +684,12 @@ def one_recycle_rule(m): Parameters ---------- - m : _type_ - _description_ + m : Pyomo.ConcreteModel + Pyomo GDP model of the CSTR superstructure. Returns ------- - _type_ + Pyomo.LogicalConstraint _description_ """ return pyo.exactly(1, m.YR) @@ -708,14 +703,14 @@ def unit_in_n_rule(m, n): Parameters ---------- - m : _type_ - _description_ - n : _type_ - _description_ + m : Pyomo.ConcreteModel + Pyomo GDP model of the CSTR superstructure. + n : int + Index of the reactor in the reactor series. Returns ------- - _type_ + Pyomo.LogicalConstraint _description_ """ if n == 1: @@ -728,21 +723,22 @@ def unit_in_n_rule(m, n): # OBJECTIVE def obj_rule(m): - """_summary_ + """ + Objective function to minimize the total reactor volume. Parameters ---------- - m : _type_ - _description_ + m : Pyomo.ConcreteModel + Pyomo GDP model of the CSTR superstructure. Returns ------- - _type_ - _description_ + Pyomo Objective + Objective function to minimize the total reactor volume. """ return sum(m.c[n] for n in m.N) - m.obj = pyo.Objective(rule=obj_rule, sense=pyo.minimize) + m.obj = pyo.Objective(rule=obj_rule, sense=pyo.minimize, doc="minimum total reactor volume") return m From faddeaf8087e008a1d72a79269d14dfa8b904684 Mon Sep 17 00:00:00 2001 From: Albert Lee Date: Sat, 18 May 2024 18:11:08 -0400 Subject: [PATCH 05/15] leftest reactor molar and volume balance. --- gdplib/cstr/gdp_reactor.py | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/gdplib/cstr/gdp_reactor.py b/gdplib/cstr/gdp_reactor.py index 305424e..e038850 100644 --- a/gdplib/cstr/gdp_reactor.py +++ b/gdplib/cstr/gdp_reactor.py @@ -122,6 +122,7 @@ def F0_Def(m, i): def unreact_mole_rule(m, i, n): """ Unreacted feed unit mole balance. + The mole balance is calculated using the inlet molar flow, the recycle molar flow, the outlet molar flow, and the reaction rate for each component. Parameters ---------- @@ -130,24 +131,26 @@ def unreact_mole_rule(m, i, n): i : float Index of the component in the reactor series. n : int - Index of the reactor in the reactor series. The reactor index starts at 1. The number increases on the left direction + Index of the reactor in the reactor series. The reactor index starts at 1. The number increases on the left direction of the reactor series. The last reactor is indexed as NT which is the feed for reagent. Returns ------- Pyomo.Constraint or Pyomo.Constraint.Skip - _description_ + The constraint for the unreacted feed unit mole balance if n is equal to NT. Otherwise, it returns Pyomo.Constraint.Skip. """ if n == NT: return m.F0[i] + m.FR[i, n] - m.F[i, n] + m.rate[i, n]*m.V[n] == 0 else: return pyo.Constraint.Skip - m.unreact_mole = pyo.Constraint(m.I, m.N, rule=unreact_mole_rule) + m.unreact_mole = pyo.Constraint(m.I, m.N, rule=unreact_mole_rule, doc="Unreacted feed unit mole balance") # Unreacted feed unit continuity def unreact_cont_rule(m, n): - """_summary_ + """ + Unreacted feed unit continuity. + The continuity is calculated using the inlet volumetric flow, the recycle flow rate, and the outlet flow rate for the reactor NT. Parameters ---------- @@ -158,15 +161,15 @@ def unreact_cont_rule(m, n): Returns ------- - _type_ - _description_ + Pyomo.Constraint or Pyomo.Constraint.Skip + The constraint for the unreacted feed unit continuity if n is equal to NT. Otherwise, it returns Pyomo.Constraint.Skip. """ if n == NT: return m.QF0 + m.QFR[n] - m.Q[n] == 0 else: return pyo.Constraint.Skip - m.unreact_cont = pyo.Constraint(m.N, rule=unreact_cont_rule) + m.unreact_cont = pyo.Constraint(m.N, rule=unreact_cont_rule, doc="Unreacted feed unit continuity") # Reactor Balances # Reactor mole balance From dc94ae153fad6dcfa4a5ba485a942142ebc2a355 Mon Sep 17 00:00:00 2001 From: Albert Lee Date: Sat, 18 May 2024 18:28:22 -0400 Subject: [PATCH 06/15] general constraint documentation --- gdplib/cstr/gdp_reactor.py | 70 +++++++++++++++++++++++--------------- 1 file changed, 42 insertions(+), 28 deletions(-) diff --git a/gdplib/cstr/gdp_reactor.py b/gdplib/cstr/gdp_reactor.py index e038850..0642019 100644 --- a/gdplib/cstr/gdp_reactor.py +++ b/gdplib/cstr/gdp_reactor.py @@ -175,7 +175,9 @@ def unreact_cont_rule(m, n): # Reactor mole balance def react_mole_rule(m, i, n): - """_summary_ + """ + Reactor mole balance. + The mole balance is calculated using the inlet molar flow, the recycle molar flow, the outlet molar flow, and the reaction rate for each component. Parameters ---------- @@ -188,20 +190,22 @@ def react_mole_rule(m, i, n): Returns ------- - _type_ - _description_ + Pyomo.Constraint or Pyomo.Constraint.Skip + The constraint for the reactor mole balance if n is different from NT. Otherwise, it returns Pyomo.Constraint.Skip. """ if n != NT: return m.F[i, n+1] + m.FR[i, n] - m.F[i, n] + m.rate[i, n]*m.V[n] == 0 else: return pyo.Constraint.Skip - m.react_mole = pyo.Constraint(m.I, m.N, rule=react_mole_rule) + m.react_mole = pyo.Constraint(m.I, m.N, rule=react_mole_rule, doc="Reactor mole balance") # Reactor continuity def react_cont_rule(m, n): - """_summary_ + """ + Reactor continuity. + The continuity is calculated using the inlet volumetric flow, the recycle flow rate, and the outlet flow rate for the reactor n. Parameters ---------- @@ -212,21 +216,23 @@ def react_cont_rule(m, n): Returns ------- - _type_ - _description_ + Pyomo.Constraint or Pyomo.Constraint.Skip + The constraint for the reactor continuity if n is different from NT. Otherwise, it returns Pyomo.Constraint.Skip. """ if n != NT: return m.Q[n+1] + m.QFR[n] - m.Q[n] == 0 else: return pyo.Constraint.Skip - m.react_cont = pyo.Constraint(m.N, rule=react_cont_rule) + m.react_cont = pyo.Constraint(m.N, rule=react_cont_rule, doc="Reactor continuity") # Splitting Point Balances # Splitting point mole balance def split_mole_rule(m, i): - """_summary_ + """ + Splitting point mole balance. + The mole balance is calculated using the product molar flow, the recycle molar flow, and the outlet molar flow for each component. Parameters ---------- @@ -237,17 +243,19 @@ def split_mole_rule(m, i): Returns ------- - _type_ - _description_ + Pyomo.Constraint + The constraint for the splitting point mole balance. """ return m.F[i, 1] - m.P[i] - m.R[i] == 0 - m.split_mole = pyo.Constraint(m.I, rule=split_mole_rule) + m.split_mole = pyo.Constraint(m.I, rule=split_mole_rule, doc="Splitting point mole balance") # Splitting point continuity def split_cont_rule(m): - """_summary_ + """ + Splitting point continuity. + The continuity is calculated using the product flow rate, the recycle flow rate, and the outlet flow rate for the reactor 1. Parameters ---------- @@ -256,17 +264,20 @@ def split_cont_rule(m): Returns ------- - _type_ - _description_ + Pyomo.Constraint + The constraint for the splitting point continuity. """ return m.Q[1] - m.QP - m.QR == 0 - m.split_cont = pyo.Constraint(rule=split_cont_rule) + m.split_cont = pyo.Constraint(rule=split_cont_rule, doc="Splitting point continuity") # Splitting point additional constraints def split_add_rule(m, i): - """_summary_ + """ + Splitting point additional constraints. + Molarity constraints over initial and final flows, read as an multiplication avoid the numerical complication. + m.P[i]/m.QP = m.F[i,1]/m.Q[1] (molarity balance) Parameters ---------- @@ -277,17 +288,19 @@ def split_add_rule(m, i): Returns ------- - _type_ - _description_ + Pyomo.Constraint + The constraint for the splitting point additional constraints. """ return m.P[i]*m.Q[1] - m.F[i, 1]*m.QP == 0 - m.split_add = pyo.Constraint(m.I, rule=split_add_rule) + m.split_add = pyo.Constraint(m.I, rule=split_add_rule, doc="Splitting point additional constraints") # Product Specification def prod_spec_rule(m): - """_summary_ + """ + Product B Specification. + The product B specification is calculated using the product flow rate and the product molar flow. Parameters ---------- @@ -296,17 +309,18 @@ def prod_spec_rule(m): Returns ------- - _type_ - _description_ + Pyomo.Constraint + The constraint for the product B specification. """ return m.QP*0.95 - m.P['B'] == 0 - m.prod_spec = pyo.Constraint(rule=prod_spec_rule) + m.prod_spec = pyo.Constraint(rule=prod_spec_rule, doc="Product B Specification") # Volume Constraint def vol_cons_rule(m, n): - """_summary_ + """ + Volume Constraint. Parameters ---------- @@ -317,15 +331,15 @@ def vol_cons_rule(m, n): Returns ------- - _type_ - _description_ + Pyomo.Constraint or Pyomo.Constraint.Skip + The constraint for the volume constraint if n is different from 1. Otherwise, it returns Pyomo.Constraint.Skip. """ if n != 1: return m.V[n] - m.V[n-1] == 0 else: return pyo.Constraint.Skip - m.vol_cons = pyo.Constraint(m.N, rule=vol_cons_rule) + m.vol_cons = pyo.Constraint(m.N, rule=vol_cons_rule, doc"Volume Constraint") # YD Disjunction block equation definition From f0c95ac8dcf5bc243d1ea35eb6afbbb6ee34c1e4 Mon Sep 17 00:00:00 2001 From: Albert Lee Date: Sun, 19 May 2024 13:45:35 -0400 Subject: [PATCH 07/15] Complete the documentation of the GDP reactor. --- gdplib/cstr/gdp_reactor.py | 241 +++++++++++++++++++++---------------- 1 file changed, 139 insertions(+), 102 deletions(-) diff --git a/gdplib/cstr/gdp_reactor.py b/gdplib/cstr/gdp_reactor.py index 0642019..b88475d 100644 --- a/gdplib/cstr/gdp_reactor.py +++ b/gdplib/cstr/gdp_reactor.py @@ -339,7 +339,7 @@ def vol_cons_rule(m, n): else: return pyo.Constraint.Skip - m.vol_cons = pyo.Constraint(m.N, rule=vol_cons_rule, doc"Volume Constraint") + m.vol_cons = pyo.Constraint(m.N, rule=vol_cons_rule, doc="Volume Constraint") # YD Disjunction block equation definition @@ -348,48 +348,55 @@ def build_cstr_equations(disjunct, n): Parameters ---------- - disjunct : _type_ - _description_ + disjunct : Pyomo.Disjunct + Pyomo Disjunct block to include the constraints for the activation of the CSTR reactor. n : int Index of the reactor in the reactor series. Returns ------- - _type_ - _description_ + None + None, the function builds the constraints for the activation of the CSTR reactor. """ m = disjunct.model() # Reaction rates calculation @disjunct.Constraint() def YPD_rate_calc(disjunct): - """_summary_ + """ + Calculate the reaction rate of A for each reactor. + The reaction rate is calculated using the kinetic constant, the outlet flow rate, the molar flow of A, and the molar flow of B. + The outlet flow is multiplied on the reaction rate to avoid numerical complications. Parameters ---------- - disjunct : _type_ - _description_ + disjunct : Pyomo.Disjunct + Pyomo Disjunct block to include the constraints for the activation of the CSTR reactor. Returns ------- - _type_ - _description_ + Pyomo.Constraint + The constraint for the calculation of the reaction rate of A for each reactor. """ return m.rate['A', n]*((m.Q[n])**m.order1)*((m.Q[n])**m.order2)+m.k*((m.F['A', n])**m.order1)*((m.F['B', n])**m.order2) == 0 # Reaction rate relation @disjunct.Constraint() def YPD_rate_rel(disjunct): - """_summary_ + """ + Reaction rate relation for defining pyomo model. + Since the chemical reaction goes from A to B, the rate of A is equal to the negative rate of B. + + A -> B Parameters ---------- - disjunct : _type_ - _description_ + disjunct : Pyomo.Disjunct + Pyomo Disjunct block to include the constraints for the activation of the CSTR reactor. Returns ------- - _type_ + Pyomo.Constraint _description_ """ return m.rate['B', n] + m.rate['A', n] == 0 @@ -397,224 +404,245 @@ def YPD_rate_rel(disjunct): # Volume activation @disjunct.Constraint() def YPD_vol_act(disjunct): - """_summary_ + """ + Volume Activation function for defining pyomo model. Parameters ---------- - disjunct : _type_ - _description_ + disjunct : Pyomo.Disjunct + Pyomo Disjunct block to include the constraints for the activation of the CSTR reactor. Returns ------- - _type_ - _description_ + Pyomo.Constraint + Activation function for the volume of the reactor. """ return m.c[n] - m.V[n] == 0 def build_bypass_equations(disjunct, n): - """_summary_ + """ + Build the constraints for the deactivation of the reactor (bypass the reactor). Parameters ---------- - disjunct : _type_ - _description_ + disjunct : Pyomo.Disjunct + Pyomo Disjunct block to include the constraints for the deactivation of the reactor (bypass the reactor). n : int Index of the reactor in the reactor series. Returns ------- - _type_ - _description_ + None + None, the function builds the constraints for the deactivation of the reactor (bypass the reactor). """ m = disjunct.model() - # FR desactivation + # FR deactivation @disjunct.Constraint(m.I) def neg_YPD_FR_desact(disjunct, i): - """_summary_ + """ + Deactivation of the recycle flow for each component in the reactor series. + There are no recycle flows when the reactor is deactivated (bypassed). Parameters ---------- - disjunct : _type_ - _description_ + disjunct : Pyomo.Disjunct + Pyomo Disjunct block to include the constraints for the deactivation of the reactor (bypass the reactor). i : float Index of the component in the reactor series. Returns ------- - _type_ - _description_ + Pyomo.Constraint + Deactivation of the recycle flow for each component in the reactor series. """ return m.FR[i, n] == 0 - # Rate desactivation + # Rate deactivation @disjunct.Constraint(m.I) def neg_YPD_rate_desact(disjunct, i): - """_summary_ + """ + Deactivate the reaction rate for each component in the reactor series. + There are no reaction rates when the reactor is deactivated (bypassed). Parameters ---------- - disjunct : _type_ - _description_ + disjunct : Pyomo.Disjunct + Pyomo Disjunct block to include the constraints for the deactivation of the reactor (bypass the reactor). i : float Index of the component in the reactor series. Returns ------- - _type_ - _description_ + Pyomo.Constraint + Deactivation of the reaction rate for each component in the reactor series. """ return m.rate[i, n] == 0 - # QFR desactivation + # QFR deactivation @disjunct.Constraint() def neg_YPD_QFR_desact(disjunct): - """_summary_ + """ + Deactivate the outlet flow rate recycle activation of the reactor. + There is no outlet flow rate recycle activation when the reactor is deactivated (bypassed). Parameters ---------- - disjunct : _type_ - _description_ + disjunct : Pyomo.Disjunct + Pyomo Disjunct block to include the constraints for the deactivation of the reactor (bypass the reactor). Returns ------- - _type_ - _description_ + Pyomo.Constraint + Deactivation of the outlet flow rate recycle activation of the reactor. """ return m.QFR[n] == 0 @disjunct.Constraint() def neg_YPD_vol_desact(disjunct): - ''' - Volume desactivation function for defining pyomo model - args: - disjunct: pyomo block with disjunct to include the constraint - n: pyomo set with reactor index - return: - return constraint - ''' + """ + Volume deactivation function for defining pyomo model. + There is no volume when the reactor is deactivated (bypassed). + + Parameters + ---------- + disjunct : Pyomo.Disjunct + Pyomo Disjunct block to include the constraints for the deactivation of the reactor (bypass the reactor). + + Returns + ------- + Pyomo.Constraint + Volume deactivation function for defining pyomo model. + """ return m.c[n] == 0 # YR Disjuction block equation definition def build_recycle_equations(disjunct, n): - """_summary_ + """ + Build the constraints for the activation of the recycle flow. Parameters ---------- - disjunct : _type_ - _description_ + disjunct : Pyomo.Disjunct + Pyomo Disjunct block to include the constraints for the activation of the reactor (recycle flow existence). n : int Index of the reactor in the reactor series. Returns ------- - _type_ - _description_ + None + None, the function builds the constraints for the activation of the recycle flow. """ m = disjunct.model() # FR activation @disjunct.Constraint(m.I) def YRD_FR_act(disjunct, i): - """_summary_ + """ + Activation of the recycle flow for each component in the reactor series. Parameters ---------- - disjunct : _type_ - _description_ + disjunct : Pyomo.Disjunct + Pyomo Disjunct block to include the constraints for the activation of the reactor (recycle flow existence). i : float Index of the component in the reactor series. Returns ------- - _type_ - _description_ + Pyomo.Constraint + Activation of the recycle flow for each component in the reactor series. """ return m.FR[i, n] - m.R[i] == 0 # QFR activation @disjunct.Constraint() def YRD_QFR_act(disjunct): - """_summary_ + """ + Activation of the outlet flow rate recycle activation of the reactor. Parameters ---------- - disjunct : _type_ - _description_ + disjunct : Pyomo.Disjunct + Pyomo Disjunct block to include the constraints for the activation of the reactor (recycle flow existence). Returns ------- - _type_ - _description_ + Pyomo.Constraint + Activation of the outlet flow rate recycle activation of the reactor. """ return m.QFR[n] - m.QR == 0 def build_no_recycle_equations(disjunct, n): - """_summary_ + """ + Build the constraints for the deactivation of the recycle flow. Parameters ---------- - disjunct : _type_ - _description_ + disjunct : Pyomo.Disjunct + Pyomo Disjunct block to include the constraints for the deactivation of the reactor (recycle flow absence). n : int Index of the reactor in the reactor series. Returns ------- - _type_ - _description_ + None + None, the function builds the constraints for the deactivation of the recycle flow. """ m = disjunct.model() - # FR desactivation + # FR deactivation @disjunct.Constraint(m.I) def neg_YRD_FR_desact(disjunct, i): - """_summary_ + """ + Deactivation of the recycle flow for each component in the reactor series. Parameters ---------- - disjunct : _type_ - _description_ + disjunct : Pyomo.Disjunct + Pyomo Disjunct block to include the constraints for the deactivation of the reactor (recycle flow absence). i : float Index of the component in the reactor series. Returns ------- - _type_ - _description_ + Pyomo.Constraint + Deactivation of the recycle flow for each component in the reactor series. """ return m.FR[i, n] == 0 - # QFR desactivation + # QFR deactivation @disjunct.Constraint() def neg_YRD_QFR_desact(disjunct): - """_summary_ + """ + Deactivation of the outlet flow rate recycle activation of the reactor. Parameters ---------- - disjunct : _type_ - _description_ + disjunct : Pyomo.Disjunct + Pyomo Disjunct block to include the constraints for the deactivation of the reactor (recycle flow absence). Returns ------- - _type_ - _description_ + Pyomo.Constraint + Deactivation of the outlet flow rate recycle activation of the reactor. """ return m.QFR[n] == 0 # Create disjunction blocks - m.YR_is_recycle = Disjunct(m.N, rule=build_recycle_equations) - m.YR_is_not_recycle = Disjunct(m.N, rule=build_no_recycle_equations) + m.YR_is_recycle = Disjunct(m.N, rule=build_recycle_equations, doc="Recycle flow in reactor n") + m.YR_is_not_recycle = Disjunct(m.N, rule=build_no_recycle_equations, doc="No recycle flow in reactor n") - m.YP_is_cstr = Disjunct(m.N, rule=build_cstr_equations) - m.YP_is_bypass = Disjunct(m.N, rule=build_bypass_equations) + m.YP_is_cstr = Disjunct(m.N, rule=build_cstr_equations, doc="CSTR reactor n") + m.YP_is_bypass = Disjunct(m.N, rule=build_bypass_equations, doc="Bypass reactor n") # Create disjunctions @m.Disjunction(m.N) def YP_is_cstr_or_bypass(m, n): - """_summary_ + """ + Build the disjunction for the activation of the CSTR reactor or bypass the reactor. Parameters ---------- @@ -626,13 +654,14 @@ def YP_is_cstr_or_bypass(m, n): Returns ------- list - _description_ + list of the disjunctions for the activation of the CSTR reactor or bypass the reactor """ return [m.YP_is_cstr[n], m.YP_is_bypass[n]] @m.Disjunction(m.N) def YR_is_recycle_or_not(m, n): - """_summary_ + """ + Build the disjunction for the existence of a recycle flow in the reactor. Parameters ---------- @@ -644,7 +673,7 @@ def YR_is_recycle_or_not(m, n): Returns ------- list - _description_ + list of the disjunctions for the existence of a recycle flow in the reactor """ return [m.YR_is_recycle[n], m.YR_is_not_recycle[n]] @@ -657,7 +686,9 @@ def YR_is_recycle_or_not(m, n): # Unit must be a CSTR to include a recycle def cstr_if_recycle_rule(m, n): - """_summary_ + """ + Build the logical constraint for the unit to be a CSTR to include a recycle. + The existence of a recycle flow implies the existence of a CSTR reactor. Parameters ---------- @@ -668,17 +699,18 @@ def cstr_if_recycle_rule(m, n): Returns ------- - _type_ - _description_ + Pyomo.LogicalConstraint + Logical constraint for the unit to be a CSTR to include a recycle. """ return m.YR[n].implies(m.YP[n]) - m.cstr_if_recycle = pyo.LogicalConstraint(m.N, rule=cstr_if_recycle_rule) + m.cstr_if_recycle = pyo.LogicalConstraint(m.N, rule=cstr_if_recycle_rule, doc="Unit must be a CSTR to include a recycle") # There is only one unreacted feed def one_unreacted_feed_rule(m): - """_summary_ + """ + Build the logical constraint for the existence of only one unreacted feed. Parameters ---------- @@ -688,16 +720,17 @@ def one_unreacted_feed_rule(m): Returns ------- Pyomo.LogicalConstraint - _description_ + Logical constraint for the existence of only one unreacted feed. """ return pyo.exactly(1, m.YF) - m.one_unreacted_feed = pyo.LogicalConstraint(rule=one_unreacted_feed_rule) + m.one_unreacted_feed = pyo.LogicalConstraint(rule=one_unreacted_feed_rule, doc="There is only one unreacted feed") # There is only one recycle stream def one_recycle_rule(m): - """_summary_ + """ + Build the logical constraint for the existence of only one recycle stream. Parameters ---------- @@ -707,16 +740,20 @@ def one_recycle_rule(m): Returns ------- Pyomo.LogicalConstraint - _description_ + Logical constraint for the existence of only one recycle stream. """ return pyo.exactly(1, m.YR) - m.one_recycle = pyo.LogicalConstraint(rule=one_recycle_rule) + m.one_recycle = pyo.LogicalConstraint(rule=one_recycle_rule, doc="There is only one recycle stream") # Unit operation in n constraint def unit_in_n_rule(m, n): - """_summary_ + """ + Build the logical constraint for the unit operation in n. + If n is equal to 1, the unit operation is a CSTR. + Otherwise, the unit operation for reactor n except reactor 1 is equivalent to the logical OR of the negation of the unreacted feed of the previous reactors and the unreacted feed of reactor n. + Reactor n is active if either the previous reactors (1 through n-1) have no unreacted feed or reactor n has unreacted feed. Parameters ---------- @@ -728,7 +765,7 @@ def unit_in_n_rule(m, n): Returns ------- Pyomo.LogicalConstraint - _description_ + Logical constraint for the unit operation in n. """ if n == 1: return m.YP[n].equivalent_to(True) From 4d3a013a71e7b1c4b68f20973f853c22e103ceed Mon Sep 17 00:00:00 2001 From: Albert Lee Date: Sun, 19 May 2024 13:47:06 -0400 Subject: [PATCH 08/15] Replaced the Spanglish function definition into English definition. --- gdplib/cstr/gdp_reactor.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/gdplib/cstr/gdp_reactor.py b/gdplib/cstr/gdp_reactor.py index b88475d..09209c9 100644 --- a/gdplib/cstr/gdp_reactor.py +++ b/gdplib/cstr/gdp_reactor.py @@ -439,7 +439,7 @@ def build_bypass_equations(disjunct, n): # FR deactivation @disjunct.Constraint(m.I) - def neg_YPD_FR_desact(disjunct, i): + def neg_YPD_FR_deactivation(disjunct, i): """ Deactivation of the recycle flow for each component in the reactor series. There are no recycle flows when the reactor is deactivated (bypassed). @@ -460,7 +460,7 @@ def neg_YPD_FR_desact(disjunct, i): # Rate deactivation @disjunct.Constraint(m.I) - def neg_YPD_rate_desact(disjunct, i): + def neg_YPD_rate_deactivation(disjunct, i): """ Deactivate the reaction rate for each component in the reactor series. There are no reaction rates when the reactor is deactivated (bypassed). @@ -481,7 +481,7 @@ def neg_YPD_rate_desact(disjunct, i): # QFR deactivation @disjunct.Constraint() - def neg_YPD_QFR_desact(disjunct): + def neg_YPD_QFR_deactivation(disjunct): """ Deactivate the outlet flow rate recycle activation of the reactor. There is no outlet flow rate recycle activation when the reactor is deactivated (bypassed). @@ -499,7 +499,7 @@ def neg_YPD_QFR_desact(disjunct): return m.QFR[n] == 0 @disjunct.Constraint() - def neg_YPD_vol_desact(disjunct): + def neg_YPD_vol_deactivation(disjunct): """ Volume deactivation function for defining pyomo model. There is no volume when the reactor is deactivated (bypassed). @@ -594,7 +594,7 @@ def build_no_recycle_equations(disjunct, n): # FR deactivation @disjunct.Constraint(m.I) - def neg_YRD_FR_desact(disjunct, i): + def neg_YRD_FR_deactivation(disjunct, i): """ Deactivation of the recycle flow for each component in the reactor series. @@ -614,7 +614,7 @@ def neg_YRD_FR_desact(disjunct, i): # QFR deactivation @disjunct.Constraint() - def neg_YRD_QFR_desact(disjunct): + def neg_YRD_QFR_deactivation(disjunct): """ Deactivation of the outlet flow rate recycle activation of the reactor. From 1734d3ab62d2aa897f6de64d057482d26820a102 Mon Sep 17 00:00:00 2001 From: Albert Lee Date: Sun, 19 May 2024 13:48:03 -0400 Subject: [PATCH 09/15] Black format --- gdplib/cstr/gdp_reactor.py | 206 ++++++++++++++++++++++++++++--------- 1 file changed, 157 insertions(+), 49 deletions(-) diff --git a/gdplib/cstr/gdp_reactor.py b/gdplib/cstr/gdp_reactor.py index 09209c9..3d3700b 100644 --- a/gdplib/cstr/gdp_reactor.py +++ b/gdplib/cstr/gdp_reactor.py @@ -10,7 +10,7 @@ def build_cstrs(NT: int = 5) -> pyo.ConcreteModel(): """ Build the CSTR superstructure model of size NT. NT is the number of reactors in series. - The CSTRs have a single 1st order auto catalytic reaction A -> B and minimizes total reactors series volume. + The CSTRs have a single 1st order auto catalytic reaction A -> B and minimizes total reactors series volume. The optimal solution should yield NT reactors with a recycle before reactor NT. Parameters @@ -32,20 +32,30 @@ def build_cstrs(NT: int = 5) -> pyo.ConcreteModel(): m.N = pyo.RangeSet(1, NT, doc='Set of units in the superstructure') # PARAMETERS - m.k = pyo.Param(initialize=2, doc="Kinetic constant [L/(mol*s)]") # Kinetic constant [L/(mol*s)] - m.order1 = pyo.Param(initialize=1, doc="Partial order of reaction 1") # Partial order of reacton 1 - m.order2 = pyo.Param(initialize=1, doc="Partial order of reaction 2") # Partial order of reaction 2 - m.QF0 = pyo.Param(initialize=1, doc="Inlet volumetric flow [L/s]") # Inlet volumetric flow [L/s] + m.k = pyo.Param( + initialize=2, doc="Kinetic constant [L/(mol*s)]" + ) # Kinetic constant [L/(mol*s)] + m.order1 = pyo.Param( + initialize=1, doc="Partial order of reaction 1" + ) # Partial order of reacton 1 + m.order2 = pyo.Param( + initialize=1, doc="Partial order of reaction 2" + ) # Partial order of reaction 2 + m.QF0 = pyo.Param( + initialize=1, doc="Inlet volumetric flow [L/s]" + ) # Inlet volumetric flow [L/s] C0_Def = {'A': 0.99, 'B': 0.01} # Initial concentration of reagents [mol/L] - m.C0 = pyo.Param(m.I, initialize=C0_Def, doc="Initial concentration of reagents [mol/L]") + m.C0 = pyo.Param( + m.I, initialize=C0_Def, doc="Initial concentration of reagents [mol/L]" + ) # Inlet molar flow [mol/s] def F0_Def(m, i): """ - Inlet molar flow [mol/s] for component i. + Inlet molar flow [mol/s] for component i. The function multiplies the initial concentration of component i by the inlet volumetric flow. Parameters @@ -60,7 +70,8 @@ def F0_Def(m, i): Pyomo.Param Inlet molar flow [mol/s] for component i """ - return m.C0[i]*m.QF0 + return m.C0[i] * m.QF0 + m.F0 = pyo.Param(m.I, initialize=F0_Def) # BOOLEAN VARIABLES @@ -78,41 +89,105 @@ def F0_Def(m, i): # Network Variables # Outlet flow rate of the superstructure unit [L/s] - m.Q = pyo.Var(m.N, initialize=0, within=pyo.NonNegativeReals, bounds=(0, 10), doc="Outlet flow rate of the superstructure unit [L/s]") + m.Q = pyo.Var( + m.N, + initialize=0, + within=pyo.NonNegativeReals, + bounds=(0, 10), + doc="Outlet flow rate of the superstructure unit [L/s]", + ) # Outlet flow rate recycle activation of the superstructure unit [L/s] - m.QFR = pyo.Var(m.N, initialize=0, - within=pyo.NonNegativeReals, bounds=(0, 10), doc="Outlet flow rate recycle activation of the superstructure unit [L/s]") + m.QFR = pyo.Var( + m.N, + initialize=0, + within=pyo.NonNegativeReals, + bounds=(0, 10), + doc="Outlet flow rate recycle activation of the superstructure unit [L/s]", + ) # Molar flow [mol/s] - m.F = pyo.Var(m.I, m.N, initialize=0, - within=pyo.NonNegativeReals, bounds=(0, 10), doc="Molar flow [mol/s]") + m.F = pyo.Var( + m.I, + m.N, + initialize=0, + within=pyo.NonNegativeReals, + bounds=(0, 10), + doc="Molar flow [mol/s]", + ) # Molar flow recycle activation [mol/s] - m.FR = pyo.Var(m.I, m.N, initialize=0, - within=pyo.NonNegativeReals, bounds=(0, 10), doc="Molar flow recycle activation [mol/s]") + m.FR = pyo.Var( + m.I, + m.N, + initialize=0, + within=pyo.NonNegativeReals, + bounds=(0, 10), + doc="Molar flow recycle activation [mol/s]", + ) # Reaction rate [mol/(L*s)] - m.rate = pyo.Var(m.I, m.N, initialize=0, within=pyo.Reals, bounds=(-10, 10), doc="Reaction rate [mol/(L*s)]") + m.rate = pyo.Var( + m.I, + m.N, + initialize=0, + within=pyo.Reals, + bounds=(-10, 10), + doc="Reaction rate [mol/(L*s)]", + ) # Reactor volume [L] - m.V = pyo.Var(m.N, initialize=0, within=pyo.NonNegativeReals, bounds=(0, 10), doc="Reactor volume [L]") + m.V = pyo.Var( + m.N, + initialize=0, + within=pyo.NonNegativeReals, + bounds=(0, 10), + doc="Reactor volume [L]", + ) # Volume activation [L] - m.c = pyo.Var(m.N, initialize=0, within=pyo.NonNegativeReals, bounds=(0, 10), doc="Volume activation [L]") + m.c = pyo.Var( + m.N, + initialize=0, + within=pyo.NonNegativeReals, + bounds=(0, 10), + doc="Volume activation [L]", + ) # Splitter Variables # Recycle flow rate [L/s] - m.QR = pyo.Var(initialize=0, within=pyo.NonNegativeReals, bounds=(0, 10), doc="Recycle flow rate [L/s]") + m.QR = pyo.Var( + initialize=0, + within=pyo.NonNegativeReals, + bounds=(0, 10), + doc="Recycle flow rate [L/s]", + ) # Product flow rate [L/s] - m.QP = pyo.Var(initialize=0, within=pyo.NonNegativeReals, bounds=(0, 10), doc="Product flow rate [L/s]") + m.QP = pyo.Var( + initialize=0, + within=pyo.NonNegativeReals, + bounds=(0, 10), + doc="Product flow rate [L/s]", + ) # Recycle molar flow [mol/s] - m.R = pyo.Var(m.I, initialize=0, within=pyo.NonNegativeReals, bounds=(0, 10), doc="Recycle molar flow [mol/s]") + m.R = pyo.Var( + m.I, + initialize=0, + within=pyo.NonNegativeReals, + bounds=(0, 10), + doc="Recycle molar flow [mol/s]", + ) # Product molar flow [mol/s] - m.P = pyo.Var(m.I, initialize=0, within=pyo.NonNegativeReals, bounds=(0, 10), doc="Product molar flow [mol/s]") + m.P = pyo.Var( + m.I, + initialize=0, + within=pyo.NonNegativeReals, + bounds=(0, 10), + doc="Product molar flow [mol/s]", + ) # CONSTRAINTS @@ -139,11 +214,13 @@ def unreact_mole_rule(m, i, n): The constraint for the unreacted feed unit mole balance if n is equal to NT. Otherwise, it returns Pyomo.Constraint.Skip. """ if n == NT: - return m.F0[i] + m.FR[i, n] - m.F[i, n] + m.rate[i, n]*m.V[n] == 0 + return m.F0[i] + m.FR[i, n] - m.F[i, n] + m.rate[i, n] * m.V[n] == 0 else: return pyo.Constraint.Skip - m.unreact_mole = pyo.Constraint(m.I, m.N, rule=unreact_mole_rule, doc="Unreacted feed unit mole balance") + m.unreact_mole = pyo.Constraint( + m.I, m.N, rule=unreact_mole_rule, doc="Unreacted feed unit mole balance" + ) # Unreacted feed unit continuity @@ -169,7 +246,9 @@ def unreact_cont_rule(m, n): else: return pyo.Constraint.Skip - m.unreact_cont = pyo.Constraint(m.N, rule=unreact_cont_rule, doc="Unreacted feed unit continuity") + m.unreact_cont = pyo.Constraint( + m.N, rule=unreact_cont_rule, doc="Unreacted feed unit continuity" + ) # Reactor Balances # Reactor mole balance @@ -194,11 +273,13 @@ def react_mole_rule(m, i, n): The constraint for the reactor mole balance if n is different from NT. Otherwise, it returns Pyomo.Constraint.Skip. """ if n != NT: - return m.F[i, n+1] + m.FR[i, n] - m.F[i, n] + m.rate[i, n]*m.V[n] == 0 + return m.F[i, n + 1] + m.FR[i, n] - m.F[i, n] + m.rate[i, n] * m.V[n] == 0 else: return pyo.Constraint.Skip - m.react_mole = pyo.Constraint(m.I, m.N, rule=react_mole_rule, doc="Reactor mole balance") + m.react_mole = pyo.Constraint( + m.I, m.N, rule=react_mole_rule, doc="Reactor mole balance" + ) # Reactor continuity @@ -220,7 +301,7 @@ def react_cont_rule(m, n): The constraint for the reactor continuity if n is different from NT. Otherwise, it returns Pyomo.Constraint.Skip. """ if n != NT: - return m.Q[n+1] + m.QFR[n] - m.Q[n] == 0 + return m.Q[n + 1] + m.QFR[n] - m.Q[n] == 0 else: return pyo.Constraint.Skip @@ -248,7 +329,9 @@ def split_mole_rule(m, i): """ return m.F[i, 1] - m.P[i] - m.R[i] == 0 - m.split_mole = pyo.Constraint(m.I, rule=split_mole_rule, doc="Splitting point mole balance") + m.split_mole = pyo.Constraint( + m.I, rule=split_mole_rule, doc="Splitting point mole balance" + ) # Splitting point continuity @@ -269,13 +352,15 @@ def split_cont_rule(m): """ return m.Q[1] - m.QP - m.QR == 0 - m.split_cont = pyo.Constraint(rule=split_cont_rule, doc="Splitting point continuity") + m.split_cont = pyo.Constraint( + rule=split_cont_rule, doc="Splitting point continuity" + ) # Splitting point additional constraints def split_add_rule(m, i): """ - Splitting point additional constraints. + Splitting point additional constraints. Molarity constraints over initial and final flows, read as an multiplication avoid the numerical complication. m.P[i]/m.QP = m.F[i,1]/m.Q[1] (molarity balance) @@ -291,9 +376,11 @@ def split_add_rule(m, i): Pyomo.Constraint The constraint for the splitting point additional constraints. """ - return m.P[i]*m.Q[1] - m.F[i, 1]*m.QP == 0 + return m.P[i] * m.Q[1] - m.F[i, 1] * m.QP == 0 - m.split_add = pyo.Constraint(m.I, rule=split_add_rule, doc="Splitting point additional constraints") + m.split_add = pyo.Constraint( + m.I, rule=split_add_rule, doc="Splitting point additional constraints" + ) # Product Specification @@ -312,7 +399,7 @@ def prod_spec_rule(m): Pyomo.Constraint The constraint for the product B specification. """ - return m.QP*0.95 - m.P['B'] == 0 + return m.QP * 0.95 - m.P['B'] == 0 m.prod_spec = pyo.Constraint(rule=prod_spec_rule, doc="Product B Specification") @@ -335,7 +422,7 @@ def vol_cons_rule(m, n): The constraint for the volume constraint if n is different from 1. Otherwise, it returns Pyomo.Constraint.Skip. """ if n != 1: - return m.V[n] - m.V[n-1] == 0 + return m.V[n] - m.V[n - 1] == 0 else: return pyo.Constraint.Skip @@ -378,7 +465,11 @@ def YPD_rate_calc(disjunct): Pyomo.Constraint The constraint for the calculation of the reaction rate of A for each reactor. """ - return m.rate['A', n]*((m.Q[n])**m.order1)*((m.Q[n])**m.order2)+m.k*((m.F['A', n])**m.order1)*((m.F['B', n])**m.order2) == 0 + return ( + m.rate['A', n] * ((m.Q[n]) ** m.order1) * ((m.Q[n]) ** m.order2) + + m.k * ((m.F['A', n]) ** m.order1) * ((m.F['B', n]) ** m.order2) + == 0 + ) # Reaction rate relation @disjunct.Constraint() @@ -387,7 +478,7 @@ def YPD_rate_rel(disjunct): Reaction rate relation for defining pyomo model. Since the chemical reaction goes from A to B, the rate of A is equal to the negative rate of B. - A -> B + A -> B Parameters ---------- @@ -508,7 +599,7 @@ def neg_YPD_vol_deactivation(disjunct): ---------- disjunct : Pyomo.Disjunct Pyomo Disjunct block to include the constraints for the deactivation of the reactor (bypass the reactor). - + Returns ------- Pyomo.Constraint @@ -604,7 +695,7 @@ def neg_YRD_FR_deactivation(disjunct, i): Pyomo Disjunct block to include the constraints for the deactivation of the reactor (recycle flow absence). i : float Index of the component in the reactor series. - + Returns ------- Pyomo.Constraint @@ -631,8 +722,12 @@ def neg_YRD_QFR_deactivation(disjunct): return m.QFR[n] == 0 # Create disjunction blocks - m.YR_is_recycle = Disjunct(m.N, rule=build_recycle_equations, doc="Recycle flow in reactor n") - m.YR_is_not_recycle = Disjunct(m.N, rule=build_no_recycle_equations, doc="No recycle flow in reactor n") + m.YR_is_recycle = Disjunct( + m.N, rule=build_recycle_equations, doc="Recycle flow in reactor n" + ) + m.YR_is_not_recycle = Disjunct( + m.N, rule=build_no_recycle_equations, doc="No recycle flow in reactor n" + ) m.YP_is_cstr = Disjunct(m.N, rule=build_cstr_equations, doc="CSTR reactor n") m.YP_is_bypass = Disjunct(m.N, rule=build_bypass_equations, doc="Bypass reactor n") @@ -704,7 +799,9 @@ def cstr_if_recycle_rule(m, n): """ return m.YR[n].implies(m.YP[n]) - m.cstr_if_recycle = pyo.LogicalConstraint(m.N, rule=cstr_if_recycle_rule, doc="Unit must be a CSTR to include a recycle") + m.cstr_if_recycle = pyo.LogicalConstraint( + m.N, rule=cstr_if_recycle_rule, doc="Unit must be a CSTR to include a recycle" + ) # There is only one unreacted feed @@ -724,7 +821,9 @@ def one_unreacted_feed_rule(m): """ return pyo.exactly(1, m.YF) - m.one_unreacted_feed = pyo.LogicalConstraint(rule=one_unreacted_feed_rule, doc="There is only one unreacted feed") + m.one_unreacted_feed = pyo.LogicalConstraint( + rule=one_unreacted_feed_rule, doc="There is only one unreacted feed" + ) # There is only one recycle stream @@ -744,14 +843,16 @@ def one_recycle_rule(m): """ return pyo.exactly(1, m.YR) - m.one_recycle = pyo.LogicalConstraint(rule=one_recycle_rule, doc="There is only one recycle stream") + m.one_recycle = pyo.LogicalConstraint( + rule=one_recycle_rule, doc="There is only one recycle stream" + ) # Unit operation in n constraint def unit_in_n_rule(m, n): """ Build the logical constraint for the unit operation in n. - If n is equal to 1, the unit operation is a CSTR. + If n is equal to 1, the unit operation is a CSTR. Otherwise, the unit operation for reactor n except reactor 1 is equivalent to the logical OR of the negation of the unreacted feed of the previous reactors and the unreacted feed of reactor n. Reactor n is active if either the previous reactors (1 through n-1) have no unreacted feed or reactor n has unreacted feed. @@ -770,7 +871,9 @@ def unit_in_n_rule(m, n): if n == 1: return m.YP[n].equivalent_to(True) else: - return m.YP[n].equivalent_to(pyo.lor(pyo.land(~m.YF[n2] for n2 in range(1, n)), m.YF[n])) + return m.YP[n].equivalent_to( + pyo.lor(pyo.land(~m.YF[n2] for n2 in range(1, n)), m.YF[n]) + ) m.unit_in_n = pyo.LogicalConstraint(m.N, rule=unit_in_n_rule) @@ -792,13 +895,18 @@ def obj_rule(m): """ return sum(m.c[n] for n in m.N) - m.obj = pyo.Objective(rule=obj_rule, sense=pyo.minimize, doc="minimum total reactor volume") + m.obj = pyo.Objective( + rule=obj_rule, sense=pyo.minimize, doc="minimum total reactor volume" + ) return m + if __name__ == "__main__": m = build_cstrs() 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) \ No newline at end of file + pyo.SolverFactory('gams').solve( + m, solver='baron', tee=True, add_options=['option optcr=1e-6;'] + ) + display(m) From 87250364e7c5ad9f03a774f1a3beef2b107f2363 Mon Sep 17 00:00:00 2001 From: Albert Lee Date: Sun, 19 May 2024 15:13:51 -0400 Subject: [PATCH 10/15] black formatted again --- gdplib/cstr/gdp_reactor.py | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/gdplib/cstr/gdp_reactor.py b/gdplib/cstr/gdp_reactor.py index 3d3700b..c010a7c 100644 --- a/gdplib/cstr/gdp_reactor.py +++ b/gdplib/cstr/gdp_reactor.py @@ -25,11 +25,11 @@ def build_cstrs(NT: int = 5) -> pyo.ConcreteModel(): """ # PYOMO MODEL - m = pyo.ConcreteModel(name='gdp_reactors') + m = pyo.ConcreteModel(name="gdp_reactors") # SETS - m.I = pyo.Set(initialize=['A', 'B'], doc='Set of components') - m.N = pyo.RangeSet(1, NT, doc='Set of units in the superstructure') + m.I = pyo.Set(initialize=["A", "B"], doc="Set of components") + m.N = pyo.RangeSet(1, NT, doc="Set of units in the superstructure") # PARAMETERS m.k = pyo.Param( @@ -44,7 +44,7 @@ def build_cstrs(NT: int = 5) -> pyo.ConcreteModel(): m.QF0 = pyo.Param( initialize=1, doc="Inlet volumetric flow [L/s]" ) # Inlet volumetric flow [L/s] - C0_Def = {'A': 0.99, 'B': 0.01} + C0_Def = {"A": 0.99, "B": 0.01} # Initial concentration of reagents [mol/L] m.C0 = pyo.Param( @@ -399,7 +399,7 @@ def prod_spec_rule(m): Pyomo.Constraint The constraint for the product B specification. """ - return m.QP * 0.95 - m.P['B'] == 0 + return m.QP * 0.95 - m.P["B"] == 0 m.prod_spec = pyo.Constraint(rule=prod_spec_rule, doc="Product B Specification") @@ -466,8 +466,8 @@ def YPD_rate_calc(disjunct): The constraint for the calculation of the reaction rate of A for each reactor. """ return ( - m.rate['A', n] * ((m.Q[n]) ** m.order1) * ((m.Q[n]) ** m.order2) - + m.k * ((m.F['A', n]) ** m.order1) * ((m.F['B', n]) ** m.order2) + m.rate["A", n] * ((m.Q[n]) ** m.order1) * ((m.Q[n]) ** m.order2) + + m.k * ((m.F["A", n]) ** m.order1) * ((m.F["B", n]) ** m.order2) == 0 ) @@ -490,7 +490,7 @@ def YPD_rate_rel(disjunct): Pyomo.Constraint _description_ """ - return m.rate['B', n] + m.rate['A', n] == 0 + return m.rate["B", n] + m.rate["A", n] == 0 # Volume activation @disjunct.Constraint() @@ -904,9 +904,9 @@ def obj_rule(m): if __name__ == "__main__": m = build_cstrs() - 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 756b4eff1ecfa63cf7dd7f21e1504698cb3728b3 Mon Sep 17 00:00:00 2001 From: Albert Lee Date: Sun, 19 May 2024 17:12:11 -0400 Subject: [PATCH 11/15] Add the documentation build_cstr_equation. --- gdplib/cstr/gdp_reactor.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/gdplib/cstr/gdp_reactor.py b/gdplib/cstr/gdp_reactor.py index c010a7c..27fb387 100644 --- a/gdplib/cstr/gdp_reactor.py +++ b/gdplib/cstr/gdp_reactor.py @@ -431,7 +431,8 @@ def vol_cons_rule(m, n): # YD Disjunction block equation definition def build_cstr_equations(disjunct, n): - """_summary_ + """ + Build the constraints for the activation of the CSTR reactor. Parameters ---------- From 995790d6db4c5b3d1012b270cfc264e89e80eb6d Mon Sep 17 00:00:00 2001 From: Albert Lee Date: Sun, 19 May 2024 17:14:19 -0400 Subject: [PATCH 12/15] Added the reaction rate relation documentation. --- gdplib/cstr/gdp_reactor.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gdplib/cstr/gdp_reactor.py b/gdplib/cstr/gdp_reactor.py index 27fb387..036fa4f 100644 --- a/gdplib/cstr/gdp_reactor.py +++ b/gdplib/cstr/gdp_reactor.py @@ -489,7 +489,7 @@ def YPD_rate_rel(disjunct): Returns ------- Pyomo.Constraint - _description_ + Reaction rate relation for defining pyomo model. """ return m.rate["B", n] + m.rate["A", n] == 0 From 33685ebfed947b552c205840a8bdbbeb98a021a5 Mon Sep 17 00:00:00 2001 From: Albert Lee Date: Sun, 19 May 2024 17:33:59 -0400 Subject: [PATCH 13/15] Added the header of the code. --- gdplib/cstr/gdp_reactor.py | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/gdplib/cstr/gdp_reactor.py b/gdplib/cstr/gdp_reactor.py index 036fa4f..7473f0d 100644 --- a/gdplib/cstr/gdp_reactor.py +++ b/gdplib/cstr/gdp_reactor.py @@ -4,7 +4,16 @@ from pyomo.core.base.misc import display from pyomo.gdp import Disjunct, Disjunction from pyomo.opt.base.solvers import SolverFactory +""" +gdp_reactor.py +This script builds a Pyomo GDP model for a CSTR superstructure with a single 1st order auto catalytic reaction A -> B. +The model minimizes the total reactors series volume. +The optimal solution should yield NT reactors with a recycle before reactor NT. + +Reference: + > Linan, D. A., Bernal, D. E., Gomez, J. M., & Ricardez-Sandoval, L. A. (2021). Optimal synthesis and design of catalytic distillation columns: A rate-based modeling approach. Chemical Engineering Science, 231, 116294. https://doi.org/10.1016/j.ces.2020.116294 +""" def build_cstrs(NT: int = 5) -> pyo.ConcreteModel(): """ From eecbcd12afbca382bbc52579f3998d2aab86eeb7 Mon Sep 17 00:00:00 2001 From: Albert Lee Date: Sun, 19 May 2024 17:35:01 -0400 Subject: [PATCH 14/15] black format --- gdplib/cstr/gdp_reactor.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/gdplib/cstr/gdp_reactor.py b/gdplib/cstr/gdp_reactor.py index 7473f0d..09cd536 100644 --- a/gdplib/cstr/gdp_reactor.py +++ b/gdplib/cstr/gdp_reactor.py @@ -4,6 +4,7 @@ from pyomo.core.base.misc import display from pyomo.gdp import Disjunct, Disjunction from pyomo.opt.base.solvers import SolverFactory + """ gdp_reactor.py @@ -15,6 +16,7 @@ > Linan, D. A., Bernal, D. E., Gomez, J. M., & Ricardez-Sandoval, L. A. (2021). Optimal synthesis and design of catalytic distillation columns: A rate-based modeling approach. Chemical Engineering Science, 231, 116294. https://doi.org/10.1016/j.ces.2020.116294 """ + def build_cstrs(NT: int = 5) -> pyo.ConcreteModel(): """ Build the CSTR superstructure model of size NT. From 0ade99df270015be2dc2733dc713e6745fa15133 Mon Sep 17 00:00:00 2001 From: Albert Lee Date: Sun, 19 May 2024 20:12:56 -0400 Subject: [PATCH 15/15] relocate the header. --- gdplib/cstr/gdp_reactor.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/gdplib/cstr/gdp_reactor.py b/gdplib/cstr/gdp_reactor.py index 09cd536..b357384 100644 --- a/gdplib/cstr/gdp_reactor.py +++ b/gdplib/cstr/gdp_reactor.py @@ -1,10 +1,3 @@ -import os -import sys -import pyomo.environ as pyo -from pyomo.core.base.misc import display -from pyomo.gdp import Disjunct, Disjunction -from pyomo.opt.base.solvers import SolverFactory - """ gdp_reactor.py @@ -16,6 +9,13 @@ > Linan, D. A., Bernal, D. E., Gomez, J. M., & Ricardez-Sandoval, L. A. (2021). Optimal synthesis and design of catalytic distillation columns: A rate-based modeling approach. Chemical Engineering Science, 231, 116294. https://doi.org/10.1016/j.ces.2020.116294 """ +import os +import sys +import pyomo.environ as pyo +from pyomo.core.base.misc import display +from pyomo.gdp import Disjunct, Disjunction +from pyomo.opt.base.solvers import SolverFactory + def build_cstrs(NT: int = 5) -> pyo.ConcreteModel(): """