updated storage autosizing model

Author: elenya-grant

h2integrate/storage/simple_storage_auto_sizing.py

@@ -1,13 +1,16 @@
 import numpy as np
 from attrs import field, define
 
-from h2integrate.core.utilities import BaseConfig, merge_shared_inputs
-from h2integrate.core.validators import range_val, range_val_or_none
-from h2integrate.core.model_baseclasses import PerformanceModelBaseClass
+from h2integrate.core.utilities import merge_shared_inputs
+from h2integrate.core.validators import range_val_or_none
+from h2integrate.storage.storage_baseclass import (
+    StoragePerformanceBase,
+    StoragePerformanceBaseConfig,
+)
 
 
 @define(kw_only=True)
-class StorageSizingModelConfig(BaseConfig):
+class StorageSizingModelConfig(StoragePerformanceBaseConfig):
     """Configuration class for the StorageAutoSizingModel.
 
     Attributes:
@@ -38,8 +41,8 @@ class StorageSizingModelConfig(BaseConfig):
     commodity: str = field(converter=(str.strip, str.lower))
     commodity_rate_units: str = field(converter=str.strip)
 
-    min_soc_fraction: float = field(validator=range_val(0, 1))
-    max_soc_fraction: float = field(validator=range_val(0, 1))
+    # min_soc_fraction: float = field(validator=range_val(0, 1))
+    # max_soc_fraction: float = field(validator=range_val(0, 1))
 
     # TODO: add in logic for having different discharge rate
     # charge_equals_discharge: bool = field(default=True)
@@ -80,7 +83,7 @@ def __attrs_post_init__(self):
             self.commodity_amount_units = f"({self.commodity_rate_units})*h"
 
 
-class StorageAutoSizingModel(PerformanceModelBaseClass):
+class StorageAutoSizingModel(StoragePerformanceBase):
     """Performance model that calculates the storage charge rate and capacity needed
     to either:
 
@@ -136,23 +139,6 @@ def setup(self):
 
         super().setup()
 
-        # Inputs
-        self.add_input(
-            f"{self.commodity}_demand",
-            val=self.config.demand_profile,
-            shape=self.n_timesteps,
-            units=self.commodity_rate_units,
-            desc=f"{self.commodity} demand profile timeseries",
-        )
-
-        self.add_input(
-            f"{self.commodity}_in",
-            val=0.0,
-            shape=self.n_timesteps,
-            units=self.commodity_rate_units,
-            desc=f"{self.commodity} input timeseries from production to storage",
-        )
-
         # Capacity outputs
         self.add_output(
             "storage_capacity",
@@ -175,87 +161,6 @@ def setup(self):
             units=self.commodity_rate_units,
         )
 
-        self.add_output(
-            "storage_duration",
-            units=f"({self.commodity_amount_units})/({self.commodity_rate_units})",
-            desc="Estimated storage duration based on max capacity and discharge rate",
-        )
-
-        # Storage performance outputs
-        self.add_output(
-            "SOC",
-            val=0.0,
-            shape=self.n_timesteps,
-            units="percent",
-            desc="State of charge of storage",
-        )
-
-        self.add_output(
-            f"storage_{self.commodity}_discharge",
-            val=0.0,
-            shape=self.n_timesteps,
-            units=self.commodity_rate_units,
-            desc=f"{self.commodity} output from storage only",
-        )
-
-        self.add_output(
-            f"storage_{self.commodity}_charge",
-            val=0.0,
-            shape=self.n_timesteps,
-            units=self.commodity_rate_units,
-            desc=f"{self.commodity} input to storage only",
-        )
-
-        self.add_output(
-            f"storage_{self.commodity}_out",
-            val=0.0,
-            shape=self.n_timesteps,
-            units=self.commodity_rate_units,
-            desc=f"{self.commodity} input and output from storage",
-        )
-
-        # Combined storage and input outputs
-        self.add_output(
-            f"unmet_{self.commodity}_demand_out",
-            val=0.0,
-            shape=self.n_timesteps,
-            units=self.commodity_rate_units,
-            desc=f"Unmet {self.commodity} demand",
-        )
-
-        self.add_output(
-            f"unused_{self.commodity}_out",
-            val=0.0,
-            shape=self.n_timesteps,
-            units=self.commodity_rate_units,
-            desc="Unused generated commodity",
-        )
-
-        self.dt_hr = self.dt / 3600  # convert from seconds to hours
-
-        # create a variable to determine whether we are using feedback control
-        # for this technology
-        using_feedback_control = False
-        # create inputs for pyomo control model
-        if "tech_to_dispatch_connections" in self.options["plant_config"]:
-            self.tech_group_name = self.pathname.split(".")
-            for _source_tech, intended_dispatch_tech in self.options["plant_config"][
-                "tech_to_dispatch_connections"
-            ]:
-                if any(intended_dispatch_tech in name for name in self.tech_group_name):
-                    self.add_discrete_input("pyomo_dispatch_solver", val=dummy_function)
-                    using_feedback_control = True
-                    break
-
-        if not using_feedback_control:
-            # using an open-loop storage controller
-            self.add_input(
-                f"{self.commodity}_set_point",
-                val=0.0,
-                shape=self.n_timesteps,
-                units=self.commodity_rate_units,
-            )
-
     def compute(self, inputs, outputs, discrete_inputs=[], discrete_outputs=[]):
         """
         Part 0: get demand profile based on user input parameters:
@@ -357,221 +262,28 @@ def compute(self, inputs, outputs, discrete_inputs=[], discrete_outputs=[]):
             [self.config.min_soc_fraction, commodity_storage_soc[0] / rated_storage_capacity]
         )
 
-        # 2. Simulate the storage performance using the `simulate()`
-        if "pyomo_dispatch_solver" in discrete_inputs:
-            # Simulate the storage with provided dispatch inputs
-            dispatch = discrete_inputs["pyomo_dispatch_solver"]
-            # kwargs are tech-specific inputs to the simulate() method
-            kwargs = {
-                "charge_rate": storage_max_fill_rate,
-                "discharge_rate": storage_max_empty_rate,
-                "storage_capacity": rated_storage_capacity,
-            }
-
-            # Make dictionary of inputs containing information to pass to the controller
-            # (such as demand profile, charge rate, and storage capacity)
-            inputs_adjusted = dict(inputs.items())
-            inputs_adjusted["storage_capacity"] = np.array([rated_storage_capacity])
-            inputs_adjusted["max_charge_rate"] = np.array([storage_max_fill_rate])
-
-            if self.config.set_demand_as_avg_commodity_in:
-                inputs_adjusted[f"{self.commodity}_demand"] = commodity_demand
-            storage_commodity_out, soc = dispatch(self.simulate, kwargs, inputs_adjusted)
-
-        else:
-            # Simulate the storage with provided inputs using dispatch commands from
-            # an open-loop controller. The commodity_set_point should come from an
-            # open-loop controller. commodity_set_point is negative when commanding
-            # storage to charge and positive when commanding storage to discharge
-            storage_commodity_out, soc = self.simulate(
-                inputs[f"{self.commodity}_set_point"],
-                storage_max_fill_rate,
-                storage_max_empty_rate,
-                rated_storage_capacity,
-            )
-
-        # soc output from simulating storage is represented as a percentage
-        storage_commodity_out = np.array(storage_commodity_out)
-
-        # 3. Calculate outputs
-
-        # calculate combined commodity out from inflow source and storage
-        # (note: storage_commodity_out is negative when charging)
-        combined_commodity_out = inputs[f"{self.commodity}_in"] + storage_commodity_out
-
-        # find the total commodity out to meet demand
-        total_commodity_out = np.minimum(commodity_demand, combined_commodity_out)
-
-        # determine how much of the inflow commodity was unused
-        unused_commodity = np.maximum(0, combined_commodity_out - commodity_demand)
-
-        # determine how much demand was not met
-        unmet_demand = np.maximum(0, commodity_demand - combined_commodity_out)
-
-        # Output the storage performance outputs
-        outputs[f"storage_{self.commodity}_charge"] = np.where(
-            storage_commodity_out < 0, storage_commodity_out, 0
-        )
-        outputs[f"storage_{self.commodity}_discharge"] = np.where(
-            storage_commodity_out > 0, storage_commodity_out, 0
-        )
-        outputs[f"unmet_{self.commodity}_demand_out"] = unmet_demand
-        outputs[f"unused_{self.commodity}_out"] = unused_commodity
-        outputs[f"storage_{self.commodity}_out"] = storage_commodity_out
-        outputs["SOC"] = soc
-
-        # commodity_out is the commodity_in - charge_storage + discharge_storage
-        outputs[f"{self.commodity}_out"] = total_commodity_out
-
         # Output the calculated storage sizes (charge rate and capacity)
         outputs["max_charge_rate"] = storage_max_fill_rate
         outputs["max_discharge_rate"] = storage_max_empty_rate
         outputs["storage_capacity"] = rated_storage_capacity
-        outputs["storage_duration"] = outputs["storage_capacity"] / outputs["max_discharge_rate"]
 
-        # The rated_commodity_production is based on the discharge rate
-        # (which is assumed equal to the charge rate)
-        outputs[f"rated_{self.commodity}_production"] = storage_max_fill_rate
+        # 2. Simulate the storage performance using the `simulate()`
 
-        # The total_commodity_produced is the sum of the commodity discharged from storage
-        outputs[f"total_{self.commodity}_produced"] = np.sum(total_commodity_out)
+        # Make dictionary of inputs containing information to pass to the controller
+        # (such as demand profile, charge rate, and storage capacity)
+        inputs_adjusted = dict(inputs.items())
+        if self.config.set_demand_as_avg_commodity_in:
+            inputs_adjusted[f"{self.commodity}_demand"] = commodity_demand
 
-        # Adjust the total_commodity_produced to a year-long simulation
-        outputs[f"annual_{self.commodity}_produced"] = outputs[
-            f"total_{self.commodity}_produced"
-        ] * (1 / self.fraction_of_year_simulated)
+        if "pyomo_dispatch_solver" in discrete_inputs:
+            inputs_adjusted["storage_capacity"] = np.array([rated_storage_capacity])
+            inputs_adjusted["max_charge_rate"] = np.array([storage_max_fill_rate])
 
-        # Capacity factor is total discharged commodity / maximum discharged commodity possible
-        outputs["capacity_factor"] = outputs[f"total_{self.commodity}_produced"] / (
-            outputs[f"rated_{self.commodity}_production"] * self.n_timesteps
+        outputs = self.run_storage(
+            storage_max_fill_rate,
+            storage_max_empty_rate,
+            rated_storage_capacity,
+            inputs_adjusted,
+            outputs,
+            discrete_inputs,
         )
-
-    def simulate(
-        self,
-        storage_dispatch_commands: list,
-        charge_rate: float,
-        discharge_rate: float,
-        storage_capacity: float,
-        sim_start_index: int = 0,
-    ):
-        """Run the storage model over a control window of ``n_control_window`` timesteps.
-
-        Iterates through ``storage_dispatch_commands`` one timestep at a time.
-        A negative command requests charging; a positive command requests
-        discharging.  Each command is clipped to the most restrictive of three
-        limits before it is applied:
-
-        1. **SOC headroom** - the remaining capacity (charge) or remaining
-           stored commodity (discharge), converted to a rate via
-           ``storage_capacity / dt_hr``.
-        2. **Hardware rate limit** - ``charge_rate`` or ``discharge_rate``,
-           divided by the corresponding efficiency so the limit is expressed
-           in pre-efficiency rate units.
-        3. **Commanded magnitude** - the absolute value of the dispatch command
-           itself (we never exceed what was asked for).
-
-        After clipping, the result is scaled by the charge or discharge
-        efficiency to obtain the actual commodity flow into or out of the
-        storage, and the SOC is updated accordingly.
-
-        This method is separated from ``compute()`` so the Pyomo dispatch
-        controller can call it directly to evaluate candidate schedules.
-
-        Args:
-            storage_dispatch_commands (array_like[float]):
-                Dispatch set-points for each timestep in ``commodity_rate_units``.
-                Negative values command charging; positive values command
-                discharging.  Length must equal ``config.n_control_window``.
-            charge_rate (float):
-                Maximum commodity input rate to storage in
-                ``commodity_rate_units`` (before charge efficiency is applied).
-            discharge_rate (float):
-                Maximum commodity output rate from storage in
-                ``commodity_rate_units`` (before discharge efficiency is applied).
-            storage_capacity (float):
-                Rated storage capacity in ``commodity_amount_units``.
-            sim_start_index (int, optional):
-                Starting index for writing into persistent output arrays.
-                Defaults to 0.
-
-        Returns:
-            tuple[np.ndarray, np.ndarray]
-                storage_commodity_out_timesteps :
-                    Commodity flow per timestep in ``commodity_rate_units``.
-                    Positive = discharge (commodity leaving storage),
-                    negative = charge (commodity entering storage).
-                soc_timesteps :
-                    State of charge at the end of each timestep, in percent
-                    (0-100).
-        """
-
-        n = len(storage_dispatch_commands)
-        storage_commodity_out_timesteps = np.zeros(n)
-        soc_timesteps = np.zeros(n)
-
-        # Early return when storage cannot operate: zero capacity or both
-        # charge and discharge rates are zero.
-        if storage_capacity <= 0 or (charge_rate <= 0 and discharge_rate <= 0):
-            soc_timesteps[:] = self.current_soc * 100.0
-            return storage_commodity_out_timesteps, soc_timesteps
-
-        # Pre-compute scalar constants to avoid repeated attribute lookups
-        # and redundant divisions inside the per-timestep loop.
-        charge_eff = self.config.charge_efficiency
-        discharge_eff = self.config.discharge_efficiency
-        soc_max = self.config.max_soc_fraction
-        soc_min = self.config.min_soc_fraction
-
-        commands = np.asarray(storage_dispatch_commands, dtype=float)
-        soc = float(self.current_soc)
-
-        for t, cmd in enumerate(commands):
-            if cmd < 0.0:
-                # --- Charging ---
-                # headroom: how much more commodity the storage can accept,
-                # expressed as a rate (commodity_rate_units).
-                headroom = (soc_max - soc) * storage_capacity / self.dt_hr
-
-                # Clip to the most restrictive limit, then apply efficiency.
-                # max(0, ...) guards against negative headroom when SOC
-                # slightly exceeds soc_max.
-                # correct headroom to not include charge_eff.
-                actual_charge = max(0.0, min(headroom / charge_eff, charge_rate, -cmd)) * charge_eff
-
-                # Update SOC (actual_charge is in post-efficiency units)
-                soc += actual_charge / storage_capacity
-
-                # Update the amount of commodity used to charge from the input stream
-                # If charge_eff<1, more commodity is pulled from the input stream than
-                # the commodity that goes into the storage.
-                storage_commodity_out_timesteps[t] = -actual_charge / charge_eff
-            else:
-                # --- Discharging ---
-                # headroom: how much commodity can still be drawn before
-                # hitting the minimum SOC, expressed as a rate.
-                headroom = (soc - soc_min) * storage_capacity / self.dt_hr
-
-                # Clip to the most restrictive limit without applied efficiency.
-                # Discharge efficiency losses occur as energy leaves storage.
-                actual_discharge = max(
-                    0.0, min(headroom, discharge_rate / discharge_eff, cmd / discharge_eff)
-                )
-
-                # Update SOC (actual_discharge is before efficiency losses are applied.)
-                soc -= actual_discharge / storage_capacity
-
-                # If discharge_eff<1, then less commodity is output from the storage
-                # than the commodity discharged from storage
-                storage_commodity_out_timesteps[t] = actual_discharge * discharge_eff
-
-            soc_timesteps[t] = soc * 100.0
-
-        # Persist the final SOC so subsequent simulate() calls (e.g. from the
-        # Pyomo controller across rolling windows) start where we left off.
-        self.current_soc = soc
-        return storage_commodity_out_timesteps, soc_timesteps
-
-
-def dummy_function():
-    # this function is required for initializing the pyomo control input and nothing else
-    pass