Merge branch 'main' into pr/zacharybinger/152

src/watertap_contrib/reflo/analysis/case_studies/KBHDP/KBHDP_ZLD.py

@@ -276,6 +276,7 @@ def add_treatment_costing(m):
     treatment.pump.costing = UnitModelCostingBlock(
         flowsheet_costing_block=treatment.costing,
     )
+
     add_ec_costing(m, treatment.EC, treatment.costing)
     add_UF_costing(m, treatment.UF, treatment.costing)
     add_ro_costing(m, treatment.RO, treatment.costing)
@@ -469,7 +470,7 @@ def display_unfixed_vars(blk, report=True):
             print(f"\t{v2.name:<40s}")
 
 
-def set_operating_conditions(m, RO_pressure=15e5):
+def set_operating_conditions(m, RO_pressure=20e5):
     treatment = m.fs.treatment
     pump_efi = 0.8  # pump efficiency [-]
     # Set inlet conditions and operating conditions for each unit

src/watertap_contrib/reflo/analysis/case_studies/KBHDP/components/EC.py

@@ -15,7 +15,7 @@
     assert_optimal_termination,
     units as pyunits,
 )
-
+from pyomo.util.calc_var_value import calculate_variable_from_constraint as cvc
 from idaes.core import FlowsheetBlock, UnitModelCostingBlock
 from idaes.core.solvers import get_solver
 
@@ -39,6 +39,20 @@
 import numpy as np
 from watertap.costing import WaterTAPCosting
 import math
+from watertap_contrib.reflo.analysis.case_studies.KBHDP.utils import (
+    check_jac,
+    calc_scale,
+)
+
+__all__ = [
+    "build_ec",
+    "set_ec_operating_conditions",
+    "init_ec",
+    "add_ec_costing",
+    "add_ec_scaling",
+    "report_EC",
+    "print_EC_costing_breakdown",
+]
 
 __all__ = [
     "build_ec",
@@ -157,12 +171,21 @@ def set_system_operating_conditions(m):
     # # initialize feed
 
 
-def set_ec_operating_conditions(m, blk):
+def set_ec_operating_conditions(m, blk, conv=5e3):
     """Set EC operating conditions"""
     # Check if the set up of the ec inputs is correct
     print(f"EC Degrees of Freedom: {degrees_of_freedom(blk.ec)}")
 
     blk.ec.load_parameters_from_database(use_default_removal=True)
+    # blk.ec.conductivity.unfix()
+    # tds_ec_conversion = conv * (pyunits.mg * pyunits.m) / (pyunits.liter * pyunits.S)
+    # blk.ec.conductivity_constr = Constraint(
+    #     expr=blk.ec.conductivity
+    #     == pyunits.convert(
+    #         blk.feed.properties[0].conc_mass_comp["tds"] / tds_ec_conversion,
+    #         to_units=pyunits.S / pyunits.m,
+    #     )
+    # )
     # blk.feed.properties[0.0].flow_mass_comp["tss"].fix(5.22e-6)
     # blk.ec.overpotential.fix(2)
     print(f"EC Degrees of Freedom: {degrees_of_freedom(blk.ec)}")
@@ -175,24 +198,28 @@ def calc_scale(value):
 
     scale_flow = calc_scale(m.fs.feed.flow_mass_comp[0, "H2O"].value)
     scale_tds = calc_scale(m.fs.feed.flow_mass_comp[0, "tds"].value)
+    scale_tss = calc_scale(m.fs.feed.flow_mass_comp[0, "tss"].value)
 
     m.fs.properties.set_default_scaling(
         "flow_mass_comp", 10**-scale_flow, index=("H2O")
     )
     m.fs.properties.set_default_scaling("flow_mass_comp", 10**-scale_tds, index=("tds"))
+    m.fs.properties.set_default_scaling("flow_mass_comp", 10**-scale_tss, index=("tss"))
     calculate_scaling_factors(m)
 
-    badly_scaled_var_list = list_badly_scaled_variables(m)
-    if len(badly_scaled_var_list) > 0:
-        [print(i[0].name, i[1]) for i in badly_scaled_var_list]
-    else:
-        print("Variables are scaled well")
-
 
 def add_ec_scaling(m, blk):
-    set_scaling_factor(blk.ec.charge_loading_rate, 1e-3)
-    set_scaling_factor(blk.ec.reactor_volume, 1e-1)
-    set_scaling_factor(blk.ec.power_required, 1e-6)
+    set_scaling_factor(blk.ec.charge_loading_rate, 1e1)
+    set_scaling_factor(blk.ec.electrode_volume, 1)
+    # set_scaling_factor(blk.ec.cell_voltage, )
+    # set_scaling_factor(blk.ec.applied_current, 1e5)
+    # set_scaling_factor(blk.ec.power_required, 1e-6)
+
+    constraint_scaling_transform(blk.ec.eq_power_required, 1e-4)
+    set_scaling_factor(blk.ec.properties_in[0.0].flow_mass_comp["H2O"], 1e1)
+    set_scaling_factor(blk.ec.properties_in[0.0].flow_mass_comp["tds"], 1)
+    set_scaling_factor(blk.ec.properties_treated[0.0].flow_mass_comp["H2O"], 1e1)
+    set_scaling_factor(blk.ec.properties_byproduct[0.0].flow_mass_comp["H2O"], 1)
 
 
 def init_system(m, solver=None):
@@ -222,19 +249,41 @@ def init_ec(m, blk, solver=None):
         solver = get_solver()
 
     optarg = solver.options
-    # assert_no_degrees_of_freedom(m)
-    _initialize(blk.feed)
+
+    blk.feed.initialize(optarg=optarg)
     propagate_state(blk.feed_to_ec)
 
-    _initialize(blk.ec)
+    cvc(blk.ec.overpotential, blk.ec.eq_overpotential)
+    cvc(blk.ec.applied_current, blk.ec.eq_applied_current)
+    cvc(blk.ec.anode_area, blk.ec.eq_electrode_area_total)
+    cvc(blk.ec.ohmic_resistance, blk.ec.eq_ohmic_resistance)
 
+    blk.ec.initialize(optarg=optarg)
     propagate_state(blk.ec_to_product)
     propagate_state(blk.ec_to_disposal)
 
 
+# def init_ec(m, blk, solver=None):
+#     """Initialize IX model"""
+
+#     if solver is None:
+#         solver = get_solver()
+
+#     optarg = solver.options
+#     # assert_no_degrees_of_freedom(m)
+#     _initialize(blk.feed)
+#     propagate_state(blk.feed_to_ec)
+
+#     _initialize(blk.ec)
+
+#     propagate_state(blk.ec_to_product)
+#     propagate_state(blk.ec_to_disposal)
+
+
 def add_system_costing(m):
     """Add system level costing components"""
     m.fs.costing = ZeroOrderCosting()
+    m.fs.costing.base_currency = pyunits.USD_2022
     add_ec_costing(m, m.fs.EC)
     calc_costing(m, m.fs.EC)
 
@@ -254,6 +303,44 @@ def calc_costing(m, blk):
     m.fs.costing.add_electricity_intensity(blk.ec.properties_treated[0].flow_vol)
 
 
+def solve(m, solver=None, tee=True, raise_on_failure=True, debug=False):
+    # ---solving---
+    if solver is None:
+        solver = get_solver()
+        solver.options["max_iter"] = 2000
+
+    print("\n--------- SOLVING ---------\n")
+
+    results = solver.solve(m, tee=tee)
+
+    if check_optimal_termination(results):
+        print("\n--------- OPTIMAL SOLVE!!! ---------\n")
+        if debug:
+            print("\n--------- CHECKING JACOBIAN ---------\n")
+            print("\n--------- TREATMENT ---------\n")
+            check_jac(m)
+
+            print("\n--------- CLOSE TO BOUNDS ---------\n")
+            print_close_to_bounds(m)
+        return results
+    msg = (
+        "The current configuration is infeasible. Please adjust the decision variables."
+    )
+    if raise_on_failure:
+        print('\n{"=======> INFEASIBLE BOUNDS <=======":^60}\n')
+        print_infeasible_bounds(m)
+        print('\n{"=======> INFEASIBLE CONSTRAINTS <=======":^60}\n')
+        print_infeasible_constraints(m)
+        print('\n{"=======> CLOSE TO BOUNDS <=======":^60}\n')
+        print_close_to_bounds(m)
+
+        raise RuntimeError(msg)
+    else:
+        print("\n--------- FAILED SOLVE!!! ---------\n")
+        print(msg)
+        assert False
+
+
 def report_EC(blk):
 
     print(f"\n\n-------------------- EC Report --------------------\n")
@@ -319,16 +406,15 @@ def breakdown_dof(blk):
     m = build_system()
     set_system_operating_conditions(m)
     set_ec_operating_conditions(m, m.fs.EC)
-    set_scaling(m, m.fs.EC)
     add_ec_scaling(m, m.fs.EC)
+    set_scaling(m, m.fs.EC)
     init_system(m)
     add_system_costing(m)
 
     solver = get_solver()
     m.fs.objective_lcow = Objective(expr=m.fs.costing.LCOW)
-    results = solver.solve(m)
+    results = solve(m, debug=True)
 
-    print(m.fs.objective_lcow())
     report_EC(m.fs.EC)
     print_EC_costing_breakdown(m.fs.EC)
-    # print(m.fs.EC.ec.display())
+    m.fs.EC.ec.conductivity.display()