Issue135 forecast uncertainty 2

README.md

@@ -63,7 +63,7 @@ API requests that interact with a running test case (those that require a ``test
 | Receive boundary condition forecast from current communication step for the given point names for the horizon and at the interval in seconds.   |  PUT ``forecast/{testid}`` with required arguments ``point_names=<list of strings>``, ``horizon=<value>``, ``interval=<value>``|
 | Receive boundary condition forecast available point names and metadata. |  GET ``forecast_points/{testid}``                              |
 | Receive current test scenario.                                         |  GET ``scenario/{testid}``                                   |
-| Set test scenario. Setting the argument ``time_period`` performs an initialization with predefined start time and warmup period and will only simulate for predefined duration. |  PUT ``scenario/{testid}`` with optional arguments ``electricity_price=<string>``, ``time_period=<string>``.  See the [Test Case](https://ibpsa.github.io/project1-boptest/testcases/index.html) page for options and documentation.|
+| Set test scenario. Setting the argument ``time_period`` performs an initialization with predefined start time and warmup period and will only simulate for predefined duration. |  PUT ``scenario/{testid}`` with optional arguments ``electricity_price=<string>``, ``time_period=<string>``, ``temperature_uncertainty=<string>``, ``solar_uncertainty=<string>``, ``seed=<int>``.  See the [Test Case](https://ibpsa.github.io/project1-boptest/testcases/index.html) page for options and documentation.|
 | Get test status as `Running` or `Queued`                                                                                    | GET ``status/{testid}``                                    |
 | Stop a queued or running test.  Needed to deploy a new test case when no more idle workers are avaiable.                                                                                              | PUT ``stop/{testid}``                                      |
 | Receive BOPTEST version.                                               |  GET ``version/{testid}``                                             |

contributors.md

@@ -30,6 +30,9 @@ Thank you to all who have provided guidance on the development of this software.
 - Christian Veje, University of Southern Denmark
 - Draguna Vrabie, Pacific Northwest National Laboratory
 - Harald Walnum, SINTEF
+- Zhe Wang, The Hong Kong University of Science and Technology
 - Michael Wetter, Lawrence Berkeley National Laboratory
 - Tao Yang, University of Southern Denmark
+- Laura Zabala, R2M Solution Spain
 - Ettore Zanetti, Lawrence Berkeley National Laboratory
+- Wanfu Zheng, The Hong Kong University of Science and Technology

forecast/error_emulator.py

@@ -0,0 +1,103 @@
+'''
+Created on Sept 5, 2022
+
+@author: Laura Zabala, Wanfu Zheng
+
+'''
+
+import numpy as np
+
+
+def mean_filter(data, window_size=3):
+    """
+    Apply mean filter to a 1D list of data.
+
+    Parameters
+    ----------
+    data: list of numbers
+        List of numbers to be filtered.
+    window_size: int
+        Size of the filtering window. Must be an odd number.
+
+    Returns
+    -------
+    filtered_data : list of numbers
+        List of filtered data.
+
+    """
+
+    if window_size % 2 == 0:
+        raise ValueError("Window size must be an odd number.")
+
+    half_window = window_size // 2
+    filtered_data = data.copy()
+
+    for i in range(half_window, len(data) - half_window):
+        if data[i] == 0:
+            continue
+
+        window_data = data[i - half_window : i + half_window + 1]
+        filtered_data[i] = sum(window_data) / len(window_data)
+
+    return filtered_data
+
+def predict_temperature_error_AR1(hp, F0, K0, F, K, mu):
+    '''
+    Generates an error for the temperature forecast with an AR model with normal distribution in the hp points of the predictions horizon.
+
+    Parameters
+    ----------
+    hp : int
+        Number of points in the prediction horizon.
+    F0 : float
+        Mean of the initial error model.
+    K0 : float
+        Standard deviation of the initial error model.
+    F : float
+        Autocorrelation factor of the AR error model, value should be between 0 and 1.
+    K : float
+        Standard deviation of the AR error model.
+    mu : float
+        Mean value of the distribution function integrated in the AR error model.
+
+
+    Returns
+    -------
+    error : 1D array
+        Array containing the error values in the hp points.
+
+    '''
+
+    error = np.zeros(hp)
+    error[0] = np.random.normal(F0, K0)
+    for i_c in range(hp - 1):
+        error[i_c + 1] = np.random.normal(
+            error[i_c] * F + mu, K
+        )
+
+    return error
+
+def predict_solar_error_AR1(hp, ag0, bg0, phi, ag, bg):
+    '''
+    Generates an error for the solar forecast based on the specified parameters using an AR model with Laplace distribution in the hp points of the predictions horizon.
+
+    Parameters
+    ----------
+    hp : int
+        Number of points in the prediction horizon.
+    ag0, bg0, phi, ag, bg : float
+        Parameters for the AR1 model.
+
+    Returns
+    -------
+    error : 1D numpy array
+        Contains the error values in the hp points.
+
+    '''
+
+    error = np.zeros(hp)
+    error[0] = np.random.laplace(ag0, bg0)
+    for i_c in range(1, hp):
+        error[i_c] = np.random.laplace(error[i_c - 1] * phi + ag, bg)
+
+    return error

forecast/forecast_uncertainty_params.json

@@ -0,0 +1,49 @@
+{
+    "temperature": {
+        "low": {
+            "F0": 0,
+            "K0": 0.6,
+            "F": 0.92,
+            "K": 0.4,
+            "mu": 0
+        },
+        "medium": {
+            "F0": 0.15,
+            "K0": 1.2,
+            "F": 0.93,
+            "K": 0.6,
+            "mu": 0
+        },
+        "high": {
+            "F0": -0.58,
+            "K0": 1.5,
+            "F": 0.95,
+            "K": 0.7,
+            "mu": -0.015
+        }
+    },
+    "solar": {
+        "low": {
+            "ag0": 4.44,
+            "bg0": 57.42,
+            "phi": 0.62,
+            "ag": 1.86,
+            "bg": 45.64
+        },
+        "medium": {
+            "ag0": 15.02,
+            "bg0": 122.6,
+            "phi": 0.63,
+            "ag": 4.44,
+            "bg": 91.97
+        },
+        "high": {
+            "ag0": 32.09,
+            "bg0": 119.94,
+            "phi": 0.67,
+            "ag": 10.63,
+            "bg": 87.44
+        }
+    }
+
+}

forecast/forecaster.py

@@ -1,21 +1,26 @@
 '''
 Created on Apr 25, 2019
 
-@author: Javier Arroyo
+@author: Javier Arroyo, Laura Zabala, and Wanfu Zheng
 
 This module contains the Forecaster class with methods to obtain
 forecast data for the test case. It relies on the data_manager object
-of the test case to provide deterministic forecast.
+of the test case to provide a deterministic forecast, and the
+error_emulator module to generate errors for uncertain forecasts.
 
 '''
+from .error_emulator import predict_temperature_error_AR1, predict_solar_error_AR1, mean_filter
+import numpy as np
+import json
+
 
 class Forecaster(object):
-    '''This class retrieves test case data forecast for its use in
-    optimal control strategies.
+    '''
+    This class retrieves test case data forecast for its use in optimal control strategies.
 
     '''
 
-    def __init__(self, testcase):
+    def __init__(self, testcase, forecast_uncertainty_params_path='forecast/forecast_uncertainty_params.json'):
         '''
         Constructor
 
@@ -24,51 +29,138 @@ def __init__(self, testcase):
         testcase: BOPTEST TestCase object
             object of an already deployed test case that
             contains the data stored from the test case run
+        forecast_uncertainty_params_path : str, optional
+            Path to the JSON file containing the uncertainty parameters.
+            Default is 'forecast/forecast_uncertainty_params.json'.
 
         '''
 
         # Point to the test case object
         self.case = testcase
+        # Load forecast uncertainty parameters
+        self.uncertainty_params = self.load_uncertainty_params(forecast_uncertainty_params_path)
 
-    def get_forecast(self,point_names, horizon=24*3600, interval=3600,
-                     category=None, plot=False):
-        '''Returns forecast of the test case data
+    def get_forecast(self, point_names, horizon=24 * 3600, interval=3600,
+                     wea_tem_dry_bul=None, wea_sol_glo_hor=None, seed=None):
+        '''
+        Retrieves forecast data for specified points over a given horizon and interval.
 
         Parameters
         ----------
         point_names : list of str
-            List of forecast point names for which to get data.
-        horizon : int, default is 86400 (one day)
-            Length of the requested forecast in seconds. If None,
-            the test case horizon will be used instead.
-        interval : int, default is 3600 (one hour)
-            resampling time interval in seconds. If None,
-            the test case interval will be used instead.
-        category : string, default is None
-            Type of data to retrieve from the test case.
-            If None it will return all available test case
-            data without filtering it by any category.
-            Possible options are 'weather', 'prices',
-            'emissions', 'occupancy', internalGains, 'setpoints'
-        plot : boolean, default is False
-            True if desired to plot the forecast
+            List of data point names for which the forecast is to be retrieved.
+        horizon : int, optional
+            Forecast horizon in seconds (default is 86400 seconds, i.e., one day).
+        interval : int, optional
+            Time interval between forecast points in seconds (default is 3600 seconds, i.e., one hour).
+        wea_tem_dry_bul : str, optional
+            Uncertainty level for outside air dry bulb temperature.  'low', 'medium', or 'high'
+            If None, defaults to no forecast error.
+            Default is None.
+        wea_sol_glo_hor : dict, optional
+            Uncertainty level for outside solar radiation.  'low', 'medium', or 'high'
+            If None, defaults to no forecast error.
+            Default is None.
+        seed : int, optional
+            Seed for the random number generator to ensure reproducibility of the stochastic forecast error.
+            If None, no seed is used.
+            Default is None.
 
         Returns
         -------
         forecast : dict
-            Dictionary with the requested forecast data
-            {<variable_name>:<variable_forecast_trajectory>}
-            where <variable_name> is a string with the variable
-            key and <variable_forecast_trajectory> is a list with
-            the forecasted values. 'time' is included as a variable
+            A dictionary containing the forecast data for the requested points with applied error models.
 
         '''
 
+        # Set uncertainty parameters to 0 if no forecast uncertainty
+        temperature_params = {
+            "F0": 0, "K0": 0, "F": 0, "K": 0, "mu": 0
+        }
+
+        solar_params = {
+            "ag0": 0, "bg0": 0, "phi": 0, "ag": 0, "bg": 0
+        }
+
+        if wea_tem_dry_bul is not None:
+            temperature_params.update(self.uncertainty_params['temperature'][wea_tem_dry_bul])
+
+        if wea_sol_glo_hor is not None:
+            solar_params.update(self.uncertainty_params['solar'][wea_sol_glo_hor])
+
         # Get the forecast
         forecast = self.case.data_manager.get_data(variables=point_names,
                                                    horizon=horizon,
-                                                   interval=interval,
-                                                   category=category,
-                                                   plot=plot)
+                                                   interval=interval)
+
+        # Add any outside dry bulb temperature error
+        if 'TDryBul' in point_names and any(temperature_params.values()):
+            if seed is not None:
+                np.random.seed(seed)
+            # error in the forecast
+            error_forecast_temp = predict_temperature_error_AR1(
+                hp=int(horizon / interval + 1),
+                F0=temperature_params["F0"],
+                K0=temperature_params["K0"],
+                F=temperature_params["F"],
+                K=temperature_params["K"],
+                mu=temperature_params["mu"]
+            )
+
+            # forecast error just added to dry bulb temperature
+            forecast['TDryBul'] = forecast['TDryBul'] - error_forecast_temp
+            forecast['TDryBul'] = forecast['TDryBul'].tolist()
+
+        # Add any outside global horizontal irradiation error
+        if 'HGloHor' in point_names and any(solar_params.values()):
+
+            original_HGloHor = np.array(forecast['HGloHor']).copy()
+            lower_bound = 0.2 * original_HGloHor
+            upper_bound = 2 * original_HGloHor
+            indices = np.where(original_HGloHor > 50)[0]
+
+            for i in range(200):
+                if seed is not None:
+                    np.random.seed(seed+i*i)
+                error_forecast_solar = predict_solar_error_AR1(
+                    int(horizon / interval + 1),
+                    solar_params["ag0"],
+                    solar_params["bg0"],
+                    solar_params["phi"],
+                    solar_params["ag"],
+                    solar_params["bg"]
+                )
+
+                forecast['HGloHor'] = original_HGloHor - error_forecast_solar
+
+                # Check if any point in forecast['HGloHor'] is out of the specified range
+                condition = np.any((forecast['HGloHor'][indices] > 2 * original_HGloHor[indices]) |
+                                   (forecast['HGloHor'][indices] < 0.2 * original_HGloHor[indices]))
+                # forecast['HGloHor']=gaussian_filter_ignoring_nans(forecast['HGloHor'])
+                forecast['HGloHor'] = mean_filter(forecast['HGloHor'])
+                forecast['HGloHor'] = np.clip(forecast['HGloHor'], lower_bound, upper_bound)
+                forecast['HGloHor'] = forecast['HGloHor'].tolist()
+                if not condition:
+                    break
 
         return forecast
+
+    def load_uncertainty_params(self, filepath):
+        '''Load the uncertainty parameters from a JSON file.
+
+        Parameters
+        ----------
+        filepath : str
+            Path to the JSON file containing the uncertainty parameters.
+
+        Returns
+        -------
+        dict
+            Uncertainty parameters loaded from the JSON file.
+
+        '''
+
+        with open(filepath, 'r') as f:
+            uncertainty_params = json.load(f)
+
+        return uncertainty_params

releasenotes.md

@@ -8,6 +8,7 @@ Released on xx/xx/xxxx.
 
 - Add note to ``README.md`` about using environment variable ``BOPTEST_TIMEOUT`` to edit the timeout period for idle workers.  This is for [#715](https://github.com/ibpsa/project1-boptest/issues/715).
 - Add note to ``README.md`` about a Julia interface implemented by [BOPTestAPI.jl](https://terion-io.github.io/BOPTestAPI.jl/stable/).  This is for [#707](https://github.com/ibpsa/project1-boptest/issues/707).
+- Add weather forecast uncertainty as new scenario options for dry bulb temperature and global horizontal irradiation.  The corresponding new scenario keys are ``temperature_uncertainty`` and ``solar_uncertainty``, which can take values ``low``, ``medium``, or ``high``.  A new scenario key ``seed`` is also added to set an integer seed for reproducible uncertainty generation.  This is for [#135](https://github.com/ibpsa/project1-boptest/issues/135).
 
 
 ## BOPTEST v0.7.0