Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

BSM2-P Effluent Metrics #1489

Closed
wants to merge 12 commits into from
Closed

BSM2-P Effluent Metrics #1489

Show file tree
Hide file tree
Changes from 2 commits
Commits
Show all changes
12 commits
Select commit Hold shift + click to select a range
0ed0d42
Add effluent metrics to prop pack
MarcusHolly Sep 12, 2024
9541a6a
Add preliminary effluent violations to the BSM2-P flowsheet
MarcusHolly Sep 12, 2024
07ea01a
Revert TSS change
MarcusHolly Sep 12, 2024
56b5145
Make BOD5_factor a parameter
MarcusHolly Sep 13, 2024
e2ac3b0
Update BSM2 GUI images
MarcusHolly Sep 13, 2024
4d05c61
Define effluent metrics as expressions (as done in ASM1) rather than …
MarcusHolly Sep 13, 2024
a09c691
Initial attempt at adding optimization to BSM2-P
MarcusHolly Sep 13, 2024
a19f1c6
Disable phosphorus effluent violation
MarcusHolly Sep 13, 2024
26df31e
Solves to Acceptable Level
MarcusHolly Sep 13, 2024
ef4f5e6
Update unit model tests
MarcusHolly Sep 13, 2024
c2ee4ac
Touch effluent metrics before scaling
MarcusHolly Sep 18, 2024
ec483ba
minor updates
MarcusHolly Sep 18, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
35 changes: 17 additions & 18 deletions watertap/flowsheets/full_water_resource_recovery_facility/BSM2.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -526,35 +526,34 @@ def setup_optimization(m, reactor_volume_equalities=False):
add_effluent_violations(m)


def add_effluent_violations(m):
# TODO: update "m" to blk; change ref to m.fs.Treated instead of CL1 effluent
Comment on lines -529 to -530
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I haven't been able to figure out why, but this change is breaking the BSM2 GUI, so I've reverted it in #1491

m.fs.TSS_max = pyo.Var(initialize=0.03, units=pyo.units.kg / pyo.units.m**3)
m.fs.TSS_max.fix()
def add_effluent_violations(blk):
blk.fs.TSS_max = pyo.Var(initialize=0.03, units=pyo.units.kg / pyo.units.m**3)
blk.fs.TSS_max.fix()

@m.fs.Constraint(m.fs.time)
@blk.fs.Constraint(blk.fs.time)
def eq_TSS_max(self, t):
return m.fs.CL1.effluent_state[0].TSS <= m.fs.TSS_max
return blk.fs.Treated.properties[0].TSS <= blk.fs.TSS_max

m.fs.COD_max = pyo.Var(initialize=0.1, units=pyo.units.kg / pyo.units.m**3)
m.fs.COD_max.fix()
blk.fs.COD_max = pyo.Var(initialize=0.1, units=pyo.units.kg / pyo.units.m**3)
blk.fs.COD_max.fix()

@m.fs.Constraint(m.fs.time)
@blk.fs.Constraint(blk.fs.time)
def eq_COD_max(self, t):
return m.fs.CL1.effluent_state[0].COD <= m.fs.COD_max
return blk.fs.Treated.properties[0].COD <= blk.fs.COD_max

m.fs.totalN_max = pyo.Var(initialize=0.018, units=pyo.units.kg / pyo.units.m**3)
m.fs.totalN_max.fix()
blk.fs.totalN_max = pyo.Var(initialize=0.018, units=pyo.units.kg / pyo.units.m**3)
blk.fs.totalN_max.fix()

@m.fs.Constraint(m.fs.time)
@blk.fs.Constraint(blk.fs.time)
def eq_totalN_max(self, t):
return m.fs.CL1.effluent_state[0].Total_N <= m.fs.totalN_max
return blk.fs.Treated.properties[0].Total_N <= blk.fs.totalN_max

m.fs.BOD5_max = pyo.Var(initialize=0.01, units=pyo.units.kg / pyo.units.m**3)
m.fs.BOD5_max.fix()
blk.fs.BOD5_max = pyo.Var(initialize=0.01, units=pyo.units.kg / pyo.units.m**3)
blk.fs.BOD5_max.fix()

@m.fs.Constraint(m.fs.time)
@blk.fs.Constraint(blk.fs.time)
def eq_BOD5_max(self, t):
return m.fs.CL1.effluent_state[0].BOD5["effluent"] <= m.fs.BOD5_max
return blk.fs.Treated.properties[0].BOD5["effluent"] <= blk.fs.BOD5_max


def add_reactor_volume_equalities(m):
Expand Down
54 changes: 53 additions & 1 deletion watertap/flowsheets/full_water_resource_recovery_facility/BSM2_P_extension.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -120,7 +120,7 @@ def main(bio_P=False):
m.fs.R6.outlet.conc_mass_comp[:, "S_O2"].unfix()
m.fs.R7.outlet.conc_mass_comp[:, "S_O2"].unfix()

# Resolve with controls in place
# Re-solve with controls in place
results = solve(m)

pyo.assert_optimal_termination(results)
Expand All @@ -131,6 +131,10 @@ def main(bio_P=False):
fail_flag=True,
)

# Re-solve with effluent violation constraints
add_effluent_violations(m)
results = solve(m)

add_costing(m)
m.fs.costing.initialize()

Expand Down Expand Up @@ -693,6 +697,54 @@ def solve(m, solver=None):
return results


def add_effluent_violations(blk):
# TODO: Revisit the max effluent concentration values

# Max value taken from Flores-Alsina Excel
blk.fs.TSS_max = pyo.Var(initialize=0.03, units=pyo.units.kg / pyo.units.m**3)
blk.fs.TSS_max.fix()

@blk.fs.Constraint(blk.fs.time)
def eq_TSS_max(self, t):
return blk.fs.Treated.properties[0].TSS <= blk.fs.TSS_max

# Max value carried over from BSM2
blk.fs.COD_max = pyo.Var(initialize=0.1, units=pyo.units.kg / pyo.units.m**3)
blk.fs.COD_max.fix()

@blk.fs.Constraint(blk.fs.time)
def eq_COD_max(self, t):
return blk.fs.Treated.properties[0].COD <= blk.fs.COD_max

# Max value taken from Flores-Alsina Excel
blk.fs.SNKj_max = pyo.Var(initialize=0.004, units=pyo.units.kg / pyo.units.m**3)
blk.fs.SNKj_max.fix()

@blk.fs.Constraint(blk.fs.time)
def eq_SNKj_max(self, t):
return blk.fs.Treated.properties[0].SNKj <= blk.fs.SNKj_max

# Max value carried over from BSM2
blk.fs.BOD5_max = pyo.Var(initialize=0.01, units=pyo.units.kg / pyo.units.m**3)
blk.fs.BOD5_max.fix()

@blk.fs.Constraint(blk.fs.time)
def eq_BOD5_max(self, t):
return blk.fs.Treated.properties[0].BOD5 <= blk.fs.BOD5_max

# Max value taken from Flores-Alsina Excel
blk.fs.total_P_max = pyo.Var(initialize=0.002, units=pyo.units.kg / pyo.units.m**3)
blk.fs.total_P_max.fix()

@blk.fs.Constraint(blk.fs.time)
def eq_total_P_max(self, t):
return (
blk.fs.Treated.properties[0].SP_organic
+ blk.fs.Treated.properties[0].SP_inorganic
<= blk.fs.total_P_max
)


def add_costing(m):
m.fs.costing = WaterTAPCosting()
m.fs.costing.base_currency = pyo.units.USD_2020
Expand Down
162 changes: 161 additions & 1 deletion watertap/property_models/unit_specific/activated_sludge/modified_asm2d_properties.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -210,6 +210,86 @@ def build(self):
doc="ISS fractional content of biomass",
)

# EQI parameters
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can you add some general refs so we know where to look for these in future?

self.i_NSF = pyo.Var(
initialize=0.03352,
units=pyo.units.dimensionless,
domain=pyo.NonNegativeReals,
doc="N content of fermentable substrate, S_F, [kg N/kg COD]",
)
self.i_NSI = pyo.Var(
initialize=0.06003,
units=pyo.units.dimensionless,
domain=pyo.NonNegativeReals,
doc="N content of inert soluble COD S_I, [kg N/kg COD]",
)
self.i_NXI = pyo.Var(
initialize=0.06003,
units=pyo.units.dimensionless,
domain=pyo.NonNegativeReals,
doc="N content of inert particulate COD X_I, [kg N/kg COD]",
)
self.i_NXS = pyo.Var(
initialize=0.03352,
units=pyo.units.dimensionless,
domain=pyo.NonNegativeReals,
doc="N content of slowly biodegradable substrate X_S, [kg N/kg COD]",
)
self.i_NBM = pyo.Var(
initialize=0.08615,
units=pyo.units.dimensionless,
domain=pyo.NonNegativeReals,
doc="N content of biomass, X_H, X_PAO, X_AUT, [kg N/kg COD]",
)
self.f_SI = pyo.Var(
initialize=0.00,
units=pyo.units.dimensionless,
domain=pyo.NonNegativeReals,
doc="Production of S_I in hydrolysis, [kg COD/kg COD]",
)
self.f_XIH = pyo.Var(
initialize=0.1,
units=pyo.units.dimensionless,
domain=pyo.NonNegativeReals,
doc="Fraction of inert COD generated in lysis of X_H, [kg COD/kg COD]",
)
self.f_XIP = pyo.Var(
initialize=0.1,
units=pyo.units.dimensionless,
domain=pyo.NonNegativeReals,
doc="Fraction of inert COD generated in lysis of X_PAO and X_PHA, [kg COD/kg COD]",
)
self.f_XIA = pyo.Var(
initialize=0.1,
units=pyo.units.dimensionless,
domain=pyo.NonNegativeReals,
doc="Fraction of inert COD generated in lysis of X_AUT, [kg COD/kg COD]",
)
self.i_PSF = pyo.Var(
initialize=0.00559,
units=pyo.units.dimensionless,
domain=pyo.NonNegativeReals,
doc="P content of fermentable substrate, S_F, [kg P/kg COD]",
)
self.i_PXI = pyo.Var(
initialize=0.00649,
units=pyo.units.dimensionless,
domain=pyo.NonNegativeReals,
doc="P content of inert particulate COD X_I, [kg P/kg COD]",
)
self.i_PXS = pyo.Var(
initialize=0.00559,
units=pyo.units.dimensionless,
domain=pyo.NonNegativeReals,
doc="P content of slowly biodegradable substrate X_S, [kg P/kg COD]",
)
self.i_PBM = pyo.Var(
initialize=0.02154,
units=pyo.units.dimensionless,
domain=pyo.NonNegativeReals,
doc="P content of biomass, X_H, X_PAO, X_AUT, [kg P/kg COD]",
)

# Fix Vars that are treated as Params
for v in self.component_objects(pyo.Var):
v.fix()
Expand Down Expand Up @@ -442,10 +522,90 @@ def _TSS(self):
)

def rule_TSS(b):
return b.TSS == b.VSS + b.ISS
return b.TSS == b.conc_mass_comp["X_H"]
MarcusHolly marked this conversation as resolved.
Show resolved Hide resolved

self.eq_TSS = pyo.Constraint(rule=rule_TSS)

def _COD(self):
self.COD = pyo.Var(
initialize=1,
domain=pyo.NonNegativeReals,
doc="Chemical oxygen demand",
units=pyo.units.kg / pyo.units.m**3,
)

def rule_COD(b):
return (
b.COD
== b.conc_mass_comp["S_F"]
+ b.conc_mass_comp["S_A"]
+ b.conc_mass_comp["S_I"]
+ b.conc_mass_comp["X_I"]
+ b.conc_mass_comp["X_S"]
+ b.conc_mass_comp["X_H"]
+ b.conc_mass_comp["X_PAO"]
+ b.conc_mass_comp["X_PHA"]
+ b.conc_mass_comp["X_AUT"]
)

self.eq_COD = pyo.Constraint(rule=rule_COD)

def rule_SNKj(b):
return b.SNKj == b.conc_mass_comp[
"S_NH4"
] + b.params.i_NSF * b.conc_mass_comp[
"S_F"
] + b.params.i_NSI * b.conc_mass_comp[
"S_I"
] + b.params.i_NXI * b.conc_mass_comp[
"X_I"
] + b.params.i_NXS * b.conc_mass_comp[
"X_S"
] + b.params.i_NBM * (
b.conc_mass_comp["X_H"]
+ b.conc_mass_comp["X_PAO"]
+ b.conc_mass_comp["X_AUT"]
)

self.eq_SNKj = pyo.Constraint(rule=rule_SNKj)

def rule_BOD5(b):
return b.BOD5 == 0.25 * (
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

From memory, I recall that the factor (0.25) could change depending on whether this is calculated for influent or effluent. Can you double-check this? I handled this in some way for ASM1 (though it might not be ideal approach).

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If I am recalling correctly, I would at least opt for using a mutable Param instead of the hard-coded 0.25. Then it would be on the user to change the Param value as needed at the flowsheet level.

Copy link
Contributor Author

@MarcusHolly MarcusHolly Sep 12, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am familiar with what you are talking about, but the c-code and the literature reference both hard-code 0.25, so I didn't see a need to do what you did in ASM1
image

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The problem here is if I want to apply to the influent, that value changes, right? Since this can be constructed on any state block and not only one that corresponds to treated effluent, we cannot just hard-code the value to what would correspond to the effluent.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can you show the same snapshot for BOD5i or whatever corresponds to influent in that file?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ah yes, the influent value is 0.65
image

b.conc_mass_comp["S_F"]
+ b.conc_mass_comp["S_A"]
+ (1 - b.params.f_SI) * b.conc_mass_comp["X_S"]
+ (1 - b.params.f_XIH) * b.conc_mass_comp["X_H"]
+ (1 - b.params.f_XIP)
* (b.conc_mass_comp["X_PAO"] + b.conc_mass_comp["X_PHA"])
+ (1 - b.params.f_XIA) * b.conc_mass_comp["X_AUT"]
)

self.eq_BOD5 = pyo.Constraint(rule=rule_BOD5)

def rule_SP_organic(b):
return b.SP_organic == b.conc_mass_comp[
"X_PP"
] + b.params.i_PSF * b.conc_mass_comp[
"S_F"
] + b.params.i_PSI * b.conc_mass_comp[
"S_I"
] + b.params.i_PXI * b.conc_mass_comp[
"X_I"
] + b.params.i_PXS * b.conc_mass_comp[
"X_S"
] + b.params.i_PBM * (
b.conc_mass_comp["X_H"]
+ b.conc_mass_comp["X_PAO"]
+ b.conc_mass_comp["X_AUT"]
)

self.eq_SP_organic = pyo.Constraint(rule=rule_SP_organic)

def rule_SP_inorganic(b):
return b.SP_inorganic == b.conc_mass_comp["S_PO4"]

self.eq_SP_inorganic = pyo.Constraint(rule=rule_SP_inorganic)

def get_material_flow_terms(self, p, j):
if j == "H2O":
return self.flow_vol * self.params.dens_mass
Expand Down
Loading