diff --git a/h2integrate/storage/storage_baseclass.py b/h2integrate/storage/storage_baseclass.py index 96997fa85..4e4d36c54 100644 --- a/h2integrate/storage/storage_baseclass.py +++ b/h2integrate/storage/storage_baseclass.py @@ -1,5 +1,6 @@ import numpy as np from attrs import field, define +from openmdao.utils import units as om_units from h2integrate.core.utilities import BaseConfig from h2integrate.core.validators import range_val @@ -46,8 +47,8 @@ class StoragePerformanceBase(PerformanceModelBaseClass): """ _time_step_bounds = ( - 3600, - 3600, + 1, + 36000, ) # (min, max) time step lengths (in seconds) compatible with this model def setup(self): @@ -185,10 +186,16 @@ def setup(self): ) self.using_feedback_control = using_feedback_control - # convert from seconds to hours - self.dt_hr = int(self.options["plant_config"]["plant"]["simulation"]["dt"]) / ( - 3600 - ) # convert from seconds to hours + # convert from seconds to hours (kept for PySAM and legacy callers) + self.dt_hr = self.dt / 3600.0 + + # dt expressed in (commodity_amount_units / commodity_rate_units), i.e. the + # timestep width in whatever time unit makes rate * dt_amount = amount. + self.dt_amount = om_units.convert_units( + self.dt, + "s", + f"({self.commodity_amount_units})/({self.commodity_rate_units})", + ) def compute(self, inputs, outputs, discrete_inputs=[], discrete_outputs=[]): """Run the storage model. @@ -296,7 +303,8 @@ def run_storage( # Performance model outputs outputs[f"rated_{self.commodity}_production"] = discharge_rate - outputs[f"total_{self.commodity}_produced"] = np.sum(storage_commodity_out) + # rate * dt_amount = commodity_amount_units (works for any commodity_rate_units) + outputs[f"total_{self.commodity}_produced"] = np.sum(storage_commodity_out) * self.dt_amount outputs[f"annual_{self.commodity}_produced"] = outputs[ f"total_{self.commodity}_produced" ] * (1 / self.fraction_of_year_simulated) @@ -306,12 +314,14 @@ def run_storage( outputs["standard_capacity_factor"] = 0.0 else: outputs["capacity_factor"] = outputs[f"total_{self.commodity}_produced"] / ( - outputs[f"rated_{self.commodity}_production"] * self.n_timesteps + outputs[f"rated_{self.commodity}_production"] * self.n_timesteps * self.dt_amount ) # standard_capacity_factor is the ratio of commodity discharged to the discharge rate - total_commodity_discharged = outputs[f"storage_{self.commodity}_discharge"].sum() + total_commodity_discharged = ( + outputs[f"storage_{self.commodity}_discharge"].sum() * self.dt_amount + ) outputs["standard_capacity_factor"] = total_commodity_discharged / ( - outputs[f"rated_{self.commodity}_production"] * self.n_timesteps + outputs[f"rated_{self.commodity}_production"] * self.n_timesteps * self.dt_amount ) return outputs @@ -402,7 +412,7 @@ def simulate( # --- Charging --- # headroom: how much more commodity the storage can accept, # expressed as a rate (commodity_rate_units). - headroom = (soc_max - soc) * storage_capacity / self.dt_hr + headroom = (soc_max - soc) * storage_capacity / self.dt_amount # charge available based on the available input commodity charge_available = commodity_available[sim_start_index + t] @@ -418,7 +428,7 @@ def simulate( ) # Update SOC (actual_charge is in post-efficiency units) - soc += actual_charge / storage_capacity + soc += actual_charge * self.dt_amount / storage_capacity # Update the amount of commodity used to charge from the input stream # If charge_eff<1, more commodity is pulled from the input stream than @@ -428,7 +438,7 @@ def simulate( # --- Discharging --- # headroom: how much commodity can still be drawn before # hitting the minimum SOC, expressed as a rate. - headroom = (soc - soc_min) * storage_capacity / self.dt_hr + headroom = (soc - soc_min) * storage_capacity / self.dt_amount # Clip to the most restrictive limit without applied efficiency. # Efficiency losses occur as energy leaves storage. @@ -437,7 +447,7 @@ def simulate( ) # Update SOC (actual_discharge is before efficiency losses are applied.) - soc -= actual_discharge / storage_capacity + soc -= actual_discharge * self.dt_amount / storage_capacity # If discharge_eff<1, then less commodity is output from the storage # than the commodity discharged from storage diff --git a/h2integrate/storage/storage_performance_model.py b/h2integrate/storage/storage_performance_model.py index 2486f1d66..d597446b7 100644 --- a/h2integrate/storage/storage_performance_model.py +++ b/h2integrate/storage/storage_performance_model.py @@ -120,7 +120,7 @@ class StoragePerformanceModel(StoragePerformanceBase): """OpenMDAO component for a storage component.""" _time_step_bounds = ( - 3600, + 1, 3600, ) # (min, max) time step lengths (in seconds) compatible with this model diff --git a/h2integrate/storage/test/test_storage_performance_model.py b/h2integrate/storage/test/test_storage_performance_model.py index bc70c7de4..2ad903400 100644 --- a/h2integrate/storage/test/test_storage_performance_model.py +++ b/h2integrate/storage/test/test_storage_performance_model.py @@ -1210,3 +1210,311 @@ def test_round_trip_efficiency_preserved_in_config(subtests): assert config_dict["round_trip_efficiency"] == round_trip_eff assert config_dict["charge_efficiency"] == pytest.approx(np.sqrt(round_trip_eff)) assert config_dict["discharge_efficiency"] == pytest.approx(np.sqrt(round_trip_eff)) + + +@pytest.fixture +def plant_config_non_hourly(n_timesteps): + plant = { + "plant": { + "plant_life": 30, + "simulation": { + "dt": 1800, + "n_timesteps": n_timesteps, + }, + }, + } + return plant + + +@pytest.mark.regression +@pytest.mark.parametrize("n_timesteps", [4]) +def test_storage_half_hourly_known_outputs(subtests, plant_config_non_hourly): + """Verify SOC, charge/discharge profiles, and scalar outputs against calculated + values for a simple scenario at dt=1800s (30-min dt). + + Scenario (4 timesteps * 30 min = 2 hours total): + t0, t1: charge at 10 kg/h - stores 5 kg each step + t2, t3: discharge at 10 kg/h — removes 5 kg each step + + With capacity=40 kg, init_soc=0.1, eff=1.0, min_soc=0.1, max_soc=1.0: + SOC[0] = 0.1 + 5/40 = 0.225 + SOC[1] = 0.225 + 5/40 = 0.35 + SOC[2] = 0.35 - 5/40 = 0.225 + SOC[3] = 0.225 - 5/40 = 0.1 + total_hydrogen_produced = (-10 - 10 + 10 + 10) * 0.5 hr = 0 kg + standard_capacity_factor = (10+10)*0.5 / (10 * 4 * 0.5) = 10/20 = 0.5 -> 50 % + """ + + model_inputs = { + "shared_parameters": { + "commodity": "hydrogen", + "commodity_rate_units": "kg/h", + }, + "performance_parameters": { + "max_capacity": 40.0, + "max_charge_rate": 10.0, + "min_soc_fraction": 0.1, + "max_soc_fraction": 1.0, + "init_soc_fraction": 0.1, + "commodity_amount_units": "kg", + "charge_equals_discharge": True, + "charge_efficiency": 1.0, + "discharge_efficiency": 1.0, + "demand_profile": 0.0, + }, + } + + commodity_in = np.array([10.0, 10.0, 0.0, 0.0]) + set_point = np.array([-10.0, -10.0, 10.0, 10.0]) + + prob = om.Problem() + prob.model.add_subsystem( + "IVC1", + om.IndepVarComp("hydrogen_in", val=commodity_in, units="kg/h"), + promotes=["*"], + ) + prob.model.add_subsystem( + "IVC2", + om.IndepVarComp("hydrogen_set_point", val=set_point, units="kg/h"), + promotes=["*"], + ) + prob.model.add_subsystem( + "storage", + StoragePerformanceModel( + plant_config=plant_config_non_hourly, + tech_config={"model_inputs": model_inputs}, + ), + promotes=["*"], + ) + prob.setup() + prob.run_model() + + with subtests.test("SOC profile matches hand-calculated values"): + expected_soc_pct = np.array([22.5, 35.0, 22.5, 10.0]) + np.testing.assert_allclose( + prob.get_val("storage.SOC", units="percent"), expected_soc_pct, rtol=1e-9 + ) + + with subtests.test("Charge profile"): + np.testing.assert_allclose( + prob.get_val("storage.storage_hydrogen_charge", units="kg/h"), + np.array([-10.0, -10.0, 0.0, 0.0]), + rtol=1e-9, + ) + + with subtests.test("Discharge profile"): + np.testing.assert_allclose( + prob.get_val("storage.storage_hydrogen_discharge", units="kg/h"), + np.array([0.0, 0.0, 10.0, 10.0]), + rtol=1e-9, + ) + + with subtests.test("total_hydrogen_produced = 0 kg (charge equals discharge)"): + assert ( + pytest.approx(prob.get_val("storage.total_hydrogen_produced", units="kg")[0], abs=1e-9) + == 0.0 + ) + + with subtests.test("standard_capacity_factor = 50 %"): + assert ( + pytest.approx( + prob.get_val("storage.standard_capacity_factor", units="percent")[0], rel=1e-9 + ) + == 50.0 + ) + + +@pytest.mark.regression +@pytest.mark.parametrize("n_timesteps", [4]) +def test_storage_half_hourly_known_outputs_kg_s(subtests, plant_config_non_hourly): + """Verify SOC, charge/discharge profiles against hand-calculated values when + commodity_rate_units='kg/s' and dt=1800s. + + This test verifies that dt_amount is correctly computed via OpenMDAO + that SOC increments are correct for a non-hourly rate unit. + + Scenario + capacity=3600 kg, charge_rate=1 kg/s, init_soc=0.1, min/max_soc=0.1/1.0, eff=1.0 + commodity_in = [1, 1, 0, 0] kg/s + set_point = [-1, -1, 1, 1] kg/s (negative=charge, positive=discharge) + + With 1800s (0.5 hr) timesteps, the expected SOC profile is: + SOC[0] = 0.1 + (1 kg/s * 1800s) / 3600 kg = 0.1 + 0.5 = 0.6 -> 60% + SOC[1] = 0.6 + (1 kg/s * 1800s) / 3600 kg = 0.6 + 0 = 1.0 -> 100% + SOC[2] = 1.0 - (1 kg/s * 1800s) / 3600 kg = 1.0 - 0.5 = 0.5 -> 50% + SOC[3] = 0.5 - (1 kg/s * 1800s) / 3600 kg = 0.5 - 0.5 = 0.1 -> 10% + + total_produced = (-1-0.8+1+0.8)*1800 = 0 kg + standard_capacity_factor = (1+0.8)*1800 / (1*4*1800) = 3240/7200 = 45% + """ + + model_inputs = { + "shared_parameters": { + "commodity": "hydrogen", + "commodity_rate_units": "kg/s", + }, + "performance_parameters": { + "max_capacity": 3600.0, + "max_charge_rate": 1.0, + "min_soc_fraction": 0.1, + "max_soc_fraction": 1.0, + "init_soc_fraction": 0.1, + "commodity_amount_units": "kg", + "charge_equals_discharge": True, + "charge_efficiency": 1.0, + "discharge_efficiency": 1.0, + "demand_profile": 0.0, + }, + } + + commodity_in = np.array([1.0, 1.0, 0.0, 0.0]) + set_point = np.array([-1.0, -1.0, 1.0, 1.0]) + + prob = om.Problem() + prob.model.add_subsystem( + "IVC1", + om.IndepVarComp("hydrogen_in", val=commodity_in, units="kg/s"), + promotes=["*"], + ) + prob.model.add_subsystem( + "IVC2", + om.IndepVarComp("hydrogen_set_point", val=set_point, units="kg/s"), + promotes=["*"], + ) + prob.model.add_subsystem( + "storage", + StoragePerformanceModel( + plant_config=plant_config_non_hourly, + tech_config={"model_inputs": model_inputs}, + ), + promotes=["*"], + ) + prob.setup() + prob.run_model() + + with subtests.test("SOC profile matches hand-calculated values"): + expected_soc_pct = np.array([60.0, 100.0, 50.0, 10.0]) + np.testing.assert_allclose( + prob.get_val("storage.SOC", units="percent"), expected_soc_pct, rtol=1e-9 + ) + + with subtests.test("Charge profile"): + np.testing.assert_allclose( + prob.get_val("storage.storage_hydrogen_charge", units="kg/s"), + np.array([-1.0, -0.8, 0.0, 0.0]), + rtol=1e-9, + ) + + with subtests.test("Discharge profile"): + np.testing.assert_allclose( + prob.get_val("storage.storage_hydrogen_discharge", units="kg/s"), + np.array([0.0, 0.0, 1.0, 0.8]), + rtol=1e-9, + ) + + with subtests.test("total_hydrogen_produced = 0 kg (charge equals discharge)"): + assert ( + pytest.approx(prob.get_val("storage.total_hydrogen_produced", units="kg")[0], abs=1e-9) + == 0.0 + ) + + with subtests.test("standard_capacity_factor = 45 %"): + assert ( + pytest.approx( + prob.get_val("storage.standard_capacity_factor", units="percent")[0], rel=1e-9 + ) + == 45.0 + ) + + +@pytest.mark.regression +def test_storage_half_hourly_kw_kwh_2hr(subtests): + model_inputs = { + "shared_parameters": { + "commodity": "electricity", + "commodity_rate_units": "kW", + }, + "performance_parameters": { + "max_capacity": 10.0, + "max_charge_rate": 10.0, + "min_soc_fraction": 0.1, + "max_soc_fraction": 1.0, + "init_soc_fraction": 0.1, + "commodity_amount_units": "kW*h", + "charge_equals_discharge": True, + "charge_efficiency": 1.0, + "discharge_efficiency": 1.0, + "demand_profile": 0.0, + }, + } + plant_config_non_hourly = { + "plant": { + "plant_life": 30, + "simulation": { + "dt": 7200, + "n_timesteps": 4, + }, + }, + } + + commodity_in = np.array([10.0, 10.0, 10.0, 10.0]) + set_point = np.array([-10.0, -10.0, 10.0, 10.0]) + + prob = om.Problem() + prob.model.add_subsystem( + "IVC1", + om.IndepVarComp("electricity_in", val=commodity_in, units="kW"), + promotes=["*"], + ) + prob.model.add_subsystem( + "IVC2", + om.IndepVarComp("electricity_set_point", val=set_point, units="kW"), + promotes=["*"], + ) + prob.model.add_subsystem( + "storage", + StoragePerformanceModel( + plant_config=plant_config_non_hourly, + tech_config={"model_inputs": model_inputs}, + ), + promotes=["*"], + ) + prob.setup() + prob.run_model() + + with subtests.test("SOC profile"): + np.testing.assert_allclose( + prob.get_val("storage.SOC", units="percent"), + np.array([100.0, 100.0, 10.0, 10.0]), + rtol=1e-9, + ) + + with subtests.test("Charge profile"): + np.testing.assert_allclose( + prob.get_val("storage.storage_electricity_charge", units="kW"), + np.array([-4.5, -0.0, 0.0, 0.0]), + rtol=1e-9, + ) + + with subtests.test("Discharge profile"): + np.testing.assert_allclose( + prob.get_val("storage.storage_electricity_discharge", units="kW"), + np.array([0.0, 0.0, 4.5, 0.0]), + rtol=1e-9, + ) + + with subtests.test("total_electricity_produced = 0 kWh"): + assert ( + pytest.approx( + prob.get_val("storage.total_electricity_produced", units="kW*h")[0], abs=1e-9 + ) + == 0.0 + ) + + with subtests.test("standard_capacity_factor"): + assert ( + pytest.approx( + prob.get_val("storage.standard_capacity_factor", units="percent")[0], rel=1e-9 + ) + == 11.25 + )