# Deep Learning-Based Multivariate Probabilistic Forecasting for Short-Term Scheduling in Power Markets

By Caspar Boersma
4445074
[e-mail]

Jan van Dorth tot Medler
4453344
[e-mail]

Click the button below to view the project code in Google Colab
[](https://colab.research.google.com/drive/1NqxkzaL0R6udsjrzEDHg3QW97QYa1AaY?usp=sharing)

## Introduction
#### Why is this paper relevant, problem description
As renewable power sources has gotten more popular in the past years, the operation of power networks became more complex. This is caused by the fluctuating nature of these renewable sources like solar and wind power. Compared to the steady generation of power with coal power plants and gas turbines, wind speed and cloud cover are highly variable, causing the power generation to vary too.

To succesfully manage the power networks, a bidding system is introduced to the power market where electricity generators and consumers can submit how much power they predict they will generate or consume the following day and for which price. Thus, to make the most profit, generators and consumers try to predict the day ahead energy price to make the most profitable biddings for the day head.

The authors of "Deep Learning-Based Multivariate Probabilistic Forecasting for Short-Term Scheduling in Power Markets"[1] propose to use a bidirectional LSTM network to predict the future price of energy. In this blogpost, this bidirectional LSTM network will be build from scratch using the parameters and the model architecture proposed in the paper. The objective here is to replicate figure 6d from the paper which visualizes the hourly day ahead price prediction for a week, as shown below:

###### *Visualization Original Results. Source: [1]*

As can be seen in this plot, the predictions do not consist of a single price per hour, but instead of a range of values, formed by different quantiles. These quantiles give an estimation of the uncertainty of the prediction, so when the uncertainty is high, it can be decided to use a different final bidding strategy.

By doing this reproduction of the before mentioned paper, we can test if the method proposed by the author is executable and if it indeed leads to the claimed results and conclusions.

## Method

### Dataset
To be able to train a model that predicts energy prices, historical data is needed. Renewable energy sources such as wind and solar energy are highly dependend on the weather. Therefore, hourly historical weather data is used in the input of the model. The weather properties taken into account are windspeed and cloud cover, arguably the most important for predicting wind and solar energy generation. The hourly historical energy price data is also used both as input and output (explained further in the section 'data processing').

In contrast to the paper, we will use historical weather data instead of weather predictions for the price prediction. This is because we were not able to obtain this weather prediction data.

The location for predicting the price data used by [1] is Mons, Belgium. Accordingly the weather and price data for this location is collected. The data used in this reproduction, is from a five year period ranging from january 5th 2015 to february 4th 2020, as the original paper uses a training set consisting of five years as well. In line with the paper we use a single month as test set, for which we chose the last month of the data set; a winter month like is used in the paper.

For the historical weather data, the World Weather Online Local History API [2] is used, which provides data every three hours (which is converted to hourly data by using linear interpolation). A script is written to automatically retrieve the weather data for the specific time frame that is needed. The API call returns JSON data (as seen below), which is converted to a more usable pandas dataframe in the code.

The weather data includes many different metrics. To keep things simple, only two arguably important metrics are chosen, namely windspeed and cloudcover, respectively for wind and solar energy generation.

"hourly": [
            {
            "time": "0",
            "tempC": "1",
            "tempF": "33",
            "windspeedMiles": "7",
            "windspeedKmph": "11",
            "winddirDegree": "301",
            "winddir16Point": "WNW",
            "weatherCode": "248",
            "precipMM": "0.0",
            "precipInches": "0.0",
            "humidity": "97",
            "visibility": "0",
            "visibilityMiles": "0",
            "pressure": "1028",
            "pressureInches": "31",
            "cloudcover": "100",
            "HeatIndexC": "1",
            "HeatIndexF": "33",
            "DewPointC": "0",
            "DewPointF": "32",
            "WindChillC": "-3",
            "WindChillF": "27",
            "WindGustMiles": "15",
            "WindGustKmph": "23",
            "FeelsLikeC": "-3",
            "FeelsLikeF": "27",
            "uvIndex": "1"
            }

###### *JSON data from weather API*

In the paper [1], they also use wind and solar generation as input data. To simplify things, we only use weather and price data. This is done under the assumption that the wind and solar generation are arguably highly dependent on certain weather metrics, namely cloud cover and windspeed. This would mean that wind and solar generation data would be a largely rendundant addition.

The historical hourly day ahead price data is downloaded from the Entsoe Transparency Platform[3] by the year, combined into a .csv file. In this dataset, there are 4 hours where no price is available. To be able to combine this price data properly with the corresponding weather data, these missing prices have to be filled up. This is done by taking the mean value of the prices an hour before and after the respective missing value. Then this final .csv file can be loaded into Python using Pandas, resulting in the following dataframe:

###### *Price data in pandas*

Now that the weather data and price data are both hourly, they can be comined. This is shown in the dataframe below.

###### *Weather and price data in pandas*

Finally the data is normalized, which is only necessary to perform on the data that is not one-hot-encoded. This is done by selecting the following subsection of the dataframe, performing the scaling operation on it with an SKLearn standardscaler, and replacing these colums in the original dataset (including the one-hot-encoded variables):

###### *Scaled data*

### Data processing
Training an accurate model requires using the right training data. In [1] it is mentioned that the BLSTM performs better, when it also has information in the sequence about the same day in the previous week. That is why each input sequence in our reproduction consists of a whole weak of training data, resulting in a sequence of `7*24 = 168` timesteps. The goal of this reproduction is then to use this input sequence to predict the 24 hours of the prediction horizon of interest:

###### *Prediction horizon of interest, source: [1]*

Note that the input sequence consists of the 7 days before 12h00. This is the moment every actor in the market has to provide their trading strategy for the following day (24 hours). Note that there is a gap of 12 hours between the last known datapoints and the start of the predicted sequence.

In order to create a training and test set, to take all data into account in the specific format from the previous paragraph, a 'sliding window' creates the samples. A single sample consists of the features of the past seven days as input (X), and the day ahead prices of a whole day as corresponding output (y). The 'sliding window' slides over the whole dataset, with steps of a single hour, as visualized in the figure below. This sliding window is used to make full use of the relatively small dataset.

Here, the full black dots reprent the price data of a whole day (consisting of 24 hours of data each) and the hollow dots represent the weather data. The one-hot-encoded days of the weak, months and days of month are not shown for visualisation purposes, as the total number of features increases to 77 when they are included. Every sample has a different color and consists of large brackets and small brackets, representing the input and corresponding output data of a single sample. Note that there is half a dot (12 hours) between the end of the input and the start of the corresponding output, which is in line with the prediction horizon figure.

###### *Sliding window over the dataset*

During the sliding, every sample is added to the new dataset in the batch dimension as seen in the figure below. This new dataset will be split into a train and test set.

###### *New dataset with batch dimension*

This operation is performed with the following code:

n_future = 24   # Number of hours we want to predict
n_inter = 12  # Number of hours between the prediction area of interest and the last known data
n_past = 168   # Number of past hours we want to use to predict the future (7 days)

for i in range(n_past, len(df_small) - n_future - n_inter + 1): 
    X.append(df_small[i - n_past:i, 0:df_small.shape[1]])
    df_inter = df_small[i + n_inter: i + n_inter + n_future, df.columns.get_loc("price") - 1]
    df_inter =df_inter.reshape(n_future, 1)
    y.append(df_inter)

The energy price is, next to the weather, also influenced by periodicity. At night, the price of energy is usually significantly lower than during daytime. This periodicity needs to be represented in the features and this is done by adding the hours, weekdays (e.g. monday), months and monthdays (e.g. 5th day of the month) to the data with one-hot encoding.

Then, the data is split into a training and a test set. The training set ranging from january 5th 2015 to januari 4th 2020 and the test set ranging from januari 5th 2020 to february 4th 2020, which is a single winter month, as was used in the original paper. For the test set we iterate back with steps of 24 hours from the end of the dataset. This is done because our test set is the last month of the dataset and we do not want overlapping data in this month, because this can't be plotted in the same way as is done in the figure that will be reproduced.

### Neural network
The neural network model used is an Bidirectional Long Short Term Memory network (BLSTM). A BLSTM is a form of a Bidirectional Recurrent Neural Network (BRRN). RNN's are commonly used for learning patterns in sequences and are useful for tasks such as language translation and time series predictions.
RNN's have the problem of vanishing gradients, which causes long distance relations to be difficult to learn. LSTM's solve this with its gates and 'conveyer belt' memory (further details below) that allow information to propagate through the network unchanged or little changed, resulting in the possibility of uninterrupted gradient flow.

A LSTM network consists of a sequence of LSTM cells. To explain the LSTM network it is wise to first examine these individual cells. They key to the operations of these cells is that they make use of several diferent gates.
These gates are:
1. Forget gate F: what previous hidden state information to use/discard
2. Input gate I: what input information to use/discard
3. Output gate O: what to output

Further more, a LSTM has a cell state memory C which can be thought of as a 'conveyor belt' [5]. The gates described above can influence what information of the cell state is let through and what is added to it.

The gates work with sigmoid functions, which have outputs between 0 and 1.
1. The forget gate looks at the previous hidden state Ht-1 and the input Xt. It applying the weights and outputs a number between 0 and 1. This number is then multiplied with the previous memory state C_t-1. A 0 value means to 'forget' the previous cell state and a 1 means to fully keep it.
2. The input gate looks at H_t-1sub> and Xt, applies weights and also outputs a number between 0 and 1. This number is then multiplied with the candidate input to decide how much of the input is added to the cell state. A tanh function is applied to the canditate input to prevent an exploding gradient (which would make backprop more difficult).
3. The output gate again looks at H_t-1 and X_t, applies weights and again outputs a number between 0 and 1. It is multiplied with the new cell state (which resulted from adding the input in the previous step), after first applying a tanh layer to the cell state, again for preventing exploding gradients. Here is decided how much of the new cell state is outputted to the new hidden layer.

All gates use learnable weights, meaning these can be finetuned by backpropagationg while training the model.
A visual representation is provided in the figure below.

###### *Source: Jan van Gemert Slides Lecture 6 CS4240 2020/2021*

Below, the formulae used for calculating outputs in a LSTM cell are described.

###### *Source: Jan van Gemert Slides Lecture 6 CS4240 2020/2021*

When multiple of these LSTM cells are chained in a sequence as visualized in the figure below, they become a LSTM network. Here , the outputs C_t-1 and H_t-1 of the previous cell are forwarded to the current cell. The current cell uses X_t, C_t-1 and H_t-1 to calculate C_t and H_t and output them to the next cell. This sequence continues for every cell.

###### *Sequence of multiple LSTM cells forming a LSTM network. Source: [5]*

In comparison to a LSTM, a BLSTM has two cells per time step, as seen in the figure below. Like a LSTM, one (the A cell) is for propagating information forward in time. But unlike a LSTM, the second cell (the A' cell) is for propagating information backwards. Both A and A' are connected to the same output.

This network architecture makes it possible to -next to forward relations- also learn backward relations. A general example with speech-recognition is that it is easier to predict a word when instead of only the past words, also the future words are taken into account. This creates more context for the word that needs to be predicted, and thus is easier predicted right. The addition of the backward layer, that introduces relations from future words to the present word, makes this possible.

In the case of energy price prediction, this leads to the network not only forming relation from the past days to the future days, but also the other way around. In other words, wereas a LSTM would for example only learn the influence of conditions 7 days ago on conditions 4 days ago, a BLSTM also learns the influence of conditions 1 day ago on conditions 4 days ago. Thus, a BLSTM should generate increased performance here.

###### *BLSTM network. Source: https://towardsdatascience.com/understanding-bidirectional-rnn-in-pytorch-5bd25a5dd66*

Piling up three BLSTM layers, results in higher level representation, making it a deeper network. The outputs of the first layer are the inputs of the second layer, etc. This is done in the code block below.

inputs = Input(shape=(X_train.shape[1], X_train.shape[2]))
lstm = Bidirectional(LSTM(15, return_sequences=True, dropout=0.5))(inputs, training = True)
lstm = Bidirectional(LSTM(15, return_sequences=True, dropout=0.5))(lstm, training = True)
lstm = Bidirectional(LSTM(15, return_sequences=False, dropout=0.5))(lstm, training = True)

In the paper [1], it is argued that the best results are obtained with a 3 layer BLSTM with around 10 to 20 neurons each. Therefore, we implemented a 3 layer BLSTM with 15 neurons each for the most accurate reproduction.

### Probabilistic predictions with Quantiles
To get a sense of how certain the model is that it predicted the right input, using point forcasts won't suffice. In this case, point forcasting would provide one predicted price value per hour, resulting in no way to inspect the certainty of the model for that hour. Probabilistic predictions, however, have the power to show a probability distribution of the price at a certain hour. If the distribution is large, this would suggest low certainty and vice versa. This capacity to know when to trust the model is why probabilistic predictions are preferred. There exist two ways to obtain the distribution, namely with a parametrized model or with an empirical funtions. For the parametrized model, the type of distribution needs to be assumed. Because that asumption is difficult to get right, emperical functions are used here.

###### *Quantile regression. Source: [1]*

Because quantiles are used here, the loss function is chosen as the pinball loss. It describes the accuracy of a quantile forecasting model.

###### *Total pinball loss. Source: [1]*

def q_loss(q,y,f):
    e = (y-f)
    return K.mean(K.maximum(q*e, (q-1)*e), axis=-1)

losses = [lambda y,f: q_loss(0.1,y,f), lambda y,f: q_loss(0.5,y,f), 
          lambda y,f: q_loss(0.9,y,f), lambda y,f: q_loss(0.01,y,f),
          lambda y,f: q_loss(0.99,y,f), lambda y,f: q_loss(0.05,y,f),
          lambda y,f: q_loss(0.95,y,f), lambda y,f: q_loss(0.25,y,f),
          lambda y,f: q_loss(0.75,y,f)]

Nine quantiles, namely 0.01, 0.05, 0.1, 0.25, 0.5, 0.75, 0.90, 0.95 and 0.99 are chosen. For each quantile, there is a separate dense layer, outputting the value of each quantile for each hour. After plotting the values of the quantiles, the distribution becomes visible.

For robust converging, Adam is used as optimizer. For regularization, dropout with p=0.5 and early stopping with a patience of 3 is used. The calculations were performed with a Google Colab GPU and training the model takes around 45 minutes with 70 epochs.

### Metrics

The pinball loss will be used as metric while training, because it results in a division into optimal quantiles.

The final performance of the model in the paper is given by means of two metrics. First the Root Mean Squared Error (RMSE) between the predicted median (50^th quantile) and the actual hourly price is determined in table I of the paper. Subsequently the hourly total quantile error of the nine quantiles is used to compare the model to one using a Gaussian to estimate the uncertainty in table III of the paper.

## Results
#### Graphical

In the original paper, the predictions of several days, concatenated with eachother, were plotted to visualize the results. In the following figure these predictions and the distances between certain quantiles can be seen.
###### *Visualization Original Results. Source: [1]*

In this blog post, tried to reproduce this figure from the original paper specifically. To clarify the value we will use to calculate the mean squared error with in the following subsection, we specifically showed the value of the median, which is the 50^th quantile, by means of a black line. The results are visualized below:

###### *Visualization of Our Results*

It can be seen that our model is much more certain about it's predictions, based on historical data, than the model from the paper. This can be seen by the lower spread between the first and 99^th quantile in our plot, compared to that of the original paper. Furthermore it can be seen that our model generates prices that are often much higher than the actual value.

#### Numerical

We calculated the mean squared error (mse) with the predicted median. This resulted in a mean squared error of €141.27 for the whole test month. The result of the original paper (table III in [1]) was a €28 for their BLSTM using quantile regression.

Using the total hourly quantile loss of the entire test month, we obtained a value of €608.21. This is a lot higher than the original paper, with a value of €28.

## Conclusion
In this reproduction we have used a sequence of historical weather and electricity price. We have managed to reproduce results using the proposed architecture from the paper [1]. However, our results of the metrics are a lot worse compared to the original paper. In the discussion we give multiple reasons why this could be the case.

## Discussion
For the discrepancies between our results and the results in the paper there might be several explainations. These can also be seen as recommendations for future reproduction.

1. As described earlier in this blog, we didn't use the PV and wind load data for price prediction. This data could have more influence on the predictions than we earlier assumed.
2. The low spread between different quantiles might have something to do with the fact that we used a sliding window when generating our training data, which can give the model a false sense of accuracy on a limited dataset like we used.
3. In this reproduction we used historical weather data instead of weather predictions as done in the paper. User weather predictions as input would, without doubt, result in substantially more accurate predictions than when using past weather data.
4. We did not tune any hyperparameters, but have used the ones from the paper as much as possible. However, the paper is often very vague in the parameters they used. Further hyperparameter tuning could improve our results as well.

## Bibliography
[1] Toubeau, J. F., Bottieau, J., Vallée, F., & De Grève, Z. (2018). Deep learning-based multivariate probabilistic forecasting for short-term scheduling in power markets. IEEE Transactions on Power Systems, 34(2), 1203-1215.
[2] https://www.worldweatheronline.com/developer/api/historical-weather-api.aspx
[3] https://transparency.entsoe.eu/transmission-domain/r2/dayAheadPrices/show
[4] https://medium.com/@raghavaggarwal0089/bi-lstm-bc3d68da8bd0
[5] https://colah.github.io/posts/2015-08-Understanding-LSTMs/

## Overview of who did what
Caspar Boersma:
- Neural network design
- Price data processing
- Predicting
- Visualization

Jan van Dorth tot Medler:
- Data collection
- Weather data processing
- General data processing
- Quantile regression