From 283830467c93bfcb865d13a2cf9e2f4686b54e70 Mon Sep 17 00:00:00 2001 From: Carolina Tristan Date: Tue, 13 Aug 2024 12:13:15 -0400 Subject: [PATCH 1/6] Revised vapor composition constraints fix ZeroDivisionError when gdp.hull reformulation is applied --- gdplib/kaibel/kaibel_solve_gdp.py | 164 ++++++++++++++++++++++++------ 1 file changed, 132 insertions(+), 32 deletions(-) diff --git a/gdplib/kaibel/kaibel_solve_gdp.py b/gdplib/kaibel/kaibel_solve_gdp.py index 7aca6f6e..5231ccbc 100644 --- a/gdplib/kaibel/kaibel_solve_gdp.py +++ b/gdplib/kaibel/kaibel_solve_gdp.py @@ -26,8 +26,6 @@ from gdplib.kaibel.kaibel_side_flash import calc_side_feed_flash -# from .kaibel_side_flash import calc_side_feed_flash - def build_model(): """ @@ -185,8 +183,8 @@ def build_model(): kc = {} # Light and heavy key components Tbounds[1] = m.Tcon # Condenser temperature [K] Tbounds[2] = m.Treb # Reboiler temperature [K] - kc[1] = m.lc - kc[2] = m.hc + kc[1] = m.lc # Light key component + kc[2] = m.hc # Heavy key component ## Heat of vaporization calculation for each component in the feed. for comp in m.comp: @@ -245,7 +243,7 @@ def build_model(): for comp in m.comp: m.dHvap[comp] = dHvapb[comp] / m.Hscale - ## Heat capacity calculation for liquid and vapor phases using Ruczika-D method for each component in the feed, section, and tray + ## Heat capacity calculation for liquid and vapor phases for each component in the feed, section, and tray for sec in m.section: for n_tray in m.tray: for comp in m.comp: @@ -269,7 +267,7 @@ def build_model(): / m.Hscale ) - ## Liquid and vapor enthalpy calculation using Ruczika-D method for each component in the feed, section, and tray + ## Liquid and vapor enthalpy calculation for each component in the feed, section, and tray for sec in m.section: for n_tray in m.tray: for comp in m.comp: @@ -282,7 +280,7 @@ def build_model(): m.hvf = {} # Vapor enthalpy for side feed [J/mol] m.F0 = {} # Side feed flowrate per component [mol/s] - ## Heat capacity in liquid and vapor phases for side feed for each component using Ruczika-D method + ## Heat capacity in liquid and vapor phases for side feed for each component. for comp in m.comp: for cp in m.cplv: m.cpdTf[comp, cp] = ( @@ -302,7 +300,7 @@ def build_model(): / m.Hscale ) - ## Side feed flowrate and liquid and vapor enthalpy calculation using Ruczika-D method for each component in the feed + ## Side feed flowrate and liquid and vapor enthalpy calculation for each component in the feed for comp in m.comp: m.F0[comp] = ( m.xfi[comp] * m.Fi @@ -317,7 +315,7 @@ def build_model(): m.P = Var( m.section, m.tray, - doc="Pressure at each potential tray in bars", + doc="Pressure at each potential tray [bar]", domain=NonNegativeReals, bounds=(m.Pcon, m.Preb), initialize=m.P0, @@ -331,11 +329,24 @@ def build_model(): initialize=m.T0, ) + m.Tr = Var( + m.section, + m.tray, + m.comp, + doc='Temperature term for vapor pressure [-]', + domain=NonNegativeReals, + bounds=lambda m, sec, tray, comp: ( + m.prop[comp, 'TC'] / m.Tup, + m.prop[comp, 'TC'] / m.Tlo, + ), # (0, None), + initialize=lambda m, sec, tray, comp: m.T0[sec, tray] / m.Tlo, + ) + m.x = Var( m.section, m.tray, m.comp, - doc="Liquid composition", + doc="Liquid composition [mol/mol]", domain=NonNegativeReals, bounds=(0, 1), initialize=m.x0, @@ -344,7 +355,7 @@ def build_model(): m.section, m.tray, m.comp, - doc="Vapor composition", + doc="Vapor composition [mol/mol]", domain=NonNegativeReals, bounds=(0, 1), initialize=m.y0, @@ -1094,7 +1105,6 @@ def _logic_proposition_main(m, sec, n_tray): ) else: return Constraint.NoConstraint - # TOCHECK: Update the logic proposition constraint for the main section with the new pyomo.gdp syntax @m.Constraint(m.candidate_trays_feed) def _logic_proposition_feed(m, n_tray): @@ -1125,7 +1135,6 @@ def _logic_proposition_feed(m, n_tray): ) else: return Constraint.NoConstraint - # TODO: Update the logic proposition constraint for the feed section with the new pyomo.gdp syntax @m.Constraint(m.candidate_trays_product) def _logic_proposition_section3(m, n_tray): @@ -1151,7 +1160,6 @@ def _logic_proposition_section3(m, n_tray): ) else: return Constraint.NoConstraint - # TODO: Update the logic proposition constraint for the product section with the new pyomo.gdp syntax @m.Constraint(m.tray) def equality_feed_product_side(m, n_tray): @@ -1175,8 +1183,6 @@ def equality_feed_product_side(m, n_tray): == m.tray_exists[3, n_tray].binary_indicator_var ) - # TODO: Update the equality constraint for the feed and product side trays with the new pyomo.gdp syntax - @m.Constraint() def _existent_minimum_numbertrays(m): """ @@ -1431,6 +1437,29 @@ def bottom_vapor_composition_summation(disj): """ return sum(m.y[1, n_tray, comp] for comp in m.comp) - 1 == m.erry[1, n_tray] + @disj.Constraint( + m.comp, + doc="Temperature term for the bottom section 1 used in the vapor composition equation", + ) + def _bottom_reduced_temperature(disj, comp): + """Calculate the temperature term for the bottom section vapor composition. + + The constraint is used to avoid division by zero in the MINLP Hull reformulation. + + Parameters + ---------- + disj : Disjunct + The disjunct object for the bottom section in the column. + comp : int + The component index. + + Returns + ------- + Constraint + The constraint equation to calculate the temperature term for the bottom section vapor composition. + """ + return m.Tr[1, n_tray, comp] * m.T[1, n_tray] == m.prop[comp, 'TC'] + @disj.Constraint(m.comp, doc="Bottom section 1 vapor composition") def bottom_vapor_composition(disj, comp): """ @@ -1450,15 +1479,16 @@ def bottom_vapor_composition(disj, comp): The equation is derived from the vapor-liquid equilibrium relationship. """ return ( - m.y[1, n_tray, comp] + m.y[1, n_tray, comp] * m.P[1, n_tray] == m.x[1, n_tray, comp] * ( m.actv[1, n_tray, comp] * ( m.prop[comp, 'PC'] * exp( - m.prop[comp, 'TC'] - / m.T[1, n_tray] + m.Tr[1, n_tray, comp] + # m.prop[comp, 'TC'] + # / m.T[1, n_tray] * ( m.prop[comp, 'vpA'] * (1 - m.T[1, n_tray] / m.prop[comp, 'TC']) @@ -1472,7 +1502,7 @@ def bottom_vapor_composition(disj, comp): ) ) ) - / m.P[1, n_tray] + # / m.P[1, n_tray] ) @disj.Constraint(m.comp, doc="Bottom section 1 liquid enthalpy") @@ -1738,6 +1768,28 @@ def feedside_vapor_composition_summation(disj): """ return sum(m.y[2, n_tray, comp] for comp in m.comp) - 1 == m.erry[2, n_tray] + @disj.Constraint( + m.comp, doc="Feed section 2 temperature term for the vapor composition equation" + ) + def _feedside_reduced_temperature(disj, comp): + """Calculate the temperature term for the feedside section vapor composition. + + The constraint is used to avoid division by zero in the MINLP Hull reformulation. + + Parameters + ---------- + disj : Disjunct + The disjunct object for the bottom section in the column. + comp : int + The component index. + + Returns + ------- + Constraint + The constraint equation to calculate the temperature term for the feedside section vapor composition. + """ + return m.Tr[2, n_tray, comp] * m.T[2, n_tray] == m.prop[comp, 'TC'] + @disj.Constraint(m.comp, doc="Feed section 2 vapor composition") def feedside_vapor_composition(disj, comp): """ @@ -1757,15 +1809,16 @@ def feedside_vapor_composition(disj, comp): The equation is derived from the vapor-liquid equilibrium relationship. """ return ( - m.y[2, n_tray, comp] + m.y[2, n_tray, comp] * m.P[2, n_tray] == m.x[2, n_tray, comp] * ( m.actv[2, n_tray, comp] * ( m.prop[comp, 'PC'] * exp( - m.prop[comp, 'TC'] - / m.T[2, n_tray] + m.Tr[2, n_tray, comp] + # m.prop[comp, 'TC'] + # / m.T[2, n_tray] * ( m.prop[comp, 'vpA'] * (1 - m.T[2, n_tray] / m.prop[comp, 'TC']) @@ -1779,7 +1832,7 @@ def feedside_vapor_composition(disj, comp): ) ) ) - / m.P[2, n_tray] + # / m.P[2, n_tray] ) @disj.Constraint(m.comp, doc="Feed section 2 liquid enthalpy") @@ -2049,6 +2102,29 @@ def productside_vapor_composition_summation(disj): """ return sum(m.y[3, n_tray, comp] for comp in m.comp) - 1 == m.erry[3, n_tray] + @disj.Constraint( + m.comp, + doc="Product section 3 temperature term for the vapor composition equation", + ) + def _feedside_reduced_temperature(disj, comp): + """Calculate the temperature term for the product section vapor composition. + + The constraint is used to avoid division by zero in the MINLP Hull reformulation. + + Parameters + ---------- + disj : Disjunct + The disjunct object for the bottom section in the column. + comp : int + The component index. + + Returns + ------- + Constraint + The constraint equation to calculate the temperature term for the product section vapor composition. + """ + return m.Tr[3, n_tray, comp] * m.T[3, n_tray] == m.prop[comp, 'TC'] + @disj.Constraint(m.comp, doc="Product section 3 vapor composition") def productside_vapor_composition(disj, comp): """ @@ -2068,15 +2144,16 @@ def productside_vapor_composition(disj, comp): The equation is derived from the vapor-liquid equilibrium relationship. """ return ( - m.y[3, n_tray, comp] + m.y[3, n_tray, comp] * m.P[3, n_tray] == m.x[3, n_tray, comp] * ( m.actv[3, n_tray, comp] * ( m.prop[comp, 'PC'] * exp( - m.prop[comp, 'TC'] - / m.T[3, n_tray] + m.Tr[3, n_tray, comp] + # m.prop[comp, 'TC'] + # / m.T[3, n_tray] * ( m.prop[comp, 'vpA'] * (1 - m.T[3, n_tray] / m.prop[comp, 'TC']) @@ -2090,7 +2167,7 @@ def productside_vapor_composition(disj, comp): ) ) ) - / m.P[3, n_tray] + # / m.P[3, n_tray] ) @disj.Constraint(m.comp, doc="Product section 3 liquid enthalpy") @@ -2368,6 +2445,28 @@ def top_vapor_composition_summation(disj): """ return sum(m.y[4, n_tray, comp] for comp in m.comp) - 1 == m.erry[4, n_tray] + @disj.Constraint( + m.comp, doc="Top section 4 temperature term for the vapor composition equation" + ) + def _feedside_reduced_temperature(disj, comp): + """Calculate the temperature term for the top section vapor composition. + + The constraint is used to avoid division by zero in the MINLP Hull reformulation. + + Parameters + ---------- + disj : Disjunct + The disjunct object for the bottom section in the column. + comp : int + The component index. + + Returns + ------- + Constraint + The constraint equation to calculate the temperature term for the top section vapor composition. + """ + return m.Tr[4, n_tray, comp] * m.T[4, n_tray] == m.prop[comp, 'TC'] + @disj.Constraint(m.comp, doc="Top scetion 4 vapor composition") def top_vapor_composition(disj, comp): """ @@ -2387,15 +2486,16 @@ def top_vapor_composition(disj, comp): The equation is derived from the vapor-liquid equilibrium relationship. """ return ( - m.y[4, n_tray, comp] + m.y[4, n_tray, comp] * m.P[4, n_tray] == m.x[4, n_tray, comp] * ( m.actv[4, n_tray, comp] * ( m.prop[comp, 'PC'] * exp( - m.prop[comp, 'TC'] - / m.T[4, n_tray] + m.Tr[4, n_tray, comp] + # m.prop[comp, 'TC'] + # / m.T[4, n_tray] * ( m.prop[comp, 'vpA'] * (1 - m.T[4, n_tray] / m.prop[comp, 'TC']) @@ -2409,7 +2509,7 @@ def top_vapor_composition(disj, comp): ) ) ) - / m.P[4, n_tray] + # / m.P[4, n_tray] ) @disj.Constraint(m.comp, doc="Top section 4 liquid enthalpy") From f798d9d793cda2a8f496a4173a3e2ee77af406e4 Mon Sep 17 00:00:00 2001 From: ZedongPeng Date: Tue, 13 Aug 2024 17:58:18 -0400 Subject: [PATCH 2/6] replace gam with gamma --- gdplib/hda/HDA_GDP_gdpopt.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/gdplib/hda/HDA_GDP_gdpopt.py b/gdplib/hda/HDA_GDP_gdpopt.py index cf817b1d..39c8c852 100644 --- a/gdplib/hda/HDA_GDP_gdpopt.py +++ b/gdplib/hda/HDA_GDP_gdpopt.py @@ -20,7 +20,7 @@ def HDA_model(): m.alpha = Param(initialize=0.3665, doc="compressor coefficient") m.compeff = Param(initialize=0.750, doc="compressor efficiency") - m.gam = Param(initialize=1.300, doc="ratio of cp to cv") + m.gamma = Param(initialize=1.300, doc="ratio of cp to cv") m.abseff = Param(initialize=0.333, doc="absorber tray efficiency") m.disteff = Param(initialize=0.5000, doc="column tray efficiency") m.uflow = Param(initialize=50, doc="upper bound - flow logicals") @@ -1312,7 +1312,7 @@ def Compelec(_m, comp_): * m.f[stream] / 60.0 * (1.0 / m.compeff) - * (m.gam / (m.gam - 1.0)) + * (m.gamma / (m.gamma - 1.0)) for (comp1, stream) in m.icomp if comp_ == comp1 ) @@ -1324,7 +1324,7 @@ def Compelec(_m, comp_): def Ratio(_m, comp_): if comp == comp_: - return m.presrat[comp_] ** (m.gam / (m.gam - 1.0)) == sum( + return m.presrat[comp_] ** (m.gamma / (m.gamma - 1.0)) == sum( m.p[stream] for (comp1, stream) in m.ocomp if comp_ == comp1 ) / sum(m.p[stream] for (comp1, stream) in m.icomp if comp1 == comp_) return Constraint.Skip @@ -2248,11 +2248,11 @@ def Valcmb(_m, valve, compon): def Valt(_m, valve): return sum( - m.t[stream] / (m.p[stream] ** ((m.gam - 1.0) / m.gam)) + m.t[stream] / (m.p[stream] ** ((m.gamma - 1.0) / m.gamma)) for (valv, stream) in m.oval if valv == valve ) == sum( - m.t[stream] / (m.p[stream] ** ((m.gam - 1.0) / m.gam)) + m.t[stream] / (m.p[stream] ** ((m.gamma - 1.0) / m.gamma)) for (valv, stream) in m.ival if valv == valve ) From a70395bcc8cbdd43b75da53b5a090478ca130a2a Mon Sep 17 00:00:00 2001 From: "David E. Bernal Neira" Date: Sun, 10 May 2026 15:22:34 -0400 Subject: [PATCH 3/6] Address PR 47 review feedback --- gdplib/hda/HDA_GDP_gdpopt.py | 26 +++- gdplib/kaibel/kaibel_solve_gdp.py | 249 +++++++++++++----------------- tests/test_pr47_regressions.py | 31 ++++ 3 files changed, 158 insertions(+), 148 deletions(-) create mode 100644 tests/test_pr47_regressions.py diff --git a/gdplib/hda/HDA_GDP_gdpopt.py b/gdplib/hda/HDA_GDP_gdpopt.py index 39c8c852..05bea29d 100644 --- a/gdplib/hda/HDA_GDP_gdpopt.py +++ b/gdplib/hda/HDA_GDP_gdpopt.py @@ -20,7 +20,7 @@ def HDA_model(): m.alpha = Param(initialize=0.3665, doc="compressor coefficient") m.compeff = Param(initialize=0.750, doc="compressor efficiency") - m.gamma = Param(initialize=1.300, doc="ratio of cp to cv") + m.heat_capacity_ratio = Param(initialize=1.300, doc="ratio of cp to cv") m.abseff = Param(initialize=0.333, doc="absorber tray efficiency") m.disteff = Param(initialize=0.5000, doc="column tray efficiency") m.uflow = Param(initialize=50, doc="upper bound - flow logicals") @@ -1312,7 +1312,7 @@ def Compelec(_m, comp_): * m.f[stream] / 60.0 * (1.0 / m.compeff) - * (m.gamma / (m.gamma - 1.0)) + * (m.heat_capacity_ratio / (m.heat_capacity_ratio - 1.0)) for (comp1, stream) in m.icomp if comp_ == comp1 ) @@ -1324,9 +1324,13 @@ def Compelec(_m, comp_): def Ratio(_m, comp_): if comp == comp_: - return m.presrat[comp_] ** (m.gamma / (m.gamma - 1.0)) == sum( + return m.presrat[comp_] ** ( + m.heat_capacity_ratio / (m.heat_capacity_ratio - 1.0) + ) == sum( m.p[stream] for (comp1, stream) in m.ocomp if comp_ == comp1 - ) / sum(m.p[stream] for (comp1, stream) in m.icomp if comp1 == comp_) + ) / sum( + m.p[stream] for (comp1, stream) in m.icomp if comp1 == comp_ + ) return Constraint.Skip b.ratio = Constraint([comp], rule=Ratio, doc='pressure ratio (out to in)') @@ -1416,7 +1420,7 @@ def Relvol(_m, dist_): return m.avevlt[dist] == sqrt(divided1 * divided2) return Constraint.Skip - b.relvol = Constraint([dist], rule=Relvol, doc='average relative volatilty') + b.relvol = Constraint([dist], rule=Relvol, doc='average relative volatility') def Undwood(_m, dist_): if dist_ == dist: @@ -2248,11 +2252,19 @@ def Valcmb(_m, valve, compon): def Valt(_m, valve): return sum( - m.t[stream] / (m.p[stream] ** ((m.gamma - 1.0) / m.gamma)) + m.t[stream] + / ( + m.p[stream] + ** ((m.heat_capacity_ratio - 1.0) / m.heat_capacity_ratio) + ) for (valv, stream) in m.oval if valv == valve ) == sum( - m.t[stream] / (m.p[stream] ** ((m.gamma - 1.0) / m.gamma)) + m.t[stream] + / ( + m.p[stream] + ** ((m.heat_capacity_ratio - 1.0) / m.heat_capacity_ratio) + ) for (valv, stream) in m.ival if valv == valve ) diff --git a/gdplib/kaibel/kaibel_solve_gdp.py b/gdplib/kaibel/kaibel_solve_gdp.py index 5231ccbc..5426d50b 100644 --- a/gdplib/kaibel/kaibel_solve_gdp.py +++ b/gdplib/kaibel/kaibel_solve_gdp.py @@ -1,7 +1,7 @@ -""" Kaibel Column model: GDP formulation. +"""Kaibel Column model: GDP formulation. The solution requires the specification of certain parameters, such as the number trays, feed location, etc., and an initialization procedure, which consists of the next three steps: -(i) a preliminary design of the separation considering a sequence of indirect continuous distillation columns (CDCs) to obtain the minimum number of stages with Fenske Equation in the function initialize_kaibel in kaibel_init.py +(i) a preliminary design of the separation considering a sequence of indirect continuous distillation columns (CDCs) to obtain the minimum number of stages with Fenske Equation in the function initialize_kaibel in kaibel_init.py (ii) flash calculation for the feed with the function calc_side_feed_flash in kaibel_side_flash.py (iii) calculation of variable bounds by solving the NLP problem. @@ -101,11 +101,10 @@ def build_model(): m.so = RangeSet(2, doc="Side product outlets") m.bounds = RangeSet(2, doc="Number of boundary condition values") - m.candidate_trays_main = Set( - initialize=m.tray - [m.con_tray, m.reb_tray], - doc="Candidate trays for top and \ - bottom sections 1 and 4", - ) + m.candidate_trays_main = Set( + initialize=m.tray - [m.con_tray, m.reb_tray], doc="Candidate trays for top and \ + bottom sections 1 and 4" + ) m.candidate_trays_feed = Set( initialize=m.tray - [m.con_tray, m.feed_tray, m.reb_tray], doc="Candidate trays for feed section 2", @@ -243,7 +242,7 @@ def build_model(): for comp in m.comp: m.dHvap[comp] = dHvapb[comp] / m.Hscale - ## Heat capacity calculation for liquid and vapor phases for each component in the feed, section, and tray + ## Heat capacity calculation using the Ruczika-D method for each component in the feed, section, and tray for sec in m.section: for n_tray in m.tray: for comp in m.comp: @@ -267,7 +266,7 @@ def build_model(): / m.Hscale ) - ## Liquid and vapor enthalpy calculation for each component in the feed, section, and tray + ## Liquid and vapor enthalpy calculation using the Ruczika-D method for each component in the feed, section, and tray for sec in m.section: for n_tray in m.tray: for comp in m.comp: @@ -280,7 +279,7 @@ def build_model(): m.hvf = {} # Vapor enthalpy for side feed [J/mol] m.F0 = {} # Side feed flowrate per component [mol/s] - ## Heat capacity in liquid and vapor phases for side feed for each component. + ## Heat capacity in liquid and vapor phases for side feed using the Ruczika-D method. for comp in m.comp: for cp in m.cplv: m.cpdTf[comp, cp] = ( @@ -300,7 +299,7 @@ def build_model(): / m.Hscale ) - ## Side feed flowrate and liquid and vapor enthalpy calculation for each component in the feed + ## Side feed flowrate and enthalpy calculation using the Ruczika-D method for each component in the feed for comp in m.comp: m.F0[comp] = ( m.xfi[comp] * m.Fi @@ -338,8 +337,8 @@ def build_model(): bounds=lambda m, sec, tray, comp: ( m.prop[comp, 'TC'] / m.Tup, m.prop[comp, 'TC'] / m.Tlo, - ), # (0, None), - initialize=lambda m, sec, tray, comp: m.T0[sec, tray] / m.Tlo, + ), + initialize=lambda m, sec, tray, comp: m.prop[comp, 'TC'] / m.T0[sec, tray], ) m.x = Var( @@ -1460,8 +1459,8 @@ def _bottom_reduced_temperature(disj, comp): """ return m.Tr[1, n_tray, comp] * m.T[1, n_tray] == m.prop[comp, 'TC'] - @disj.Constraint(m.comp, doc="Bottom section 1 vapor composition") - def bottom_vapor_composition(disj, comp): + @disj.Constraint(m.comp, doc="Bottom section 1 vapor composition") + def bottom_vapor_composition(disj, comp): """ Vapor composition for the bottom section in the column. @@ -1478,32 +1477,24 @@ def bottom_vapor_composition(disj, comp): The constraint expression that enforces the vapor composition for the bottom section in the column. The equation is derived from the vapor-liquid equilibrium relationship. """ - return ( - m.y[1, n_tray, comp] * m.P[1, n_tray] - == m.x[1, n_tray, comp] - * ( - m.actv[1, n_tray, comp] - * ( - m.prop[comp, 'PC'] - * exp( - m.Tr[1, n_tray, comp] - # m.prop[comp, 'TC'] - # / m.T[1, n_tray] - * ( - m.prop[comp, 'vpA'] - * (1 - m.T[1, n_tray] / m.prop[comp, 'TC']) - + m.prop[comp, 'vpB'] - * (1 - m.T[1, n_tray] / m.prop[comp, 'TC']) ** 1.5 - + m.prop[comp, 'vpC'] - * (1 - m.T[1, n_tray] / m.prop[comp, 'TC']) ** 3 - + m.prop[comp, 'vpD'] - * (1 - m.T[1, n_tray] / m.prop[comp, 'TC']) ** 6 - ) - ) - ) - ) - # / m.P[1, n_tray] - ) + return m.y[1, n_tray, comp] * m.P[1, n_tray] == m.x[1, n_tray, comp] * ( + m.actv[1, n_tray, comp] + * ( + m.prop[comp, 'PC'] + * exp( + m.Tr[1, n_tray, comp] + * ( + m.prop[comp, 'vpA'] * (1 - m.T[1, n_tray] / m.prop[comp, 'TC']) + + m.prop[comp, 'vpB'] + * (1 - m.T[1, n_tray] / m.prop[comp, 'TC']) ** 1.5 + + m.prop[comp, 'vpC'] + * (1 - m.T[1, n_tray] / m.prop[comp, 'TC']) ** 3 + + m.prop[comp, 'vpD'] + * (1 - m.T[1, n_tray] / m.prop[comp, 'TC']) ** 6 + ) + ) + ) + ) @disj.Constraint(m.comp, doc="Bottom section 1 liquid enthalpy") def bottom_liquid_enthalpy(disj, comp): @@ -1768,10 +1759,10 @@ def feedside_vapor_composition_summation(disj): """ return sum(m.y[2, n_tray, comp] for comp in m.comp) - 1 == m.erry[2, n_tray] - @disj.Constraint( - m.comp, doc="Feed section 2 temperature term for the vapor composition equation" - ) - def _feedside_reduced_temperature(disj, comp): + @disj.Constraint( + m.comp, doc="Feed section 2 temperature term for the vapor composition equation" + ) + def _feedside_reduced_temperature(disj, comp): """Calculate the temperature term for the feedside section vapor composition. The constraint is used to avoid division by zero in the MINLP Hull reformulation. @@ -1779,7 +1770,7 @@ def _feedside_reduced_temperature(disj, comp): Parameters ---------- disj : Disjunct - The disjunct object for the bottom section in the column. + The disjunct object for the feed side section in the column. comp : int The component index. @@ -1790,8 +1781,8 @@ def _feedside_reduced_temperature(disj, comp): """ return m.Tr[2, n_tray, comp] * m.T[2, n_tray] == m.prop[comp, 'TC'] - @disj.Constraint(m.comp, doc="Feed section 2 vapor composition") - def feedside_vapor_composition(disj, comp): + @disj.Constraint(m.comp, doc="Feed section 2 vapor composition") + def feedside_vapor_composition(disj, comp): """ Vapor composition for the feed side section in the column. @@ -1808,32 +1799,24 @@ def feedside_vapor_composition(disj, comp): The constraint expression that enforces the vapor composition for the feed side section in the column. The equation is derived from the vapor-liquid equilibrium relationship. """ - return ( - m.y[2, n_tray, comp] * m.P[2, n_tray] - == m.x[2, n_tray, comp] - * ( - m.actv[2, n_tray, comp] - * ( - m.prop[comp, 'PC'] - * exp( - m.Tr[2, n_tray, comp] - # m.prop[comp, 'TC'] - # / m.T[2, n_tray] - * ( - m.prop[comp, 'vpA'] - * (1 - m.T[2, n_tray] / m.prop[comp, 'TC']) - + m.prop[comp, 'vpB'] - * (1 - m.T[2, n_tray] / m.prop[comp, 'TC']) ** 1.5 - + m.prop[comp, 'vpC'] - * (1 - m.T[2, n_tray] / m.prop[comp, 'TC']) ** 3 - + m.prop[comp, 'vpD'] - * (1 - m.T[2, n_tray] / m.prop[comp, 'TC']) ** 6 - ) - ) - ) - ) - # / m.P[2, n_tray] - ) + return m.y[2, n_tray, comp] * m.P[2, n_tray] == m.x[2, n_tray, comp] * ( + m.actv[2, n_tray, comp] + * ( + m.prop[comp, 'PC'] + * exp( + m.Tr[2, n_tray, comp] + * ( + m.prop[comp, 'vpA'] * (1 - m.T[2, n_tray] / m.prop[comp, 'TC']) + + m.prop[comp, 'vpB'] + * (1 - m.T[2, n_tray] / m.prop[comp, 'TC']) ** 1.5 + + m.prop[comp, 'vpC'] + * (1 - m.T[2, n_tray] / m.prop[comp, 'TC']) ** 3 + + m.prop[comp, 'vpD'] + * (1 - m.T[2, n_tray] / m.prop[comp, 'TC']) ** 6 + ) + ) + ) + ) @disj.Constraint(m.comp, doc="Feed section 2 liquid enthalpy") def feedside_liquid_enthalpy(disj, comp): @@ -2102,11 +2085,11 @@ def productside_vapor_composition_summation(disj): """ return sum(m.y[3, n_tray, comp] for comp in m.comp) - 1 == m.erry[3, n_tray] - @disj.Constraint( - m.comp, - doc="Product section 3 temperature term for the vapor composition equation", - ) - def _feedside_reduced_temperature(disj, comp): + @disj.Constraint( + m.comp, + doc="Product section 3 temperature term for the vapor composition equation", + ) + def _productside_reduced_temperature(disj, comp): """Calculate the temperature term for the product section vapor composition. The constraint is used to avoid division by zero in the MINLP Hull reformulation. @@ -2114,7 +2097,7 @@ def _feedside_reduced_temperature(disj, comp): Parameters ---------- disj : Disjunct - The disjunct object for the bottom section in the column. + The disjunct object for the product side section in the column. comp : int The component index. @@ -2125,8 +2108,8 @@ def _feedside_reduced_temperature(disj, comp): """ return m.Tr[3, n_tray, comp] * m.T[3, n_tray] == m.prop[comp, 'TC'] - @disj.Constraint(m.comp, doc="Product section 3 vapor composition") - def productside_vapor_composition(disj, comp): + @disj.Constraint(m.comp, doc="Product section 3 vapor composition") + def productside_vapor_composition(disj, comp): """ Vapor composition for the product side section in the column. @@ -2143,32 +2126,24 @@ def productside_vapor_composition(disj, comp): The constraint expression that enforces the vapor composition for the product side section in the column. The equation is derived from the vapor-liquid equilibrium relationship. """ - return ( - m.y[3, n_tray, comp] * m.P[3, n_tray] - == m.x[3, n_tray, comp] - * ( - m.actv[3, n_tray, comp] - * ( - m.prop[comp, 'PC'] - * exp( - m.Tr[3, n_tray, comp] - # m.prop[comp, 'TC'] - # / m.T[3, n_tray] - * ( - m.prop[comp, 'vpA'] - * (1 - m.T[3, n_tray] / m.prop[comp, 'TC']) - + m.prop[comp, 'vpB'] - * (1 - m.T[3, n_tray] / m.prop[comp, 'TC']) ** 1.5 - + m.prop[comp, 'vpC'] - * (1 - m.T[3, n_tray] / m.prop[comp, 'TC']) ** 3 - + m.prop[comp, 'vpD'] - * (1 - m.T[3, n_tray] / m.prop[comp, 'TC']) ** 6 - ) - ) - ) - ) - # / m.P[3, n_tray] - ) + return m.y[3, n_tray, comp] * m.P[3, n_tray] == m.x[3, n_tray, comp] * ( + m.actv[3, n_tray, comp] + * ( + m.prop[comp, 'PC'] + * exp( + m.Tr[3, n_tray, comp] + * ( + m.prop[comp, 'vpA'] * (1 - m.T[3, n_tray] / m.prop[comp, 'TC']) + + m.prop[comp, 'vpB'] + * (1 - m.T[3, n_tray] / m.prop[comp, 'TC']) ** 1.5 + + m.prop[comp, 'vpC'] + * (1 - m.T[3, n_tray] / m.prop[comp, 'TC']) ** 3 + + m.prop[comp, 'vpD'] + * (1 - m.T[3, n_tray] / m.prop[comp, 'TC']) ** 6 + ) + ) + ) + ) @disj.Constraint(m.comp, doc="Product section 3 liquid enthalpy") def productside_liquid_enthalpy(disj, comp): @@ -2310,7 +2285,7 @@ def _top_mass_percomponent_balances(disj, comp): m.D[comp] if n_tray == m.con_tray else 0 ) == 0 - @disj.Constraint(doc="Top scetion 4 energy balances") + @disj.Constraint(doc="Top section 4 energy balances") def _top_energy_balances(disj): """ Energy balances for the top section in the column. @@ -2445,10 +2420,10 @@ def top_vapor_composition_summation(disj): """ return sum(m.y[4, n_tray, comp] for comp in m.comp) - 1 == m.erry[4, n_tray] - @disj.Constraint( - m.comp, doc="Top section 4 temperature term for the vapor composition equation" - ) - def _feedside_reduced_temperature(disj, comp): + @disj.Constraint( + m.comp, doc="Top section 4 temperature term for the vapor composition equation" + ) + def _top_reduced_temperature(disj, comp): """Calculate the temperature term for the top section vapor composition. The constraint is used to avoid division by zero in the MINLP Hull reformulation. @@ -2456,7 +2431,7 @@ def _feedside_reduced_temperature(disj, comp): Parameters ---------- disj : Disjunct - The disjunct object for the bottom section in the column. + The disjunct object for the top section in the column. comp : int The component index. @@ -2467,7 +2442,7 @@ def _feedside_reduced_temperature(disj, comp): """ return m.Tr[4, n_tray, comp] * m.T[4, n_tray] == m.prop[comp, 'TC'] - @disj.Constraint(m.comp, doc="Top scetion 4 vapor composition") + @disj.Constraint(m.comp, doc="Top section 4 vapor composition") def top_vapor_composition(disj, comp): """ Vapor composition for the top section in the column. @@ -2485,32 +2460,24 @@ def top_vapor_composition(disj, comp): The constraint expression that enforces the vapor composition for the top section in the column. The equation is derived from the vapor-liquid equilibrium relationship. """ - return ( - m.y[4, n_tray, comp] * m.P[4, n_tray] - == m.x[4, n_tray, comp] - * ( - m.actv[4, n_tray, comp] - * ( - m.prop[comp, 'PC'] - * exp( - m.Tr[4, n_tray, comp] - # m.prop[comp, 'TC'] - # / m.T[4, n_tray] - * ( - m.prop[comp, 'vpA'] - * (1 - m.T[4, n_tray] / m.prop[comp, 'TC']) - + m.prop[comp, 'vpB'] - * (1 - m.T[4, n_tray] / m.prop[comp, 'TC']) ** 1.5 - + m.prop[comp, 'vpC'] - * (1 - m.T[4, n_tray] / m.prop[comp, 'TC']) ** 3 - + m.prop[comp, 'vpD'] - * (1 - m.T[4, n_tray] / m.prop[comp, 'TC']) ** 6 - ) - ) - ) - ) - # / m.P[4, n_tray] - ) + return m.y[4, n_tray, comp] * m.P[4, n_tray] == m.x[4, n_tray, comp] * ( + m.actv[4, n_tray, comp] + * ( + m.prop[comp, 'PC'] + * exp( + m.Tr[4, n_tray, comp] + * ( + m.prop[comp, 'vpA'] * (1 - m.T[4, n_tray] / m.prop[comp, 'TC']) + + m.prop[comp, 'vpB'] + * (1 - m.T[4, n_tray] / m.prop[comp, 'TC']) ** 1.5 + + m.prop[comp, 'vpC'] + * (1 - m.T[4, n_tray] / m.prop[comp, 'TC']) ** 3 + + m.prop[comp, 'vpD'] + * (1 - m.T[4, n_tray] / m.prop[comp, 'TC']) ** 6 + ) + ) + ) + ) @disj.Constraint(m.comp, doc="Top section 4 liquid enthalpy") def top_liquid_enthalpy(disj, comp): diff --git a/tests/test_pr47_regressions.py b/tests/test_pr47_regressions.py new file mode 100644 index 00000000..fea5d8a5 --- /dev/null +++ b/tests/test_pr47_regressions.py @@ -0,0 +1,31 @@ +import pytest +from pyomo.environ import TransformationFactory, value + +from gdplib.hda import HDA_model +from gdplib.kaibel import build_model as build_kaibel_model + + +def test_hda_model_constructs_with_scalar_heat_capacity_ratio(): + model = HDA_model() + + assert value(model.heat_capacity_ratio) == 1.3 + assert model.gamma.is_indexed() + + +def test_kaibel_reduced_temperature_initial_values_are_within_bounds(): + model = build_kaibel_model() + + out_of_bounds = [ + var.name + for var in model.Tr.values() + if value(var) < var.lb or value(var) > var.ub + ] + + assert out_of_bounds == [] + + +@pytest.mark.parametrize('transformation', ['gdp.hull', 'gdp.bigm']) +def test_kaibel_gdp_transformation_constructs(transformation): + model = build_kaibel_model() + + TransformationFactory(transformation).apply_to(model) From c19329c0795b1d3ac81ede1cd97bbe36daa0ee62 Mon Sep 17 00:00:00 2001 From: "David E. Bernal Neira" Date: Sun, 10 May 2026 15:44:42 -0400 Subject: [PATCH 4/6] Fix repository lint checks --- gdplib/biofuel/model.py | 2 +- gdplib/gdp_col/initialize.py | 6 +- gdplib/kaibel/kaibel_init.py | 48 +++++------ gdplib/kaibel/kaibel_prop.py | 2 +- gdplib/kaibel/kaibel_side_flash.py | 2 +- gdplib/kaibel/main_gdp.py | 6 +- gdplib/logical/positioning.py | 18 ++--- gdplib/logical/spectralog.py | 20 ++--- gdplib/methanol/methanol.py | 1 - gdplib/mod_hens/modular_discrete.py | 2 +- .../modular_discrete_single_module.py | 80 ++++++++++--------- gdplib/stranded_gas/model.py | 6 +- generate_model_size_report.py | 1 - setup.py | 2 +- 14 files changed, 91 insertions(+), 105 deletions(-) diff --git a/gdplib/biofuel/model.py b/gdplib/biofuel/model.py index 84820f42..75e48918 100644 --- a/gdplib/biofuel/model.py +++ b/gdplib/biofuel/model.py @@ -4,7 +4,7 @@ The model enforces constraints to ensure that raw material supplies do not exceed available amounts, product shipments meet market demands exactly, and production at each site matches outgoing shipments and available resources. It also optimizes transportation costs by managing both variable and fixed costs associated with active transportation routes. -The disjunctions in the model define the operational modes for facility sites (modular, conventional, or inactive) and the activity status of supply and product routes (active or inactive). +The disjunctions in the model define the operational modes for facility sites (modular, conventional, or inactive) and the activity status of supply and product routes (active or inactive). These elements allow the model to simulate different operational scenarios and strategic decisions, optimizing the network's layout and logistics based on economic and market conditions. The objective of the model is to optimize the network layout and production allocation to minimize total costs, which include setup and teardown of facilities, production costs, and transportation costs. diff --git a/gdplib/gdp_col/initialize.py b/gdplib/gdp_col/initialize.py index a5b0140c..2382f639 100644 --- a/gdplib/gdp_col/initialize.py +++ b/gdplib/gdp_col/initialize.py @@ -72,9 +72,9 @@ def set_value_if_not_fixed(var, val): for i in range(1, num_active_trays) ] for tray in range(2, num_active_trays): - indx = new_indices[tray - 1] - lower = floor(indx) - frac_above = indx - lower + index = new_indices[tray - 1] + lower = floor(index) + frac_above = index - lower # Take linear combination of values tray_indexed_data.loc[tray] = ( tray_indexed_data.loc[lower] * (1 - frac_above) diff --git a/gdplib/kaibel/kaibel_init.py b/gdplib/kaibel/kaibel_init.py index 5a510be4..03590951 100644 --- a/gdplib/kaibel/kaibel_init.py +++ b/gdplib/kaibel/kaibel_init.py @@ -1,28 +1,28 @@ """ - Calculation of the theoretical minimum number of trays and initial - temperature values. - (written by E. Soraya Rawlings, esoraya@rwlngs.net) - - The separation of four components require a sequence of at least three distillation - columns. Here, we calculate the minimum number of theoretical trays for the three - columns. The sequence is shown in Figure 2. - - COLUMN 1 COLUMN 2 COLUMN 3 - ----- ---- ----- - | | | | | | - ----- | A ----- | ----- | - | |<---> B -- | |<----> A -- | |<---> A - | | C | | | B | | | - A | | | | | | | | - B | | | | | | | | - C --->| | -->| | -->| | - D | | | | | | - | | | | | | - | |<- | |<- | |<- - ----- | ----- | ----- | - | | | | | | - -------> D -------> C -------> B - Figure 2. Sequence of columns for the separation of a quaternary mixture + Calculation of the theoretical minimum number of trays and initial + temperature values. + (written by E. Soraya Rawlings, esoraya@rwlngs.net) + +The separation of four components require a sequence of at least three distillation +columns. Here, we calculate the minimum number of theoretical trays for the three +columns. The sequence is shown in Figure 2. + + COLUMN 1 COLUMN 2 COLUMN 3 + ----- ---- ----- + | | | | | | + ----- | A ----- | ----- | + | |<---> B -- | |<----> A -- | |<---> A + | | C | | | B | | | + A | | | | | | | | + B | | | | | | | | + C --->| | -->| | -->| | + D | | | | | | + | | | | | | + | |<- | |<- | |<- + ----- | ----- | ----- | + | | | | | | + -------> D -------> C -------> B + Figure 2. Sequence of columns for the separation of a quaternary mixture """ from __future__ import division diff --git a/gdplib/kaibel/kaibel_prop.py b/gdplib/kaibel/kaibel_prop.py index 3cae0f13..afae66f0 100644 --- a/gdplib/kaibel/kaibel_prop.py +++ b/gdplib/kaibel/kaibel_prop.py @@ -1,4 +1,4 @@ -""" Properties of the system """ +"""Properties of the system""" from pyomo.environ import ConcreteModel diff --git a/gdplib/kaibel/kaibel_side_flash.py b/gdplib/kaibel/kaibel_side_flash.py index dfe38f4e..7d7af0ca 100644 --- a/gdplib/kaibel/kaibel_side_flash.py +++ b/gdplib/kaibel/kaibel_side_flash.py @@ -1,4 +1,4 @@ -""" Side feed flash """ +"""Side feed flash""" from pyomo.environ import ( ConcreteModel, diff --git a/gdplib/kaibel/main_gdp.py b/gdplib/kaibel/main_gdp.py index 3a66121c..0b817c90 100644 --- a/gdplib/kaibel/main_gdp.py +++ b/gdplib/kaibel/main_gdp.py @@ -135,8 +135,7 @@ def intro_message(m): """ - print( - """ + print(""" If you use this model and/or initialization strategy, you may cite the following: Rawlings, ES; Chen, Q; Grossmann, IE; Caballero, JA. Kaibel Column: Modeling, @@ -145,8 +144,7 @@ def intro_message(m): DOI: https://doi.org/10.1016/j.compchemeng.2019.03.006 - """ - ) + """) def display_results(m): diff --git a/gdplib/logical/positioning.py b/gdplib/logical/positioning.py index 8047c213..e345bd88 100644 --- a/gdplib/logical/positioning.py +++ b/gdplib/logical/positioning.py @@ -114,8 +114,7 @@ def build_model(): # Bounds for the locations location_bounds = {1: (2, 4.5), 2: (0, 8.0), 3: (3, 9.0), 4: (0, 5.0), 5: (4, 10)} - ideal_points_data = StringIO( - """ + ideal_points_data = StringIO(""" 1 2 3 4 5 1 2.26 5.15 4.03 1.74 4.74 2 5.51 9.01 3.84 1.47 9.92 @@ -142,8 +141,7 @@ def build_model(): 23 1.37 0.54 1.55 5.56 5.85 24 8.79 5.04 4.83 6.94 0.38 25 2.66 4.19 6.49 8.04 1.66 - """ - ) + """) ideal_points_table = pd.read_csv( ideal_points_data, delimiter=r'\s+' ) # ideal points for each consumer and product @@ -157,8 +155,7 @@ def build_model(): doc='Ideal points for each consumer and product', ) - weights_data = StringIO( - """ + weights_data = StringIO(""" 1 2 3 4 5 1 9.57 2.74 9.75 3.96 8.67 2 8.38 3.93 5.18 5.2 7.82 @@ -185,8 +182,7 @@ def build_model(): 23 1.47 5.71 6.95 1.42 3.49 24 5.4 3.12 5.37 6.1 3.71 25 6.32 0.81 6.12 6.73 7.93 - """ - ) + """) weights_table = pd.read_csv(weights_data, delimiter=r'\s+') weights_dict = { (k[0], int(k[1])): v for k, v in weights_table.stack().to_dict().items() @@ -198,8 +194,7 @@ def build_model(): doc='Weights for each consumer and product', ) - existing_products_data = StringIO( - """ + existing_products_data = StringIO(""" 1 2 3 4 5 1 0.62 5.06 7.82 0.22 4.42 2 5.21 2.66 9.54 5.03 8.01 @@ -211,8 +206,7 @@ def build_model(): 8 8.35 3.79 1.19 1.96 5.88 9 6.44 0.17 9.93 6.8 9.75 10 6.49 1.92 0.05 4.89 6.43 - """ - ) + """) existing_products_table = pd.read_csv(existing_products_data, delimiter=r'\s+') existing_products_dict = { (k[0], int(k[1])): v diff --git a/gdplib/logical/spectralog.py b/gdplib/logical/spectralog.py index ea0c7f35..0aab0860 100644 --- a/gdplib/logical/spectralog.py +++ b/gdplib/logical/spectralog.py @@ -6,7 +6,7 @@ [1] Vecchietti, A., & Grossmann, I. E. (1997). LOGMIP: a disjunctive 0-1 nonlinear optimizer for process systems models. Computers & chemical engineering, 21, S427-S432. https://doi.org/10.1016/S0098-1354(97)87539-4 [2] Brink, A., & Westerlund, T. (1995). The joint problem of model structure determination and parameter estimation in quantitative IR spectroscopy. Chemometrics and intelligent laboratory systems, 29(1), 29-36. https://doi.org/10.1016/0169-7439(95)00033-3 - + Optimal value: 12.0893 """ @@ -40,8 +40,7 @@ def build_model(): [2] Brink, A., & Westerlund, T. (1995). The joint problem of model structure determination and parameter estimation in quantitative IR spectroscopy. Chemometrics and intelligent laboratory systems, 29(1), 29-36. https://doi.org/10.1016/0169-7439(95)00033-3 """ # Matrix of absorbance values across different wave numbers (rows) and spectra numbers (columns) - spectroscopic_data = StringIO( - """ + spectroscopic_data = StringIO(""" 1 2 3 4 5 6 7 8 1 0.0003 0.0764 0.0318 0.0007 0.0534 0.0773 0.0536 0.0320 2 0.0007 0.0003 0.0004 0.0009 0.0005 0.0009 0.0005 0.0003 @@ -53,8 +52,7 @@ def build_model(): 8 0.0507 0.0361 0.0433 0.0635 0.0048 0.0891 0.0213 0.0310 9 0.0905 0.0600 0.0754 0.1098 0.0038 0.1443 0.0420 0.0574 10 0.0016 0.0209 0.0063 0.0010 0.0132 0.0203 0.0139 0.0057 - """ - ) + """) # Note: this could come from an external data file spectroscopic_data_table = pd.read_csv(spectroscopic_data, delimiter=r'\s+') flat_spectro_data = spectroscopic_data_table.stack() @@ -64,28 +62,24 @@ def build_model(): # Measured concentration data for each compound(row) and spectra number(column) # Units for concentration for each component 1, 2, and 3 are ppm, ppm, and % for CO, NO, and CO2, respectively. - c_data = StringIO( - """ + c_data = StringIO(""" 1 2 3 4 5 6 7 8 1 502 204 353 702 0 1016 104 204 2 97 351 351 351 700 0 201 97 3 0 22 8 0 14 22 14 8 - """ - ) + """) c_data_table = pd.read_csv(c_data, delimiter=r'\s+') c_data_dict = { (k[0], int(k[1])): v for k, v in c_data_table.stack().to_dict().items() } # Covariance matrix; It is assumed to be known that it is equal to the identity matrix at first problem iteration - r_data = StringIO( - """ + r_data = StringIO(""" 1 2 3 1 1 0 0 2 0 1 0 3 0 0 1 - """ - ) + """) r_data_table = pd.read_csv(r_data, delimiter=r'\s+') r_data_dict = { (k[0], int(k[1])): v for k, v in r_data_table.stack().to_dict().items() diff --git a/gdplib/methanol/methanol.py b/gdplib/methanol/methanol.py index fe859182..225b880b 100644 --- a/gdplib/methanol/methanol.py +++ b/gdplib/methanol/methanol.py @@ -4,7 +4,6 @@ import logging import pyomo.gdp as gdp - assert sys.version_info.major == 3 assert sys.version_info.minor >= 6 diff --git a/gdplib/mod_hens/modular_discrete.py b/gdplib/mod_hens/modular_discrete.py index 6c70f1f7..a7a4450a 100644 --- a/gdplib/mod_hens/modular_discrete.py +++ b/gdplib/mod_hens/modular_discrete.py @@ -81,7 +81,7 @@ def build_modular_option(cafaro_approx, num_stages): """ m = build_model(cafaro_approx, num_stages) - # Optimize for the least cost configuration across all stages and matche + # Optimize for the least cost configuration across all stages and matches for hot, cold in m.valid_matches: lowest_price = float('inf') for size in sorted(m.possible_sizes, reverse=True): diff --git a/gdplib/mod_hens/modular_discrete_single_module.py b/gdplib/mod_hens/modular_discrete_single_module.py index 0fbe9800..3eb7819b 100644 --- a/gdplib/mod_hens/modular_discrete_single_module.py +++ b/gdplib/mod_hens/modular_discrete_single_module.py @@ -64,7 +64,7 @@ def build_model(use_cafaro_approximation, num_stages): configs = [(i + 1, size) for i in range(m.max_num_modules[size])] configurations_list += configs - # Map of config indx: (# modules, module size) + # Map of config index: (# modules, module size) m.configurations_map = {(k + 1): v for k, v in enumerate(configurations_list)} m.module_index_set = RangeSet(len(configurations_list)) @@ -79,7 +79,7 @@ def build_model(use_cafaro_approximation, num_stages): ) @m.Param(m.module_index_set, doc="Area of each configuration") - def module_area(m, indx): + def module_area(m, index): """ Calculates the total area of a module configuration based on the number of modules and the size of each module. @@ -87,7 +87,7 @@ def module_area(m, indx): ---------- m : Pyomo.ConcreteModel A Pyomo concrete model representing the heat exchange model modified to support discretization of area simplifying the nonlinear expressions, specialized to the case of allowing only a single exchanger module type (size). - indx : int + index : int Index of the module configuration in the model. Returns @@ -95,7 +95,7 @@ def module_area(m, indx): Pyomo.Parameter The total area of the configuration corresponding to the given index. """ - num_modules, size = m.configurations_map[indx] + num_modules, size = m.configurations_map[index] return num_modules * size @m.Param( @@ -103,7 +103,7 @@ def module_area(m, indx): m.module_index_set, doc="Area cost for each modular configuration.", ) - def modular_size_cost(m, hot, cold, indx): + def modular_size_cost(m, hot, cold, index): """ Determines the cost associated with a specific modular configuration, taking into account the number of modules and their individual sizes. @@ -115,7 +115,7 @@ def modular_size_cost(m, hot, cold, indx): The index for the hot stream involved in the heat exchanger. cold : str The index for the cold stream involved in the heat exchanger. - indx : int + index : int Index of the module configuration in the model. Returns @@ -123,7 +123,7 @@ def modular_size_cost(m, hot, cold, indx): Pyomo.Parameter Cost associated with the specified modular configuration. """ - num_modules, size = m.configurations_map[indx] + num_modules, size = m.configurations_map[index] return num_modules * m.module_area_cost[hot, cold, size] @m.Param( @@ -131,7 +131,7 @@ def modular_size_cost(m, hot, cold, indx): m.module_index_set, doc="Fixed cost for each modular exchanger size.", ) - def modular_fixed_cost(m, hot, cold, indx): + def modular_fixed_cost(m, hot, cold, index): """ Computes the fixed cost for a modular exchanger configuration, factoring in the number of modules and the set fixed cost per unit. @@ -141,7 +141,7 @@ def modular_fixed_cost(m, hot, cold, indx): A Pyomo concrete model representing the heat exchange model modified for the case of allowing only a single exchanger module type (size). cold : str The index for the cold stream involved in the heat exchanger. - indx : int + index : int Index of the module configuration in the model. Returns @@ -149,7 +149,7 @@ def modular_fixed_cost(m, hot, cold, indx): Pyomo.Parameter Fixed cost for the given modular exchanger configuration. """ - num_modules, size = m.configurations_map[indx] + num_modules, size = m.configurations_map[index] return num_modules * m.module_fixed_unit_cost m.LMTD_discretize = Var( @@ -166,8 +166,8 @@ def modular_fixed_cost(m, hot, cold, indx): disj = m.exchanger_exists[hot, cold, stg] disj.choose_one_config = Constraint( expr=sum( - m.module_config_active[hot, cold, stg, indx] - for indx in m.module_index_set + m.module_config_active[hot, cold, stg, index] + for index in m.module_index_set ) == 1, doc="Enforce a single active configuration per exchanger per stage.", @@ -176,20 +176,20 @@ def modular_fixed_cost(m, hot, cold, indx): disj.exchanger_area_cost = Constraint( expr=m.exchanger_area_cost[stg, hot, cold] * 1e-3 == sum( - m.modular_size_cost[hot, cold, indx] + m.modular_size_cost[hot, cold, index] * 1e-3 - * m.module_config_active[hot, cold, stg, indx] - for indx in m.module_index_set + * m.module_config_active[hot, cold, stg, index] + for index in m.module_index_set ), doc="Compute total area cost from active configurations.", ) disj.exchanger_fixed_cost = Constraint( expr=m.exchanger_fixed_cost[stg, hot, cold] == sum( - m.modular_fixed_cost[hot, cold, indx] + m.modular_fixed_cost[hot, cold, index] * 1e-3 - * m.module_config_active[hot, cold, stg, indx] - for indx in m.module_index_set + * m.module_config_active[hot, cold, stg, index] + for index in m.module_index_set ), doc="Sum fixed costs of active configurations for total investment.", ) @@ -197,8 +197,8 @@ def modular_fixed_cost(m, hot, cold, indx): disj.discretize_area = Constraint( expr=m.exchanger_area[stg, hot, cold] == sum( - m.module_area[indx] * m.module_config_active[hot, cold, stg, indx] - for indx in m.module_index_set + m.module_area[index] * m.module_config_active[hot, cold, stg, index] + for index in m.module_index_set ), doc="Match exchanger area with sum of active configuration areas.", ) @@ -206,13 +206,13 @@ def modular_fixed_cost(m, hot, cold, indx): disj.discretized_LMTD = Constraint( expr=m.LMTD[hot, cold, stg] == sum( - m.LMTD_discretize[hot, cold, stg, indx] for indx in m.module_index_set + m.LMTD_discretize[hot, cold, stg, index] for index in m.module_index_set ), doc="Aggregate LMTD from active configurations for thermal modeling.", ) @disj.Constraint(m.module_index_set) - def discretized_LMTD_LB(disj, indx): + def discretized_LMTD_LB(disj, index): """ Sets the lower bound on the discretized Log Mean Temperature Difference (LMTD) for each module configuration. @@ -220,7 +220,7 @@ def discretized_LMTD_LB(disj, indx): ---------- disj : Pyomo.Disjunct The disjunct object representing a specific module configuration. - indx : int + index : int Index of the module configuration in the model. Returns @@ -229,11 +229,12 @@ def discretized_LMTD_LB(disj, indx): A constraint ensuring that the discretized LMTD respects the specified lower bound for active configurations. """ return ( - m.LMTD[hot, cold, stg].lb * m.module_config_active[hot, cold, stg, indx] - ) <= m.LMTD_discretize[hot, cold, stg, indx] + m.LMTD[hot, cold, stg].lb + * m.module_config_active[hot, cold, stg, index] + ) <= m.LMTD_discretize[hot, cold, stg, index] @disj.Constraint(m.module_index_set) - def discretized_LMTD_UB(disj, indx): + def discretized_LMTD_UB(disj, index): """ Sets the upper bound on the discretized Log Mean Temperature Difference (LMTD) for each module configuration. @@ -241,7 +242,7 @@ def discretized_LMTD_UB(disj, indx): ---------- disj : Pyomo.Disjunct The disjunct object representing a specific module configuration. - indx : int + index : int Index of the module configuration in the model. Returns @@ -249,15 +250,16 @@ def discretized_LMTD_UB(disj, indx): Pyomo.Constraint A constraint ensuring that the discretized LMTD does not exceed the specified upper bound for active configurations. """ - return m.LMTD_discretize[hot, cold, stg, indx] <= ( - m.LMTD[hot, cold, stg].ub * m.module_config_active[hot, cold, stg, indx] + return m.LMTD_discretize[hot, cold, stg, index] <= ( + m.LMTD[hot, cold, stg].ub + * m.module_config_active[hot, cold, stg, index] ) disj.exchanger_required_area = Constraint( expr=m.U[hot, cold] * sum( - m.module_area[indx] * m.LMTD_discretize[hot, cold, stg, indx] - for indx in m.module_index_set + m.module_area[index] * m.LMTD_discretize[hot, cold, stg, index] + for index in m.module_index_set ) >= m.heat_exchanged[hot, cold, stg], doc="Ensures sufficient heat transfer capacity for required heat exchange.", @@ -282,7 +284,7 @@ def module_type(disj, size): """ @disj.Constraint(m.valid_matches, m.stages, m.module_index_set) - def no_other_module_types(_, hot, cold, stg, indx): + def no_other_module_types(_, hot, cold, stg, index): """ Ensures only modules of the selected size are active, deactivating other sizes. @@ -296,7 +298,7 @@ def no_other_module_types(_, hot, cold, stg, indx): The index for the cold stream involved in the heat exchanger. stg : int The index for the stage involved in the heat exchanger. - indx : int + index : int Index of the module configuration in the model. Returns @@ -304,17 +306,17 @@ def no_other_module_types(_, hot, cold, stg, indx): Pyomo.Constraint A constraint expression that ensures only modules of the specified size are active, effectively disabling other module sizes for the current configuration. """ - # num_modules, size = configurations_map[indx] - if m.configurations_map[indx][1] != size: - return m.module_config_active[hot, cold, stg, indx] == 0 + # num_modules, size = configurations_map[index] + if m.configurations_map[index][1] != size: + return m.module_config_active[hot, cold, stg, index] == 0 else: return Constraint.NoConstraint # disj.no_other_module_types = Constraint( # expr=sum( - # m.module_config_active[hot, cold, stg, indx] - # for indx in m.module_index_set - # if m.configurations_map[indx][1] != size + # m.module_config_active[hot, cold, stg, index] + # for index in m.module_index_set + # if m.configurations_map[index][1] != size # ) # == 0, # doc="Deactivates non-selected module sizes.", diff --git a/gdplib/stranded_gas/model.py b/gdplib/stranded_gas/model.py index adc6720b..d7dabd53 100644 --- a/gdplib/stranded_gas/model.py +++ b/gdplib/stranded_gas/model.py @@ -1,10 +1,10 @@ """ model.py -Pyomo ConcreteModel for optimizing a modular stranded gas processing network. +Pyomo ConcreteModel for optimizing a modular stranded gas processing network. -The model is designed to convert stranded gas into gasoline using a modular and intensified GTL process. +The model is designed to convert stranded gas into gasoline using a modular and intensified GTL process. It incorporates the economic dynamics of module investments, gas processing, and product transportation. -Constraints manage the balance of gas supply and consumption, module availability and movement, and production capacities at potential sites. +Constraints manage the balance of gas supply and consumption, module availability and movement, and production capacities at potential sites. Disjunctions delineate operational scenarios, such as the existence or absence of pipelines and the activation status of sites, enabling dynamic and flexible system configuration. The objective function aims to maximize the network's net profit by optimizing revenue from gasoline sales while minimizing operational and capital expenditures across the network. diff --git a/generate_model_size_report.py b/generate_model_size_report.py index 96ddac08..f0978d2f 100644 --- a/generate_model_size_report.py +++ b/generate_model_size_report.py @@ -3,7 +3,6 @@ from pyomo.util.model_size import build_model_size_report import pandas as pd - if __name__ == "__main__": instance_list = [ # "batch_processing", diff --git a/setup.py b/setup.py index bd40d617..5a158b60 100644 --- a/setup.py +++ b/setup.py @@ -2,12 +2,12 @@ """ GDPlib open source model library for Generalized Disjunctive Programming """ + import sys from logging import warning from setuptools import setup, find_packages - kwargs = dict( name='gdplib', packages=find_packages(), From 8aceb546948f4c1bbe1e780e4db9db0ea8373612 Mon Sep 17 00:00:00 2001 From: "David E. Bernal Neira" Date: Sun, 10 May 2026 15:52:07 -0400 Subject: [PATCH 5/6] Fix benchmark stdout restoration --- benchmark.py | 18 +++++------ tests/test_benchmark.py | 67 +++++++++++++++++++++++++++++++++++++++++ 2 files changed, 75 insertions(+), 10 deletions(-) create mode 100644 tests/test_benchmark.py diff --git a/benchmark.py b/benchmark.py index 1bb9df84..ce347af0 100644 --- a/benchmark.py +++ b/benchmark.py @@ -1,7 +1,7 @@ import os import json import time -import sys +from contextlib import redirect_stdout from datetime import datetime from importlib import import_module from pyomo.environ import * @@ -28,14 +28,13 @@ def benchmark(model, strategy, timelimit, result_dir, subsolver="scip"): None """ model = model.clone() - stdout = sys.stdout + log_path = os.path.join(result_dir, f"{strategy}_{subsolver}.log") + json_path = os.path.join(result_dir, f"{strategy}_{subsolver}.json") if strategy in ["gdp.bigm", "gdp.hull"]: transformation_start_time = time.time() TransformationFactory(strategy).apply_to(model) transformation_end_time = time.time() - with open( - result_dir + "/" + strategy + "_" + subsolver + ".log", "w" - ) as sys.stdout: + with open(log_path, "w") as log_file, redirect_stdout(log_file): results = SolverFactory(subsolver).solve( model, tee=True, timelimit=timelimit ) @@ -50,9 +49,7 @@ def benchmark(model, strategy, timelimit, result_dir, subsolver="scip"): "gdpopt.lbb", "gdpopt.ric", ]: - with open( - result_dir + "/" + strategy + "_" + subsolver + ".log", "w" - ) as sys.stdout: + with open(log_path, "w") as log_file, redirect_stdout(log_file): results = SolverFactory(strategy).solve( model, tee=True, @@ -63,9 +60,10 @@ def benchmark(model, strategy, timelimit, result_dir, subsolver="scip"): time_limit=timelimit, ) print(results) + else: + raise ValueError(f"Unknown benchmark strategy: {strategy}") - sys.stdout = stdout - with open(result_dir + "/" + strategy + "_" + subsolver + ".json", "w") as f: + with open(json_path, "w") as f: json.dump(results.json_repn(), f) return None diff --git a/tests/test_benchmark.py b/tests/test_benchmark.py new file mode 100644 index 00000000..6de32018 --- /dev/null +++ b/tests/test_benchmark.py @@ -0,0 +1,67 @@ +import sys +from importlib.util import module_from_spec, spec_from_file_location +from pathlib import Path + +import pytest + +benchmark_path = Path(__file__).resolve().parents[1] / "benchmark.py" +benchmark_spec = spec_from_file_location("benchmark_module", benchmark_path) +assert benchmark_spec is not None +assert benchmark_spec.loader is not None +benchmark_module = module_from_spec(benchmark_spec) +benchmark_spec.loader.exec_module(benchmark_module) + + +class DummyModel: + def clone(self): + return self + + +class DummyTransformation: + def apply_to(self, model): + return None + + +class FailingSolver: + def solve(self, *args, **kwargs): + raise RuntimeError("solver failed") + + +def test_benchmark_restores_stdout_when_direct_solver_raises(tmp_path, monkeypatch): + original_stdout = sys.stdout + monkeypatch.setattr( + benchmark_module, + "TransformationFactory", + lambda strategy: DummyTransformation(), + ) + monkeypatch.setattr( + benchmark_module, "SolverFactory", lambda subsolver: FailingSolver() + ) + + with pytest.raises(RuntimeError, match="solver failed"): + benchmark_module.benchmark(DummyModel(), "gdp.bigm", 1, str(tmp_path)) + + assert sys.stdout is original_stdout + assert (tmp_path / "gdp.bigm_scip.log").exists() + assert not (tmp_path / "gdp.bigm_scip.json").exists() + + +def test_benchmark_restores_stdout_when_gdpopt_solver_raises(tmp_path, monkeypatch): + original_stdout = sys.stdout + monkeypatch.setattr( + benchmark_module, "SolverFactory", lambda strategy: FailingSolver() + ) + + with pytest.raises(RuntimeError, match="solver failed"): + benchmark_module.benchmark(DummyModel(), "gdpopt.loa", 1, str(tmp_path)) + + assert sys.stdout is original_stdout + assert (tmp_path / "gdpopt.loa_scip.log").exists() + assert not (tmp_path / "gdpopt.loa_scip.json").exists() + + +def test_benchmark_rejects_unknown_strategy(tmp_path): + with pytest.raises( + ValueError, match="Unknown benchmark strategy: no.such.strategy" + ): + benchmark_module.benchmark(DummyModel(), "no.such.strategy", 1, str(tmp_path)) From 90e289d2adb37ed931bef78009f1f404a1a57698 Mon Sep 17 00:00:00 2001 From: "David E. Bernal Neira" Date: Sun, 10 May 2026 16:23:57 -0400 Subject: [PATCH 6/6] Address PR 47 review follow-ups --- .github/workflows/typos.toml | 1 + gdplib/hda/HDA_GDP_gdpopt.py | 24 +- gdplib/kaibel/README.md | 5 +- gdplib/kaibel/kaibel_init.py | 798 ++--- gdplib/kaibel/kaibel_prop.py | 2 +- gdplib/kaibel/kaibel_solve_gdp.py | 5250 ++++++++++++++--------------- gdplib/kaibel/main_gdp.py | 574 ++-- tests/test_pr47_regressions.py | 63 +- 8 files changed, 3385 insertions(+), 3332 deletions(-) diff --git a/.github/workflows/typos.toml b/.github/workflows/typos.toml index a7e28f03..6fcec49c 100644 --- a/.github/workflows/typos.toml +++ b/.github/workflows/typos.toml @@ -3,3 +3,4 @@ hda = "hda" HDA = "HDA" equil = "equil" +gam = "gam" diff --git a/gdplib/hda/HDA_GDP_gdpopt.py b/gdplib/hda/HDA_GDP_gdpopt.py index 05bea29d..a4cfef1a 100644 --- a/gdplib/hda/HDA_GDP_gdpopt.py +++ b/gdplib/hda/HDA_GDP_gdpopt.py @@ -20,7 +20,7 @@ def HDA_model(): m.alpha = Param(initialize=0.3665, doc="compressor coefficient") m.compeff = Param(initialize=0.750, doc="compressor efficiency") - m.heat_capacity_ratio = Param(initialize=1.300, doc="ratio of cp to cv") + m.gam = Param(initialize=1.300, doc="ratio of cp to cv") m.abseff = Param(initialize=0.333, doc="absorber tray efficiency") m.disteff = Param(initialize=0.5000, doc="column tray efficiency") m.uflow = Param(initialize=50, doc="upper bound - flow logicals") @@ -1312,7 +1312,7 @@ def Compelec(_m, comp_): * m.f[stream] / 60.0 * (1.0 / m.compeff) - * (m.heat_capacity_ratio / (m.heat_capacity_ratio - 1.0)) + * (m.gam / (m.gam - 1.0)) for (comp1, stream) in m.icomp if comp_ == comp1 ) @@ -1324,13 +1324,9 @@ def Compelec(_m, comp_): def Ratio(_m, comp_): if comp == comp_: - return m.presrat[comp_] ** ( - m.heat_capacity_ratio / (m.heat_capacity_ratio - 1.0) - ) == sum( + return m.presrat[comp_] ** (m.gam / (m.gam - 1.0)) == sum( m.p[stream] for (comp1, stream) in m.ocomp if comp_ == comp1 - ) / sum( - m.p[stream] for (comp1, stream) in m.icomp if comp1 == comp_ - ) + ) / sum(m.p[stream] for (comp1, stream) in m.icomp if comp1 == comp_) return Constraint.Skip b.ratio = Constraint([comp], rule=Ratio, doc='pressure ratio (out to in)') @@ -2252,19 +2248,11 @@ def Valcmb(_m, valve, compon): def Valt(_m, valve): return sum( - m.t[stream] - / ( - m.p[stream] - ** ((m.heat_capacity_ratio - 1.0) / m.heat_capacity_ratio) - ) + m.t[stream] / (m.p[stream] ** ((m.gam - 1.0) / m.gam)) for (valv, stream) in m.oval if valv == valve ) == sum( - m.t[stream] - / ( - m.p[stream] - ** ((m.heat_capacity_ratio - 1.0) / m.heat_capacity_ratio) - ) + m.t[stream] / (m.p[stream] ** ((m.gam - 1.0) / m.gam)) for (valv, stream) in m.ival if valv == valve ) diff --git a/gdplib/kaibel/README.md b/gdplib/kaibel/README.md index c3f480d9..d3e5cb5c 100644 --- a/gdplib/kaibel/README.md +++ b/gdplib/kaibel/README.md @@ -6,6 +6,10 @@ Source paper: > Rawlings, E. S., Chen, Q., Grossmann, I. E., & Caballero, J. A. (2019). Kaibel Column: Modeling, optimization, and conceptual design of multi-product dividing wall columns. *Computers and Chemical Engineering*, 125, 31–39. https://doi.org/10.1016/j.compchemeng.2019.03.006 +Heat-capacity and enthalpy calculations use the Ruzicka-Domalski group-additivity heat-capacity method: + +> Ruzicka, V. Jr., & Domalski, E. S. Estimation of the heat capacities of organic liquids as a function of temperature using group additivity. *Journal of Physical and Chemical Reference Data*, 22(3). https://www.osti.gov/biblio/6251468 + ## Problem Details ### Solution @@ -30,4 +34,3 @@ Best known objective value: 99,709 - ``disjtn``: disjunctions - diff --git a/gdplib/kaibel/kaibel_init.py b/gdplib/kaibel/kaibel_init.py index 03590951..c1131b5f 100644 --- a/gdplib/kaibel/kaibel_init.py +++ b/gdplib/kaibel/kaibel_init.py @@ -1,399 +1,399 @@ -""" - Calculation of the theoretical minimum number of trays and initial - temperature values. - (written by E. Soraya Rawlings, esoraya@rwlngs.net) - -The separation of four components require a sequence of at least three distillation -columns. Here, we calculate the minimum number of theoretical trays for the three -columns. The sequence is shown in Figure 2. - - COLUMN 1 COLUMN 2 COLUMN 3 - ----- ---- ----- - | | | | | | - ----- | A ----- | ----- | - | |<---> B -- | |<----> A -- | |<---> A - | | C | | | B | | | - A | | | | | | | | - B | | | | | | | | - C --->| | -->| | -->| | - D | | | | | | - | | | | | | - | |<- | |<- | |<- - ----- | ----- | ----- | - | | | | | | - -------> D -------> C -------> B - Figure 2. Sequence of columns for the separation of a quaternary mixture -""" - -from __future__ import division - -from pyomo.environ import ( - exp, - log10, - minimize, - NonNegativeReals, - Objective, - RangeSet, - SolverFactory, - value, - Var, -) - -from gdplib.kaibel.kaibel_prop import get_model_with_properties - -# from .kaibel_prop import get_model_with_properties - - -def initialize_kaibel(): - """Initialize the Kaibel optimization model. - - This function initializes the Kaibel optimization model by setting up the operating conditions, - initial liquid compositions, and creating the necessary variables and constraints. - - Returns - ------- - None - """ - - ## Get the model with properties from kaibel_prop.py - m = get_model_with_properties() - - ## Operating conditions - m.Preb = 1.2 # Reboiler pressure in bar - m.Pcon = 1.05 # Condenser pressure in bar - m.Pf = 1.02 - - Pnmin = {} # Pressure in bars - Pnmin[1] = m.Preb # Reboiler pressure in bars - Pnmin[3] = m.Pcon # Distillate pressure in bars - Pnmin[2] = m.Pf # Side feed pressure in bars - - xi_nmin = {} # Initial liquid composition: first number = column and - # second number = 1 reboiler, 2 side feed, and - # 3 for condenser - - ## Column 1 - c_c1 = 4 # Components in Column 1 - lc_c1 = 3 # Light component in Column 1 - hc_c1 = 4 # Heavy component in Column 1 - inter1_c1 = 1 # Intermediate component in Column 1 - inter2_c1 = 2 # Intermediate component in Column 1 - - xi_nmin[1, 1, hc_c1] = 0.999 - xi_nmin[1, 1, lc_c1] = (1 - xi_nmin[1, 1, hc_c1]) / (c_c1 - 1) - xi_nmin[1, 1, inter1_c1] = (1 - xi_nmin[1, 1, hc_c1]) / (c_c1 - 1) - xi_nmin[1, 1, inter2_c1] = (1 - xi_nmin[1, 1, hc_c1]) / (c_c1 - 1) - xi_nmin[1, 3, lc_c1] = 0.33 - xi_nmin[1, 3, inter1_c1] = 0.33 - xi_nmin[1, 3, inter2_c1] = 0.33 - xi_nmin[1, 3, hc_c1] = 1 - ( - xi_nmin[1, 3, lc_c1] + xi_nmin[1, 3, inter1_c1] + xi_nmin[1, 3, inter2_c1] - ) - xi_nmin[1, 2, lc_c1] = 1 / c_c1 - xi_nmin[1, 2, inter1_c1] = 1 / c_c1 - xi_nmin[1, 2, inter2_c1] = 1 / c_c1 - xi_nmin[1, 2, hc_c1] = 1 / c_c1 - - ## Column 2 - c_c2 = 3 # Light components in Column 2 - lc_c2 = 2 # Light component in Column 2 - hc_c2 = 3 # Heavy component in Column 2 - inter_c2 = 1 # Intermediate component in Column 2 - - xi_nmin[2, 1, hc_c2] = 0.999 - xi_nmin[2, 1, lc_c2] = (1 - xi_nmin[2, 1, hc_c2]) / (c_c2 - 1) - xi_nmin[2, 1, inter_c2] = (1 - xi_nmin[2, 1, hc_c2]) / (c_c2 - 1) - xi_nmin[2, 3, lc_c2] = 0.499 - xi_nmin[2, 3, inter_c2] = 0.499 - xi_nmin[2, 3, hc_c2] = 1 - (xi_nmin[2, 3, lc_c2] + xi_nmin[2, 3, inter_c2]) - xi_nmin[2, 2, lc_c2] = 1 / c_c2 - xi_nmin[2, 2, inter_c2] = 1 / c_c2 - xi_nmin[2, 2, hc_c2] = 1 / c_c2 - - ## Column 3 - c_c3 = 2 # Components in Column 3 - lc_c3 = 1 # Light component in Column 3 - hc_c3 = 2 # Heavy component in Column 3 - - xi_nmin[3, 1, hc_c3] = 0.999 - xi_nmin[3, 1, lc_c3] = 1 - xi_nmin[3, 1, hc_c3] - xi_nmin[3, 3, lc_c3] = 0.999 - xi_nmin[3, 3, hc_c3] = 1 - xi_nmin[3, 3, lc_c3] - xi_nmin[3, 2, lc_c3] = 0.50 - xi_nmin[3, 2, hc_c3] = 0.50 - - #### - - mn = m.clone() # Clone the model to add the initialization code - - mn.name = "Initialization Code" - - mn.cols = RangeSet(3, doc='Number of columns ') - mn.sec = RangeSet(3, doc='Sections in column: 1 reb, 2 side feed, 3 cond') - mn.nc1 = RangeSet(c_c1, doc='Number of components in Column 1') - mn.nc2 = RangeSet(c_c2, doc='Number of components in Column 2') - mn.nc3 = RangeSet(c_c3, doc='Number of components in Column 3') - - mn.Tnmin = Var( - mn.cols, - mn.sec, - doc='Temperature in K', - bounds=(0, 500), - domain=NonNegativeReals, - ) - mn.Tr1nmin = Var( - mn.cols, - mn.sec, - mn.nc1, - doc='Temperature term for vapor pressure', - domain=NonNegativeReals, - bounds=(0, None), - ) - mn.Tr2nmin = Var( - mn.cols, - mn.sec, - mn.nc2, - doc='Temperature term for vapor pressure', - domain=NonNegativeReals, - bounds=(0, None), - ) - mn.Tr3nmin = Var( - mn.cols, - mn.sec, - mn.nc3, - doc='Temperature term for vapor pressure', - domain=NonNegativeReals, - bounds=(0, None), - ) - - @mn.Constraint(mn.cols, mn.sec, mn.nc1, doc="Temperature term for vapor pressure") - def _column1_reduced_temperature(mn, col, sec, nc): - """Calculate the reduced temperature for column 1. - - This function calculates the reduced temperature for column 1 based on the given parameters using the Peng-Robinson equation of state. - - Parameters - ---------- - mn : Pyomo ConcreteModel - The optimization model - col : int - The column index - sec : int - The section index - nc : int - The component index in column 1 - - Returns - ------- - Constraint - The constraint statement to calculate the reduced temperature. - """ - return mn.Tr1nmin[col, sec, nc] * mn.Tnmin[col, sec] == mn.prop[nc, 'TC'] - - @mn.Constraint(mn.cols, mn.sec, mn.nc2, doc="Temperature term for vapor pressure") - def _column2_reduced_temperature(mn, col, sec, nc): - """Calculate the reduced temperature for column 2. - - This function calculates the reduced temperature for column 2 based on the given parameters using the Peng-Robinson equation of state. - - Parameters - ---------- - mn : Pyomo ConcreteModel - The optimization model - col : int - The column index - sec : int - The section index - nc : int - The component index in column 2 - - Returns - ------- - Constraint - The constraint equation to calculate the reduced temperature - """ - return mn.Tr2nmin[col, sec, nc] * mn.Tnmin[col, sec] == mn.prop[nc, 'TC'] - - @mn.Constraint(mn.cols, mn.sec, mn.nc3, doc="Temperature term for vapor pressure") - def _column3_reduced_temperature(mn, col, sec, nc): - """Calculate the reduced temperature for column 3. - - This function calculates the reduced temperature for column 3 based on the given parameters. - - Parameters - ---------- - mn : Pyomo ConcreteModel - The optimization model - col : int - The column index - sec : int - The section index - nc : int - The component index in column 3 - - Returns - ------- - Constraint - The constraint equation to calculate the reduced temperature in column 3 - """ - return mn.Tr3nmin[col, sec, nc] * mn.Tnmin[col, sec] == mn.prop[nc, 'TC'] - - @mn.Constraint(mn.cols, mn.sec, doc="Boiling point temperature") - def _equilibrium_equation(mn, col, sec): - """Equilibrium equations for a given column and section. - - Parameters - ---------- - mn : Pyomo ConcreteModel - The optimization model object with properties - col : int - The column index - sec : int - The section index - - Returns - ------- - Constraint - The constraint equation to calculate the boiling point temperature using the Peng-Robinson equation of state - """ - if col == 1: - a = mn.Tr1nmin - b = mn.nc1 - elif col == 2: - a = mn.Tr2nmin - b = mn.nc2 - elif col == 3: - a = mn.Tr3nmin - b = mn.nc3 - return ( - sum( - xi_nmin[col, sec, nc] - * mn.prop[nc, 'PC'] - * exp( - a[col, sec, nc] - * ( - mn.prop[nc, 'vpA'] - * (1 - mn.Tnmin[col, sec] / mn.prop[nc, 'TC']) - + mn.prop[nc, 'vpB'] - * (1 - mn.Tnmin[col, sec] / mn.prop[nc, 'TC']) ** 1.5 - + mn.prop[nc, 'vpC'] - * (1 - mn.Tnmin[col, sec] / mn.prop[nc, 'TC']) ** 3 - + mn.prop[nc, 'vpD'] - * (1 - mn.Tnmin[col, sec] / mn.prop[nc, 'TC']) ** 6 - ) - ) - / Pnmin[sec] - for nc in b - ) - == 1 - ) - - mn.OBJ = Objective(expr=1, sense=minimize) - - #### - - SolverFactory('ipopt').solve(mn) - - yc = {} # Vapor composition - kl = {} # Light key component - kh = {} # Heavy key component - alpha = {} # Relative volatility of kl - ter = {} # Term to calculate the minimum number of trays - Nmin = {} # Minimum number of stages - Nminopt = {} # Total optimal minimum number of trays - Nfeed = {} # Side feed optimal location using Kirkbride's method: - # 1 = number of trays in rectifying section and - # 2 = number of trays in stripping section - side_feed = {} # Side feed location - av_alpha = {} # Average relative volatilities - xi_lhc = {} # Liquid composition in columns - rel = mn.Bdes / mn.Ddes # Ratio between products flowrates - ln = {} # Light component for the different columns - hn = {} # Heavy component for the different columns - ln[1] = lc_c1 - ln[2] = lc_c2 - ln[3] = lc_c3 - hn[1] = hc_c1 - hn[2] = hc_c2 - hn[3] = hc_c3 - - for col in mn.cols: - if col == 1: - b = mn.nc1 - elif col == 2: - b = mn.nc2 - else: - b = mn.nc3 - # For each component in the column and section calculate the vapor composition with the Peng-Robinson equation of state - for sec in mn.sec: - for nc in b: - yc[col, sec, nc] = ( - xi_nmin[col, sec, nc] - * mn.prop[nc, 'PC'] - * exp( - mn.prop[nc, 'TC'] - / value(mn.Tnmin[col, sec]) - * ( - mn.prop[nc, 'vpA'] - * (1 - value(mn.Tnmin[col, sec]) / mn.prop[nc, 'TC']) - + mn.prop[nc, 'vpB'] - * (1 - value(mn.Tnmin[col, sec]) / mn.prop[nc, 'TC']) ** 1.5 - + mn.prop[nc, 'vpC'] - * (1 - value(mn.Tnmin[col, sec]) / mn.prop[nc, 'TC']) ** 3 - + mn.prop[nc, 'vpD'] - * (1 - value(mn.Tnmin[col, sec]) / mn.prop[nc, 'TC']) ** 6 - ) - ) - ) / Pnmin[ - sec - ] # Vapor composition in the different sections for the different components in the columns - - for col in mn.cols: - # Calculate the relative volatility of the light and heavy components in the different sections for the different columns - xi_lhc[col, 4] = ( - xi_nmin[col, 1, ln[col]] / xi_nmin[col, 3, hn[col]] - ) # Liquid composition in the different sections with the initial liquid composition of the components in the different sections and columns and ln and hn which are the light and heavy components in the different columns - for sec in mn.sec: - kl[col, sec] = ( - yc[col, sec, ln[col]] / xi_nmin[col, sec, ln[col]] - ) # Light component in the different sections - kh[col, sec] = ( - yc[col, sec, hn[col]] / xi_nmin[col, sec, hn[col]] - ) # Heavy component in the different sections - xi_lhc[col, sec] = ( - xi_nmin[col, sec, hn[col]] / xi_nmin[col, sec, ln[col]] - ) # Liquid composition in the different sections - alpha[col, sec] = ( - kl[col, sec] / kh[col, sec] - ) # Relative volatility in the different sections - - for col in mn.cols: - # Calculate the average relative volatilities and the minimum number of trays with Fenske's and Kirkbride's method - av_alpha[col] = (alpha[col, 1] * alpha[col, 2] * alpha[col, 3]) ** ( - 1 / 3 - ) # Average relative volatilities calculated with the relative volatilities of the components in the three sections - Nmin[col] = log10((1 / xi_lhc[col, 3]) * xi_lhc[col, 1]) / log10( - av_alpha[col] - ) # Minimum number of trays calculated with Fenske's method - ter[col] = ( - rel * xi_lhc[col, 2] * (xi_lhc[col, 4] ** 2) - ) ** 0.206 # Term to calculate the minimum number of trays with Kirkbride's method - # Side feed optimal location using Kirkbride's method - Nfeed[1, col] = ( - ter[col] * Nmin[col] / (1 + ter[col]) - ) # Number of trays in rectifying section - Nfeed[2, col] = Nfeed[1, col] / ter[col] # Number of trays in stripping section - side_feed[col] = Nfeed[2, col] # Side feed location - - m.Nmintot = sum(Nmin[col] for col in mn.cols) # Total minimum number of trays - m.Knmin = int(m.Nmintot) + 1 # Total optimal minimum number of trays - - m.TB0 = value(mn.Tnmin[1, 1]) # Reboiler temperature in K in column 1 - m.Tf0 = value(mn.Tnmin[1, 2]) # Side feed temperature in K in column 1 - m.TD0 = value(mn.Tnmin[2, 3]) # Distillate temperature in K in column 2 - - return m - - -if __name__ == "__main__": - initialize_kaibel() +""" + Calculation of the theoretical minimum number of trays and initial + temperature values. + (written by E. Soraya Rawlings, esoraya@rwlngs.net) + +The separation of four components require a sequence of at least three distillation +columns. Here, we calculate the minimum number of theoretical trays for the three +columns. The sequence is shown in Figure 2. + + COLUMN 1 COLUMN 2 COLUMN 3 + ----- ---- ----- + | | | | | | + ----- | A ----- | ----- | + | |<---> B -- | |<----> A -- | |<---> A + | | C | | | B | | | + A | | | | | | | | + B | | | | | | | | + C --->| | -->| | -->| | + D | | | | | | + | | | | | | + | |<- | |<- | |<- + ----- | ----- | ----- | + | | | | | | + -------> D -------> C -------> B + Figure 2. Sequence of columns for the separation of a quaternary mixture +""" + +from __future__ import division + +from pyomo.environ import ( + exp, + log10, + minimize, + NonNegativeReals, + Objective, + RangeSet, + SolverFactory, + value, + Var, +) + +from gdplib.kaibel.kaibel_prop import get_model_with_properties + +# from .kaibel_prop import get_model_with_properties + + +def initialize_kaibel(): + """Initialize the Kaibel optimization model. + + This function initializes the Kaibel optimization model by setting up the operating conditions, + initial liquid compositions, and creating the necessary variables and constraints. + + Returns + ------- + None + """ + + ## Get the model with properties from kaibel_prop.py + m = get_model_with_properties() + + ## Operating conditions + m.Preb = 1.2 # Reboiler pressure in bar + m.Pcon = 1.05 # Condenser pressure in bar + m.Pf = 1.02 + + Pnmin = {} # Pressure in bars + Pnmin[1] = m.Preb # Reboiler pressure in bars + Pnmin[3] = m.Pcon # Distillate pressure in bars + Pnmin[2] = m.Pf # Side feed pressure in bars + + xi_nmin = {} # Initial liquid composition: first number = column and + # second number = 1 reboiler, 2 side feed, and + # 3 for condenser + + ## Column 1 + c_c1 = 4 # Components in Column 1 + lc_c1 = 3 # Light component in Column 1 + hc_c1 = 4 # Heavy component in Column 1 + inter1_c1 = 1 # Intermediate component in Column 1 + inter2_c1 = 2 # Intermediate component in Column 1 + + xi_nmin[1, 1, hc_c1] = 0.999 + xi_nmin[1, 1, lc_c1] = (1 - xi_nmin[1, 1, hc_c1]) / (c_c1 - 1) + xi_nmin[1, 1, inter1_c1] = (1 - xi_nmin[1, 1, hc_c1]) / (c_c1 - 1) + xi_nmin[1, 1, inter2_c1] = (1 - xi_nmin[1, 1, hc_c1]) / (c_c1 - 1) + xi_nmin[1, 3, lc_c1] = 0.33 + xi_nmin[1, 3, inter1_c1] = 0.33 + xi_nmin[1, 3, inter2_c1] = 0.33 + xi_nmin[1, 3, hc_c1] = 1 - ( + xi_nmin[1, 3, lc_c1] + xi_nmin[1, 3, inter1_c1] + xi_nmin[1, 3, inter2_c1] + ) + xi_nmin[1, 2, lc_c1] = 1 / c_c1 + xi_nmin[1, 2, inter1_c1] = 1 / c_c1 + xi_nmin[1, 2, inter2_c1] = 1 / c_c1 + xi_nmin[1, 2, hc_c1] = 1 / c_c1 + + ## Column 2 + c_c2 = 3 # Light components in Column 2 + lc_c2 = 2 # Light component in Column 2 + hc_c2 = 3 # Heavy component in Column 2 + inter_c2 = 1 # Intermediate component in Column 2 + + xi_nmin[2, 1, hc_c2] = 0.999 + xi_nmin[2, 1, lc_c2] = (1 - xi_nmin[2, 1, hc_c2]) / (c_c2 - 1) + xi_nmin[2, 1, inter_c2] = (1 - xi_nmin[2, 1, hc_c2]) / (c_c2 - 1) + xi_nmin[2, 3, lc_c2] = 0.499 + xi_nmin[2, 3, inter_c2] = 0.499 + xi_nmin[2, 3, hc_c2] = 1 - (xi_nmin[2, 3, lc_c2] + xi_nmin[2, 3, inter_c2]) + xi_nmin[2, 2, lc_c2] = 1 / c_c2 + xi_nmin[2, 2, inter_c2] = 1 / c_c2 + xi_nmin[2, 2, hc_c2] = 1 / c_c2 + + ## Column 3 + c_c3 = 2 # Components in Column 3 + lc_c3 = 1 # Light component in Column 3 + hc_c3 = 2 # Heavy component in Column 3 + + xi_nmin[3, 1, hc_c3] = 0.999 + xi_nmin[3, 1, lc_c3] = 1 - xi_nmin[3, 1, hc_c3] + xi_nmin[3, 3, lc_c3] = 0.999 + xi_nmin[3, 3, hc_c3] = 1 - xi_nmin[3, 3, lc_c3] + xi_nmin[3, 2, lc_c3] = 0.50 + xi_nmin[3, 2, hc_c3] = 0.50 + + #### + + mn = m.clone() # Clone the model to add the initialization code + + mn.name = "Initialization Code" + + mn.cols = RangeSet(3, doc='Number of columns ') + mn.sec = RangeSet(3, doc='Sections in column: 1 reb, 2 side feed, 3 cond') + mn.nc1 = RangeSet(c_c1, doc='Number of components in Column 1') + mn.nc2 = RangeSet(c_c2, doc='Number of components in Column 2') + mn.nc3 = RangeSet(c_c3, doc='Number of components in Column 3') + + mn.Tnmin = Var( + mn.cols, + mn.sec, + doc='Temperature in K', + bounds=(0, 500), + domain=NonNegativeReals, + ) + mn.Tr1nmin = Var( + mn.cols, + mn.sec, + mn.nc1, + doc='Temperature term for vapor pressure', + domain=NonNegativeReals, + bounds=(0, None), + ) + mn.Tr2nmin = Var( + mn.cols, + mn.sec, + mn.nc2, + doc='Temperature term for vapor pressure', + domain=NonNegativeReals, + bounds=(0, None), + ) + mn.Tr3nmin = Var( + mn.cols, + mn.sec, + mn.nc3, + doc='Temperature term for vapor pressure', + domain=NonNegativeReals, + bounds=(0, None), + ) + + @mn.Constraint(mn.cols, mn.sec, mn.nc1, doc="Temperature term for vapor pressure") + def _column1_reduced_temperature(mn, col, sec, nc): + """Calculate the reduced temperature for column 1. + + This function calculates the reduced temperature for column 1 based on the given parameters using the Peng-Robinson equation of state. + + Parameters + ---------- + mn : Pyomo ConcreteModel + The optimization model + col : int + The column index + sec : int + The section index + nc : int + The component index in column 1 + + Returns + ------- + Constraint + The constraint statement to calculate the reduced temperature. + """ + return mn.Tr1nmin[col, sec, nc] * mn.Tnmin[col, sec] == mn.prop[nc, 'TC'] + + @mn.Constraint(mn.cols, mn.sec, mn.nc2, doc="Temperature term for vapor pressure") + def _column2_reduced_temperature(mn, col, sec, nc): + """Calculate the reduced temperature for column 2. + + This function calculates the reduced temperature for column 2 based on the given parameters using the Peng-Robinson equation of state. + + Parameters + ---------- + mn : Pyomo ConcreteModel + The optimization model + col : int + The column index + sec : int + The section index + nc : int + The component index in column 2 + + Returns + ------- + Constraint + The constraint equation to calculate the reduced temperature + """ + return mn.Tr2nmin[col, sec, nc] * mn.Tnmin[col, sec] == mn.prop[nc, 'TC'] + + @mn.Constraint(mn.cols, mn.sec, mn.nc3, doc="Temperature term for vapor pressure") + def _column3_reduced_temperature(mn, col, sec, nc): + """Calculate the reduced temperature for column 3. + + This function calculates the reduced temperature for column 3 based on the given parameters. + + Parameters + ---------- + mn : Pyomo ConcreteModel + The optimization model + col : int + The column index + sec : int + The section index + nc : int + The component index in column 3 + + Returns + ------- + Constraint + The constraint equation to calculate the reduced temperature in column 3 + """ + return mn.Tr3nmin[col, sec, nc] * mn.Tnmin[col, sec] == mn.prop[nc, 'TC'] + + @mn.Constraint(mn.cols, mn.sec, doc="Boiling point temperature") + def _equilibrium_equation(mn, col, sec): + """Equilibrium equations for a given column and section. + + Parameters + ---------- + mn : Pyomo ConcreteModel + The optimization model object with properties + col : int + The column index + sec : int + The section index + + Returns + ------- + Constraint + The constraint equation to calculate the boiling point temperature using the Peng-Robinson equation of state + """ + if col == 1: + a = mn.Tr1nmin + b = mn.nc1 + elif col == 2: + a = mn.Tr2nmin + b = mn.nc2 + elif col == 3: + a = mn.Tr3nmin + b = mn.nc3 + return ( + sum( + xi_nmin[col, sec, nc] + * mn.prop[nc, 'PC'] + * exp( + a[col, sec, nc] + * ( + mn.prop[nc, 'vpA'] + * (1 - mn.Tnmin[col, sec] / mn.prop[nc, 'TC']) + + mn.prop[nc, 'vpB'] + * (1 - mn.Tnmin[col, sec] / mn.prop[nc, 'TC']) ** 1.5 + + mn.prop[nc, 'vpC'] + * (1 - mn.Tnmin[col, sec] / mn.prop[nc, 'TC']) ** 3 + + mn.prop[nc, 'vpD'] + * (1 - mn.Tnmin[col, sec] / mn.prop[nc, 'TC']) ** 6 + ) + ) + / Pnmin[sec] + for nc in b + ) + == 1 + ) + + mn.OBJ = Objective(expr=1, sense=minimize) + + #### + + SolverFactory('ipopt').solve(mn) + + yc = {} # Vapor composition + kl = {} # Light key component + kh = {} # Heavy key component + alpha = {} # Relative volatility of kl + ter = {} # Term to calculate the minimum number of trays + Nmin = {} # Minimum number of stages + Nminopt = {} # Total optimal minimum number of trays + Nfeed = {} # Side feed optimal location using Kirkbride's method: + # 1 = number of trays in rectifying section and + # 2 = number of trays in stripping section + side_feed = {} # Side feed location + av_alpha = {} # Average relative volatilities + xi_lhc = {} # Liquid composition in columns + rel = mn.Bdes / mn.Ddes # Ratio between products flowrates + ln = {} # Light component for the different columns + hn = {} # Heavy component for the different columns + ln[1] = lc_c1 + ln[2] = lc_c2 + ln[3] = lc_c3 + hn[1] = hc_c1 + hn[2] = hc_c2 + hn[3] = hc_c3 + + for col in mn.cols: + if col == 1: + b = mn.nc1 + elif col == 2: + b = mn.nc2 + else: + b = mn.nc3 + # For each component in the column and section calculate the vapor composition with the Peng-Robinson equation of state + for sec in mn.sec: + for nc in b: + yc[col, sec, nc] = ( + xi_nmin[col, sec, nc] + * mn.prop[nc, 'PC'] + * exp( + mn.prop[nc, 'TC'] + / value(mn.Tnmin[col, sec]) + * ( + mn.prop[nc, 'vpA'] + * (1 - value(mn.Tnmin[col, sec]) / mn.prop[nc, 'TC']) + + mn.prop[nc, 'vpB'] + * (1 - value(mn.Tnmin[col, sec]) / mn.prop[nc, 'TC']) ** 1.5 + + mn.prop[nc, 'vpC'] + * (1 - value(mn.Tnmin[col, sec]) / mn.prop[nc, 'TC']) ** 3 + + mn.prop[nc, 'vpD'] + * (1 - value(mn.Tnmin[col, sec]) / mn.prop[nc, 'TC']) ** 6 + ) + ) + ) / Pnmin[ + sec + ] # Vapor composition in the different sections for the different components in the columns + + for col in mn.cols: + # Calculate the relative volatility of the light and heavy components in the different sections for the different columns + xi_lhc[col, 4] = ( + xi_nmin[col, 1, ln[col]] / xi_nmin[col, 3, hn[col]] + ) # Liquid composition in the different sections with the initial liquid composition of the components in the different sections and columns and ln and hn which are the light and heavy components in the different columns + for sec in mn.sec: + kl[col, sec] = ( + yc[col, sec, ln[col]] / xi_nmin[col, sec, ln[col]] + ) # Light component in the different sections + kh[col, sec] = ( + yc[col, sec, hn[col]] / xi_nmin[col, sec, hn[col]] + ) # Heavy component in the different sections + xi_lhc[col, sec] = ( + xi_nmin[col, sec, hn[col]] / xi_nmin[col, sec, ln[col]] + ) # Liquid composition in the different sections + alpha[col, sec] = ( + kl[col, sec] / kh[col, sec] + ) # Relative volatility in the different sections + + for col in mn.cols: + # Calculate the average relative volatilities and the minimum number of trays with Fenske's and Kirkbride's method + av_alpha[col] = (alpha[col, 1] * alpha[col, 2] * alpha[col, 3]) ** ( + 1 / 3 + ) # Average relative volatilities calculated with the relative volatilities of the components in the three sections + Nmin[col] = log10((1 / xi_lhc[col, 3]) * xi_lhc[col, 1]) / log10( + av_alpha[col] + ) # Minimum number of trays calculated with Fenske's method + ter[col] = ( + rel * xi_lhc[col, 2] * (xi_lhc[col, 4] ** 2) + ) ** 0.206 # Term to calculate the minimum number of trays with Kirkbride's method + # Side feed optimal location using Kirkbride's method + Nfeed[1, col] = ( + ter[col] * Nmin[col] / (1 + ter[col]) + ) # Number of trays in rectifying section + Nfeed[2, col] = Nfeed[1, col] / ter[col] # Number of trays in stripping section + side_feed[col] = Nfeed[2, col] # Side feed location + + m.Nmintot = sum(Nmin[col] for col in mn.cols) # Total minimum number of trays + m.Knmin = int(m.Nmintot) + 1 # Total optimal minimum number of trays + + m.TB0 = value(mn.Tnmin[1, 1]) # Reboiler temperature in K in column 1 + m.Tf0 = value(mn.Tnmin[1, 2]) # Side feed temperature in K in column 1 + m.TD0 = value(mn.Tnmin[2, 3]) # Distillate temperature in K in column 2 + + return m + + +if __name__ == "__main__": + initialize_kaibel() diff --git a/gdplib/kaibel/kaibel_prop.py b/gdplib/kaibel/kaibel_prop.py index afae66f0..cbc51a30 100644 --- a/gdplib/kaibel/kaibel_prop.py +++ b/gdplib/kaibel/kaibel_prop.py @@ -90,7 +90,7 @@ def get_model_with_properties(): # ------------------------------------------------------------------ m.prop = {} # Properties of components: - cpL = {} # Ruczika-D method for liquid heat capacity calculation + cpL = {} # Ruzicka-Domalski method for liquid heat capacity calculation # (Reference A, page 6.20) sumA = {} sumB = {} diff --git a/gdplib/kaibel/kaibel_solve_gdp.py b/gdplib/kaibel/kaibel_solve_gdp.py index 5426d50b..385a99a7 100644 --- a/gdplib/kaibel/kaibel_solve_gdp.py +++ b/gdplib/kaibel/kaibel_solve_gdp.py @@ -1,1482 +1,1482 @@ """Kaibel Column model: GDP formulation. - -The solution requires the specification of certain parameters, such as the number trays, feed location, etc., and an initialization procedure, which consists of the next three steps: + +The solution requires the specification of certain parameters, such as the number trays, feed location, etc., and an initialization procedure, which consists of the next three steps: (i) a preliminary design of the separation considering a sequence of indirect continuous distillation columns (CDCs) to obtain the minimum number of stages with Fenske Equation in the function initialize_kaibel in kaibel_init.py -(ii) flash calculation for the feed with the function calc_side_feed_flash in kaibel_side_flash.py -(iii) calculation of variable bounds by solving the NLP problem. - -After the initialization, the GDP model is built. -""" - -from math import copysign - -from pyomo.environ import ( - Constraint, - exp, - minimize, - NonNegativeReals, - Objective, - RangeSet, - Set, - Var, -) -from pyomo.gdp import Disjunct - -from gdplib.kaibel.kaibel_init import initialize_kaibel - -from gdplib.kaibel.kaibel_side_flash import calc_side_feed_flash - - -def build_model(): - """ - Build the GDP Kaibel Column model. - It combines the initialization of the model and the flash calculation for the side feed before the GDP formulation. - - Returns - ------- - ConcreteModel - The constructed GDP Kaibel Column model. - """ - - # Calculation of the theoretical minimum number of trays (Knmin) and initial temperature values (TB0, Tf0, TD0). - m = initialize_kaibel() - - # Side feed init. Returns side feed vapor composition yfi and vapor fraction q_init - m = calc_side_feed_flash(m) - - m.name = "GDP Kaibel Column" - - #### Calculated initial values - m.Treb = m.TB0 + 5 # Reboiler temperature [K] - m.Tbot = m.TB0 # Bottom-most tray temperature [K] - m.Ttop = m.TD0 # Top-most tray temperature [K] - m.Tcon = m.TD0 - 5 # Condenser temperature [K] - - m.dv0 = {} # Initial vapor distributor value - m.dl0 = {} # Initial liquid distributor value - m.dv0[2] = 0.516 - m.dv0[3] = 1 - m.dv0[2] - m.dl0[2] = 0.36 - m.dl0[3] = 1 - m.dl0[2] - - #### Calculated upper and lower bounds - m.min_tray = m.Knmin * 0.8 # Lower bound on number of trays - m.Tlo = m.Tcon - 20 # Temperature lower bound - m.Tup = m.Treb + 20 # Temperature upper bound - - m.flow_max = 1e3 # Flowrates upper bound [mol/s] - m.Qmax = 60 # Heat loads upper bound [J/s] - - #### Column tray details - m.num_trays = m.np # Trays per section. np = 25 - m.min_num_trays = 10 # Minimum number of trays per section - m.num_total = m.np * 3 # Total number of trays - m.feed_tray = 12 # Side feed tray - m.sideout1_tray = 8 # Side outlet 1 tray - m.sideout2_tray = 17 # Side outlet 2 tray - m.reb_tray = 1 # Reboiler tray. Dividing wall starting tray - m.con_tray = m.num_trays # Condenser tray. Dividing wall ending tray - - # ------------------------------------------------------------------ - - # Beginning of model - - # ------------------------------------------------------------------ - - ## Sets - m.section = RangeSet( - 4, doc="Column sections:1=top, 2=feed side, 3=prod side, 4=bot" - ) - m.section_main = Set(initialize=[1, 4], doc="Main sections of the column") - - m.tray = RangeSet(m.np, doc="Potential trays in each section") - m.tray_total = RangeSet(m.num_total, doc="Total trays in the column") - m.tray_below_feed = RangeSet(m.feed_tray, doc="Trays below feed") - m.tray_below_so1 = RangeSet(m.sideout1_tray, doc="Trays below side outlet 1") - m.tray_below_so2 = RangeSet(m.sideout2_tray, doc="Trays below side outlet 2") - - m.comp = RangeSet(4, doc="Components") - m.dw = RangeSet(2, 3, doc="Dividing wall sections") - m.cplv = RangeSet(2, doc="Heat capacity: 1=liquid, 2=vapor") - m.so = RangeSet(2, doc="Side product outlets") - m.bounds = RangeSet(2, doc="Number of boundary condition values") - +(ii) flash calculation for the feed with the function calc_side_feed_flash in kaibel_side_flash.py +(iii) calculation of variable bounds by solving the NLP problem. + +After the initialization, the GDP model is built. +""" + +from math import copysign + +from pyomo.environ import ( + Constraint, + exp, + minimize, + NonNegativeReals, + Objective, + RangeSet, + Set, + Var, +) +from pyomo.gdp import Disjunct + +from gdplib.kaibel.kaibel_init import initialize_kaibel + +from gdplib.kaibel.kaibel_side_flash import calc_side_feed_flash + + +def build_model(): + """ + Build the GDP Kaibel Column model. + It combines the initialization of the model and the flash calculation for the side feed before the GDP formulation. + + Returns + ------- + ConcreteModel + The constructed GDP Kaibel Column model. + """ + + # Calculation of the theoretical minimum number of trays (Knmin) and initial temperature values (TB0, Tf0, TD0). + m = initialize_kaibel() + + # Side feed init. Returns side feed vapor composition yfi and vapor fraction q_init + m = calc_side_feed_flash(m) + + m.name = "GDP Kaibel Column" + + #### Calculated initial values + m.Treb = m.TB0 + 5 # Reboiler temperature [K] + m.Tbot = m.TB0 # Bottom-most tray temperature [K] + m.Ttop = m.TD0 # Top-most tray temperature [K] + m.Tcon = m.TD0 - 5 # Condenser temperature [K] + + m.dv0 = {} # Initial vapor distributor value + m.dl0 = {} # Initial liquid distributor value + m.dv0[2] = 0.516 + m.dv0[3] = 1 - m.dv0[2] + m.dl0[2] = 0.36 + m.dl0[3] = 1 - m.dl0[2] + + #### Calculated upper and lower bounds + m.min_tray = m.Knmin * 0.8 # Lower bound on number of trays + m.Tlo = m.Tcon - 20 # Temperature lower bound + m.Tup = m.Treb + 20 # Temperature upper bound + + m.flow_max = 1e3 # Flowrates upper bound [mol/s] + m.Qmax = 60 # Heat loads upper bound [J/s] + + #### Column tray details + m.num_trays = m.np # Trays per section. np = 25 + m.min_num_trays = 10 # Minimum number of trays per section + m.num_total = m.np * 3 # Total number of trays + m.feed_tray = 12 # Side feed tray + m.sideout1_tray = 8 # Side outlet 1 tray + m.sideout2_tray = 17 # Side outlet 2 tray + m.reb_tray = 1 # Reboiler tray. Dividing wall starting tray + m.con_tray = m.num_trays # Condenser tray. Dividing wall ending tray + + # ------------------------------------------------------------------ + + # Beginning of model + + # ------------------------------------------------------------------ + + ## Sets + m.section = RangeSet( + 4, doc="Column sections:1=top, 2=feed side, 3=prod side, 4=bot" + ) + m.section_main = Set(initialize=[1, 4], doc="Main sections of the column") + + m.tray = RangeSet(m.np, doc="Potential trays in each section") + m.tray_total = RangeSet(m.num_total, doc="Total trays in the column") + m.tray_below_feed = RangeSet(m.feed_tray, doc="Trays below feed") + m.tray_below_so1 = RangeSet(m.sideout1_tray, doc="Trays below side outlet 1") + m.tray_below_so2 = RangeSet(m.sideout2_tray, doc="Trays below side outlet 2") + + m.comp = RangeSet(4, doc="Components") + m.dw = RangeSet(2, 3, doc="Dividing wall sections") + m.cplv = RangeSet(2, doc="Heat capacity: 1=liquid, 2=vapor") + m.so = RangeSet(2, doc="Side product outlets") + m.bounds = RangeSet(2, doc="Number of boundary condition values") + m.candidate_trays_main = Set( initialize=m.tray - [m.con_tray, m.reb_tray], doc="Candidate trays for top and \ bottom sections 1 and 4" ) - m.candidate_trays_feed = Set( - initialize=m.tray - [m.con_tray, m.feed_tray, m.reb_tray], - doc="Candidate trays for feed section 2", - ) - m.candidate_trays_product = Set( - initialize=m.tray - [m.con_tray, m.sideout1_tray, m.sideout2_tray, m.reb_tray], - doc="Candidate trays for product section 3", - ) - - ## Calculation of initial values - m.dHvap = {} # Heat of vaporization [J/mol] - - m.P0 = {} # Initial pressure [bar] - m.T0 = {} # Initial temperature [K] - m.L0 = {} # Initial individual liquid flowrate [mol/s] - m.V0 = {} # Initial individual vapor flowrate [mol/s] - m.Vtotal0 = {} # Initial total vapor flowrate [mol/s] - m.Ltotal0 = {} # Initial liquid flowrate [mol/s] - m.x0 = {} # Initial liquid composition - m.y0 = {} # Initial vapor composition - m.actv0 = {} # Initial activity coefficients - m.cpdT0 = {} # Initial heat capacity for liquid and vapor phases [J/mol/K] - m.hl0 = {} # Initial liquid enthalpy [J/mol] - m.hv0 = {} # Initial vapor enthalpy [J/mol] - m.Pi = m.Preb # Initial given pressure value [bar] - m.Ti = {} # Initial known temperature values [K] - - ## Initial values for pressure, temperature, flowrates, composition, and enthalpy - for sec in m.section: - for n_tray in m.tray: - m.P0[sec, n_tray] = m.Pi - - for sec in m.section: - for n_tray in m.tray: - for comp in m.comp: - m.L0[sec, n_tray, comp] = m.Li - m.V0[sec, n_tray, comp] = m.Vi - - for sec in m.section: - for n_tray in m.tray: - m.Ltotal0[sec, n_tray] = sum(m.L0[sec, n_tray, comp] for comp in m.comp) - m.Vtotal0[sec, n_tray] = sum(m.V0[sec, n_tray, comp] for comp in m.comp) - - for n_tray in m.tray_total: - if n_tray == m.reb_tray: - m.Ti[n_tray] = m.Treb - elif n_tray == m.num_total: - m.Ti[n_tray] = m.Tcon - else: - m.Ti[n_tray] = m.Tbot + (m.Ttop - m.Tbot) * (n_tray - 2) / (m.num_total - 3) - - for n_tray in m.tray_total: - if n_tray <= m.num_trays: - m.T0[1, n_tray] = m.Ti[n_tray] - elif n_tray >= m.num_trays and n_tray <= m.num_trays * 2: - m.T0[2, n_tray - m.num_trays] = m.Ti[n_tray] - m.T0[3, n_tray - m.num_trays] = m.Ti[n_tray] - elif n_tray >= m.num_trays * 2: - m.T0[4, n_tray - m.num_trays * 2] = m.Ti[n_tray] - - ## Initial vapor and liquid composition of the feed and activity coefficients - for sec in m.section: - for n_tray in m.tray: - for comp in m.comp: - m.x0[sec, n_tray, comp] = m.xfi[comp] - m.actv0[sec, n_tray, comp] = 1 - m.y0[sec, n_tray, comp] = m.xfi[comp] - - ## Assigns the enthalpy boundary values, heat capacity, heat of vaporization calculation, temperature bounds, and light and heavy key components. - hlb = {} # Liquid enthalpy [J/mol] - hvb = {} # Vapor enthalpy [J/mol] - cpb = {} # Heact capacity [J/mol/K] - dHvapb = {} # Heat of vaporization [J/mol] - Tbounds = {} # Temperature bounds [K] - kc = {} # Light and heavy key components - Tbounds[1] = m.Tcon # Condenser temperature [K] - Tbounds[2] = m.Treb # Reboiler temperature [K] - kc[1] = m.lc # Light key component - kc[2] = m.hc # Heavy key component - - ## Heat of vaporization calculation for each component in the feed. - for comp in m.comp: - dHvapb[comp] = -( - m.Rgas - * m.prop[comp, 'TC'] - * ( - m.prop[comp, 'vpA'] * (1 - m.Tref / m.prop[comp, 'TC']) - + m.prop[comp, 'vpB'] * (1 - m.Tref / m.prop[comp, 'TC']) ** 1.5 - + m.prop[comp, 'vpC'] * (1 - m.Tref / m.prop[comp, 'TC']) ** 3 - + m.prop[comp, 'vpD'] * (1 - m.Tref / m.prop[comp, 'TC']) ** 6 - ) - + m.Rgas - * m.Tref - * ( - m.prop[comp, 'vpA'] - + 1.5 * m.prop[comp, 'vpB'] * (1 - m.Tref / m.prop[comp, 'TC']) ** 0.5 - + 3 * m.prop[comp, 'vpC'] * (1 - m.Tref / m.prop[comp, 'TC']) ** 2 - + 6 * m.prop[comp, 'vpD'] * (1 - m.Tref / m.prop[comp, 'TC']) ** 5 - ) - ) - - ## Boundary values for heat capacity and enthalpy of liquid and vapor phases for light and heavy key components in the feed. - for b in m.bounds: - for cp in m.cplv: - cpb[b, cp] = m.cpc[cp] * ( - (Tbounds[b] - m.Tref) * m.prop[kc[b], 'cpA', cp] - + (Tbounds[b] ** 2 - m.Tref**2) - * m.prop[kc[b], 'cpB', cp] - * m.cpc2['A', cp] - / 2 - + (Tbounds[b] ** 3 - m.Tref**3) - * m.prop[kc[b], 'cpC', cp] - * m.cpc2['B', cp] - / 3 - + (Tbounds[b] ** 4 - m.Tref**4) * m.prop[kc[b], 'cpD', cp] / 4 - ) - hlb[b] = cpb[b, 1] - hvb[b] = cpb[b, 2] + dHvapb[b] - - m.hllo = ( - (1 - copysign(0.2, hlb[1])) * hlb[1] / m.Hscale - ) # Liquid enthalpy lower bound - m.hlup = ( - (1 + copysign(0.2, hlb[2])) * hlb[2] / m.Hscale - ) # Liquid enthalpy upper bound - m.hvlo = ( - (1 - copysign(0.2, hvb[1])) * hvb[1] / m.Hscale - ) # Vapor enthalpy lower bound - m.hvup = ( - (1 + copysign(0.2, hvb[2])) * hvb[2] / m.Hscale - ) # Vapor enthalpy upper bound - # copysign is a function that returns the first argument with the sign of the second argument - - ## Heat of vaporization for each component in the feed scaled by Hscale - for comp in m.comp: - m.dHvap[comp] = dHvapb[comp] / m.Hscale - - ## Heat capacity calculation using the Ruczika-D method for each component in the feed, section, and tray - for sec in m.section: - for n_tray in m.tray: - for comp in m.comp: - for cp in m.cplv: - m.cpdT0[sec, n_tray, comp, cp] = ( - m.cpc[cp] - * ( - (m.T0[sec, n_tray] - m.Tref) * m.prop[comp, 'cpA', cp] - + (m.T0[sec, n_tray] ** 2 - m.Tref**2) - * m.prop[comp, 'cpB', cp] - * m.cpc2['A', cp] - / 2 - + (m.T0[sec, n_tray] ** 3 - m.Tref**3) - * m.prop[comp, 'cpC', cp] - * m.cpc2['B', cp] - / 3 - + (m.T0[sec, n_tray] ** 4 - m.Tref**4) - * m.prop[comp, 'cpD', cp] - / 4 - ) - / m.Hscale - ) - - ## Liquid and vapor enthalpy calculation using the Ruczika-D method for each component in the feed, section, and tray - for sec in m.section: - for n_tray in m.tray: - for comp in m.comp: - m.hl0[sec, n_tray, comp] = m.cpdT0[sec, n_tray, comp, 1] - m.hv0[sec, n_tray, comp] = m.cpdT0[sec, n_tray, comp, 2] + m.dHvap[comp] - - #### Side feed - m.cpdTf = {} # Heat capacity for side feed [J/mol K] - m.hlf = {} # Liquid enthalpy for side feed [J/mol] - m.hvf = {} # Vapor enthalpy for side feed [J/mol] - m.F0 = {} # Side feed flowrate per component [mol/s] - - ## Heat capacity in liquid and vapor phases for side feed using the Ruczika-D method. - for comp in m.comp: - for cp in m.cplv: - m.cpdTf[comp, cp] = ( - m.cpc[cp] - * ( - (m.Tf - m.Tref) * m.prop[comp, 'cpA', cp] - + (m.Tf**2 - m.Tref**2) - * m.prop[comp, 'cpB', cp] - * m.cpc2['A', cp] - / 2 - + (m.Tf**3 - m.Tref**3) - * m.prop[comp, 'cpC', cp] - * m.cpc2['B', cp] - / 3 - + (m.Tf**4 - m.Tref**4) * m.prop[comp, 'cpD', cp] / 4 - ) - / m.Hscale - ) - - ## Side feed flowrate and enthalpy calculation using the Ruczika-D method for each component in the feed - for comp in m.comp: - m.F0[comp] = ( - m.xfi[comp] * m.Fi - ) # Side feed flowrate per component computed from the feed composition and flowrate Fi - m.hlf[comp] = m.cpdTf[ - comp, 1 - ] # Liquid enthalpy for side feed computed from the heat capacity for side feed and liquid phase - m.hvf[comp] = ( - m.cpdTf[comp, 2] + m.dHvap[comp] - ) # Vapor enthalpy for side feed computed from the heat capacity for side feed and vapor phase and heat of vaporization - - m.P = Var( - m.section, - m.tray, - doc="Pressure at each potential tray [bar]", - domain=NonNegativeReals, - bounds=(m.Pcon, m.Preb), - initialize=m.P0, - ) - m.T = Var( - m.section, - m.tray, - doc="Temperature at each potential tray [K]", - domain=NonNegativeReals, - bounds=(m.Tlo, m.Tup), - initialize=m.T0, - ) - - m.Tr = Var( - m.section, - m.tray, - m.comp, - doc='Temperature term for vapor pressure [-]', - domain=NonNegativeReals, - bounds=lambda m, sec, tray, comp: ( - m.prop[comp, 'TC'] / m.Tup, - m.prop[comp, 'TC'] / m.Tlo, + m.candidate_trays_feed = Set( + initialize=m.tray - [m.con_tray, m.feed_tray, m.reb_tray], + doc="Candidate trays for feed section 2", + ) + m.candidate_trays_product = Set( + initialize=m.tray - [m.con_tray, m.sideout1_tray, m.sideout2_tray, m.reb_tray], + doc="Candidate trays for product section 3", + ) + + ## Calculation of initial values + m.dHvap = {} # Heat of vaporization [J/mol] + + m.P0 = {} # Initial pressure [bar] + m.T0 = {} # Initial temperature [K] + m.L0 = {} # Initial individual liquid flowrate [mol/s] + m.V0 = {} # Initial individual vapor flowrate [mol/s] + m.Vtotal0 = {} # Initial total vapor flowrate [mol/s] + m.Ltotal0 = {} # Initial liquid flowrate [mol/s] + m.x0 = {} # Initial liquid composition + m.y0 = {} # Initial vapor composition + m.actv0 = {} # Initial activity coefficients + m.cpdT0 = {} # Initial heat capacity for liquid and vapor phases [J/mol/K] + m.hl0 = {} # Initial liquid enthalpy [J/mol] + m.hv0 = {} # Initial vapor enthalpy [J/mol] + m.Pi = m.Preb # Initial given pressure value [bar] + m.Ti = {} # Initial known temperature values [K] + + ## Initial values for pressure, temperature, flowrates, composition, and enthalpy + for sec in m.section: + for n_tray in m.tray: + m.P0[sec, n_tray] = m.Pi + + for sec in m.section: + for n_tray in m.tray: + for comp in m.comp: + m.L0[sec, n_tray, comp] = m.Li + m.V0[sec, n_tray, comp] = m.Vi + + for sec in m.section: + for n_tray in m.tray: + m.Ltotal0[sec, n_tray] = sum(m.L0[sec, n_tray, comp] for comp in m.comp) + m.Vtotal0[sec, n_tray] = sum(m.V0[sec, n_tray, comp] for comp in m.comp) + + for n_tray in m.tray_total: + if n_tray == m.reb_tray: + m.Ti[n_tray] = m.Treb + elif n_tray == m.num_total: + m.Ti[n_tray] = m.Tcon + else: + m.Ti[n_tray] = m.Tbot + (m.Ttop - m.Tbot) * (n_tray - 2) / (m.num_total - 3) + + for n_tray in m.tray_total: + if n_tray <= m.num_trays: + m.T0[1, n_tray] = m.Ti[n_tray] + elif n_tray >= m.num_trays and n_tray <= m.num_trays * 2: + m.T0[2, n_tray - m.num_trays] = m.Ti[n_tray] + m.T0[3, n_tray - m.num_trays] = m.Ti[n_tray] + elif n_tray >= m.num_trays * 2: + m.T0[4, n_tray - m.num_trays * 2] = m.Ti[n_tray] + + ## Initial vapor and liquid composition of the feed and activity coefficients + for sec in m.section: + for n_tray in m.tray: + for comp in m.comp: + m.x0[sec, n_tray, comp] = m.xfi[comp] + m.actv0[sec, n_tray, comp] = 1 + m.y0[sec, n_tray, comp] = m.xfi[comp] + + ## Assigns the enthalpy boundary values, heat capacity, heat of vaporization calculation, temperature bounds, and light and heavy key components. + hlb = {} # Liquid enthalpy [J/mol] + hvb = {} # Vapor enthalpy [J/mol] + cpb = {} # Heact capacity [J/mol/K] + dHvapb = {} # Heat of vaporization [J/mol] + Tbounds = {} # Temperature bounds [K] + kc = {} # Light and heavy key components + Tbounds[1] = m.Tcon # Condenser temperature [K] + Tbounds[2] = m.Treb # Reboiler temperature [K] + kc[1] = m.lc # Light key component + kc[2] = m.hc # Heavy key component + + ## Heat of vaporization calculation for each component in the feed. + for comp in m.comp: + dHvapb[comp] = -( + m.Rgas + * m.prop[comp, 'TC'] + * ( + m.prop[comp, 'vpA'] * (1 - m.Tref / m.prop[comp, 'TC']) + + m.prop[comp, 'vpB'] * (1 - m.Tref / m.prop[comp, 'TC']) ** 1.5 + + m.prop[comp, 'vpC'] * (1 - m.Tref / m.prop[comp, 'TC']) ** 3 + + m.prop[comp, 'vpD'] * (1 - m.Tref / m.prop[comp, 'TC']) ** 6 + ) + + m.Rgas + * m.Tref + * ( + m.prop[comp, 'vpA'] + + 1.5 * m.prop[comp, 'vpB'] * (1 - m.Tref / m.prop[comp, 'TC']) ** 0.5 + + 3 * m.prop[comp, 'vpC'] * (1 - m.Tref / m.prop[comp, 'TC']) ** 2 + + 6 * m.prop[comp, 'vpD'] * (1 - m.Tref / m.prop[comp, 'TC']) ** 5 + ) + ) + + ## Boundary values for heat capacity and enthalpy of liquid and vapor phases for light and heavy key components in the feed. + for b in m.bounds: + for cp in m.cplv: + cpb[b, cp] = m.cpc[cp] * ( + (Tbounds[b] - m.Tref) * m.prop[kc[b], 'cpA', cp] + + (Tbounds[b] ** 2 - m.Tref**2) + * m.prop[kc[b], 'cpB', cp] + * m.cpc2['A', cp] + / 2 + + (Tbounds[b] ** 3 - m.Tref**3) + * m.prop[kc[b], 'cpC', cp] + * m.cpc2['B', cp] + / 3 + + (Tbounds[b] ** 4 - m.Tref**4) * m.prop[kc[b], 'cpD', cp] / 4 + ) + hlb[b] = cpb[b, 1] + hvb[b] = cpb[b, 2] + dHvapb[b] + + m.hllo = ( + (1 - copysign(0.2, hlb[1])) * hlb[1] / m.Hscale + ) # Liquid enthalpy lower bound + m.hlup = ( + (1 + copysign(0.2, hlb[2])) * hlb[2] / m.Hscale + ) # Liquid enthalpy upper bound + m.hvlo = ( + (1 - copysign(0.2, hvb[1])) * hvb[1] / m.Hscale + ) # Vapor enthalpy lower bound + m.hvup = ( + (1 + copysign(0.2, hvb[2])) * hvb[2] / m.Hscale + ) # Vapor enthalpy upper bound + # copysign is a function that returns the first argument with the sign of the second argument + + ## Heat of vaporization for each component in the feed scaled by Hscale + for comp in m.comp: + m.dHvap[comp] = dHvapb[comp] / m.Hscale + + ## Heat capacity calculation using the Ruzicka-Domalski method for each component in the feed, section, and tray + for sec in m.section: + for n_tray in m.tray: + for comp in m.comp: + for cp in m.cplv: + m.cpdT0[sec, n_tray, comp, cp] = ( + m.cpc[cp] + * ( + (m.T0[sec, n_tray] - m.Tref) * m.prop[comp, 'cpA', cp] + + (m.T0[sec, n_tray] ** 2 - m.Tref**2) + * m.prop[comp, 'cpB', cp] + * m.cpc2['A', cp] + / 2 + + (m.T0[sec, n_tray] ** 3 - m.Tref**3) + * m.prop[comp, 'cpC', cp] + * m.cpc2['B', cp] + / 3 + + (m.T0[sec, n_tray] ** 4 - m.Tref**4) + * m.prop[comp, 'cpD', cp] + / 4 + ) + / m.Hscale + ) + + ## Liquid and vapor enthalpy calculation using the Ruzicka-Domalski method for each component in the feed, section, and tray + for sec in m.section: + for n_tray in m.tray: + for comp in m.comp: + m.hl0[sec, n_tray, comp] = m.cpdT0[sec, n_tray, comp, 1] + m.hv0[sec, n_tray, comp] = m.cpdT0[sec, n_tray, comp, 2] + m.dHvap[comp] + + #### Side feed + m.cpdTf = {} # Heat capacity for side feed [J/mol K] + m.hlf = {} # Liquid enthalpy for side feed [J/mol] + m.hvf = {} # Vapor enthalpy for side feed [J/mol] + m.F0 = {} # Side feed flowrate per component [mol/s] + + ## Heat capacity in liquid and vapor phases for side feed using the Ruzicka-Domalski method. + for comp in m.comp: + for cp in m.cplv: + m.cpdTf[comp, cp] = ( + m.cpc[cp] + * ( + (m.Tf - m.Tref) * m.prop[comp, 'cpA', cp] + + (m.Tf**2 - m.Tref**2) + * m.prop[comp, 'cpB', cp] + * m.cpc2['A', cp] + / 2 + + (m.Tf**3 - m.Tref**3) + * m.prop[comp, 'cpC', cp] + * m.cpc2['B', cp] + / 3 + + (m.Tf**4 - m.Tref**4) * m.prop[comp, 'cpD', cp] / 4 + ) + / m.Hscale + ) + + ## Side feed flowrate and enthalpy calculation using the Ruzicka-Domalski method for each component in the feed + for comp in m.comp: + m.F0[comp] = ( + m.xfi[comp] * m.Fi + ) # Side feed flowrate per component computed from the feed composition and flowrate Fi + m.hlf[comp] = m.cpdTf[ + comp, 1 + ] # Liquid enthalpy for side feed computed from the heat capacity for side feed and liquid phase + m.hvf[comp] = ( + m.cpdTf[comp, 2] + m.dHvap[comp] + ) # Vapor enthalpy for side feed computed from the heat capacity for side feed and vapor phase and heat of vaporization + + m.P = Var( + m.section, + m.tray, + doc="Pressure at each potential tray [bar]", + domain=NonNegativeReals, + bounds=(m.Pcon, m.Preb), + initialize=m.P0, + ) + m.T = Var( + m.section, + m.tray, + doc="Temperature at each potential tray [K]", + domain=NonNegativeReals, + bounds=(m.Tlo, m.Tup), + initialize=m.T0, + ) + + m.Tr = Var( + m.section, + m.tray, + m.comp, + doc='Temperature term for vapor pressure [-]', + domain=NonNegativeReals, + bounds=lambda m, sec, tray, comp: ( + m.prop[comp, 'TC'] / m.Tup, + m.prop[comp, 'TC'] / m.Tlo, ), initialize=lambda m, sec, tray, comp: m.prop[comp, 'TC'] / m.T0[sec, tray], - ) - - m.x = Var( - m.section, - m.tray, - m.comp, - doc="Liquid composition [mol/mol]", - domain=NonNegativeReals, - bounds=(0, 1), - initialize=m.x0, - ) - m.y = Var( - m.section, - m.tray, - m.comp, - doc="Vapor composition [mol/mol]", - domain=NonNegativeReals, - bounds=(0, 1), - initialize=m.y0, - ) - - m.dl = Var( - m.dw, - doc="Liquid distributor in the dividing wall sections", - bounds=(0.2, 0.8), - initialize=m.dl0, - ) - m.dv = Var( - m.dw, - doc="Vapor distributor in the dividing wall sections", - bounds=(0, 1), - domain=NonNegativeReals, - initialize=m.dv0, - ) - - m.V = Var( - m.section, - m.tray, - m.comp, - doc="Vapor flowrate [mol/s]", - domain=NonNegativeReals, - bounds=(0, m.flow_max), - initialize=m.V0, - ) - m.L = Var( - m.section, - m.tray, - m.comp, - doc="Liquid flowrate [mol/s]", - domain=NonNegativeReals, - bounds=(0, m.flow_max), - initialize=m.L0, - ) - m.Vtotal = Var( - m.section, - m.tray, - doc="Total vapor flowrate [mol/s]", - domain=NonNegativeReals, - bounds=(0, m.flow_max), - initialize=m.Vtotal0, - ) - m.Ltotal = Var( - m.section, - m.tray, - doc="Total liquid flowrate [mol/s]", - domain=NonNegativeReals, - bounds=(0, m.flow_max), - initialize=m.Ltotal0, - ) - - m.D = Var( - m.comp, - doc="Distillate flowrate [mol/s]", - domain=NonNegativeReals, - bounds=(0, m.flow_max), - initialize=m.Ddes, - ) - m.B = Var( - m.comp, - doc="Bottoms flowrate [mol/s]", - domain=NonNegativeReals, - bounds=(0, m.flow_max), - initialize=m.Bdes, - ) - m.S = Var( - m.so, - m.comp, - doc="Product 2 and 3 flowrates [mol/s]", - domain=NonNegativeReals, - bounds=(0, m.flow_max), - initialize=m.Sdes, - ) - m.Dtotal = Var( - doc="Distillate flowrate [mol/s]", - domain=NonNegativeReals, - bounds=(0, m.flow_max), - initialize=m.Ddes, - ) - m.Btotal = Var( - doc="Bottoms flowrate [mol/s]", - domain=NonNegativeReals, - bounds=(0, m.flow_max), - initialize=m.Bdes, - ) - m.Stotal = Var( - m.so, - doc="Total product 2 and 3 side flowrate [mol/s]", - domain=NonNegativeReals, - bounds=(0, m.flow_max), - initialize=m.Sdes, - ) - - m.hl = Var( - m.section, - m.tray, - m.comp, - doc='Liquid enthalpy [J/mol]', - bounds=(m.hllo, m.hlup), - initialize=m.hl0, - ) - m.hv = Var( - m.section, - m.tray, - m.comp, - doc='Vapor enthalpy [J/mol]', - bounds=(m.hvlo, m.hvup), - initialize=m.hv0, - ) - m.Qreb = Var( - doc="Reboiler heat duty [J/s]", - domain=NonNegativeReals, - bounds=(0, m.Qmax), - initialize=1, - ) - m.Qcon = Var( - doc="Condenser heat duty [J/s]", - domain=NonNegativeReals, - bounds=(0, m.Qmax), - initialize=1, - ) - - m.rr = Var( - doc="Internal reflux ratio in the column", - domain=NonNegativeReals, - bounds=(0.7, 1), - initialize=m.rr0, - ) - m.bu = Var( - doc="Boilup rate in the reboiler", - domain=NonNegativeReals, - bounds=(0.7, 1), - initialize=m.bu0, - ) - - m.F = Var( - m.comp, - doc="Side feed flowrate [mol/s]", - domain=NonNegativeReals, - bounds=(0, 50), - initialize=m.F0, - ) - m.q = Var( - doc="Vapor fraction in side feed", - domain=NonNegativeReals, - bounds=(0, 1), - initialize=1, - ) - - m.actv = Var( - m.section, - m.tray, - m.comp, - doc="Liquid activity coefficient", - domain=NonNegativeReals, - bounds=(0, 10), - initialize=m.actv0, - ) - - m.errx = Var( - m.section, - m.tray, - doc="Error in liquid composition [mol/mol]", - bounds=(-1e-3, 1e-3), - initialize=0, - ) - m.erry = Var( - m.section, - m.tray, - doc="Error in vapor composition [mol/mol]", - bounds=(-1e-3, 1e-3), - initialize=0, - ) - m.slack = Var( - m.section, - m.tray, - m.comp, - doc="Slack variable", - bounds=(-1e-8, 1e-8), - initialize=0, - ) - - m.tray_exists = Disjunct( - m.section, - m.tray, - doc="Disjunct that enforce the existence of each tray", - rule=_build_tray_equations, - ) - m.tray_absent = Disjunct( - m.section, - m.tray, - doc="Disjunct that enforce the absence of each tray", - rule=_build_pass_through_eqns, - ) - - @m.Disjunction( - m.section, m.tray, doc="Disjunction between whether each tray exists or not" - ) - def tray_exists_or_not(m, sec, n_tray): - """ - Disjunction between whether each tray exists or not. - - Parameters - ---------- - m : Pyomo ConcreteModel - The Pyomo model object. - sec : int - The section index. - n_tray : int - The tray index. - - Returns - ------- - Disjunction - The disjunction between whether each tray exists or not. - """ - return [m.tray_exists[sec, n_tray], m.tray_absent[sec, n_tray]] - - @m.Constraint(m.section_main) - def minimum_trays_main(m, sec): - """ - Constraint that ensures the minimum number of trays in the main section. - - Parameters - ---------- - m : Pyomo ConcreteModel - The model object for the GDP Kaibel Column. - sec : Set - The section index. - - Returns - ------- - Constraint - A constraint expression that enforces the minimum number of trays in the main section to be greater than or equal to the minimum number of trays. - """ - return ( - sum( - m.tray_exists[sec, n_tray].binary_indicator_var - for n_tray in m.candidate_trays_main - ) - + 1 - >= m.min_num_trays - ) - - @m.Constraint() - def minimum_trays_feed(m): - """ - Constraint function that ensures the minimum number of trays in the feed section is met. - - Parameters - ---------- - m : Pyomo ConcreteModel - The Pyomo model object. - - Returns - ------- - Constraint - The constraint expression that enforces the minimum number of trays is greater than or equal to the minimum number of trays. - """ - return ( - sum( - m.tray_exists[2, n_tray].binary_indicator_var - for n_tray in m.candidate_trays_feed - ) - + 1 - >= m.min_num_trays - ) - - # TOCHECK: pyomo.GDP Syntax - - @m.Constraint() - def minimum_trays_product(m): - """ - Constraint function to calculate the minimum number of trays in the product section. - - Parameters - ---------- - m : Pyomo ConcreteModel - The optimization model. - - Returns - ------- - Constraint - The constraint expression that enforces the minimum number of trays is greater than or equal to the minimum number of trays. - """ - return ( - sum( - m.tray_exists[3, n_tray].binary_indicator_var - for n_tray in m.candidate_trays_product - ) - + 1 - >= m.min_num_trays - ) - - ## Fixed trays - enforce_tray_exists(m, 1, 1) # reboiler - enforce_tray_exists(m, 1, m.num_trays) # vapor distributor - enforce_tray_exists(m, 2, 1) # dividing wall starting tray - enforce_tray_exists(m, 2, m.feed_tray) # feed tray - enforce_tray_exists(m, 2, m.num_trays) # dividing wall ending tray - enforce_tray_exists(m, 3, 1) # dividing wall starting tray - enforce_tray_exists(m, 3, m.sideout1_tray) # side outlet 1 for product 3 - enforce_tray_exists(m, 3, m.sideout2_tray) # side outlet 2 for product 2 - enforce_tray_exists(m, 3, m.num_trays) # dividing wall ending tray - enforce_tray_exists(m, 4, 1) # liquid distributor - enforce_tray_exists(m, 4, m.num_trays) # condenser - - #### Global constraints - @m.Constraint( - m.dw, m.tray, doc="Monotonic temperature in the dividing wall sections" - ) - def monotonic_temperature(m, sec, n_tray): - """This function returns a constraint object representing the monotonic temperature constraint. - - The monotonic temperature constraint ensures that the temperature on each tray in the distillation column - is less than or equal to the temperature on the top tray. - - Parameters - ---------- - m : Pyomo ConcreteModel - The Pyomo model object. - sec : Set - The set of sections in the dividing wall sections. - n_tray : Set - The set of trays in the distillation column. - - Returns - ------- - Constraint - The monotonic temperature constraint specifying that the temperature on each tray is less than or equal to the temperature on the top tray of section 1 which is the condenser. - """ - return m.T[sec, n_tray] <= m.T[1, m.num_trays] - - @m.Constraint(doc="Liquid distributor") - def liquid_distributor(m): - """Defines the liquid distributor constraint. - - This constraint ensures that the sum of the liquid distributors in all sections is equal to 1. - - Parameters - ---------- - m : Pyomo ConcreteModel - The optimization model. - - Returns - ------- - Constraint - The liquid distributor constraint that enforces the sum of the liquid flow rates in all sections is equal to 1. - - """ - return sum(m.dl[sec] for sec in m.dw) - 1 == 0 - - @m.Constraint(doc="Vapor distributor") - def vapor_distributor(m): - """ - Add a constraint to ensure that the sum of the vapor distributors is equal to 1. - - Parameters - ---------- - m : Pyomo ConcreteModel - The Pyomo model object. - - Returns - ------- - Constraint - The vapor distributor constraint. - """ - return sum(m.dv[sec] for sec in m.dw) - 1 == 0 - - @m.Constraint(doc="Reboiler composition specification") - def heavy_product(m): - """ - Reboiler composition specification for the heavy component in the feed. - - Parameters - ---------- - m : Pyomo ConcreteModel - The optimization model. - - Returns - ------- - Constraint - The constraint that enforces the reboiler composition is greater than or equal to the specified composition xspechc final liquid composition for butanol, the heavy component in the feed. - """ - return m.x[1, m.reb_tray, m.hc] >= m.xspec_hc - - @m.Constraint(doc="Condenser composition specification") - def light_product(m): - """ - Condenser composition specification for the light component in the feed. - - Parameters - ---------- - m : Model - The optimization model. - - Returns - ------- - Constraint - The constraint that enforces the condenser composition is greater than or equal to the specified final liquid composition for ethanol, xspeclc , the light component in the feed. - """ - return m.x[4, m.con_tray, m.lc] >= m.xspec_lc - - @m.Constraint(doc="Side outlet 1 final liquid composition") - def intermediate1_product(m): - ''' - This constraint ensures that the intermediate 1 final liquid composition is greater than or equal to the specified composition xspec_inter1, which is the final liquid composition for ethanol. - - Parameters - ---------- - m : Pyomo ConcreteModel. - The optimization model. - - Returns - ------- - Constraint - The constraint that enforces the intermediate 1 final liquid composition is greater than or equal to the specified composition xspec_inter1, which is the final liquid composition for ethanol. - ''' - return m.x[3, m.sideout1_tray, 3] >= m.xspec_inter3 - - @m.Constraint(doc="Side outlet 2 final liquid composition") - def intermediate2_product(m): - """ - This constraint ensures that the intermediate 2 final liquid composition is greater than or equal to the specified composition xspec_inter2, which is the final liquid composition for butanol. - - Parameters - ---------- - m : Pyomo ConcreteModel - The optimization model. - - Returns - ------- - Constraint - The constraint that enforces the intermediate 2 final liquid composition is greater than or equal to the specified composition xspec_inter2, which is the final liquid composition for butanol. - """ - return m.x[3, m.sideout2_tray, 2] >= m.xspec_inter2 - - @m.Constraint(doc="Reboiler flowrate") - def _heavy_product_flow(m): - """ - Reboiler flowrate constraint that ensures the reboiler flowrate is greater than or equal to the specified flowrate Bdes, which is the flowrate of butanol. - - Parameters - ---------- - m : Pyomo ConcreteModel - The optimization model. - - Returns - ------- - Constraint - The constraint that enforces the reboiler flowrate is greater than or equal to the specified flowrate Bdes, which is the flowrate of butanol. - """ - return m.Btotal >= m.Bdes - - @m.Constraint(doc="Condenser flowrate") - def _light_product_flow(m): - """ - Condenser flowrate constraint that ensures the condenser flowrate is greater than or equal to the specified flowrate Ddes, which is the flowrate of ethanol. - - Parameters - ---------- - m : Pyomo ConcreteModel - The optimization model. - - Returns - ------- - Constraint - The constraint that enforces the condenser flowrate is greater than or equal to the specified flowrate Ddes, which is the flowrate of ethanol. - """ - return m.Dtotal >= m.Ddes - - @m.Constraint(m.so, doc="Intermediate flowrate") - def _intermediate_product_flow(m, so): - """ - Intermediate flowrate constraint that ensures the intermediate flowrate is greater than or equal to the specified flowrate Sdes, which is the flowrate of the intermediate side product 2 and 3. - - Parameters - ---------- - m : Pyomo ConcreteModel - The optimization model. - so : int - The side product outlet index. - - Returns - ------- - Constraint - The constraint that enforces the intermediate flowrate is greater than or equal to the specified flowrate Sdes, which is the flowrate of the intermediate side product 2 and 3. - """ - return m.Stotal[so] >= m.Sdes - - @m.Constraint(doc="Internal boilup ratio, V/L") - def _internal_boilup_ratio(m): - """ - Internal boilup ratio constraint that ensures the internal boilup ratio is equal to the boilup rate times the liquid flowrate on the reboiler tray is equal to the vapor flowrate on the tray above the reboiler tray. - - Parameters - ---------- - m : Pyomo ConcreteModel - The optimization model. - - Returns - ------- - Constraint - The constraint that enforces the boilup rate times the liquid flowrate on the reboiler tray is equal to the vapor flowrate on the tray above the reboiler tray. - """ - return m.bu * m.Ltotal[1, m.reb_tray + 1] == m.Vtotal[1, m.reb_tray] - - @m.Constraint(doc="Internal reflux ratio, L/V") - def internal_reflux_ratio(m): - """ - Internal reflux ratio constraint that ensures the internal reflux ratio is equal to the reflux rate times the vapor flowrate on the tray above the condenser tray is equal to the liquid flowrate on the condenser tray. - - Parameters - ---------- - m : Pyomo ConcreteModel - The optimization model. - - Returns - ------- - Constraint - The constraint that enforces the reflux rate times the vapor flowrate on the tray above the condenser tray is equal to the liquid flowrate on the condenser tray. - """ - return m.rr * m.Vtotal[4, m.con_tray - 1] == m.Ltotal[4, m.con_tray] - - @m.Constraint(doc="External boilup ratio relation with bottoms") - def _external_boilup_ratio(m): - """ - External boilup ratio constraint that ensures the external boilup ratio times the liquid flowrate on the reboiler tray is equal to the bottoms flowrate. - - Parameters - ---------- - m : Pyomo ConcreteModel - The optimization model. - - Returns - ------- - Constraint - The constraint that enforces the external boilup ratio times the liquid flowrate on the reboiler tray is equal to the bottoms flowrate. - """ - return m.Btotal == (1 - m.bu) * m.Ltotal[1, m.reb_tray + 1] - - @m.Constraint(doc="External reflux ratio relation with distillate") - def _external_reflux_ratio(m): - """ - External reflux ratio constraint that ensures the external reflux ratio times the vapor flowrate on the tray above the condenser tray is equal to the distillate flowrate. - - Parameters - ---------- - m : Pyomo ConcreteModel - The optimization model. - - Returns - ------- - Constraint - The constraint that enforces the external reflux ratio times the vapor flowrate on the tray above the condenser tray is equal to the distillate flowrate. - """ - return m.Dtotal == (1 - m.rr) * m.Vtotal[4, m.con_tray - 1] - - @m.Constraint(m.section, m.tray, doc="Total vapor flowrate") - def _total_vapor_flowrate(m, sec, n_tray): - """ - Constraint that ensures the total vapor flowrate is equal to the sum of the vapor flowrates of each component on each tray. - - Parameters - ---------- - m : Pyomo ConcreteModel - The optimization model. - sec : int - The section index. - n_tray : int - The tray index. - - Returns - ------- - Constraint - The constraint that enforces the total vapor flowrate is equal to the sum of the vapor flowrates of each component on each tray on each section. - """ - return sum(m.V[sec, n_tray, comp] for comp in m.comp) == m.Vtotal[sec, n_tray] - - @m.Constraint(m.section, m.tray, doc="Total liquid flowrate") - def _total_liquid_flowrate(m, sec, n_tray): - """ - Constraint that ensures the total liquid flowrate is equal to the sum of the liquid flowrates of each component on each tray on each section. - - Parameters - ---------- - m : Pyomo ConcreteModel - The optimization model. - sec : int - The section index. - n_tray : int - The tray index. - - Returns - ------- - Constraint - The constraint that enforces the total liquid flowrate is equal to the sum of the liquid flowrates of each component on each tray on each section. - """ - return sum(m.L[sec, n_tray, comp] for comp in m.comp) == m.Ltotal[sec, n_tray] - - @m.Constraint(m.comp, doc="Bottoms and liquid relation") - def bottoms_equality(m, comp): - """ - Constraint that ensures the bottoms flowrate is equal to the liquid flowrate of each component on the reboiler tray. - - Parameters - ---------- - m : Pyomo ConcreteModel - The optimization model. - comp : int - The component index. - - Returns - ------- - Constraint - The constraint that enforces the bottoms flowrate is equal to the liquid flowrate of each component on the reboiler tray. - """ - return m.B[comp] == m.L[1, m.reb_tray, comp] - - @m.Constraint(m.comp) - def condenser_total(m, comp): - """ - Constraint that ensures the distillate flowrate in the condenser is null. - - Parameters - ---------- - m : Pyomo ConcreteModel - The optimization model. - comp : int - The component index. - Returns - ------- - Constraint - The constraint that enforces the distillate flowrate is equal to zero in the condenser. - """ - return m.V[4, m.con_tray, comp] == 0 - - @m.Constraint() - def total_bottoms_product(m): - """ - Constraint that ensures the total bottoms flowrate is equal to the sum of the bottoms flowrates of each component. - - Parameters - ---------- - m : Pyomo ConcreteModel - The optimization model. - - Returns - ------- - Constraint - The constraint that enforces the total bottoms flowrate is equal to the sum of the bottoms flowrates of each component. - """ - return sum(m.B[comp] for comp in m.comp) == m.Btotal - - @m.Constraint() - def total_distillate_product(m): - """ - Constraint that ensures the total distillate flowrate is equal to the sum of the distillate flowrates of each component. - - Parameters - ---------- - m : Pyomo ConcreteModel - The optimization model. - - Returns - ------- - Constraint - The constraint that enforces the total distillate flowrate is equal to the sum of the distillate flowrates of each component. - """ - return sum(m.D[comp] for comp in m.comp) == m.Dtotal - - @m.Constraint(m.so) - def total_side_product(m, so): - """ - Constraint that ensures the total side product flowrate is equal to the sum of the side product flowrates of each component. - - Parameters - ---------- - m : Pyomo ConcreteModel - The optimization model. - so : int - The side product index, 2 or 3 for the intermediate side products. - - Returns - ------- - Constraint - The constraint that enforces the total side product flowrate is equal to the sum of the side product flowrates of each component. - """ - return sum(m.S[so, comp] for comp in m.comp) == m.Stotal[so] - - # Considers the number of existent trays and operating costs (condenser and reboiler heat duties) in the column. To ensure equal weights to the capital and operating costs, the number of existent trays is multiplied by a weight coefficient of 1000. - - m.obj = Objective( - expr=(m.Qcon + m.Qreb) * m.Hscale - + 1e3 - * ( - sum( - sum( - m.tray_exists[sec, n_tray].binary_indicator_var - for n_tray in m.candidate_trays_main - ) - for sec in m.section_main - ) - + sum( - m.tray_exists[2, n_tray].binary_indicator_var - for n_tray in m.candidate_trays_feed - ) - + sum( - m.tray_exists[3, n_tray].binary_indicator_var - for n_tray in m.candidate_trays_product - ) - + 1 - ), - sense=minimize, - doc="Objective function to minimize the operating costs and number of existent trays in the column", - ) - - @m.Constraint( - m.section_main, m.candidate_trays_main, doc="Logic proposition for main section" - ) - def _logic_proposition_main(m, sec, n_tray): - """ - Apply a logic proposition constraint to the main section and candidate trays to specify the order of trays in the column is from bottom to top provided the condition is met. - - Parameters - ---------- - m : Pyomo ConcreteModel - The optimization model. - sec : int - The section index. - n_tray : int - The tray index. - - Returns - ------- - Constraint or NoConstraint - The constraint expression or NoConstraint if the condition is not met. - """ - - if n_tray > m.reb_tray and (n_tray + 1) < m.num_trays: - return ( - m.tray_exists[sec, n_tray].binary_indicator_var - <= m.tray_exists[sec, n_tray + 1].binary_indicator_var - ) - else: - return Constraint.NoConstraint - - @m.Constraint(m.candidate_trays_feed) - def _logic_proposition_feed(m, n_tray): - """ - Apply a logic proposition constraint to the feed section and candidate trays to specify the order of trays in the column is from bottom to top provided the condition is met. - - Parameters - ---------- - m : Pyomo ConcreteModel - The optimization model. - n_tray : int - The tray index. - - Returns - ------- - Constraint or NoConstraint - The constraint expression or NoConstraint if the condition is not met. - """ - if n_tray > m.reb_tray and (n_tray + 1) < m.feed_tray: - return ( - m.tray_exists[2, n_tray].binary_indicator_var - <= m.tray_exists[2, n_tray + 1].binary_indicator_var - ) - elif n_tray > m.feed_tray and (n_tray + 1) < m.con_tray: - return ( - m.tray_exists[2, n_tray + 1].binary_indicator_var - <= m.tray_exists[2, n_tray].binary_indicator_var - ) - else: - return Constraint.NoConstraint - - @m.Constraint(m.candidate_trays_product) - def _logic_proposition_section3(m, n_tray): - """ - Apply a logic proposition constraint to the product section and candidate trays to specify the order of trays in the column is from bottom to top provided the condition is met. - - Parameters - ---------- - m : Pyomo ConcreteModel - The optimization model. - n_tray : int - The tray index. - - Returns - ------- - Constraint or NoConstraint - The constraint expression or NoConstraint if the condition is not met. - """ - if n_tray > 1 and (n_tray + 1) < m.num_trays: - return ( - m.tray_exists[3, n_tray].binary_indicator_var - <= m.tray_exists[3, n_tray + 1].binary_indicator_var - ) - else: - return Constraint.NoConstraint - - @m.Constraint(m.tray) - def equality_feed_product_side(m, n_tray): - """ - Constraint that enforces the equality of the binary indicator variables for the feed and product side trays. - - Parameters - ---------- - m : Pyomo ConcreteModel - The optimization model. - n_tray : int - The tray index. - - Returns - ------- - Constraint - The constraint expression that enforces the equality of the binary indicator variables for the feed and product side trays. - """ - return ( - m.tray_exists[2, n_tray].binary_indicator_var - == m.tray_exists[3, n_tray].binary_indicator_var - ) - - @m.Constraint() - def _existent_minimum_numbertrays(m): - """ - Constraint that enforces the minimum number of trays in the column. - - Parameters - ---------- - m : Pyomo ConcreteModel - The optimization model. - - Returns - ------- - Constraint - The constraint expression that enforces the minimum number of trays in each section and each tray to be greater than or equal to the minimum number of trays. - """ - return sum( - sum(m.tray_exists[sec, n_tray].binary_indicator_var for n_tray in m.tray) - for sec in m.section - ) - sum( - m.tray_exists[3, n_tray].binary_indicator_var for n_tray in m.tray - ) >= int( - m.min_tray - ) - - return m - - -def enforce_tray_exists(m, sec, n_tray): - """ - Enforce the existence of a tray in the column. - - Parameters - ---------- - m : Pyomo ConcreteModel - The optimization model. - sec : int - The section index. - n_tray : int - The tray index. - """ - m.tray_exists[sec, n_tray].indicator_var.fix(True) - m.tray_absent[sec, n_tray].deactivate() - - -def _build_tray_equations(m, sec, n_tray): - """ - Build the equations for the tray in the column as a function of the section when the tray exists. - Points to the appropriate function to build the equations for the section in the column. - - Parameters - ---------- - m : Pyomo ConcreteModel - The optimization model. - sec : int - The section index. - n_tray : int - The tray index. - - Returns - ------- - None - """ - build_function = { - 1: _build_bottom_equations, - 2: _build_feed_side_equations, - 3: _build_product_side_equations, - 4: _build_top_equations, - } - build_function[sec](m, n_tray) - - -def _build_bottom_equations(disj, n_tray): - """ - Build the equations for the bottom section in the column. - - Parameters - ---------- - disj : Disjunct - The disjunct object for the bottom section in the column. - n_tray : int - The tray index. - - Returns - ------- - None - """ - m = disj.model() - - @disj.Constraint(m.comp, doc="Bottom section 1 mass per component balances") - def _bottom_mass_percomponent_balances(disj, comp): - """ - Mass per component balances for the bottom section in the column. - - Parameters - ---------- - disj : Disjunct - The disjunct object for the bottom section in the column. - comp : int - The component index. - - Returns - ------- - Constraint - The constraint expression that enforces the mass balance per component in the bottom section of the column. - """ - return (m.L[1, n_tray + 1, comp] if n_tray < m.num_trays else 0) + ( - m.L[2, 1, comp] if n_tray == m.num_trays else 0 - ) + (m.L[3, 1, comp] if n_tray == m.num_trays else 0) - ( - m.L[1, n_tray, comp] if n_tray > m.reb_tray else 0 - ) + ( - m.V[1, n_tray - 1, comp] if n_tray > m.reb_tray else 0 - ) - ( - m.V[1, n_tray, comp] * m.dv[2] if n_tray == m.num_trays else 0 - ) - ( - m.V[1, n_tray, comp] * m.dv[3] if n_tray == m.num_trays else 0 - ) - ( - m.V[1, n_tray, comp] if n_tray < m.num_trays else 0 - ) - ( - m.B[comp] if n_tray == m.reb_tray else 0 - ) == m.slack[ - 1, n_tray, comp - ] - - @disj.Constraint(doc="Bottom section 1 energy balances") - def _bottom_energy_balances(disj): - """ - Energy balances for the bottom section in the column. - - Parameters - ---------- - disj : Disjunct - The disjunct object for the bottom section in the column. - - Returns - ------- - Constraint - The constraint expression that enforces the energy balance for the bottom section in the column. - """ - return ( - sum( - ( - m.L[1, n_tray + 1, comp] * m.hl[1, n_tray + 1, comp] - if n_tray < m.num_trays - else 0 - ) - + (m.L[2, 1, comp] * m.hl[2, 1, comp] if n_tray == m.num_trays else 0) - + (m.L[3, 1, comp] * m.hl[3, 1, comp] if n_tray == m.num_trays else 0) - - ( - m.L[1, n_tray, comp] * m.hl[1, n_tray, comp] - if n_tray > m.reb_tray - else 0 - ) - + ( - m.V[1, n_tray - 1, comp] * m.hv[1, n_tray - 1, comp] - if n_tray > m.reb_tray - else 0 - ) - - ( - m.V[1, n_tray, comp] * m.dv[2] * m.hv[1, n_tray, comp] - if n_tray == m.num_trays - else 0 - ) - - ( - m.V[1, n_tray, comp] * m.dv[3] * m.hv[1, n_tray, comp] - if n_tray == m.num_trays - else 0 - ) - - ( - m.V[1, n_tray, comp] * m.hv[1, n_tray, comp] - if n_tray < m.num_trays - else 0 - ) - - (m.B[comp] * m.hl[1, n_tray, comp] if n_tray == m.reb_tray else 0) - for comp in m.comp - ) - * m.Qscale - + (m.Qreb if n_tray == m.reb_tray else 0) - == 0 - ) - - @disj.Constraint(m.comp, doc="Bottom section 1 liquid flowrate per component") - def _bottom_liquid_percomponent(disj, comp): - """ - Liquid flowrate per component in the bottom section of the column. - - Parameters - ---------- - disj : Disjunct - The disjunct object for the bottom section in the column. - comp : int - The component index. - - Returns - ------- - Constraint - The constraint expression that enforces the liquid flowrate per component in the bottom section of the column. - """ - return m.L[1, n_tray, comp] == m.Ltotal[1, n_tray] * m.x[1, n_tray, comp] - - @disj.Constraint(m.comp, doc="Bottom section 1 vapor flowrate per component") - def _bottom_vapor_percomponent(disj, comp): - """ - Vapor flowrate per component in the bottom section of the column. - - Parameters - ---------- - disj : Disjunct - The disjunct object for the bottom section in the column. - comp : int - The component index. - - Returns - ------- - Constraint - The constraint expression that enforces the vapor flowrate per component in the bottom section of the column. - """ - return m.V[1, n_tray, comp] == m.Vtotal[1, n_tray] * m.y[1, n_tray, comp] - - @disj.Constraint(doc="Bottom section 1 liquid composition equilibrium summation") - def bottom_liquid_composition_summation(disj): - """ - Liquid composition equilibrium summation for the bottom section in the column. - - Parameters - ---------- - disj : Disjunct - The disjunct object for the bottom section in the column. - - Returns - ------- - Constraint - The constraint expression that enforces the liquid composition equilibrium summation for the bottom section in the column. - It ensures the sum of the liquid compositions is equal to 1 plus the error in the liquid composition. - """ - return sum(m.x[1, n_tray, comp] for comp in m.comp) - 1 == m.errx[1, n_tray] - - @disj.Constraint(doc="Bottom section 1 vapor composition equilibrium summation") - def bottom_vapor_composition_summation(disj): - """ - Vapor composition equilibrium summation for the bottom section in the column. - - Parameters - ---------- - disj : Disjunct - The disjunct object for the bottom section in the column. - - Returns - ------- - Constraint - The constraint expression that enforces the vapor composition equilibrium summation for the bottom section in the column. - It ensures the sum of the vapor compositions is equal to 1 plus the error in the vapor composition. - """ - return sum(m.y[1, n_tray, comp] for comp in m.comp) - 1 == m.erry[1, n_tray] - - @disj.Constraint( - m.comp, - doc="Temperature term for the bottom section 1 used in the vapor composition equation", - ) - def _bottom_reduced_temperature(disj, comp): - """Calculate the temperature term for the bottom section vapor composition. - - The constraint is used to avoid division by zero in the MINLP Hull reformulation. - - Parameters - ---------- - disj : Disjunct - The disjunct object for the bottom section in the column. - comp : int - The component index. - - Returns - ------- - Constraint - The constraint equation to calculate the temperature term for the bottom section vapor composition. - """ - return m.Tr[1, n_tray, comp] * m.T[1, n_tray] == m.prop[comp, 'TC'] - + ) + + m.x = Var( + m.section, + m.tray, + m.comp, + doc="Liquid composition [mol/mol]", + domain=NonNegativeReals, + bounds=(0, 1), + initialize=m.x0, + ) + m.y = Var( + m.section, + m.tray, + m.comp, + doc="Vapor composition [mol/mol]", + domain=NonNegativeReals, + bounds=(0, 1), + initialize=m.y0, + ) + + m.dl = Var( + m.dw, + doc="Liquid distributor in the dividing wall sections", + bounds=(0.2, 0.8), + initialize=m.dl0, + ) + m.dv = Var( + m.dw, + doc="Vapor distributor in the dividing wall sections", + bounds=(0, 1), + domain=NonNegativeReals, + initialize=m.dv0, + ) + + m.V = Var( + m.section, + m.tray, + m.comp, + doc="Vapor flowrate [mol/s]", + domain=NonNegativeReals, + bounds=(0, m.flow_max), + initialize=m.V0, + ) + m.L = Var( + m.section, + m.tray, + m.comp, + doc="Liquid flowrate [mol/s]", + domain=NonNegativeReals, + bounds=(0, m.flow_max), + initialize=m.L0, + ) + m.Vtotal = Var( + m.section, + m.tray, + doc="Total vapor flowrate [mol/s]", + domain=NonNegativeReals, + bounds=(0, m.flow_max), + initialize=m.Vtotal0, + ) + m.Ltotal = Var( + m.section, + m.tray, + doc="Total liquid flowrate [mol/s]", + domain=NonNegativeReals, + bounds=(0, m.flow_max), + initialize=m.Ltotal0, + ) + + m.D = Var( + m.comp, + doc="Distillate flowrate [mol/s]", + domain=NonNegativeReals, + bounds=(0, m.flow_max), + initialize=m.Ddes, + ) + m.B = Var( + m.comp, + doc="Bottoms flowrate [mol/s]", + domain=NonNegativeReals, + bounds=(0, m.flow_max), + initialize=m.Bdes, + ) + m.S = Var( + m.so, + m.comp, + doc="Product 2 and 3 flowrates [mol/s]", + domain=NonNegativeReals, + bounds=(0, m.flow_max), + initialize=m.Sdes, + ) + m.Dtotal = Var( + doc="Distillate flowrate [mol/s]", + domain=NonNegativeReals, + bounds=(0, m.flow_max), + initialize=m.Ddes, + ) + m.Btotal = Var( + doc="Bottoms flowrate [mol/s]", + domain=NonNegativeReals, + bounds=(0, m.flow_max), + initialize=m.Bdes, + ) + m.Stotal = Var( + m.so, + doc="Total product 2 and 3 side flowrate [mol/s]", + domain=NonNegativeReals, + bounds=(0, m.flow_max), + initialize=m.Sdes, + ) + + m.hl = Var( + m.section, + m.tray, + m.comp, + doc='Liquid enthalpy [J/mol]', + bounds=(m.hllo, m.hlup), + initialize=m.hl0, + ) + m.hv = Var( + m.section, + m.tray, + m.comp, + doc='Vapor enthalpy [J/mol]', + bounds=(m.hvlo, m.hvup), + initialize=m.hv0, + ) + m.Qreb = Var( + doc="Reboiler heat duty [J/s]", + domain=NonNegativeReals, + bounds=(0, m.Qmax), + initialize=1, + ) + m.Qcon = Var( + doc="Condenser heat duty [J/s]", + domain=NonNegativeReals, + bounds=(0, m.Qmax), + initialize=1, + ) + + m.rr = Var( + doc="Internal reflux ratio in the column", + domain=NonNegativeReals, + bounds=(0.7, 1), + initialize=m.rr0, + ) + m.bu = Var( + doc="Boilup rate in the reboiler", + domain=NonNegativeReals, + bounds=(0.7, 1), + initialize=m.bu0, + ) + + m.F = Var( + m.comp, + doc="Side feed flowrate [mol/s]", + domain=NonNegativeReals, + bounds=(0, 50), + initialize=m.F0, + ) + m.q = Var( + doc="Vapor fraction in side feed", + domain=NonNegativeReals, + bounds=(0, 1), + initialize=1, + ) + + m.actv = Var( + m.section, + m.tray, + m.comp, + doc="Liquid activity coefficient", + domain=NonNegativeReals, + bounds=(0, 10), + initialize=m.actv0, + ) + + m.errx = Var( + m.section, + m.tray, + doc="Error in liquid composition [mol/mol]", + bounds=(-1e-3, 1e-3), + initialize=0, + ) + m.erry = Var( + m.section, + m.tray, + doc="Error in vapor composition [mol/mol]", + bounds=(-1e-3, 1e-3), + initialize=0, + ) + m.slack = Var( + m.section, + m.tray, + m.comp, + doc="Slack variable", + bounds=(-1e-8, 1e-8), + initialize=0, + ) + + m.tray_exists = Disjunct( + m.section, + m.tray, + doc="Disjunct that enforce the existence of each tray", + rule=_build_tray_equations, + ) + m.tray_absent = Disjunct( + m.section, + m.tray, + doc="Disjunct that enforce the absence of each tray", + rule=_build_pass_through_eqns, + ) + + @m.Disjunction( + m.section, m.tray, doc="Disjunction between whether each tray exists or not" + ) + def tray_exists_or_not(m, sec, n_tray): + """ + Disjunction between whether each tray exists or not. + + Parameters + ---------- + m : Pyomo ConcreteModel + The Pyomo model object. + sec : int + The section index. + n_tray : int + The tray index. + + Returns + ------- + Disjunction + The disjunction between whether each tray exists or not. + """ + return [m.tray_exists[sec, n_tray], m.tray_absent[sec, n_tray]] + + @m.Constraint(m.section_main) + def minimum_trays_main(m, sec): + """ + Constraint that ensures the minimum number of trays in the main section. + + Parameters + ---------- + m : Pyomo ConcreteModel + The model object for the GDP Kaibel Column. + sec : Set + The section index. + + Returns + ------- + Constraint + A constraint expression that enforces the minimum number of trays in the main section to be greater than or equal to the minimum number of trays. + """ + return ( + sum( + m.tray_exists[sec, n_tray].binary_indicator_var + for n_tray in m.candidate_trays_main + ) + + 1 + >= m.min_num_trays + ) + + @m.Constraint() + def minimum_trays_feed(m): + """ + Constraint function that ensures the minimum number of trays in the feed section is met. + + Parameters + ---------- + m : Pyomo ConcreteModel + The Pyomo model object. + + Returns + ------- + Constraint + The constraint expression that enforces the minimum number of trays is greater than or equal to the minimum number of trays. + """ + return ( + sum( + m.tray_exists[2, n_tray].binary_indicator_var + for n_tray in m.candidate_trays_feed + ) + + 1 + >= m.min_num_trays + ) + + # TOCHECK: pyomo.GDP Syntax + + @m.Constraint() + def minimum_trays_product(m): + """ + Constraint function to calculate the minimum number of trays in the product section. + + Parameters + ---------- + m : Pyomo ConcreteModel + The optimization model. + + Returns + ------- + Constraint + The constraint expression that enforces the minimum number of trays is greater than or equal to the minimum number of trays. + """ + return ( + sum( + m.tray_exists[3, n_tray].binary_indicator_var + for n_tray in m.candidate_trays_product + ) + + 1 + >= m.min_num_trays + ) + + ## Fixed trays + enforce_tray_exists(m, 1, 1) # reboiler + enforce_tray_exists(m, 1, m.num_trays) # vapor distributor + enforce_tray_exists(m, 2, 1) # dividing wall starting tray + enforce_tray_exists(m, 2, m.feed_tray) # feed tray + enforce_tray_exists(m, 2, m.num_trays) # dividing wall ending tray + enforce_tray_exists(m, 3, 1) # dividing wall starting tray + enforce_tray_exists(m, 3, m.sideout1_tray) # side outlet 1 for product 3 + enforce_tray_exists(m, 3, m.sideout2_tray) # side outlet 2 for product 2 + enforce_tray_exists(m, 3, m.num_trays) # dividing wall ending tray + enforce_tray_exists(m, 4, 1) # liquid distributor + enforce_tray_exists(m, 4, m.num_trays) # condenser + + #### Global constraints + @m.Constraint( + m.dw, m.tray, doc="Monotonic temperature in the dividing wall sections" + ) + def monotonic_temperature(m, sec, n_tray): + """This function returns a constraint object representing the monotonic temperature constraint. + + The monotonic temperature constraint ensures that the temperature on each tray in the distillation column + is less than or equal to the temperature on the top tray. + + Parameters + ---------- + m : Pyomo ConcreteModel + The Pyomo model object. + sec : Set + The set of sections in the dividing wall sections. + n_tray : Set + The set of trays in the distillation column. + + Returns + ------- + Constraint + The monotonic temperature constraint specifying that the temperature on each tray is less than or equal to the temperature on the top tray of section 1 which is the condenser. + """ + return m.T[sec, n_tray] <= m.T[1, m.num_trays] + + @m.Constraint(doc="Liquid distributor") + def liquid_distributor(m): + """Defines the liquid distributor constraint. + + This constraint ensures that the sum of the liquid distributors in all sections is equal to 1. + + Parameters + ---------- + m : Pyomo ConcreteModel + The optimization model. + + Returns + ------- + Constraint + The liquid distributor constraint that enforces the sum of the liquid flow rates in all sections is equal to 1. + + """ + return sum(m.dl[sec] for sec in m.dw) - 1 == 0 + + @m.Constraint(doc="Vapor distributor") + def vapor_distributor(m): + """ + Add a constraint to ensure that the sum of the vapor distributors is equal to 1. + + Parameters + ---------- + m : Pyomo ConcreteModel + The Pyomo model object. + + Returns + ------- + Constraint + The vapor distributor constraint. + """ + return sum(m.dv[sec] for sec in m.dw) - 1 == 0 + + @m.Constraint(doc="Reboiler composition specification") + def heavy_product(m): + """ + Reboiler composition specification for the heavy component in the feed. + + Parameters + ---------- + m : Pyomo ConcreteModel + The optimization model. + + Returns + ------- + Constraint + The constraint that enforces the reboiler composition is greater than or equal to the specified composition xspechc final liquid composition for butanol, the heavy component in the feed. + """ + return m.x[1, m.reb_tray, m.hc] >= m.xspec_hc + + @m.Constraint(doc="Condenser composition specification") + def light_product(m): + """ + Condenser composition specification for the light component in the feed. + + Parameters + ---------- + m : Model + The optimization model. + + Returns + ------- + Constraint + The constraint that enforces the condenser composition is greater than or equal to the specified final liquid composition for ethanol, xspeclc , the light component in the feed. + """ + return m.x[4, m.con_tray, m.lc] >= m.xspec_lc + + @m.Constraint(doc="Side outlet 1 final liquid composition") + def intermediate1_product(m): + ''' + This constraint ensures that the intermediate 1 final liquid composition is greater than or equal to the specified composition xspec_inter1, which is the final liquid composition for ethanol. + + Parameters + ---------- + m : Pyomo ConcreteModel. + The optimization model. + + Returns + ------- + Constraint + The constraint that enforces the intermediate 1 final liquid composition is greater than or equal to the specified composition xspec_inter1, which is the final liquid composition for ethanol. + ''' + return m.x[3, m.sideout1_tray, 3] >= m.xspec_inter3 + + @m.Constraint(doc="Side outlet 2 final liquid composition") + def intermediate2_product(m): + """ + This constraint ensures that the intermediate 2 final liquid composition is greater than or equal to the specified composition xspec_inter2, which is the final liquid composition for butanol. + + Parameters + ---------- + m : Pyomo ConcreteModel + The optimization model. + + Returns + ------- + Constraint + The constraint that enforces the intermediate 2 final liquid composition is greater than or equal to the specified composition xspec_inter2, which is the final liquid composition for butanol. + """ + return m.x[3, m.sideout2_tray, 2] >= m.xspec_inter2 + + @m.Constraint(doc="Reboiler flowrate") + def _heavy_product_flow(m): + """ + Reboiler flowrate constraint that ensures the reboiler flowrate is greater than or equal to the specified flowrate Bdes, which is the flowrate of butanol. + + Parameters + ---------- + m : Pyomo ConcreteModel + The optimization model. + + Returns + ------- + Constraint + The constraint that enforces the reboiler flowrate is greater than or equal to the specified flowrate Bdes, which is the flowrate of butanol. + """ + return m.Btotal >= m.Bdes + + @m.Constraint(doc="Condenser flowrate") + def _light_product_flow(m): + """ + Condenser flowrate constraint that ensures the condenser flowrate is greater than or equal to the specified flowrate Ddes, which is the flowrate of ethanol. + + Parameters + ---------- + m : Pyomo ConcreteModel + The optimization model. + + Returns + ------- + Constraint + The constraint that enforces the condenser flowrate is greater than or equal to the specified flowrate Ddes, which is the flowrate of ethanol. + """ + return m.Dtotal >= m.Ddes + + @m.Constraint(m.so, doc="Intermediate flowrate") + def _intermediate_product_flow(m, so): + """ + Intermediate flowrate constraint that ensures the intermediate flowrate is greater than or equal to the specified flowrate Sdes, which is the flowrate of the intermediate side product 2 and 3. + + Parameters + ---------- + m : Pyomo ConcreteModel + The optimization model. + so : int + The side product outlet index. + + Returns + ------- + Constraint + The constraint that enforces the intermediate flowrate is greater than or equal to the specified flowrate Sdes, which is the flowrate of the intermediate side product 2 and 3. + """ + return m.Stotal[so] >= m.Sdes + + @m.Constraint(doc="Internal boilup ratio, V/L") + def _internal_boilup_ratio(m): + """ + Internal boilup ratio constraint that ensures the internal boilup ratio is equal to the boilup rate times the liquid flowrate on the reboiler tray is equal to the vapor flowrate on the tray above the reboiler tray. + + Parameters + ---------- + m : Pyomo ConcreteModel + The optimization model. + + Returns + ------- + Constraint + The constraint that enforces the boilup rate times the liquid flowrate on the reboiler tray is equal to the vapor flowrate on the tray above the reboiler tray. + """ + return m.bu * m.Ltotal[1, m.reb_tray + 1] == m.Vtotal[1, m.reb_tray] + + @m.Constraint(doc="Internal reflux ratio, L/V") + def internal_reflux_ratio(m): + """ + Internal reflux ratio constraint that ensures the internal reflux ratio is equal to the reflux rate times the vapor flowrate on the tray above the condenser tray is equal to the liquid flowrate on the condenser tray. + + Parameters + ---------- + m : Pyomo ConcreteModel + The optimization model. + + Returns + ------- + Constraint + The constraint that enforces the reflux rate times the vapor flowrate on the tray above the condenser tray is equal to the liquid flowrate on the condenser tray. + """ + return m.rr * m.Vtotal[4, m.con_tray - 1] == m.Ltotal[4, m.con_tray] + + @m.Constraint(doc="External boilup ratio relation with bottoms") + def _external_boilup_ratio(m): + """ + External boilup ratio constraint that ensures the external boilup ratio times the liquid flowrate on the reboiler tray is equal to the bottoms flowrate. + + Parameters + ---------- + m : Pyomo ConcreteModel + The optimization model. + + Returns + ------- + Constraint + The constraint that enforces the external boilup ratio times the liquid flowrate on the reboiler tray is equal to the bottoms flowrate. + """ + return m.Btotal == (1 - m.bu) * m.Ltotal[1, m.reb_tray + 1] + + @m.Constraint(doc="External reflux ratio relation with distillate") + def _external_reflux_ratio(m): + """ + External reflux ratio constraint that ensures the external reflux ratio times the vapor flowrate on the tray above the condenser tray is equal to the distillate flowrate. + + Parameters + ---------- + m : Pyomo ConcreteModel + The optimization model. + + Returns + ------- + Constraint + The constraint that enforces the external reflux ratio times the vapor flowrate on the tray above the condenser tray is equal to the distillate flowrate. + """ + return m.Dtotal == (1 - m.rr) * m.Vtotal[4, m.con_tray - 1] + + @m.Constraint(m.section, m.tray, doc="Total vapor flowrate") + def _total_vapor_flowrate(m, sec, n_tray): + """ + Constraint that ensures the total vapor flowrate is equal to the sum of the vapor flowrates of each component on each tray. + + Parameters + ---------- + m : Pyomo ConcreteModel + The optimization model. + sec : int + The section index. + n_tray : int + The tray index. + + Returns + ------- + Constraint + The constraint that enforces the total vapor flowrate is equal to the sum of the vapor flowrates of each component on each tray on each section. + """ + return sum(m.V[sec, n_tray, comp] for comp in m.comp) == m.Vtotal[sec, n_tray] + + @m.Constraint(m.section, m.tray, doc="Total liquid flowrate") + def _total_liquid_flowrate(m, sec, n_tray): + """ + Constraint that ensures the total liquid flowrate is equal to the sum of the liquid flowrates of each component on each tray on each section. + + Parameters + ---------- + m : Pyomo ConcreteModel + The optimization model. + sec : int + The section index. + n_tray : int + The tray index. + + Returns + ------- + Constraint + The constraint that enforces the total liquid flowrate is equal to the sum of the liquid flowrates of each component on each tray on each section. + """ + return sum(m.L[sec, n_tray, comp] for comp in m.comp) == m.Ltotal[sec, n_tray] + + @m.Constraint(m.comp, doc="Bottoms and liquid relation") + def bottoms_equality(m, comp): + """ + Constraint that ensures the bottoms flowrate is equal to the liquid flowrate of each component on the reboiler tray. + + Parameters + ---------- + m : Pyomo ConcreteModel + The optimization model. + comp : int + The component index. + + Returns + ------- + Constraint + The constraint that enforces the bottoms flowrate is equal to the liquid flowrate of each component on the reboiler tray. + """ + return m.B[comp] == m.L[1, m.reb_tray, comp] + + @m.Constraint(m.comp) + def condenser_total(m, comp): + """ + Constraint that ensures the distillate flowrate in the condenser is null. + + Parameters + ---------- + m : Pyomo ConcreteModel + The optimization model. + comp : int + The component index. + Returns + ------- + Constraint + The constraint that enforces the distillate flowrate is equal to zero in the condenser. + """ + return m.V[4, m.con_tray, comp] == 0 + + @m.Constraint() + def total_bottoms_product(m): + """ + Constraint that ensures the total bottoms flowrate is equal to the sum of the bottoms flowrates of each component. + + Parameters + ---------- + m : Pyomo ConcreteModel + The optimization model. + + Returns + ------- + Constraint + The constraint that enforces the total bottoms flowrate is equal to the sum of the bottoms flowrates of each component. + """ + return sum(m.B[comp] for comp in m.comp) == m.Btotal + + @m.Constraint() + def total_distillate_product(m): + """ + Constraint that ensures the total distillate flowrate is equal to the sum of the distillate flowrates of each component. + + Parameters + ---------- + m : Pyomo ConcreteModel + The optimization model. + + Returns + ------- + Constraint + The constraint that enforces the total distillate flowrate is equal to the sum of the distillate flowrates of each component. + """ + return sum(m.D[comp] for comp in m.comp) == m.Dtotal + + @m.Constraint(m.so) + def total_side_product(m, so): + """ + Constraint that ensures the total side product flowrate is equal to the sum of the side product flowrates of each component. + + Parameters + ---------- + m : Pyomo ConcreteModel + The optimization model. + so : int + The side product index, 2 or 3 for the intermediate side products. + + Returns + ------- + Constraint + The constraint that enforces the total side product flowrate is equal to the sum of the side product flowrates of each component. + """ + return sum(m.S[so, comp] for comp in m.comp) == m.Stotal[so] + + # Considers the number of existent trays and operating costs (condenser and reboiler heat duties) in the column. To ensure equal weights to the capital and operating costs, the number of existent trays is multiplied by a weight coefficient of 1000. + + m.obj = Objective( + expr=(m.Qcon + m.Qreb) * m.Hscale + + 1e3 + * ( + sum( + sum( + m.tray_exists[sec, n_tray].binary_indicator_var + for n_tray in m.candidate_trays_main + ) + for sec in m.section_main + ) + + sum( + m.tray_exists[2, n_tray].binary_indicator_var + for n_tray in m.candidate_trays_feed + ) + + sum( + m.tray_exists[3, n_tray].binary_indicator_var + for n_tray in m.candidate_trays_product + ) + + 1 + ), + sense=minimize, + doc="Objective function to minimize the operating costs and number of existent trays in the column", + ) + + @m.Constraint( + m.section_main, m.candidate_trays_main, doc="Logic proposition for main section" + ) + def _logic_proposition_main(m, sec, n_tray): + """ + Apply a logic proposition constraint to the main section and candidate trays to specify the order of trays in the column is from bottom to top provided the condition is met. + + Parameters + ---------- + m : Pyomo ConcreteModel + The optimization model. + sec : int + The section index. + n_tray : int + The tray index. + + Returns + ------- + Constraint or NoConstraint + The constraint expression or NoConstraint if the condition is not met. + """ + + if n_tray > m.reb_tray and (n_tray + 1) < m.num_trays: + return ( + m.tray_exists[sec, n_tray].binary_indicator_var + <= m.tray_exists[sec, n_tray + 1].binary_indicator_var + ) + else: + return Constraint.NoConstraint + + @m.Constraint(m.candidate_trays_feed) + def _logic_proposition_feed(m, n_tray): + """ + Apply a logic proposition constraint to the feed section and candidate trays to specify the order of trays in the column is from bottom to top provided the condition is met. + + Parameters + ---------- + m : Pyomo ConcreteModel + The optimization model. + n_tray : int + The tray index. + + Returns + ------- + Constraint or NoConstraint + The constraint expression or NoConstraint if the condition is not met. + """ + if n_tray > m.reb_tray and (n_tray + 1) < m.feed_tray: + return ( + m.tray_exists[2, n_tray].binary_indicator_var + <= m.tray_exists[2, n_tray + 1].binary_indicator_var + ) + elif n_tray > m.feed_tray and (n_tray + 1) < m.con_tray: + return ( + m.tray_exists[2, n_tray + 1].binary_indicator_var + <= m.tray_exists[2, n_tray].binary_indicator_var + ) + else: + return Constraint.NoConstraint + + @m.Constraint(m.candidate_trays_product) + def _logic_proposition_section3(m, n_tray): + """ + Apply a logic proposition constraint to the product section and candidate trays to specify the order of trays in the column is from bottom to top provided the condition is met. + + Parameters + ---------- + m : Pyomo ConcreteModel + The optimization model. + n_tray : int + The tray index. + + Returns + ------- + Constraint or NoConstraint + The constraint expression or NoConstraint if the condition is not met. + """ + if n_tray > 1 and (n_tray + 1) < m.num_trays: + return ( + m.tray_exists[3, n_tray].binary_indicator_var + <= m.tray_exists[3, n_tray + 1].binary_indicator_var + ) + else: + return Constraint.NoConstraint + + @m.Constraint(m.tray) + def equality_feed_product_side(m, n_tray): + """ + Constraint that enforces the equality of the binary indicator variables for the feed and product side trays. + + Parameters + ---------- + m : Pyomo ConcreteModel + The optimization model. + n_tray : int + The tray index. + + Returns + ------- + Constraint + The constraint expression that enforces the equality of the binary indicator variables for the feed and product side trays. + """ + return ( + m.tray_exists[2, n_tray].binary_indicator_var + == m.tray_exists[3, n_tray].binary_indicator_var + ) + + @m.Constraint() + def _existent_minimum_numbertrays(m): + """ + Constraint that enforces the minimum number of trays in the column. + + Parameters + ---------- + m : Pyomo ConcreteModel + The optimization model. + + Returns + ------- + Constraint + The constraint expression that enforces the minimum number of trays in each section and each tray to be greater than or equal to the minimum number of trays. + """ + return sum( + sum(m.tray_exists[sec, n_tray].binary_indicator_var for n_tray in m.tray) + for sec in m.section + ) - sum( + m.tray_exists[3, n_tray].binary_indicator_var for n_tray in m.tray + ) >= int( + m.min_tray + ) + + return m + + +def enforce_tray_exists(m, sec, n_tray): + """ + Enforce the existence of a tray in the column. + + Parameters + ---------- + m : Pyomo ConcreteModel + The optimization model. + sec : int + The section index. + n_tray : int + The tray index. + """ + m.tray_exists[sec, n_tray].indicator_var.fix(True) + m.tray_absent[sec, n_tray].deactivate() + + +def _build_tray_equations(m, sec, n_tray): + """ + Build the equations for the tray in the column as a function of the section when the tray exists. + Points to the appropriate function to build the equations for the section in the column. + + Parameters + ---------- + m : Pyomo ConcreteModel + The optimization model. + sec : int + The section index. + n_tray : int + The tray index. + + Returns + ------- + None + """ + build_function = { + 1: _build_bottom_equations, + 2: _build_feed_side_equations, + 3: _build_product_side_equations, + 4: _build_top_equations, + } + build_function[sec](m, n_tray) + + +def _build_bottom_equations(disj, n_tray): + """ + Build the equations for the bottom section in the column. + + Parameters + ---------- + disj : Disjunct + The disjunct object for the bottom section in the column. + n_tray : int + The tray index. + + Returns + ------- + None + """ + m = disj.model() + + @disj.Constraint(m.comp, doc="Bottom section 1 mass per component balances") + def _bottom_mass_percomponent_balances(disj, comp): + """ + Mass per component balances for the bottom section in the column. + + Parameters + ---------- + disj : Disjunct + The disjunct object for the bottom section in the column. + comp : int + The component index. + + Returns + ------- + Constraint + The constraint expression that enforces the mass balance per component in the bottom section of the column. + """ + return (m.L[1, n_tray + 1, comp] if n_tray < m.num_trays else 0) + ( + m.L[2, 1, comp] if n_tray == m.num_trays else 0 + ) + (m.L[3, 1, comp] if n_tray == m.num_trays else 0) - ( + m.L[1, n_tray, comp] if n_tray > m.reb_tray else 0 + ) + ( + m.V[1, n_tray - 1, comp] if n_tray > m.reb_tray else 0 + ) - ( + m.V[1, n_tray, comp] * m.dv[2] if n_tray == m.num_trays else 0 + ) - ( + m.V[1, n_tray, comp] * m.dv[3] if n_tray == m.num_trays else 0 + ) - ( + m.V[1, n_tray, comp] if n_tray < m.num_trays else 0 + ) - ( + m.B[comp] if n_tray == m.reb_tray else 0 + ) == m.slack[ + 1, n_tray, comp + ] + + @disj.Constraint(doc="Bottom section 1 energy balances") + def _bottom_energy_balances(disj): + """ + Energy balances for the bottom section in the column. + + Parameters + ---------- + disj : Disjunct + The disjunct object for the bottom section in the column. + + Returns + ------- + Constraint + The constraint expression that enforces the energy balance for the bottom section in the column. + """ + return ( + sum( + ( + m.L[1, n_tray + 1, comp] * m.hl[1, n_tray + 1, comp] + if n_tray < m.num_trays + else 0 + ) + + (m.L[2, 1, comp] * m.hl[2, 1, comp] if n_tray == m.num_trays else 0) + + (m.L[3, 1, comp] * m.hl[3, 1, comp] if n_tray == m.num_trays else 0) + - ( + m.L[1, n_tray, comp] * m.hl[1, n_tray, comp] + if n_tray > m.reb_tray + else 0 + ) + + ( + m.V[1, n_tray - 1, comp] * m.hv[1, n_tray - 1, comp] + if n_tray > m.reb_tray + else 0 + ) + - ( + m.V[1, n_tray, comp] * m.dv[2] * m.hv[1, n_tray, comp] + if n_tray == m.num_trays + else 0 + ) + - ( + m.V[1, n_tray, comp] * m.dv[3] * m.hv[1, n_tray, comp] + if n_tray == m.num_trays + else 0 + ) + - ( + m.V[1, n_tray, comp] * m.hv[1, n_tray, comp] + if n_tray < m.num_trays + else 0 + ) + - (m.B[comp] * m.hl[1, n_tray, comp] if n_tray == m.reb_tray else 0) + for comp in m.comp + ) + * m.Qscale + + (m.Qreb if n_tray == m.reb_tray else 0) + == 0 + ) + + @disj.Constraint(m.comp, doc="Bottom section 1 liquid flowrate per component") + def _bottom_liquid_percomponent(disj, comp): + """ + Liquid flowrate per component in the bottom section of the column. + + Parameters + ---------- + disj : Disjunct + The disjunct object for the bottom section in the column. + comp : int + The component index. + + Returns + ------- + Constraint + The constraint expression that enforces the liquid flowrate per component in the bottom section of the column. + """ + return m.L[1, n_tray, comp] == m.Ltotal[1, n_tray] * m.x[1, n_tray, comp] + + @disj.Constraint(m.comp, doc="Bottom section 1 vapor flowrate per component") + def _bottom_vapor_percomponent(disj, comp): + """ + Vapor flowrate per component in the bottom section of the column. + + Parameters + ---------- + disj : Disjunct + The disjunct object for the bottom section in the column. + comp : int + The component index. + + Returns + ------- + Constraint + The constraint expression that enforces the vapor flowrate per component in the bottom section of the column. + """ + return m.V[1, n_tray, comp] == m.Vtotal[1, n_tray] * m.y[1, n_tray, comp] + + @disj.Constraint(doc="Bottom section 1 liquid composition equilibrium summation") + def bottom_liquid_composition_summation(disj): + """ + Liquid composition equilibrium summation for the bottom section in the column. + + Parameters + ---------- + disj : Disjunct + The disjunct object for the bottom section in the column. + + Returns + ------- + Constraint + The constraint expression that enforces the liquid composition equilibrium summation for the bottom section in the column. + It ensures the sum of the liquid compositions is equal to 1 plus the error in the liquid composition. + """ + return sum(m.x[1, n_tray, comp] for comp in m.comp) - 1 == m.errx[1, n_tray] + + @disj.Constraint(doc="Bottom section 1 vapor composition equilibrium summation") + def bottom_vapor_composition_summation(disj): + """ + Vapor composition equilibrium summation for the bottom section in the column. + + Parameters + ---------- + disj : Disjunct + The disjunct object for the bottom section in the column. + + Returns + ------- + Constraint + The constraint expression that enforces the vapor composition equilibrium summation for the bottom section in the column. + It ensures the sum of the vapor compositions is equal to 1 plus the error in the vapor composition. + """ + return sum(m.y[1, n_tray, comp] for comp in m.comp) - 1 == m.erry[1, n_tray] + + @disj.Constraint( + m.comp, + doc="Temperature term for the bottom section 1 used in the vapor composition equation", + ) + def _bottom_reduced_temperature(disj, comp): + """Calculate the temperature term for the bottom section vapor composition. + + The constraint is used to avoid division by zero in the MINLP Hull reformulation. + + Parameters + ---------- + disj : Disjunct + The disjunct object for the bottom section in the column. + comp : int + The component index. + + Returns + ------- + Constraint + The constraint equation to calculate the temperature term for the bottom section vapor composition. + """ + return m.Tr[1, n_tray, comp] * m.T[1, n_tray] == m.prop[comp, 'TC'] + @disj.Constraint(m.comp, doc="Bottom section 1 vapor composition") def bottom_vapor_composition(disj, comp): - """ - Vapor composition for the bottom section in the column. - - Parameters - ---------- - disj : Disjunct - The disjunct object for the bottom section in the column. - comp : int - The component index. - - Returns - ------- - Constraint - The constraint expression that enforces the vapor composition for the bottom section in the column. - The equation is derived from the vapor-liquid equilibrium relationship. - """ + """ + Vapor composition for the bottom section in the column. + + Parameters + ---------- + disj : Disjunct + The disjunct object for the bottom section in the column. + comp : int + The component index. + + Returns + ------- + Constraint + The constraint expression that enforces the vapor composition for the bottom section in the column. + The equation is derived from the vapor-liquid equilibrium relationship. + """ return m.y[1, n_tray, comp] * m.P[1, n_tray] == m.x[1, n_tray, comp] * ( m.actv[1, n_tray, comp] * ( @@ -1495,310 +1495,310 @@ def bottom_vapor_composition(disj, comp): ) ) ) - - @disj.Constraint(m.comp, doc="Bottom section 1 liquid enthalpy") - def bottom_liquid_enthalpy(disj, comp): - """ - Liquid enthalpy for the bottom section in the column. - - Parameters - ---------- - disj : Disjunct - The disjunct object for the bottom section in the column. - comp : int - The component index. - - Returns - ------- - Constraint - The constraint expression that enforces the liquid enthalpy for the bottom section in the column. - """ - return m.hl[1, n_tray, comp] == ( - m.cpc[1] - * ( - (m.T[1, n_tray] - m.Tref) * m.prop[comp, 'cpA', 1] - + (m.T[1, n_tray] ** 2 - m.Tref**2) - * m.prop[comp, 'cpB', 1] - * m.cpc2['A', 1] - / 2 - + (m.T[1, n_tray] ** 3 - m.Tref**3) - * m.prop[comp, 'cpC', 1] - * m.cpc2['B', 1] - / 3 - + (m.T[1, n_tray] ** 4 - m.Tref**4) * m.prop[comp, 'cpD', 1] / 4 - ) - / m.Hscale - ) - - @disj.Constraint(m.comp, doc="Bottom section 1 vapor enthalpy") - def bottom_vapor_enthalpy(disj, comp): - """ - Vapor enthalpy for the bottom section in the column. - - Parameters - ---------- - disj : Disjunct - The disjunct object for the bottom section in the column. - comp : int - The component index. - - Returns - ------- - Constraint - The constraint expression that enforces the vapor enthalpy for the bottom section in the column. - """ - return m.hv[1, n_tray, comp] == ( - m.cpc[2] - * ( - (m.T[1, n_tray] - m.Tref) * m.prop[comp, 'cpA', 2] - + (m.T[1, n_tray] ** 2 - m.Tref**2) - * m.prop[comp, 'cpB', 2] - * m.cpc2['A', 2] - / 2 - + (m.T[1, n_tray] ** 3 - m.Tref**3) - * m.prop[comp, 'cpC', 2] - * m.cpc2['B', 2] - / 3 - + (m.T[1, n_tray] ** 4 - m.Tref**4) * m.prop[comp, 'cpD', 2] / 4 - ) - / m.Hscale - + m.dHvap[comp] - ) - - @disj.Constraint(m.comp, doc="Bottom section 1 liquid activity coefficient") - def bottom_activity_coefficient(disj, comp): - """ - Liquid activity coefficient for the bottom section in the column equal to 1. - - Parameters - ---------- - disj : Disjunct - The disjunct object for the bottom section in the column. - comp : int - The component index. - - Returns - ------- - Constraint - The constraint expression that enforces the liquid activity coefficient for the bottom section for each component and tray in the column to be equal to 1. - """ - return m.actv[1, n_tray, comp] == 1 - - -def _build_feed_side_equations(disj, n_tray): - """ - Build the equations for the feed side section in the column. - - Parameters - ---------- - disj : Disjunct - The disjunct object for the feed side section in the column. - n_tray : int - The tray index. - - Returns - ------- - None - """ - m = disj.model() - - @disj.Constraint(m.comp, doc="Feed section 2 mass per component balances") - def _feedside_masspercomponent_balances(disj, comp): - """ - Mass per component balances for the feed side section in the column. - - Parameters - ---------- - disj : Disjunct - The disjunct object for the feed side section in the column. - comp : int - The component index. - - Returns - ------- - Constraint - The constraint expression that enforces the mass balance per component in the feed side section of the column. - """ - return (m.L[2, n_tray + 1, comp] if n_tray < m.num_trays else 0) + ( - m.L[4, 1, comp] * m.dl[2] if n_tray == m.num_trays else 0 - ) - m.L[2, n_tray, comp] + ( - m.V[1, m.num_trays, comp] * m.dv[2] if n_tray == 1 else 0 - ) + ( - m.V[2, n_tray - 1, comp] if n_tray > 1 else 0 - ) - m.V[ - 2, n_tray, comp - ] + ( - m.F[comp] if n_tray == m.feed_tray else 0 - ) == 0 - - @disj.Constraint(doc="Feed section 2 energy balances") - def _feedside_energy_balances(disj): - """ - Energy balances for the feed side section in the column. - - Parameters - ---------- - disj : Disjunct - The disjunct object for the feed side section in the column. - - Returns - ------- - Constraint - The constraint expression that enforces the energy balance for the feed side section in the column. - """ - return ( - sum( - ( - m.L[2, n_tray + 1, comp] * m.hl[2, n_tray + 1, comp] - if n_tray < m.num_trays - else 0 - ) - + ( - m.L[4, 1, comp] * m.dl[2] * m.hl[4, 1, comp] - if n_tray == m.num_trays - else 0 - ) - - m.L[2, n_tray, comp] * m.hl[2, n_tray, comp] - + ( - m.V[1, m.num_trays, comp] * m.dv[2] * m.hv[1, m.num_trays, comp] - if n_tray == 1 - else 0 - ) - + ( - m.V[2, n_tray - 1, comp] * m.hv[2, n_tray - 1, comp] - if n_tray > 1 - else 0 - ) - - m.V[2, n_tray, comp] * m.hv[2, n_tray, comp] - for comp in m.comp - ) - * m.Qscale - + sum( - ( - m.F[comp] * (m.hlf[comp] * (1 - m.q) + m.hvf[comp] * m.q) - if n_tray == m.feed_tray - else 0 - ) - for comp in m.comp - ) - * m.Qscale - == 0 - ) - - @disj.Constraint(m.comp, doc="Feed section 2 liquid flowrate per component") - def _feedside_liquid_percomponent(disj, comp): - """ - Liquid flowrate per component in the feed side section of the column. - - Parameters - ---------- - disj : Disjunct - The disjunct object for the feed side section in the column. - comp : int - The component index. - - Returns - ------- - Constraint - The constraint expression that enforces the liquid flowrate per component in the feed side section of the column is equal to the total liquid flowrate times the liquid composition. - """ - return m.L[2, n_tray, comp] == m.Ltotal[2, n_tray] * m.x[2, n_tray, comp] - - @disj.Constraint(m.comp, doc="Feed section 2 vapor flowrate per component") - def _feedside_vapor_percomponent(disj, comp): - """ - Vapor flowrate per component in the feed side section of the column. - - Parameters - ---------- - disj : Disjunct - The disjunct object for the feed side section in the column. - comp : int - The component index. - - Returns - ------- - Constraint - The constraint expression that enforces the vapor flowrate per component in the feed side section of the column is equal to the total vapor flowrate times the vapor composition. - """ - return m.V[2, n_tray, comp] == m.Vtotal[2, n_tray] * m.y[2, n_tray, comp] - - @disj.Constraint(doc="Feed section 2 liquid composition equilibrium summation") - def feedside_liquid_composition_summation(disj): - """ - Liquid composition equilibrium summation for the feed side section in the column. - - Parameters - ---------- - disj : Disjunct - The disjunct object for the feed side section in the column. - - Returns - ------- - Constraint - The constraint expression that enforces the liquid composition equilibrium summation for the feed side section in the column. - It ensures the sum of the liquid compositions is equal to 1 plus the error in the liquid composition. - """ - return sum(m.x[2, n_tray, comp] for comp in m.comp) - 1 == m.errx[2, n_tray] - - @disj.Constraint(doc="Feed section 2 vapor composition equilibrium summation") - def feedside_vapor_composition_summation(disj): - """ - Vapor composition equilibrium summation for the feed side section in the column. - - Parameters - ---------- - disj : Disjunct - The disjunct object for the feed side section in the column. - - Returns - ------- - Constraint - The constraint expression that enforces the vapor composition equilibrium summation for the feed side section in the column. - It ensures the sum of the vapor compositions is equal to 1 plus the error in the vapor composition. - """ - return sum(m.y[2, n_tray, comp] for comp in m.comp) - 1 == m.erry[2, n_tray] - + + @disj.Constraint(m.comp, doc="Bottom section 1 liquid enthalpy") + def bottom_liquid_enthalpy(disj, comp): + """ + Liquid enthalpy for the bottom section in the column. + + Parameters + ---------- + disj : Disjunct + The disjunct object for the bottom section in the column. + comp : int + The component index. + + Returns + ------- + Constraint + The constraint expression that enforces the liquid enthalpy for the bottom section in the column. + """ + return m.hl[1, n_tray, comp] == ( + m.cpc[1] + * ( + (m.T[1, n_tray] - m.Tref) * m.prop[comp, 'cpA', 1] + + (m.T[1, n_tray] ** 2 - m.Tref**2) + * m.prop[comp, 'cpB', 1] + * m.cpc2['A', 1] + / 2 + + (m.T[1, n_tray] ** 3 - m.Tref**3) + * m.prop[comp, 'cpC', 1] + * m.cpc2['B', 1] + / 3 + + (m.T[1, n_tray] ** 4 - m.Tref**4) * m.prop[comp, 'cpD', 1] / 4 + ) + / m.Hscale + ) + + @disj.Constraint(m.comp, doc="Bottom section 1 vapor enthalpy") + def bottom_vapor_enthalpy(disj, comp): + """ + Vapor enthalpy for the bottom section in the column. + + Parameters + ---------- + disj : Disjunct + The disjunct object for the bottom section in the column. + comp : int + The component index. + + Returns + ------- + Constraint + The constraint expression that enforces the vapor enthalpy for the bottom section in the column. + """ + return m.hv[1, n_tray, comp] == ( + m.cpc[2] + * ( + (m.T[1, n_tray] - m.Tref) * m.prop[comp, 'cpA', 2] + + (m.T[1, n_tray] ** 2 - m.Tref**2) + * m.prop[comp, 'cpB', 2] + * m.cpc2['A', 2] + / 2 + + (m.T[1, n_tray] ** 3 - m.Tref**3) + * m.prop[comp, 'cpC', 2] + * m.cpc2['B', 2] + / 3 + + (m.T[1, n_tray] ** 4 - m.Tref**4) * m.prop[comp, 'cpD', 2] / 4 + ) + / m.Hscale + + m.dHvap[comp] + ) + + @disj.Constraint(m.comp, doc="Bottom section 1 liquid activity coefficient") + def bottom_activity_coefficient(disj, comp): + """ + Liquid activity coefficient for the bottom section in the column equal to 1. + + Parameters + ---------- + disj : Disjunct + The disjunct object for the bottom section in the column. + comp : int + The component index. + + Returns + ------- + Constraint + The constraint expression that enforces the liquid activity coefficient for the bottom section for each component and tray in the column to be equal to 1. + """ + return m.actv[1, n_tray, comp] == 1 + + +def _build_feed_side_equations(disj, n_tray): + """ + Build the equations for the feed side section in the column. + + Parameters + ---------- + disj : Disjunct + The disjunct object for the feed side section in the column. + n_tray : int + The tray index. + + Returns + ------- + None + """ + m = disj.model() + + @disj.Constraint(m.comp, doc="Feed section 2 mass per component balances") + def _feedside_masspercomponent_balances(disj, comp): + """ + Mass per component balances for the feed side section in the column. + + Parameters + ---------- + disj : Disjunct + The disjunct object for the feed side section in the column. + comp : int + The component index. + + Returns + ------- + Constraint + The constraint expression that enforces the mass balance per component in the feed side section of the column. + """ + return (m.L[2, n_tray + 1, comp] if n_tray < m.num_trays else 0) + ( + m.L[4, 1, comp] * m.dl[2] if n_tray == m.num_trays else 0 + ) - m.L[2, n_tray, comp] + ( + m.V[1, m.num_trays, comp] * m.dv[2] if n_tray == 1 else 0 + ) + ( + m.V[2, n_tray - 1, comp] if n_tray > 1 else 0 + ) - m.V[ + 2, n_tray, comp + ] + ( + m.F[comp] if n_tray == m.feed_tray else 0 + ) == 0 + + @disj.Constraint(doc="Feed section 2 energy balances") + def _feedside_energy_balances(disj): + """ + Energy balances for the feed side section in the column. + + Parameters + ---------- + disj : Disjunct + The disjunct object for the feed side section in the column. + + Returns + ------- + Constraint + The constraint expression that enforces the energy balance for the feed side section in the column. + """ + return ( + sum( + ( + m.L[2, n_tray + 1, comp] * m.hl[2, n_tray + 1, comp] + if n_tray < m.num_trays + else 0 + ) + + ( + m.L[4, 1, comp] * m.dl[2] * m.hl[4, 1, comp] + if n_tray == m.num_trays + else 0 + ) + - m.L[2, n_tray, comp] * m.hl[2, n_tray, comp] + + ( + m.V[1, m.num_trays, comp] * m.dv[2] * m.hv[1, m.num_trays, comp] + if n_tray == 1 + else 0 + ) + + ( + m.V[2, n_tray - 1, comp] * m.hv[2, n_tray - 1, comp] + if n_tray > 1 + else 0 + ) + - m.V[2, n_tray, comp] * m.hv[2, n_tray, comp] + for comp in m.comp + ) + * m.Qscale + + sum( + ( + m.F[comp] * (m.hlf[comp] * (1 - m.q) + m.hvf[comp] * m.q) + if n_tray == m.feed_tray + else 0 + ) + for comp in m.comp + ) + * m.Qscale + == 0 + ) + + @disj.Constraint(m.comp, doc="Feed section 2 liquid flowrate per component") + def _feedside_liquid_percomponent(disj, comp): + """ + Liquid flowrate per component in the feed side section of the column. + + Parameters + ---------- + disj : Disjunct + The disjunct object for the feed side section in the column. + comp : int + The component index. + + Returns + ------- + Constraint + The constraint expression that enforces the liquid flowrate per component in the feed side section of the column is equal to the total liquid flowrate times the liquid composition. + """ + return m.L[2, n_tray, comp] == m.Ltotal[2, n_tray] * m.x[2, n_tray, comp] + + @disj.Constraint(m.comp, doc="Feed section 2 vapor flowrate per component") + def _feedside_vapor_percomponent(disj, comp): + """ + Vapor flowrate per component in the feed side section of the column. + + Parameters + ---------- + disj : Disjunct + The disjunct object for the feed side section in the column. + comp : int + The component index. + + Returns + ------- + Constraint + The constraint expression that enforces the vapor flowrate per component in the feed side section of the column is equal to the total vapor flowrate times the vapor composition. + """ + return m.V[2, n_tray, comp] == m.Vtotal[2, n_tray] * m.y[2, n_tray, comp] + + @disj.Constraint(doc="Feed section 2 liquid composition equilibrium summation") + def feedside_liquid_composition_summation(disj): + """ + Liquid composition equilibrium summation for the feed side section in the column. + + Parameters + ---------- + disj : Disjunct + The disjunct object for the feed side section in the column. + + Returns + ------- + Constraint + The constraint expression that enforces the liquid composition equilibrium summation for the feed side section in the column. + It ensures the sum of the liquid compositions is equal to 1 plus the error in the liquid composition. + """ + return sum(m.x[2, n_tray, comp] for comp in m.comp) - 1 == m.errx[2, n_tray] + + @disj.Constraint(doc="Feed section 2 vapor composition equilibrium summation") + def feedside_vapor_composition_summation(disj): + """ + Vapor composition equilibrium summation for the feed side section in the column. + + Parameters + ---------- + disj : Disjunct + The disjunct object for the feed side section in the column. + + Returns + ------- + Constraint + The constraint expression that enforces the vapor composition equilibrium summation for the feed side section in the column. + It ensures the sum of the vapor compositions is equal to 1 plus the error in the vapor composition. + """ + return sum(m.y[2, n_tray, comp] for comp in m.comp) - 1 == m.erry[2, n_tray] + @disj.Constraint( m.comp, doc="Feed section 2 temperature term for the vapor composition equation" ) def _feedside_reduced_temperature(disj, comp): - """Calculate the temperature term for the feedside section vapor composition. - - The constraint is used to avoid division by zero in the MINLP Hull reformulation. - - Parameters - ---------- - disj : Disjunct + """Calculate the temperature term for the feedside section vapor composition. + + The constraint is used to avoid division by zero in the MINLP Hull reformulation. + + Parameters + ---------- + disj : Disjunct The disjunct object for the feed side section in the column. - comp : int - The component index. - - Returns - ------- - Constraint - The constraint equation to calculate the temperature term for the feedside section vapor composition. - """ - return m.Tr[2, n_tray, comp] * m.T[2, n_tray] == m.prop[comp, 'TC'] - + comp : int + The component index. + + Returns + ------- + Constraint + The constraint equation to calculate the temperature term for the feedside section vapor composition. + """ + return m.Tr[2, n_tray, comp] * m.T[2, n_tray] == m.prop[comp, 'TC'] + @disj.Constraint(m.comp, doc="Feed section 2 vapor composition") def feedside_vapor_composition(disj, comp): - """ - Vapor composition for the feed side section in the column. - - Parameters - ---------- - disj : Disjunct - The disjunct object for the feed side section in the column. - comp : int - The component index. - - Returns - ------- - Constraint - The constraint expression that enforces the vapor composition for the feed side section in the column. - The equation is derived from the vapor-liquid equilibrium relationship. - """ + """ + Vapor composition for the feed side section in the column. + + Parameters + ---------- + disj : Disjunct + The disjunct object for the feed side section in the column. + comp : int + The component index. + + Returns + ------- + Constraint + The constraint expression that enforces the vapor composition for the feed side section in the column. + The equation is derived from the vapor-liquid equilibrium relationship. + """ return m.y[2, n_tray, comp] * m.P[2, n_tray] == m.x[2, n_tray, comp] * ( m.actv[2, n_tray, comp] * ( @@ -1817,315 +1817,315 @@ def feedside_vapor_composition(disj, comp): ) ) ) - - @disj.Constraint(m.comp, doc="Feed section 2 liquid enthalpy") - def feedside_liquid_enthalpy(disj, comp): - """ - Liquid enthalpy for the feed side section in the column. - - Parameters - ---------- - disj : Disjunct - The disjunct object for the feed side section in the column. - comp : int - The component index. - - Returns - ------- - Constraint - The constraint expression that enforces the liquid enthalpy for the feed side section in the column. - """ - return m.hl[2, n_tray, comp] == ( - m.cpc[1] - * ( - (m.T[2, n_tray] - m.Tref) * m.prop[comp, 'cpA', 1] - + (m.T[2, n_tray] ** 2 - m.Tref**2) - * m.prop[comp, 'cpB', 1] - * m.cpc2['A', 1] - / 2 - + (m.T[2, n_tray] ** 3 - m.Tref**3) - * m.prop[comp, 'cpC', 1] - * m.cpc2['B', 1] - / 3 - + (m.T[2, n_tray] ** 4 - m.Tref**4) * m.prop[comp, 'cpD', 1] / 4 - ) - / m.Hscale - ) - - @disj.Constraint(m.comp, doc="Feed section 2 vapor enthalpy") - def feedside_vapor_enthalpy(disj, comp): - """ - Vapor enthalpy for the feed side section in the column. - - Parameters - ---------- - disj : Disjunct - The disjunct object for the feed side section in the column. - comp : int - The component index. - - Returns - ------- - Constraint - The constraint expression that enforces the vapor enthalpy for the feed side section in the column. - """ - return m.hv[2, n_tray, comp] == ( - m.cpc[2] - * ( - (m.T[2, n_tray] - m.Tref) * m.prop[comp, 'cpA', 2] - + (m.T[2, n_tray] ** 2 - m.Tref**2) - * m.prop[comp, 'cpB', 2] - * m.cpc2['A', 2] - / 2 - + (m.T[2, n_tray] ** 3 - m.Tref**3) - * m.prop[comp, 'cpC', 2] - * m.cpc2['B', 2] - / 3 - + (m.T[2, n_tray] ** 4 - m.Tref**4) * m.prop[comp, 'cpD', 2] / 4 - ) - / m.Hscale - + m.dHvap[comp] - ) - - @disj.Constraint(m.comp, doc="Feed section 2 liquid activity coefficient") - def feedside_activity_coefficient(disj, comp): - """ - Liquid activity coefficient for the feed side section in the column equal to 1. - - Parameters - ---------- - disj : Disjunct - The disjunct object for the feed side section in the column. - comp : int - The component index. - - Returns - ------- - Constraint - The constraint expression that enforces the liquid activity coefficient for the feed side section for each component and tray in the column to be equal to 1. - This is an assumption for the feed side section, since the feed is assumed to be ideal. - """ - return m.actv[2, n_tray, comp] == 1 - - -def _build_product_side_equations(disj, n_tray): - """ - Build the equations for the product side section in the column. - - Parameters - ---------- - disj : Disjunct - The disjunct object for the product side section in the column. - n_tray : int - The tray index. - - Returns - ------- - None - """ - m = disj.model() - - @disj.Constraint(m.comp, doc="Product section 3 mass per component balances") - def _productside_masspercomponent_balances(disj, comp): - """ - Mass per component balances for the product side section in the column. - - Parameters - ---------- - disj : Disjunct - The disjunct object for the product side section in the column. - comp : int - The component index. - - Returns - ------- - Constraint - The constraint expression that enforces the mass balance per component in the product side section of the column. - """ - return (m.L[3, n_tray + 1, comp] if n_tray < m.num_trays else 0) + ( - m.L[4, 1, comp] * m.dl[3] if n_tray == m.num_trays else 0 - ) - m.L[3, n_tray, comp] + ( - m.V[1, m.num_trays, comp] * m.dv[3] if n_tray == 1 else 0 - ) + ( - m.V[3, n_tray - 1, comp] if n_tray > 1 else 0 - ) - m.V[ - 3, n_tray, comp - ] - ( - m.S[1, comp] if n_tray == m.sideout1_tray else 0 - ) - ( - m.S[2, comp] if n_tray == m.sideout2_tray else 0 - ) == 0 - - @disj.Constraint(doc="Product section 3 energy balances") - def _productside_energy_balances(disj): - """ - Energy balances for the product side section in the column. - - Parameters - ---------- - disj : Disjunct - The disjunct object for the product side section in the column. - - Returns - ------- - Constraint - The constraint expression that enforces the energy balance for the product side section in the column. - """ - return ( - sum( - ( - m.L[3, n_tray + 1, comp] * m.hl[3, n_tray + 1, comp] - if n_tray < m.num_trays - else 0 - ) - + ( - m.L[4, 1, comp] * m.dl[3] * m.hl[4, 1, comp] - if n_tray == m.num_trays - else 0 - ) - - m.L[3, n_tray, comp] * m.hl[3, n_tray, comp] - + ( - m.V[1, m.num_trays, comp] * m.dv[3] * m.hv[1, m.num_trays, comp] - if n_tray == 1 - else 0 - ) - + ( - m.V[3, n_tray - 1, comp] * m.hv[3, n_tray - 1, comp] - if n_tray > 1 - else 0 - ) - - m.V[3, n_tray, comp] * m.hv[3, n_tray, comp] - - ( - m.S[1, comp] * m.hl[3, n_tray, comp] - if n_tray == m.sideout1_tray - else 0 - ) - - ( - m.S[2, comp] * m.hl[3, n_tray, comp] - if n_tray == m.sideout2_tray - else 0 - ) - for comp in m.comp - ) - * m.Qscale - == 0 - ) - - @disj.Constraint(m.comp, doc="Product section 3 liquid flowrate per component") - def _productside_liquid_percomponent(disj, comp): - """ - Liquid flowrate per component in the product side section of the column. - - Parameters - ---------- - disj : Disjunct - The disjunct object for the product side section in the column. - comp : int - The component index. - - Returns - ------- - Constraint - The constraint expression that enforces the liquid flowrate per component in the product side section of the column is equal to the total liquid flowrate times the liquid composition. - """ - return m.L[3, n_tray, comp] == m.Ltotal[3, n_tray] * m.x[3, n_tray, comp] - - @disj.Constraint(m.comp, doc="Product section 3 vapor flowrate per component") - def _productside_vapor_percomponent(disj, comp): - """ - Vapor flowrate per component in the product side section of the column. - - Parameters - ---------- - disj : Disjunct - The disjunct object for the product side section in the column. - comp : int - The component index. - - Returns - ------- - Constraint - The constraint expression that enforces the vapor flowrate per component in the product side section of the column is equal to the total vapor flowrate times the vapor composition. - """ - return m.V[3, n_tray, comp] == m.Vtotal[3, n_tray] * m.y[3, n_tray, comp] - - @disj.Constraint(doc="Product section 3 liquid composition equilibrium summation") - def productside_liquid_composition_summation(disj): - """ - Liquid composition equilibrium summation for the product side section in the column. - - Parameters - ---------- - disj : Disjunct - The disjunct object for the product side section in the column. - - Returns - ------- - Constraint - The constraint expression that enforces the liquid composition equilibrium summation for the product side section in the column. - It ensures the sum of the liquid compositions is equal to 1 plus the error in the liquid composition. - """ - return sum(m.x[3, n_tray, comp] for comp in m.comp) - 1 == m.errx[3, n_tray] - - @disj.Constraint(doc="Product section 3 vapor composition equilibrium summation") - def productside_vapor_composition_summation(disj): - """ - Vapor composition equilibrium summation for the product side section in the column. - - Parameters - ---------- - disj : Disjunct - The disjunct object for the product side section in the column. - - Returns - ------- - Constraint - The constraint expression that enforces the vapor composition equilibrium summation for the product side section in the column. - It ensures the sum of the vapor compositions is equal to 1 plus the error in the vapor composition. - """ - return sum(m.y[3, n_tray, comp] for comp in m.comp) - 1 == m.erry[3, n_tray] - + + @disj.Constraint(m.comp, doc="Feed section 2 liquid enthalpy") + def feedside_liquid_enthalpy(disj, comp): + """ + Liquid enthalpy for the feed side section in the column. + + Parameters + ---------- + disj : Disjunct + The disjunct object for the feed side section in the column. + comp : int + The component index. + + Returns + ------- + Constraint + The constraint expression that enforces the liquid enthalpy for the feed side section in the column. + """ + return m.hl[2, n_tray, comp] == ( + m.cpc[1] + * ( + (m.T[2, n_tray] - m.Tref) * m.prop[comp, 'cpA', 1] + + (m.T[2, n_tray] ** 2 - m.Tref**2) + * m.prop[comp, 'cpB', 1] + * m.cpc2['A', 1] + / 2 + + (m.T[2, n_tray] ** 3 - m.Tref**3) + * m.prop[comp, 'cpC', 1] + * m.cpc2['B', 1] + / 3 + + (m.T[2, n_tray] ** 4 - m.Tref**4) * m.prop[comp, 'cpD', 1] / 4 + ) + / m.Hscale + ) + + @disj.Constraint(m.comp, doc="Feed section 2 vapor enthalpy") + def feedside_vapor_enthalpy(disj, comp): + """ + Vapor enthalpy for the feed side section in the column. + + Parameters + ---------- + disj : Disjunct + The disjunct object for the feed side section in the column. + comp : int + The component index. + + Returns + ------- + Constraint + The constraint expression that enforces the vapor enthalpy for the feed side section in the column. + """ + return m.hv[2, n_tray, comp] == ( + m.cpc[2] + * ( + (m.T[2, n_tray] - m.Tref) * m.prop[comp, 'cpA', 2] + + (m.T[2, n_tray] ** 2 - m.Tref**2) + * m.prop[comp, 'cpB', 2] + * m.cpc2['A', 2] + / 2 + + (m.T[2, n_tray] ** 3 - m.Tref**3) + * m.prop[comp, 'cpC', 2] + * m.cpc2['B', 2] + / 3 + + (m.T[2, n_tray] ** 4 - m.Tref**4) * m.prop[comp, 'cpD', 2] / 4 + ) + / m.Hscale + + m.dHvap[comp] + ) + + @disj.Constraint(m.comp, doc="Feed section 2 liquid activity coefficient") + def feedside_activity_coefficient(disj, comp): + """ + Liquid activity coefficient for the feed side section in the column equal to 1. + + Parameters + ---------- + disj : Disjunct + The disjunct object for the feed side section in the column. + comp : int + The component index. + + Returns + ------- + Constraint + The constraint expression that enforces the liquid activity coefficient for the feed side section for each component and tray in the column to be equal to 1. + This is an assumption for the feed side section, since the feed is assumed to be ideal. + """ + return m.actv[2, n_tray, comp] == 1 + + +def _build_product_side_equations(disj, n_tray): + """ + Build the equations for the product side section in the column. + + Parameters + ---------- + disj : Disjunct + The disjunct object for the product side section in the column. + n_tray : int + The tray index. + + Returns + ------- + None + """ + m = disj.model() + + @disj.Constraint(m.comp, doc="Product section 3 mass per component balances") + def _productside_masspercomponent_balances(disj, comp): + """ + Mass per component balances for the product side section in the column. + + Parameters + ---------- + disj : Disjunct + The disjunct object for the product side section in the column. + comp : int + The component index. + + Returns + ------- + Constraint + The constraint expression that enforces the mass balance per component in the product side section of the column. + """ + return (m.L[3, n_tray + 1, comp] if n_tray < m.num_trays else 0) + ( + m.L[4, 1, comp] * m.dl[3] if n_tray == m.num_trays else 0 + ) - m.L[3, n_tray, comp] + ( + m.V[1, m.num_trays, comp] * m.dv[3] if n_tray == 1 else 0 + ) + ( + m.V[3, n_tray - 1, comp] if n_tray > 1 else 0 + ) - m.V[ + 3, n_tray, comp + ] - ( + m.S[1, comp] if n_tray == m.sideout1_tray else 0 + ) - ( + m.S[2, comp] if n_tray == m.sideout2_tray else 0 + ) == 0 + + @disj.Constraint(doc="Product section 3 energy balances") + def _productside_energy_balances(disj): + """ + Energy balances for the product side section in the column. + + Parameters + ---------- + disj : Disjunct + The disjunct object for the product side section in the column. + + Returns + ------- + Constraint + The constraint expression that enforces the energy balance for the product side section in the column. + """ + return ( + sum( + ( + m.L[3, n_tray + 1, comp] * m.hl[3, n_tray + 1, comp] + if n_tray < m.num_trays + else 0 + ) + + ( + m.L[4, 1, comp] * m.dl[3] * m.hl[4, 1, comp] + if n_tray == m.num_trays + else 0 + ) + - m.L[3, n_tray, comp] * m.hl[3, n_tray, comp] + + ( + m.V[1, m.num_trays, comp] * m.dv[3] * m.hv[1, m.num_trays, comp] + if n_tray == 1 + else 0 + ) + + ( + m.V[3, n_tray - 1, comp] * m.hv[3, n_tray - 1, comp] + if n_tray > 1 + else 0 + ) + - m.V[3, n_tray, comp] * m.hv[3, n_tray, comp] + - ( + m.S[1, comp] * m.hl[3, n_tray, comp] + if n_tray == m.sideout1_tray + else 0 + ) + - ( + m.S[2, comp] * m.hl[3, n_tray, comp] + if n_tray == m.sideout2_tray + else 0 + ) + for comp in m.comp + ) + * m.Qscale + == 0 + ) + + @disj.Constraint(m.comp, doc="Product section 3 liquid flowrate per component") + def _productside_liquid_percomponent(disj, comp): + """ + Liquid flowrate per component in the product side section of the column. + + Parameters + ---------- + disj : Disjunct + The disjunct object for the product side section in the column. + comp : int + The component index. + + Returns + ------- + Constraint + The constraint expression that enforces the liquid flowrate per component in the product side section of the column is equal to the total liquid flowrate times the liquid composition. + """ + return m.L[3, n_tray, comp] == m.Ltotal[3, n_tray] * m.x[3, n_tray, comp] + + @disj.Constraint(m.comp, doc="Product section 3 vapor flowrate per component") + def _productside_vapor_percomponent(disj, comp): + """ + Vapor flowrate per component in the product side section of the column. + + Parameters + ---------- + disj : Disjunct + The disjunct object for the product side section in the column. + comp : int + The component index. + + Returns + ------- + Constraint + The constraint expression that enforces the vapor flowrate per component in the product side section of the column is equal to the total vapor flowrate times the vapor composition. + """ + return m.V[3, n_tray, comp] == m.Vtotal[3, n_tray] * m.y[3, n_tray, comp] + + @disj.Constraint(doc="Product section 3 liquid composition equilibrium summation") + def productside_liquid_composition_summation(disj): + """ + Liquid composition equilibrium summation for the product side section in the column. + + Parameters + ---------- + disj : Disjunct + The disjunct object for the product side section in the column. + + Returns + ------- + Constraint + The constraint expression that enforces the liquid composition equilibrium summation for the product side section in the column. + It ensures the sum of the liquid compositions is equal to 1 plus the error in the liquid composition. + """ + return sum(m.x[3, n_tray, comp] for comp in m.comp) - 1 == m.errx[3, n_tray] + + @disj.Constraint(doc="Product section 3 vapor composition equilibrium summation") + def productside_vapor_composition_summation(disj): + """ + Vapor composition equilibrium summation for the product side section in the column. + + Parameters + ---------- + disj : Disjunct + The disjunct object for the product side section in the column. + + Returns + ------- + Constraint + The constraint expression that enforces the vapor composition equilibrium summation for the product side section in the column. + It ensures the sum of the vapor compositions is equal to 1 plus the error in the vapor composition. + """ + return sum(m.y[3, n_tray, comp] for comp in m.comp) - 1 == m.erry[3, n_tray] + @disj.Constraint( m.comp, doc="Product section 3 temperature term for the vapor composition equation", ) def _productside_reduced_temperature(disj, comp): - """Calculate the temperature term for the product section vapor composition. - - The constraint is used to avoid division by zero in the MINLP Hull reformulation. - - Parameters - ---------- - disj : Disjunct + """Calculate the temperature term for the product section vapor composition. + + The constraint is used to avoid division by zero in the MINLP Hull reformulation. + + Parameters + ---------- + disj : Disjunct The disjunct object for the product side section in the column. - comp : int - The component index. - - Returns - ------- - Constraint - The constraint equation to calculate the temperature term for the product section vapor composition. - """ - return m.Tr[3, n_tray, comp] * m.T[3, n_tray] == m.prop[comp, 'TC'] - + comp : int + The component index. + + Returns + ------- + Constraint + The constraint equation to calculate the temperature term for the product section vapor composition. + """ + return m.Tr[3, n_tray, comp] * m.T[3, n_tray] == m.prop[comp, 'TC'] + @disj.Constraint(m.comp, doc="Product section 3 vapor composition") def productside_vapor_composition(disj, comp): - """ - Vapor composition for the product side section in the column. - - Parameters - ---------- - disj : Disjunct - The disjunct object for the product side section in the column. - comp : int - The component index. - - Returns - ------- - Constraint - The constraint expression that enforces the vapor composition for the product side section in the column. - The equation is derived from the vapor-liquid equilibrium relationship. - """ + """ + Vapor composition for the product side section in the column. + + Parameters + ---------- + disj : Disjunct + The disjunct object for the product side section in the column. + comp : int + The component index. + + Returns + ------- + Constraint + The constraint expression that enforces the vapor composition for the product side section in the column. + The equation is derived from the vapor-liquid equilibrium relationship. + """ return m.y[3, n_tray, comp] * m.P[3, n_tray] == m.x[3, n_tray, comp] * ( m.actv[3, n_tray, comp] * ( @@ -2144,322 +2144,322 @@ def productside_vapor_composition(disj, comp): ) ) ) - - @disj.Constraint(m.comp, doc="Product section 3 liquid enthalpy") - def productside_liquid_enthalpy(disj, comp): - """ - Liquid enthalpy for the product side section in the column. - - Parameters - ---------- - disj : Disjunct - The disjunct object for the product side section in the column. - comp : int - The component index. - - Returns - ------- - Constraint - The constraint expression that enforces the liquid enthalpy for the product side section in the column. - """ - return m.hl[3, n_tray, comp] == ( - m.cpc[1] - * ( - (m.T[3, n_tray] - m.Tref) * m.prop[comp, 'cpA', 1] - + (m.T[3, n_tray] ** 2 - m.Tref**2) - * m.prop[comp, 'cpB', 1] - * m.cpc2['A', 1] - / 2 - + (m.T[3, n_tray] ** 3 - m.Tref**3) - * m.prop[comp, 'cpC', 1] - * m.cpc2['B', 1] - / 3 - + (m.T[3, n_tray] ** 4 - m.Tref**4) * m.prop[comp, 'cpD', 1] / 4 - ) - / m.Hscale - ) - - @disj.Constraint(m.comp, doc="Product section 3 vapor enthalpy") - def productside_vapor_enthalpy(disj, comp): - """ - Vapor enthalpy for the product side section in the column. - - Parameters - ---------- - disj : Disjunct - The disjunct object for the product side section in the column. - comp : int - The component index. - - Returns - ------- - Constraint - The constraint expression that enforces the vapor enthalpy for the product side section in the column. - """ - return m.hv[3, n_tray, comp] == ( - m.cpc[2] - * ( - (m.T[3, n_tray] - m.Tref) * m.prop[comp, 'cpA', 2] - + (m.T[3, n_tray] ** 2 - m.Tref**2) - * m.prop[comp, 'cpB', 2] - * m.cpc2['A', 2] - / 2 - + (m.T[3, n_tray] ** 3 - m.Tref**3) - * m.prop[comp, 'cpC', 2] - * m.cpc2['B', 2] - / 3 - + (m.T[3, n_tray] ** 4 - m.Tref**4) * m.prop[comp, 'cpD', 2] / 4 - ) - / m.Hscale - + m.dHvap[comp] - ) - - @disj.Constraint(m.comp, doc="Product section 3 liquid activity coefficient") - def productside_activity_coefficient(disj, comp): - """ - Liquid activity coefficient for the product side section in the column equal to 1. - - Parameters - ---------- - disj : Disjunct - The disjunct object for the product side section in the column. - comp : int - The component index. - - Returns - ------- - Constraint - The constraint expression that enforces the liquid activity coefficient for the product side section for each component and tray in the column to be equal to 1. - This is an assumption for the product side section, since the product is assumed to be ideal. - """ - return m.actv[3, n_tray, comp] == 1 - - -def _build_top_equations(disj, n_tray): - """ - Build the equations for the top section in the column. - - Parameters - ---------- - disj : Disjunct - The disjunct object for the top section in the column. - n_tray : int - The tray index. - - Returns - ------- - None - """ - m = disj.model() - - @disj.Constraint(m.comp, doc="Top section 4 mass per component balances") - def _top_mass_percomponent_balances(disj, comp): - """ - Mass per component balances for the top section in the column. - - Parameters - ---------- - disj : Disjunct - The disjunct object for the top section in the column. - comp : int - The component index. - - Returns - ------- - Constraint - The constraint expression that enforces the mass balance per component in the top section of the column. - """ - return (m.L[4, n_tray + 1, comp] if n_tray < m.con_tray else 0) - ( - m.L[4, n_tray, comp] * m.dl[2] if n_tray == 1 else 0 - ) - (m.L[4, n_tray, comp] * m.dl[3] if n_tray == 1 else 0) - ( - m.L[4, n_tray, comp] if n_tray > 1 else 0 - ) + ( - m.V[2, m.num_trays, comp] if n_tray == 1 else 0 - ) + ( - m.V[3, m.num_trays, comp] if n_tray == 1 else 0 - ) + ( - m.V[4, n_tray - 1, comp] if n_tray > 1 else 0 - ) - ( - m.V[4, n_tray, comp] if n_tray < m.con_tray else 0 - ) - ( - m.D[comp] if n_tray == m.con_tray else 0 - ) == 0 - + + @disj.Constraint(m.comp, doc="Product section 3 liquid enthalpy") + def productside_liquid_enthalpy(disj, comp): + """ + Liquid enthalpy for the product side section in the column. + + Parameters + ---------- + disj : Disjunct + The disjunct object for the product side section in the column. + comp : int + The component index. + + Returns + ------- + Constraint + The constraint expression that enforces the liquid enthalpy for the product side section in the column. + """ + return m.hl[3, n_tray, comp] == ( + m.cpc[1] + * ( + (m.T[3, n_tray] - m.Tref) * m.prop[comp, 'cpA', 1] + + (m.T[3, n_tray] ** 2 - m.Tref**2) + * m.prop[comp, 'cpB', 1] + * m.cpc2['A', 1] + / 2 + + (m.T[3, n_tray] ** 3 - m.Tref**3) + * m.prop[comp, 'cpC', 1] + * m.cpc2['B', 1] + / 3 + + (m.T[3, n_tray] ** 4 - m.Tref**4) * m.prop[comp, 'cpD', 1] / 4 + ) + / m.Hscale + ) + + @disj.Constraint(m.comp, doc="Product section 3 vapor enthalpy") + def productside_vapor_enthalpy(disj, comp): + """ + Vapor enthalpy for the product side section in the column. + + Parameters + ---------- + disj : Disjunct + The disjunct object for the product side section in the column. + comp : int + The component index. + + Returns + ------- + Constraint + The constraint expression that enforces the vapor enthalpy for the product side section in the column. + """ + return m.hv[3, n_tray, comp] == ( + m.cpc[2] + * ( + (m.T[3, n_tray] - m.Tref) * m.prop[comp, 'cpA', 2] + + (m.T[3, n_tray] ** 2 - m.Tref**2) + * m.prop[comp, 'cpB', 2] + * m.cpc2['A', 2] + / 2 + + (m.T[3, n_tray] ** 3 - m.Tref**3) + * m.prop[comp, 'cpC', 2] + * m.cpc2['B', 2] + / 3 + + (m.T[3, n_tray] ** 4 - m.Tref**4) * m.prop[comp, 'cpD', 2] / 4 + ) + / m.Hscale + + m.dHvap[comp] + ) + + @disj.Constraint(m.comp, doc="Product section 3 liquid activity coefficient") + def productside_activity_coefficient(disj, comp): + """ + Liquid activity coefficient for the product side section in the column equal to 1. + + Parameters + ---------- + disj : Disjunct + The disjunct object for the product side section in the column. + comp : int + The component index. + + Returns + ------- + Constraint + The constraint expression that enforces the liquid activity coefficient for the product side section for each component and tray in the column to be equal to 1. + This is an assumption for the product side section, since the product is assumed to be ideal. + """ + return m.actv[3, n_tray, comp] == 1 + + +def _build_top_equations(disj, n_tray): + """ + Build the equations for the top section in the column. + + Parameters + ---------- + disj : Disjunct + The disjunct object for the top section in the column. + n_tray : int + The tray index. + + Returns + ------- + None + """ + m = disj.model() + + @disj.Constraint(m.comp, doc="Top section 4 mass per component balances") + def _top_mass_percomponent_balances(disj, comp): + """ + Mass per component balances for the top section in the column. + + Parameters + ---------- + disj : Disjunct + The disjunct object for the top section in the column. + comp : int + The component index. + + Returns + ------- + Constraint + The constraint expression that enforces the mass balance per component in the top section of the column. + """ + return (m.L[4, n_tray + 1, comp] if n_tray < m.con_tray else 0) - ( + m.L[4, n_tray, comp] * m.dl[2] if n_tray == 1 else 0 + ) - (m.L[4, n_tray, comp] * m.dl[3] if n_tray == 1 else 0) - ( + m.L[4, n_tray, comp] if n_tray > 1 else 0 + ) + ( + m.V[2, m.num_trays, comp] if n_tray == 1 else 0 + ) + ( + m.V[3, m.num_trays, comp] if n_tray == 1 else 0 + ) + ( + m.V[4, n_tray - 1, comp] if n_tray > 1 else 0 + ) - ( + m.V[4, n_tray, comp] if n_tray < m.con_tray else 0 + ) - ( + m.D[comp] if n_tray == m.con_tray else 0 + ) == 0 + @disj.Constraint(doc="Top section 4 energy balances") - def _top_energy_balances(disj): - """ - Energy balances for the top section in the column. - - Parameters - ---------- - disj : Disjunct - The disjunct object for the top section in the column. - - Returns - ------- - Constraint - The constraint expression that enforces the energy balance for the top section in the column. - """ - return ( - sum( - ( - m.L[4, n_tray + 1, comp] * m.hl[4, n_tray + 1, comp] - if n_tray < m.con_tray - else 0 - ) - - ( - m.L[4, n_tray, comp] * m.dl[2] * m.hl[4, n_tray, comp] - if n_tray == 1 - else 0 - ) - - ( - m.L[4, n_tray, comp] * m.dl[3] * m.hl[4, n_tray, comp] - if n_tray == 1 - else 0 - ) - - (m.L[4, n_tray, comp] * m.hl[4, n_tray, comp] if n_tray > 1 else 0) - + ( - m.V[2, m.num_trays, comp] * m.hv[2, m.num_trays, comp] - if n_tray == 1 - else 0 - ) - + ( - m.V[3, m.num_trays, comp] * m.hv[3, m.num_trays, comp] - if n_tray == 1 - else 0 - ) - + ( - m.V[4, n_tray - 1, comp] * m.hv[4, n_tray - 1, comp] - if n_tray > 1 - else 0 - ) - - ( - m.V[4, n_tray, comp] * m.hv[4, n_tray, comp] - if n_tray < m.con_tray - else 0 - ) - - (m.D[comp] * m.hl[4, n_tray, comp] if n_tray == m.con_tray else 0) - for comp in m.comp - ) - * m.Qscale - - (m.Qcon if n_tray == m.con_tray else 0) - == 0 - ) - - @disj.Constraint(m.comp, doc="Top section 4 liquid flowrate per component") - def _top_liquid_percomponent(disj, comp): - """ - Liquid flowrate per component in the top section of the column. - - Parameters - ---------- - disj : Disjunct - The disjunct object for the top section in the column. - comp : int - The component index. - - Returns - ------- - Constraint - The constraint expression that enforces the liquid flowrate per component in the top section of the column is equal to the total liquid flowrate times the liquid composition. - """ - return m.L[4, n_tray, comp] == m.Ltotal[4, n_tray] * m.x[4, n_tray, comp] - - @disj.Constraint(m.comp, doc="Top section 4 vapor flowrate per component") - def _top_vapor_percomponent(disj, comp): - """ - Vapor flowrate per component in the top section of the column. - - Parameters - ---------- - disj : Disjunct - The disjunct object for the top section in the column. - comp : int - The component index. - - Returns - ------- - Constraint - The constraint expression that enforces the vapor flowrate per component in the top section of the column is equal to the total vapor flowrate times the vapor composition. - """ - return m.V[4, n_tray, comp] == m.Vtotal[4, n_tray] * m.y[4, n_tray, comp] - - @disj.Constraint(doc="Top section 4 liquid composition equilibrium summation") - def top_liquid_composition_summation(disj): - """ - Liquid composition equilibrium summation for the top section in the column. - - Parameters - ---------- - disj : Disjunct - The disjunct object for the top section in the column. - - Returns - ------- - Constraint - The constraint expression that enforces the liquid composition equilibrium summation for the top section in the column. - It ensures the sum of the liquid compositions is equal to 1 plus the error in the liquid composition. - """ - return sum(m.x[4, n_tray, comp] for comp in m.comp) - 1 == m.errx[4, n_tray] - - @disj.Constraint(doc="Top section 4 vapor composition equilibrium summation") - def top_vapor_composition_summation(disj): - """ - Vapor composition equilibrium summation for the top section in the column. - - Parameters - ---------- - disj : Disjunct - The disjunct object for the top section in the column. - - Returns - ------- - Constraint - The constraint expression that enforces the vapor composition equilibrium summation for the top section in the column. - It ensures the sum of the vapor compositions is equal to 1 plus the error in the vapor composition. - """ - return sum(m.y[4, n_tray, comp] for comp in m.comp) - 1 == m.erry[4, n_tray] - + def _top_energy_balances(disj): + """ + Energy balances for the top section in the column. + + Parameters + ---------- + disj : Disjunct + The disjunct object for the top section in the column. + + Returns + ------- + Constraint + The constraint expression that enforces the energy balance for the top section in the column. + """ + return ( + sum( + ( + m.L[4, n_tray + 1, comp] * m.hl[4, n_tray + 1, comp] + if n_tray < m.con_tray + else 0 + ) + - ( + m.L[4, n_tray, comp] * m.dl[2] * m.hl[4, n_tray, comp] + if n_tray == 1 + else 0 + ) + - ( + m.L[4, n_tray, comp] * m.dl[3] * m.hl[4, n_tray, comp] + if n_tray == 1 + else 0 + ) + - (m.L[4, n_tray, comp] * m.hl[4, n_tray, comp] if n_tray > 1 else 0) + + ( + m.V[2, m.num_trays, comp] * m.hv[2, m.num_trays, comp] + if n_tray == 1 + else 0 + ) + + ( + m.V[3, m.num_trays, comp] * m.hv[3, m.num_trays, comp] + if n_tray == 1 + else 0 + ) + + ( + m.V[4, n_tray - 1, comp] * m.hv[4, n_tray - 1, comp] + if n_tray > 1 + else 0 + ) + - ( + m.V[4, n_tray, comp] * m.hv[4, n_tray, comp] + if n_tray < m.con_tray + else 0 + ) + - (m.D[comp] * m.hl[4, n_tray, comp] if n_tray == m.con_tray else 0) + for comp in m.comp + ) + * m.Qscale + - (m.Qcon if n_tray == m.con_tray else 0) + == 0 + ) + + @disj.Constraint(m.comp, doc="Top section 4 liquid flowrate per component") + def _top_liquid_percomponent(disj, comp): + """ + Liquid flowrate per component in the top section of the column. + + Parameters + ---------- + disj : Disjunct + The disjunct object for the top section in the column. + comp : int + The component index. + + Returns + ------- + Constraint + The constraint expression that enforces the liquid flowrate per component in the top section of the column is equal to the total liquid flowrate times the liquid composition. + """ + return m.L[4, n_tray, comp] == m.Ltotal[4, n_tray] * m.x[4, n_tray, comp] + + @disj.Constraint(m.comp, doc="Top section 4 vapor flowrate per component") + def _top_vapor_percomponent(disj, comp): + """ + Vapor flowrate per component in the top section of the column. + + Parameters + ---------- + disj : Disjunct + The disjunct object for the top section in the column. + comp : int + The component index. + + Returns + ------- + Constraint + The constraint expression that enforces the vapor flowrate per component in the top section of the column is equal to the total vapor flowrate times the vapor composition. + """ + return m.V[4, n_tray, comp] == m.Vtotal[4, n_tray] * m.y[4, n_tray, comp] + + @disj.Constraint(doc="Top section 4 liquid composition equilibrium summation") + def top_liquid_composition_summation(disj): + """ + Liquid composition equilibrium summation for the top section in the column. + + Parameters + ---------- + disj : Disjunct + The disjunct object for the top section in the column. + + Returns + ------- + Constraint + The constraint expression that enforces the liquid composition equilibrium summation for the top section in the column. + It ensures the sum of the liquid compositions is equal to 1 plus the error in the liquid composition. + """ + return sum(m.x[4, n_tray, comp] for comp in m.comp) - 1 == m.errx[4, n_tray] + + @disj.Constraint(doc="Top section 4 vapor composition equilibrium summation") + def top_vapor_composition_summation(disj): + """ + Vapor composition equilibrium summation for the top section in the column. + + Parameters + ---------- + disj : Disjunct + The disjunct object for the top section in the column. + + Returns + ------- + Constraint + The constraint expression that enforces the vapor composition equilibrium summation for the top section in the column. + It ensures the sum of the vapor compositions is equal to 1 plus the error in the vapor composition. + """ + return sum(m.y[4, n_tray, comp] for comp in m.comp) - 1 == m.erry[4, n_tray] + @disj.Constraint( m.comp, doc="Top section 4 temperature term for the vapor composition equation" ) def _top_reduced_temperature(disj, comp): - """Calculate the temperature term for the top section vapor composition. - - The constraint is used to avoid division by zero in the MINLP Hull reformulation. - - Parameters - ---------- - disj : Disjunct + """Calculate the temperature term for the top section vapor composition. + + The constraint is used to avoid division by zero in the MINLP Hull reformulation. + + Parameters + ---------- + disj : Disjunct The disjunct object for the top section in the column. - comp : int - The component index. - - Returns - ------- - Constraint - The constraint equation to calculate the temperature term for the top section vapor composition. - """ - return m.Tr[4, n_tray, comp] * m.T[4, n_tray] == m.prop[comp, 'TC'] - + comp : int + The component index. + + Returns + ------- + Constraint + The constraint equation to calculate the temperature term for the top section vapor composition. + """ + return m.Tr[4, n_tray, comp] * m.T[4, n_tray] == m.prop[comp, 'TC'] + @disj.Constraint(m.comp, doc="Top section 4 vapor composition") - def top_vapor_composition(disj, comp): - """ - Vapor composition for the top section in the column. - - Parameters - ---------- - disj : Disjunct - The disjunct object for the top section in the column. - comp : int - The component index. - - Returns - ------- - Constraint - The constraint expression that enforces the vapor composition for the top section in the column. - The equation is derived from the vapor-liquid equilibrium relationship. - """ + def top_vapor_composition(disj, comp): + """ + Vapor composition for the top section in the column. + + Parameters + ---------- + disj : Disjunct + The disjunct object for the top section in the column. + comp : int + The component index. + + Returns + ------- + Constraint + The constraint expression that enforces the vapor composition for the top section in the column. + The equation is derived from the vapor-liquid equilibrium relationship. + """ return m.y[4, n_tray, comp] * m.P[4, n_tray] == m.x[4, n_tray, comp] * ( m.actv[4, n_tray, comp] * ( @@ -2478,252 +2478,252 @@ def top_vapor_composition(disj, comp): ) ) ) - - @disj.Constraint(m.comp, doc="Top section 4 liquid enthalpy") - def top_liquid_enthalpy(disj, comp): - """ - Liquid enthalpy for the top section in the column. - - Parameters - ---------- - disj : Disjunct - The disjunct object for the top section in the column. - comp : int - The component index. - - Returns - ------- - Constraint - The constraint expression that enforces the liquid enthalpy for the top section in the column. - """ - return m.hl[4, n_tray, comp] == ( - m.cpc[1] - * ( - (m.T[4, n_tray] - m.Tref) * m.prop[comp, 'cpA', 1] - + (m.T[4, n_tray] ** 2 - m.Tref**2) - * m.prop[comp, 'cpB', 1] - * m.cpc2['A', 1] - / 2 - + (m.T[4, n_tray] ** 3 - m.Tref**3) - * m.prop[comp, 'cpC', 1] - * m.cpc2['B', 1] - / 3 - + (m.T[4, n_tray] ** 4 - m.Tref**4) * m.prop[comp, 'cpD', 1] / 4 - ) - / m.Hscale - ) - - @disj.Constraint(m.comp, doc="Top section 4 vapor enthalpy") - def top_vapor_enthalpy(disj, comp): - """ - Vapor enthalpy for the top section in the column. - - Parameters - ---------- - disj : Disjunct - The disjunct object for the top section in the column. - comp : int - The component index. - - Returns - ------- - Constraint - The constraint expression that enforces the vapor enthalpy for the top section in the column. - """ - return m.hv[4, n_tray, comp] == ( - m.cpc[2] - * ( - (m.T[4, n_tray] - m.Tref) * m.prop[comp, 'cpA', 2] - + (m.T[4, n_tray] ** 2 - m.Tref**2) - * m.prop[comp, 'cpB', 2] - * m.cpc2['A', 2] - / 2 - + (m.T[4, n_tray] ** 3 - m.Tref**3) - * m.prop[comp, 'cpC', 2] - * m.cpc2['B', 2] - / 3 - + (m.T[4, n_tray] ** 4 - m.Tref**4) * m.prop[comp, 'cpD', 2] / 4 - ) - / m.Hscale - + m.dHvap[comp] - ) - - @disj.Constraint(m.comp, doc="Top section 4 liquid activity coefficient") - def top_activity_coefficient(disj, comp): - """ - Liquid activity coefficient for the top section in the column equal to 1. - - Parameters - ---------- - disj : Disjunct - The disjunct object for the top section in the column. - comp : int - The component index. - - Returns - ------- - Constraint - The constraint expression that enforces the liquid activity coefficient for the top section for each component and tray in the column to be equal to 1. - This is an assumption for the top section, since the product is assumed to be ideal. - """ - return m.actv[4, n_tray, comp] == 1 - - -def _build_pass_through_eqns(disj, sec, n_tray): - """ - Build the equations for the pass through section in the column when a given tray in the disjunct is not active if it is the first or last tray. - - Parameters - ---------- - disj : Disjunct - The disjunct object for the pass through section in the column. - sec : int - The section index. - n_tray : int - The tray index. - - Returns - ------- - None - """ - m = disj.model() - - # If the tray is the first or last tray, then the liquid and vapor flowrates, compositions, enthalpies, and temperature are passed through. - if n_tray == 1 or n_tray == m.num_trays: - return - - @disj.Constraint(m.comp, doc="Pass through liquid flowrate") - def pass_through_liquid_flowrate(disj, comp): - """ - Pass through liquid flowrate for the given tray in the column. - The constraint enforces the liquid flowrate for the given tray is equal to the liquid flowrate for the tray above it. - - Parameters - ---------- - disj : Disjunct - The disjunct object for the pass through when the tray is inactive. - comp : int - The component index. - - Returns - ------- - Constraint - The constraint expression that enforces the liquid flowrate for the given tray is equal to the liquid flowrate for the tray above it. - """ - return m.L[sec, n_tray, comp] == m.L[sec, n_tray + 1, comp] - - @disj.Constraint(m.comp, doc="Pass through vapor flowrate") - def pass_through_vapor_flowrate(disj, comp): - """ - Pass through vapor flowrate for the given tray in the column. - - Parameters - ---------- - disj : Disjunct - The disjunct object for the pass through when the tray is inactive. - comp : int - The component index. - - Returns - ------- - Constraint - The constraint expression that enforces the vapor flowrate for the given tray is equal to the vapor flowrate for the tray below it. - """ - return m.V[sec, n_tray, comp] == m.V[sec, n_tray - 1, comp] - - @disj.Constraint(m.comp, doc="Pass through liquid composition") - def pass_through_liquid_composition(disj, comp): - """ - Pass through liquid composition for the given tray in the column. - - Parameters - ---------- - disj : Disjunct - The disjunct object for the pass through when the tray is inactive. - comp : int - The component index. - - Returns - ------- - Constraint - The constraint expression that enforces the liquid composition for the given tray is equal to the liquid composition for the tray above it. - """ - return m.x[sec, n_tray, comp] == m.x[sec, n_tray + 1, comp] - - @disj.Constraint(m.comp, doc="Pass through vapor composition") - def pass_through_vapor_composition(disj, comp): - """ - Pass through vapor composition for the given tray in the column. - - Parameters - ---------- - disj : Disjunct - The disjunct object for the pass through when the tray is inactive. - comp : int - The component index. - - Returns - ------- - Constraint - The constraint expression that enforces the vapor composition for the given tray is equal to the vapor composition for the tray below it. - """ - return m.y[sec, n_tray, comp] == m.y[sec, n_tray + 1, comp] - - @disj.Constraint(m.comp, doc="Pass through liquid enthalpy") - def pass_through_liquid_enthalpy(disj, comp): - """ - Pass through liquid enthalpy for the given tray in the column. - - Parameters - ---------- - disj : Disjunct - The disjunct object for the pass through when the tray is inactive. - comp : int - The component index. - - Returns - ------- - Constraint - The constraint expression that enforces the liquid enthalpy for the given tray is equal to the liquid enthalpy for the tray above it. - """ - return m.hl[sec, n_tray, comp] == m.hl[sec, n_tray + 1, comp] - - @disj.Constraint(m.comp, doc="Pass through vapor enthalpy") - def pass_through_vapor_enthalpy(disj, comp): - """ - Pass through vapor enthalpy for the given tray in the column. - - Parameters - ---------- - disj : Disjunct - The disjunct object for the pass through when the tray is inactive. - comp : int - The component index. - - Returns - ------- - Constraint - The constraint expression that enforces the vapor enthalpy for the given tray is equal to the vapor enthalpy for the tray below it. - """ - return m.hv[sec, n_tray, comp] == m.hv[sec, n_tray - 1, comp] - - @disj.Constraint(doc="Pass through temperature") - def pass_through_temperature(disj): - """ - Pass through temperature for the given tray in the column. - - Parameters - ---------- - disj : Disjunct - The disjunct object for the pass through when the tray is inactive. - - Returns - ------- - Constraint - The constraint expression that enforces the temperature for the given tray is equal to the temperature for the tray below it. - """ - return m.T[sec, n_tray] == m.T[sec, n_tray - 1] - - -if __name__ == "__main__": - model = build_model() + + @disj.Constraint(m.comp, doc="Top section 4 liquid enthalpy") + def top_liquid_enthalpy(disj, comp): + """ + Liquid enthalpy for the top section in the column. + + Parameters + ---------- + disj : Disjunct + The disjunct object for the top section in the column. + comp : int + The component index. + + Returns + ------- + Constraint + The constraint expression that enforces the liquid enthalpy for the top section in the column. + """ + return m.hl[4, n_tray, comp] == ( + m.cpc[1] + * ( + (m.T[4, n_tray] - m.Tref) * m.prop[comp, 'cpA', 1] + + (m.T[4, n_tray] ** 2 - m.Tref**2) + * m.prop[comp, 'cpB', 1] + * m.cpc2['A', 1] + / 2 + + (m.T[4, n_tray] ** 3 - m.Tref**3) + * m.prop[comp, 'cpC', 1] + * m.cpc2['B', 1] + / 3 + + (m.T[4, n_tray] ** 4 - m.Tref**4) * m.prop[comp, 'cpD', 1] / 4 + ) + / m.Hscale + ) + + @disj.Constraint(m.comp, doc="Top section 4 vapor enthalpy") + def top_vapor_enthalpy(disj, comp): + """ + Vapor enthalpy for the top section in the column. + + Parameters + ---------- + disj : Disjunct + The disjunct object for the top section in the column. + comp : int + The component index. + + Returns + ------- + Constraint + The constraint expression that enforces the vapor enthalpy for the top section in the column. + """ + return m.hv[4, n_tray, comp] == ( + m.cpc[2] + * ( + (m.T[4, n_tray] - m.Tref) * m.prop[comp, 'cpA', 2] + + (m.T[4, n_tray] ** 2 - m.Tref**2) + * m.prop[comp, 'cpB', 2] + * m.cpc2['A', 2] + / 2 + + (m.T[4, n_tray] ** 3 - m.Tref**3) + * m.prop[comp, 'cpC', 2] + * m.cpc2['B', 2] + / 3 + + (m.T[4, n_tray] ** 4 - m.Tref**4) * m.prop[comp, 'cpD', 2] / 4 + ) + / m.Hscale + + m.dHvap[comp] + ) + + @disj.Constraint(m.comp, doc="Top section 4 liquid activity coefficient") + def top_activity_coefficient(disj, comp): + """ + Liquid activity coefficient for the top section in the column equal to 1. + + Parameters + ---------- + disj : Disjunct + The disjunct object for the top section in the column. + comp : int + The component index. + + Returns + ------- + Constraint + The constraint expression that enforces the liquid activity coefficient for the top section for each component and tray in the column to be equal to 1. + This is an assumption for the top section, since the product is assumed to be ideal. + """ + return m.actv[4, n_tray, comp] == 1 + + +def _build_pass_through_eqns(disj, sec, n_tray): + """ + Build the equations for the pass through section in the column when a given tray in the disjunct is not active if it is the first or last tray. + + Parameters + ---------- + disj : Disjunct + The disjunct object for the pass through section in the column. + sec : int + The section index. + n_tray : int + The tray index. + + Returns + ------- + None + """ + m = disj.model() + + # If the tray is the first or last tray, then the liquid and vapor flowrates, compositions, enthalpies, and temperature are passed through. + if n_tray == 1 or n_tray == m.num_trays: + return + + @disj.Constraint(m.comp, doc="Pass through liquid flowrate") + def pass_through_liquid_flowrate(disj, comp): + """ + Pass through liquid flowrate for the given tray in the column. + The constraint enforces the liquid flowrate for the given tray is equal to the liquid flowrate for the tray above it. + + Parameters + ---------- + disj : Disjunct + The disjunct object for the pass through when the tray is inactive. + comp : int + The component index. + + Returns + ------- + Constraint + The constraint expression that enforces the liquid flowrate for the given tray is equal to the liquid flowrate for the tray above it. + """ + return m.L[sec, n_tray, comp] == m.L[sec, n_tray + 1, comp] + + @disj.Constraint(m.comp, doc="Pass through vapor flowrate") + def pass_through_vapor_flowrate(disj, comp): + """ + Pass through vapor flowrate for the given tray in the column. + + Parameters + ---------- + disj : Disjunct + The disjunct object for the pass through when the tray is inactive. + comp : int + The component index. + + Returns + ------- + Constraint + The constraint expression that enforces the vapor flowrate for the given tray is equal to the vapor flowrate for the tray below it. + """ + return m.V[sec, n_tray, comp] == m.V[sec, n_tray - 1, comp] + + @disj.Constraint(m.comp, doc="Pass through liquid composition") + def pass_through_liquid_composition(disj, comp): + """ + Pass through liquid composition for the given tray in the column. + + Parameters + ---------- + disj : Disjunct + The disjunct object for the pass through when the tray is inactive. + comp : int + The component index. + + Returns + ------- + Constraint + The constraint expression that enforces the liquid composition for the given tray is equal to the liquid composition for the tray above it. + """ + return m.x[sec, n_tray, comp] == m.x[sec, n_tray + 1, comp] + + @disj.Constraint(m.comp, doc="Pass through vapor composition") + def pass_through_vapor_composition(disj, comp): + """ + Pass through vapor composition for the given tray in the column. + + Parameters + ---------- + disj : Disjunct + The disjunct object for the pass through when the tray is inactive. + comp : int + The component index. + + Returns + ------- + Constraint + The constraint expression that enforces the vapor composition for the given tray is equal to the vapor composition for the tray below it. + """ + return m.y[sec, n_tray, comp] == m.y[sec, n_tray + 1, comp] + + @disj.Constraint(m.comp, doc="Pass through liquid enthalpy") + def pass_through_liquid_enthalpy(disj, comp): + """ + Pass through liquid enthalpy for the given tray in the column. + + Parameters + ---------- + disj : Disjunct + The disjunct object for the pass through when the tray is inactive. + comp : int + The component index. + + Returns + ------- + Constraint + The constraint expression that enforces the liquid enthalpy for the given tray is equal to the liquid enthalpy for the tray above it. + """ + return m.hl[sec, n_tray, comp] == m.hl[sec, n_tray + 1, comp] + + @disj.Constraint(m.comp, doc="Pass through vapor enthalpy") + def pass_through_vapor_enthalpy(disj, comp): + """ + Pass through vapor enthalpy for the given tray in the column. + + Parameters + ---------- + disj : Disjunct + The disjunct object for the pass through when the tray is inactive. + comp : int + The component index. + + Returns + ------- + Constraint + The constraint expression that enforces the vapor enthalpy for the given tray is equal to the vapor enthalpy for the tray below it. + """ + return m.hv[sec, n_tray, comp] == m.hv[sec, n_tray - 1, comp] + + @disj.Constraint(doc="Pass through temperature") + def pass_through_temperature(disj): + """ + Pass through temperature for the given tray in the column. + + Parameters + ---------- + disj : Disjunct + The disjunct object for the pass through when the tray is inactive. + + Returns + ------- + Constraint + The constraint expression that enforces the temperature for the given tray is equal to the temperature for the tray below it. + """ + return m.T[sec, n_tray] == m.T[sec, n_tray - 1] + + +if __name__ == "__main__": + model = build_model() diff --git a/gdplib/kaibel/main_gdp.py b/gdplib/kaibel/main_gdp.py index 0b817c90..1a164dcf 100644 --- a/gdplib/kaibel/main_gdp.py +++ b/gdplib/kaibel/main_gdp.py @@ -1,287 +1,287 @@ -""" - KAIBEL COLUMN: Modeling, Optimization, and - Conceptual Design of Multi-product Dividing Wall Columns - (written by E. Soraya Rawlings, esoraya@rwlngs.net, - in collaboration with Qi Chen) - -This is a dividing wall distillation column design problem to determine the optimal minimum number of trays and the optimal location of side streams for the separation of a quaternary mixture: - 1 = methanol - 2 = ethanol - 3 = propanol - 4 = butanol -while minimizing its capital and operating costs. - - The scheme of the Kaibel Column is shown in Figure 1: - ____ - --|Cond|--- - | ---- | - -------------- | - | sect 4 |<------> D mostly 1 - | ---------- | - | | - | ----- |----- | - | | |-------> S2 - Fj ------>| sect | sect | - | 2 | 3 | - | | |-------> S1 - | ----- |----- | - | | - | ---------- | - | sect 1 |<------> B mostly 4 - -------------- | - | ____ | - ---|Reb |--- - ---- - Figure 1. Kaibel Column scheme - -Permanent trays: -- Reboiler and vapor distributor in the bottom section (sect 1) -- Liquid distributor and condenser in the top section (sect 4) -- Side feed tray for the feed side and dividing wall starting and and ening tray in the feed section (sect 2). -- Side product trays and dividing wall starting and ending tray in the product section (sect 3). - -The trays in each section are counted from bottom to top, being tray 1 the bottom tray in each section and tray np the top tray in each section, where np is a specified upper bound for the number of possible trays for each section. -Each section has the same number of possible trays. - -Six degrees of freedom: the reflux ratio, the product outlets (bottom, intermediate, and distillate product flowrates), and the liquid and vapor flowrates between the two sections of the dividing wall, controlled by a liquid and vapor distributor on the top and bottom section of the column, respectively. -including also the vapor and liquid flowrate and the energy consumption in the reboiler and condenser. -The vapor distributor is fixed and remain constant during the column operation. - -Source paper: -Rawlings, E. S., Chen, Q., Grossmann, I. E., & Caballero, J. A. (2019). Kaibel Column: Modeling, optimization, and conceptual design of multi-product dividing wall columns. *Computers and Chemical Engineering*, 125, 31–39. https://doi.org/10.1016/j.compchemeng.2019.03.006 - -""" - -from math import fabs - -import matplotlib.pyplot as plt -from pyomo.environ import * - -# from kaibel_solve_gdp import build_model -from gdplib.kaibel.kaibel_solve_gdp import build_model - - -def main(): - """ - This is the main function that executes the optimization process. - - It builds the model, fixes certain variables, sets initial values for tray existence or absence, - solves the model using the 'gdpopt' solver, and displays the results. - - Returns: - None - """ - m = build_model() - - # Fixing variables - m.F[1].fix(50) # feed flowrate in mol/s - m.F[2].fix(50) - m.F[3].fix(50) - m.F[4].fix(50) - m.q.fix( - m.q_init - ) # vapor fraction q_init from the feed set in the build_model function - m.dv[2].fix(0.394299) # vapor distributor in the feed section - - for sec in m.section: - for n_tray in m.tray: - m.P[sec, n_tray].fix(m.Preb) - - ## Initial values for the tray existence or absence - for n_tray in m.candidate_trays_main: - for sec in m.section_main: - m.tray_exists[sec, n_tray].indicator_var.set_value(1) - m.tray_absent[sec, n_tray].indicator_var.set_value(0) - for n_tray in m.candidate_trays_feed: - m.tray_exists[2, n_tray].indicator_var.set_value(1) - m.tray_absent[2, n_tray].indicator_var.set_value(0) - for n_tray in m.candidate_trays_product: - m.tray_exists[3, n_tray].indicator_var.set_value(1) - m.tray_absent[3, n_tray].indicator_var.set_value(0) - - intro_message(m) - - results = SolverFactory('gdpopt').solve( - m, - strategy='LOA', - tee=True, - time_limit=3600, - mip_solver='gams', - mip_solver_args=dict(solver='cplex'), - ) - - m.calc_nt = sum( - sum(m.tray_exists[sec, n_tray].indicator_var.value for n_tray in m.tray) - for sec in m.section - ) - sum(m.tray_exists[3, n_tray].indicator_var.value for n_tray in m.tray) - m.dw_start = ( - sum(m.tray_exists[1, n_tray].indicator_var.value for n_tray in m.tray) + 1 - ) - m.dw_end = sum( - m.tray_exists[1, n_tray].indicator_var.value for n_tray in m.tray - ) + sum(m.tray_exists[2, n_tray].indicator_var.value for n_tray in m.tray) - - display_results(m) - - print(' ', results) - print(' Solver Status: ', results.solver.status) - print(' Termination condition: ', results.solver.termination_condition) - - -def intro_message(m): - """ - Display the introduction message. - - - - """ - print(""" - -If you use this model and/or initialization strategy, you may cite the following: -Rawlings, ES; Chen, Q; Grossmann, IE; Caballero, JA. Kaibel Column: Modeling, -optimization, and conceptual design of multi-product dividing wall columns. -Comp. and Chem. Eng., 2019, 125, 31-39. -DOI: https://doi.org/10.1016/j.compchemeng.2019.03.006 - - - """) - - -def display_results(m): - """ - Display the results of the optimization process. - - Parameters - ---------- - m : Pyomo ConcreteModel - The Pyomo model object containing the results of the optimization process. - """ - print('') - print('Components:') - print('1 methanol') - print('2 ethanol') - print('3 butanol') - print('4 propanol') - print(' ') - print('Objective: %s' % value(m.obj)) - print('Trays: %s' % value(m.calc_nt)) - print('Dividing_wall_start: %s' % value(m.dw_start)) - print('Dividing_wall_end: %s' % value(m.dw_end)) - print(' ') - print( - 'Qreb: {: >3.0f}kW B_1: {: > 2.0f} B_2: {: >2.0f} B_3: {: >2.0f} B_4: {: >2.0f} Btotal: {: >2.0f}'.format( - value(m.Qreb / m.Qscale), - value(m.B[1]), - value(m.B[2]), - value(m.B[3]), - value(m.B[4]), - value(m.Btotal), - ) - ) - print( - 'Qcon: {: >2.0f}kW D_1: {: >2.0f} D_2: {: >2.0f} D_3: {: >2.0f} D_4: {: >2.0f} Dtotal: {: >2.0f}'.format( - value(m.Qcon / m.Qscale), - value(m.D[1]), - value(m.D[2]), - value(m.D[3]), - value(m.D[4]), - value(m.Dtotal), - ) - ) - print(' ') - print('Reflux: {: >3.4f}'.format(value(m.rr))) - print('Reboil: {: >3.4f} '.format(value(m.bu))) - print(' ') - print('Flowrates[mol/s]') - print( - 'F_1: {: > 3.0f} F_2: {: >2.0f} F_3: {: >2.0f} F_4: {: >2.0f} Ftotal: {: >2.0f}'.format( - value(m.F[1]), - value(m.F[2]), - value(m.F[3]), - value(m.F[4]), - sum(value(m.F[comp]) for comp in m.comp), - ) - ) - print( - 'S1_1: {: > 1.0f} S1_2: {: >2.0f} S1_3: {: >2.0f} S1_4: {: >2.0f} S1total: {: >2.0f}'.format( - value(m.S[1, 1]), - value(m.S[1, 2]), - value(m.S[1, 3]), - value(m.S[1, 4]), - sum(value(m.S[1, comp]) for comp in m.comp), - ) - ) - print( - 'S2_1: {: > 1.0f} S2_2: {: >2.0f} S2_3: {: >2.0f} S2_4: {: >2.0f} S2total: {: >2.0f}'.format( - value(m.S[2, 1]), - value(m.S[2, 2]), - value(m.S[2, 3]), - value(m.S[2, 4]), - sum(value(m.S[2, comp]) for comp in m.comp), - ) - ) - print(' ') - print('Distributors:') - print('dl[2]: {: >3.4f} dl[3]: {: >3.4f}'.format(value(m.dl[2]), value(m.dl[3]))) - print('dv[2]: {: >3.4f} dv[3]: {: >3.4f}'.format(value(m.dv[2]), value(m.dv[3]))) - print(' ') - print(' ') - print(' ') - print(' Kaibel Column Sections ') - print('__________________________________________') - print(' ') - print(' Tray Bottom Feed ') - print('__________________________________________') - for t in reversed(list(m.tray)): - print( - '[{: >2.0f}] {: >9.0g} {: >18.0g} F:{: >3.0f} '.format( - t, - ( - fabs(value(m.tray_exists[1, t].indicator_var)) - if t in m.candidate_trays_main - else 1 - ), - ( - fabs(value(m.tray_exists[2, t].indicator_var)) - if t in m.candidate_trays_feed - else 1 - ), - sum(value(m.F[comp]) for comp in m.comp) if t == m.feed_tray else 0, - ) - ) - print(' ') - print('__________________________________________') - print(' ') - print(' Product Top ') - print('__________________________________________') - for t in reversed(list(m.tray)): - print( - '[{: >2.0f}] {: >9.0g} S1:{: >2.0f} S2:{: >2.0f} {: >8.0g}'.format( - t, - ( - fabs(value(m.tray_exists[3, t].indicator_var)) - if t in m.candidate_trays_product - else 1 - ), - ( - sum(value(m.S[1, comp]) for comp in m.comp) - if t == m.sideout1_tray - else 0 - ), - ( - sum(value(m.S[2, comp]) for comp in m.comp) - if t == m.sideout2_tray - else 0 - ), - ( - fabs(value(m.tray_exists[4, t].indicator_var)) - if t in m.candidate_trays_main - else 1 - ), - ) - ) - print(' 1 = trays exists, 0 = absent tray') - - -if __name__ == "__main__": - main() +""" + KAIBEL COLUMN: Modeling, Optimization, and + Conceptual Design of Multi-product Dividing Wall Columns + (written by E. Soraya Rawlings, esoraya@rwlngs.net, + in collaboration with Qi Chen) + +This is a dividing wall distillation column design problem to determine the optimal minimum number of trays and the optimal location of side streams for the separation of a quaternary mixture: + 1 = methanol + 2 = ethanol + 3 = propanol + 4 = butanol +while minimizing its capital and operating costs. + + The scheme of the Kaibel Column is shown in Figure 1: + ____ + --|Cond|--- + | ---- | + -------------- | + | sect 4 |<------> D mostly 1 + | ---------- | + | | + | ----- |----- | + | | |-------> S2 + Fj ------>| sect | sect | + | 2 | 3 | + | | |-------> S1 + | ----- |----- | + | | + | ---------- | + | sect 1 |<------> B mostly 4 + -------------- | + | ____ | + ---|Reb |--- + ---- + Figure 1. Kaibel Column scheme + +Permanent trays: +- Reboiler and vapor distributor in the bottom section (sect 1) +- Liquid distributor and condenser in the top section (sect 4) +- Side feed tray for the feed side and dividing wall starting and and ening tray in the feed section (sect 2). +- Side product trays and dividing wall starting and ending tray in the product section (sect 3). + +The trays in each section are counted from bottom to top, being tray 1 the bottom tray in each section and tray np the top tray in each section, where np is a specified upper bound for the number of possible trays for each section. +Each section has the same number of possible trays. + +Six degrees of freedom: the reflux ratio, the product outlets (bottom, intermediate, and distillate product flowrates), and the liquid and vapor flowrates between the two sections of the dividing wall, controlled by a liquid and vapor distributor on the top and bottom section of the column, respectively. +including also the vapor and liquid flowrate and the energy consumption in the reboiler and condenser. +The vapor distributor is fixed and remain constant during the column operation. + +Source paper: +Rawlings, E. S., Chen, Q., Grossmann, I. E., & Caballero, J. A. (2019). Kaibel Column: Modeling, optimization, and conceptual design of multi-product dividing wall columns. *Computers and Chemical Engineering*, 125, 31–39. https://doi.org/10.1016/j.compchemeng.2019.03.006 + +""" + +from math import fabs + +import matplotlib.pyplot as plt +from pyomo.environ import * + +# from kaibel_solve_gdp import build_model +from gdplib.kaibel.kaibel_solve_gdp import build_model + + +def main(): + """ + This is the main function that executes the optimization process. + + It builds the model, fixes certain variables, sets initial values for tray existence or absence, + solves the model using the 'gdpopt' solver, and displays the results. + + Returns: + None + """ + m = build_model() + + # Fixing variables + m.F[1].fix(50) # feed flowrate in mol/s + m.F[2].fix(50) + m.F[3].fix(50) + m.F[4].fix(50) + m.q.fix( + m.q_init + ) # vapor fraction q_init from the feed set in the build_model function + m.dv[2].fix(0.394299) # vapor distributor in the feed section + + for sec in m.section: + for n_tray in m.tray: + m.P[sec, n_tray].fix(m.Preb) + + ## Initial values for the tray existence or absence + for n_tray in m.candidate_trays_main: + for sec in m.section_main: + m.tray_exists[sec, n_tray].indicator_var.set_value(1) + m.tray_absent[sec, n_tray].indicator_var.set_value(0) + for n_tray in m.candidate_trays_feed: + m.tray_exists[2, n_tray].indicator_var.set_value(1) + m.tray_absent[2, n_tray].indicator_var.set_value(0) + for n_tray in m.candidate_trays_product: + m.tray_exists[3, n_tray].indicator_var.set_value(1) + m.tray_absent[3, n_tray].indicator_var.set_value(0) + + intro_message(m) + + results = SolverFactory('gdpopt').solve( + m, + strategy='LOA', + tee=True, + time_limit=3600, + mip_solver='gams', + mip_solver_args=dict(solver='cplex'), + ) + + m.calc_nt = sum( + sum(m.tray_exists[sec, n_tray].indicator_var.value for n_tray in m.tray) + for sec in m.section + ) - sum(m.tray_exists[3, n_tray].indicator_var.value for n_tray in m.tray) + m.dw_start = ( + sum(m.tray_exists[1, n_tray].indicator_var.value for n_tray in m.tray) + 1 + ) + m.dw_end = sum( + m.tray_exists[1, n_tray].indicator_var.value for n_tray in m.tray + ) + sum(m.tray_exists[2, n_tray].indicator_var.value for n_tray in m.tray) + + display_results(m) + + print(' ', results) + print(' Solver Status: ', results.solver.status) + print(' Termination condition: ', results.solver.termination_condition) + + +def intro_message(m): + """ + Display the introduction message. + + + + """ + print(""" + +If you use this model and/or initialization strategy, you may cite the following: +Rawlings, ES; Chen, Q; Grossmann, IE; Caballero, JA. Kaibel Column: Modeling, +optimization, and conceptual design of multi-product dividing wall columns. +Comp. and Chem. Eng., 2019, 125, 31-39. +DOI: https://doi.org/10.1016/j.compchemeng.2019.03.006 + + + """) + + +def display_results(m): + """ + Display the results of the optimization process. + + Parameters + ---------- + m : Pyomo ConcreteModel + The Pyomo model object containing the results of the optimization process. + """ + print('') + print('Components:') + print('1 methanol') + print('2 ethanol') + print('3 butanol') + print('4 propanol') + print(' ') + print('Objective: %s' % value(m.obj)) + print('Trays: %s' % value(m.calc_nt)) + print('Dividing_wall_start: %s' % value(m.dw_start)) + print('Dividing_wall_end: %s' % value(m.dw_end)) + print(' ') + print( + 'Qreb: {: >3.0f}kW B_1: {: > 2.0f} B_2: {: >2.0f} B_3: {: >2.0f} B_4: {: >2.0f} Btotal: {: >2.0f}'.format( + value(m.Qreb / m.Qscale), + value(m.B[1]), + value(m.B[2]), + value(m.B[3]), + value(m.B[4]), + value(m.Btotal), + ) + ) + print( + 'Qcon: {: >2.0f}kW D_1: {: >2.0f} D_2: {: >2.0f} D_3: {: >2.0f} D_4: {: >2.0f} Dtotal: {: >2.0f}'.format( + value(m.Qcon / m.Qscale), + value(m.D[1]), + value(m.D[2]), + value(m.D[3]), + value(m.D[4]), + value(m.Dtotal), + ) + ) + print(' ') + print('Reflux: {: >3.4f}'.format(value(m.rr))) + print('Reboil: {: >3.4f} '.format(value(m.bu))) + print(' ') + print('Flowrates[mol/s]') + print( + 'F_1: {: > 3.0f} F_2: {: >2.0f} F_3: {: >2.0f} F_4: {: >2.0f} Ftotal: {: >2.0f}'.format( + value(m.F[1]), + value(m.F[2]), + value(m.F[3]), + value(m.F[4]), + sum(value(m.F[comp]) for comp in m.comp), + ) + ) + print( + 'S1_1: {: > 1.0f} S1_2: {: >2.0f} S1_3: {: >2.0f} S1_4: {: >2.0f} S1total: {: >2.0f}'.format( + value(m.S[1, 1]), + value(m.S[1, 2]), + value(m.S[1, 3]), + value(m.S[1, 4]), + sum(value(m.S[1, comp]) for comp in m.comp), + ) + ) + print( + 'S2_1: {: > 1.0f} S2_2: {: >2.0f} S2_3: {: >2.0f} S2_4: {: >2.0f} S2total: {: >2.0f}'.format( + value(m.S[2, 1]), + value(m.S[2, 2]), + value(m.S[2, 3]), + value(m.S[2, 4]), + sum(value(m.S[2, comp]) for comp in m.comp), + ) + ) + print(' ') + print('Distributors:') + print('dl[2]: {: >3.4f} dl[3]: {: >3.4f}'.format(value(m.dl[2]), value(m.dl[3]))) + print('dv[2]: {: >3.4f} dv[3]: {: >3.4f}'.format(value(m.dv[2]), value(m.dv[3]))) + print(' ') + print(' ') + print(' ') + print(' Kaibel Column Sections ') + print('__________________________________________') + print(' ') + print(' Tray Bottom Feed ') + print('__________________________________________') + for t in reversed(list(m.tray)): + print( + '[{: >2.0f}] {: >9.0g} {: >18.0g} F:{: >3.0f} '.format( + t, + ( + fabs(value(m.tray_exists[1, t].indicator_var)) + if t in m.candidate_trays_main + else 1 + ), + ( + fabs(value(m.tray_exists[2, t].indicator_var)) + if t in m.candidate_trays_feed + else 1 + ), + sum(value(m.F[comp]) for comp in m.comp) if t == m.feed_tray else 0, + ) + ) + print(' ') + print('__________________________________________') + print(' ') + print(' Product Top ') + print('__________________________________________') + for t in reversed(list(m.tray)): + print( + '[{: >2.0f}] {: >9.0g} S1:{: >2.0f} S2:{: >2.0f} {: >8.0g}'.format( + t, + ( + fabs(value(m.tray_exists[3, t].indicator_var)) + if t in m.candidate_trays_product + else 1 + ), + ( + sum(value(m.S[1, comp]) for comp in m.comp) + if t == m.sideout1_tray + else 0 + ), + ( + sum(value(m.S[2, comp]) for comp in m.comp) + if t == m.sideout2_tray + else 0 + ), + ( + fabs(value(m.tray_exists[4, t].indicator_var)) + if t in m.candidate_trays_main + else 1 + ), + ) + ) + print(' 1 = trays exists, 0 = absent tray') + + +if __name__ == "__main__": + main() diff --git a/tests/test_pr47_regressions.py b/tests/test_pr47_regressions.py index fea5d8a5..45855997 100644 --- a/tests/test_pr47_regressions.py +++ b/tests/test_pr47_regressions.py @@ -1,3 +1,5 @@ +from math import exp + import pytest from pyomo.environ import TransformationFactory, value @@ -8,7 +10,7 @@ def test_hda_model_constructs_with_scalar_heat_capacity_ratio(): model = HDA_model() - assert value(model.heat_capacity_ratio) == 1.3 + assert value(model.gam) == 1.3 assert model.gamma.is_indexed() @@ -29,3 +31,62 @@ def test_kaibel_gdp_transformation_constructs(transformation): model = build_kaibel_model() TransformationFactory(transformation).apply_to(model) + + +def _original_vapor_composition_residual(model, section, tray, component): + temperature = value(model.T[section, tray]) + reduced_temperature = value(model.prop[component, 'TC']) / temperature + vapor_pressure_exponent = reduced_temperature * ( + value(model.prop[component, 'vpA']) + * (1 - temperature / value(model.prop[component, 'TC'])) + + value(model.prop[component, 'vpB']) + * (1 - temperature / value(model.prop[component, 'TC'])) ** 1.5 + + value(model.prop[component, 'vpC']) + * (1 - temperature / value(model.prop[component, 'TC'])) ** 3 + + value(model.prop[component, 'vpD']) + * (1 - temperature / value(model.prop[component, 'TC'])) ** 6 + ) + vapor_pressure = value(model.prop[component, 'PC']) * exp(vapor_pressure_exponent) + original_rhs = ( + value(model.x[section, tray, component]) + * value(model.actv[section, tray, component]) + * vapor_pressure + / value(model.P[section, tray]) + ) + return value(model.y[section, tray, component]) - original_rhs + + +@pytest.mark.parametrize( + ('section', 'reduced_temperature_name', 'vapor_composition_name'), + [ + (1, '_bottom_reduced_temperature', 'bottom_vapor_composition'), + (2, '_feedside_reduced_temperature', 'feedside_vapor_composition'), + (3, '_productside_reduced_temperature', 'productside_vapor_composition'), + (4, '_top_reduced_temperature', 'top_vapor_composition'), + ], +) +def test_kaibel_reduced_temperature_reformulation_matches_original_algebra( + section, reduced_temperature_name, vapor_composition_name +): + model = build_kaibel_model() + tray = 1 + + for component in model.comp: + temperature = value(model.T[section, tray]) + model.Tr[section, tray, component].set_value( + value(model.prop[component, 'TC']) / temperature + ) + + disjunct = model.tray_exists[section, tray] + reduced_temperature = getattr(disjunct, reduced_temperature_name)[component] + vapor_composition = getattr(disjunct, vapor_composition_name)[component] + scaled_reformulation_residual = value(vapor_composition.body) / value( + model.P[section, tray] + ) + + assert value(reduced_temperature.body) == pytest.approx( + value(reduced_temperature.lower) + ) + assert scaled_reformulation_residual == pytest.approx( + _original_vapor_composition_residual(model, section, tray, component) + )