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/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..75c0d6d04 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: @@ -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/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/run_battery_degradation.py b/examples/98_battery_degradation/run_battery_degradation.py new file mode 100644 index 000000000..906858ab6 --- /dev/null +++ b/examples/98_battery_degradation/run_battery_degradation.py @@ -0,0 +1,313 @@ +""" +Example 98: Battery degradation through the standard H2Integrate interface. + +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 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) + 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) + +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 + +import numpy as np +import pandas as pd +import matplotlib + + +matplotlib.use("Agg") +import matplotlib.pyplot as plt + +from h2integrate import H2IntegrateModel + + +HERE = Path(__file__).parent + +# --------------------------------------------------------------------------- +# Build and set up the H2Integrate model (standard interface) +# --------------------------------------------------------------------------- +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 + +# --------------------------------------------------------------------------- +# 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. +# --------------------------------------------------------------------------- +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), +) + +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)") + +# --------------------------------------------------------------------------- +# Drive the battery via its set-point (passthrough controller -> command value) +# and run the model. +# --------------------------------------------------------------------------- +model.prob.set_val("battery.electricity_set_point", set_point, units="kW") +model.run() + +# --------------------------------------------------------------------------- +# 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. +# --------------------------------------------------------------------------- +index = pd.date_range("2026-01-01", periods=n_steps, freq=f"{int(dt)}s") +df = pd.DataFrame( + { + "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, +) + +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 + +# --------------------------------------------------------------------------- +# Plot 1: AC & DC power (first 7 days) +# --------------------------------------------------------------------------- +fig, ax = plt.subplots(figsize=(12, 3)) +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}")) +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 3: Battery temperature (first 7 days) +# --------------------------------------------------------------------------- +fig, ax = plt.subplots(figsize=(12, 3)) +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 4a: Terminal voltage (first 7 days) +# --------------------------------------------------------------------------- +fig, ax = plt.subplots(figsize=(12, 3)) +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 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]") +fig.tight_layout() +fig.savefig(HERE / "plot_voltage.png", dpi=150) +print("Saved plot_voltage.png") + +# --------------------------------------------------------------------------- +# 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 [%]", +) +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.axhline(20, color="steelblue", linewidth=0.8, linestyle="--", alpha=0.6) +ax1.set_title( + 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") +fig.tight_layout() +fig.savefig(HERE / "plot_degradation.png", dpi=150) +print("Saved plot_degradation.png") + +# --------------------------------------------------------------------------- +# Plot 6: Losses and heat (daily mean) +# --------------------------------------------------------------------------- +fig, ax = plt.subplots(figsize=(12, 3)) +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() +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/examples/98_battery_degradation/tech_config.yaml b/examples/98_battery_degradation/tech_config.yaml new file mode 100644 index 000000000..df28739cc --- /dev/null +++ b/examples/98_battery_degradation/tech_config.yaml @@ -0,0 +1,48 @@ +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 + 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 + 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) + # 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 + battery_temperature_c: 25.0 + converter_efficiency: 0.96 + 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 + 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 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/supported_models.py b/h2integrate/core/supported_models.py index d80d5e1ec..950760994 100644 --- a/h2integrate/core/supported_models.py +++ b/h2integrate/core/supported_models.py @@ -147,6 +147,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/core/test/test_framework.py b/h2integrate/core/test/test_framework.py index f901111b8..5464269f7 100644 --- a/h2integrate/core/test/test_framework.py +++ b/h2integrate/core/test/test_framework.py @@ -358,20 +358,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/__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/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 diff --git a/h2integrate/storage/battery/battery_performance.py b/h2integrate/storage/battery/battery_performance.py new file mode 100644 index 000000000..dd489a92c --- /dev/null +++ b/h2integrate/storage/battery/battery_performance.py @@ -0,0 +1,556 @@ +""" +# 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 +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.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, + 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 +from h2integrate.storage.storage_baseclass import ( + StoragePerformanceBase, + StoragePerformanceBaseConfig, +) + + +class LFP280Ah(SonyLFP): + """280 Ah / 3.2 V prismatic LFP cell scaled from the SonyLFP OCV/resistance curves. + + 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__() + 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 __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 * 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: + 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 * 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, + ) -> float: + delta_fec = half_cycle.full_equivalent_cycles + if delta_fec == 0.0: + return 0.0 + 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: + 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 * 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 + + +@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)) + + 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) + 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) + + 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 = ( + 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"), + 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 + + super().setup() + + self.add_discrete_input( + "solar_resource_data", + val={}, + 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", + ) + + # 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=[]): + """Run the storage performance model.""" + self.current_soc = self.config.init_soc_fraction + + 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] + + # H2I dispatch command: positive = discharge, negative = charge (commodity_rate_units) + power_profile = inputs[f"{self.commodity}_command_value"] + + # --------------------------------------------------------------------------- + # Battery pack + inverter (fixed Megapack-style topology, from config) + # --------------------------------------------------------------------------- + cell = LFP280Ah() + + # TODO check sizing + battery = Battery( + cell=cell, + circuit=(self.config.series_count, self.config.parallel_count), + initial_states={ + "start_soc": self.config.init_soc_fraction, + "start_T": self.config.battery_temperature_c, + }, + degradation=DegradationModel( + 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), + ), + ) + + 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, + ) + + # --------------------------------------------------------------------------- + # Thermal model: constant ambient, battery registered as thermal node + # #TODO check battery temp ambient + # --------------------------------------------------------------------------- + thermal = AmbientThermalModel( + T_ambient=self.config.battery_temperature_c, components=[battery] + ) + + # --------------------------------------------------------------------------- + # Simulation loop + # --------------------------------------------------------------------------- + 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): + # 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 + + ############# + + # 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). + 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"] = 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) + + # --- 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: + # 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 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..7f0641c82 --- /dev/null +++ b/h2integrate/storage/battery/test/test_battery_performance.py @@ -0,0 +1 @@ +# TODO degradation: tests for battery performance 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",