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.

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

BSM2-P Effluent Metrics #1489

Closed
35 changes: 17 additions & 18 deletions watertap/flowsheets/full_water_resource_recovery_facility/BSM2.py
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
Original file line number Diff line number Diff line change
Expand Up @@ -89,11 +89,13 @@
cost_primary_clarifier,
)

from idaes.core.util.model_diagnostics import DiagnosticsToolbox

# Set up logger
_log = idaeslog.getLogger(__name__)


def main(bio_P=False):
def main(bio_P=False, reactor_volume_equalities=False):
m = build(bio_P=bio_P)
set_operating_conditions(m)

Expand All @@ -120,7 +122,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,12 +133,16 @@ def main(bio_P=False):
fail_flag=True,
)

# Re-solve with effluent violation constraints
setup_optimization(m, reactor_volume_equalities=reactor_volume_equalities)
results = solve(m)

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

interval_initializer(m.fs.costing)

assert_degrees_of_freedom(m, 0)
# assert_degrees_of_freedom(m, 0)

results = solve(m)
pyo.assert_optimal_termination(results)
Expand Down Expand Up @@ -552,6 +558,14 @@ def scale_variables(m):
scale_variables(m)
iscale.calculate_scaling_factors(m)

# Touch treated properties
m.fs.Treated.properties[0].TSS
m.fs.Treated.properties[0].COD
m.fs.Treated.properties[0].BOD5
m.fs.Treated.properties[0].SNKj
m.fs.Treated.properties[0].SP_organic
m.fs.Treated.properties[0].SP_inorganic
Copy link
Contributor

Choose a reason for hiding this comment

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

You should touch these before the solve.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

These are touched before the solve. I placed it at the end of set_operating_conditions. Do you mean these should be touched right before the solve? How would that functionally be any different than the current placement?

Copy link
Contributor

Choose a reason for hiding this comment

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

My mistake- I misread the GH diff and it looked as if these lines were after the solve. That said, I wonder if it makes more sense to touch before applying scaling. Right now, I think you made each an expression so it might not matter. If they were reverted back to vars, then it would matter.



def initialize_system(m, bio_P=False, solver=None):
# Initialize flowsheet
Expand Down Expand Up @@ -693,6 +707,104 @@ def solve(m, solver=None):
return results


def setup_optimization(m, reactor_volume_equalities=False):

for i in ["R1", "R2", "R3", "R4", "R5", "R6", "R7"]:
reactor = getattr(m.fs, i)
reactor.volume.unfix()
reactor.volume.setlb(1)
# reactor.volume.setub(2000)
if reactor_volume_equalities:
add_reactor_volume_equalities(m)

m.fs.R5.outlet.conc_mass_comp[:, "S_O2"].unfix()
m.fs.R5.outlet.conc_mass_comp[:, "S_O2"].setub(1e-2)

m.fs.R6.outlet.conc_mass_comp[:, "S_O2"].unfix()
m.fs.R6.outlet.conc_mass_comp[:, "S_O2"].setub(1e-2)

m.fs.R7.outlet.conc_mass_comp[:, "S_O2"].unfix()
m.fs.R7.outlet.conc_mass_comp[:, "S_O2"].setub(1e-2)

# m.fs.R5.injection[:, :, :].unfix()
# m.fs.R6.injection[:, :, :].unfix()
# m.fs.R7.injection[:, :, :].unfix()

# Unfix fraction of outflow from reactor 7 that goes to recycle
m.fs.SP1.split_fraction[:, "underflow"].unfix()
# m.fs.SP1.split_fraction[:, "underflow"].setlb(0.45)
m.fs.SP2.split_fraction[:, "recycle"].unfix()

add_effluent_violations(m)


def add_reactor_volume_equalities(m):
# TODO: These constraints were applied for initial optimization of AS reactor volumes; otherwise, volumes drive towards lower bound. Revisit
@m.fs.Constraint(m.fs.time)
def Vol_1(self, t):
return m.fs.R1.volume[0] == m.fs.R2.volume[0]

@m.fs.Constraint(m.fs.time)
def Vol_2(self, t):
return m.fs.R3.volume[0] == m.fs.R4.volume[0]

@m.fs.Constraint(m.fs.time)
def Vol_3(self, t):
return m.fs.R5.volume[0] == m.fs.R6.volume[0]

@m.fs.Constraint(m.fs.time)
def Vol_4(self, t):
return m.fs.R7.volume[0] >= m.fs.R6.volume[0] * 0.5


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["effluent"] <= 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 Expand Up @@ -872,6 +984,20 @@ def display_performance_metrics(m):
pyo.units.get_units(m.fs.AD.liquid_phase.properties_in[0].flow_vol),
)

TSS = pyo.value(m.fs.Treated.properties[0].TSS)
COD = pyo.value(m.fs.Treated.properties[0].COD)
BOD = pyo.value(m.fs.Treated.properties[0].BOD5["effluent"])
SNKj = pyo.value(m.fs.Treated.properties[0].SNKj)
sp_org = pyo.value(m.fs.Treated.properties[0].SP_organic)
sp_inorg = pyo.value(m.fs.Treated.properties[0].SP_inorganic)

print(f"TSS Concentration: {TSS}")
print(f"COD Concentration: {COD}")
print(f"BOD Concentration: {BOD}")
print(f"SNKj Concentration: {SNKj}")
print(f"SP Organic Concentration: {sp_org}")
print(f"SP Inorganic Concentration: {sp_inorg}")


if __name__ == "__main__":
# This method builds and runs a steady state activated sludge flowsheet.
Expand Down
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Original file line number Diff line number Diff line change
Expand Up @@ -149,7 +149,7 @@ def build(self):
domain=pyo.PositiveReals,
doc="Conversion factor applied for TSS calculation",
)
self.BOD5_factor = pyo.Var(
self.BOD5_factor = pyo.Param(
["raw", "effluent"],
initialize={"raw": 0.65, "effluent": 0.25},
units=pyo.units.dimensionless,
Expand Down
Loading
Loading