Costing Unification: move from factor_total_investment to TIC (wa…

…tertap-org#1175)

* resolving TIC / factor_total_investment difference

* fixing idiosyncrasy with caoh2 unit cost

* standardizing cost_factor; adding aggregate_direct_captial_cost
Expression

* fix CSTR costing test

* make test_lssro_paper_analysis more robust when model is infeasible

* fix bad oaro_multi test

* reconciling ion exchange costing

* setting custom options to stabilize flowsheet_softening_two_stage

* make LCOW an Expression

* fixing issue with multiple choice costing block

---------

Co-authored-by: Bernard Knueven <bknueven@el2.ib0.cm.hpc.nrel.gov>

tutorials/unit_model_customization_example.ipynb

@@ -708,7 +708,7 @@
     "    pressure.append(m.fs.RO.inlet.pressure[0].value/1e5)\n",
     "    sh_avg=np.average([m.fs.RO.feed_side.N_Sh_comp[t,x,j].value for (t,x,j) in m.fs.RO.feed_side.N_Sh_comp])\n",
     "    avg_sh.append(sh_avg)\n",
-    "    lcow.append(m.fs.costing.LCOW.value)\n",
+    "    lcow.append(value(m.fs.costing.LCOW))\n",
     "    print(\"Solved multiplier {}\".format(adj))"
    ]
   },
@@ -818,15 +818,15 @@
     "    result = solver.solve(m, tee=False)\n",
     "    assert_optimal_termination(result)\n",
     "    pressures.append(m.fs.RO.inlet.pressure[0].value/1e5)    \n",
-    "    lcow_fixed_mem_cost.append(m.fs.costing.LCOW.value)\n",
+    "    lcow_fixed_mem_cost.append(value(m.fs.costing.LCOW))\n",
     "    mem_cost_fixed.append(m.fs.costing.reverse_osmosis.membrane_cost.value)\n",
     "    # second unfix our membrane cost and activate variable cost constraint and solve\n",
     "    m.fs.costing.reverse_osmosis.membrane_cost.unfix()\n",
     "    m.fs.RO_cost_pressure_constraint.activate()\n",
     "    assert degrees_of_freedom(m) == 0\n",
     "    result = solver.solve(m, tee=False)\n",
     "    assert_optimal_termination(result)\n",
-    "    lcow_variable_mem_cost.append(m.fs.costing.LCOW.value)\n",
+    "    lcow_variable_mem_cost.append(value(m.fs.costing.LCOW))\n",
     "    mem_cost_variable.append(m.fs.costing.reverse_osmosis.membrane_cost.value)\n",
     "\n",
     "    print(\"Solved con {}\".format(con))"

watertap/core/zero_order_base.py

@@ -329,22 +329,6 @@ def _get_performance_contents(self, time_point=0):
 
     # -------------------------------------------------------------------------
     # Unit operation costing methods
-    @staticmethod
-    def _add_cost_factor(blk, factor):
-        if factor == "TPEC":
-            blk.cost_factor = pyo.Expression(
-                expr=blk.config.flowsheet_costing_block.TPEC
-            )
-        elif factor == "TIC":
-            blk.cost_factor = pyo.Expression(
-                expr=blk.config.flowsheet_costing_block.TIC
-            )
-        else:
-            blk.cost_factor = pyo.Expression(expr=1.0)
-        blk.direct_capital_cost = pyo.Expression(
-            expr=blk.capital_cost / blk.cost_factor
-        )
-
     @staticmethod
     def _get_unit_cost_method(blk):
         """
@@ -462,7 +446,7 @@ def _general_power_law_form(
 
         expr *= number_of_parallel_units
 
-        blk.unit_model._add_cost_factor(blk, factor)
+        blk.costing_package.add_cost_factor(blk, factor)
 
         blk.capital_cost_constraint = pyo.Constraint(
             expr=blk.capital_cost == blk.cost_factor * expr

watertap/costing/costing_base.py

@@ -37,10 +37,6 @@ class WaterTAPCostingBlockData(FlowsheetCostingBlockData):
         CSTR: cost_cstr,
     }
 
-    def build(self):
-        super().build()
-        self._registered_LCOWs = {}
-
     def add_LCOW(self, flow_rate, name="LCOW"):
         """
         Add Levelized Cost of Water (LCOW) to costing block.
@@ -50,29 +46,22 @@ def add_LCOW(self, flow_rate, name="LCOW"):
             name (optional) - name for the LCOW variable (default: LCOW)
         """
 
-        LCOW = pyo.Var(
-            doc=f"Levelized Cost of Water based on flow {flow_rate.name}",
-            units=self.base_currency / pyo.units.m**3,
-        )
-        self.add_component(name, LCOW)
-
-        LCOW_constraint = pyo.Constraint(
-            expr=LCOW
-            == (
-                self.total_capital_cost * self.capital_recovery_factor
-                + self.total_operating_cost
-            )
-            / (
-                pyo.units.convert(
-                    flow_rate, to_units=pyo.units.m**3 / self.base_period
+        self.add_component(
+            name,
+            pyo.Expression(
+                expr=(
+                    self.total_capital_cost * self.capital_recovery_factor
+                    + self.total_operating_cost
                 )
-                * self.utilization_factor
+                / (
+                    pyo.units.convert(
+                        flow_rate, to_units=pyo.units.m**3 / self.base_period
+                    )
+                    * self.utilization_factor
+                ),
+                doc=f"Levelized Cost of Water based on flow {flow_rate.name}",
             ),
-            doc=f"Constraint for Levelized Cost of Water based on flow {flow_rate.name}",
         )
-        self.add_component(name + "_constraint", LCOW_constraint)
-
-        self._registered_LCOWs[name] = (LCOW, LCOW_constraint)
 
     def add_specific_energy_consumption(
         self, flow_rate, name="specific_energy_consumption"
@@ -206,19 +195,31 @@ def _build_common_global_params(self):
         )
 
         self.TPEC = pyo.Var(
-            initialize=3.4,
+            initialize=3.4 * (2.0 / 1.65),
             doc="Total Purchased Equipment Cost (TPEC)",
             units=pyo.units.dimensionless,
         )
 
         self.TIC = pyo.Var(
-            initialize=1.65,
+            initialize=2.0,
             doc="Total Installed Cost (TIC)",
             units=pyo.units.dimensionless,
         )
 
         self.fix_all_vars()
 
+    @staticmethod
+    def add_cost_factor(blk, factor):
+        if factor == "TPEC":
+            blk.cost_factor = pyo.Expression(expr=blk.costing_package.TPEC)
+        elif factor == "TIC":
+            blk.cost_factor = pyo.Expression(expr=blk.costing_package.TIC)
+        else:
+            blk.cost_factor = pyo.Expression(expr=1.0)
+        blk.direct_capital_cost = pyo.Expression(
+            expr=blk.capital_cost / blk.cost_factor
+        )
+
     def _get_costing_method_for(self, unit_model):
         """
         Allow the unit model to register its default costing method,
@@ -229,6 +230,44 @@ def _get_costing_method_for(self, unit_model):
             return unit_model.default_costing_method
         return super()._get_costing_method_for(unit_model)
 
+    def aggregate_costs(self):
+        """
+        This method aggregates costs from all the unit models and flows
+        registered with this FlowsheetCostingBlock and creates aggregate
+        variables for these on the FlowsheetCostingBlock that can be used for
+        further process-wide costing calculations.
+
+        The following costing variables are aggregated from all the registered
+        UnitModelCostingBlocks (if they exist):
+
+        * capital_cost,
+        * direct_capital_cost,
+        * fixed_operating_cost, and
+        * variable_operating_cost
+
+        Additionally, aggregate flow variables are created for all registered
+        flow types along with aggregate costs associated with each of these.
+
+        Args:
+            None
+        """
+        super().aggregate_costs()
+        c_units = self.base_currency
+
+        @self.Expression(doc="Aggregation Expression for direct capital cost")
+        def aggregate_direct_capital_cost(blk):
+            e = 0
+            for u in self._registered_unit_costing:
+                # Allow for units that might only have a subset of cost Vars
+                if hasattr(u, "direct_capital_cost"):
+                    e += pyo.units.convert(u.direct_capital_cost, to_units=c_units)
+                elif hasattr(u, "capital_cost"):
+                    raise RuntimeError(
+                        f"WaterTAP models with a capital_cost must also supply a direct_capital_cost. Found unit {u.unit_model} with `capital_cost` but no `direct_capital_cost`."
+                    )
+
+            return e
+
     def register_flow_type(self, flow_type, cost):
         """
         This method allows users to register new material and utility flows

watertap/costing/tests/test_costing_interop.py

@@ -96,15 +96,18 @@ def add_wt_costing(m):
 def test_hrcs_case_1575_wtcosting():
     hrcs.add_costing = add_wt_costing
 
-    m, results = simple_main()
-    pyo.assert_optimal_termination(results)
+    try:
+        m, results = simple_main()
+        pyo.assert_optimal_termination(results)
 
-    # check costing -- baseline is 0.02003276
-    assert pyo.value(m.fs.costing.LCOW) == pytest.approx(
-        0.02605311, rel=1e-3
-    )  # in $/m**3
-
-    hrcs.add_costing = hrcs_original_costing
+        # check costing -- baseline is 0.02003276
+        assert pyo.value(m.fs.costing.LCOW) == pytest.approx(
+            0.02419106, rel=1e-3
+        )  # in $/m**3
+    except:
+        raise
+    finally:
+        hrcs.add_costing = hrcs_original_costing
 
 
 def add_zo_costing(m):
@@ -167,12 +170,15 @@ def add_zo_costing(m):
 def test_hrcs_case_1575_zocosting():
     hrcs.add_costing = add_zo_costing
 
-    m, results = simple_main()
-    pyo.assert_optimal_termination(results)
-
-    # check costing -- baseline is 0.02003276
-    assert pyo.value(m.fs.costing.LCOW) == pytest.approx(
-        0.02087999, rel=1e-3
-    )  # in $/m**3
-
-    hrcs.add_costing = hrcs_original_costing
+    try:
+        m, results = simple_main()
+        pyo.assert_optimal_termination(results)
+
+        # check costing -- baseline is 0.02003276
+        assert pyo.value(m.fs.costing.LCOW) == pytest.approx(
+            0.02172723, rel=1e-3
+        )  # in $/m**3
+    except:
+        raise
+    finally:
+        hrcs.add_costing = hrcs_original_costing

watertap/costing/tests/test_multiple_choice_costing_block.py

@@ -238,7 +238,7 @@ def test_multiple_choice_costing_block():
         m.fs.RO2.costing.costing_blocks["high_pressure"].capital_cost
     )
 
-    assert m.fs.costing.total_capital_cost.value == 2 * (
+    assert m.fs.costing.total_capital_cost.value == (
         m.fs.RO.costing.costing_blocks["normal_pressure"].capital_cost.value
         + m.fs.RO2.costing.costing_blocks["high_pressure"].capital_cost.value
     )
@@ -250,21 +250,21 @@ def test_multiple_choice_costing_block():
     )
 
     # need to re-initialize
-    assert m.fs.costing.total_capital_cost.value == 2 * (
+    assert m.fs.costing.total_capital_cost.value == (
         m.fs.RO.costing.costing_blocks["normal_pressure"].capital_cost.value
         + m.fs.RO2.costing.costing_blocks["high_pressure"].capital_cost.value
     )
 
     m.fs.costing.initialize()
-    assert m.fs.costing.total_capital_cost.value == 2 * (
+    assert m.fs.costing.total_capital_cost.value == (
         m.fs.RO.costing.costing_blocks["high_pressure"].capital_cost.value
         + m.fs.RO2.costing.costing_blocks["high_pressure"].capital_cost.value
     )
     assert m.fs.costing.aggregate_variable_operating_cost.value == 42
 
     m.fs.RO3.costing.select_costing_block("high_pressure")
     m.fs.costing.initialize()
-    assert m.fs.costing.total_capital_cost.value == 2 * (
+    assert m.fs.costing.total_capital_cost.value == (
         m.fs.RO.costing.costing_blocks["high_pressure"].capital_cost.value
         + m.fs.RO2.costing.costing_blocks["high_pressure"].capital_cost.value
         + m.fs.RO3.costing.costing_blocks["high_pressure"].capital_cost.value

watertap/costing/tests/test_util.py

@@ -136,6 +136,7 @@ def build_dummy_cost_rectifier(blk):
 )
 def dummy_cost_rectifier(blk):
     cost_rectifier(blk)
+    blk.costing_package.add_cost_factor(blk, "TIC")
     blk.capital_cost_constraint = pyo.Constraint(
         expr=blk.capital_cost == blk.capital_cost_rectifier
     )

watertap/costing/tests/test_zero_order_costing.py

@@ -426,8 +426,7 @@ def test_process_costing(self, model):
     def test_add_LCOW(self, model):
         model.fs.costing.add_LCOW(model.fs.unit1.properties_in[0].flow_vol)
 
-        assert isinstance(model.fs.costing.LCOW, Var)
-        assert isinstance(model.fs.costing.LCOW_constraint, Constraint)
+        assert isinstance(model.fs.costing.LCOW, Expression)
 
         assert_units_consistent(model.fs)
         assert degrees_of_freedom(model.fs) == 0

watertap/costing/unit_models/anaerobic_digestor.py

@@ -83,10 +83,11 @@ def cost_anaerobic_digestor_capital(
         flow_in / blk.reference_flow, to_units=pyo.units.dimensionless
     )
 
-    print(f"base_currency: {blk.costing_package.base_currency}")
+    blk.costing_package.add_cost_factor(blk, "TIC")
     blk.capital_cost_constraint = pyo.Constraint(
         expr=blk.capital_cost
-        == pyo.units.convert(
+        == blk.cost_factor
+        * pyo.units.convert(
             blk.capital_a_parameter * sizing_term**blk.capital_b_parameter,
             to_units=blk.costing_package.base_currency,
         )

watertap/costing/unit_models/clarifier.py

@@ -84,14 +84,16 @@ def cost_circular_clarifier(blk, cost_electricity_flow=True):
     Circular clarifier costing method [1]
     """
     make_capital_cost_var(blk)
+    blk.costing_package.add_cost_factor(blk, "TIC")
 
     surface_area = pyo.units.convert(
         blk.unit_model.surface_area, to_units=pyo.units.ft**2
     )
 
     blk.capital_cost_constraint = pyo.Constraint(
         expr=blk.capital_cost
-        == pyo.units.convert(
+        == blk.cost_factor
+        * pyo.units.convert(
             blk.costing_package.circular.construction_a_parameter * surface_area**2
             + blk.costing_package.circular.construction_b_parameter * surface_area
             + blk.costing_package.circular.construction_c_parameter,
@@ -140,14 +142,16 @@ def cost_rectangular_clarifier(blk, cost_electricity_flow=True):
     Rectangular clarifier costing method [1]
     """
     make_capital_cost_var(blk)
+    blk.costing_package.add_cost_factor(blk, "TIC")
 
     surface_area = pyo.units.convert(
         blk.unit_model.surface_area, to_units=pyo.units.ft**2
     )
 
     blk.capital_cost_constraint = pyo.Constraint(
         expr=blk.capital_cost
-        == pyo.units.convert(
+        == blk.cost_factor
+        * pyo.units.convert(
             blk.costing_package.rectangular.construction_a_parameter * surface_area**2
             + blk.costing_package.rectangular.construction_b_parameter * surface_area
             + blk.costing_package.rectangular.construction_c_parameter,
@@ -190,14 +194,16 @@ def cost_primary_clarifier(blk, cost_electricity_flow=True):
     Primary clarifier costing method [2]
     """
     make_capital_cost_var(blk)
+    blk.costing_package.add_cost_factor(blk, "TIC")
 
     t0 = blk.flowsheet().time.first()
     flow_in = pyo.units.convert(
         blk.unit_model.inlet.flow_vol[t0], to_units=pyo.units.gallon / pyo.units.day
     )
     blk.capital_cost_constraint = pyo.Constraint(
         expr=blk.capital_cost
-        == pyo.units.convert(
+        == blk.cost_factor
+        * pyo.units.convert(
             blk.costing_package.primary.capital_a_parameter
             * pyo.units.convert(
                 flow_in / (1e6 * pyo.units.gallon / pyo.units.day),

watertap/costing/unit_models/compressor.py

@@ -34,9 +34,11 @@ def build_compressor_cost_param_block(blk):
 )
 def cost_compressor(blk, cost_electricity_flow=True):
     make_capital_cost_var(blk)
+    blk.costing_package.add_cost_factor(blk, "TIC")
     blk.capital_cost_constraint = pyo.Constraint(
         expr=blk.capital_cost
-        == pyo.units.convert(
+        == blk.cost_factor
+        * pyo.units.convert(
             blk.costing_package.compressor.unit_cost
             * blk.unit_model.control_volume.properties_in[0].flow_mass_phase_comp[
                 "Vap", "H2O"