diff --git a/CHANGELOG.md b/CHANGELOG.md index 1b8f207d5..155ff73fb 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -37,6 +37,7 @@ - Remove `pytest-subtests` as it's incorporated into pytest as of v9, and is an archived project. - Added [Ard](https://github.com/NLRWindSystems/Ard) as a combined performance and cost model - Added `PerformanceModelBaseClass` and standardized outputs of converter performance models +- Updated storage control strategies to output `commodity_set_point` and generic storage performance models to output `commodity_out` and input `commodity_set_point`. ## 0.5.1 [December 18, 2025] diff --git a/examples/test/test_all_examples.py b/examples/test/test_all_examples.py index 05ef94564..f55632960 100644 --- a/examples/test/test_all_examples.py +++ b/examples/test/test_all_examples.py @@ -258,7 +258,15 @@ def test_ammonia_synloop_example(subtests): assert pytest.approx(model.prob.get_val("ammonia.CapEx"), rel=1e-6) == 1.15173753e09 with subtests.test("Check ammonia OpEx"): - assert pytest.approx(model.prob.get_val("ammonia.OpEx"), rel=1e-4) == 25737370.661763854 + assert pytest.approx(model.prob.get_val("ammonia.OpEx")[0], rel=1e-4) == 25414748.989416014 + + with subtests.test("Check ammonia production"): + assert ( + pytest.approx( + model.prob.get_val("ammonia.annual_ammonia_produced", units="t/yr").mean(), rel=1e-4 + ) + == 406333.161 + ) with subtests.test("Check total adjusted CapEx"): assert ( @@ -273,7 +281,7 @@ def test_ammonia_synloop_example(subtests): pytest.approx( model.prob.get_val("finance_subgroup_nh3.total_opex_adjusted")[0], rel=1e-6 ) - == 79744581.00552343 + == 79421959.33317558 ) with subtests.test("Check LCOH"): @@ -285,7 +293,7 @@ def test_ammonia_synloop_example(subtests): with subtests.test("Check LCOA"): assert ( pytest.approx(model.prob.get_val("finance_subgroup_nh3.LCOA")[0], rel=1e-6) - == 1.2310335361130984 + == 1.1022714567388747 ) diff --git a/h2integrate/control/control_strategies/converters/demand_openloop_controller.py b/h2integrate/control/control_strategies/converters/demand_openloop_controller.py index b9125e7ec..faa67d4a3 100644 --- a/h2integrate/control/control_strategies/converters/demand_openloop_controller.py +++ b/h2integrate/control/control_strategies/converters/demand_openloop_controller.py @@ -70,6 +70,6 @@ def compute(self, inputs, outputs): ) # Calculate actual output based on demand met and curtailment - outputs[f"{commodity}_out"] = ( + outputs[f"{commodity}_set_point"] = ( inputs[f"{commodity}_in"] - outputs[f"{commodity}_unused_commodity"] ) diff --git a/h2integrate/control/control_strategies/converters/flexible_demand_openloop_controller.py b/h2integrate/control/control_strategies/converters/flexible_demand_openloop_controller.py index 398c89f85..43d7923a5 100644 --- a/h2integrate/control/control_strategies/converters/flexible_demand_openloop_controller.py +++ b/h2integrate/control/control_strategies/converters/flexible_demand_openloop_controller.py @@ -297,6 +297,6 @@ def compute(self, inputs, outputs): ) # Calculate actual output based on demand met and curtailment - outputs[f"{commodity}_out"] = ( + outputs[f"{commodity}_set_point"] = ( inputs[f"{commodity}_in"] - outputs[f"{commodity}_unused_commodity"] ) diff --git a/h2integrate/control/control_strategies/demand_openloop_controller.py b/h2integrate/control/control_strategies/demand_openloop_controller.py index f0724d922..bbab59687 100644 --- a/h2integrate/control/control_strategies/demand_openloop_controller.py +++ b/h2integrate/control/control_strategies/demand_openloop_controller.py @@ -97,7 +97,7 @@ def setup(self): ) self.add_output( - f"{commodity}_out", + f"{commodity}_set_point", val=0.0, shape=(n_timesteps), units=self.config.commodity_units, diff --git a/h2integrate/control/control_strategies/passthrough_openloop_controller.py b/h2integrate/control/control_strategies/passthrough_openloop_controller.py index 9f3261449..9dafa6082 100644 --- a/h2integrate/control/control_strategies/passthrough_openloop_controller.py +++ b/h2integrate/control/control_strategies/passthrough_openloop_controller.py @@ -35,7 +35,7 @@ def setup(self): ) self.add_output( - f"{self.config.commodity_name}_out", + f"{self.config.commodity_name}_set_point", copy_shape=f"{self.config.commodity_name}_in", units=self.config.commodity_units, desc=f"{self.config.commodity_name} output timeseries from plant after storage", @@ -53,7 +53,9 @@ def compute(self, inputs, outputs): """ # Assign the input to the output - outputs[f"{self.config.commodity_name}_out"] = inputs[f"{self.config.commodity_name}_in"] + outputs[f"{self.config.commodity_name}_set_point"] = inputs[ + f"{self.config.commodity_name}_in" + ] def setup_partials(self): """ @@ -74,7 +76,7 @@ def setup_partials(self): # Declare partials sparsely for all elements as an identity matrix # (diagonal elements are 1.0, others are 0.0) self.declare_partials( - of=f"{self.config.commodity_name}_out", + of=f"{self.config.commodity_name}_set_point", wrt=f"{self.config.commodity_name}_in", rows=np.arange(size), cols=np.arange(size), diff --git a/h2integrate/control/control_strategies/storage/demand_openloop_controller.py b/h2integrate/control/control_strategies/storage/demand_openloop_controller.py index 64cac74f0..a07853f0d 100644 --- a/h2integrate/control/control_strategies/storage/demand_openloop_controller.py +++ b/h2integrate/control/control_strategies/storage/demand_openloop_controller.py @@ -279,7 +279,7 @@ def compute(self, inputs, outputs): # initialize outputs soc_array = outputs[f"{commodity}_soc"] unused_commodity_array = outputs[f"{commodity}_unused_commodity"] - output_array = outputs[f"{commodity}_out"] + output_array = outputs[f"{commodity}_set_point"] unmet_demand_array = outputs[f"{commodity}_unmet_demand"] # Loop through each time step @@ -338,7 +338,7 @@ def compute(self, inputs, outputs): # Record the missed load at the current time step unmet_demand_array[t] = max(0.0, (demand_t - output_array[t])) - outputs[f"{commodity}_out"] = output_array + outputs[f"{commodity}_set_point"] = output_array # Return the SOC outputs[f"{commodity}_soc"] = soc_array diff --git a/h2integrate/control/test/test_openloop_controllers.py b/h2integrate/control/test/test_openloop_controllers.py index 68f154bd1..013b14c70 100644 --- a/h2integrate/control/test/test_openloop_controllers.py +++ b/h2integrate/control/test/test_openloop_controllers.py @@ -72,7 +72,7 @@ def test_pass_through_controller(subtests): # Run the test with subtests.test("Check output"): - assert pytest.approx(prob.get_val("hydrogen_out"), rel=1e-3) == np.arange(10) + assert pytest.approx(prob.get_val("hydrogen_set_point"), rel=1e-3) == np.arange(10) # Run the test with subtests.test("Check derivatives"): @@ -80,7 +80,7 @@ def test_pass_through_controller(subtests): assert_check_totals( prob.check_totals( of=[ - "hydrogen_out", + "hydrogen_set_point", ], wrt=[ "hydrogen_in", @@ -149,7 +149,7 @@ def test_storage_demand_controller(subtests): # Run the test with subtests.test("Check output"): assert pytest.approx([0.5, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]) == prob.get_val( - "hydrogen_out" + "hydrogen_set_point" ) with subtests.test("Check curtailment"): @@ -241,7 +241,9 @@ def set_up_and_run_problem(config): # Run the test with subtests.test("Check output"): - assert pytest.approx(prob_ioe.get_val("hydrogen_out")) == prob_rte.get_val("hydrogen_out") + assert pytest.approx(prob_ioe.get_val("hydrogen_set_point")) == prob_rte.get_val( + "hydrogen_set_point" + ) with subtests.test("Check curtailment"): assert pytest.approx(prob_ioe.get_val("hydrogen_unused_commodity")) == prob_rte.get_val( @@ -324,7 +326,7 @@ def test_generic_storage_demand_controller(subtests): # # Run the test with subtests.test("Check output"): assert pytest.approx([0.5, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]) == prob.get_val( - "hydrogen_out" + "hydrogen_set_point" ) with subtests.test("Check curtailment"): @@ -396,7 +398,7 @@ def test_demand_converter_controller(subtests): # # Run the test with subtests.test("Check output"): assert pytest.approx([0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 5.0, 5.0, 5.0, 5.0]) == prob.get_val( - "hydrogen_out" + "hydrogen_set_point" ) with subtests.test("Check curtailment"): diff --git a/h2integrate/core/model_baseclasses.py b/h2integrate/core/model_baseclasses.py index 721423527..60226ce08 100644 --- a/h2integrate/core/model_baseclasses.py +++ b/h2integrate/core/model_baseclasses.py @@ -30,19 +30,21 @@ def setup(self): # n_timesteps is number of timesteps in a simulation self.n_timesteps = self.options["plant_config"]["plant"]["simulation"]["n_timesteps"] + # dt is seconds per timestep self.dt = self.options["plant_config"]["plant"]["simulation"]["dt"] + # plant_life is number of years the plant is expected to operate for self.plant_life = int(self.options["plant_config"]["plant"]["plant_life"]) - hours_per_year = 8760 + # hours simulated is the number of hours in a simulation hours_simulated = (self.dt / 3600) * self.n_timesteps + # fraction_of_year_simulated is the ratio of simulation length to length of year # and may be used to estimate annual performance from simulation performance + hours_per_year = 8760 self.fraction_of_year_simulated = hours_simulated / hours_per_year - self.set_outputs() - def set_outputs(self): # Check that the required attributes have been instantiated required = ("commodity", "commodity_rate_units", "commodity_amount_units") missing = [el for el in required if not hasattr(self, el)] diff --git a/h2integrate/postprocess/test/test_sql_timeseries_to_csv.py b/h2integrate/postprocess/test/test_sql_timeseries_to_csv.py index e9c7a4893..1dd80538e 100644 --- a/h2integrate/postprocess/test/test_sql_timeseries_to_csv.py +++ b/h2integrate/postprocess/test/test_sql_timeseries_to_csv.py @@ -35,7 +35,7 @@ def test_save_csv_all_results(subtests, run_example_02_sql_fpath): res = save_case_timeseries_as_csv(run_example_02_sql_fpath, save_to_file=True) with subtests.test("Check number of columns"): - assert len(res.columns.to_list()) == 35 + assert len(res.columns.to_list()) == 36 with subtests.test("Check number of rows"): assert len(res) == 8760 diff --git a/h2integrate/storage/simple_generic_storage.py b/h2integrate/storage/simple_generic_storage.py index b378e5a9d..f67471200 100644 --- a/h2integrate/storage/simple_generic_storage.py +++ b/h2integrate/storage/simple_generic_storage.py @@ -1,7 +1,7 @@ -import openmdao.api as om from attrs import field, define from h2integrate.core.utilities import BaseConfig, merge_shared_inputs +from h2integrate.core.model_baseclasses import PerformanceModelBaseClass @define(kw_only=True) @@ -10,26 +10,48 @@ class SimpleGenericStorageConfig(BaseConfig): commodity_units: str = field() # TODO: update to commodity_rate_units -class SimpleGenericStorage(om.ExplicitComponent): - """ - Simple generic storage model. +class SimpleGenericStorage(PerformanceModelBaseClass): """ + Simple generic storage model that acts as a pass-through component. + + Note: this storage performance model is intended to be used with the + `DemandOpenLoopStorageController` controller. - def initialize(self): - self.options.declare("tech_config", types=dict) - self.options.declare("plant_config", types=dict) - self.options.declare("driver_config", types=dict) + """ def setup(self): - n_timesteps = self.options["plant_config"]["plant"]["simulation"]["n_timesteps"] self.config = SimpleGenericStorageConfig.from_dict( merge_shared_inputs(self.options["tech_config"]["model_inputs"], "performance"), strict=False, additional_cls_name=self.__class__.__name__, ) - commodity_name = self.config.commodity_name - commodity_units = self.config.commodity_units - self.add_input(f"{commodity_name}_in", val=0.0, shape=n_timesteps, units=commodity_units) + self.commodity = self.config.commodity_name + self.commodity_rate_units = self.config.commodity_units + self.commodity_amount_units = f"({self.commodity_rate_units})*h" + super().setup() + self.add_input( + f"{self.commodity}_set_point", + val=0.0, + shape=self.n_timesteps, + units=self.commodity_rate_units, + ) def compute(self, inputs, outputs): - pass + # Pass the commodity_out as the commodity_set_point + outputs[f"{self.commodity}_out"] = inputs[f"{self.commodity}_set_point"] + + # Estimate the rated commodity production as the maximum value from the commodity_set_point + outputs[f"rated_{self.commodity}_production"] = inputs[f"{self.commodity}_set_point"].max() + + # Calculate the total and annual commodity produced + outputs[f"total_{self.commodity}_produced"] = outputs[f"{self.commodity}_out"].sum() + outputs[f"annual_{self.commodity}_produced"] = outputs[ + f"total_{self.commodity}_produced" + ] * (1 / self.fraction_of_year_simulated) + + # Calculate the maximum commodity production over the simulation + max_production = ( + inputs[f"{self.commodity}_set_point"].max() * self.n_timesteps * (self.dt / 3600) + ) + + outputs["capacity_factor"] = outputs[f"total_{self.commodity}_produced"] / max_production diff --git a/h2integrate/storage/simple_storage_auto_sizing.py b/h2integrate/storage/simple_storage_auto_sizing.py index 7d6fb32db..4b7b7287b 100644 --- a/h2integrate/storage/simple_storage_auto_sizing.py +++ b/h2integrate/storage/simple_storage_auto_sizing.py @@ -1,8 +1,10 @@ +from copy import deepcopy + import numpy as np -import openmdao.api as om from attrs import field, define from h2integrate.core.utilities import BaseConfig, merge_shared_inputs +from h2integrate.core.model_baseclasses import PerformanceModelBaseClass @define(kw_only=True) @@ -22,30 +24,49 @@ class StorageSizingModelConfig(BaseConfig): demand_profile: int | float | list = field(default=0.0) -class StorageAutoSizingModel(om.ExplicitComponent): +class StorageAutoSizingModel(PerformanceModelBaseClass): """Performance model that calculates the storage charge rate and capacity needed to either: - 1. supply the comodity at a constant rate based on the commodity production profile or + 1. supply the commodity at a constant rate based on the commodity production profile or 2. try to meet the commodity demand with the given commodity production profile. + Then simulates performance of a basic storage component using the charge rate and + capacity calculated. + + Note: this storage performance model is intended to be used with the + `PassThroughOpenLoopController` controller and is not compatible with the + `DemandOpenLoopStorageController` controller. + Inputs: - {commodity_name}_in (float): Input commodity flow timeseries (e.g., hydrogen production). + {commodity}_in (float): Input commodity flow timeseries (e.g., hydrogen production) + used to estimate the demand if `commodity_demand_profile` is zero. - Units: Defined in `commodity_units` (e.g., "kg/h"). - {commodity_name}_demand_profile (float): Demand profile of commodity. + {commodity}_set_point (float): Input commodity flow timeseries (e.g., hydrogen production) + used as the available input commodity to meet the demand. + {commodity}_demand_profile (float): Demand profile of commodity. - Units: Defined in `commodity_units` (e.g., "kg/h"). + Outputs: max_capacity (float): Maximum storage capacity of the commodity. - Units: in non-rate units, e.g., "kg" if `commodity_units` is "kg/h" max_charge_rate (float): Maximum rate at which the commodity can be charged - Units: Defined in `commodity_units` (e.g., "kg/h"). - """ + Assumed to also be the discharge rate. + {commodity}_out (np.ndarray): the commodity used to meet demand from the available + input commodity and storage component. Defined in `commodity_units`. + total_{commodity}_produced (float): sum of commodity discharged from storage over + the simulation. Defined in `commodity_units*h` + rated_{commodity}_production (float): maximum commodity that could be discharged + in a timestep. Defined in `commodity_units` + annual_{commodity}_produced (np.ndarray): total commodity discharged per year. + Defined in `commodity_units*h/year` + capacity_factor (np.ndarray): ratio of commodity discharged to the maximum + commodity that could be discharged over the simulation. + Defined as a ratio (units of `unitless`) - def initialize(self): - self.options.declare("driver_config", types=dict) - self.options.declare("plant_config", types=dict) - self.options.declare("tech_config", types=dict) + """ def setup(self): self.config = StorageSizingModelConfig.from_dict( @@ -54,25 +75,32 @@ def setup(self): additional_cls_name=self.__class__.__name__, ) - super().setup() - - n_timesteps = self.options["plant_config"]["plant"]["simulation"]["n_timesteps"] + self.commodity = self.config.commodity_name + self.commodity_rate_units = self.config.commodity_units + self.commodity_amount_units = f"({self.commodity_rate_units})*h" - commodity_name = self.config.commodity_name + super().setup() self.add_input( - f"{commodity_name}_demand_profile", + f"{self.commodity}_demand_profile", units=f"{self.config.commodity_units}", val=self.config.demand_profile, - shape=n_timesteps, - desc=f"{commodity_name} demand profile timeseries", + shape=self.n_timesteps, + desc=f"{self.commodity} demand profile timeseries", ) self.add_input( - f"{commodity_name}_in", + f"{self.commodity}_in", shape_by_conn=True, units=f"{self.config.commodity_units}", - desc=f"{commodity_name} input timeseries from production to storage", + desc=f"{self.commodity} input timeseries from production to storage", + ) + + self.add_input( + f"{self.commodity}_set_point", + shape_by_conn=True, + units=f"{self.config.commodity_units}", + desc=f"{self.commodity} input set point from controller", ) self.add_output( @@ -90,38 +118,104 @@ def setup(self): ) def compute(self, inputs, outputs): - commodity_name = self.config.commodity_name - storage_max_fill_rate = np.max(inputs[f"{commodity_name}_in"]) + # Step 1: Auto-size the storage to meet the demand + + # Auto-size the fill rate as the max of the input commodity + storage_max_fill_rate = np.max(inputs[f"{self.commodity}_in"]) - ########### get storage size ########### - if np.sum(inputs[f"{commodity_name}_demand_profile"]) > 0: - commodity_demand = inputs[f"{commodity_name}_demand_profile"] + # Set the demand profile + if np.sum(inputs[f"{self.commodity}_demand_profile"]) > 0: + commodity_demand = inputs[f"{self.commodity}_demand_profile"] else: - commodity_demand = np.mean( - inputs[f"{commodity_name}_in"] + # If the commodity_demand_profile is zero, use the average + # commodity_in as the demand + commodity_demand = np.mean(inputs[f"{self.commodity}_in"]) * np.ones( + self.n_timesteps ) # TODO: update demand based on end-use needs - commodity_production = inputs[f"{commodity_name}_in"] + # The commodity_set_point is the production set by the controller + commodity_production = inputs[f"{self.commodity}_set_point"] # TODO: SOC is just an absolute value and is not a percentage. Ideally would calculate as shortfall in future. + # Size the storage capacity to meet the demand as much as possible commodity_storage_soc = [] for j in range(len(commodity_production)): if j == 0: - commodity_storage_soc.append(commodity_production[j] - commodity_demand) + commodity_storage_soc.append(commodity_production[j] - commodity_demand[j]) else: commodity_storage_soc.append( - commodity_storage_soc[j - 1] + commodity_production[j] - commodity_demand + commodity_storage_soc[j - 1] + commodity_production[j] - commodity_demand[j] ) minimum_soc = np.min(commodity_storage_soc) - # adjust soc so it's not negative. + # Adjust soc so it's not negative. if minimum_soc < 0: commodity_storage_soc = [x + np.abs(minimum_soc) for x in commodity_storage_soc] + # Calculate the maximum hydrogen storage capacity needed to meet the demand commodity_storage_capacity_kg = np.max(commodity_storage_soc) - np.min( commodity_storage_soc ) + # Step 2: Simulate the storage performance based on the sizes calculated + + # Initialize output arrays of charge and discharge + discharge_storage = np.zeros(self.n_timesteps) + charge_storage = np.zeros(self.n_timesteps) + output_array = np.zeros(self.n_timesteps) + + # Initialize state-of-charge value as the soc at t=0 + soc = deepcopy(commodity_storage_soc[0]) + + # Simulate a basic storage component + for t, demand_t in enumerate(commodity_demand): + input_flow = commodity_production[t] + available_charge = float(commodity_storage_capacity_kg - soc) + available_discharge = float(soc) + + # If demand is greater than the input, discharge storage + if demand_t > input_flow: + # Discharge storage to meet demand. + discharge_needed = demand_t - input_flow + discharge = min(discharge_needed, available_discharge, storage_max_fill_rate) + # Update SOC + soc -= discharge + + discharge_storage[t] = discharge + output_array[t] = input_flow + discharge + + # If input is greater than the demand, charge storage + else: + # Charge storage with unused input + unused_input = input_flow - demand_t + charge = min(unused_input, available_charge, storage_max_fill_rate) + # Update SOC + soc += charge + + charge_storage[t] = charge + output_array[t] = demand_t + + # Output the storage sizes (charge rate and capacity) outputs["max_charge_rate"] = storage_max_fill_rate outputs["max_capacity"] = commodity_storage_capacity_kg + + # commodity_out is the commodity_set_point - charge_storage + discharge_storage + outputs[f"{self.commodity}_out"] = output_array + + # The rated_commodity_production is based on the discharge rate + # (which is assumed equal to the charge rate) + outputs[f"rated_{self.commodity}_production"] = storage_max_fill_rate + + # The total_commodity_produced is the sum of the commodity discharged from storage + outputs[f"total_{self.commodity}_produced"] = discharge_storage.sum() + # Adjust the total_commodity_produced to a year-long simulation + outputs[f"annual_{self.commodity}_produced"] = outputs[ + f"total_{self.commodity}_produced" + ] * (1 / self.fraction_of_year_simulated) + + # The maximum production is based on the charge/discharge rate + max_production = storage_max_fill_rate * self.n_timesteps * (self.dt / 3600) + + # Capacity factor is total discharged commodity / maximum discharged commodity possible + outputs["capacity_factor"] = outputs[f"total_{self.commodity}_produced"] / max_production