BSM2-P Effluent Metrics

watertap/flowsheets/full_water_resource_recovery_facility/BSM2.py

@@ -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
-    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):

watertap/flowsheets/full_water_resource_recovery_facility/BSM2_P_extension.py

@@ -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)
 
@@ -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)
@@ -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)
@@ -526,6 +532,14 @@ def set_operating_conditions(m):
     m.fs.thickener.hydraulic_retention_time.fix(86400 * pyo.units.s)
     m.fs.thickener.diameter.fix(10 * pyo.units.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
+
     def scale_variables(m):
         for var in m.fs.component_data_objects(pyo.Var, descend_into=True):
             if "flow_vol" in var.name:
@@ -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
@@ -872,10 +984,74 @@ def display_performance_metrics(m):
         pyo.units.get_units(m.fs.AD.liquid_phase.properties_in[0].flow_vol),
     )
 
+    print("---- Feed Metrics----")
+    print(
+        "Feed TSS concentration",
+        pyo.value(m.fs.FeedWater.properties[0].TSS),
+        pyo.units.get_units(m.fs.FeedWater.properties[0].TSS),
+    )
+    print(
+        "Feed COD concentration",
+        pyo.value(m.fs.FeedWater.properties[0].COD),
+        pyo.units.get_units(m.fs.FeedWater.properties[0].COD),
+    )
+    print(
+        "BOD5 concentration",
+        pyo.value(m.fs.FeedWater.properties[0].BOD5["effluent"]),
+        pyo.units.get_units(m.fs.FeedWater.properties[0].BOD5["effluent"]),
+    )
+    print(
+        "SNKj concentration",
+        pyo.value(m.fs.FeedWater.properties[0].SNKj),
+        pyo.units.get_units(m.fs.FeedWater.properties[0].SNKj),
+    )
+    print(
+        "Organic phosphorus concentration",
+        pyo.value(m.fs.FeedWater.properties[0].SP_organic),
+        pyo.units.get_units(m.fs.FeedWater.properties[0].SP_organic),
+    )
+    print(
+        "Inorganic phosphorus concentration",
+        pyo.value(m.fs.FeedWater.properties[0].SP_inorganic),
+        pyo.units.get_units(m.fs.FeedWater.properties[0].SP_inorganic),
+    )
+
+    print("---- Effluent Metrics----")
+    print(
+        "TSS concentration",
+        pyo.value(m.fs.Treated.properties[0].TSS),
+        pyo.units.get_units(m.fs.Treated.properties[0].TSS),
+    )
+    print(
+        "COD concentration",
+        pyo.value(m.fs.Treated.properties[0].COD),
+        pyo.units.get_units(m.fs.Treated.properties[0].COD),
+    )
+    print(
+        "BOD5 concentration",
+        pyo.value(m.fs.Treated.properties[0].BOD5["effluent"]),
+        pyo.units.get_units(m.fs.Treated.properties[0].BOD5["effluent"]),
+    )
+    print(
+        "SNKj concentration",
+        pyo.value(m.fs.Treated.properties[0].SNKj),
+        pyo.units.get_units(m.fs.Treated.properties[0].SNKj),
+    )
+    print(
+        "Organic phosphorus concentration",
+        pyo.value(m.fs.Treated.properties[0].SP_organic),
+        pyo.units.get_units(m.fs.Treated.properties[0].SP_organic),
+    )
+    print(
+        "Inorganic phosphorus concentration",
+        pyo.value(m.fs.Treated.properties[0].SP_inorganic),
+        pyo.units.get_units(m.fs.Treated.properties[0].SP_inorganic),
+    )
+
 
 if __name__ == "__main__":
     # This method builds and runs a steady state activated sludge flowsheet.
-    m, results = main(bio_P=False)
+    m, results = main(bio_P=True)
 
     stream_table = create_stream_table_dataframe(
         {

watertap/property_models/unit_specific/activated_sludge/asm1_properties.py

@@ -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,