Issue135 forecast uncertainty

.gitignore

@@ -15,3 +15,7 @@ testcase2/doc/build/
 parser/*.fmu
 parser/*.mo
 parser/*.json
+
+xx.ipynb
+connect_boptest.py
+\!

forecast/error_emulator.py

@@ -0,0 +1,95 @@
+'''
+Created on Sept 5, 2022
+
+@author: Laura Zabala
+
+'''
+
+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 to be filtered.
+    - window_size: Size of the filtering window. Must be an odd number.
+
+    Returns:
+    - 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,13 +1,16 @@
 '''
 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.
 
 '''
+from .error_emulator import predict_temperature_error_AR1, predict_solar_error_AR1, mean_filter
+import numpy as np
+
 
 class Forecaster(object):
     '''This class retrieves test case data forecast for its use in
@@ -30,45 +33,146 @@ def __init__(self, testcase):
         # Point to the test case object
         self.case = testcase
 
-    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 : dict, optional
+            Parameters for the AR1 model to simulate forecast error in dry bulb temperature:
+                - F0, K0, F, K, mu : parameters used in the AR1 model.
+            If None, defaults to a dictionary with all parameters set to zero, simulating no forecast error.
+        wea_sol_glo_hor : dict, optional
+            Parameters for the AR1 model to simulate forecast error in global horizontal solar irradiation:
+                - ag0, bg0, phi, ag, bg : parameters used in the AR1 model.
+            If None, defaults to a dictionary with all parameters set to zero, simulating no forecast error.
+        seed : int, optional
+            Seed for the random number generator to ensure reproducibility of the stochastic forecast error.
 
         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.
+=======
+                     weather_temperature_dry_bulb=None, weather_solar_global_horizontal=None, seed=None,
+                     category=None, plot=False):
 
-        '''
+        if weather_temperature_dry_bulb is None:
+            weather_temperature_dry_bulb = {
+                "F0": 0, "K0": 0, "F": 0, "K": 0, "mu": 0
+            }
 
+        if weather_solar_global_horizontal is None:
+            weather_solar_global_horizontal = {
+                "ag0": 0, "bg0": 0, "phi": 0, "ag": 0, "bg": 0
+            }
         # Get the forecast
         forecast = self.case.data_manager.get_data(variables=point_names,
                                                    horizon=horizon,
-                                                   interval=interval,
-                                                   category=category,
-                                                   plot=plot)
+                                                   interval=interval)
+
+        if 'TDryBul' in point_names and any(wea_tem_dry_bul.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=wea_tem_dry_bul["F0"],
+                K0=wea_tem_dry_bul["K0"],
+                F=wea_tem_dry_bul["F"],
+                K=wea_tem_dry_bul["K"],
+                mu=wea_tem_dry_bul["mu"]
+            )
+
+            # forecast error just added to dry bulb temperature
+            forecast['TDryBul'] = forecast['TDryBul'] - error_forecast_temp
+            forecast['TDryBul'] = forecast['TDryBul'].tolist()
+        if 'HGloHor' in point_names and any(wea_sol_glo_hor.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),
+                    wea_sol_glo_hor["ag0"],
+                    wea_sol_glo_hor["bg0"],
+                    wea_sol_glo_hor["phi"],
+                    wea_sol_glo_hor["ag"],
+                    wea_sol_glo_hor["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
+
+        if 'TDryBul' in point_names and any(weather_temperature_dry_bulb.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=weather_temperature_dry_bulb["F0"],
+                K0=weather_temperature_dry_bulb["K0"],
+                F=weather_temperature_dry_bulb["F"],
+                K=weather_temperature_dry_bulb["K"],
+                mu=weather_temperature_dry_bulb["mu"]
+            )
+
+            # forecast error just added to dry bulb temperature
+            forecast['TDryBul'] = forecast['TDryBul'] - error_forecast_temp
+            forecast['TDryBul'] = forecast['TDryBul'].tolist()
+        if 'HGloHor' in point_names and any(weather_solar_global_horizontal.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),
+                    weather_solar_global_horizontal["ag0"],
+                    weather_solar_global_horizontal["bg0"],
+                    weather_solar_global_horizontal["phi"],
+                    weather_solar_global_horizontal["ag"],
+                    weather_solar_global_horizontal["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

restapi.py

@@ -82,12 +82,23 @@ def handle_validation_error(self, error, bundle_errors):
 parser_scenario = reqparse.RequestParser(argument_class=CustomArgument)
 parser_scenario.add_argument('electricity_price', type=str)
 parser_scenario.add_argument('time_period', type=str)
+parser_scenario.add_argument('temperature_uncertainty', type=str)
+parser_scenario.add_argument('solar_uncertainty', type=str)
+parser_scenario.add_argument('seed', type=int)
 # ``forecast`` interface
 parser_forecast_points = reqparse.RequestParser(argument_class=CustomArgument)
 parser_forecast_points.add_argument('point_names', type=list, action='append', required=True)
-forecast_parameters = ['horizon', 'interval']
-for arg in forecast_parameters:
-    parser_forecast_points.add_argument(arg, required=True)
+
+forecast_parameters = ['horizon', 'interval', 'temperature_uncertainty', 'solar_uncertainty']
+
+# Add required parameters
+parser_forecast_points.add_argument('horizon', required=True)
+parser_forecast_points.add_argument('interval', required=True)
+
+# Add optional uncertainty parameters
+parser_forecast_points.add_argument('temperature_uncertainty', required=False)
+parser_forecast_points.add_argument('solar_uncertainty', required=False)
+
 # ``results`` interface
 results_var = reqparse.RequestParser(argument_class=CustomArgument)
 results_var.add_argument('point_names', type=list, action='append', required=True)
@@ -200,10 +211,19 @@ def put(self):
         args = parser_forecast_points.parse_args()
         horizon = args['horizon']
         interval = args['interval']
+
+        # Extract optional uncertainty parameters if provided
+        temperature_uncertainty = args.get('temperature_uncertainty', None)
+        solar_uncertainty = args.get('solar_uncertainty', None)
+
         point_names = []
         for point_name in args['point_names']:
             point_names.append(''.join(point_name))
-        status, message, payload = case.get_forecast(point_names, horizon, interval)
+
+        # Modify the get_forecast function call to handle the new parameters
+        status, message, payload = case.get_forecast(point_names, horizon, interval, temperature_uncertainty,
+                                                     solar_uncertainty)
+
         return construct(status, message, payload)
 
 class Scenario(Resource):
@@ -276,3 +296,4 @@ def post(self):
 
 if __name__ == '__main__':
     app.run(host='0.0.0.0')
+