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

Swap LCOW LCOT in REFLOSystemCosting #168

Open
wants to merge 9 commits into
base: main
Choose a base branch
from
Original file line number Diff line number Diff line change
Expand Up @@ -192,7 +192,6 @@ def set_ec_operating_conditions(m, blk, conv=5e3):


def set_scaling(m, blk):

def calc_scale(value):
return math.floor(math.log(value, 10))

Expand Down
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
#################################################################################
# WaterTAP Copyright (c) 2020-2024, The Regents of the University of California,
# WaterTAP Copyright (c) 2020-2025, The Regents of the University of California,
# through Lawrence Berkeley National Laboratory, Oak Ridge National Laboratory,
# National Renewable Energy Laboratory, and National Energy Technology
# Laboratory (subject to receipt of any required approvals from the U.S. Dept.
Expand Down Expand Up @@ -492,7 +492,6 @@ def test_default_build(self, default_build):


class TestElectricityGenOnlyWithHeat:

@pytest.fixture(scope="class")
def energy_gen_only_with_heat(self):

Expand Down Expand Up @@ -629,7 +628,6 @@ def test_optimize_frac_from_grid(self):


class TestElectricityGenOnlyNoHeat:

@pytest.fixture(scope="class")
def elec_gen_only_no_heat(self):

Expand Down Expand Up @@ -899,7 +897,6 @@ def test_optimize_frac_from_grid(self):


class TestElectricityAndHeatGen:

@pytest.fixture(scope="class")
def heat_and_elec_gen(self):

Expand Down Expand Up @@ -948,7 +945,7 @@ def test_build(self, heat_and_elec_gen):
assert hasattr(m.fs.energy.costing, "LCOE")
assert hasattr(m.fs.energy.costing, "LCOH")
# treatment metrics on TreatmentCosting
assert not hasattr(m.fs.treatment.costing, "LCOT")
assert hasattr(m.fs.treatment.costing, "LCOT")
assert hasattr(m.fs.treatment.costing, "SEEC")
assert hasattr(m.fs.treatment.costing, "STEC")

Expand Down Expand Up @@ -1094,6 +1091,116 @@ def test_no_energy_treatment_block():
m.fs.costing = REFLOSystemCosting()


@pytest.mark.component
def test_add_LCOW_and_LCOT_to_treatment_and_system_costing():

m = ConcreteModel()
m.fs = FlowsheetBlock(dynamic=False)
m.fs.properties = SeawaterParameterBlock()

#### TREATMENT BLOCK
m.fs.treatment = Block()
m.fs.treatment.costing = TreatmentCosting()
m.fs.treatment.costing.base_currency = pyunits.USD_2011

m.fs.treatment.unit = DummyTreatmentUnit(property_package=m.fs.properties)
m.fs.treatment.unit.costing = UnitModelCostingBlock(
flowsheet_costing_block=m.fs.treatment.costing
)

m.fs.treatment.unit.design_var_a.fix()
m.fs.treatment.unit.design_var_b.fix()
m.fs.treatment.unit.electricity_consumption.fix(110)
m.fs.treatment.unit.heat_consumption.fix(250)

m.fs.treatment.costing.cost_process()
# add_LCOW and add_LCOT to TreatmentCosting before adding to REFLOSystemCosting.
# At this point, LCOW and LCOT are parameters on TreatmentCosting.
m.fs.treatment.costing.add_LCOW(
m.fs.treatment.unit.properties[0].flow_vol_phase["Liq"]
)
m.fs.treatment.costing.add_LCOT(
m.fs.treatment.unit.properties[0].flow_vol_phase["Liq"]
)
LCOT_id = id(m.fs.treatment.costing.LCOT)
LCOW_id = id(m.fs.treatment.costing.LCOW)

#### ENERGY BLOCK
m.fs.energy = Block()
m.fs.energy.costing = EnergyCosting()
m.fs.energy.costing.base_currency = pyunits.USD_2011
m.fs.energy.heat_unit = DummyHeatUnit()
m.fs.energy.elec_unit = DummyElectricityUnit()
m.fs.energy.heat_unit.costing = UnitModelCostingBlock(
flowsheet_costing_block=m.fs.energy.costing
)
m.fs.energy.elec_unit.costing = UnitModelCostingBlock(
flowsheet_costing_block=m.fs.energy.costing
)
m.fs.energy.heat_unit.heat.fix(50)
m.fs.energy.elec_unit.electricity.fix(100)
m.fs.energy.costing.cost_process()

#### SYSTEM COSTING
m.fs.costing = REFLOSystemCosting()
m.fs.costing.cost_process()

#### SCALING
m.fs.properties.set_default_scaling(
"flow_mass_phase_comp", 1e-1, index=("Liq", "H2O")
)
m.fs.properties.set_default_scaling(
"flow_mass_phase_comp", 1e-1, index=("Liq", "TDS")
)
calculate_scaling_factors(m)

#### INITIALIZE
m.fs.treatment.unit.properties.calculate_state(
var_args={
("flow_vol_phase", "Liq"): 0.4381,
("conc_mass_phase_comp", ("Liq", "TDS")): 20,
("temperature", None): 293,
("pressure", None): 101325,
},
hold_state=True,
)

assert degrees_of_freedom(m) == 0

m.fs.treatment.unit.initialize()
m.fs.treatment.costing.initialize()
m.fs.energy.costing.initialize()
m.fs.costing.initialize()
m.fs.costing.add_LCOE()
m.fs.costing.add_LCOH()
m.fs.costing.add_specific_electric_energy_consumption(
m.fs.treatment.unit.properties[0].flow_vol_phase["Liq"]
)
m.fs.costing.add_specific_thermal_energy_consumption(
m.fs.treatment.unit.properties[0].flow_vol_phase["Liq"]
)
m.fs.costing.add_LCOW(m.fs.treatment.unit.properties[0].flow_vol_phase["Liq"])
m.fs.costing.add_LCOT(m.fs.treatment.unit.properties[0].flow_vol_phase["Liq"])

results = solver.solve(m)
assert_optimal_termination(results)

# check that we have created new objects and
# not just references to LCOW and LCOT
assert id(m.fs.costing.LCOT) != LCOT_id
assert id(m.fs.costing.LCOW) != LCOW_id
# check that LCOW doesn't exist anymore on treatment costing
assert not hasattr(m.fs.treatment.costing, "LCOW")
assert hasattr(m.fs.treatment.costing, "_LCOW_flow_rate")
assert (
m.fs.treatment.costing._LCOW_flow_rate
is m.fs.treatment.unit.properties[0].flow_vol_phase["Liq"]
)
assert m.fs.costing.LCOT is m.fs.treatment.costing.LCOT
assert m.fs.costing.LCOE is m.fs.energy.costing.LCOE
assert m.fs.costing.LCOH is m.fs.energy.costing.LCOH


@pytest.mark.component
def test_add_LCOW_to_energy_costing():

Expand Down
3 changes: 2 additions & 1 deletion src/watertap_contrib/reflo/costing/units/lt_med_surrogate.py
Original file line number Diff line number Diff line change
Expand Up @@ -159,7 +159,8 @@ def cost_lt_med_surrogate(blk):
dist.flow_vol_phase["Liq"], to_units=pyo.units.m**3 / pyo.units.day
)
blk.capacity_dimensionless = pyo.units.convert(
blk.capacity * pyo.units.day * pyo.units.m**-3, to_units=pyo.units.dimensionless
blk.capacity * pyo.units.day * pyo.units.m**-3,
to_units=pyo.units.dimensionless,
)

blk.annual_dist_production = pyo.units.convert(
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -150,7 +150,8 @@ def cost_med_tvc_surrogate(blk):
dist.flow_vol_phase["Liq"], to_units=pyo.units.m**3 / pyo.units.day
)
blk.capacity_dimensionless = pyo.units.convert(
blk.capacity * pyo.units.day * pyo.units.m**-3, to_units=pyo.units.dimensionless
blk.capacity * pyo.units.day * pyo.units.m**-3,
to_units=pyo.units.dimensionless,
)

blk.annual_dist_production = pyo.units.convert(
Expand Down
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
#################################################################################
# WaterTAP Copyright (c) 2020-2024, The Regents of the University of California,
# WaterTAP Copyright (c) 2020-2025, The Regents of the University of California,
# through Lawrence Berkeley National Laboratory, Oak Ridge National Laboratory,
# National Renewable Energy Laboratory, and National Energy Technology
# Laboratory (subject to receipt of any required approvals from the U.S. Dept.
Expand Down Expand Up @@ -168,6 +168,18 @@ def add_specific_thermal_energy_consumption(
specific_thermal_energy_consumption_constraint,
)

def add_LCOW(self, flow_rate, name="LCOW"):
super().add_LCOW(flow_rate, name=name)
# NOTE: Add reference to flow rate here so that if a user tries to add_LCOW
# on a REFLOSystemCosting block after add_LCOW is called on TreatmentCosting,
# the parameter can be renamed from LCOW to LCOT and use the same flow rate.
add_object_reference(self, "_LCOW_flow_rate", flow_rate)

def add_LCOT(self, flow_rate):
Copy link
Contributor

Choose a reason for hiding this comment

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

Why do we need this separate function? You can just use add_LCOW with name of "LCOT" and point to relevant flowrate, right?


if not hasattr(self, "LCOT"):
self.add_LCOW(flow_rate, name="LCOT")


@declare_process_block_class("EnergyCosting")
class EnergyCostingData(REFLOCostingData):
Expand Down Expand Up @@ -233,7 +245,6 @@ def build_process_costs(self):
super().build_process_costs()

def build_LCOE_params(self):

def rule_yearly_electricity_production(b, y):
if y == 0:
return b.yearly_electricity_production[y] == pyo.units.convert(
Expand Down Expand Up @@ -269,7 +280,6 @@ def rule_lifetime_electricity_production(b):
set_scaling_factor(self.lifetime_electricity_production, 1e-4)

def build_LCOH_params(self):

def rule_yearly_heat_production(b, y):
if y == 0:
return b.yearly_heat_production[y] == pyo.units.convert(
Expand Down Expand Up @@ -728,10 +738,10 @@ def initialize_build(self):

super().initialize_build()

if hasattr(self, "LCOT"):
if hasattr(self, "LCOW"):
calculate_variable_from_constraint(
self.LCOT,
self.LCOT_constraint,
self.LCOW,
self.LCOW_constraint,
)

if hasattr(self, "LCOE"):
Expand Down Expand Up @@ -810,22 +820,37 @@ def build_process_costs(self):
"""
pass

def add_LCOT(self, flow_rate):
def add_LCOW(self, flow_rate):
"""
Add Levelized Cost of Treatment (LCOT) to costing block.
Add Levelized Cost of Water (LCOW) to costing block.
Args:
flow_rate - flow rate of water (volumetric) to be used in
calculating LCOT
calculating LCOW
"""
treat_cost = self._get_treatment_cost_block()

LCOT = pyo.Var(
doc=f"Levelized Cost of Treatment based on flow {flow_rate.name}",
if hasattr(treat_cost, "LCOW"):
# NOTE: If a user has called add_LCOW on a TreatmentCosting block
# prior to calling add_LCOW on REFLOSystemCosting, that parameter
# is recreated here using the same flow rate and renamed to 'LCOT'.
# This is done to ensure there aren't two parameters named 'LCOW' on
# TreatmentCosting and REFLOSystemCosting that are calculated with
# different parameters.
warning_msg = "When adding LCOW to REFLOSystemCosting, 'LCOW' was "
warning_msg += "found on the TreatmentCosting block and renamed 'LCOT' "
warning_msg += "to avoid ambiguity."
_log.warning(warning_msg)
treat_cost.del_component("LCOW")
treat_cost.add_LCOT(treat_cost._LCOW_flow_rate)
Comment on lines +842 to +844
Copy link
Contributor

Choose a reason for hiding this comment

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

I am not sure we should be doing this. Typically, we'd be more restrictive and just raise an error, telling the user that LCOW already exists and to specify a different name. Also not clear when just looking at the diff why we have a second add_LCOW definition.

Copy link
Contributor

Choose a reason for hiding this comment

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

All this monkey business I do here is to make a consideration for when the user add_LCOW using the treatment costing package (e.g., m.fs.treatment.costing.add_LCOW(flow_rate)) before doing add_LCOW with REFLOSystemCosting (e.g., m.fs.costing.add_LCOW(flow_rate)) . In that scenario, you would have two different identically named variables (LCOW) that are calculated differently.

This seemed problematic to me and potentially confusing to the user. Especially as the behavior we wanted to enforce is to have is swapping LCOW and LCOT (LCOT really just being a metric we sort of made up that is focused purely on the treatment side and to distinguish from LCOW, the levelized metric we are all familiar with). Thus, this is a way to make sure the metrics are where we want them to be. The _LCOW_flow_rate attribute is for the (probably extreme) edge case where the user wants the two levelized metrics calculated on the basis of different flow rates. (I may have really over coddled the user on this one.)


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("LCOT", LCOT)
self.add_component("LCOW", LCOW)

LCOT_constraint = pyo.Constraint(
expr=LCOT
LCOW_constraint = pyo.Constraint(
expr=LCOW
== (
self.total_capital_cost * self.capital_recovery_factor
+ self.total_operating_cost
Expand All @@ -834,9 +859,31 @@ def add_LCOT(self, flow_rate):
pyo.units.convert(flow_rate, to_units=pyo.units.m**3 / self.base_period)
* self.utilization_factor
),
doc=f"Constraint for Levelized Cost of Treatment based on flow {flow_rate.name}",
doc=f"Constraint for Levelized Cost of Water based on flow {flow_rate.name}",
)
self.add_component("LCOT_constraint", LCOT_constraint)
self.add_component("LCOW_constraint", LCOW_constraint)

def add_LCOT(self, flow_rate):
"""
Add Levelized Cost of Treatment (LCOT) to costing block.
"""

treat_cost = self._get_treatment_cost_block()

if hasattr(treat_cost, "LCOT"):
# NOTE: If LCOT exists already on TreatmentCosting, it is deleted to
# ensure the LCOT being added here uses the intended flow_rate.
warning_msg = "When adding LCOT to REFLOSystemCosting, 'LCOT' was "
warning_msg += "found on the TreatmentCosting block. That parameter was "
warning_msg += f"deleted and recreated using {flow_rate.name}."
_log.warning(warning_msg)
treat_cost.del_component("LCOT")
treat_cost.add_LCOT(flow_rate)

else:
treat_cost.add_LCOT(flow_rate)

add_object_reference(self, "LCOT", getattr(treat_cost, "LCOT"))

def add_LCOE(self):
"""
Expand All @@ -849,18 +896,6 @@ def add_LCOE(self):

add_object_reference(self, "LCOE", energy_cost.LCOE)

def add_LCOW(self, flow_rate, name="LCOW"):
"""
Add Levelized Cost of Water (LCOW) to costing block.
"""

treat_cost = self._get_treatment_cost_block()

if not hasattr(treat_cost, "LCOW"):
treat_cost.add_LCOW(flow_rate, name="LCOW")

add_object_reference(self, name, getattr(treat_cost, name))

def add_LCOH(self):
"""
Add Levelized Cost of Heat (LCOH) to costing block.
Expand Down
9 changes: 5 additions & 4 deletions src/watertap_contrib/reflo/unit_models/solar_still.py
Original file line number Diff line number Diff line change
Expand Up @@ -192,10 +192,11 @@ def build(self):
f" g/L TDS with {self.config.water_yield_calculation_args['initial_water_depth']} m initial water depth."
)

daily_water_yield_mass, num_zld_cycles_per_year = (
self.calculate_daily_water_yield(
**self.config.water_yield_calculation_args
)
(
daily_water_yield_mass,
num_zld_cycles_per_year,
) = self.calculate_daily_water_yield(
**self.config.water_yield_calculation_args
)

unit_log.info(
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -392,7 +392,6 @@ def test_costing(self, dwi_frame):


class TestDeepWellInjection_SimpleCosting:

@pytest.fixture(scope="class")
def dwi_frame(self):

Expand Down Expand Up @@ -492,7 +491,6 @@ def test_costing_as_opex(self, dwi_frame):


class TestDeepWellInjection_10000ft:

@pytest.fixture(scope="class")
def dwi_10000_frame(self):
m = build_dwi_10000()
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -82,7 +82,6 @@ def create_input_arrays(
temperature_col=None, # column from weather_data to use as temperature input data
wind_velocity_col=None, # column from weather_data to use as wind velocity input data
):

def generate_continuous_day_series():
# Days in each month for non-leap year
days_in_month = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]
Expand Down Expand Up @@ -163,9 +162,11 @@ def get_solar_still_daily_water_yield(
if len(blk.weather_data) > 8760:
blk.weather_data = blk.weather_data.loc[:8760]

blk.ambient_temp_by_hr, blk.irradiance_by_hr, blk.wind_vel_by_hr = (
create_input_arrays(blk, **kwargs)
)
(
blk.ambient_temp_by_hr,
blk.irradiance_by_hr,
blk.wind_vel_by_hr,
) = create_input_arrays(blk, **kwargs)

len_data_hr = len(blk.irradiance_by_hr)

Expand Down