Repository: NREL/HOPP

Issue:
Add load following dispatch example
<!--
IMPORTANT NOTES

Thank you for taking the time to suggest a feature. Before submitting,
please reread your description to ensure that other readers can reasonably
understand the motivation and proposed solution.

This form is written in GitHub's Markdown format. For a reference on this type
of syntax, see GitHub's documentation:
https://docs.github.com/en/get-started/writing-on-github/getting-started-with-writing-and-formatting-on-github/basic-writing-and-formatting-syntax

This template contains guidance for your submission within the < ! - -, - - > blocks.
These are comments in HTML syntax and will not appear in the submission.
Be sure to use the "Preview" feature on GitHub to ensure your submission is formatted as intended.

When including code snippets, please paste the text itself and wrap the code block with
ticks (see the other character on the tilde ~ key in a US keyboard) to format it as code.
For example, Python code should be wrapped in ticks like this:
```python
def a_func():
    return 1

a = 1
b = a_func()
print(a + b)
```
This is preferred over screen shots since it is searchable and others can copy/paste
the text to run it.
-->

# Add Load Following Dispatch Example
<!--
High level description of the feature request including motivation and background.
-->
#262 Introduces a new feature: heuristic load following dispatch. Create an example notebook in the style of HOPP v2.0+ examples.

## Proposed solution
<!--
Here's an opportunity to prototype a feature. Please include pseudocode, images, or
any other visual aids to communicate the idea.
-->
There's an existing [example](https://github.com/NREL/HOPP/blob/00536414b520023c8af71ba4bbf296594b652cbf/examples/simulate_hybrid_load_following.py) that has become out of date. Rework this into a notebook.

Discussion:
@bayc: @genevievestarke Do we have an example for this?@elenya-grant: completed in PR #391 

Solution patch:
diff --git a/hopp/simulation/hybrid_simulation.py b/hopp/simulation/hybrid_simulation.py
index dcad3619..457a9c76 100644
--- a/hopp/simulation/hybrid_simulation.py
+++ b/hopp/simulation/hybrid_simulation.py
@@ -660,7 +660,14 @@ class HybridSimulation(BaseClass):
         # 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}.")
diff --git a/hopp/simulation/technologies/dispatch/hybrid_dispatch_builder_solver.py b/hopp/simulation/technologies/dispatch/hybrid_dispatch_builder_solver.py
index 1c336553..cbc5d694 100644
--- a/hopp/simulation/technologies/dispatch/hybrid_dispatch_builder_solver.py
+++ b/hopp/simulation/technologies/dispatch/hybrid_dispatch_builder_solver.py
@@ -509,7 +509,22 @@ class HybridDispatchBuilderSolver:
             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:
diff --git a/hopp/simulation/technologies/dispatch/hybrid_dispatch_options.py b/hopp/simulation/technologies/dispatch/hybrid_dispatch_options.py
index b6a5fb4e..69c26f1d 100644
--- a/hopp/simulation/technologies/dispatch/hybrid_dispatch_options.py
+++ b/hopp/simulation/technologies/dispatch/hybrid_dispatch_options.py
@@ -5,7 +5,8 @@ from hopp.simulation.technologies.dispatch.power_storage import (
     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 @@ class HybridDispatchOptions:
         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 @@ class HybridDispatchOptions:
             '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:
diff --git a/hopp/simulation/technologies/dispatch/power_storage/__init__.py b/hopp/simulation/technologies/dispatch/power_storage/__init__.py
index 15951a51..73592b77 100644
--- a/hopp/simulation/technologies/dispatch/power_storage/__init__.py
+++ b/hopp/simulation/technologies/dispatch/power_storage/__init__.py
@@ -4,3 +4,4 @@ from hopp.simulation.technologies.dispatch.power_storage.one_cycle_battery_dispa
 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
diff --git a/hopp/simulation/technologies/dispatch/power_storage/heuristic_load_following_dispatch.py b/hopp/simulation/technologies/dispatch/power_storage/heuristic_load_following_dispatch.py
new file mode 100644
index 00000000..2d6ae049
--- /dev/null
+++ b/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        
\ No newline at end of file
diff --git a/hopp/simulation/technologies/dispatch/power_storage/simple_battery_dispatch_heuristic.py b/hopp/simulation/technologies/dispatch/power_storage/simple_battery_dispatch_heuristic.py
index 09077aac..8c576069 100644
--- a/hopp/simulation/technologies/dispatch/power_storage/simple_battery_dispatch_heuristic.py
+++ b/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 @@ class SimpleBatteryDispatchHeuristic(SimpleBatteryDispatch):
                  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 @@ class SimpleBatteryDispatchHeuristic(SimpleBatteryDispatch):
             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 @@ class SimpleBatteryDispatchHeuristic(SimpleBatteryDispatch):
                                                                                        / 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):
diff --git a/hopp/simulation/technologies/grid.py b/hopp/simulation/technologies/grid.py
index 70105211..a080061a 100644
--- a/hopp/simulation/technologies/grid.py
+++ b/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.base import BaseClass
 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 @@ class Grid(PowerSource):
         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 @@ class Grid(PowerSource):
                 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 @@ class Grid(PowerSource):
             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?
diff --git a/hopp/tools/dispatch/plot_tools.py b/hopp/tools/dispatch/plot_tools.py
index 7581eee8..3fcf119a 100644
--- a/hopp/tools/dispatch/plot_tools.py
+++ b/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:


Test patch:
diff --git a/tests/hopp/test_dispatch.py b/tests/hopp/test_dispatch.py
index 026179b7..b42c0c41 100644
--- a/tests/hopp/test_dispatch.py
+++ b/tests/hopp/test_dispatch.py
@@ -6,7 +6,7 @@ from pyomo.opt import TerminationCondition
 from pyomo.util.check_units import assert_units_consistent
 
 from hopp.simulation import HoppInterface
-from hopp.simulation.technologies.sites import SiteInfo
+from hopp.simulation.technologies.sites import SiteInfo, flatirons_site
 from hopp.simulation.technologies.financial.custom_financial_model import CustomFinancialModel
 from hopp.simulation.technologies.wind.wind_plant import WindPlant, WindConfig
 from hopp.simulation.technologies.pv.pv_plant import PVPlant, PVConfig
@@ -992,3 +992,35 @@ def test_simple_battery_dispatch_lifecycle_limit(site):
     assert sum(battery.dispatch.discharge_power) > 0.0
     assert (sum(battery.dispatch.charge_power) * battery.dispatch.round_trip_efficiency / 100.0
             == pytest.approx(sum(battery.dispatch.discharge_power)))
+
+
+def test_hybrid_dispatch_baseload_heuristic_and_analysis(site):
+
+    desired_schedule = 8760*[20]
+
+    desired_schedule_site = SiteInfo(flatirons_site,
+                                     desired_schedule=desired_schedule)
+    wind_solar_battery = {key: technologies[key] for key in ('pv', 'wind', 'battery')}
+
+    dispatch_options = {'battery_dispatch': 'load_following_heuristic',
+                        'use_higher_hours': True, 
+                        'higher_hours': {'min_regulation_hours': 4, 'min_regulation_power': 5000}}
+
+    hopp_config = {
+        "site": desired_schedule_site,
+        "technologies": wind_solar_battery,
+        "config": {
+            "dispatch_options": dispatch_options
+        }
+    }
+    hopp_config["technologies"]["grid"] = {
+        "interconnect_kw": interconnect_mw * 1000
+    }
+    hi = HoppInterface(hopp_config)
+    hi.simulate(1)
+
+    hybrid_plant = hi.system
+
+    assert hybrid_plant.grid.time_load_met == pytest.approx(91.9, 1e-2)
+    assert hybrid_plant.grid.capacity_factor_load == pytest.approx(94.45, 1e-2)
+    assert hybrid_plant.grid.total_number_hours == pytest.approx(3732, 1e-2)
\ No newline at end of file