From 553742ee95b9624ee5eed74287f92af436144cbf Mon Sep 17 00:00:00 2001 From: Jared Thomas Date: Thu, 28 May 2026 12:55:19 -0600 Subject: [PATCH 01/18] create battery model files for including degradation --- h2integrate/core/supported_models.py | 1 + h2integrate/storage/battery/__init__.py | 1 + .../storage/battery/battery_performance.py | 164 ++++++++++++++++++ .../battery/test/test_battery_performance.py | 1 + 4 files changed, 167 insertions(+) create mode 100644 h2integrate/storage/battery/battery_performance.py create mode 100644 h2integrate/storage/battery/test/test_battery_performance.py diff --git a/h2integrate/core/supported_models.py b/h2integrate/core/supported_models.py index a22ae427e..429d60bae 100644 --- a/h2integrate/core/supported_models.py +++ b/h2integrate/core/supported_models.py @@ -141,6 +141,7 @@ def copy(self): "GenericSummerPerformanceModel": "transporters:GenericSummerPerformanceModel", # Storage "PySAMBatteryPerformanceModel": "storage.battery:PySAMBatteryPerformanceModel", + "BatteryPerformanceModel": "storage.battery:BatteryPerformanceModel", "StoragePerformanceModel": "storage:StoragePerformanceModel", "StorageAutoSizingModel": "storage:StorageAutoSizingModel", "LinedRockCavernStorageCostModel": "storage.hydrogen:LinedRockCavernStorageCostModel", diff --git a/h2integrate/storage/battery/__init__.py b/h2integrate/storage/battery/__init__.py index 1a32933c3..432e095a4 100644 --- a/h2integrate/storage/battery/__init__.py +++ b/h2integrate/storage/battery/__init__.py @@ -1,2 +1,3 @@ from h2integrate.storage.battery.pysam_battery import PySAMBatteryPerformanceModel from h2integrate.storage.battery.atb_battery_cost import ATBBatteryCostModel +from h2integrate.storage.battery.battery_performance import BatteryPerformanceModel diff --git a/h2integrate/storage/battery/battery_performance.py b/h2integrate/storage/battery/battery_performance.py new file mode 100644 index 000000000..d2b676a59 --- /dev/null +++ b/h2integrate/storage/battery/battery_performance.py @@ -0,0 +1,164 @@ +import numpy as np +from attrs import field, define + +from h2integrate.core.utilities import merge_shared_inputs +from h2integrate.core.validators import gt_zero, range_val, range_val_or_none +from h2integrate.storage.storage_baseclass import ( + StoragePerformanceBase, + StoragePerformanceBaseConfig, +) + + +@define(kw_only=True) +class BatteryPerformanceModelConfig(StoragePerformanceBaseConfig): + """Configuration class for storage performance models. + + This class defines configuration parameters for simulating storage + performance with the Pyomo controllers. It includes + specifications such as capacity, charge rate, state-of-charge limits, + and charge/discharge efficiencies. + + Attributes: + commodity (str): name of commodity + commodity_rate_units (str): Units of the commodity (e.g., "kg/h"). + demand_profile (int | float | list): Demand values for each timestep, in + the same units as `commodity_rate_units`. May be a scalar for constant + demand or a list/array for time-varying demand. + max_capacity (float): Maximum storage energy capacity in commodity_amount_units. + Must be greater than zero. + max_charge_rate (float): Rated commodity capacity of the storage in commodity_rate_units. + Must be greater than zero. + min_soc_fraction (float): Minimum allowable state of charge as a fraction (0 to 1). + max_soc_fraction (float): Maximum allowable state of charge as a fraction (0 to 1). + init_soc_fraction (float): Initial state of charge as a fraction (0 to 1). + commodity_amount_units (str | None, optional): Units of the commodity as an amount + (i.e., kW*h or kg). If not provided, defaults to commodity_rate_units*h. + max_discharge_rate (float | None, optional): Maximum rate at which the commodity can be + discharged (in units per time step, e.g., "kg/time step"). This rate does not include + the discharge_efficiency. Only required if `charge_equals_discharge` is False. + charge_equals_discharge (bool, optional): If True, set the max_discharge_rate equal to the + max_charge_rate. If False, specify the max_discharge_rate as a value different than + the max_charge_rate. Defaults to True. + charge_efficiency (float | None, optional): Efficiency of charging the storage, represented + as a decimal between 0 and 1 (e.g., 0.9 for 90% efficiency). Optional if + `round_trip_efficiency` is provided. + discharge_efficiency (float | None, optional): Efficiency of discharging the storage, + represented as a decimal between 0 and 1 (e.g., 0.9 for 90% efficiency). Optional if + `round_trip_efficiency` is provided. + round_trip_efficiency (float | None, optional): Combined efficiency of charging and + discharging the storage, represented as a decimal between 0 and 1 (e.g., 0.81 for + 81% efficiency). Optional if `charge_efficiency` and `discharge_efficiency` are + provided. + + """ + + commodity: str = field() + commodity_rate_units: str = field() + + max_capacity: float = field(validator=gt_zero) + max_charge_rate: float = field(validator=gt_zero) + + init_soc_fraction: float = field(validator=range_val(0, 1)) + + commodity_amount_units: str = field(default=None) + max_discharge_rate: float | None = field(default=None) + charge_equals_discharge: bool = field(default=True) + + charge_efficiency: float | None = field(default=None, validator=range_val_or_none(0, 1)) + discharge_efficiency: float | None = field(default=None, validator=range_val_or_none(0, 1)) + round_trip_efficiency: float | None = field(default=None, validator=range_val_or_none(0, 1)) + + # TODO degredation: add additional parameters for degradation here + + def __attrs_post_init__(self): + """ + Post-initialization logic to validate and calculate efficiencies. + + Ensures that either `charge_efficiency` and `discharge_efficiency` are provided, + or `round_trip_efficiency` is provided. If `round_trip_efficiency` is provided, + it calculates `charge_efficiency` and `discharge_efficiency` as the square root + of `round_trip_efficiency`. + """ + if (self.round_trip_efficiency is not None) and ( + self.charge_efficiency is None and self.discharge_efficiency is None + ): + # Calculate charge and discharge efficiencies from round-trip efficiency + self.charge_efficiency = np.sqrt(self.round_trip_efficiency) + self.discharge_efficiency = np.sqrt(self.round_trip_efficiency) + + if self.charge_efficiency is None or self.discharge_efficiency is None: + raise ValueError( + "Exactly one of the following sets of parameters must be set: (a) " + "`round_trip_efficiency`, or (b) both `charge_efficiency` " + "and `discharge_efficiency`." + ) + + if self.charge_equals_discharge: + if ( + self.max_discharge_rate is not None + and self.max_discharge_rate != self.max_charge_rate + ): + msg = ( + "Max discharge rate does not equal max charge rate but charge_equals_discharge " + f"is True. Discharge rate is {self.max_discharge_rate} and charge rate " + f"is {self.max_charge_rate}." + ) + raise ValueError(msg) + + self.max_discharge_rate = self.max_charge_rate + + if not self.charge_equals_discharge and self.max_discharge_rate is None: + msg = ( + "max_discharge_rate is required when charge_equals_discharge is False. " + "Please input the discharge rate using the key `max_discharge_rate`." + ) + raise ValueError(msg) + + if self.commodity_amount_units is None: + self.commodity_amount_units = f"({self.commodity_rate_units})*h" + + +class BatteryPerformanceModel(StoragePerformanceBase): + """OpenMDAO component for a storage component.""" + + _time_step_bounds = ( + 1, + 3600, + ) # (min, max) time step lengths (in seconds) compatible with this model + + def setup(self): + self.config = BatteryPerformanceModelConfig.from_dict( + merge_shared_inputs(self.options["tech_config"]["model_inputs"], "performance"), + strict=False, + additional_cls_name=self.__class__.__name__, + ) + + self.commodity = self.config.commodity + self.commodity_rate_units = self.config.commodity_rate_units + self.commodity_amount_units = self.config.commodity_amount_units + + # TODO degredation: adjustments for degradation + + super().setup() + + def compute(self, inputs, outputs, discrete_inputs=[], discrete_outputs=[]): + """Run the storage performance model.""" + self.current_soc = self.config.init_soc_fraction + + charge_rate = inputs["max_charge_rate"][0] + if "max_discharge_rate" in inputs: + discharge_rate = inputs["max_discharge_rate"][0] + else: + discharge_rate = inputs["max_charge_rate"][0] + storage_capacity = inputs["storage_capacity"][0] + + # TODO degredation: adjust compute method for degradation as needed + outputs = self.run_storage( + charge_rate, discharge_rate, storage_capacity, inputs, outputs, discrete_inputs + ) + + # TODO degredation: add degradation method here + def degradation( + self, + ): + return diff --git a/h2integrate/storage/battery/test/test_battery_performance.py b/h2integrate/storage/battery/test/test_battery_performance.py new file mode 100644 index 000000000..83808572c --- /dev/null +++ b/h2integrate/storage/battery/test/test_battery_performance.py @@ -0,0 +1 @@ +# TODO degredation: tests for battery performance From bdcfcc29889c04e0ef6b3a03f53b8243953369c0 Mon Sep 17 00:00:00 2001 From: Jared Thomas Date: Thu, 28 May 2026 14:52:30 -0600 Subject: [PATCH 02/18] spelling correction --- h2integrate/storage/battery/battery_performance.py | 8 ++++---- .../storage/battery/test/test_battery_performance.py | 2 +- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/h2integrate/storage/battery/battery_performance.py b/h2integrate/storage/battery/battery_performance.py index d2b676a59..40cd5a353 100644 --- a/h2integrate/storage/battery/battery_performance.py +++ b/h2integrate/storage/battery/battery_performance.py @@ -68,7 +68,7 @@ class BatteryPerformanceModelConfig(StoragePerformanceBaseConfig): discharge_efficiency: float | None = field(default=None, validator=range_val_or_none(0, 1)) round_trip_efficiency: float | None = field(default=None, validator=range_val_or_none(0, 1)) - # TODO degredation: add additional parameters for degradation here + # TODO degradation: add additional parameters for degradation here def __attrs_post_init__(self): """ @@ -137,7 +137,7 @@ def setup(self): self.commodity_rate_units = self.config.commodity_rate_units self.commodity_amount_units = self.config.commodity_amount_units - # TODO degredation: adjustments for degradation + # TODO degradation: adjustments for degradation super().setup() @@ -152,12 +152,12 @@ def compute(self, inputs, outputs, discrete_inputs=[], discrete_outputs=[]): discharge_rate = inputs["max_charge_rate"][0] storage_capacity = inputs["storage_capacity"][0] - # TODO degredation: adjust compute method for degradation as needed + # TODO degradation: adjust compute method for degradation as needed outputs = self.run_storage( charge_rate, discharge_rate, storage_capacity, inputs, outputs, discrete_inputs ) - # TODO degredation: add degradation method here + # TODO degradation: add degradation method here def degradation( self, ): diff --git a/h2integrate/storage/battery/test/test_battery_performance.py b/h2integrate/storage/battery/test/test_battery_performance.py index 83808572c..7f0641c82 100644 --- a/h2integrate/storage/battery/test/test_battery_performance.py +++ b/h2integrate/storage/battery/test/test_battery_performance.py @@ -1 +1 @@ -# TODO degredation: tests for battery performance +# TODO degradation: tests for battery performance From 3c48f03d50d1fb7739a4933189a28c94e0876a73 Mon Sep 17 00:00:00 2001 From: Jared Thomas Date: Thu, 28 May 2026 17:38:45 -0600 Subject: [PATCH 03/18] make solar resource data available to the battery performance model --- examples/24_solar_battery_grid/plant_config.yaml | 4 +++- examples/24_solar_battery_grid/tech_config.yaml | 2 +- h2integrate/storage/battery/battery_performance.py | 8 ++++++++ 3 files changed, 12 insertions(+), 2 deletions(-) diff --git a/examples/24_solar_battery_grid/plant_config.yaml b/examples/24_solar_battery_grid/plant_config.yaml index 3fed5d46d..114bccce4 100644 --- a/examples/24_solar_battery_grid/plant_config.yaml +++ b/examples/24_solar_battery_grid/plant_config.yaml @@ -34,8 +34,10 @@ technology_interconnections: - [solar, combiner, electricity, cable] - [grid_buy, combiner, electricity, cable] resource_to_tech_connections: - # connect the wind resource to the wind technology + # connect the solar resource to the pv solar technology - [site.solar_resource, solar, solar_resource_data] + # connect the solar resource to the battery technology + - [site.solar_resource, battery, solar_resource_data] plant: plant_life: 30 simulation: diff --git a/examples/24_solar_battery_grid/tech_config.yaml b/examples/24_solar_battery_grid/tech_config.yaml index 969cfbf63..d16665938 100644 --- a/examples/24_solar_battery_grid/tech_config.yaml +++ b/examples/24_solar_battery_grid/tech_config.yaml @@ -31,7 +31,7 @@ technologies: cost_year: 2024 battery: performance_model: - model: StoragePerformanceModel + model: BatteryPerformanceModel cost_model: model: ATBBatteryCostModel control_strategy: diff --git a/h2integrate/storage/battery/battery_performance.py b/h2integrate/storage/battery/battery_performance.py index 40cd5a353..8344241b5 100644 --- a/h2integrate/storage/battery/battery_performance.py +++ b/h2integrate/storage/battery/battery_performance.py @@ -137,6 +137,12 @@ def setup(self): self.commodity_rate_units = self.config.commodity_rate_units self.commodity_amount_units = self.config.commodity_amount_units + self.add_discrete_input( + "solar_resource_data", + val={}, + desc="Solar resource data dictionary", + ) + # TODO degradation: adjustments for degradation super().setup() @@ -153,6 +159,7 @@ def compute(self, inputs, outputs, discrete_inputs=[], discrete_outputs=[]): storage_capacity = inputs["storage_capacity"][0] # TODO degradation: adjust compute method for degradation as needed + self.degradation(discrete_inputs["solar_resource_data"]) outputs = self.run_storage( charge_rate, discharge_rate, storage_capacity, inputs, outputs, discrete_inputs ) @@ -160,5 +167,6 @@ def compute(self, inputs, outputs, discrete_inputs=[], discrete_outputs=[]): # TODO degradation: add degradation method here def degradation( self, + solar_resource_data, ): return From f62b25f6df635c954b0aec52687190ea73fadece Mon Sep 17 00:00:00 2001 From: Jared Thomas Date: Thu, 28 May 2026 17:46:08 -0600 Subject: [PATCH 04/18] add doc string comment --- h2integrate/storage/battery/battery_performance.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/h2integrate/storage/battery/battery_performance.py b/h2integrate/storage/battery/battery_performance.py index 8344241b5..d5ffd168d 100644 --- a/h2integrate/storage/battery/battery_performance.py +++ b/h2integrate/storage/battery/battery_performance.py @@ -169,4 +169,10 @@ def degradation( self, solar_resource_data, ): + """_summary_ + + Args: + solar_resource_data (_type_): a dictionary of hourly (or by timestep) solar resource + information including irradiance and temperature + """ return From 37b813a44b326b0066826141ca2773e156a8efe7 Mon Sep 17 00:00:00 2001 From: Jared Thomas Date: Fri, 29 May 2026 09:15:38 -0600 Subject: [PATCH 05/18] use plm example 33 for degradation instead of example 24 --- examples/24_solar_battery_grid/plant_config.yaml | 4 +--- examples/24_solar_battery_grid/tech_config.yaml | 2 +- examples/33_peak_load_management/plant_config.yaml | 10 ++++++++++ examples/33_peak_load_management/tech_config.yaml | 2 +- 4 files changed, 13 insertions(+), 5 deletions(-) diff --git a/examples/24_solar_battery_grid/plant_config.yaml b/examples/24_solar_battery_grid/plant_config.yaml index 114bccce4..3fed5d46d 100644 --- a/examples/24_solar_battery_grid/plant_config.yaml +++ b/examples/24_solar_battery_grid/plant_config.yaml @@ -34,10 +34,8 @@ technology_interconnections: - [solar, combiner, electricity, cable] - [grid_buy, combiner, electricity, cable] resource_to_tech_connections: - # connect the solar resource to the pv solar technology + # connect the wind resource to the wind technology - [site.solar_resource, solar, solar_resource_data] - # connect the solar resource to the battery technology - - [site.solar_resource, battery, solar_resource_data] plant: plant_life: 30 simulation: diff --git a/examples/24_solar_battery_grid/tech_config.yaml b/examples/24_solar_battery_grid/tech_config.yaml index d16665938..969cfbf63 100644 --- a/examples/24_solar_battery_grid/tech_config.yaml +++ b/examples/24_solar_battery_grid/tech_config.yaml @@ -31,7 +31,7 @@ technologies: cost_year: 2024 battery: performance_model: - model: BatteryPerformanceModel + model: StoragePerformanceModel cost_model: model: ATBBatteryCostModel control_strategy: diff --git a/examples/33_peak_load_management/plant_config.yaml b/examples/33_peak_load_management/plant_config.yaml index 93be1cecf..26ef6faad 100644 --- a/examples/33_peak_load_management/plant_config.yaml +++ b/examples/33_peak_load_management/plant_config.yaml @@ -2,14 +2,24 @@ name: plant_config description: Demonstrates multivariable streams with a gas combiner plant: plant_life: 30 + latitude: 30.6617 + longitude: -101.7096 simulation: n_timesteps: 8760 dt: 3600 timezone: -6 start_time: 2025/07/01 00:00:00 + resources: + solar_resource: + resource_model: GOESAggregatedSolarAPI + resource_parameters: + resource_year: 2025 technology_interconnections: - [grid_buy, battery, electricity, cable] # include battery charge/discharge in the load - [battery, electrical_load_demand, [electricity_out, electricity_in]] # buy power from the grid to fulfill demand including to accommodate battery operation - [electrical_load_demand, grid_buy, [unmet_electricity_demand_out, electricity_set_point]] +resource_to_tech_connections: + # connect the solar resource to the battery technology + - [site.solar_resource, battery, solar_resource_data] diff --git a/examples/33_peak_load_management/tech_config.yaml b/examples/33_peak_load_management/tech_config.yaml index e1a9ff88f..9f247d33b 100644 --- a/examples/33_peak_load_management/tech_config.yaml +++ b/examples/33_peak_load_management/tech_config.yaml @@ -3,7 +3,7 @@ description: This plant charges a battery from the grid to reduce peak demand technologies: battery: performance_model: - model: StoragePerformanceModel + model: BatteryPerformanceModel cost_model: model: ATBBatteryCostModel control_strategy: From 185c6d67dfb061bd0d03a43d78a089d6d72852c7 Mon Sep 17 00:00:00 2001 From: Jared Thomas Date: Thu, 4 Jun 2026 10:56:00 -0600 Subject: [PATCH 06/18] add performance param example --- examples/33_peak_load_management/tech_config.yaml | 3 +++ h2integrate/storage/battery/battery_performance.py | 1 + 2 files changed, 4 insertions(+) diff --git a/examples/33_peak_load_management/tech_config.yaml b/examples/33_peak_load_management/tech_config.yaml index 9f247d33b..75c0d6d04 100644 --- a/examples/33_peak_load_management/tech_config.yaml +++ b/examples/33_peak_load_management/tech_config.yaml @@ -22,6 +22,9 @@ technologies: charge_efficiency: 1.0 # percent as decimal discharge_efficiency: 1.0 # percent as decimal demand_profile: !include demand_profiles/demand_profile.yaml + performance_parameters: + cop: 0.3 + # TODO: add remaining performance parameters here control_parameters: demand_profile_upstream: !include demand_profiles/demand_profile_upstream.yaml # demand used to define when the supervisor commands battery dispatch. This may represent an upstream load. dispatch_priority_demand_profile: demand_profile_upstream # demand profile the controller prioritizes when deciding dispatch diff --git a/h2integrate/storage/battery/battery_performance.py b/h2integrate/storage/battery/battery_performance.py index d5ffd168d..971fb380e 100644 --- a/h2integrate/storage/battery/battery_performance.py +++ b/h2integrate/storage/battery/battery_performance.py @@ -69,6 +69,7 @@ class BatteryPerformanceModelConfig(StoragePerformanceBaseConfig): round_trip_efficiency: float | None = field(default=None, validator=range_val_or_none(0, 1)) # TODO degradation: add additional parameters for degradation here + cop: float = field(validator=gt_zero) def __attrs_post_init__(self): """ From fd86452cd4a018684502f5836a48e19fd0957bfd Mon Sep 17 00:00:00 2001 From: Jared Thomas Date: Fri, 12 Jun 2026 11:02:56 -0600 Subject: [PATCH 07/18] add battery auxiliary power output --- h2integrate/storage/battery/battery_performance.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/h2integrate/storage/battery/battery_performance.py b/h2integrate/storage/battery/battery_performance.py index 971fb380e..ed526b134 100644 --- a/h2integrate/storage/battery/battery_performance.py +++ b/h2integrate/storage/battery/battery_performance.py @@ -144,6 +144,12 @@ def setup(self): desc="Solar resource data dictionary", ) + self.add_output( + f"{self.commodity}_auxiliary_demand", + shape=self.n_timesteps, + desc="Electricity demand for running battery auxiliary systems", + ) + # TODO degradation: adjustments for degradation super().setup() From 33ddd07e1e371d37adf1ed926bec08ec9f1b4c5a Mon Sep 17 00:00:00 2001 From: Jared Thomas Date: Fri, 12 Jun 2026 11:38:47 -0600 Subject: [PATCH 08/18] add commented code for assigning auxiliary power --- h2integrate/storage/battery/battery_performance.py | 1 + 1 file changed, 1 insertion(+) diff --git a/h2integrate/storage/battery/battery_performance.py b/h2integrate/storage/battery/battery_performance.py index ed526b134..04db3e9f3 100644 --- a/h2integrate/storage/battery/battery_performance.py +++ b/h2integrate/storage/battery/battery_performance.py @@ -167,6 +167,7 @@ def compute(self, inputs, outputs, discrete_inputs=[], discrete_outputs=[]): # TODO degradation: adjust compute method for degradation as needed self.degradation(discrete_inputs["solar_resource_data"]) + # outputs[f"{self.commodity}_auxiliary_demand"] = outputs = self.run_storage( charge_rate, discharge_rate, storage_capacity, inputs, outputs, discrete_inputs ) From 879221f446a541b433c6dd4124121dac47a0f79b Mon Sep 17 00:00:00 2001 From: Jared Thomas Date: Wed, 15 Jul 2026 08:18:49 -0600 Subject: [PATCH 09/18] wip: integrate simses --- .../storage/battery/battery_performance.py | 125 +++++++++++++++--- 1 file changed, 108 insertions(+), 17 deletions(-) diff --git a/h2integrate/storage/battery/battery_performance.py b/h2integrate/storage/battery/battery_performance.py index 04db3e9f3..ca1ffff28 100644 --- a/h2integrate/storage/battery/battery_performance.py +++ b/h2integrate/storage/battery/battery_performance.py @@ -1,5 +1,10 @@ import numpy as np +import pandas as pd from attrs import field, define +from openmdao.utils import units as om_units +from simses.battery.battery import Battery +from simses.degradation.state import DegradationState +from simses.model.cell.sony_lfp import SonyLFP from h2integrate.core.utilities import merge_shared_inputs from h2integrate.core.validators import gt_zero, range_val, range_val_or_none @@ -9,6 +14,19 @@ ) +class LFP280Ah(SonyLFP): + """280 Ah / 3.2 V prismatic LFP cell scaled from the SonyLFP OCV/resistance curves.""" + + _SCALE = 3.0 / 280.0 # resistance scales inversely with capacity + + def __init__(self): + super().__init__() + self.electrical.nominal_capacity = 280.0 # Ah + + def internal_resistance(self, state): + return super().internal_resistance(state) * self._SCALE + + @define(kw_only=True) class BatteryPerformanceModelConfig(StoragePerformanceBaseConfig): """Configuration class for storage performance models. @@ -158,29 +176,102 @@ def compute(self, inputs, outputs, discrete_inputs=[], discrete_outputs=[]): """Run the storage performance model.""" self.current_soc = self.config.init_soc_fraction - charge_rate = inputs["max_charge_rate"][0] + inputs["max_charge_rate"][0] if "max_discharge_rate" in inputs: discharge_rate = inputs["max_discharge_rate"][0] else: discharge_rate = inputs["max_charge_rate"][0] storage_capacity = inputs["storage_capacity"][0] - # TODO degradation: adjust compute method for degradation as needed - self.degradation(discrete_inputs["solar_resource_data"]) - # outputs[f"{self.commodity}_auxiliary_demand"] = - outputs = self.run_storage( - charge_rate, discharge_rate, storage_capacity, inputs, outputs, discrete_inputs + ### from Ankit + cell = LFP280Ah() + + battery = Battery( + cell=cell, + circuit=(239, 18), + initial_states={"start_soc": 0.5, "start_T": 25.0}, + degradation=LFP280Ah.default_degradation_model( + initial_soc=0.5, + initial_state=DegradationState(qloss_cal=1e-4), + ), + ) + + summary = pd.Series( + { + "nominal_capacity [Ah]": battery.nominal_capacity, + "nominal_voltage [V]": battery.nominal_voltage, + "nominal_energy [kWh]": battery.nominal_energy_capacity / 1e3, + "max_charge_current [A]": battery.max_charge_current, + "max_discharge_current [A]": battery.max_discharge_current, + "thermal_capacity [kJ/K]": battery.thermal_capacity / 1e3, + } + ) + print("Battery summary:") + print(summary.to_string()) + print() + + dt = self.dt_hr * 60.0 + + # power_profile = inputs[f"{self.commodity}_in"] + power_profile = inputs[f"{self.commodity}_set_point"] + + log = { + k: np.empty(self.n_timesteps) + for k in ["soc", "v", "i", "T", "loss", "heat", "soh_Q", "soh_R", "power"] + } + for i, p in enumerate(power_profile): + battery.step(float(p), dt) + for k in log: + log[k][i] = getattr(battery.state, k) + + index = pd.date_range("2025-01-01", periods=self.n_timesteps, freq=f"{int(dt)}s") + df_bat = pd.DataFrame(log, index=index) + print("\nFirst rows:") + print(df_bat.head().to_string()) + + ############# + + # Populate all OpenMDAO outputs defined in this class and its parent classes, + # pulling the time-series results from the simses battery run (``df_bat``). + + soc_ts = df_bat["soc"].to_numpy() + # Convert battery power (W) into the desired commodity rate units. + # Sign convention: positive = discharge (out of storage), negative = charge. + power_ts = om_units.convert_units( + df_bat["power"].to_numpy(), "W", self.commodity_rate_units ) - # TODO degradation: add degradation method here - def degradation( - self, - solar_resource_data, - ): - """_summary_ + # --- BatteryPerformanceModel outputs --- + outputs[f"{self.commodity}_auxiliary_demand"] = np.zeros(self.n_timesteps) - Args: - solar_resource_data (_type_): a dictionary of hourly (or by timestep) solar resource - information including irradiance and temperature - """ - return + # --- StoragePerformanceBase outputs --- + outputs["storage_duration"] = ( + storage_capacity / discharge_rate if discharge_rate > 0 else 0.0 + ) + outputs["SOC"] = soc_ts * 100.0 # fraction -> percent + outputs[f"storage_{self.commodity}_charge"] = np.where(power_ts < 0, power_ts, 0.0) + outputs[f"storage_{self.commodity}_discharge"] = np.where(power_ts > 0, power_ts, 0.0) + + # --- PerformanceModelBaseClass outputs --- + outputs[f"{self.commodity}_out"] = power_ts + outputs[f"rated_{self.commodity}_production"] = discharge_rate + outputs[f"total_{self.commodity}_produced"] = np.sum(power_ts) * self.dt_amount + outputs[f"annual_{self.commodity}_produced"] = outputs[ + f"total_{self.commodity}_produced" + ] * (1 / self.fraction_of_year_simulated) + outputs["replacement_schedule"] = np.zeros(self.plant_life) + outputs["operational_life"] = self.plant_life + + if discharge_rate <= 0: + outputs["capacity_factor"] = 0.0 + outputs["standard_capacity_factor"] = 0.0 + else: + outputs["capacity_factor"] = outputs[f"total_{self.commodity}_produced"] / ( + discharge_rate * self.n_timesteps * self.dt_amount + ) + total_commodity_discharged = ( + outputs[f"storage_{self.commodity}_discharge"].sum() * self.dt_amount + ) + outputs["standard_capacity_factor"] = total_commodity_discharged / ( + discharge_rate * self.n_timesteps * self.dt_amount + ) From f4a23b9ff1460b883c0021b1d5edaae2df9164e6 Mon Sep 17 00:00:00 2001 From: Jared Thomas Date: Fri, 17 Jul 2026 11:19:45 -0600 Subject: [PATCH 10/18] progress towards revised approach --- .../storage/battery/battery_performance.py | 155 +++++++++++++++--- 1 file changed, 134 insertions(+), 21 deletions(-) diff --git a/h2integrate/storage/battery/battery_performance.py b/h2integrate/storage/battery/battery_performance.py index ca1ffff28..77515d594 100644 --- a/h2integrate/storage/battery/battery_performance.py +++ b/h2integrate/storage/battery/battery_performance.py @@ -1,10 +1,40 @@ +import math + import numpy as np import pandas as pd from attrs import field, define from openmdao.utils import units as om_units +from simses.degradation import DegradationModel +from simses.battery.state import BatteryState from simses.battery.battery import Battery +from simses.thermal.ambient import AmbientThermalModel from simses.degradation.state import DegradationState from simses.model.cell.sony_lfp import SonyLFP +from simses.degradation.cycle_detector import HalfCycle +from simses.model.degradation.sony_lfp_cyclic import ( + A_RINC, + B_RINC, + C_RINC as CYC_C_RINC, + D_RINC as CYC_D_RINC, + A_QLOSS, + B_QLOSS, + C_QLOSS as CYC_C_QLOSS, + D_QLOSS as CYC_D_QLOSS, + SonyLFPCyclicDegradation, +) +from simses.model.degradation.sony_lfp_calendar import ( + T_REF, + C_RINC as CAL_C_RINC, + D_RINC as CAL_D_RINC, + C_QLOSS as CAL_C_QLOSS, + D_QLOSS as CAL_D_QLOSS, + EA_RINC, + EA_QLOSS, + K_REF_RINC, + K_REF_QLOSS, + R, + SonyLFPCalendarDegradation, +) from h2integrate.core.utilities import merge_shared_inputs from h2integrate.core.validators import gt_zero, range_val, range_val_or_none @@ -22,11 +52,76 @@ class LFP280Ah(SonyLFP): def __init__(self): super().__init__() self.electrical.nominal_capacity = 280.0 # Ah + # Thermal properties for a large-format prismatic cell (vs. the 70 g 26650 reference). + # mass=1.5 kg, h=23 W/m²K gives C_th≈6.5 MJ/K, R_th≈1.73 mK/W, τ≈3.1 h + # → ΔT ≈ 10 °C at end of 2-hour C/2 discharge. + self.thermal.mass = 3.0 # kg per cell + self.thermal.convection_coefficient = 23.0 # W/m²K def internal_resistance(self, state): return super().internal_resistance(state) * self._SCALE +class ScaledLFPCalendarDegradation(SonyLFPCalendarDegradation): + def update_capacity( + self, state: BatteryState, dt: float, accumulated_qloss: float, _DEG_SCALE + ) -> float: + if dt == 0.0: + return 0.0 + T_K = state.T + 273.15 + T_REF_K = T_REF + 273.15 + k_T_q = (K_REF_QLOSS * _DEG_SCALE) * math.exp(-EA_QLOSS / R * (1.0 / T_K - 1.0 / T_REF_K)) + k_soc_q = CAL_C_QLOSS * (state.soc - 0.5) ** 3 + CAL_D_QLOSS + stress_q = k_T_q * k_soc_q + if stress_q > 0.0: + virtual_time = (accumulated_qloss / stress_q) ** 2 + delta_q = stress_q * math.sqrt(virtual_time + dt) - accumulated_qloss + else: + delta_q = 0.0 + return delta_q + + def update_resistance(self, state: BatteryState, dt: float, _DEG_SCALE) -> float: + if dt == 0.0: + return 0.0 + T_K = state.T + 273.15 + T_REF_K = T_REF + 273.15 + k_T_r = (K_REF_RINC * _DEG_SCALE) * math.exp(-EA_RINC / R * (1.0 / T_K - 1.0 / T_REF_K)) + k_soc_r = CAL_C_RINC * (state.soc - 0.5) ** 2 + CAL_D_RINC + return k_T_r * k_soc_r * dt + + +class ScaledLFPCyclicDegradation(SonyLFPCyclicDegradation): + def update_capacity( + self, + state: BatteryState, + half_cycle: HalfCycle, + accumulated_qloss: float, + _DEG_SCALE: float, + ) -> float: + delta_fec = half_cycle.full_equivalent_cycles + if delta_fec == 0.0: + return 0.0 + k_crate_q = (A_QLOSS * _DEG_SCALE) * half_cycle.c_rate + (B_QLOSS * _DEG_SCALE) + k_dod_q = CYC_C_QLOSS * (half_cycle.depth_of_discharge - 0.6) ** 3 + CYC_D_QLOSS + stress_q = k_crate_q * k_dod_q + if stress_q > 0.0: + virtual_fec = (accumulated_qloss * 100.0 / stress_q) ** 2 + delta_q = stress_q * math.sqrt(virtual_fec + delta_fec) / 100.0 - accumulated_qloss + else: + delta_q = 0.0 + return delta_q + + def update_resistance( + self, state: BatteryState, half_cycle: HalfCycle, _DEG_SCALE: float + ) -> float: + delta_fec = half_cycle.full_equivalent_cycles + if delta_fec == 0.0: + return 0.0 + k_crate_r = (A_RINC * _DEG_SCALE) * half_cycle.c_rate + (B_RINC * _DEG_SCALE) + k_dod_r = CYC_C_RINC * (half_cycle.depth_of_discharge - 0.5) ** 3 + CYC_D_RINC + return k_crate_r * k_dod_r * delta_fec / 100.0 + + @define(kw_only=True) class BatteryPerformanceModelConfig(StoragePerformanceBaseConfig): """Configuration class for storage performance models. @@ -86,6 +181,8 @@ class BatteryPerformanceModelConfig(StoragePerformanceBaseConfig): discharge_efficiency: float | None = field(default=None, validator=range_val_or_none(0, 1)) round_trip_efficiency: float | None = field(default=None, validator=range_val_or_none(0, 1)) + _DEG_SCALE: float = field(default=0.4, validator=range_val(0, 1)) + # TODO degradation: add additional parameters for degradation here cop: float = field(validator=gt_zero) @@ -183,15 +280,25 @@ def compute(self, inputs, outputs, discrete_inputs=[], discrete_outputs=[]): discharge_rate = inputs["max_charge_rate"][0] storage_capacity = inputs["storage_capacity"][0] + power_profile = inputs[f"{self.commodity}_command_value"] ### from Ankit + + # --------------------------------------------------------------------------- + # Battery pack: 239s x 18p -> 764.8 V * 5040 Ah ~ 3855 kWh + # --------------------------------------------------------------------------- cell = LFP280Ah() battery = Battery( cell=cell, - circuit=(239, 18), - initial_states={"start_soc": 0.5, "start_T": 25.0}, - degradation=LFP280Ah.default_degradation_model( - initial_soc=0.5, + circuit=(239, 18), # TODO update to be based on provided battery power and energy + initial_states={ + "start_soc": self.config.init_soc_fraction, + "start_T": 25.0, + }, # TODO should be user inputs + degradation=DegradationModel( + calendar=ScaledLFPCalendarDegradation(), + cyclic=ScaledLFPCyclicDegradation(), + initial_soc=self.config.init_soc_fraction, initial_state=DegradationState(qloss_cal=1e-4), ), ) @@ -210,36 +317,42 @@ def compute(self, inputs, outputs, discrete_inputs=[], discrete_outputs=[]): print(summary.to_string()) print() - dt = self.dt_hr * 60.0 + # --------------------------------------------------------------------------- + # Thermal model: constant 25 °C ambient, battery registered as thermal node + # --------------------------------------------------------------------------- + thermal = AmbientThermalModel(T_ambient=25.0, components=[battery]) - # power_profile = inputs[f"{self.commodity}_in"] - power_profile = inputs[f"{self.commodity}_set_point"] + # --------------------------------------------------------------------------- + # Simulation loop + # --------------------------------------------------------------------------- + keys = ["soc", "v", "i", "T", "loss", "heat", "soh_Q", "soh_R", "power"] + log = {k: np.empty(self.n_timesteps) for k in keys} - log = { - k: np.empty(self.n_timesteps) - for k in ["soc", "v", "i", "T", "loss", "heat", "soh_Q", "soh_R", "power"] - } for i, p in enumerate(power_profile): - battery.step(float(p), dt) - for k in log: + battery.step(float(p), self.dt) + thermal.step(self.dt) # update battery temperature after each step + for k in keys: log[k][i] = getattr(battery.state, k) - index = pd.date_range("2025-01-01", periods=self.n_timesteps, freq=f"{int(dt)}s") - df_bat = pd.DataFrame(log, index=index) - print("\nFirst rows:") - print(df_bat.head().to_string()) + index = pd.date_range("2026-01-01", periods=self.n_timesteps, freq=f"{int(self.dt)}s") + df = pd.DataFrame(log, index=index) + print("\nFirst rows:") + print(df.head().to_string()) + # print( + # f"\nFinal SOH_Q : {df['soh_Q'].iloc[-1]:.4f} ({(1 - df['soh_Q'].iloc[-1]) + # * 100:.2f} % capacity fade)" + # ) + print(f"Final SOH_R : {df['soh_R'].iloc[-1]:.4f}") ############# # Populate all OpenMDAO outputs defined in this class and its parent classes, # pulling the time-series results from the simses battery run (``df_bat``). - soc_ts = df_bat["soc"].to_numpy() + soc_ts = df["soc"].to_numpy() # Convert battery power (W) into the desired commodity rate units. # Sign convention: positive = discharge (out of storage), negative = charge. - power_ts = om_units.convert_units( - df_bat["power"].to_numpy(), "W", self.commodity_rate_units - ) + power_ts = om_units.convert_units(df["power"].to_numpy(), "W", self.commodity_rate_units) # --- BatteryPerformanceModel outputs --- outputs[f"{self.commodity}_auxiliary_demand"] = np.zeros(self.n_timesteps) From 77c8310e39b1286db7350efd62a55067d314348b Mon Sep 17 00:00:00 2001 From: Jared Thomas Date: Mon, 20 Jul 2026 15:45:03 -0600 Subject: [PATCH 11/18] including batt deg example --- .../run_battery_degradation.py | 322 ++++++++++++++++++ 1 file changed, 322 insertions(+) create mode 100644 examples/98_battery_degradation/run_battery_degradation.py diff --git a/examples/98_battery_degradation/run_battery_degradation.py b/examples/98_battery_degradation/run_battery_degradation.py new file mode 100644 index 000000000..507305f2f --- /dev/null +++ b/examples/98_battery_degradation/run_battery_degradation.py @@ -0,0 +1,322 @@ +import math + +import matplotlib + + +matplotlib.use("Agg") +from pathlib import Path + +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt +from tqdm import tqdm + + +HERE = Path(__file__).parent + +from simses.degradation import DegradationModel +from simses.battery.state import BatteryState +from simses.battery.battery import Battery +from simses.thermal.ambient import AmbientThermalModel +from simses.degradation.state import DegradationState +from simses.model.cell.sony_lfp import SonyLFP +from simses.degradation.cycle_detector import HalfCycle +from simses.model.degradation.sony_lfp_cyclic import ( + A_RINC, + B_RINC, + C_RINC as CYC_C_RINC, + D_RINC as CYC_D_RINC, + A_QLOSS, + B_QLOSS, + C_QLOSS as CYC_C_QLOSS, + D_QLOSS as CYC_D_QLOSS, + SonyLFPCyclicDegradation, +) +from simses.model.degradation.sony_lfp_calendar import ( + T_REF, + C_RINC as CAL_C_RINC, + D_RINC as CAL_D_RINC, + C_QLOSS as CAL_C_QLOSS, + D_QLOSS as CAL_D_QLOSS, + EA_RINC, + EA_QLOSS, + K_REF_RINC, + K_REF_QLOSS, + R, + SonyLFPCalendarDegradation, +) + + +# --------------------------------------------------------------------------- +# Large-format 280 Ah LFP cell +# --------------------------------------------------------------------------- + + +class LFP280Ah(SonyLFP): + """280 Ah / 3.2 V prismatic LFP cell scaled from the SonyLFP OCV/resistance curves.""" + + _SCALE = 3.0 / 280.0 # resistance scales inversely with capacity + + def __init__(self): + super().__init__() + self.electrical.nominal_capacity = 280.0 # Ah + # Thermal properties for a large-format prismatic cell (vs. the 70 g 26650 reference). + # mass=1.5 kg, h=23 W/m²K gives C_th≈6.5 MJ/K, R_th≈1.73 mK/W, τ≈3.1 h + # → ΔT ≈ 10 °C at end of 2-hour C/2 discharge. + self.thermal.mass = 3.0 # kg per cell + self.thermal.convection_coefficient = 23.0 # W/m²K + + def internal_resistance(self, state): + return super().internal_resistance(state) * self._SCALE + + +# --------------------------------------------------------------------------- +# Scaled degradation models targeting ~1-2 % capacity fade over 365 days. +# The Sony LFP parametrization is for a 3 Ah cylindrical cell; large-format +# prismatic LFP cells degrade roughly 5x slower per calendar/cycle unit. +# A uniform scale factor of 0.2 brings the total fade into the 1-2 % range. +# --------------------------------------------------------------------------- +_DEG_SCALE = 0.4 + + +class ScaledLFPCalendarDegradation(SonyLFPCalendarDegradation): + def update_capacity(self, state: BatteryState, dt: float, accumulated_qloss: float) -> float: + if dt == 0.0: + return 0.0 + T_K = state.T + 273.15 + T_REF_K = T_REF + 273.15 + k_T_q = (K_REF_QLOSS * _DEG_SCALE) * math.exp(-EA_QLOSS / R * (1.0 / T_K - 1.0 / T_REF_K)) + k_soc_q = CAL_C_QLOSS * (state.soc - 0.5) ** 3 + CAL_D_QLOSS + stress_q = k_T_q * k_soc_q + if stress_q > 0.0: + virtual_time = (accumulated_qloss / stress_q) ** 2 + delta_q = stress_q * math.sqrt(virtual_time + dt) - accumulated_qloss + else: + delta_q = 0.0 + return delta_q + + def update_resistance(self, state: BatteryState, dt: float) -> float: + if dt == 0.0: + return 0.0 + T_K = state.T + 273.15 + T_REF_K = T_REF + 273.15 + k_T_r = (K_REF_RINC * _DEG_SCALE) * math.exp(-EA_RINC / R * (1.0 / T_K - 1.0 / T_REF_K)) + k_soc_r = CAL_C_RINC * (state.soc - 0.5) ** 2 + CAL_D_RINC + return k_T_r * k_soc_r * dt + + +class ScaledLFPCyclicDegradation(SonyLFPCyclicDegradation): + def update_capacity( + self, state: BatteryState, half_cycle: HalfCycle, accumulated_qloss: float + ) -> float: + delta_fec = half_cycle.full_equivalent_cycles + if delta_fec == 0.0: + return 0.0 + k_crate_q = (A_QLOSS * _DEG_SCALE) * half_cycle.c_rate + (B_QLOSS * _DEG_SCALE) + k_dod_q = CYC_C_QLOSS * (half_cycle.depth_of_discharge - 0.6) ** 3 + CYC_D_QLOSS + stress_q = k_crate_q * k_dod_q + if stress_q > 0.0: + virtual_fec = (accumulated_qloss * 100.0 / stress_q) ** 2 + delta_q = stress_q * math.sqrt(virtual_fec + delta_fec) / 100.0 - accumulated_qloss + else: + delta_q = 0.0 + return delta_q + + def update_resistance(self, state: BatteryState, half_cycle: HalfCycle) -> float: + delta_fec = half_cycle.full_equivalent_cycles + if delta_fec == 0.0: + return 0.0 + k_crate_r = (A_RINC * _DEG_SCALE) * half_cycle.c_rate + (B_RINC * _DEG_SCALE) + k_dod_r = CYC_C_RINC * (half_cycle.depth_of_discharge - 0.5) ** 3 + CYC_D_RINC + return k_crate_r * k_dod_r * delta_fec / 100.0 + + +# --------------------------------------------------------------------------- +# Battery pack: 239s x 18p -> 764.8 V * 5040 Ah ~ 3855 kWh +# --------------------------------------------------------------------------- +cell = LFP280Ah() + +battery = Battery( + cell=cell, + circuit=(239, 18), + initial_states={"start_soc": 1.0, "start_T": 25.0}, + degradation=DegradationModel( + calendar=ScaledLFPCalendarDegradation(), + cyclic=ScaledLFPCyclicDegradation(), + initial_soc=1.0, + initial_state=DegradationState(qloss_cal=1e-4), + ), +) + +summary = pd.Series( + { + "nominal_capacity [Ah]": battery.nominal_capacity, + "nominal_voltage [V]": battery.nominal_voltage, + "nominal_energy [kWh]": battery.nominal_energy_capacity / 1e3, + "max_charge_current [A]": battery.max_charge_current, + "max_discharge_current [A]": battery.max_discharge_current, + "thermal_capacity [kJ/K]": battery.thermal_capacity / 1e3, + } +) +print("Battery summary:") +print(summary.to_string()) +print() + +# --------------------------------------------------------------------------- +# 15-year power demand profile (dt = 15 min = 900 s) +# +# Each day (96 steps): +# 8 steps discharge at C/2 (~-1.927 MW) +# 40 steps charge at C/10 (~+385 kW) +# 48 steps rest 0 W +# --------------------------------------------------------------------------- +dt = 900.0 +STEPS_PER_DAY = 96 +N_YEARS = 15 +N_DAYS = N_YEARS * 365 +n_steps = N_DAYS * STEPS_PER_DAY + +P_nom = battery.nominal_energy_capacity # Wh +P_disch = -0.5 * P_nom # C/2 discharge +P_chg = +0.1 * P_nom # C/10 charge + +day_profile = np.concatenate( + [ + np.full(8, P_disch), + np.full(40, P_chg), + np.full(48, 0.0), + ] +) +power_profile = np.tile(day_profile, N_DAYS) + +print(f"Simulation: {N_YEARS} years ({N_DAYS} days), {n_steps:,} steps, dt={int(dt)} s") +print(f"Discharge: {P_disch/1e6:.3f} MW (C/2)") +print(f"Charge: {P_chg/1e3:.1f} kW (C/10)") +print() + +# --------------------------------------------------------------------------- +# Thermal model: constant 25 °C ambient, battery registered as thermal node +# --------------------------------------------------------------------------- +thermal = AmbientThermalModel(T_ambient=25.0, components=[battery]) + +# --------------------------------------------------------------------------- +# Simulation loop +# --------------------------------------------------------------------------- +keys = ["soc", "v", "i", "T", "loss", "heat", "soh_Q", "soh_R", "power"] +log = {k: np.empty(n_steps) for k in keys} + +for i, p in enumerate(tqdm(power_profile, desc=f"Simulating {N_YEARS} years")): + battery.step(float(p), dt) + thermal.step(dt) # update battery temperature after each step + for k in keys: + log[k][i] = getattr(battery.state, k) + +index = pd.date_range("2026-01-01", periods=n_steps, freq=f"{int(dt)}s") +df = pd.DataFrame(log, index=index) + +print("\nFirst rows:") +print(df.head().to_string()) +print( + f"\nFinal SOH_Q : {df['soh_Q'].iloc[-1]:.4f} \ + ({(1 - df['soh_Q'].iloc[-1]) * 100:.2f} % capacity fade)" +) +print(f"Final SOH_R : {df['soh_R'].iloc[-1]:.4f}") + +# --------------------------------------------------------------------------- +# Plots +# --------------------------------------------------------------------------- + +# Plot 1: Power demand (first 7 days) +fig, ax = plt.subplots(figsize=(12, 3)) +df["power"].iloc[: 7 * STEPS_PER_DAY].plot(ax=ax) +ax.set_title("Power demand profile -- first 7 days") +ax.set_ylabel("Power [W]") +ax.axhline(0, color="gray", linewidth=0.5) +ax.yaxis.set_major_formatter(plt.FuncFormatter(lambda x, _: f"{x/1e6:.2f} MW")) +fig.tight_layout() +fig.savefig(HERE / "plot_power.png", dpi=150) +print("Saved plot_power.png") + +# Plot 2: SOC (hourly mean) +fig, ax = plt.subplots(figsize=(12, 3)) +df["soc"].resample("1h").mean().plot(ax=ax, title="State of Charge -- hourly mean") +ax.set_ylabel("SOC [p.u.]") +fig.tight_layout() +fig.savefig(HERE / "plot_soc.png", dpi=150) +print("Saved plot_soc.png") + +# Plot 2b: Battery temperature -- first 7 days (15-min resolution) +fig, ax = plt.subplots(figsize=(12, 3)) +df["T"].iloc[: 7 * STEPS_PER_DAY].plot(ax=ax, title="Battery temperature -- first 7 days") +ax.set_ylabel("Temperature [°C]") +ax.axhline(25.0, color="gray", linewidth=0.8, linestyle="--", label="T_ambient = 25 °C") +ax.legend() +fig.tight_layout() +fig.savefig(HERE / "plot_temperature_7days.png", dpi=150) +print("Saved plot_temperature_7days.png") + +# Plot 3a: Terminal voltage -- first 7 days (raw 15-min data) +fig, ax = plt.subplots(figsize=(12, 3)) +df["v"].iloc[: 7 * STEPS_PER_DAY].plot(ax=ax, title="Terminal voltage -- first 7 days [V]") +ax.set_ylabel("Voltage [V]") +fig.tight_layout() +fig.savefig(HERE / "plot_voltage_7days.png", dpi=150) +print("Saved plot_voltage_7days.png") + +# Plot 3b: Terminal voltage (daily mean, full year) +fig, ax = plt.subplots(figsize=(12, 3)) +df["v"].resample("1D").mean().plot(ax=ax, title="Terminal voltage -- daily mean [V]") +ax.set_ylabel("Voltage [V]") +fig.tight_layout() +fig.savefig(HERE / "plot_voltage.png", dpi=150) +print("Saved plot_voltage.png") + +# Plot 4: Degradation over 15 years +df_monthly = df[["soh_Q", "soh_R"]].resample("ME").last() +capacity_fade_pct = (1 - df_monthly["soh_Q"]) * 100 +resistance_growth_pct = (df_monthly["soh_R"] - 1) * 100 + +fig, ax1 = plt.subplots(figsize=(12, 4)) +ax2 = ax1.twinx() + +ax1.plot( + capacity_fade_pct.index, capacity_fade_pct.values, color="steelblue", label="Capacity fade [%]" +) +ax2.plot( + resistance_growth_pct.index, + resistance_growth_pct.values, + color="darkorange", + label="Resistance growth [%]", +) + +ax1.set_ylabel("Capacity fade [%]", color="steelblue") +ax2.set_ylabel("Resistance growth [%]", color="darkorange") +ax1.tick_params(axis="y", labelcolor="steelblue") +ax2.tick_params(axis="y", labelcolor="darkorange") +ax1.set_title( + f"Battery degradation over {N_YEARS} years " + f"(final: {capacity_fade_pct.iloc[-1]:.1f} % capacity fade, " + f"{resistance_growth_pct.iloc[-1]:.1f} % resistance growth)" +) +ax1.set_xlabel("Date") + +lines1, labels1 = ax1.get_legend_handles_labels() +lines2, labels2 = ax2.get_legend_handles_labels() +ax1.legend(lines1 + lines2, labels1 + labels2, loc="upper left") + +# Mark end-of-life reference (20 % capacity fade) +ax1.axhline(20, color="steelblue", linewidth=0.8, linestyle="--", alpha=0.6, label="EoL (20 %)") + +fig.tight_layout() +fig.savefig(HERE / "plot_degradation_15yr.png", dpi=150) +print("Saved plot_degradation_15yr.png") + +# Plot 5: Losses and heat (daily mean) +fig, ax = plt.subplots(figsize=(12, 3)) +df[["loss", "heat"]].resample("1D").mean().plot( + ax=ax, title="Battery losses and heat -- daily mean [W]" +) +ax.set_ylabel("Power [W]") +fig.tight_layout() +fig.savefig(HERE / "plot_losses.png", dpi=150) +print("Saved plot_losses.png") From f3c380a048137359ed84fc0cd1519a67f7905633 Mon Sep 17 00:00:00 2001 From: Jared Thomas Date: Thu, 13 Aug 2026 14:09:38 -0600 Subject: [PATCH 12/18] update with converter Co-authored-by: verank <9517606+verank@users.noreply.github.com> --- .../storage/battery/battery_performance.py | 136 ++++++++++-------- 1 file changed, 80 insertions(+), 56 deletions(-) diff --git a/h2integrate/storage/battery/battery_performance.py b/h2integrate/storage/battery/battery_performance.py index 77515d594..99ef3c6d0 100644 --- a/h2integrate/storage/battery/battery_performance.py +++ b/h2integrate/storage/battery/battery_performance.py @@ -1,7 +1,6 @@ import math import numpy as np -import pandas as pd from attrs import field, define from openmdao.utils import units as om_units from simses.degradation import DegradationModel @@ -9,8 +8,10 @@ from simses.battery.battery import Battery from simses.thermal.ambient import AmbientThermalModel from simses.degradation.state import DegradationState +from simses.converter.converter import Converter from simses.model.cell.sony_lfp import SonyLFP from simses.degradation.cycle_detector import HalfCycle +from simses.model.converter.fix_efficiency import FixedEfficiency from simses.model.degradation.sony_lfp_cyclic import ( A_RINC, B_RINC, @@ -45,9 +46,15 @@ class LFP280Ah(SonyLFP): - """280 Ah / 3.2 V prismatic LFP cell scaled from the SonyLFP OCV/resistance curves.""" + """280 Ah / 3.2 V prismatic LFP cell scaled from the SonyLFP OCV/resistance curves. - _SCALE = 3.0 / 280.0 # resistance scales inversely with capacity + Resistance scaling: + Step 1 - capacity scaling: R scales as 1/Q, so the first factor is 3/280. + Step 2 - design correction for large-format prismatic multi-tab cells. + Combined: _SCALE = 0.003888, matching about 0.18 mOhm at SOC=0.5, T=25 C. + """ + + _SCALE = 0.18e-3 / ((0.044767041 + 0.047827935) / 2) def __init__(self): super().__init__() @@ -63,14 +70,18 @@ def internal_resistance(self, state): class ScaledLFPCalendarDegradation(SonyLFPCalendarDegradation): - def update_capacity( - self, state: BatteryState, dt: float, accumulated_qloss: float, _DEG_SCALE - ) -> float: + def __init__(self, deg_scale: float): + super().__init__() + self._deg_scale = deg_scale + + def update_capacity(self, state: BatteryState, dt: float, accumulated_qloss: float) -> float: if dt == 0.0: return 0.0 T_K = state.T + 273.15 T_REF_K = T_REF + 273.15 - k_T_q = (K_REF_QLOSS * _DEG_SCALE) * math.exp(-EA_QLOSS / R * (1.0 / T_K - 1.0 / T_REF_K)) + k_T_q = (K_REF_QLOSS * self._deg_scale) * math.exp( + -EA_QLOSS / R * (1.0 / T_K - 1.0 / T_REF_K) + ) k_soc_q = CAL_C_QLOSS * (state.soc - 0.5) ** 3 + CAL_D_QLOSS stress_q = k_T_q * k_soc_q if stress_q > 0.0: @@ -80,28 +91,33 @@ def update_capacity( delta_q = 0.0 return delta_q - def update_resistance(self, state: BatteryState, dt: float, _DEG_SCALE) -> float: + def update_resistance(self, state: BatteryState, dt: float) -> float: if dt == 0.0: return 0.0 T_K = state.T + 273.15 T_REF_K = T_REF + 273.15 - k_T_r = (K_REF_RINC * _DEG_SCALE) * math.exp(-EA_RINC / R * (1.0 / T_K - 1.0 / T_REF_K)) + k_T_r = (K_REF_RINC * self._deg_scale) * math.exp( + -EA_RINC / R * (1.0 / T_K - 1.0 / T_REF_K) + ) k_soc_r = CAL_C_RINC * (state.soc - 0.5) ** 2 + CAL_D_RINC return k_T_r * k_soc_r * dt class ScaledLFPCyclicDegradation(SonyLFPCyclicDegradation): + def __init__(self, deg_scale: float): + super().__init__() + self._deg_scale = deg_scale + def update_capacity( self, state: BatteryState, half_cycle: HalfCycle, accumulated_qloss: float, - _DEG_SCALE: float, ) -> float: delta_fec = half_cycle.full_equivalent_cycles if delta_fec == 0.0: return 0.0 - k_crate_q = (A_QLOSS * _DEG_SCALE) * half_cycle.c_rate + (B_QLOSS * _DEG_SCALE) + k_crate_q = (A_QLOSS * self._deg_scale) * half_cycle.c_rate + (B_QLOSS * self._deg_scale) k_dod_q = CYC_C_QLOSS * (half_cycle.depth_of_discharge - 0.6) ** 3 + CYC_D_QLOSS stress_q = k_crate_q * k_dod_q if stress_q > 0.0: @@ -111,13 +127,11 @@ def update_capacity( delta_q = 0.0 return delta_q - def update_resistance( - self, state: BatteryState, half_cycle: HalfCycle, _DEG_SCALE: float - ) -> float: + def update_resistance(self, state: BatteryState, half_cycle: HalfCycle) -> float: delta_fec = half_cycle.full_equivalent_cycles if delta_fec == 0.0: return 0.0 - k_crate_r = (A_RINC * _DEG_SCALE) * half_cycle.c_rate + (B_RINC * _DEG_SCALE) + k_crate_r = (A_RINC * self._deg_scale) * half_cycle.c_rate + (B_RINC * self._deg_scale) k_dod_r = CYC_C_RINC * (half_cycle.depth_of_discharge - 0.5) ** 3 + CYC_D_RINC return k_crate_r * k_dod_r * delta_fec / 100.0 @@ -181,7 +195,13 @@ class BatteryPerformanceModelConfig(StoragePerformanceBaseConfig): discharge_efficiency: float | None = field(default=None, validator=range_val_or_none(0, 1)) round_trip_efficiency: float | None = field(default=None, validator=range_val_or_none(0, 1)) - _DEG_SCALE: float = field(default=0.4, validator=range_val(0, 1)) + _DEG_SCALE: float = field(default=0.7056, validator=range_val(0, 1)) + # TODO convert from power and energy ratings (see math in chat) + series_count: int = field(default=336, converter=int, validator=gt_zero) + parallel_count: int = field(default=16, converter=int, validator=gt_zero) + battery_temperature_c: float = field(default=25.0) + converter_efficiency: float = field(default=0.96, validator=range_val(0, 1)) + converter_max_power: float = field(default=2400.0, validator=gt_zero) # TODO degradation: add additional parameters for degradation here cop: float = field(validator=gt_zero) @@ -238,10 +258,16 @@ class BatteryPerformanceModel(StoragePerformanceBase): """OpenMDAO component for a storage component.""" _time_step_bounds = ( - 1, + 60, 3600, ) # (min, max) time step lengths (in seconds) compatible with this model + def initialize(self): + super().initialize() + self.commodity = "electricity" + self.commodity_rate_units = "kW" + self.commodity_amount_units = "kW*h" + def setup(self): self.config = BatteryPerformanceModelConfig.from_dict( merge_shared_inputs(self.options["tech_config"]["model_inputs"], "performance"), @@ -280,81 +306,79 @@ def compute(self, inputs, outputs, discrete_inputs=[], discrete_outputs=[]): discharge_rate = inputs["max_charge_rate"][0] storage_capacity = inputs["storage_capacity"][0] + # H2I dispatch command: positive = discharge, negative = charge (commodity_rate_units) power_profile = inputs[f"{self.commodity}_command_value"] + ### from Ankit # --------------------------------------------------------------------------- - # Battery pack: 239s x 18p -> 764.8 V * 5040 Ah ~ 3855 kWh + # Battery pack + inverter (fixed Megapack-style topology, from config) # --------------------------------------------------------------------------- cell = LFP280Ah() + # TODO check sizing battery = Battery( cell=cell, - circuit=(239, 18), # TODO update to be based on provided battery power and energy + circuit=(self.config.series_count, self.config.parallel_count), initial_states={ "start_soc": self.config.init_soc_fraction, - "start_T": 25.0, - }, # TODO should be user inputs + "start_T": self.config.battery_temperature_c, + }, degradation=DegradationModel( - calendar=ScaledLFPCalendarDegradation(), - cyclic=ScaledLFPCyclicDegradation(), + calendar=ScaledLFPCalendarDegradation(self.config._DEG_SCALE), + cyclic=ScaledLFPCyclicDegradation(self.config._DEG_SCALE), initial_soc=self.config.init_soc_fraction, initial_state=DegradationState(qloss_cal=1e-4), ), ) - summary = pd.Series( - { - "nominal_capacity [Ah]": battery.nominal_capacity, - "nominal_voltage [V]": battery.nominal_voltage, - "nominal_energy [kWh]": battery.nominal_energy_capacity / 1e3, - "max_charge_current [A]": battery.max_charge_current, - "max_discharge_current [A]": battery.max_discharge_current, - "thermal_capacity [kJ/K]": battery.thermal_capacity / 1e3, - } + converter = Converter( + loss_model=FixedEfficiency(self.config.converter_efficiency), + max_power=om_units.convert_units( + self.config.converter_max_power, self.commodity_rate_units, "W" + ), + storage=battery, ) - print("Battery summary:") - print(summary.to_string()) - print() # --------------------------------------------------------------------------- - # Thermal model: constant 25 °C ambient, battery registered as thermal node + # Thermal model: constant ambient, battery registered as thermal node + # #TODO check battery temp ambient # --------------------------------------------------------------------------- - thermal = AmbientThermalModel(T_ambient=25.0, components=[battery]) + thermal = AmbientThermalModel( + T_ambient=self.config.battery_temperature_c, components=[battery] + ) # --------------------------------------------------------------------------- # Simulation loop # --------------------------------------------------------------------------- - keys = ["soc", "v", "i", "T", "loss", "heat", "soh_Q", "soh_R", "power"] + keys = ["soc", "v", "i", "T", "loss", "heat", "soh_Q", "soh_R"] log = {k: np.empty(self.n_timesteps) for k in keys} + power_ac = np.empty(self.n_timesteps) + power_dc = np.empty(self.n_timesteps) + conv_loss = np.empty(self.n_timesteps) for i, p in enumerate(power_profile): - battery.step(float(p), self.dt) + # H2I sign (+discharge) -> SimSES AC sign (+charge) + converter.step( + -om_units.convert_units(float(p), self.commodity_rate_units, "W"), self.dt + ) thermal.step(self.dt) # update battery temperature after each step for k in keys: log[k][i] = getattr(battery.state, k) + power_ac[i] = converter.state.power # AC power (W), positive = charge + power_dc[i] = battery.state.power # AC power (W), positive = charge + conv_loss[i] = converter.state.loss - index = pd.date_range("2026-01-01", periods=self.n_timesteps, freq=f"{int(self.dt)}s") - df = pd.DataFrame(log, index=index) - - print("\nFirst rows:") - print(df.head().to_string()) - # print( - # f"\nFinal SOH_Q : {df['soh_Q'].iloc[-1]:.4f} ({(1 - df['soh_Q'].iloc[-1]) - # * 100:.2f} % capacity fade)" - # ) - print(f"Final SOH_R : {df['soh_R'].iloc[-1]:.4f}") ############# - # Populate all OpenMDAO outputs defined in this class and its parent classes, - # pulling the time-series results from the simses battery run (``df_bat``). - - soc_ts = df["soc"].to_numpy() - # Convert battery power (W) into the desired commodity rate units. - # Sign convention: positive = discharge (out of storage), negative = charge. - power_ts = om_units.convert_units(df["power"].to_numpy(), "W", self.commodity_rate_units) + # Populate all OpenMDAO outputs defined in this class and its parent classes. + # Convert SimSES AC power (W, +charge) back to H2I convention + # (commodity_rate_units, +discharge). + soc_ts = log["soc"] + power_ts = -om_units.convert_units(power_ac, "W", self.commodity_rate_units) # --- BatteryPerformanceModel outputs --- + # TODO calc aux power outputs[f"{self.commodity}_auxiliary_demand"] = np.zeros(self.n_timesteps) # --- StoragePerformanceBase outputs --- From dcf89cb62bc33484713f0f2d268f8a211ba7ce3e Mon Sep 17 00:00:00 2001 From: Jared Thomas Date: Fri, 14 Aug 2026 08:58:37 -0600 Subject: [PATCH 13/18] working example running battery degradation in h2i interface --- environment.yml | 2 +- .../run_battery_degradation.py | 359 +++++++----------- .../storage/battery/battery_performance.py | 20 +- 3 files changed, 151 insertions(+), 230 deletions(-) diff --git a/environment.yml b/environment.yml index 5ca9dbed3..0995d8b00 100644 --- a/environment.yml +++ b/environment.yml @@ -1,7 +1,7 @@ name: h2integrate channels: [conda-forge, defaults] dependencies: - - python=3.11 + - python=3.13 - glpk - coin-or-cbc>=2.10.12 # - wisdem # NOTE: uncomment if installing with Ard diff --git a/examples/98_battery_degradation/run_battery_degradation.py b/examples/98_battery_degradation/run_battery_degradation.py index 507305f2f..53e4a3017 100644 --- a/examples/98_battery_degradation/run_battery_degradation.py +++ b/examples/98_battery_degradation/run_battery_degradation.py @@ -1,269 +1,173 @@ -import math +""" +Example 98: Battery degradation through the standard H2Integrate interface. + +This example runs the SimSES-backed ``BatteryPerformanceModel`` through the standard +``H2IntegrateModel`` interface (YAML-configured plant/tech/driver). The battery has no +control strategy, so a passthrough controller forwards the ``electricity_set_point`` +profile to the battery as its dispatch command (positive = discharge, negative = charge). + +A repeating daily charge/discharge cycle is applied for the full simulation horizon and +the battery internal timeseries (voltage, temperature, state-of-health, losses) are read +from ``BatteryPerformanceModel.results`` to produce the reference plots: + + 1. AC and DC power (first 7 days) + 2. State of charge (hourly mean) + 3. Battery temperature (first 7 days) + 4. Terminal voltage (first 7 days and daily mean) + 5. Capacity fade and resistance growth over the horizon + 6. Battery/converter losses and heat (daily mean) + +H2Integrate constrains a single simulation to exactly one year +(``n_timesteps * dt == 31_536_000 s``), so this example runs a 1-year degradation +study. The battery pack topology and degradation scaling live in ``tech_config.yaml``. +""" -import matplotlib - - -matplotlib.use("Agg") from pathlib import Path import numpy as np import pandas as pd -import matplotlib.pyplot as plt -from tqdm import tqdm - - -HERE = Path(__file__).parent - -from simses.degradation import DegradationModel -from simses.battery.state import BatteryState -from simses.battery.battery import Battery -from simses.thermal.ambient import AmbientThermalModel -from simses.degradation.state import DegradationState -from simses.model.cell.sony_lfp import SonyLFP -from simses.degradation.cycle_detector import HalfCycle -from simses.model.degradation.sony_lfp_cyclic import ( - A_RINC, - B_RINC, - C_RINC as CYC_C_RINC, - D_RINC as CYC_D_RINC, - A_QLOSS, - B_QLOSS, - C_QLOSS as CYC_C_QLOSS, - D_QLOSS as CYC_D_QLOSS, - SonyLFPCyclicDegradation, -) -from simses.model.degradation.sony_lfp_calendar import ( - T_REF, - C_RINC as CAL_C_RINC, - D_RINC as CAL_D_RINC, - C_QLOSS as CAL_C_QLOSS, - D_QLOSS as CAL_D_QLOSS, - EA_RINC, - EA_QLOSS, - K_REF_RINC, - K_REF_QLOSS, - R, - SonyLFPCalendarDegradation, -) - - -# --------------------------------------------------------------------------- -# Large-format 280 Ah LFP cell -# --------------------------------------------------------------------------- - +import matplotlib -class LFP280Ah(SonyLFP): - """280 Ah / 3.2 V prismatic LFP cell scaled from the SonyLFP OCV/resistance curves.""" - _SCALE = 3.0 / 280.0 # resistance scales inversely with capacity +matplotlib.use("Agg") +import matplotlib.pyplot as plt - def __init__(self): - super().__init__() - self.electrical.nominal_capacity = 280.0 # Ah - # Thermal properties for a large-format prismatic cell (vs. the 70 g 26650 reference). - # mass=1.5 kg, h=23 W/m²K gives C_th≈6.5 MJ/K, R_th≈1.73 mK/W, τ≈3.1 h - # → ΔT ≈ 10 °C at end of 2-hour C/2 discharge. - self.thermal.mass = 3.0 # kg per cell - self.thermal.convection_coefficient = 23.0 # W/m²K +from h2integrate import H2IntegrateModel +from h2integrate.storage.battery.battery_performance import BatteryPerformanceModel - def internal_resistance(self, state): - return super().internal_resistance(state) * self._SCALE +HERE = Path(__file__).parent # --------------------------------------------------------------------------- -# Scaled degradation models targeting ~1-2 % capacity fade over 365 days. -# The Sony LFP parametrization is for a 3 Ah cylindrical cell; large-format -# prismatic LFP cells degrade roughly 5x slower per calendar/cycle unit. -# A uniform scale factor of 0.2 brings the total fade into the 1-2 % range. +# Build and set up the H2Integrate model (standard interface) # --------------------------------------------------------------------------- -_DEG_SCALE = 0.4 - - -class ScaledLFPCalendarDegradation(SonyLFPCalendarDegradation): - def update_capacity(self, state: BatteryState, dt: float, accumulated_qloss: float) -> float: - if dt == 0.0: - return 0.0 - T_K = state.T + 273.15 - T_REF_K = T_REF + 273.15 - k_T_q = (K_REF_QLOSS * _DEG_SCALE) * math.exp(-EA_QLOSS / R * (1.0 / T_K - 1.0 / T_REF_K)) - k_soc_q = CAL_C_QLOSS * (state.soc - 0.5) ** 3 + CAL_D_QLOSS - stress_q = k_T_q * k_soc_q - if stress_q > 0.0: - virtual_time = (accumulated_qloss / stress_q) ** 2 - delta_q = stress_q * math.sqrt(virtual_time + dt) - accumulated_qloss - else: - delta_q = 0.0 - return delta_q - - def update_resistance(self, state: BatteryState, dt: float) -> float: - if dt == 0.0: - return 0.0 - T_K = state.T + 273.15 - T_REF_K = T_REF + 273.15 - k_T_r = (K_REF_RINC * _DEG_SCALE) * math.exp(-EA_RINC / R * (1.0 / T_K - 1.0 / T_REF_K)) - k_soc_r = CAL_C_RINC * (state.soc - 0.5) ** 2 + CAL_D_RINC - return k_T_r * k_soc_r * dt - - -class ScaledLFPCyclicDegradation(SonyLFPCyclicDegradation): - def update_capacity( - self, state: BatteryState, half_cycle: HalfCycle, accumulated_qloss: float - ) -> float: - delta_fec = half_cycle.full_equivalent_cycles - if delta_fec == 0.0: - return 0.0 - k_crate_q = (A_QLOSS * _DEG_SCALE) * half_cycle.c_rate + (B_QLOSS * _DEG_SCALE) - k_dod_q = CYC_C_QLOSS * (half_cycle.depth_of_discharge - 0.6) ** 3 + CYC_D_QLOSS - stress_q = k_crate_q * k_dod_q - if stress_q > 0.0: - virtual_fec = (accumulated_qloss * 100.0 / stress_q) ** 2 - delta_q = stress_q * math.sqrt(virtual_fec + delta_fec) / 100.0 - accumulated_qloss - else: - delta_q = 0.0 - return delta_q - - def update_resistance(self, state: BatteryState, half_cycle: HalfCycle) -> float: - delta_fec = half_cycle.full_equivalent_cycles - if delta_fec == 0.0: - return 0.0 - k_crate_r = (A_RINC * _DEG_SCALE) * half_cycle.c_rate + (B_RINC * _DEG_SCALE) - k_dod_r = CYC_C_RINC * (half_cycle.depth_of_discharge - 0.5) ** 3 + CYC_D_RINC - return k_crate_r * k_dod_r * delta_fec / 100.0 - +model = H2IntegrateModel(HERE / "98_battery_degradation.yaml") +model.setup() + +# Simulation horizon (read from the plant config) +sim = model.plant_config["plant"]["simulation"] +dt = float(sim["dt"]) +n_steps = int(sim["n_timesteps"]) +steps_per_day = round(86400 / dt) + +# Fixed Megapack-style pack (read from the tech config) +shared = model.technology_config["technologies"]["battery"]["model_inputs"]["shared_parameters"] +series_count = shared["series_count"] +parallel_count = shared["parallel_count"] +cell_nominal_voltage = 3.2 # V +cell_capacity_ah = 280.0 # Ah +usable_fraction = shared["max_soc_fraction"] - shared["min_soc_fraction"] + +nominal_energy_kwh = series_count * parallel_count * cell_nominal_voltage * cell_capacity_ah / 1e3 +usable_energy_kwh = nominal_energy_kwh * usable_fraction # --------------------------------------------------------------------------- -# Battery pack: 239s x 18p -> 764.8 V * 5040 Ah ~ 3855 kWh +# Daily set-point profile (H2I convention: + = discharge, - = charge, kW) +# 0 - 2 h discharge at C/2 of usable energy +# 2 - 12 h charge at C/10 of usable energy +# 12 - 24 h rest +# Defined in continuous time so it is independent of the timestep. # --------------------------------------------------------------------------- -cell = LFP280Ah() - -battery = Battery( - cell=cell, - circuit=(239, 18), - initial_states={"start_soc": 1.0, "start_T": 25.0}, - degradation=DegradationModel( - calendar=ScaledLFPCalendarDegradation(), - cyclic=ScaledLFPCyclicDegradation(), - initial_soc=1.0, - initial_state=DegradationState(qloss_cal=1e-4), - ), +p_discharge = 0.5 * usable_energy_kwh # kW, + = discharge +p_charge = -0.1 * usable_energy_kwh # kW, - = charge + +t_mid = (np.arange(n_steps) + 0.5) * dt # midpoint time [s] from sim start +t_day = t_mid % 86400 # position within the current day +set_point = np.where( + t_day < 2 * 3600, + p_discharge, + np.where(t_day < 12 * 3600, p_charge, 0.0), ) -summary = pd.Series( - { - "nominal_capacity [Ah]": battery.nominal_capacity, - "nominal_voltage [V]": battery.nominal_voltage, - "nominal_energy [kWh]": battery.nominal_energy_capacity / 1e3, - "max_charge_current [A]": battery.max_charge_current, - "max_discharge_current [A]": battery.max_discharge_current, - "thermal_capacity [kJ/K]": battery.thermal_capacity / 1e3, - } -) -print("Battery summary:") -print(summary.to_string()) -print() +n_years = n_steps / (365 * steps_per_day) +print(f"Simulation: {n_years:.2f} years, {n_steps:,} steps, dt={int(dt)} s") +print(f"Nominal energy: {nominal_energy_kwh:.0f} kWh, usable: {usable_energy_kwh:.0f} kWh") +print(f"Discharge: {p_discharge / 1e3:.3f} MW (C/2) Charge: {p_charge / 1e3:.3f} MW (C/10)") # --------------------------------------------------------------------------- -# 15-year power demand profile (dt = 15 min = 900 s) -# -# Each day (96 steps): -# 8 steps discharge at C/2 (~-1.927 MW) -# 40 steps charge at C/10 (~+385 kW) -# 48 steps rest 0 W +# Drive the battery via its set-point (passthrough controller -> command value) +# and run the model. # --------------------------------------------------------------------------- -dt = 900.0 -STEPS_PER_DAY = 96 -N_YEARS = 15 -N_DAYS = N_YEARS * 365 -n_steps = N_DAYS * STEPS_PER_DAY - -P_nom = battery.nominal_energy_capacity # Wh -P_disch = -0.5 * P_nom # C/2 discharge -P_chg = +0.1 * P_nom # C/10 charge +model.prob.set_val("battery.electricity_set_point", set_point, units="kW") +model.run() -day_profile = np.concatenate( - [ - np.full(8, P_disch), - np.full(40, P_chg), - np.full(48, 0.0), - ] -) -power_profile = np.tile(day_profile, N_DAYS) - -print(f"Simulation: {N_YEARS} years ({N_DAYS} days), {n_steps:,} steps, dt={int(dt)} s") -print(f"Discharge: {P_disch/1e6:.3f} MW (C/2)") -print(f"Charge: {P_chg/1e3:.1f} kW (C/10)") -print() - -# --------------------------------------------------------------------------- -# Thermal model: constant 25 °C ambient, battery registered as thermal node # --------------------------------------------------------------------------- -thermal = AmbientThermalModel(T_ambient=25.0, components=[battery]) - -# --------------------------------------------------------------------------- -# Simulation loop +# Retrieve the battery internal timeseries and assemble a DataFrame # --------------------------------------------------------------------------- -keys = ["soc", "v", "i", "T", "loss", "heat", "soh_Q", "soh_R", "power"] -log = {k: np.empty(n_steps) for k in keys} - -for i, p in enumerate(tqdm(power_profile, desc=f"Simulating {N_YEARS} years")): - battery.step(float(p), dt) - thermal.step(dt) # update battery temperature after each step - for k in keys: - log[k][i] = getattr(battery.state, k) +battery = next(model.prob.model.system_iter(typ=BatteryPerformanceModel)) +res = battery.results index = pd.date_range("2026-01-01", periods=n_steps, freq=f"{int(dt)}s") -df = pd.DataFrame(log, index=index) - -print("\nFirst rows:") -print(df.head().to_string()) -print( - f"\nFinal SOH_Q : {df['soh_Q'].iloc[-1]:.4f} \ - ({(1 - df['soh_Q'].iloc[-1]) * 100:.2f} % capacity fade)" +df = pd.DataFrame( + { + "soc": res["soc"], + "v": res["voltage"], + "T": res["temperature"], + "loss": res["battery_loss"], + "heat": res["battery_heat"], + "soh_Q": res["soh_capacity"], + "soh_R": res["soh_resistance"], + "power_ac": res["power_ac"], + "power_dc": res["power_dc"], + "conv_loss": res["converter_loss"], + }, + index=index, ) -print(f"Final SOH_R : {df['soh_R'].iloc[-1]:.4f}") + +fade_pct = (1 - df["soh_Q"].iloc[-1]) * 100 +growth_pct = (df["soh_R"].iloc[-1] - 1) * 100 +print(f"Final capacity fade: {fade_pct:.2f} % resistance growth: {growth_pct:.2f} %") + +week = 7 * steps_per_day # --------------------------------------------------------------------------- -# Plots +# Plot 1: AC & DC power (first 7 days) # --------------------------------------------------------------------------- - -# Plot 1: Power demand (first 7 days) fig, ax = plt.subplots(figsize=(12, 3)) -df["power"].iloc[: 7 * STEPS_PER_DAY].plot(ax=ax) -ax.set_title("Power demand profile -- first 7 days") -ax.set_ylabel("Power [W]") +df[["power_ac", "power_dc"]].iloc[:week].plot(ax=ax) +ax.set_title("AC and DC power -- first 7 days") +ax.set_ylabel("Power [MW]") ax.axhline(0, color="gray", linewidth=0.5) -ax.yaxis.set_major_formatter(plt.FuncFormatter(lambda x, _: f"{x/1e6:.2f} MW")) +ax.yaxis.set_major_formatter(plt.FuncFormatter(lambda x, _: f"{x / 1e6:.2f}")) fig.tight_layout() fig.savefig(HERE / "plot_power.png", dpi=150) print("Saved plot_power.png") +# --------------------------------------------------------------------------- # Plot 2: SOC (hourly mean) +# --------------------------------------------------------------------------- fig, ax = plt.subplots(figsize=(12, 3)) -df["soc"].resample("1h").mean().plot(ax=ax, title="State of Charge -- hourly mean") +df["soc"].resample("1h").mean().plot(ax=ax, title="State of charge -- hourly mean") ax.set_ylabel("SOC [p.u.]") fig.tight_layout() fig.savefig(HERE / "plot_soc.png", dpi=150) print("Saved plot_soc.png") -# Plot 2b: Battery temperature -- first 7 days (15-min resolution) +# --------------------------------------------------------------------------- +# Plot 3: Battery temperature (first 7 days) +# --------------------------------------------------------------------------- fig, ax = plt.subplots(figsize=(12, 3)) -df["T"].iloc[: 7 * STEPS_PER_DAY].plot(ax=ax, title="Battery temperature -- first 7 days") -ax.set_ylabel("Temperature [°C]") -ax.axhline(25.0, color="gray", linewidth=0.8, linestyle="--", label="T_ambient = 25 °C") +df["T"].iloc[:week].plot(ax=ax, title="Battery temperature -- first 7 days") +ax.set_ylabel("Temperature [degC]") +ax.axhline(25.0, color="gray", linewidth=0.8, linestyle="--", label="T_ambient = 25 degC") ax.legend() fig.tight_layout() fig.savefig(HERE / "plot_temperature_7days.png", dpi=150) print("Saved plot_temperature_7days.png") -# Plot 3a: Terminal voltage -- first 7 days (raw 15-min data) +# --------------------------------------------------------------------------- +# Plot 4a: Terminal voltage (first 7 days) +# --------------------------------------------------------------------------- fig, ax = plt.subplots(figsize=(12, 3)) -df["v"].iloc[: 7 * STEPS_PER_DAY].plot(ax=ax, title="Terminal voltage -- first 7 days [V]") +df["v"].iloc[:week].plot(ax=ax, title="Terminal voltage -- first 7 days [V]") ax.set_ylabel("Voltage [V]") fig.tight_layout() fig.savefig(HERE / "plot_voltage_7days.png", dpi=150) print("Saved plot_voltage_7days.png") -# Plot 3b: Terminal voltage (daily mean, full year) +# --------------------------------------------------------------------------- +# Plot 4b: Terminal voltage (daily mean) +# --------------------------------------------------------------------------- fig, ax = plt.subplots(figsize=(12, 3)) df["v"].resample("1D").mean().plot(ax=ax, title="Terminal voltage -- daily mean [V]") ax.set_ylabel("Voltage [V]") @@ -271,16 +175,20 @@ def update_resistance(self, state: BatteryState, half_cycle: HalfCycle) -> float fig.savefig(HERE / "plot_voltage.png", dpi=150) print("Saved plot_voltage.png") -# Plot 4: Degradation over 15 years +# --------------------------------------------------------------------------- +# Plot 5: Degradation over the simulation horizon +# --------------------------------------------------------------------------- df_monthly = df[["soh_Q", "soh_R"]].resample("ME").last() capacity_fade_pct = (1 - df_monthly["soh_Q"]) * 100 resistance_growth_pct = (df_monthly["soh_R"] - 1) * 100 fig, ax1 = plt.subplots(figsize=(12, 4)) ax2 = ax1.twinx() - ax1.plot( - capacity_fade_pct.index, capacity_fade_pct.values, color="steelblue", label="Capacity fade [%]" + capacity_fade_pct.index, + capacity_fade_pct.values, + color="steelblue", + label="Capacity fade [%]", ) ax2.plot( resistance_growth_pct.index, @@ -288,33 +196,30 @@ def update_resistance(self, state: BatteryState, half_cycle: HalfCycle) -> float color="darkorange", label="Resistance growth [%]", ) - ax1.set_ylabel("Capacity fade [%]", color="steelblue") ax2.set_ylabel("Resistance growth [%]", color="darkorange") ax1.tick_params(axis="y", labelcolor="steelblue") ax2.tick_params(axis="y", labelcolor="darkorange") +ax1.axhline(20, color="steelblue", linewidth=0.8, linestyle="--", alpha=0.6) ax1.set_title( - f"Battery degradation over {N_YEARS} years " - f"(final: {capacity_fade_pct.iloc[-1]:.1f} % capacity fade, " + f"Battery degradation over {n_years:.1f} years " + f"(final: {capacity_fade_pct.iloc[-1]:.1f} % fade, " f"{resistance_growth_pct.iloc[-1]:.1f} % resistance growth)" ) ax1.set_xlabel("Date") - lines1, labels1 = ax1.get_legend_handles_labels() lines2, labels2 = ax2.get_legend_handles_labels() ax1.legend(lines1 + lines2, labels1 + labels2, loc="upper left") - -# Mark end-of-life reference (20 % capacity fade) -ax1.axhline(20, color="steelblue", linewidth=0.8, linestyle="--", alpha=0.6, label="EoL (20 %)") - fig.tight_layout() -fig.savefig(HERE / "plot_degradation_15yr.png", dpi=150) -print("Saved plot_degradation_15yr.png") +fig.savefig(HERE / "plot_degradation.png", dpi=150) +print("Saved plot_degradation.png") -# Plot 5: Losses and heat (daily mean) +# --------------------------------------------------------------------------- +# Plot 6: Losses and heat (daily mean) +# --------------------------------------------------------------------------- fig, ax = plt.subplots(figsize=(12, 3)) -df[["loss", "heat"]].resample("1D").mean().plot( - ax=ax, title="Battery losses and heat -- daily mean [W]" +df[["loss", "heat", "conv_loss"]].resample("1D").mean().plot( + ax=ax, title="Battery losses, heat and converter losses -- daily mean [W]" ) ax.set_ylabel("Power [W]") fig.tight_layout() diff --git a/h2integrate/storage/battery/battery_performance.py b/h2integrate/storage/battery/battery_performance.py index 99ef3c6d0..f64384fd7 100644 --- a/h2integrate/storage/battery/battery_performance.py +++ b/h2integrate/storage/battery/battery_performance.py @@ -279,6 +279,8 @@ def setup(self): self.commodity_rate_units = self.config.commodity_rate_units self.commodity_amount_units = self.config.commodity_amount_units + super().setup() + self.add_discrete_input( "solar_resource_data", val={}, @@ -293,8 +295,6 @@ def setup(self): # TODO degradation: adjustments for degradation - super().setup() - def compute(self, inputs, outputs, discrete_inputs=[], discrete_outputs=[]): """Run the storage performance model.""" self.current_soc = self.config.init_soc_fraction @@ -371,6 +371,22 @@ def compute(self, inputs, outputs, discrete_inputs=[], discrete_outputs=[]): ############# + # Store the full SimSES timeseries for downstream diagnostics/plotting + # (e.g. example 98 degradation, temperature, voltage, and loss plots). + self.results = { + "soc": log["soc"], + "voltage": log["v"], + "current": log["i"], + "temperature": log["T"], + "battery_loss": log["loss"], + "battery_heat": log["heat"], + "soh_capacity": log["soh_Q"], + "soh_resistance": log["soh_R"], + "power_ac": power_ac, + "power_dc": power_dc, + "converter_loss": conv_loss, + } + # Populate all OpenMDAO outputs defined in this class and its parent classes. # Convert SimSES AC power (W, +charge) back to H2I convention # (commodity_rate_units, +discharge). From 9353e5fdac376599a1875793a52dd1d776a6ac98 Mon Sep 17 00:00:00 2001 From: Jared Thomas Date: Fri, 14 Aug 2026 11:52:40 -0600 Subject: [PATCH 14/18] recreate full script study by Ankit in H2I --- .../run_battery_degradation.py | 126 +++++++++++-- .../storage/battery/battery_performance.py | 178 +++++++++++++++--- 2 files changed, 258 insertions(+), 46 deletions(-) diff --git a/examples/98_battery_degradation/run_battery_degradation.py b/examples/98_battery_degradation/run_battery_degradation.py index 53e4a3017..906858ab6 100644 --- a/examples/98_battery_degradation/run_battery_degradation.py +++ b/examples/98_battery_degradation/run_battery_degradation.py @@ -1,14 +1,15 @@ """ Example 98: Battery degradation through the standard H2Integrate interface. -This example runs the SimSES-backed ``BatteryPerformanceModel`` through the standard +This example reproduces the SimSES battery-degradation reference study using the standard ``H2IntegrateModel`` interface (YAML-configured plant/tech/driver). The battery has no control strategy, so a passthrough controller forwards the ``electricity_set_point`` profile to the battery as its dispatch command (positive = discharge, negative = charge). A repeating daily charge/discharge cycle is applied for the full simulation horizon and the battery internal timeseries (voltage, temperature, state-of-health, losses) are read -from ``BatteryPerformanceModel.results`` to produce the reference plots: +from the battery component's OpenMDAO outputs (via ``model.prob.get_val``) to produce the +reference plots: 1. AC and DC power (first 7 days) 2. State of charge (hourly mean) @@ -17,9 +18,11 @@ 5. Capacity fade and resistance growth over the horizon 6. Battery/converter losses and heat (daily mean) -H2Integrate constrains a single simulation to exactly one year -(``n_timesteps * dt == 31_536_000 s``), so this example runs a 1-year degradation -study. The battery pack topology and degradation scaling live in ``tech_config.yaml``. +Matching the partner reference, this runs a 15-year degradation study at 15-minute +resolution (``n_timesteps = 525600``, ``dt = 900`` s in ``plant_config.yaml``). The +horizon is fully configurable: ``n_timesteps * dt`` may cover any positive duration up to +the plant life. The battery pack topology and degradation scaling live in +``tech_config.yaml``. """ from pathlib import Path @@ -33,7 +36,6 @@ import matplotlib.pyplot as plt from h2integrate import H2IntegrateModel -from h2integrate.storage.battery.battery_performance import BatteryPerformanceModel HERE = Path(__file__).parent @@ -92,24 +94,25 @@ model.run() # --------------------------------------------------------------------------- -# Retrieve the battery internal timeseries and assemble a DataFrame +# Retrieve the battery internal timeseries and assemble a DataFrame. +# Each SimSES timeseries is exposed as an OpenMDAO output of the battery +# component and is read through the standard H2Integrate ``get_val`` interface +# (outputs are promoted to the "battery" tech group), so the model class does +# not need to be imported into this script. # --------------------------------------------------------------------------- -battery = next(model.prob.model.system_iter(typ=BatteryPerformanceModel)) -res = battery.results - index = pd.date_range("2026-01-01", periods=n_steps, freq=f"{int(dt)}s") df = pd.DataFrame( { - "soc": res["soc"], - "v": res["voltage"], - "T": res["temperature"], - "loss": res["battery_loss"], - "heat": res["battery_heat"], - "soh_Q": res["soh_capacity"], - "soh_R": res["soh_resistance"], - "power_ac": res["power_ac"], - "power_dc": res["power_dc"], - "conv_loss": res["converter_loss"], + "soc": model.prob.get_val("battery.SOC", units="unitless"), + "v": model.prob.get_val("battery.voltage", units="V"), + "T": model.prob.get_val("battery.temperature", units="degC"), + "loss": model.prob.get_val("battery.battery_loss", units="W"), + "heat": model.prob.get_val("battery.battery_heat", units="W"), + "soh_Q": model.prob.get_val("battery.soh_capacity"), + "soh_R": model.prob.get_val("battery.soh_resistance"), + "power_ac": model.prob.get_val("battery.power_ac", units="W"), + "power_dc": model.prob.get_val("battery.power_dc", units="W"), + "conv_loss": model.prob.get_val("battery.converter_loss", units="W"), }, index=index, ) @@ -225,3 +228,86 @@ fig.tight_layout() fig.savefig(HERE / "plot_losses.png", dpi=150) print("Saved plot_losses.png") + +# --------------------------------------------------------------------------- +# Finance results (ProFAST LCOE). +# The battery's per-year, degradation-aware capacity factor feeds ProFAST's +# multi-year cash-flow analysis over the full 30-year plant life. When the +# simulated capacity SOH reaches the end-of-life threshold (``eol_soh_capacity`` +# in tech_config.yaml), the battery is replaced and the capacity-factor cycle +# repeats from a fresh unit; each replacement is charged as new CapEx in ProFAST. +# --------------------------------------------------------------------------- +lcoe = model.prob.get_val("finance_subgroup_battery.LCOE", units="USD/(MW*h)")[0] +capacity_factor_by_year = model.prob.get_val("battery.capacity_factor") +replacement_schedule = model.prob.get_val("battery.replacement_schedule") +replacement_years = [int(y) for y in np.flatnonzero(replacement_schedule) + 1] +print(f"LCOE (levelized cost of storage): {lcoe:.2f} $/MWh") +print( + "Discharge capacity factor - year 1: " + f"{capacity_factor_by_year[0] * 100:.2f} %, " + f"year 15: {capacity_factor_by_year[14] * 100:.2f} %, " + f"year 30: {capacity_factor_by_year[29] * 100:.2f} %" +) +if replacement_years: + print(f"Battery reaches EOL SOH and is replaced in plant year(s): {replacement_years}") +else: + print("Battery does not reach EOL SOH within the plant life (no replacement scheduled).") + +# --------------------------------------------------------------------------- +# Plot 7: Per-year capacity factor used by the finance model (ProFAST) +# This is the exact per-year capacity factor fed to ProFAST: actual values over +# the simulated years, scaled with the projected SOH beyond the simulation, and +# reset to a fresh battery at each end-of-life replacement. +# --------------------------------------------------------------------------- +plant_years = np.arange(1, len(capacity_factor_by_year) + 1) +fig, ax = plt.subplots(figsize=(12, 3)) +ax.step(plant_years, capacity_factor_by_year * 100, where="mid", color="teal") +for i, yr in enumerate(replacement_years): + ax.axvline( + yr, + color="firebrick", + linestyle="--", + linewidth=0.8, + label="Battery replacement" if i == 0 else None, + ) +ax.set_title("Per-year discharge capacity factor used by ProFAST") +ax.set_xlabel("Plant year") +ax.set_ylabel("Capacity factor [%]") +ax.set_xlim(1, len(capacity_factor_by_year)) +if replacement_years: + ax.legend(loc="lower left") +fig.tight_layout() +fig.savefig(HERE / "plot_capacity_factor.png", dpi=150) +print("Saved plot_capacity_factor.png") + +# --------------------------------------------------------------------------- +# Plot 8: State of health over the full simulation +# --------------------------------------------------------------------------- +perf_params = model.technology_config["technologies"]["battery"]["model_inputs"][ + "performance_parameters" +] +eol_soh = float(perf_params.get("eol_soh_capacity", 0.8)) + +fig, ax1 = plt.subplots(figsize=(12, 3)) +ax2 = ax1.twinx() +ax1.plot(df.index, df["soh_Q"] * 100, color="steelblue", label="Capacity SOH [%]") +ax2.plot(df.index, df["soh_R"], color="darkorange", label="Resistance SOH [x nominal]") +ax1.axhline( + eol_soh * 100, + color="firebrick", + linestyle="--", + linewidth=0.8, + label=f"EOL SOH = {eol_soh * 100:.0f} %", +) +ax1.set_ylabel("Capacity SOH [%]", color="steelblue") +ax2.set_ylabel("Resistance SOH [x nominal]", color="darkorange") +ax1.tick_params(axis="y", labelcolor="steelblue") +ax2.tick_params(axis="y", labelcolor="darkorange") +ax1.set_xlabel("Date") +ax1.set_title("Battery state of health over the full simulation") +lines1, labels1 = ax1.get_legend_handles_labels() +lines2, labels2 = ax2.get_legend_handles_labels() +ax1.legend(lines1 + lines2, labels1 + labels2, loc="lower left") +fig.tight_layout() +fig.savefig(HERE / "plot_soh.png", dpi=150) +print("Saved plot_soh.png") diff --git a/h2integrate/storage/battery/battery_performance.py b/h2integrate/storage/battery/battery_performance.py index f64384fd7..b8534c4c3 100644 --- a/h2integrate/storage/battery/battery_performance.py +++ b/h2integrate/storage/battery/battery_performance.py @@ -1,6 +1,15 @@ +""" +# NOTE: ``simses.battery`` must be imported before ``simses.degradation`` to avoid a +# circular import within simses (>=2.1.1): importing ``simses.degradation`` first leaves +# ``simses.degradation.calendar`` partially initialized when ``simses.battery.cell`` pulls +# in ``simses.degradation.degradation``. Importing the battery package first fully loads +# both sub-packages in a safe order. (Plain ``import`` sorts above ``from`` imports.) +""" + import math import numpy as np +import simses.battery # noqa: F401 (import-order side effect; see note above) from attrs import field, define from openmdao.utils import units as om_units from simses.degradation import DegradationModel @@ -196,6 +205,7 @@ class BatteryPerformanceModelConfig(StoragePerformanceBaseConfig): round_trip_efficiency: float | None = field(default=None, validator=range_val_or_none(0, 1)) _DEG_SCALE: float = field(default=0.7056, validator=range_val(0, 1)) + eol_soh_capacity: float = field(default=0.8, validator=range_val(0, 1)) # TODO convert from power and energy ratings (see math in chat) series_count: int = field(default=336, converter=int, validator=gt_zero) parallel_count: int = field(default=16, converter=int, validator=gt_zero) @@ -293,6 +303,47 @@ def setup(self): desc="Electricity demand for running battery auxiliary systems", ) + # Internal SimSES timeseries exposed as OpenMDAO outputs (one per quantity) for + # downstream diagnostics/plotting. + self.add_output( + "voltage", shape=self.n_timesteps, units="V", desc="Battery terminal voltage" + ) + self.add_output("current", shape=self.n_timesteps, units="A", desc="Battery current") + self.add_output( + "temperature", shape=self.n_timesteps, units="degC", desc="Battery temperature" + ) + self.add_output( + "battery_loss", shape=self.n_timesteps, units="W", desc="Battery internal loss" + ) + self.add_output( + "battery_heat", shape=self.n_timesteps, units="W", desc="Battery heat generation" + ) + self.add_output( + "soh_capacity", + shape=self.n_timesteps, + units="unitless", + desc="State of health, capacity (fraction of nominal capacity)", + ) + self.add_output( + "soh_resistance", + shape=self.n_timesteps, + units="unitless", + desc="State of health, resistance (multiple of nominal resistance)", + ) + self.add_output( + "power_ac", + shape=self.n_timesteps, + units="W", + desc="AC-side power (positive = charge)", + ) + self.add_output( + "power_dc", + shape=self.n_timesteps, + units="W", + desc="DC-side power (positive = charge)", + ) + self.add_output("converter_loss", shape=self.n_timesteps, units="W", desc="Converter loss") + # TODO degradation: adjustments for degradation def compute(self, inputs, outputs, discrete_inputs=[], discrete_outputs=[]): @@ -309,8 +360,6 @@ def compute(self, inputs, outputs, discrete_inputs=[], discrete_outputs=[]): # H2I dispatch command: positive = discharge, negative = charge (commodity_rate_units) power_profile = inputs[f"{self.commodity}_command_value"] - ### from Ankit - # --------------------------------------------------------------------------- # Battery pack + inverter (fixed Megapack-style topology, from config) # --------------------------------------------------------------------------- @@ -371,37 +420,30 @@ def compute(self, inputs, outputs, discrete_inputs=[], discrete_outputs=[]): ############# - # Store the full SimSES timeseries for downstream diagnostics/plotting - # (e.g. example 98 degradation, temperature, voltage, and loss plots). - self.results = { - "soc": log["soc"], - "voltage": log["v"], - "current": log["i"], - "temperature": log["T"], - "battery_loss": log["loss"], - "battery_heat": log["heat"], - "soh_capacity": log["soh_Q"], - "soh_resistance": log["soh_R"], - "power_ac": power_ac, - "power_dc": power_dc, - "converter_loss": conv_loss, - } - # Populate all OpenMDAO outputs defined in this class and its parent classes. # Convert SimSES AC power (W, +charge) back to H2I convention # (commodity_rate_units, +discharge). - soc_ts = log["soc"] power_ts = -om_units.convert_units(power_ac, "W", self.commodity_rate_units) # --- BatteryPerformanceModel outputs --- # TODO calc aux power outputs[f"{self.commodity}_auxiliary_demand"] = np.zeros(self.n_timesteps) + outputs["voltage"] = log["v"] + outputs["current"] = log["i"] + outputs["temperature"] = log["T"] + outputs["battery_loss"] = log["loss"] + outputs["battery_heat"] = log["heat"] + outputs["soh_capacity"] = log["soh_Q"] + outputs["soh_resistance"] = log["soh_R"] + outputs["power_ac"] = power_ac + outputs["power_dc"] = power_dc + outputs["converter_loss"] = conv_loss # --- StoragePerformanceBase outputs --- outputs["storage_duration"] = ( storage_capacity / discharge_rate if discharge_rate > 0 else 0.0 ) - outputs["SOC"] = soc_ts * 100.0 # fraction -> percent + outputs["SOC"] = log["soc"] * 100.0 # fraction -> percent outputs[f"storage_{self.commodity}_charge"] = np.where(power_ts < 0, power_ts, 0.0) outputs[f"storage_{self.commodity}_discharge"] = np.where(power_ts > 0, power_ts, 0.0) @@ -419,12 +461,96 @@ def compute(self, inputs, outputs, discrete_inputs=[], discrete_outputs=[]): outputs["capacity_factor"] = 0.0 outputs["standard_capacity_factor"] = 0.0 else: - outputs["capacity_factor"] = outputs[f"total_{self.commodity}_produced"] / ( - discharge_rate * self.n_timesteps * self.dt_amount - ) - total_commodity_discharged = ( - outputs[f"storage_{self.commodity}_discharge"].sum() * self.dt_amount - ) + # Gross discharge timeseries (commodity_rate_units, discharge only). + discharge_ts = outputs[f"storage_{self.commodity}_discharge"] + total_commodity_discharged = discharge_ts.sum() * self.dt_amount + + # Scalar average discharge capacity factor over the whole simulation. outputs["standard_capacity_factor"] = total_commodity_discharged / ( discharge_rate * self.n_timesteps * self.dt_amount ) + + # Per-year discharge capacity factor and year-end capacity state-of-health over + # the simulated horizon. The simulation may span whole years plus an optional + # partial trailing year; each simulated year gets its own capacity factor and + # end-of-year SOH. + steps_per_year = round(31_536_000 / self.dt) # timesteps in one year + n_sim_years = math.ceil(self.n_timesteps / steps_per_year) + soh_capacity_ts = log["soh_Q"] + sim_cf = np.zeros(n_sim_years) + sim_soh_year_end = np.zeros(n_sim_years) + for year in range(n_sim_years): + start = year * steps_per_year + end = min(start + steps_per_year, self.n_timesteps) + segment_hours = (end - start) * self.dt_amount + sim_cf[year] = (discharge_ts[start:end].sum() * self.dt_amount) / ( + discharge_rate * segment_hours + ) + sim_soh_year_end[year] = soh_capacity_ts[end - 1] + + # Annual capacity-SOH degradation rate used to project SOH beyond the simulated + # horizon (i.e. once the simulated years are exhausted before the battery hits + # end-of-life): + # - Less than one year simulated: extrapolate the average degradation over the + # whole simulation to a per-year rate. + # - One year or more simulated: use the degradation over the last full + # simulated year. + years_simulated = self.n_timesteps / steps_per_year + soh_start = soh_capacity_ts[0] + if years_simulated < 1.0: + annual_deg_rate = (soh_start - sim_soh_year_end[-1]) / years_simulated + else: + n_full_years = int(self.n_timesteps // steps_per_year) + idx_after = n_full_years * steps_per_year - 1 + idx_before = (n_full_years - 1) * steps_per_year - 1 + soh_before = soh_capacity_ts[idx_before] if idx_before >= 0 else soh_start + annual_deg_rate = soh_before - soh_capacity_ts[idx_after] + annual_deg_rate = max(float(annual_deg_rate), 0.0) + + # Build one battery-life cycle of per-year year-end SOH and capacity factor, + # long enough to cover the whole plant life. Within the simulated years the + # actual per-year values are used. Beyond the simulated horizon the SOH keeps + # degrading at annual_deg_rate and the capacity factor is scaled down in + # proportion to the declining SOH (relative to the last simulated year), so the + # capacity factor tracks degradation rather than being held constant. + if years_simulated < 1.0: + # A sub-year simulation never completes a full year, so project every year + # from the start-of-life SOH at the extrapolated annual rate. + cycle_soh_end = soh_start - annual_deg_rate * (np.arange(self.plant_life) + 1) + else: + cycle_soh_end = np.array( + [ + sim_soh_year_end[y] + if y < n_sim_years + else sim_soh_year_end[-1] - annual_deg_rate * (y - (n_sim_years - 1)) + for y in range(self.plant_life) + ] + ) + + soh_ref = sim_soh_year_end[-1] + cycle_cf = np.array( + [ + sim_cf[y] + if y < n_sim_years + else sim_cf[-1] * max(cycle_soh_end[y], 0.0) / soh_ref + for y in range(self.plant_life) + ] + ) + + # Walk the plant life. When the projected year-end SOH reaches the user-specified + # end-of-life threshold, the battery is replaced (fresh unit) at the start of the + # following year and the degradation / capacity-factor cycle restarts. + eol_soh = self.config.eol_soh_capacity + cf_per_year = np.zeros(self.plant_life) + replacement_schedule = np.zeros(self.plant_life) + cycle_year = 0 + for plant_year in range(self.plant_life): + cf_per_year[plant_year] = cycle_cf[cycle_year] + if cycle_soh_end[cycle_year] <= eol_soh: + if plant_year + 1 < self.plant_life: + replacement_schedule[plant_year + 1] = 1.0 + cycle_year = 0 + else: + cycle_year += 1 + outputs["capacity_factor"] = cf_per_year + outputs["replacement_schedule"] = replacement_schedule From 9356b9cf777f73f0732d29b876637ac4b6ed0b59 Mon Sep 17 00:00:00 2001 From: Jared Thomas Date: Fri, 14 Aug 2026 11:54:11 -0600 Subject: [PATCH 15/18] framework changes to allow for non-annual simulation --- h2integrate/core/inputs/validation.py | 32 +++++++++++++++---- h2integrate/core/test/test_framework.py | 30 +++++++++++------ .../storage/battery/atb_battery_cost.py | 12 ++++--- 3 files changed, 53 insertions(+), 21 deletions(-) diff --git a/h2integrate/core/inputs/validation.py b/h2integrate/core/inputs/validation.py index 5f0b7afa4..78091bc69 100644 --- a/h2integrate/core/inputs/validation.py +++ b/h2integrate/core/inputs/validation.py @@ -152,15 +152,33 @@ def load_tech_yaml(finput): def load_plant_yaml(finput): plant_config = _validate(finput, fschema_plant) - n_timesteps = plant_config["plant"]["simulation"]["n_timesteps"] - dt = plant_config["plant"]["simulation"]["dt"] - - if int(n_timesteps) * int(dt) != 31536000: # seconds in simulation must be seconds/year + n_timesteps = int(plant_config["plant"]["simulation"]["n_timesteps"]) + dt = int(plant_config["plant"]["simulation"]["dt"]) + plant_life = int(plant_config["plant"]["plant_life"]) + + seconds_per_year = 31_536_000 # 8760 h/year * 3600 s/h + seconds_simulated = n_timesteps * dt + years_simulated = seconds_simulated / seconds_per_year + + # The simulation may cover any positive duration (a fraction of a year, a single + # year, or multiple years). Performance/cost/finance models annualize their results + # using ``fraction_of_year_simulated`` (see the model base classes), so an arbitrary + # horizon is supported as long as it is positive and does not exceed the plant life. + if seconds_simulated <= 0: msg = ( - "H2Integrate does not currently support simulations that are less than or " - "greater than 1-year. Please ensure that " + "The simulation horizon must be positive. Please ensure that " "plant_config['plant']['simulation']['n_timesteps'] times " - "plant_config['plant']['simulation']['dt'] equals 31536000 (s)." + "plant_config['plant']['simulation']['dt'] is greater than 0 (s)." + ) + raise ValueError(msg) + + if years_simulated > plant_life: + msg = ( + "H2Integrate does not support simulations that are longer than the plant " + f"life. The configured simulation covers {years_simulated:.4g} years " + f"(n_timesteps={n_timesteps} * dt={dt} s = {seconds_simulated} s), but " + f"plant_config['plant']['plant_life'] is {plant_life} years. Either shorten " + "the simulation horizon or increase plant_life." ) raise ValueError(msg) diff --git a/h2integrate/core/test/test_framework.py b/h2integrate/core/test/test_framework.py index ac463eaeb..4d606d868 100644 --- a/h2integrate/core/test/test_framework.py +++ b/h2integrate/core/test/test_framework.py @@ -362,20 +362,32 @@ def test_unsupported_simulation_parameters(temp_dir): plant_config_data_dt = load_plant_yaml(temp_plant_config_dt) # docs fencepost end: DO NOT REMOVE - # Modify the n_timesteps entry for the temp_plant_config_ntimesteps - plant_config_data_ntimesteps["plant"]["simulation"]["n_timesteps"] = 8759 - # Modify the dt entry for the temp_plant_config_dt - plant_config_data_dt["plant"]["simulation"]["dt"] = 3601 + plant_life = int(plant_config_data_ntimesteps["plant"]["plant_life"]) - # Save the modified plant_configs YAML back + # Sub-year and multi-year horizons are now supported as long as they are positive + # and do not exceed the plant life. A ~half-year horizon should load without error. + plant_config_data_ntimesteps["plant"]["simulation"]["n_timesteps"] = 4380 # 0.5 year at 1 h with temp_plant_config_ntimesteps.open("w") as f: yaml.safe_dump(plant_config_data_ntimesteps, f) + load_plant_yaml(temp_plant_config_ntimesteps) + + # A multi-year horizon (2 years of hourly data) should also load without error. + plant_config_data_dt["plant"]["simulation"]["n_timesteps"] = 2 * 8760 with temp_plant_config_dt.open("w") as f: yaml.safe_dump(plant_config_data_dt, f) - - # check that error is thrown when loading config with invalid number of timesteps - with pytest.raises(ValueError, match="greater than 1-year"): - load_plant_yaml(plant_config_data_ntimesteps) + load_plant_yaml(temp_plant_config_dt) + + # A horizon longer than the plant life is not supported and must raise. + over_life = deepcopy(plant_config_data_dt) + over_life["plant"]["simulation"]["n_timesteps"] = (plant_life + 1) * 8760 + with pytest.raises(ValueError, match="longer than the plant"): + load_plant_yaml(over_life) + + # A non-positive horizon is invalid and must raise. + non_positive = deepcopy(plant_config_data_dt) + non_positive["plant"]["simulation"]["dt"] = 0 + with pytest.raises(ValueError, match="must be positive"): + load_plant_yaml(non_positive) @pytest.mark.unit diff --git a/h2integrate/storage/battery/atb_battery_cost.py b/h2integrate/storage/battery/atb_battery_cost.py index f63e943fe..321c0516e 100644 --- a/h2integrate/storage/battery/atb_battery_cost.py +++ b/h2integrate/storage/battery/atb_battery_cost.py @@ -1,4 +1,5 @@ from attrs import field, define +from numpy import inf from openmdao.utils import units from h2integrate.core.utilities import merge_shared_inputs @@ -62,9 +63,11 @@ class ATBBatteryCostModel(CostModelBaseClass): """ _time_step_bounds = ( - 3600, - 3600, - ) # (min, max) time step lengths (in seconds) compatible with this model + 1e-6, + inf, + ) # (min, max) time step lengths (in seconds) compatible with this model. The ATB + # cost model is time-step independent (it only uses storage capacity and charge + # rate), so it accepts any sub-hourly-to-hourly time step. def setup(self): self.config = ATBBatteryCostConfig.from_dict( @@ -105,8 +108,7 @@ def compute(self, inputs, outputs, discrete_inputs, discrete_outputs): storage_duration_hrs = max_capacity_kWh / max_charge_rate_kW if max_charge_rate_kW < 0: msg = ( - f"max_charge_rate cannot be less than zero and has value of " - f"{max_charge_rate_kW} kW" + f"max_charge_rate cannot be less than zero and has value of {max_charge_rate_kW} kW" ) raise UserWarning(msg) # CapEx equation from Cell E29 From 2ba40b624f0f1213ee0f111ee67612db6f9af9eb Mon Sep 17 00:00:00 2001 From: Jared Thomas Date: Fri, 14 Aug 2026 11:58:56 -0600 Subject: [PATCH 16/18] h2i files for running the example --- .../98_battery_degradation.yaml | 5 ++ .../98_battery_degradation/driver_config.yaml | 5 ++ .../98_battery_degradation/plant_config.yaml | 50 +++++++++++++++++++ .../98_battery_degradation/tech_config.yaml | 44 ++++++++++++++++ 4 files changed, 104 insertions(+) create mode 100644 examples/98_battery_degradation/98_battery_degradation.yaml create mode 100644 examples/98_battery_degradation/driver_config.yaml create mode 100644 examples/98_battery_degradation/plant_config.yaml create mode 100644 examples/98_battery_degradation/tech_config.yaml diff --git a/examples/98_battery_degradation/98_battery_degradation.yaml b/examples/98_battery_degradation/98_battery_degradation.yaml new file mode 100644 index 000000000..2aff4f2bc --- /dev/null +++ b/examples/98_battery_degradation/98_battery_degradation.yaml @@ -0,0 +1,5 @@ +name: H2Integrate_config +system_summary: Battery degradation study driven by a fixed daily charge/discharge cycle +driver_config: driver_config.yaml +technology_config: tech_config.yaml +plant_config: plant_config.yaml diff --git a/examples/98_battery_degradation/driver_config.yaml b/examples/98_battery_degradation/driver_config.yaml new file mode 100644 index 000000000..4b6d68338 --- /dev/null +++ b/examples/98_battery_degradation/driver_config.yaml @@ -0,0 +1,5 @@ +name: driver_config +description: Driver configuration for the battery degradation example +general: + folder_output: outputs + create_om_reports: false diff --git a/examples/98_battery_degradation/plant_config.yaml b/examples/98_battery_degradation/plant_config.yaml new file mode 100644 index 000000000..4a2392219 --- /dev/null +++ b/examples/98_battery_degradation/plant_config.yaml @@ -0,0 +1,50 @@ +name: plant_config +description: > + Single-battery plant used to exercise the SimSES-backed BatteryPerformanceModel + over a long horizon for degradation analysis. The battery is driven directly by a + set-point profile (no generation or demand technologies are connected). +plant: + plant_life: 30 + simulation: + # 15-year degradation study at 15-min resolution (matches the partner reference + # script): 15 years * 365 days * 96 steps/day = 525600 steps. n_timesteps * dt may + # cover any horizon up to plant_life years. + n_timesteps: 525600 + dt: 900 +finance_parameters: + finance_groups: + profast_model: + commodity: electricity + finance_model: ProFastLCO + model_inputs: + params: + analysis_start_year: 2025 + installation_time: 12 + inflation_rate: 0.025 + discount_rate: 0.09 + debt_equity_ratio: 2.62 + property_tax_and_insurance: 0.015 + total_income_tax_rate: 0.257 + capital_gains_tax_rate: 0.15 + sales_tax_rate: 0.045 + debt_interest_rate: 0.06 + debt_type: Revolving debt + loan_period_if_used: 0 + cash_onhand_months: 1 + admin_expense: 0.01 + capital_items: + depr_type: MACRS + depr_period: 5 + refurb: [0.] + cost_adjustment_parameters: + cost_year_adjustment_inflation: 0.025 + target_dollar_year: 2022 + finance_subgroups: + # The battery is the sole technology and the commodity stream: its per-year + # (degradation-aware) capacity factor drives the ProFAST cash-flow / LCOE analysis + # over the full 30-year plant life. + battery: + commodity: electricity + commodity_stream: battery + finance_groups: [profast_model] + technologies: [battery] diff --git a/examples/98_battery_degradation/tech_config.yaml b/examples/98_battery_degradation/tech_config.yaml new file mode 100644 index 000000000..0cd855c3a --- /dev/null +++ b/examples/98_battery_degradation/tech_config.yaml @@ -0,0 +1,44 @@ +name: technology_config +description: > + Single battery technology using the SimSES-backed BatteryPerformanceModel. No + control strategy is defined, so a passthrough controller forwards the + electricity_set_point profile to the battery as its dispatch command + (positive = discharge, negative = charge). +technologies: + battery: + performance_model: + model: BatteryPerformanceModel + cost_model: + model: ATBBatteryCostModel + model_inputs: + shared_parameters: + commodity: electricity + commodity_rate_units: kW + max_charge_rate: 2000.0 # kW + max_discharge_rate: 2000.0 # kW + charge_equals_discharge: true + max_capacity: 3854.0 # kWh (usable, 10-90 % of the 4817 kWh nominal pack) + max_soc_fraction: 0.9 # fraction (0-1) + min_soc_fraction: 0.1 # fraction (0-1) + init_soc_fraction: 0.9 # fraction (0-1) + charge_efficiency: 1.0 # fraction (0-1) + discharge_efficiency: 1.0 # fraction (0-1) + demand_profile: 0.0 # kW (unused; battery is driven by electricity_set_point) + series_count: 336 + parallel_count: 16 + battery_temperature_c: 25.0 + converter_efficiency: 0.96 + converter_max_power: 2400.0 # kW + performance_parameters: + cop: 3.0 + eol_soh_capacity: 0.8 # capacity SOH at which the battery is replaced (0.8 = 20 % fade) + cost_parameters: + cost_year: 2022 + energy_capex: 300.0 # $/kWh (ATB utility-scale battery, energy component) + power_capex: 100.0 # $/kW (ATB utility-scale battery, power component) + opex_fraction: 0.025 # annual fixed O&M as a fraction of total system CapEx + financial_parameters: + capital_items: + # Each end-of-life replacement (from the model's replacement_schedule) costs a + # full new battery. Depreciation defaults come from the plant-level capital_items. + replacement_cost_percent: 1.0 From 4d777623eb5f4ec42bacb5bf9489d8b7ec1f2105 Mon Sep 17 00:00:00 2001 From: Jared Thomas Date: Fri, 14 Aug 2026 12:05:17 -0600 Subject: [PATCH 17/18] expose all model inputs in example tech config --- examples/98_battery_degradation/tech_config.yaml | 4 ++++ h2integrate/storage/battery/battery_performance.py | 6 +++--- 2 files changed, 7 insertions(+), 3 deletions(-) diff --git a/examples/98_battery_degradation/tech_config.yaml b/examples/98_battery_degradation/tech_config.yaml index 0cd855c3a..df28739cc 100644 --- a/examples/98_battery_degradation/tech_config.yaml +++ b/examples/98_battery_degradation/tech_config.yaml @@ -14,6 +14,7 @@ technologies: shared_parameters: commodity: electricity commodity_rate_units: kW + commodity_amount_units: kW*h # amount units for the commodity (defaults to *h) max_charge_rate: 2000.0 # kW max_discharge_rate: 2000.0 # kW charge_equals_discharge: true @@ -23,6 +24,8 @@ technologies: init_soc_fraction: 0.9 # fraction (0-1) charge_efficiency: 1.0 # fraction (0-1) discharge_efficiency: 1.0 # fraction (0-1) + # Alternative to charge_efficiency + discharge_efficiency (mutually exclusive): + # round_trip_efficiency: 1.0 # fraction (0-1); sets charge/discharge eff to its sqrt demand_profile: 0.0 # kW (unused; battery is driven by electricity_set_point) series_count: 336 parallel_count: 16 @@ -31,6 +34,7 @@ technologies: converter_max_power: 2400.0 # kW performance_parameters: cop: 3.0 + deg_scale: 0.7056 # degradation scaling (calibrated: ~20 % capacity fade after 15 years) eol_soh_capacity: 0.8 # capacity SOH at which the battery is replaced (0.8 = 20 % fade) cost_parameters: cost_year: 2022 diff --git a/h2integrate/storage/battery/battery_performance.py b/h2integrate/storage/battery/battery_performance.py index b8534c4c3..dd489a92c 100644 --- a/h2integrate/storage/battery/battery_performance.py +++ b/h2integrate/storage/battery/battery_performance.py @@ -204,7 +204,7 @@ class BatteryPerformanceModelConfig(StoragePerformanceBaseConfig): discharge_efficiency: float | None = field(default=None, validator=range_val_or_none(0, 1)) round_trip_efficiency: float | None = field(default=None, validator=range_val_or_none(0, 1)) - _DEG_SCALE: float = field(default=0.7056, validator=range_val(0, 1)) + deg_scale: float = field(default=0.7056, validator=range_val(0, 1)) eol_soh_capacity: float = field(default=0.8, validator=range_val(0, 1)) # TODO convert from power and energy ratings (see math in chat) series_count: int = field(default=336, converter=int, validator=gt_zero) @@ -374,8 +374,8 @@ def compute(self, inputs, outputs, discrete_inputs=[], discrete_outputs=[]): "start_T": self.config.battery_temperature_c, }, degradation=DegradationModel( - calendar=ScaledLFPCalendarDegradation(self.config._DEG_SCALE), - cyclic=ScaledLFPCyclicDegradation(self.config._DEG_SCALE), + calendar=ScaledLFPCalendarDegradation(self.config.deg_scale), + cyclic=ScaledLFPCyclicDegradation(self.config.deg_scale), initial_soc=self.config.init_soc_fraction, initial_state=DegradationState(qloss_cal=1e-4), ), From 7e93b02cae1b2fe06c042efa64ff8b875f92c239 Mon Sep 17 00:00:00 2001 From: Jared Thomas Date: Fri, 14 Aug 2026 12:10:01 -0600 Subject: [PATCH 18/18] update pyproject.toml --- pyproject.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/pyproject.toml b/pyproject.toml index 26941689e..7d4cc620e 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -57,6 +57,7 @@ dependencies = [ "rich>=13.7.0", "ruamel.yaml", "scipy", + "simses>=2.1.1", "orbit-nrel", "floris", "hopp>=3.3.0",