Merge pull request #262 from camirmas/load_following_heuristic

Load following heuristic

hopp/simulation/hybrid_simulation.py

@@ -660,7 +660,14 @@ def simulate_power(self, project_life: int = 25, lifetime_sim=False):
         # Consolidate grid generation by copying over power and storage generation information
         if self.battery:
             self.grid.generation_profile_wo_battery = total_gen_before_battery
-        self.grid.simulate_grid_connection(hybrid_size_kw, total_gen, project_life, lifetime_sim, total_gen_max_feasible_year1)
+        self.grid.simulate_grid_connection(
+            hybrid_size_kw, 
+            total_gen, 
+            project_life, 
+            lifetime_sim,
+            total_gen_max_feasible_year1,
+            self.dispatch_builder.options
+        )
         self.grid.hybrid_nominal_capacity = hybrid_nominal_capacity
         self.grid.total_gen_max_feasible_year1 = total_gen_max_feasible_year1
         logger.info(f"Hybrid Peformance Simulation Complete. AEPs are {self.annual_energies}.")

hopp/simulation/technologies/dispatch/hybrid_dispatch_builder_solver.py

@@ -509,7 +509,22 @@ def battery_heuristic(self):
             prices = self.power_sources['grid'].dispatch.electricity_sell_price
             self.power_sources['battery'].dispatch.prices = prices
 
-        self.power_sources['battery'].dispatch.set_fixed_dispatch(tot_gen, grid_limit)
+        if  'load_following' in self.options.battery_dispatch:
+            # TODO: Look into how to define a system as load following or not in the config file
+            required_keys = ['desired_load']
+            if self.site.follow_desired_schedule:
+                # Get difference between baseload demand and power generation and control scenario variables
+                load_value = self.site.desired_schedule
+                load_difference =  [(load_value[x] - tot_gen[x]) for x in range(len(tot_gen))]
+                self.power_sources['battery'].dispatch.load_difference = load_difference
+            else:
+                raise ValueError(type(self).__name__ + " requires the following : desired_schedule")
+            # Adding goal_power for the simple battery heuristic method for power setpoint tracking 
+            goal_power = [load_value]*self.options.n_look_ahead_periods
+            ### Note: the inputs grid_limit and goal_power are in MW ###
+            self.power_sources['battery'].dispatch.set_fixed_dispatch(tot_gen, grid_limit, load_value)
+        else:
+            self.power_sources['battery'].dispatch.set_fixed_dispatch(tot_gen, grid_limit)
 
     @property
     def pyomo_model(self) -> pyomo.ConcreteModel:

hopp/simulation/technologies/dispatch/hybrid_dispatch_options.py

@@ -5,7 +5,8 @@
     SimpleBatteryDispatchHeuristic,
     SimpleBatteryDispatch,
     NonConvexLinearVoltageBatteryDispatch,
-    ConvexLinearVoltageBatteryDispatch
+    ConvexLinearVoltageBatteryDispatch,
+    HeuristicLoadFollowingDispatch,
 )
 
 
@@ -52,6 +53,10 @@ class HybridDispatchOptions:
 
             - **clustering_divisions** (dict, default={}): Custom number of averaging periods for classification metrics for data clustering. If empty, default values will be used.
 
+            - **use_higher_hours** bool (default = False): if True, the simulation will run extra hours analysis (must be used with load following)
+
+            - **higher_hours** (dict, default = {}): Higher hour count parameters: the value of power that must be available above the schedule and the number of hours in a row   
+
     """
     def __init__(self, dispatch_options: dict = None):
         self.solver: str = 'cbc'
@@ -74,6 +79,9 @@ def __init__(self, dispatch_options: dict = None):
         self.clustering_weights: dict = {}
         self.clustering_divisions: dict = {}
 
+        self.use_higher_hours: bool = False
+        self.higher_hours: dict = {}
+
         if dispatch_options is not None:
             for key, value in dispatch_options.items():
                 if hasattr(self, key):
@@ -103,7 +111,8 @@ def __init__(self, dispatch_options: dict = None):
             'heuristic': SimpleBatteryDispatchHeuristic,
             'simple': SimpleBatteryDispatch,
             'non_convex_LV': NonConvexLinearVoltageBatteryDispatch,
-            'convex_LV': ConvexLinearVoltageBatteryDispatch}
+            'convex_LV': ConvexLinearVoltageBatteryDispatch,
+            'load_following_heuristic': HeuristicLoadFollowingDispatch}
         if self.battery_dispatch in self._battery_dispatch_model_options:
             self.battery_dispatch_class = self._battery_dispatch_model_options[self.battery_dispatch]
             if 'heuristic' in self.battery_dispatch:

hopp/simulation/technologies/dispatch/power_storage/__init__.py

@@ -4,3 +4,4 @@
 from hopp.simulation.technologies.dispatch.power_storage.power_storage_dispatch import PowerStorageDispatch
 from hopp.simulation.technologies.dispatch.power_storage.simple_battery_dispatch import SimpleBatteryDispatch
 from hopp.simulation.technologies.dispatch.power_storage.simple_battery_dispatch_heuristic import SimpleBatteryDispatchHeuristic
+from hopp.simulation.technologies.dispatch.power_storage.heuristic_load_following_dispatch import HeuristicLoadFollowingDispatch

hopp/simulation/technologies/dispatch/power_storage/heuristic_load_following_dispatch.py

@@ -0,0 +1,62 @@
+from typing import Optional, List
+
+import pyomo.environ as pyomo
+from pyomo.environ import units as u
+import PySAM.BatteryStateful as BatteryModel
+import PySAM.Singleowner as Singleowner
+
+from hopp.simulation.technologies.dispatch.power_storage.simple_battery_dispatch_heuristic import SimpleBatteryDispatchHeuristic
+
+
+class HeuristicLoadFollowingDispatch(SimpleBatteryDispatchHeuristic):
+    """Operates the battery based on heuristic rules to meet the demand profile based power available from power generation profiles and 
+        power demand profile.
+
+    Currently, enforces available generation and grid limit assuming no battery charging from grid
+    """
+    def __init__(self,
+                 pyomo_model: pyomo.ConcreteModel,
+                 index_set: pyomo.Set,
+                 system_model: BatteryModel.BatteryStateful,
+                 financial_model: Singleowner.Singleowner,
+                 fixed_dispatch: Optional[List] = None,
+                 block_set_name: str = 'heuristic_load_following_battery',
+                 dispatch_options: Optional[dict] = None):
+        """
+
+        Args:
+            fixed_dispatch: list of normalized values [-1, 1] (Charging (-), Discharging (+))
+        """
+        super().__init__(
+            pyomo_model,
+            index_set,
+            system_model,
+            financial_model,
+            fixed_dispatch,
+            block_set_name,
+            dispatch_options
+        )
+
+    def set_fixed_dispatch(self, gen: list, grid_limit: list, goal_power: list):
+        """Sets charge and discharge power of battery dispatch using fixed_dispatch attribute and enforces available
+        generation and grid limits.
+
+        """
+        self.check_gen_grid_limit(gen, grid_limit)
+        self._set_power_fraction_limits(gen, grid_limit)
+        self._heuristic_method(gen, goal_power)
+        self._fix_dispatch_model_variables()
+
+    def _heuristic_method(self, gen, goal_power):
+        """ Enforces battery power fraction limits and sets _fixed_dispatch attribute
+            Sets the _fixed_dispatch based on goal_power and gen (power genration profile)
+        """
+        for t in self.blocks.index_set():
+            fd = (goal_power[t] - gen[t]) / self.maximum_power
+            if fd > 0.0:    # Discharging
+                if fd > self.max_discharge_fraction[t]:
+                    fd = self.max_discharge_fraction[t]
+            elif fd < 0.0:  # Charging
+                if -fd > self.max_charge_fraction[t]:
+                    fd = -self.max_charge_fraction[t]
+            self._fixed_dispatch[t] = fd        

hopp/simulation/technologies/dispatch/power_storage/simple_battery_dispatch_heuristic.py

@@ -1,3 +1,5 @@
+from typing import Optional, List, Dict
+
 import pyomo.environ as pyomo
 from pyomo.environ import units as u
 
@@ -17,9 +19,9 @@ def __init__(self,
                  index_set: pyomo.Set,
                  system_model: BatteryModel.BatteryStateful,
                  financial_model: Singleowner.Singleowner,
-                 fixed_dispatch: list = None,
+                 fixed_dispatch: Optional[List] = None,
                  block_set_name: str = 'heuristic_battery',
-                 dispatch_options: dict = None):
+                 dispatch_options: Optional[Dict] = None):
         """
 
         :param fixed_dispatch: list of normalized values [-1, 1] (Charging (-), Discharging (+))
@@ -60,8 +62,13 @@ def check_gen_grid_limit(self, gen: list, grid_limit: list):
             raise ValueError("grid_limit must be the same length as fixed_dispatch.")
 
     def _set_power_fraction_limits(self, gen: list, grid_limit: list):
-        """Set battery charge and discharge power fraction limits based on available generation and grid capacity,
-        respectively.
+        """
+        Set battery charge and discharge power fraction limits based on
+        available generation and grid capacity, respectively.
+
+        Args:
+            gen: generation Blocks
+            grid_limit: grid capacity
 
         NOTE: This method assumes that battery cannot be charged by the grid.
         """
@@ -71,28 +78,52 @@ def _set_power_fraction_limits(self, gen: list, grid_limit: list):
                                                                                        / self.maximum_power)
 
     @staticmethod
-    def enforce_power_fraction_simple_bounds(power_fraction) -> float:
-        """ Enforces simple bounds (0,1) for battery power fractions."""
-        if power_fraction > 1.0:
-            power_fraction = 1.0
+    def enforce_power_fraction_simple_bounds(power_fraction: float) -> float:
+        """
+        Enforces simple bounds (0, .9) for battery power fractions.
+        
+        Args:
+            power_fraction: power fraction from heuristic method
+
+        Returns:
+            bounded power fraction
+        """
+        if power_fraction > 0.9:
+            power_fraction = 0.9
         elif power_fraction < 0.0:
             power_fraction = 0.0
         return power_fraction
 
-    def update_soc(self, power_fraction, soc0) -> float:
+    def update_soc(self, power_fraction: float, soc0: float) -> float:
+        """
+        Updates SOC based on power fraction threshold (0.1).
+        
+        Args:
+            power_fraction: power fraction from heuristic method. Below threshold
+                is charging, above is discharging
+            soc0: initial SOC
+        
+        Returns:
+            Updated SOC.
+        """
         if power_fraction > 0.0:
             discharge_power = power_fraction * self.maximum_power
             soc = soc0 - self.time_duration[0] * (1/(self.discharge_efficiency/100.) * discharge_power) / self.capacity
         elif power_fraction < 0.0:
-            charge_power = - power_fraction * self.maximum_power
+            charge_power = -power_fraction * self.maximum_power
             soc = soc0 + self.time_duration[0] * (self.charge_efficiency / 100. * charge_power) / self.capacity
         else:
             soc = soc0
-        soc = max(0, min(1, soc))
+
+        min_soc = self._system_model.value("minimum_SOC") / 100
+        max_soc = self._system_model.value("maximum_SOC") / 100
+
+        soc = max(min_soc, min(max_soc, soc))
+
         return soc
 
     def _heuristic_method(self, _):
-        """ Does specific heuristic method to fix battery dispatch."""
+        """Does specific heuristic method to fix battery dispatch."""
         self._enforce_power_fraction_limits()
 
     def _enforce_power_fraction_limits(self):

hopp/simulation/technologies/grid.py

@@ -1,4 +1,4 @@
-from typing import Iterable, List, Sequence, Optional, Union
+from typing import Iterable, List, Sequence, Optional, Union, TYPE_CHECKING
 
 import numpy as np
 from attrs import define, field
@@ -11,7 +11,10 @@
 from hopp.simulation.technologies.financial import FinancialModelType, CustomFinancialModel
 from hopp.type_dec import NDArrayFloat
 from hopp.utilities.validators import gt_zero
+from hopp.utilities.log import hybrid_logger as logger
 
+if TYPE_CHECKING:
+    from hopp.simulation.technologies.dispatch.hybrid_dispatch_options import HybridDispatchOptions
 
 @define
 class GridConfig(BaseClass):
@@ -91,7 +94,8 @@ def simulate_grid_connection(
         total_gen: Union[List[float], NDArrayFloat], 
         project_life: int, 
         lifetime_sim: bool, 
-        total_gen_max_feasible_year1: Union[List[float], NDArrayFloat]
+        total_gen_max_feasible_year1: Union[List[float], NDArrayFloat],
+        dispatch_options: Optional["HybridDispatchOptions"] = None
     ):
         """
         Sets up and simulates hybrid system grid connection. Additionally,
@@ -108,6 +112,8 @@ def simulate_grid_connection(
                 data is repeated
             total_gen_max_feasible_year1: Maximum generation profile of the hybrid
                 system (for capacity payments) [kWh]
+            dispatch_options: Hybrid dispatch options class, deliminates if the higher
+                power analysis for frequency regulation is run
 
         """
         if self.site.follow_desired_schedule:
@@ -125,8 +131,65 @@ def simulate_grid_connection(
             self.schedule_curtailed = np.array([gen - schedule if gen > schedule else 0. for (gen, schedule) in
                                             zip(total_gen, lifetime_schedule)])
             self.schedule_curtailed_percentage = sum(self.schedule_curtailed)/sum(lifetime_schedule)
+
+            # NOTE: This is currently only happening for load following, would be good to make it more general
+            #           i.e. so that this analysis can be used when load following isn't being used (without storage)
+            #           for comparison 
+            N_hybrid = len(self.generation_profile)
+
+            final_power_production = total_gen
+            schedule = [x for x in lifetime_schedule]
+            hybrid_power = [(final_power_production[x] - (schedule[x]*0.95)) for x in range(len(final_power_production))]
+
+            load_met = len([i for i in hybrid_power if i  >= 0])
+            self.time_load_met = 100 * load_met/N_hybrid
+
+            final_power_array = np.array(final_power_production)
+            power_met = np.where(final_power_array > schedule, schedule, final_power_array)
+            self.capacity_factor_load = np.sum(power_met) / np.sum(schedule) * 100
+
+            logger.info('Percent of time firm power requirement is met: ', np.round(self.time_load_met,2))
+            logger.info('Percent total firm power requirement is satisfied: ', np.round(self.capacity_factor_load,2))
+
+            ERS_keys = ['min_regulation_hours', 'min_regulation_power']
+            if dispatch_options is not None and dispatch_options.use_higher_hours:
+                """
+                Frequency regulation analysis for providing essential reliability services (ERS) availability operating case:
+                        Finds how many hours (in the group specified group size above the specified minimum
+                        power requirement) that the system has available to extra power that could be used to 
+                        provide ERS
+                Args:
+                    :param dispatch_options: need additional ERS arguments
+                                        'min_regulation_hours': minimum size of hours in a group to be considered for ERS (>= 1)
+                                        'min_regulation_power': minimum power available over the whole group of hours to be 
+                                                considered for ERS (> 0, in kW)
+
+                :returns: total_number_hours
+
+                """
+
+                # Performing frequency regulation analysis:
+                #    finding how many groups of hours satisfiy the ERS minimum power requirement
+                min_regulation_hours = dispatch_options.higher_hours['min_regulation_hours']
+                min_regulation_power = dispatch_options.higher_hours['min_regulation_power']
+
+                frequency_power_array = np.array(hybrid_power)
+                frequency_test = np.where(frequency_power_array > min_regulation_power, frequency_power_array, 0)
+                mask = (frequency_test!=0).astype(int)
+                padded_mask = np.pad(mask,(1,), "constant")
+                edge_mask = padded_mask[1:] - padded_mask[:-1]  # finding the difference between each array value
+
+                group_starts = np.where(edge_mask == 1)[0]
+                group_stops = np.where(edge_mask == -1)[0]
+
+                # Find groups and drop groups that are too small
+                groups = [group for group in zip(group_starts,group_stops) if ((group[1]-group[0]) >= min_regulation_hours)]
+                group_lengths = [len(final_power_production[group[0]:group[1]]) for group in groups]
+                self.total_number_hours = sum(group_lengths)
+
+                logger.info('Total number of hours available for ERS: ', np.round(self.total_number_hours,2))
         else:
-            self.generation_profile = list(total_gen)
+            self.generation_profile = total_gen
 
         self.total_gen_max_feasible_year1 = np.array(total_gen_max_feasible_year1)
         self.system_capacity_kw = hybrid_size_kw  # TODO: Should this be interconnection limit?

hopp/tools/dispatch/plot_tools.py

@@ -303,6 +303,13 @@ def plot_generation_profile(hybrid: HybridSimulation,
     ax.xaxis.set_ticks(list(range(start, end, hybrid.site.n_periods_per_day)))
     plt.grid()
     ax1 = plt.gca()
+
+    # Load following, if applicable
+    if hybrid.site.follow_desired_schedule:
+        desired_load = [p for p in hybrid.site.desired_schedule[time_slice]]
+        ax1.plot(time, desired_load, 'b--', label='Desired Load')
+        ax1.set_ylabel('Desired Load', fontsize=font_size)
+
     ax1.legend(fontsize=font_size-2, loc='upper left')
     ax1.set_ylabel('Power (MW)', fontsize=font_size)
 
@@ -335,6 +342,9 @@ def plot_generation_profile(hybrid: HybridSimulation,
     plt.xlabel('Time (hours)', fontsize=font_size)
     plt.title('Net Generation', fontsize=font_size)
 
+
+    plt.xlabel('Time (hours)', fontsize=font_size)
+    plt.title('Net Generation', fontsize=font_size)
     plt.tight_layout()
 
     if plot_filename is not None: