From fffd81cab9f6df71830f4699abed74fd1ebfb8fe Mon Sep 17 00:00:00 2001 From: Carolina Tristan Date: Tue, 3 Sep 2024 07:33:00 -0400 Subject: [PATCH 01/21] Add water_network instance --- gdplib/water_network/T.csv | 5 + gdplib/water_network/TU.csv | 6 + gdplib/water_network/__init__ .py | 18 + gdplib/water_network/feed.csv | 6 + gdplib/water_network/wnd.py | 787 ++++++++++++++++++++++++++++++ 5 files changed, 822 insertions(+) create mode 100644 gdplib/water_network/T.csv create mode 100644 gdplib/water_network/TU.csv create mode 100644 gdplib/water_network/__init__ .py create mode 100644 gdplib/water_network/feed.csv create mode 100644 gdplib/water_network/wnd.py diff --git a/gdplib/water_network/T.csv b/gdplib/water_network/T.csv new file mode 100644 index 0000000..55a948f --- /dev/null +++ b/gdplib/water_network/T.csv @@ -0,0 +1,5 @@ +cont,T +A,3 +B,3 +C,3 +D,3 diff --git a/gdplib/water_network/TU.csv b/gdplib/water_network/TU.csv new file mode 100644 index 0000000..a7594ea --- /dev/null +++ b/gdplib/water_network/TU.csv @@ -0,0 +1,6 @@ +TU,A,B,C,D,L,theta,beta +t1,0.95,0.95,0.6,0.0,5,1500,8000 +t2,0.08,0.8,0.6,0.0,3,1000,8000 +t3,0.0,0.6,0.95,0.95,4,4000,8000 +t4,0.0,0.6,0.8,0.85,3,3000,8000 +t5,0.0,0.6,0.85,0.8,1,3000,8000 diff --git a/gdplib/water_network/__init__ .py b/gdplib/water_network/__init__ .py new file mode 100644 index 0000000..f77f67b --- /dev/null +++ b/gdplib/water_network/__init__ .py @@ -0,0 +1,18 @@ +# from .wn_pwl import build_model_Piecewise as _pwl +# from .wn_q import build_model_Quadratic as _q + +# from .wn_minlp import build_model_MINLP as _original +from gdplib.water_network.wnd import build_model + +# def build_model(case="quadratic"): +# if case == "quadratic": +# return _q() +# elif case == "piecewise": +# return _pwl() +# else: +# raise ValueError(f"Invalid case: {case}") + + +__all__ = ['build_model'] +# if __name__ == "__main__": +# build_model diff --git a/gdplib/water_network/feed.csv b/gdplib/water_network/feed.csv new file mode 100644 index 0000000..88c91cc --- /dev/null +++ b/gdplib/water_network/feed.csv @@ -0,0 +1,6 @@ +FSU,A,B,C,D,flow_rate +fs1,1.0,1.0,1.0,1.0,10 +fs2,1.3333333333333333,0.0,0.0,0.6666666666666666,15 +fs3,2.0,2.0,0.0,4.0,5 +fs4,1.5,0.5,0.5,0.0,20 +fs5,2.0,2.0,0.0,0.0,10 diff --git a/gdplib/water_network/wnd.py b/gdplib/water_network/wnd.py new file mode 100644 index 0000000..9761866 --- /dev/null +++ b/gdplib/water_network/wnd.py @@ -0,0 +1,787 @@ +""" +Water Network Design +-------------------- + +Water networks involve water-using process and treatment units, offering several water integration alternatives that reduce freshwater consumption and wastewater generation while minimizing the total network cost subject to a specified discharge limit. + +In the Water Treatment Network (WTN) design problem, given is a set of water streams with known concentrations of contaminants and flow rate. +The objective is to find the set of treatment units and interconnections that minimize the cost of the WTN while satisfying maximum concentrations of contaminants in the reclaimed outlet stream. +The WTN superstructure consists of a set of treatment units, contaminated feed streams carrying a set of contaminants, and a discharge unit. +The fouled feed waters can be allocated to one or more treatment units or disposed of in the sink unit. Upon treatment, the reclaimed streams can be recycled, forwarded to other treatment units, or discharged into the sink unit. + +The mass balances are defined in terms of total flows and contaminants concentration. +Nonconvexities arise from bilinear terms “flows times concentration” in the mixers mass balances and concave investment cost functions of treatment units. + +The instance incorporates two approximations of the concave cost term (piecewise linear and quadratic) to reformulate the GDP model into a bilinear quadratic one. + +The general model description can be summarized as follows: +Min Cost of Treatment Units +s.t. +Physical Constraints: +(a) Mass balance around each splitter +(b) Mass balance around each mixer +(c) Mass balance around each treatment unit +Performance Constraints: +(d) Contaminant composition of the purified stream less or equal than a given limit +for each contaminant. +Logic Constraints: +(e) Treatment units not chosen have their inlet flow set to zero +(f) Every treatment unit chosen must have a minimum flow + +Assumptions: +(i) The performance of the treatment units only depends on the total flow entering the unit and its composition. +(ii) The flow of contaminants leaving the unit is a linear function of the inlet flow of contaminants. + +Case study +---------- +The WTN comprises five inlet streams with four contaminants and four treatment units. +The contaminant concentration and flow rate of the feed streams, contaminant recovery rates, minimum flow rate and cost coefficients of the treatment units, and the upper limit on the molar flow of contaminant j in the purified stream, are reported in (Ruiz and Grossmann, 2009). + +References +---------- +Ruiz J., Grossmann IE. Water Treatment Network Design. 2009 Available from CyberInfrastructure for [MINLP](), a collaboration of Carnegie Mellon University and IBM +at: www.minlp.org/library/problem/index.php?i=24 + +Ruiz, J. P., & Grossmann, I. E. (2011). Using redundancy to strengthen the relaxation for the global optimization of MINLP problems. Computers & Chemical Engineering, 35(12), 2729–2740. https://doi.org/10.1016/J.COMPCHEMENG.2011.01.035 + +""" + +# Importing the required libraries +import pyomo.environ as pyo +import pandas as pd +import numpy as np +from scipy.optimize import curve_fit +from pint import UnitRegistry + +ureg = UnitRegistry() +Q_ = ureg.Quantity + +# import pyomo.gdp as gdp +import os +import re + +# from pyomo.contrib.preprocessing.plugins import init_vars + +from pyomo.core.expr.logical_expr import * + +wnd_dir = os.path.dirname(os.path.realpath(__file__)) + +# Data for the problem +# The feed.csv file contains the contaminant concentration and flow rate of the feed streams. +# The TU.csv file contains the contaminant recovery rates, minimum flow rate, and cost coefficients of the treatment units. +# The T.csv file contains the upper limit on the molar flow of contaminant j in the purified stream. + +feed = pd.read_csv(os.path.join(wnd_dir, "feed.csv"), index_col=0) +TU = pd.read_csv(os.path.join(wnd_dir, "TU.csv"), index_col=0) +T = pd.read_csv(os.path.join(wnd_dir, "T.csv"), index_col=0) + +# GDP model + +# TU mass balances in the disjunct. + + +def build_model(approximation='quadratic'): + """ + Builds a Pyomo ConcreteModel for Water Network Design. + Generates a Pyomo model for the water network design problem with the specified approximation for the capital cost term of the active treatment units. + + Returns + ------- + pyo.ConcreteModel + A Pyomo ConcreteModel object for the water network design. + + """ + m = pyo.ConcreteModel('Water Network Design') + + # ============================================================================ + # Sets + # ============================================================================= + + m.contaminant = pyo.Set(doc='Set of contaminants', initialize=T.index.tolist()) + + m.FSU = pyo.Set(doc='Set of feed splitters', initialize=feed.index.tolist()) + + m.feed = pyo.Set( + doc='Set of feedstreams', + initialize=['f' + str(nf) for nf in pyo.RangeSet(len(m.FSU))], + ) + + m.TU = pyo.Set(doc='Set of treatment units, TU', initialize=TU.index.tolist()) + + m.inTU = pyo.Set( + doc="Inlet TU Port", + initialize=[ + 'mt' + str(ntu) for ntu in pyo.RangeSet(len(m.TU)) + ], # pyo.RangeSet(m.nr)] + ) + + m.outTU = pyo.Set( + doc="Outlet TU Port", + initialize=[ + 'st' + str(ntu) for ntu in pyo.RangeSet(len(m.TU)) + ], # pyo.RangeSet(m.nr)] + ) + + m.TUport = pyo.Set(doc="Inlet and Outlet TU Ports", initialize=m.inTU | m.outTU) + + m.units = pyo.Set( + doc="Superstructure Units", + initialize=m.feed | m.FSU | m.TU | m.TUport | ['dm'] | ['sink'], + ) + + m.splitters = pyo.Set( + doc="Set of splitters", within=m.units, initialize=m.FSU | m.outTU + ) + + m.mixers = pyo.Set(doc="Set of mixers", within=m.units, initialize=['dm'] | m.inTU) + + m.MU_TU_streams = pyo.Set( + doc="MU to TU 1-1 port pairing", + initialize=m.inTU * m.TU, + filter=lambda _, x, y: re.findall(r'\d+', x) == re.findall(r'\d+', y), + ) + + m.TU_SU_streams = pyo.Set( + doc="TU to SU 1-1 port pairing", + initialize=m.TU * m.outTU, + filter=lambda _, x, y: re.findall(r'\d+', x) == re.findall(r'\d+', y), + ) + + m.TU_streams = pyo.Set( + doc="Set of feasible TU streams", + initialize=m.MU_TU_streams | m.outTU * m.inTU | m.TU_SU_streams, + ) + + m.feed_streams = pyo.Set( + doc="Feed to FSU 1-1 port pairing", + initialize=m.feed * m.FSU, + filter=lambda _, x, y: re.findall(r'\d+', x) == re.findall(r'\d+', y), + ) + + m.streams = pyo.Set( + doc="Set of feasible streams", + initialize=m.feed_streams + | m.FSU * ['dm'] + | m.FSU * m.inTU + | m.TU_streams + | m.outTU * ['dm'] + | [('dm', 'sink')], + # initialize=m.feed_streams|m.FSU*['dm']|m.FSU*m.inTU|m.outTU*m.inTU|m.outTU*['dm']|[('dm','sink')] + ) + + m.from_splitters = pyo.Set( + doc='Streams from splitters', + within=m.streams, + initialize=m.streams, + filter=lambda _, x, y: x in m.splitters or (x, y) == ('dm', 'sink'), + ) + m.to_splitters = pyo.Set( + doc='Streams to splitters', + within=m.streams, + initialize=m.streams, + filter=lambda _, x, y: y in m.splitters or (x, y) in m.feed_streams, + ) + + # ============================================================================= + # Treatment Units Parameters + # ============================================================================= + + m.removal_ratio = pyo.Param( + m.contaminant, + m.TU, + doc='Removal ratio for contaminant j in treatment unit t', + within=pyo.NonNegativeReals, + initialize=lambda _, j, k: TU.loc[k, j], + ) + m.theta = pyo.Param( + m.TU, + within=pyo.NonNegativeReals, + doc='Cost parameter theta for unit TU', + initialize=lambda _, k: TU.loc[k, 'theta'], + ) + m.beta = pyo.Param( + m.TU, + within=pyo.NonNegativeReals, + doc='Cost parameter beta for unit TU', + initialize=lambda _, k: TU.loc[k, 'beta'], + ) + m.gamma = pyo.Param( + m.TU, + within=pyo.NonNegativeReals, + doc='Cost parameter gamma for unit TU', + default=0, + initialize=0, + ) + + # ============================================================================= + # Variables + # ============================================================================= + + m.conc = pyo.Var( + m.contaminant, + m.streams - m.MU_TU_streams - m.TU_SU_streams, + doc="Concentration", + domain=pyo.NonNegativeReals, + bounds=(0, 100), + # bounds=lambda _, j, i, k: (feed[j].min(),feed[j].max()), + initialize=lambda _, j, i, k: feed.loc[i, j] if i in m.FSU else None, + ) + + # Concentration of component j in feedstream i + for j, i, k in m.contaminant * (m.FSU * ['dm'] | m.FSU * m.inTU | m.feed_streams): + # print(j,i,k) + if i in m.feed: + m.conc[j, i, k].fix(feed.loc[k, j]) + else: + m.conc[j, i, k].fix(feed.loc[i, j]) + + m.flow = pyo.Var( + m.streams - m.MU_TU_streams - m.TU_SU_streams, + doc="Flowrate", + domain=pyo.NonNegativeReals, + bounds=lambda _, i, k: ( + (None, feed.loc[i, 'flow_rate']) + if i in m.FSU + # else (None,100), + else (0, feed['flow_rate'].sum()) + ), + ) + + # Fix the flow rates of the feed streams + [m.flow[i, k].fix(feed.loc[k, 'flow_rate']) for i, k in m.feed_streams] + # Fix the flow rate from the mixer to the sink + m.flow['dm', 'sink'].fix(feed['flow_rate'].sum()) + + m.costTU = pyo.Var( + m.TU, + domain=pyo.NonNegativeReals, + doc='CTUk cost of TU', + # bounds=lambda _,tu: (0, m.beta[tu] * feed['flow_rate'].sum() + m.gamma[tu] + m.theta[tu] * 41) # Q = optimal flow t1 + # bounds=lambda _,tu: (0, m.beta[tu] * feed['flow_rate'].sum() + m.gamma[tu] + m.theta[tu] * feed['flow_rate'].sum()) + bounds=lambda _, tu: ( + 0, + m.beta[tu] * feed['flow_rate'].sum() + + m.gamma[tu] + + m.theta[tu] * feed['flow_rate'].sum() ** 0.7, + ), + # initialize=lambda _,tu: m.beta[tu] * m.flow_into[mt] + m.gamma[tu] + m.theta[tu] * m.flow_into[mt]**0.7 for mt,t in m.MU_TU_streams if tu==t + ) + + # ============================================================================= + # Mass balances + # ============================================================================= + + # @m.Expression(m.units) + @m.Expression(m.units - m.TU - m.outTU) + def _flow_into(m, option): + """ + Expression for the flow into a unit. + + Parameters + ---------- + m : ConcreteModel + The Pyomo model. + option : str + The unit except the treatment units and the outlet treatment units. + + Returns + ------- + Expression + The flow into the unit, which is the sum of the flow from the feed splitters and the flow from the mixers. + """ + return sum(m.flow[src, sink] for src, sink in m.streams if sink == option) + + # @m.Expression(m.units, m.contaminant) + @m.Expression(m.units - m.TU - m.outTU, m.contaminant) + def _conc_into(m, option, j): + """ + Expression for the mass balance of the contaminant in the unit to compute the concentration into the unit. + + Parameters + ---------- + m : ConcreteModel + The Pyomo model. + option : str + The units except the treatment units and the treatment units' outlet. + j : str + The contaminant. + + Returns + ------- + Expression + The concentration into the unit, which is the sum of the flow from the feed splitters and the flow from the mixers multiplied by the concentration of the contaminant. + """ + if option in m.mixers: + return sum( + m.flow[src, sink] * m.conc[j, src, sink] + for src, sink in m.streams + if sink == option + ) + return pyo.Expression.Skip + + # @m.Expression(m.units) + @m.Expression(m.units - m.TU - m.inTU) + def _flow_out_from(m, option): + """ + Expression for to compute the flow rate leaving the unit. + + Parameters + ---------- + m : ConcreteModel + The Pyomo model. + option : str + The units except the treatment units and the treatment units' inlet. + + Returns + ------- + Expression + The flow rate leaving the unit, which is the sum of the flow to the discharge mixers and the splitters. + """ + return sum(m.flow[src, sink] for src, sink in m.streams if src == option) + + # @m.Expression(m.units, m.contaminant) + @m.Expression(m.units - m.inTU, m.contaminant) + def _conc_out_from(m, option, j): + """The mass balance of the contaminant in the unit to compute the concentration out of the unit. + + Parameters + ---------- + m : ConcreteModel + The Pyomo model. + option : str + The unit except the treatment units and the treatment units' inlet. + j : str + The contaminant. + + Returns + ------- + Expression + The concentration out of the unit, which is the sum of the flow from the feed splitters and the flow from the mixers multiplied by the concentration of the contaminant. + """ + if option in m.mixers: + return sum( + m.flow[src, sink] * m.conc[j, src, sink] + for src, sink in m.streams + if src == option + ) + return pyo.Expression.Skip + + # Global constraints + + # Adding the mass balances for the splitters and mixers + m.mixer_balances = pyo.ConstraintList() + # for mu in m.mixers: + for mu in m.mixers - m.TU - m.inTU: + # Flowrate balance + [m.mixer_balances.add(m._flow_into[mu] == m._flow_out_from[mu])] + # Concentration balance + [ + m.mixer_balances.add(m._conc_into[mu, j] == m._conc_out_from[mu, j]) + for j in m.contaminant + ] + # m.mixer_balances_global.add(sum(m._conc_into[mu,sol] for sol in m.SOL) == sum(m._conc_out_from[mu,sol] for sol in m.SOL)) + + m.splitter_balances = pyo.ConstraintList() + + for su in m.splitters - m.TU - m.outTU: + # Flowrate balance + [m.splitter_balances.add(m._flow_into[su] == m._flow_out_from[su])] + + for src1, sink1 in m.from_splitters - m.TU_SU_streams: + for src2, sink2 in m.to_splitters - m.TU_SU_streams: + # Concentration balance + [ + m.splitter_balances.add( + m.conc[j, src2, sink2] == m.conc[j, src1, sink1] + ) + for j in m.contaminant + if src1 == sink2 + ] + + @m.Constraint(m.contaminant) + def outlet_contaminant_limit(m, j): + """Constraint: Outlet contaminant limit. + The constraint ensures that the concentration of the contaminant in the outlet stream is less than or equal to the specified limit. + + Parameters + ---------- + m : ConcreteModel + The Pyomo model. + j : str + The contaminant. + + Returns + ------- + Constraint + The constraint that the concentration of the contaminant in the outlet stream is less than or equal to the specified limit. + """ + return m._conc_into['dm', j] <= T.loc[j].values[0] * 10 + + @m.Disjunct(m.TU) + def unit_exists(disj, unit): + '''Disjunct: Unit exists. + This disjunct specifies that the unit exists. + ''' + pass + + @m.Disjunct(m.TU) + def unit_absent(no_unit, unit): + '''Disjunct: Unit absent. + This disjunct specifies that the unit is absent.''' + + @no_unit.Constraint(m.inTU) + def _no_flow_in(disj, mt): + '''Constraint: No flow enters the unit when the unit is inactive. + + This constraint ensures that no flow enters the unit when the unit is inactive. + Flow from the feed splitters and mixers to the unit is set to zero when the unit is inactive. + + Parameters + ---------- + disj : Disjunct + The disjunct for the inactive unit. + mt : str + The inlet TU port. + + Returns + ------- + Constraint + The constraint that no flow enters the unit when the unit is inactive. + ''' + if (mt, unit) in m.MU_TU_streams: + return m._flow_into[mt] == 0 + return pyo.Constraint.Skip + + @no_unit.Constraint(doc='Cost inactive TU') + def _no_cost(disj): + '''Constraint: Cost of inactive unit. + This constraint ensures that the cost of the inactive unit is zero. + + Returns + ------- + Constraint + The constraint enforcing the cost of the inactive unit to be zero. + ''' + return m.costTU[unit] == 0 + + pass + + @m.Disjunction(m.TU) + def unit_exists_or_not(m, unit): + '''Disjunction: Unit exists or not. + This disjunctiont specifies if the treatment unit exists or does not exist. + + Parameters + ---------- + m : ConcreteModel + The Pyomo model. + unit : str + The treatment unit. + + Returns + ------- + Disjunction + The disjunction specifying if the treatment unit exists or does not exist. + ''' + return [m.unit_exists[unit], m.unit_absent[unit]] + + # Yunit is a Boolean variable that indicates if the unit is active. + m.Yunit = pyo.BooleanVar( + m.TU, doc="Boolean variable for existence of a treatment unit" + ) + for unit in m.TU: + m.Yunit[unit].associate_binary_var( + m.unit_exists[unit].indicator_var.get_associated_binary() + ) + + # Logical constraints that ensure that at least one unit is active + # m.atleast_oneRU = pyo.LogicalConstraint(expr=atleast(1, m.Yunit)) + + for unit in m.TU: + ue = m.unit_exists[unit] + ue.streams = pyo.Set( + doc="Streams in active TU", + initialize=m.TU_streams, + filter=lambda _, x, y: x == unit or y == unit, + ) + ue.MU_TU_streams = pyo.Set( + # doc='Mixer to treatment unit streams', + doc="MU to TU 1-1 port pairing", + initialize=m.inTU * m.TU, + # within=ue.streams, + filter=lambda _, x, y: re.findall(r'\d+', x) == re.findall(r'\d+', y) + and y == unit, + ) + + ue.flow = pyo.Var( + ue.streams, + doc="TU streams flowrate", + domain=pyo.NonNegativeReals, + # bounds=lambda _,i,k:(TU.loc[unit,'L'],100) + bounds=lambda _, i, k: (TU.loc[unit, 'L'], feed['flow_rate'].sum()), + # initialize=lambda _,i,k: + ) + + ue.conc = pyo.Var( + m.contaminant, + ue.streams, + doc="TU streams concentration", + domain=pyo.NonNegativeReals, + bounds=lambda _, j, i, k: (0, 100) if i == unit else (0, 4), + # bounds=lambda _, j, i, k: (feed[j].min(),100) if i==unit else (feed[j].min(),4) + # bounds=lambda _, j, i, k: (feed[j].min(),feed[j].max()), + # initialize= lambda _, j, i, k: feed.loc[i,j] if i in m.FSU else None + ) + + # Adding the mass balances for the active treatment units + ue.balances_con = pyo.ConstraintList(doc='TU Material balances') + # The flowrate at the inlet of the treatment unit is equal to the flowrate at the outlet of the treatment unit. + [ + ue.balances_con.add(ue.flow[mt, unit] == ue.flow[unit, st]) + for mt, t in ue.streams + if t == unit + for t, st in ue.streams + if t == unit + ] + # The concentration of the contaminant at the inlet of the treatment unit is equal to the concentration of the contaminant at the outlet of the treatment unit times the removal ratio. + [ + ue.balances_con.add( + ue.conc[j, unit, st] + == (1 - m.removal_ratio[j, t]) * ue.conc[j, mt, unit] + ) + for mt, t in ue.streams + if t == unit + for t, st in ue.streams + if t == unit + for j in m.contaminant + ] + # Treatment unit's mixer mass balance on the flowrate. + [ + ue.balances_con.add(m._flow_into[mt] == ue.flow[mt, unit]) + for mt, t in ue.streams + if t == unit + ] + # Treatment unit's mixer mass balance on the concentration of contaminants. + [ + ue.balances_con.add( + m._conc_into[mt, j] == ue.conc[j, mt, unit] * ue.flow[mt, unit] + ) + for mt, t in ue.streams + if t == unit + for j in m.contaminant + ] + # Treatment unit's splitter mass balance on the flowrate. + [ + ue.balances_con.add(ue.flow[unit, st] == m._flow_out_from[st]) + for t, st in ue.streams + if t == unit + ] + # Treatment unit's splitter mass balance on the concentration of contaminants. + [ + ue.balances_con.add(ue.conc[j, src2, sink2] == m.conc[j, src1, sink1]) + for src1, sink1 in m.from_splitters + for src2, sink2 in ue.streams + if src2 == unit and src1 == sink2 + for j in m.contaminant + ] + # Setting inlet flowrate bounds for the active treatment units. + ue.flow_bound = pyo.ConstraintList(doc='Flowrate bounds to/from active RU') + [ + ue.flow_bound.add( + (ue.flow[mt, unit].lb, m._flow_into[mt], ue.flow[mt, unit].ub) + ) + for mt, t in m.MU_TU_streams + if t == unit + ] + + if approximation == 'quadratic': + # New variable for potential term in capital cost. Z = sum(Q)**0.7, Z >= 0 + ue.cost_var = pyo.Var( + m.TU, + domain=pyo.NonNegativeReals, + initialize=lambda _, unit: TU.loc[unit, 'L'] ** 0.7, + bounds=(0, 100**0.7), + doc="New var for potential term in capital cost", + ) + + def _quadratic_curve_fit(lb, ub, x): + """This function fits a quadratic curve to the function Z^(1/0.7). + + Parameters + ---------- + lb : float + The lower bound of the flow rate. + ub : float + The upper bound of the flow rate. + x : variable + The flow rate. + + Returns + ------- + Expression + The quadratic curve fit to the function Z^(1/0.7). + """ + # Equally spaced points between the lower and upper bounds of the flow rate + z = np.linspace(lb, ub, 100) + + def _func(x, a, b): + """This function computes the quadratic curve fit. + The curve fit is a quadratic function of the form a*x + b*x^2. + + Parameters + ---------- + x : variable + The flow rate. + a : float + The coefficient a. + b : float + The coefficient b. + + Returns + ------- + Expression + The quadratic curve fit for the function Z^(1/0.7). + """ + return a * x + b * x**2 + + # popt is the optimal values for the coefficients a and b, pcov is the covariance matrix + # curve_fit is used to fit the quadratic curve to the function Z^(1/0.7) given the flow rate and the quadratic expression _func. + popt, pcov = curve_fit(_func, z, z ** (1 / 0.7)) + + return popt[0] * x + popt[1] * x**2 + + # Z^(1/0.7)=sum(Q) + @ue.Constraint(doc='New var potential term in capital cost cstr.') + def _cost_nv(ue): + """Constraint: New variable potential term in capital cost. + The constraint ensures that the new variable potential term in the capital cost is equal to the flow rate entering the treatment unit. + + Parameters + ---------- + ue : Disjunct + The disjunct for the active treatment unit. + + Returns + ------- + Constraint + The constraint that the new variable potential term in the capital cost is equal to the flow rate. + """ + for mt, t in ue.streams: + if t == unit: + return ( + _quadratic_curve_fit( + 0, ue.cost_var[unit].ub, ue.cost_var[unit] + ) + == ue.flow[mt, unit] + ) + + elif approximation == "piecewise": + # New variable for potential term in capital cost. Z = sum(Q)**0.7, Z >= 0 + ue.cost_var = pyo.Var( + ue.MU_TU_streams, + domain=pyo.NonNegativeReals, + initialize=lambda _, mt, unit: TU.loc[unit, 'L'] ** 0.7, + bounds=lambda _, mt, unit: ( + TU.loc[unit, 'L'] ** 0.7, + feed['flow_rate'].sum() ** 0.7, + ), + # bounds= lambda _,mt,unit: (TU.loc[unit,'L']**0.7,100**0.7), + doc="New var for potential term in capital cost", + ) + + # to avoid warnings, we set breakpoints at or beyond the bounds + PieceCnt = 100 + bpts = [] + for mt, t in ue.streams: + if t == unit: + Topx = ue.flow[mt, unit].ub + for i in range(PieceCnt + 2): + bpts.append(float((i * Topx) / PieceCnt)) + + # def _breakpoints(lb, ub): + # x = np.linspace(lb, ub, 100) + # my_pwlf = pwlf.PiecewiseLinFit(x, x**0.7, degree=1) + # # my_pwlf.fit(8) + # my_pwlf.fitfast(10, pop=100) + # return my_pwlf.fit_breaks.tolist() + + def _func(model, i, j, xp): + """This function provides the expression for the piecewise linear approximation of the capital cost using the Picewise class. + + Parameters + ---------- + model : ConcreteModel + The Pyomo model. + i : str + j : str + xp : variable + The flow rate. + + Returns + ------- + Expression + The expression for the piecewise linear approximation of the capital cost. + """ + # we not need i, j, but it are passed as the index for the constraint + return xp**0.7 + + # Piecewise is a class that provides a piecewise linear approximation of a function using the INC representation. + # The INC representation is a piecewise linear approximation of a function using the incremental form. + ue.ComputeObj = pyo.Piecewise( + ue.MU_TU_streams, + ue.cost_var, + ue.flow, + pw_pts=bpts, + pw_constr_type='EQ', + f_rule=_func, + # ue.MU_TU_streams, ue.cost_var, ue.flow, pw_pts=_breakpoints(0,Topx), pw_constr_type='EQ', f_rule=func, + # pw_repn='BIGM_BIN' + pw_repn='INC', + ) + + # for i,j in ue.MU_TU_streams: + # ue.ComputeObj[i,j].SOS2_y.setub(1.) + + # Setting bounds for the INC_delta variable which is a binary variable that indicates the piecewise linear segment. + for i, j in ue.MU_TU_streams: + ue.ComputeObj[i, j].INC_delta.setub(1) + ue.ComputeObj[i, j].INC_delta.setlb(0) + + @ue.Constraint(doc='Cost active TU') + def costTU(ue): + """Constraint: Cost of active treatment unit. + The constraint ensures that the cost of the active treatment unit is equal to the sum of an investment cost which is proportional to the total flow to 0.7 exponent and an operating cost which is proportional to the flow + + Parameters + ---------- + ue : Disjunct + The disjunct for the active treatment unit. + + Returns + ------- + Constraint + The constraint that the cost of the active treatment unit is equal to the sum of an investment cost and an operating cost. + """ + for mt, t in ue.streams: + if approximation == 'quadratic': + new_var = ue.cost_var[unit] + elif approximation == 'piecewise': + new_var = ue.cost_var[mt, unit] + if t == unit: + return ( + m.costTU[unit] + == m.beta[unit] * ue.flow[mt, unit] + + m.gamma[unit] + + m.theta[unit] * new_var + ) + + # Objective function: minimize the total cost of the treatment units + m.obj = pyo.Objective(expr=sum(m.costTU[k] for k in m.TU), sense=pyo.minimize) + + # init_vars.InitMidpoint().apply_to(m) + + return m + + +# if __name__ == "__main__": +# build_model() From f784ca0459ac8ea80906f8f9bb1dfced480389c9 Mon Sep 17 00:00:00 2001 From: Carolina Tristan Date: Tue, 3 Sep 2024 15:44:17 -0400 Subject: [PATCH 02/21] Upload README.md --- gdplib/README.md | 62 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 62 insertions(+) create mode 100644 gdplib/README.md diff --git a/gdplib/README.md b/gdplib/README.md new file mode 100644 index 0000000..8cfa313 --- /dev/null +++ b/gdplib/README.md @@ -0,0 +1,62 @@ +## Water Network Design + +In the Water Treatment Network (WTN) design problem, given is a set of water streams with known concentrations of contaminants and flow rate. +The objective is to find the set of treatment units and interconnections that minimize the cost of the WTN while satisfying maximum concentrations of contaminants in the reclaimed outlet stream. +The WTN superstructure consists of a set of treatment units, contaminated feed streams carrying a set of contaminants, and a discharge unit. +The fouled feed waters can be allocated to one or more treatment units or disposed of in the sink unit. Upon treatment, the reclaimed streams can be recycled, forwarded to other treatment units, or discharged into the sink unit. + +The mass balances are defined in terms of total flows and contaminants concentration. +Nonconvexities arise from bilinear terms “flows times concentration” in the mixers mass balances and concave investment cost functions of treatment units. + +The instance incorporates two approximations of the concave cost term (piecewise linear and quadratic) to reformulate the GDP model into a bilinear quadratic one. + +The general model description can be summarized as follows: +Min Cost of Treatment Units +s.t. +Physical Constraints: +(a) Mass balance around each splitter +(b) Mass balance around each mixer +(c) Mass balance around each treatment unit +Performance Constraints: +(d) Contaminant composition of the purified stream less or equal than a given limit +for each contaminant. +Logic Constraints: +(e) Treatment units not chosen have their inlet flow set to zero +(f) Every treatment unit chosen must have a minimum flow + +Assumptions: +(i) The performance of the treatment units only depends on the total flow entering the unit and its composition. +(ii) The flow of contaminants leaving the unit is a linear function of the inlet flow of contaminants. + +### Case Study + +The WTN comprises five inlet streams with four contaminants and four treatment units. +The contaminant concentration and flow rate of the feed streams, contaminant recovery rates, minimum flow rate and cost coefficients of the treatment units, and the upper limit on the molar flow of contaminant in the purified stream, are reported in [1]. + +### Solution + +Best known objective value: 348,340 + +### Size + +Number of reactors in series is 5. + +| Problem | vars | Bool | bin | int | cont | cons | nl | disj | disjtn | +| ----------------- | ---- | ---- | --- | --- | ---- | ---- | -- | ---- | ------ | +| `water_network` | | | | | | | | | | + +- ``vars``: variables +- ``Bool``: Boolean variables +- ``bin``: binary variables +- ``int``: integer variables +- ``cont``: continuous variables +- ``cons``: constraints +- ``nl``: nonlinear constraints +- ``disj``: disjuncts +- ``disjtn``: disjunctions + +### References + +> [1] Ruiz J., Grossmann I. E. Water Treatment Network Design. 2009 Available from CyberInfrastructure for [MINLP](www.minlp.org), a collaboration of Carnegie Mellon University and IBM at: www.minlp.org/library/problem/index.php?i=24 +> +> [2] Ruiz, J., & Grossmann, I. E. (2011). Using redundancy to strengthen the relaxation for the global optimization of MINLP problems. Computers & Chemical Engineering, 35(12), 2729–2740. https://doi.org/10.1016/J.COMPCHEMENG.2011.01.035 From 8b5f14aa341b9779094771a1e760001fce4f698d Mon Sep 17 00:00:00 2001 From: Carolina Tristan Date: Wed, 4 Sep 2024 07:23:35 -0400 Subject: [PATCH 03/21] Remove water_network module and update dependencies --- benchmark.py | 1 + gdplib/__init__.py | 1 + gdplib/water_network/__init__ .py | 18 ------------------ gdplib/water_network/__init__.py | 3 +++ generate_model_size_report.py | 1 + 5 files changed, 6 insertions(+), 18 deletions(-) delete mode 100644 gdplib/water_network/__init__ .py create mode 100644 gdplib/water_network/__init__.py diff --git a/benchmark.py b/benchmark.py index 1bb9df8..5f02dfc 100644 --- a/benchmark.py +++ b/benchmark.py @@ -86,6 +86,7 @@ def benchmark(model, strategy, timelimit, result_dir, subsolver="scip"): # "modprodnet", # "stranded_gas", # "syngas", + # "water_network", ] strategy_list = [ "gdp.bigm", diff --git a/gdplib/__init__.py b/gdplib/__init__.py index a3cb501..dca67c2 100644 --- a/gdplib/__init__.py +++ b/gdplib/__init__.py @@ -12,3 +12,4 @@ import gdplib.disease_model import gdplib.med_term_purchasing import gdplib.syngas +import gdplib.water_network diff --git a/gdplib/water_network/__init__ .py b/gdplib/water_network/__init__ .py deleted file mode 100644 index f77f67b..0000000 --- a/gdplib/water_network/__init__ .py +++ /dev/null @@ -1,18 +0,0 @@ -# from .wn_pwl import build_model_Piecewise as _pwl -# from .wn_q import build_model_Quadratic as _q - -# from .wn_minlp import build_model_MINLP as _original -from gdplib.water_network.wnd import build_model - -# def build_model(case="quadratic"): -# if case == "quadratic": -# return _q() -# elif case == "piecewise": -# return _pwl() -# else: -# raise ValueError(f"Invalid case: {case}") - - -__all__ = ['build_model'] -# if __name__ == "__main__": -# build_model diff --git a/gdplib/water_network/__init__.py b/gdplib/water_network/__init__.py new file mode 100644 index 0000000..d3ac558 --- /dev/null +++ b/gdplib/water_network/__init__.py @@ -0,0 +1,3 @@ +from .wnd import build_model + +__all__ = ['build_model'] diff --git a/generate_model_size_report.py b/generate_model_size_report.py index 96ddac0..5995147 100644 --- a/generate_model_size_report.py +++ b/generate_model_size_report.py @@ -20,6 +20,7 @@ # "modprodnet", # "stranded_gas", # "syngas", + # "water_network" ] current_time = datetime.now().strftime("%Y-%m-%d_%H-%M-%S") timelimit = 600 From 435b80cfac6ada3a9419a6e26f69ea8a5858ba92 Mon Sep 17 00:00:00 2001 From: Carolina Tristan Date: Mon, 9 Sep 2024 20:58:32 -0400 Subject: [PATCH 04/21] Add `ue` in typos.toml to prevent typo flags. --- .github/workflows/typos.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/workflows/typos.toml b/.github/workflows/typos.toml index a7e28f0..53a49eb 100644 --- a/.github/workflows/typos.toml +++ b/.github/workflows/typos.toml @@ -3,3 +3,4 @@ hda = "hda" HDA = "HDA" equil = "equil" +ue = "ue" From b69965c9e89aa0dcc7a4919b317adda7f65a80a4 Mon Sep 17 00:00:00 2001 From: Carolina Tristan Date: Mon, 9 Sep 2024 21:02:45 -0400 Subject: [PATCH 05/21] Update README.md with best known objective value and component sizes --- gdplib/README.md | 27 +++++++++++---------------- 1 file changed, 11 insertions(+), 16 deletions(-) diff --git a/gdplib/README.md b/gdplib/README.md index 8cfa313..2be0b37 100644 --- a/gdplib/README.md +++ b/gdplib/README.md @@ -35,25 +35,20 @@ The contaminant concentration and flow rate of the feed streams, contaminant rec ### Solution -Best known objective value: 348,340 +Best known objective value: $ 348,340 ### Size -Number of reactors in series is 5. - -| Problem | vars | Bool | bin | int | cont | cons | nl | disj | disjtn | -| ----------------- | ---- | ---- | --- | --- | ---- | ---- | -- | ---- | ------ | -| `water_network` | | | | | | | | | | - -- ``vars``: variables -- ``Bool``: Boolean variables -- ``bin``: binary variables -- ``int``: integer variables -- ``cont``: continuous variables -- ``cons``: constraints -- ``nl``: nonlinear constraints -- ``disj``: disjuncts -- ``disjtn``: disjunctions +| Component | pwl | quadratic | +| :-------------------- | :--: | :-------: | +| variables | 1405 | 420 | +| binary_variables | 510 | 10 | +| integer_variables | 0 | 0 | +| continuous_variables | 895 | 410 | +| disjunctions | 5 | 5 | +| disjuncts | 10 | 10 | +| constraints | 1339 | 334 | +| nonlinear_constraints | 28 | 33 | ### References From 6ea4d16d40928f79c63c4f7ec88fdd03c92d2c1f Mon Sep 17 00:00:00 2001 From: Carolina Tristan Date: Mon, 9 Sep 2024 21:16:20 -0400 Subject: [PATCH 06/21] Shifted `README.md` to the `water_network` directory --- gdplib/water_network/README.md | 57 ++++++++++++++++++++++++++++++++++ 1 file changed, 57 insertions(+) create mode 100644 gdplib/water_network/README.md diff --git a/gdplib/water_network/README.md b/gdplib/water_network/README.md new file mode 100644 index 0000000..2be0b37 --- /dev/null +++ b/gdplib/water_network/README.md @@ -0,0 +1,57 @@ +## Water Network Design + +In the Water Treatment Network (WTN) design problem, given is a set of water streams with known concentrations of contaminants and flow rate. +The objective is to find the set of treatment units and interconnections that minimize the cost of the WTN while satisfying maximum concentrations of contaminants in the reclaimed outlet stream. +The WTN superstructure consists of a set of treatment units, contaminated feed streams carrying a set of contaminants, and a discharge unit. +The fouled feed waters can be allocated to one or more treatment units or disposed of in the sink unit. Upon treatment, the reclaimed streams can be recycled, forwarded to other treatment units, or discharged into the sink unit. + +The mass balances are defined in terms of total flows and contaminants concentration. +Nonconvexities arise from bilinear terms “flows times concentration” in the mixers mass balances and concave investment cost functions of treatment units. + +The instance incorporates two approximations of the concave cost term (piecewise linear and quadratic) to reformulate the GDP model into a bilinear quadratic one. + +The general model description can be summarized as follows: +Min Cost of Treatment Units +s.t. +Physical Constraints: +(a) Mass balance around each splitter +(b) Mass balance around each mixer +(c) Mass balance around each treatment unit +Performance Constraints: +(d) Contaminant composition of the purified stream less or equal than a given limit +for each contaminant. +Logic Constraints: +(e) Treatment units not chosen have their inlet flow set to zero +(f) Every treatment unit chosen must have a minimum flow + +Assumptions: +(i) The performance of the treatment units only depends on the total flow entering the unit and its composition. +(ii) The flow of contaminants leaving the unit is a linear function of the inlet flow of contaminants. + +### Case Study + +The WTN comprises five inlet streams with four contaminants and four treatment units. +The contaminant concentration and flow rate of the feed streams, contaminant recovery rates, minimum flow rate and cost coefficients of the treatment units, and the upper limit on the molar flow of contaminant in the purified stream, are reported in [1]. + +### Solution + +Best known objective value: $ 348,340 + +### Size + +| Component | pwl | quadratic | +| :-------------------- | :--: | :-------: | +| variables | 1405 | 420 | +| binary_variables | 510 | 10 | +| integer_variables | 0 | 0 | +| continuous_variables | 895 | 410 | +| disjunctions | 5 | 5 | +| disjuncts | 10 | 10 | +| constraints | 1339 | 334 | +| nonlinear_constraints | 28 | 33 | + +### References + +> [1] Ruiz J., Grossmann I. E. Water Treatment Network Design. 2009 Available from CyberInfrastructure for [MINLP](www.minlp.org), a collaboration of Carnegie Mellon University and IBM at: www.minlp.org/library/problem/index.php?i=24 +> +> [2] Ruiz, J., & Grossmann, I. E. (2011). Using redundancy to strengthen the relaxation for the global optimization of MINLP problems. Computers & Chemical Engineering, 35(12), 2729–2740. https://doi.org/10.1016/J.COMPCHEMENG.2011.01.035 From eae3d937884afba408143d0b0ee3cce07a74464f Mon Sep 17 00:00:00 2001 From: Carolina Tristan Date: Mon, 9 Sep 2024 21:18:58 -0400 Subject: [PATCH 07/21] Delet `README.md` from `gdplib` directory. --- gdplib/README.md | 57 ------------------------------------------------ 1 file changed, 57 deletions(-) delete mode 100644 gdplib/README.md diff --git a/gdplib/README.md b/gdplib/README.md deleted file mode 100644 index 2be0b37..0000000 --- a/gdplib/README.md +++ /dev/null @@ -1,57 +0,0 @@ -## Water Network Design - -In the Water Treatment Network (WTN) design problem, given is a set of water streams with known concentrations of contaminants and flow rate. -The objective is to find the set of treatment units and interconnections that minimize the cost of the WTN while satisfying maximum concentrations of contaminants in the reclaimed outlet stream. -The WTN superstructure consists of a set of treatment units, contaminated feed streams carrying a set of contaminants, and a discharge unit. -The fouled feed waters can be allocated to one or more treatment units or disposed of in the sink unit. Upon treatment, the reclaimed streams can be recycled, forwarded to other treatment units, or discharged into the sink unit. - -The mass balances are defined in terms of total flows and contaminants concentration. -Nonconvexities arise from bilinear terms “flows times concentration” in the mixers mass balances and concave investment cost functions of treatment units. - -The instance incorporates two approximations of the concave cost term (piecewise linear and quadratic) to reformulate the GDP model into a bilinear quadratic one. - -The general model description can be summarized as follows: -Min Cost of Treatment Units -s.t. -Physical Constraints: -(a) Mass balance around each splitter -(b) Mass balance around each mixer -(c) Mass balance around each treatment unit -Performance Constraints: -(d) Contaminant composition of the purified stream less or equal than a given limit -for each contaminant. -Logic Constraints: -(e) Treatment units not chosen have their inlet flow set to zero -(f) Every treatment unit chosen must have a minimum flow - -Assumptions: -(i) The performance of the treatment units only depends on the total flow entering the unit and its composition. -(ii) The flow of contaminants leaving the unit is a linear function of the inlet flow of contaminants. - -### Case Study - -The WTN comprises five inlet streams with four contaminants and four treatment units. -The contaminant concentration and flow rate of the feed streams, contaminant recovery rates, minimum flow rate and cost coefficients of the treatment units, and the upper limit on the molar flow of contaminant in the purified stream, are reported in [1]. - -### Solution - -Best known objective value: $ 348,340 - -### Size - -| Component | pwl | quadratic | -| :-------------------- | :--: | :-------: | -| variables | 1405 | 420 | -| binary_variables | 510 | 10 | -| integer_variables | 0 | 0 | -| continuous_variables | 895 | 410 | -| disjunctions | 5 | 5 | -| disjuncts | 10 | 10 | -| constraints | 1339 | 334 | -| nonlinear_constraints | 28 | 33 | - -### References - -> [1] Ruiz J., Grossmann I. E. Water Treatment Network Design. 2009 Available from CyberInfrastructure for [MINLP](www.minlp.org), a collaboration of Carnegie Mellon University and IBM at: www.minlp.org/library/problem/index.php?i=24 -> -> [2] Ruiz, J., & Grossmann, I. E. (2011). Using redundancy to strengthen the relaxation for the global optimization of MINLP problems. Computers & Chemical Engineering, 35(12), 2729–2740. https://doi.org/10.1016/J.COMPCHEMENG.2011.01.035 From ac2397beb228f83ec4029b409dc0cc7a2875d22b Mon Sep 17 00:00:00 2001 From: Carolina Tristan Date: Mon, 9 Sep 2024 21:22:43 -0400 Subject: [PATCH 08/21] Remove unused import statement in wnd.py --- gdplib/water_network/wnd.py | 1 - 1 file changed, 1 deletion(-) diff --git a/gdplib/water_network/wnd.py b/gdplib/water_network/wnd.py index 9761866..b1c286b 100644 --- a/gdplib/water_network/wnd.py +++ b/gdplib/water_network/wnd.py @@ -56,7 +56,6 @@ ureg = UnitRegistry() Q_ = ureg.Quantity -# import pyomo.gdp as gdp import os import re From d6b9cba56798ba8beb6f39596a1e5bab380d88e1 Mon Sep 17 00:00:00 2001 From: Carolina Tristan Date: Mon, 9 Sep 2024 21:24:57 -0400 Subject: [PATCH 09/21] Refactor import statements in wnd.py --- gdplib/water_network/wnd.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/gdplib/water_network/wnd.py b/gdplib/water_network/wnd.py index b1c286b..70c75a5 100644 --- a/gdplib/water_network/wnd.py +++ b/gdplib/water_network/wnd.py @@ -50,14 +50,14 @@ import pyomo.environ as pyo import pandas as pd import numpy as np +import os +import re from scipy.optimize import curve_fit from pint import UnitRegistry ureg = UnitRegistry() Q_ = ureg.Quantity -import os -import re # from pyomo.contrib.preprocessing.plugins import init_vars From 88c7b6d5dbb464242cdd4515be698fdfd6daa223 Mon Sep 17 00:00:00 2001 From: Carolina Tristan Date: Mon, 9 Sep 2024 21:46:59 -0400 Subject: [PATCH 10/21] Remove comments in `wnd.py` --- gdplib/water_network/wnd.py | 45 +++++-------------------------------- 1 file changed, 5 insertions(+), 40 deletions(-) diff --git a/gdplib/water_network/wnd.py b/gdplib/water_network/wnd.py index 70c75a5..ac05460 100644 --- a/gdplib/water_network/wnd.py +++ b/gdplib/water_network/wnd.py @@ -58,9 +58,6 @@ ureg = UnitRegistry() Q_ = ureg.Quantity - -# from pyomo.contrib.preprocessing.plugins import init_vars - from pyomo.core.expr.logical_expr import * wnd_dir = os.path.dirname(os.path.realpath(__file__)) @@ -109,16 +106,12 @@ def build_model(approximation='quadratic'): m.inTU = pyo.Set( doc="Inlet TU Port", - initialize=[ - 'mt' + str(ntu) for ntu in pyo.RangeSet(len(m.TU)) - ], # pyo.RangeSet(m.nr)] + initialize=['mt' + str(ntu) for ntu in pyo.RangeSet(len(m.TU))], ) m.outTU = pyo.Set( doc="Outlet TU Port", - initialize=[ - 'st' + str(ntu) for ntu in pyo.RangeSet(len(m.TU)) - ], # pyo.RangeSet(m.nr)] + initialize=['st' + str(ntu) for ntu in pyo.RangeSet(len(m.TU))], ) m.TUport = pyo.Set(doc="Inlet and Outlet TU Ports", initialize=m.inTU | m.outTU) @@ -165,7 +158,6 @@ def build_model(approximation='quadratic'): | m.TU_streams | m.outTU * ['dm'] | [('dm', 'sink')], - # initialize=m.feed_streams|m.FSU*['dm']|m.FSU*m.inTU|m.outTU*m.inTU|m.outTU*['dm']|[('dm','sink')] ) m.from_splitters = pyo.Set( @@ -222,7 +214,6 @@ def build_model(approximation='quadratic'): doc="Concentration", domain=pyo.NonNegativeReals, bounds=(0, 100), - # bounds=lambda _, j, i, k: (feed[j].min(),feed[j].max()), initialize=lambda _, j, i, k: feed.loc[i, j] if i in m.FSU else None, ) @@ -241,7 +232,7 @@ def build_model(approximation='quadratic'): bounds=lambda _, i, k: ( (None, feed.loc[i, 'flow_rate']) if i in m.FSU - # else (None,100), + # else (None,100), # Upper bound for the flow rate from Ruiz and Grossmann (2009) is 100 else (0, feed['flow_rate'].sum()) ), ) @@ -255,15 +246,12 @@ def build_model(approximation='quadratic'): m.TU, domain=pyo.NonNegativeReals, doc='CTUk cost of TU', - # bounds=lambda _,tu: (0, m.beta[tu] * feed['flow_rate'].sum() + m.gamma[tu] + m.theta[tu] * 41) # Q = optimal flow t1 - # bounds=lambda _,tu: (0, m.beta[tu] * feed['flow_rate'].sum() + m.gamma[tu] + m.theta[tu] * feed['flow_rate'].sum()) bounds=lambda _, tu: ( 0, m.beta[tu] * feed['flow_rate'].sum() + m.gamma[tu] + m.theta[tu] * feed['flow_rate'].sum() ** 0.7, ), - # initialize=lambda _,tu: m.beta[tu] * m.flow_into[mt] + m.gamma[tu] + m.theta[tu] * m.flow_into[mt]**0.7 for mt,t in m.MU_TU_streams if tu==t ) # ============================================================================= @@ -290,7 +278,6 @@ def _flow_into(m, option): """ return sum(m.flow[src, sink] for src, sink in m.streams if sink == option) - # @m.Expression(m.units, m.contaminant) @m.Expression(m.units - m.TU - m.outTU, m.contaminant) def _conc_into(m, option, j): """ @@ -318,7 +305,6 @@ def _conc_into(m, option, j): ) return pyo.Expression.Skip - # @m.Expression(m.units) @m.Expression(m.units - m.TU - m.inTU) def _flow_out_from(m, option): """ @@ -338,7 +324,6 @@ def _flow_out_from(m, option): """ return sum(m.flow[src, sink] for src, sink in m.streams if src == option) - # @m.Expression(m.units, m.contaminant) @m.Expression(m.units - m.inTU, m.contaminant) def _conc_out_from(m, option, j): """The mass balance of the contaminant in the unit to compute the concentration out of the unit. @@ -378,7 +363,6 @@ def _conc_out_from(m, option, j): m.mixer_balances.add(m._conc_into[mu, j] == m._conc_out_from[mu, j]) for j in m.contaminant ] - # m.mixer_balances_global.add(sum(m._conc_into[mu,sol] for sol in m.SOL) == sum(m._conc_out_from[mu,sol] for sol in m.SOL)) m.splitter_balances = pyo.ConstraintList() @@ -504,10 +488,8 @@ def unit_exists_or_not(m, unit): filter=lambda _, x, y: x == unit or y == unit, ) ue.MU_TU_streams = pyo.Set( - # doc='Mixer to treatment unit streams', doc="MU to TU 1-1 port pairing", initialize=m.inTU * m.TU, - # within=ue.streams, filter=lambda _, x, y: re.findall(r'\d+', x) == re.findall(r'\d+', y) and y == unit, ) @@ -516,9 +498,8 @@ def unit_exists_or_not(m, unit): ue.streams, doc="TU streams flowrate", domain=pyo.NonNegativeReals, - # bounds=lambda _,i,k:(TU.loc[unit,'L'],100) + # bounds=lambda _,i,k:(TU.loc[unit,'L'],100) # Upper bound for the flow rate from Ruiz and Grossmann (2009) is 100 bounds=lambda _, i, k: (TU.loc[unit, 'L'], feed['flow_rate'].sum()), - # initialize=lambda _,i,k: ) ue.conc = pyo.Var( @@ -684,7 +665,7 @@ def _cost_nv(ue): TU.loc[unit, 'L'] ** 0.7, feed['flow_rate'].sum() ** 0.7, ), - # bounds= lambda _,mt,unit: (TU.loc[unit,'L']**0.7,100**0.7), + # bounds= lambda _,mt,unit: (TU.loc[unit,'L']**0.7,100**0.7), # Upper bound for the flow rate from Ruiz and Grossmann (2009) is 100 doc="New var for potential term in capital cost", ) @@ -697,13 +678,6 @@ def _cost_nv(ue): for i in range(PieceCnt + 2): bpts.append(float((i * Topx) / PieceCnt)) - # def _breakpoints(lb, ub): - # x = np.linspace(lb, ub, 100) - # my_pwlf = pwlf.PiecewiseLinFit(x, x**0.7, degree=1) - # # my_pwlf.fit(8) - # my_pwlf.fitfast(10, pop=100) - # return my_pwlf.fit_breaks.tolist() - def _func(model, i, j, xp): """This function provides the expression for the piecewise linear approximation of the capital cost using the Picewise class. @@ -733,14 +707,9 @@ def _func(model, i, j, xp): pw_pts=bpts, pw_constr_type='EQ', f_rule=_func, - # ue.MU_TU_streams, ue.cost_var, ue.flow, pw_pts=_breakpoints(0,Topx), pw_constr_type='EQ', f_rule=func, - # pw_repn='BIGM_BIN' pw_repn='INC', ) - # for i,j in ue.MU_TU_streams: - # ue.ComputeObj[i,j].SOS2_y.setub(1.) - # Setting bounds for the INC_delta variable which is a binary variable that indicates the piecewise linear segment. for i, j in ue.MU_TU_streams: ue.ComputeObj[i, j].INC_delta.setub(1) @@ -780,7 +749,3 @@ def costTU(ue): # init_vars.InitMidpoint().apply_to(m) return m - - -# if __name__ == "__main__": -# build_model() From cd21c072320d24266a8b3588279cd3d5cc950e98 Mon Sep 17 00:00:00 2001 From: Carolina Tristan Date: Tue, 10 Sep 2024 11:27:22 -0400 Subject: [PATCH 11/21] Refactor README.md in water_network directory --- gdplib/water_network/README.md | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/gdplib/water_network/README.md b/gdplib/water_network/README.md index 2be0b37..bd40c1a 100644 --- a/gdplib/water_network/README.md +++ b/gdplib/water_network/README.md @@ -9,6 +9,12 @@ The mass balances are defined in terms of total flows and contaminants concentra Nonconvexities arise from bilinear terms “flows times concentration” in the mixers mass balances and concave investment cost functions of treatment units. The instance incorporates two approximations of the concave cost term (piecewise linear and quadratic) to reformulate the GDP model into a bilinear quadratic one. +The user can create each instance like this: + +``` +build_model(approximation='quadratic') +build_model(approximation='piecewise') +``` The general model description can be summarized as follows: Min Cost of Treatment Units From 6e3333c3bea6d08141406a69d94ec4dc3fec5dcc Mon Sep 17 00:00:00 2001 From: Carolina Tristan Date: Mon, 16 Sep 2024 11:51:57 -0400 Subject: [PATCH 12/21] Refactor capital cost approximation in `wnd.py` --- gdplib/water_network/wnd.py | 66 +++++++++++++++++++++++++++++++++++-- 1 file changed, 63 insertions(+), 3 deletions(-) diff --git a/gdplib/water_network/wnd.py b/gdplib/water_network/wnd.py index ac05460..36caa48 100644 --- a/gdplib/water_network/wnd.py +++ b/gdplib/water_network/wnd.py @@ -12,7 +12,12 @@ The mass balances are defined in terms of total flows and contaminants concentration. Nonconvexities arise from bilinear terms “flows times concentration” in the mixers mass balances and concave investment cost functions of treatment units. -The instance incorporates two approximations of the concave cost term (piecewise linear and quadratic) to reformulate the GDP model into a bilinear quadratic one. +The instance incorporates two approximations of the concave cost term (piecewise linear and quadratic) to reformulate the GDP model (approximation = 'none') into a bilinear quadratic one. +The user can create each instance like this: + +build_model(approximation='none') +build_model(approximation='quadratic') +build_model(approximation='piecewise') The general model description can be summarized as follows: Min Cost of Treatment Units @@ -87,6 +92,7 @@ def build_model(approximation='quadratic'): A Pyomo ConcreteModel object for the water network design. """ + print(f"Using {approximation} approximation for the capital cost term.") m = pyo.ConcreteModel('Water Network Design') # ============================================================================ @@ -574,6 +580,8 @@ def unit_exists_or_not(m, unit): if t == unit ] + # Approximation of the concave capital cost term for the active treatment units in the objective function. + if approximation == 'quadratic': # New variable for potential term in capital cost. Z = sum(Q)**0.7, Z >= 0 ue.cost_var = pyo.Var( @@ -628,7 +636,7 @@ def _func(x, a, b): # curve_fit is used to fit the quadratic curve to the function Z^(1/0.7) given the flow rate and the quadratic expression _func. popt, pcov = curve_fit(_func, z, z ** (1 / 0.7)) - return popt[0] * x + popt[1] * x**2 + return _func(x, *popt) # Z^(1/0.7)=sum(Q) @ue.Constraint(doc='New var potential term in capital cost cstr.') @@ -655,6 +663,51 @@ def _cost_nv(ue): == ue.flow[mt, unit] ) + elif approximation == "quadratic2": + + def _func2(x, a, b, c): + """This function computes the quadratic curve fit. + The curve fit is a quadratic function of the form a + b*x + c*x^2 given the flow rate and the coefficients a, b, and c obtained from the curve_fit function. + + Parameters + ---------- + x : variable + The flow rate. + a : float + The coefficient a. + b : float + The coefficient b. + c : float + The coefficient c. + + Returns + ------- + Expression + The quadratic curve fit for the function q^0.7. + """ + return a + b * x + c * x**2 + + def _g(x): + """This function provides the expression for the quadratic curve fit of the capital cost using the curve_fit function. + + Parameters + ---------- + x : variable + The flow rate. + + Returns + ------- + Expression + The quadratic curve fit of the capital cost + """ + # Equally spaced points between 0 and 100 + q = np.linspace(0, 100, 100) + # curve_fit is used to fit the quadratic curve to the function q^0.7 given the flow rate and the quadratic expression func2. + # popt is the optimal values for the coefficients a, b, and c, pcov is the covariance matrix + # func is a quadratic function of the form a + b*x + c*x^2. + popt, pcov = curve_fit(_func2, q, q**0.7) + return _func2(x, *popt) + elif approximation == "piecewise": # New variable for potential term in capital cost. Z = sum(Q)**0.7, Z >= 0 ue.cost_var = pyo.Var( @@ -718,7 +771,10 @@ def _func(model, i, j, xp): @ue.Constraint(doc='Cost active TU') def costTU(ue): """Constraint: Cost of active treatment unit. - The constraint ensures that the cost of the active treatment unit is equal to the sum of an investment cost which is proportional to the total flow to 0.7 exponent and an operating cost which is proportional to the flow + The constraint ensures that the cost of the active treatment unit is equal to the sum of an investment cost which is proportional to the total flow to 0.7 exponent and an operating cost which is proportional to the flow. + If approximation is quadratic, the investment cost is approximated by a quadratic function of the flow rate. + If approximation is piecewise, the investment cost is approximated by a piecewise linear function of the flow rate. + If approximation is none, the investment cost is equal to the flow rate to the 0.7 exponent, the original concave function. Parameters ---------- @@ -733,8 +789,12 @@ def costTU(ue): for mt, t in ue.streams: if approximation == 'quadratic': new_var = ue.cost_var[unit] + elif approximation == 'quadratic2': + new_var = _g(ue.flow[mt, unit]) elif approximation == 'piecewise': new_var = ue.cost_var[mt, unit] + elif approximation == 'none': + new_var = ue.flow[mt, unit] ** 0.7 if t == unit: return ( m.costTU[unit] From ac4acaeca23b8faafaaea9fc1d3eb67b9c07d482 Mon Sep 17 00:00:00 2001 From: Carolina Tristan Date: Mon, 16 Sep 2024 11:59:23 -0400 Subject: [PATCH 13/21] Refactor README.md in `water_network` directory --- gdplib/water_network/README.md | 25 +++++++++++++------------ 1 file changed, 13 insertions(+), 12 deletions(-) diff --git a/gdplib/water_network/README.md b/gdplib/water_network/README.md index bd40c1a..146f163 100644 --- a/gdplib/water_network/README.md +++ b/gdplib/water_network/README.md @@ -8,10 +8,11 @@ The fouled feed waters can be allocated to one or more treatment units or dispos The mass balances are defined in terms of total flows and contaminants concentration. Nonconvexities arise from bilinear terms “flows times concentration” in the mixers mass balances and concave investment cost functions of treatment units. -The instance incorporates two approximations of the concave cost term (piecewise linear and quadratic) to reformulate the GDP model into a bilinear quadratic one. +The instance incorporates two approximations of the concave cost term (piecewise linear and quadratic) to reformulate the orginal GDP model into a bilinear quadratic one. The user can create each instance like this: ``` +build_model(approximation='none') build_model(approximation='quadratic') build_model(approximation='piecewise') ``` @@ -41,20 +42,20 @@ The contaminant concentration and flow rate of the feed streams, contaminant rec ### Solution -Best known objective value: $ 348,340 +Best known objective value: $348,340 ### Size -| Component | pwl | quadratic | -| :-------------------- | :--: | :-------: | -| variables | 1405 | 420 | -| binary_variables | 510 | 10 | -| integer_variables | 0 | 0 | -| continuous_variables | 895 | 410 | -| disjunctions | 5 | 5 | -| disjuncts | 10 | 10 | -| constraints | 1339 | 334 | -| nonlinear_constraints | 28 | 33 | +| Component | original | pwl | quadratic | +| :-------------------- | -------- | :--: | :-------: | +| variables | 395 | 1405 | 420 | +| binary_variables | 10 | 510 | 10 | +| integer_variables | 0 | 0 | 0 | +| continuous_variables | 385 | 895 | 410 | +| disjunctions | 5 | 5 | 5 | +| disjuncts | 10 | 10 | 10 | +| constraints | 329 | 1339 | 334 | +| nonlinear_constraints | 33 | 28 | 33 | ### References From d9a3031e163869545a748874b4913f32f2c07ee2 Mon Sep 17 00:00:00 2001 From: Carolina Tristan Date: Mon, 16 Sep 2024 12:03:50 -0400 Subject: [PATCH 14/21] Fixed typo in README.md --- gdplib/water_network/README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gdplib/water_network/README.md b/gdplib/water_network/README.md index 146f163..60d673d 100644 --- a/gdplib/water_network/README.md +++ b/gdplib/water_network/README.md @@ -8,7 +8,7 @@ The fouled feed waters can be allocated to one or more treatment units or dispos The mass balances are defined in terms of total flows and contaminants concentration. Nonconvexities arise from bilinear terms “flows times concentration” in the mixers mass balances and concave investment cost functions of treatment units. -The instance incorporates two approximations of the concave cost term (piecewise linear and quadratic) to reformulate the orginal GDP model into a bilinear quadratic one. +The instance incorporates two approximations of the concave cost term (piecewise linear and quadratic) to reformulate the original GDP model into a bilinear quadratic one. The user can create each instance like this: ``` From 13d7a7c63f8331f800b8ab012090547335013dfb Mon Sep 17 00:00:00 2001 From: Carolina Tristan Date: Wed, 9 Oct 2024 15:53:38 -0400 Subject: [PATCH 15/21] Refactor references in README.md --- gdplib/water_network/README.md | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/gdplib/water_network/README.md b/gdplib/water_network/README.md index 60d673d..4c34f55 100644 --- a/gdplib/water_network/README.md +++ b/gdplib/water_network/README.md @@ -38,7 +38,7 @@ Assumptions: ### Case Study The WTN comprises five inlet streams with four contaminants and four treatment units. -The contaminant concentration and flow rate of the feed streams, contaminant recovery rates, minimum flow rate and cost coefficients of the treatment units, and the upper limit on the molar flow of contaminant in the purified stream, are reported in [1]. +The contaminant concentration and flow rate of the feed streams, contaminant recovery rates, minimum flow rate and cost coefficients of the treatment units, and the upper limit on the molar flow of contaminant in the purified stream, are reported in [2]. ### Solution @@ -59,6 +59,8 @@ Best known objective value: $348,340 ### References -> [1] Ruiz J., Grossmann I. E. Water Treatment Network Design. 2009 Available from CyberInfrastructure for [MINLP](www.minlp.org), a collaboration of Carnegie Mellon University and IBM at: www.minlp.org/library/problem/index.php?i=24 +> [1] Tristán C., Fallanza M., Ibáñez R., Grossmann I. E., and Bernal Neira D. E. (2024). Global Optimization via Quadratic Disjunctive Programming for Water Networks Design with Energy Recovery. Computer Aided Chemical Engineering, 53, 2161–2166. https://doi.org/10.1016/B978-0-443-28824-1.50361-6 > -> [2] Ruiz, J., & Grossmann, I. E. (2011). Using redundancy to strengthen the relaxation for the global optimization of MINLP problems. Computers & Chemical Engineering, 35(12), 2729–2740. https://doi.org/10.1016/J.COMPCHEMENG.2011.01.035 +> [2] Ruiz J., and Grossmann I. E. Water Treatment Network Design. 2009 Available from CyberInfrastructure for [MINLP](www.minlp.org), a collaboration of Carnegie Mellon University and IBM at: www.minlp.org/library/problem/index.php?i=24 +> +> [3] Ruiz, J., and Grossmann, I. E. (2011). Using redundancy to strengthen the relaxation for the global optimization of MINLP problems. Computers & Chemical Engineering, 35(12), 2729–2740. https://doi.org/10.1016/J.COMPCHEMENG.2011.01.035 From 325890c0de97ddf059d2ac575a0d56c74fceb4df Mon Sep 17 00:00:00 2001 From: Carolina Tristan Date: Thu, 10 Oct 2024 12:36:09 -0400 Subject: [PATCH 16/21] Remove commented `bounds` in `ue.conc` --- gdplib/water_network/wnd.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/gdplib/water_network/wnd.py b/gdplib/water_network/wnd.py index 36caa48..922b40f 100644 --- a/gdplib/water_network/wnd.py +++ b/gdplib/water_network/wnd.py @@ -514,8 +514,6 @@ def unit_exists_or_not(m, unit): doc="TU streams concentration", domain=pyo.NonNegativeReals, bounds=lambda _, j, i, k: (0, 100) if i == unit else (0, 4), - # bounds=lambda _, j, i, k: (feed[j].min(),100) if i==unit else (feed[j].min(),4) - # bounds=lambda _, j, i, k: (feed[j].min(),feed[j].max()), # initialize= lambda _, j, i, k: feed.loc[i,j] if i in m.FSU else None ) From 493365fdbc6e9253d5e1caf963388ab73258d9b0 Mon Sep 17 00:00:00 2001 From: Carolina Tristan Date: Thu, 10 Oct 2024 12:42:33 -0400 Subject: [PATCH 17/21] Update `ue` to `unit_exists` in `wnd.py` for better clarity. --- gdplib/water_network/wnd.py | 125 ++++++++++++++++++++---------------- 1 file changed, 70 insertions(+), 55 deletions(-) diff --git a/gdplib/water_network/wnd.py b/gdplib/water_network/wnd.py index 922b40f..465dc35 100644 --- a/gdplib/water_network/wnd.py +++ b/gdplib/water_network/wnd.py @@ -487,30 +487,30 @@ def unit_exists_or_not(m, unit): # m.atleast_oneRU = pyo.LogicalConstraint(expr=atleast(1, m.Yunit)) for unit in m.TU: - ue = m.unit_exists[unit] - ue.streams = pyo.Set( + unit_exists = m.unit_exists[unit] + unit_exists.streams = pyo.Set( doc="Streams in active TU", initialize=m.TU_streams, filter=lambda _, x, y: x == unit or y == unit, ) - ue.MU_TU_streams = pyo.Set( + unit_exists.MU_TU_streams = pyo.Set( doc="MU to TU 1-1 port pairing", initialize=m.inTU * m.TU, filter=lambda _, x, y: re.findall(r'\d+', x) == re.findall(r'\d+', y) and y == unit, ) - ue.flow = pyo.Var( - ue.streams, + unit_exists.flow = pyo.Var( + unit_exists.streams, doc="TU streams flowrate", domain=pyo.NonNegativeReals, # bounds=lambda _,i,k:(TU.loc[unit,'L'],100) # Upper bound for the flow rate from Ruiz and Grossmann (2009) is 100 bounds=lambda _, i, k: (TU.loc[unit, 'L'], feed['flow_rate'].sum()), ) - ue.conc = pyo.Var( + unit_exists.conc = pyo.Var( m.contaminant, - ue.streams, + unit_exists.streams, doc="TU streams concentration", domain=pyo.NonNegativeReals, bounds=lambda _, j, i, k: (0, 100) if i == unit else (0, 4), @@ -518,61 +518,74 @@ def unit_exists_or_not(m, unit): ) # Adding the mass balances for the active treatment units - ue.balances_con = pyo.ConstraintList(doc='TU Material balances') + unit_exists.balances_con = pyo.ConstraintList(doc='TU Material balances') # The flowrate at the inlet of the treatment unit is equal to the flowrate at the outlet of the treatment unit. [ - ue.balances_con.add(ue.flow[mt, unit] == ue.flow[unit, st]) - for mt, t in ue.streams + unit_exists.balances_con.add( + unit_exists.flow[mt, unit] == unit_exists.flow[unit, st] + ) + for mt, t in unit_exists.streams if t == unit - for t, st in ue.streams + for t, st in unit_exists.streams if t == unit ] # The concentration of the contaminant at the inlet of the treatment unit is equal to the concentration of the contaminant at the outlet of the treatment unit times the removal ratio. [ - ue.balances_con.add( - ue.conc[j, unit, st] - == (1 - m.removal_ratio[j, t]) * ue.conc[j, mt, unit] + unit_exists.balances_con.add( + unit_exists.conc[j, unit, st] + == (1 - m.removal_ratio[j, t]) * unit_exists.conc[j, mt, unit] ) - for mt, t in ue.streams + for mt, t in unit_exists.streams if t == unit - for t, st in ue.streams + for t, st in unit_exists.streams if t == unit for j in m.contaminant ] # Treatment unit's mixer mass balance on the flowrate. [ - ue.balances_con.add(m._flow_into[mt] == ue.flow[mt, unit]) - for mt, t in ue.streams + unit_exists.balances_con.add(m._flow_into[mt] == unit_exists.flow[mt, unit]) + for mt, t in unit_exists.streams if t == unit ] # Treatment unit's mixer mass balance on the concentration of contaminants. [ - ue.balances_con.add( - m._conc_into[mt, j] == ue.conc[j, mt, unit] * ue.flow[mt, unit] + unit_exists.balances_con.add( + m._conc_into[mt, j] + == unit_exists.conc[j, mt, unit] * unit_exists.flow[mt, unit] ) - for mt, t in ue.streams + for mt, t in unit_exists.streams if t == unit for j in m.contaminant ] # Treatment unit's splitter mass balance on the flowrate. [ - ue.balances_con.add(ue.flow[unit, st] == m._flow_out_from[st]) - for t, st in ue.streams + unit_exists.balances_con.add( + unit_exists.flow[unit, st] == m._flow_out_from[st] + ) + for t, st in unit_exists.streams if t == unit ] # Treatment unit's splitter mass balance on the concentration of contaminants. [ - ue.balances_con.add(ue.conc[j, src2, sink2] == m.conc[j, src1, sink1]) + unit_exists.balances_con.add( + unit_exists.conc[j, src2, sink2] == m.conc[j, src1, sink1] + ) for src1, sink1 in m.from_splitters - for src2, sink2 in ue.streams + for src2, sink2 in unit_exists.streams if src2 == unit and src1 == sink2 for j in m.contaminant ] # Setting inlet flowrate bounds for the active treatment units. - ue.flow_bound = pyo.ConstraintList(doc='Flowrate bounds to/from active RU') + unit_exists.flow_bound = pyo.ConstraintList( + doc='Flowrate bounds to/from active RU' + ) [ - ue.flow_bound.add( - (ue.flow[mt, unit].lb, m._flow_into[mt], ue.flow[mt, unit].ub) + unit_exists.flow_bound.add( + ( + unit_exists.flow[mt, unit].lb, + m._flow_into[mt], + unit_exists.flow[mt, unit].ub, + ) ) for mt, t in m.MU_TU_streams if t == unit @@ -582,7 +595,7 @@ def unit_exists_or_not(m, unit): if approximation == 'quadratic': # New variable for potential term in capital cost. Z = sum(Q)**0.7, Z >= 0 - ue.cost_var = pyo.Var( + unit_exists.cost_var = pyo.Var( m.TU, domain=pyo.NonNegativeReals, initialize=lambda _, unit: TU.loc[unit, 'L'] ** 0.7, @@ -637,14 +650,14 @@ def _func(x, a, b): return _func(x, *popt) # Z^(1/0.7)=sum(Q) - @ue.Constraint(doc='New var potential term in capital cost cstr.') - def _cost_nv(ue): + @unit_exists.Constraint(doc='New var potential term in capital cost cstr.') + def _cost_nv(unit_exists): """Constraint: New variable potential term in capital cost. The constraint ensures that the new variable potential term in the capital cost is equal to the flow rate entering the treatment unit. Parameters ---------- - ue : Disjunct + unit_exists : Disjunct The disjunct for the active treatment unit. Returns @@ -652,13 +665,15 @@ def _cost_nv(ue): Constraint The constraint that the new variable potential term in the capital cost is equal to the flow rate. """ - for mt, t in ue.streams: + for mt, t in unit_exists.streams: if t == unit: return ( _quadratic_curve_fit( - 0, ue.cost_var[unit].ub, ue.cost_var[unit] + 0, + unit_exists.cost_var[unit].ub, + unit_exists.cost_var[unit], ) - == ue.flow[mt, unit] + == unit_exists.flow[mt, unit] ) elif approximation == "quadratic2": @@ -708,8 +723,8 @@ def _g(x): elif approximation == "piecewise": # New variable for potential term in capital cost. Z = sum(Q)**0.7, Z >= 0 - ue.cost_var = pyo.Var( - ue.MU_TU_streams, + unit_exists.cost_var = pyo.Var( + unit_exists.MU_TU_streams, domain=pyo.NonNegativeReals, initialize=lambda _, mt, unit: TU.loc[unit, 'L'] ** 0.7, bounds=lambda _, mt, unit: ( @@ -723,9 +738,9 @@ def _g(x): # to avoid warnings, we set breakpoints at or beyond the bounds PieceCnt = 100 bpts = [] - for mt, t in ue.streams: + for mt, t in unit_exists.streams: if t == unit: - Topx = ue.flow[mt, unit].ub + Topx = unit_exists.flow[mt, unit].ub for i in range(PieceCnt + 2): bpts.append(float((i * Topx) / PieceCnt)) @@ -751,10 +766,10 @@ def _func(model, i, j, xp): # Piecewise is a class that provides a piecewise linear approximation of a function using the INC representation. # The INC representation is a piecewise linear approximation of a function using the incremental form. - ue.ComputeObj = pyo.Piecewise( - ue.MU_TU_streams, - ue.cost_var, - ue.flow, + unit_exists.ComputeObj = pyo.Piecewise( + unit_exists.MU_TU_streams, + unit_exists.cost_var, + unit_exists.flow, pw_pts=bpts, pw_constr_type='EQ', f_rule=_func, @@ -762,12 +777,12 @@ def _func(model, i, j, xp): ) # Setting bounds for the INC_delta variable which is a binary variable that indicates the piecewise linear segment. - for i, j in ue.MU_TU_streams: - ue.ComputeObj[i, j].INC_delta.setub(1) - ue.ComputeObj[i, j].INC_delta.setlb(0) + for i, j in unit_exists.MU_TU_streams: + unit_exists.ComputeObj[i, j].INC_delta.setub(1) + unit_exists.ComputeObj[i, j].INC_delta.setlb(0) - @ue.Constraint(doc='Cost active TU') - def costTU(ue): + @unit_exists.Constraint(doc='Cost active TU') + def costTU(unit_exists): """Constraint: Cost of active treatment unit. The constraint ensures that the cost of the active treatment unit is equal to the sum of an investment cost which is proportional to the total flow to 0.7 exponent and an operating cost which is proportional to the flow. If approximation is quadratic, the investment cost is approximated by a quadratic function of the flow rate. @@ -776,7 +791,7 @@ def costTU(ue): Parameters ---------- - ue : Disjunct + unit_exists : Disjunct The disjunct for the active treatment unit. Returns @@ -784,19 +799,19 @@ def costTU(ue): Constraint The constraint that the cost of the active treatment unit is equal to the sum of an investment cost and an operating cost. """ - for mt, t in ue.streams: + for mt, t in unit_exists.streams: if approximation == 'quadratic': - new_var = ue.cost_var[unit] + new_var = unit_exists.cost_var[unit] elif approximation == 'quadratic2': - new_var = _g(ue.flow[mt, unit]) + new_var = _g(unit_exists.flow[mt, unit]) elif approximation == 'piecewise': - new_var = ue.cost_var[mt, unit] + new_var = unit_exists.cost_var[mt, unit] elif approximation == 'none': - new_var = ue.flow[mt, unit] ** 0.7 + new_var = unit_exists.flow[mt, unit] ** 0.7 if t == unit: return ( m.costTU[unit] - == m.beta[unit] * ue.flow[mt, unit] + == m.beta[unit] * unit_exists.flow[mt, unit] + m.gamma[unit] + m.theta[unit] * new_var ) From 3d68a252cd6f6f562b1f045b2d0216b88c0fb576 Mon Sep 17 00:00:00 2001 From: Carolina Tristan Date: Thu, 10 Oct 2024 18:45:56 -0400 Subject: [PATCH 18/21] Rename `quadratic2` and `quadratic` in `wnd.py` --- gdplib/water_network/wnd.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/gdplib/water_network/wnd.py b/gdplib/water_network/wnd.py index 465dc35..249c15e 100644 --- a/gdplib/water_network/wnd.py +++ b/gdplib/water_network/wnd.py @@ -593,7 +593,7 @@ def unit_exists_or_not(m, unit): # Approximation of the concave capital cost term for the active treatment units in the objective function. - if approximation == 'quadratic': + if approximation == 'quadratic_zero_origin': # New variable for potential term in capital cost. Z = sum(Q)**0.7, Z >= 0 unit_exists.cost_var = pyo.Var( m.TU, @@ -676,7 +676,7 @@ def _cost_nv(unit_exists): == unit_exists.flow[mt, unit] ) - elif approximation == "quadratic2": + elif approximation == "quadratic_nonzero_origin": def _func2(x, a, b, c): """This function computes the quadratic curve fit. From ab9db0c7265fea7d52c759c71e288a42f3835dee Mon Sep 17 00:00:00 2001 From: Carolina Tristan Date: Thu, 10 Oct 2024 18:47:46 -0400 Subject: [PATCH 19/21] Refactor quadratic approximations in README.md --- gdplib/water_network/README.md | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/gdplib/water_network/README.md b/gdplib/water_network/README.md index 4c34f55..1257dd4 100644 --- a/gdplib/water_network/README.md +++ b/gdplib/water_network/README.md @@ -8,12 +8,22 @@ The fouled feed waters can be allocated to one or more treatment units or dispos The mass balances are defined in terms of total flows and contaminants concentration. Nonconvexities arise from bilinear terms “flows times concentration” in the mixers mass balances and concave investment cost functions of treatment units. -The instance incorporates two approximations of the concave cost term (piecewise linear and quadratic) to reformulate the original GDP model into a bilinear quadratic one. +The instance incorporates three approximations of the concave cost term when the treatment unit exist, including one piecewise linear and two quadratic approximations, to reformulate the original GDP model into a bilinear quadratic problem [1]. + + +The quadratic reformulations are as follows: + +* `quadratic_zero_origin` which reads as $f(x) = a \ x + b \ x^2$. Absence of flow leads to no costs for the pump, which is what we expect. Optimal solution: $346,654 with a relative error with respect to the best known objective value of 0.5% +* `quadratic_nonzero_origin` takes the form of $f(x) = a + b \ x + b \ x^2$. This constraint leads to a non-zero pump expense in the absence of flow. Optimal solution: $349,521 with a relative error with respect to the best known objective value of 0.3%. + +The two quadratic approximations are both effective in capturing pump costs, but `quadratic_nonzero_origin` provides a better fit for the active treatment unit's flowrate range. + The user can create each instance like this: ``` build_model(approximation='none') -build_model(approximation='quadratic') +build_model(approximation='quadratic_zero_origin') +build_model(approximation='quadratic_nonzero_origin') build_model(approximation='piecewise') ``` From da82210c0d23624222ae1809f020c37506c9e240 Mon Sep 17 00:00:00 2001 From: Carolina Tristan Date: Thu, 10 Oct 2024 19:04:15 -0400 Subject: [PATCH 20/21] Refactor references in `wnd.py` --- gdplib/water_network/wnd.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/gdplib/water_network/wnd.py b/gdplib/water_network/wnd.py index 249c15e..81f0b2d 100644 --- a/gdplib/water_network/wnd.py +++ b/gdplib/water_network/wnd.py @@ -44,6 +44,8 @@ References ---------- +Tristán C., Fallanza M., Ibáñez R., Grossmann I. E., and Bernal Neira D. E. (2024). Global Optimization via Quadratic Disjunctive Programming for Water Networks Design with Energy Recovery. Computer Aided Chemical Engineering, 53, 2161–2166. https://doi.org/10.1016/B978-0-443-28824-1.50361-6 + Ruiz J., Grossmann IE. Water Treatment Network Design. 2009 Available from CyberInfrastructure for [MINLP](), a collaboration of Carnegie Mellon University and IBM at: www.minlp.org/library/problem/index.php?i=24 From 6b68964b8542e90925bc0fb3f17c12d7ac598649 Mon Sep 17 00:00:00 2001 From: Carolina Tristan Date: Mon, 4 Nov 2024 08:07:16 -0500 Subject: [PATCH 21/21] Add comment to clarify variable initialization in `build_model` function --- gdplib/water_network/wnd.py | 1 + 1 file changed, 1 insertion(+) diff --git a/gdplib/water_network/wnd.py b/gdplib/water_network/wnd.py index 81f0b2d..6471dc2 100644 --- a/gdplib/water_network/wnd.py +++ b/gdplib/water_network/wnd.py @@ -821,6 +821,7 @@ def costTU(unit_exists): # Objective function: minimize the total cost of the treatment units m.obj = pyo.Objective(expr=sum(m.costTU[k] for k in m.TU), sense=pyo.minimize) + # Initialization of the variables to midpoint between the bounds to avoid numerical issues if no initial values are provided. # init_vars.InitMidpoint().apply_to(m) return m