**Author Contributions:** Conceptualization, R.C., J.K. and M.S.; methodology, R.C.; data curation and data analysis, R.C.; supervision, J.K.; resources, M.S.; writing—original draft preparation, R.C., J.K. and M.S.; writing—review and editing, J.K. and M.S.; project administration, J.K.; funding acquisition, J.K. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data are available within this document.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **Abbreviations**

The following abbreviations are used in this manuscript:

#### **References**

### *Article* **Offshore Wind Power Forecasting—A New Hyperparameter Optimisation Algorithm for Deep Learning Models**

**Shahram Hanifi 1,\*, Saeid Lotfian 2,\*, Hossein Zare-Behtash 1 and Andrea Cammarano 1**

**Abstract:** The main obstacle against the penetration of wind power into the power grid is its high variability in terms of wind speed fluctuations. Accurate power forecasting, while making maintenance more efficient, leads to the profit maximisation of power traders, whether for a wind turbine or a wind farm. Machine learning (ML) models are recognised as an accurate and fast method of wind power prediction, but their accuracy depends on the selection of the correct hyperparameters. The incorrect choice of hyperparameters will make it impossible to extract the maximum performance of the ML models, which is attributed to the weakness of the forecasting models. This paper uses a novel optimisation algorithm to tune the long short-term memory (LSTM) model for short-term wind power forecasting. The proposed method improves the power prediction accuracy and accelerates the optimisation process. Historical power data of an offshore wind turbine in Scotland is utilised to validate the proposed method and compare its outcome with regular ML models tuned by grid search. The results revealed the significant effect of the optimisation algorithm on the forecasting models' performance, with improvements of the RMSE of 7.89, 5.9, and 2.65 percent, compared to the persistence and conventional grid search-tuned Auto-Regressive Integrated Moving Average (ARIMA) and LSTM models.

**Keywords:** auto-regressive integrated moving average (ARIMA); long short-term memory (LSTM); Optuna; isolation forest (IF); elliptic envelope (EE); one-class support vector machine (OCSVM)

#### **1. Introduction**

Undoubtedly, to accelerate economic growth, power production through renewable energy sources needs to increase because conventional methods such as using fossil fuels have irreparable consequences, including pollution, climate change, and the depletion of the ozone layer [1].

In recent decades, various renewable energies, such as wind, solar, waves, etc., have received increasing attention. Among all these energies, wind power has played the most important role in replacing fossil fuels [2]. As reported by the World Wind Energy Council, the installed global capacity of wind energy in the world in 2021 has reached 837 GW, with an increase of 92 GW compared to 2020 [3]. Figure 1 shows the global wind power installed capacity increment over the past 21 years [3]. In this figure, the blue columns represent the capacity of installed wind power on land, while the red columns represent the offshore installed wind energy.

One main obstacle hindering the increase of wind power penetration into the power grid is the production uncertainty due to fluctuations in wind speed [1]. Therefore, adequate planning in electricity distribution to meet consumers' demand, determining the best time for operation and maintenance, and the fairest pricing on the market requires accurate wind power forecasting in the upcoming time steps.

**Citation:** Hanifi, S.; Lotfian, S.; Zare-Behtash, H.; Cammarano, A. Offshore Wind Power Forecasting—A New Hyperparameter Optimisation Algorithm for Deep Learning Models. *Energies* **2022**, *15*, 6919. https:// doi.org/10.3390/en15196919

Academic Editors: Paweł Piotrowski, Grzegorz Dudek and Dariusz Baczy ´nski

Received: 2 September 2022 Accepted: 16 September 2022 Published: 21 September 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

**Figure 1.** Global wind power installed capacity increment during the last 21 years [3].

Hanifi et al. [1] categorised wind power forecasting into three main methods, including physical, statistical, and hybrid approaches. Physical methods utilise numerical weather prediction (NWP) data, wind turbine geographic descriptions, and weather information to predict wind power [1]. These methods are computationally complex and very sensitive to initial information [2]. On the other hand, statistical methods work based on building an accurate mapping between input variables (such as NWP data, historical data, etc.) and target variables (wind speed or wind power). These methods include two main approaches: time-series-based methods and machine learning (ML) approaches [1]. Timeseries-based methods can predict wind speed or wind power based on the history of the predicted variable itself. They can recognise the concealed random features of wind speed and are used for very short-term (minutes to a few hours) forecasting. The Auto-Regressive Integrated Moving Average (ARIMA) model proposed by Box–Jenkins [4] is one of the common statistical methods which is used in various research. For example, in Western Australia, Yatiyana et al. [5] applied the ARIMA model for wind speed and direction forecasting. They proved that their proposed model could predict wind speed and direction with a maximum of 5% and 16% error, respectively. Firat et al. [6] proposed an autoregressive (AR) wind speed prediction model for a wind farm in the Netherlands. They used six years of hourly wind speed and achieved a high accuracy for 2–14 h ahead. In another study, De Felice et al. [7] applied 14 months of temperature readings in Italy to train an ARIMA model for electricity demand prediction. Their proposed method demonstrated higher accuracy, particularly in hot locations, compared with persistence methods. Duran et al. [8] proposed a method to combine AR and exogenous variable (ARX) models to predict the wind power generation in a wind farm located in Spain up to one day in advance. They used different model orders and training periods to prove that the application of the AR models presents lower errors than a persistent model. Kavasseri et al. [9] examined the application of fractional ARIMA models to predict wind farm hourly average wind speed for one- and two-day-ahead time horizons. The results of the predictions showed a 42% improvement compared to persistent methods. Later, the predicted wind speeds were applied to the power curve of an operating wind turbine to predict the relevant wind powers. In another study, Torres et al. [10] used the ARMA and the persistence model for hourly average wind speed forecasting up to 10 h ahead. The ARMA model demonstrated a better performance compared to the persistence method, with a 12% to 20% lower root mean square error (RMSE) when forecasting 10 h in advance.

ML methods such as neural networks (NNs) can establish deductive models by learning dependencies between input and output variables. These methods are easy to create, do not require further geographic information, and can predict over longer timeframes. One of the common ML methods is the LSTM model, which can address the long-term dependency issues [11], which is important in forecasting time-series with long input sequences [12].

LSTM is variously used in research for wind power prediction. For instance, Zhang et al. [13] proposed an LSTM wind power forecasting model for three wind turbines of a wind farm in China. They utilised three months of wind speed and historical power data and achieved the highest forecasting accuracy in a one-to-five time-steps ahead compared to the radial basis function (RBF) and deep belief network (DBN). Fu et al. [14] demonstrated LSTM and gated recurrent unit (GRU) for a one-to-four step-ahead forecasting of a 3 MW wind turbine in China, based on the first three-month dataset of 2014, with a resolution of 15 min. The comparison with ARIMA and support vector machine (SVM) methods showed the superiority of their proposed methods. Cali and Sharma [15] proposed an LSTM-based model with one hidden layer for 1 to 24 h ahead of wind power forecasting. The model was trained with 9-month data and evaluated in the last three months of 2016. They used nine combinations of input data, including wind speed at various levels, wind direction, temperature, and surface pressure. They demonstrated that temperature, wind speed, and direction positively impacted model performance; however, adding surface pressure to the input features led to worse performance.

As well as the training data, ML models' accuracy strongly depends on the adequate selection of their parameters and hyperparameters. The parameters of ML models (e.g., the weights of each neuron) are determined during the training process of the algorithm. In contrast, hyperparameters are not directly learnt by the learning algorithm and need to be specified outside the training process. The main role of the hyperparameters is to control the capacity of the models in learning dependencies. They also prevent overfitting and improve the generalisation of the algorithm. Hyperparameter optimisation or tuning improves forecasting accuracy and reduces models' complexity [16].

The literature's most common hyperparameter tuning methods are the grid search and random search. Grid search can be used for simple models with a few parameters. The calculation will be extremely time-consuming by increasing the number of parameters and expanding the space of the possible configurations [17]. Therefore, researchers usually consider a narrow range of hyperparameters during the grid search [16]. On the other hand, a random search algorithm looks randomly for a set of combinations rather than searching for better results.

Both these search methods generate all candidate combinations of hyperparameters upfront and then evaluate them in parallel. Based on the evaluation of all combinations, the best hyperparameters can be selected. Trying all possible combinations is very costly; as a result, it is vital to develop advanced techniques to intelligently select which hyperparameters to assess and then decide where to sample next after evaluating their quality.

The advanced optimisation of the ML-based time-series forecasting models for wind turbine-related predictions remained untouched. However, a few studies have proposed methods for optimising ARIMA and LSTM models within other applications than wind power forecasting. For example, Al-Douri et al. [18] designed a genetic algorithm (GA) to find the best parameters of an ARIMA model for the better cost prediction of used fans in Swedish road tunnels, and provided results which proved a significant improvement in data forecasting. In another study, F. Shahid et al. [19] employed GA to optimise the window size and neuron numbers of LSTM layers. This approach improved the power prediction accuracy of wind farms in Europe by up to about 30% compared to existing methods such as support vector regressors.

As the review of the literature indicates, several examples use linear and nonlinear regression models for challenges related to predicting wind power. Each study provides the use of one model type or a comparison of various model types. Nevertheless, without the tuning and selection of the hyperparameters, it is not possible to obtain their maximum benefit [16]. This advanced tuning method plays an important role when the hyperparameter search space grows exponentially, and the use of exhaustive grid search becomes extremely time-consuming.

This paper proposes a framework for developing accurate and robust ML models for wind power forecasting. The framework outlines the model development procedure from data engineering to precision evaluation and fine-tuning. Furthermore, an advanced algorithm is utilised to optimise wind power forecasting models to reduce time calculation costs, as well as to improve accuracy. For the case study, two ML models were selected: the LSTM model, which is proven to have remarkable prediction performance on time-seriesbased models, and ARIMA, a traditional model, for the purpose of benchmarking.

The novelty of this work lies in developing a short-term wind power forecasting model through an intelligent application of the long short-term memory (LSTM) model, while a new optimisation algorithm tunes its main hyperparameters. In addition, the distinguished aspects of the methodology are summarised, based on importance, as follows:

The rest of this paper is organised as follows: Section 2 discusses the optimisation process, the forecasting models, and the studied supervisory control and data acquisition (SCADA) data. This section includes the steps taken for preprocessing, resampling, and outlier treatment. Section 3 presents the results of the trained, optimised LSTM model in terms of model accuracy and the time cost compared to other prediction methods. Finally, Section 4 summarises the paper's contributions.

#### **2. Methodology**

The proposed procedure of this study is illustrated in Figure 2. At the beginning of this study, three required features, including the time stamps, wind speeds, and active wind powers, are selected to improve the computational time. At the next step, negative power values are removed or replaced. This data preprocessing is followed by resampling the dataset and removing outliers in three different ways. After finishing the data preprocessing and providing proper data for forecasting, data predictability and stationarity are assessed as two important specifications for accurate power forecasting. Afterwards, three different approaches are employed for forecasting, and their best performance is gained by the selection of their most appropriate hyperparameters.

**Figure 2.** Diagram of applied methodology.

#### *2.1. ARIMA Model*

In this study, the standard approach of the Box–Jenkins method [20] was traced for the ARIMA model development. The ARIMA model is a widely used set of statistical models for analysing and predicting time-series data [21]. This model can be expressed as [22]:

$$X\_t = \phi\_1 X\_{t-1} + \phi\_2 X\_{t-2} + \dots + \phi\_p X\_{t-p} + \varepsilon\_t - \theta\_1 \varepsilon\_{t-1} - \theta\_2 \varepsilon\_{t-2} - \dots - \theta\_q \varepsilon\_{t-q} \tag{1}$$

While *φt* and *θt* are coefficients, *p*, *q*, and *d* are the lag number of observations in the model, the order of moving average, and the degree of difference, respectively. Degree of difference (*d*) values greater than 0 imply that the data has been nonstationary but has become stationary after some degree of difference.

The ARIMA model combines the AR, moving average (MA), and the Integrated (I) components, which denotes the data substitution with the value of the difference between its values and the preceding values [23]. The forecasting accuracy of the ARIMA model depends on selecting the most appropriate combination of *p*, *d*, and *q*. Normally, for small data sets, the autocorrelation function (ACF) and partial autocorrelation function (PACF) can be used to determine which AR or MA component should be selected in the ARIMA model [24].

These two factors, which can be graphically plotted, are widely used elements in analysing and predicting time-series. They highlight the relationship between an observation and the observations' value at prior time steps. The difference between ACF and PACF is that, in PACF, while assessing the relationship between observation of two time steps, the relationships of the intervening observations are removed. Figure 3a,b show the observations' ACF and PACF plots. An appropriate ARIMA model can be selected based on the simple explanations in Table 1 [9], and the value of *d* (degree of difference) depends on the number of differencing until the data is stationary.

**Figure 3.** ACF (**a**) and PACF (**b**) plots for generated power of LDT. The blue points represent the value of autocorrelation and partial autocorrelation of different time lags.

The ARIMA model forecasting steps after resampling and outlier treatment can be seen in Figure 4. The first step is assessing the stationarity of the time-series. Stationary is one of the assumptions during time-series modelling, which shows the consistency of the summary statistics of the observations.

**Figure 4.** Flowchart of ARIMA and LSTM wind power forecasting models.

When a time-series is stationary, it means that the statistical properties of the timeseries (such as mean, variance, and autocorrelation) do not change over time. This property can be violated by having any trend, seasonality, and other time-dependent structures. There are two main methods for the stationarity assessment of time-series, the visualisation approach and the augmented Dickey–Fuller (ADF) test. The visualisation method uses graphs to show whether the standard deviation changes over time. On the other hand, the ADF method is a statistical significance test that compares the *p*-value with the critical values and does hypothesis testing. This test makes the stationarity of data clear at different levels of confidence.

Regarding the data used in this study, due to the high number of observations and wide dispersion, it is not possible to check stationarity through the visualisation method. Therefore, in this study, the ADF method was used.

The ADF test's execution provides a *p*-value which, by comparing it with a threshold (such as 5% or 1%), can identify the stationarity of the data. Nonstationary data in this step need to be changed to stationary by methods such as differencing. After ensuring the timeseries is stationary, a persistence method as a baseline is created. Then, through a detailed grid search, the best hyperparameters for the ARIMA forecasting for each preprocessed data were found. The last step is ARIMA forecasting and comparing its error with the error of the persistence method.

#### *2.2. LSTM Model*

The recurrent neural network (RNN) is a model in which the connection of its units creates cycles. RNN has a high ability to represent all dynamics. However, its effectiveness is affected by the limitations of the learning process. The main limitation of gradient-based methods that use back propagation is their path integral time-dependence on assigned weight [13]. When the time lag between the input signal and the target signal increases to more than 5–10 time-steps, the normal RNN loses the learning ability, and the backpropagation error either vanishes or explodes. This error elimination raises the question of whether normal RNNs can show practical benefits for feed-forward networks. To address this problem, the LSTM has been developed based on memory cells. The LSTM consists of a recurrently attached linear unit known as the constant error carousel (CEC). CECs, by keeping the local error backflow constant, mitigate the gradient's vanishing problem [25]. They can be trained by adjusting both the back propagation over time and the real-time recurrent learning algorithm [26]. Figure 5 shows the typical structure of the LSTM.

As can be seen, there are three gate units in a basic LSTM cell, including the input, output, and forget gates. The gate activation vectors of *it*, *ot* and *ft* for input, output, and forget gates, respectively, are calculated in Equations (2)–(4).

$$\dot{a}\_{l} = \sigma\_{l} (\mathsf{W}\_{l}\mathsf{x}\_{l} + L\_{l}\mathsf{h}\_{t-1} + b\_{i}) \tag{2}$$

$$\boldsymbol{\sigma}\_{t} = \boldsymbol{\sigma}\_{l} (\mathcal{W}\_{o}\mathbf{x}\_{t} + \mathcal{U}\_{o}\mathbf{h}\_{t-1} + \mathbf{b}\_{o}) \tag{3}$$

$$f\_t = \sigma\_l \left(\mathcal{W}\_f \mathbf{x}\_t + \mathcal{U}\_f h\_{t-1} + b\_f\right) \tag{4}$$

In these equations, *Wi*, *Wo*, *Wf Ui*, *Uo*, and *Uf* represent the assigned weights, and *bi*, *bo*, and *bf* represent the biases in conjunction with relevant activation functions *σl*. In addition, *xt* is the neuron input at time step t, and the cell state vector at time step t − 1 is *ht*−1. As shown in Equation (5), the next evaluated value of the state *S*%*t* can be calculated based on the relevant activation function *σs*.

$$\mathcal{S}\_t = \sigma\_s(\mathcal{W}\_s \mathbf{x}\_t + \mathcal{U}\_s \mathbf{h}\_{t-1} + \mathbf{b}\_s) \tag{5}$$

In Equation (7), the newly assessed value of *S*%*t* and the prior cell state *St*−1 are used to calculate cell state *St*, which by itself will be used with the output gate control signal *ot* and the activation function *σlh* to obtain the overall output *ht* according to Equation (8).

$$S\_t = f\_t S\_{t-1} + i\_t \dot{S}\_t \tag{6}$$

$$h\_t = o\_t \sigma\_{lh} \left( \mathcal{S}\_t \right) \tag{7}$$

**Figure 5.** Typical structure of the LSTM.

As can be seen in Equations (6) and (7), the output *ht* is dependent on the state *St* of the LSTM cell and the activation function *σlh* that is usually tanh (x). The state *St* depends on the state of the prior step *St*−1 as well as the new value of the state *S*%*t*.

In accordance with all the relations mentioned above, the function of the LSTM model can be concluded as:

Specifying the best LSTM model for wind power forecasting requires the determination of the neural network's best combination of hyperparameters. LSTMs have five main hyperparameters, including the number of lag observations as inputs of the model, the quantity of LSTM units for the hidden layer, the model exposure frequency to the whole training dataset, the number of samples inside an epoch in each weight updating, and finally, the used difference order for making nonstationary data stationary.

#### *2.3. Grid Search for ARIMA and LSTM Models*

ARIMA model factors (i.e., *p*, *d*, and *q*) can be estimated through iterative trial and error by revising the ACF and PACF plot. This part of defining the ARIMA forecasting model can be very challenging and time-consuming, leading to prediction errors. As a result, researchers attempt to find these hyperparameters using an automatic grid search approach. Similar to the ARIMA model, specifying the best LSTM model for wind power forecasting requires the determination of the best combination of hyperparameters in this neural network. This study also specified a grid of the LSTM parameters to iterate. An LSTM model is created based on each combination, and its forecasting accuracy is assessed by calculating its RMSE.

#### *2.4. Persistence Method*

It is vital to create a baseline for any time-series prediction approach. As a reference, for comparing all modelling approaches, this baseline can show how well a model makes predictions. Models which perform worse than the performance level of the baseline can be ignored.

Benchmarks for forecasting problems need to be very simple to train, fast to implement, and repeatable. The persistence model is one of the most commonly used references for wind speed and power prediction (short-term forecasting methods in particular). Based on the definition of this method, wind power in the future will be equivalent to the generated power in the present [27], as given by Equation (8):

$$
\hat{P}\_{t+k/t} = P\_t \tag{8}
$$

where *Pt* is the measured wind power at time *t* and *P*ˆ *t*+*k*/*t* is the predicted wind power for the future time *k*. This model performs better than most short-term physical and statistical forecasting methods. Therefore, it is still widely used in very short-term prediction [28]. This research uses the persistence model to compare the performance of the ARIMA and LSTM models for different datasets.

#### *2.5. Hyperparameter Optimisation with Optuna*

This study uses the Optuna optimisation method to optimise the forecasting models. Optuna is an open-source optimisation software with several advantages over the other optimisation frameworks [29]. Other optimisation tools usually differ depending on the algorithm used to select the parameters. For example, GPyOpt and Spearmint [30] apply Gaussian processes, SMAC [31] employs random forests, and Hyperopt [32] uses a treestructured Parzen estimator (TPE). These methods have three main drawbacks. Firstly, they need the parameter search space to be statically defined by the user, a process that is extremely hard for large-scale experiments with many possible parameters. Furthermore, they do not have an efficient pruning strategy for high-performance optimisation when accessing limited resources. In addition, they cannot handle large-scale experiments with minimal setup requirements. On the other hand, Optuna, with a define-by-run design, enables the user to create the search space dynamically. This optimisation framework is an open-source, easy-to-set-up package that benefits effective sampling and pruning algorithms [29]. Optuna optimises the model through minimising/maximising an objective function (here, the RMSE of the forecasted wind power rather than the real generated values) that assumes a group of hyperparameters as input and returns its validation core. The optimisation process is called a study, and each objective function's evaluation is called a trial [29].

At the beginning of the optimisation, the user is asked to provide the search space for the dynamic generation of the hyperparameters for each trial. Then, the model builds the objective function by interacting with the trial object. After this step, the next hyperparameter selection is based on the history of previously evaluated trials. This algorithm optimises ML models in two steps. First, a search strategy determines a set of parameters to be examined, and second, a performance assessment strategy known as a pruning algorithm excludes the improper parameters based on the estimation of the value of the currently investigated parameters [29].

Since the initial prediction accuracy assessment of the ARIMA and LSTM models (both tuned by grid search) highlighted the better performance of the LSTM model compared to ARIMA, it was decided to apply the optimisation framework only to the LSTM model.

In this way, the hyperparameter ranges of the LSTM model increased from what was examined in its grid search to wider ranges, as shown in Table 2. In other words, the hyperparameter combinations increased from 48 combinations to more than a million combinations.

**Table 2.** Hyperparameter ranges and their total combinations in LSTM–grid search and LSTM– Optuna methods.

#### *2.6. Wind Power Dataset*

The source SCADA data are measured at a 1 Hz frequency from the Levenmouth Demonstration Turbine (LDT), an offshore wind turbine which is located just 50 m from the coast at Leven, a seaside town in Fife, Scotland [33]. This wind turbine was acquired by the Offshore Renewable Energy (ORE) Catapult in 2015, while its construction was completed by Samsung in October 2013 [34].

ORE Catapult's wind turbine is a three-bladed upwind turbine installed on a jacket structure [25]. The turbine is ranked to work at 7 MW, but to decrease the noise, it is limited to operating at the highest power of 6.5 MW [33]. This turbine's rotor diameter is 171.2 m, and its hub height is 110.6 m. Each blade of this turbine measures 83.5 m and weighs 30 tons. The defined cut-in speed for this turbine is 3.5 m/s, which means its electricity generation will start when wind speeds reach this speed. It will shut down if the wind is blowing too hard (roughly 25 m/s) so to prevent equipment damage. Its operating temperature is between −10 ◦C to +25 ◦C, and it has been designed to work for 25 years [35]. Figure 6 shows the configuration and main parameters of the LDT.

**Figure 6.** Main parameters and schematic of Levenmouth wind turbine [35].

#### *2.7. Feature Selection*

This study recorded the SCADA datasets for five months, from 1 January 2019 to 31 May 2019, at a 1 Hz frequency (with one-second intervals). Each timestamp in this time-series data includes 574 different observations, including the generated power, wind speed at different levels, blade pitch angle, nacelle orientation, etc. At the beginning of the data processing, a feature selection was carried out to decrease the size of the dataset to reduce the computation time by excluding unnecessary variables. This process was vital to making this study possible. All variables except the time stamp, wind speed, and active power were removed at this stage, which was useless in the ARIMA and univariate LSTM forecasting methods. Keeping the wind speed variable was vital in this project, as it verified the accuracy of generated power. For example, failure to generate power when high wind speeds were recorded was recognised as a stop in power generation due to reasons such as maintenance. After removing the redundant information, observations of wind speed and active power were plotted as shown in Figure 7a,b.

**Figure 7.** Wind speed observations (**a**), wind active power observations (**b**).

The histograms of this dataset for wind speed and active power are presented in Figure 8a,b, and Table 3 shows their statistical descriptions.

**Figure 8.** Histogram of active power (**a**) and wind speed (**b**).

**Table 3.** Statistical descriptions of the SCADA datasets.

#### *2.8. Obvious Outlier Removal*

An initial assessment of Figure 7b specified that a large part of the recorded generated power at the end of this time-series (May 2019) equals zero. Usually, the generated power of a turbine can be zero when no wind is blown. However, the evaluation of Figure 7a shows a continuous wind blowing with fluctuations similar to previous months. Therefore, it is speculated that the turbine was out of production during this period. Based on this assumption, it was decided that this month (May 2019) should be removed entirely from the dataset. The time-series after this omission was reduced to four months, from 1 January 2019 to 30 April 2019. A closer look at the active power, as shown in Figure 9, revealed another obvious error in the SCADA data, the existence of negative values. Negative values are values of which there is no practical meaning in wind power generation. Shen et al. [36] believe that these values represent time stamps when turbine blades do not rotate, but the turbine's control system needs electricity [36]. These values need to be eliminated along with the corresponding parameters of the same timestamp for better forecasting results [25]. Since the elimination of these negative values disrupts the time continuity of the time-series, and can possibly lead to errors in wind power prediction, at this stage it was decided to create and assess three types of datasets based on different actions against negative values. Assessment of the impact of these actions on forecasting accuracy became another goal of this study.

**Figure 9.** Wind power observations (only power values under 1000 kW are shown). The dotted red line indicates the power value of zero (the boundary of negative/positive values).

These three preprocessing methods against the negative values are:

#### *2.9. Resampling*

The effect of wind turbulence as one of the obstacles to increasing the wind energy penetration in energy markets is more significant in horizontal axis wind turbines. This is because the wind speed and direction change rapidly after hitting swept blade rotors. Therefore, the amount of wind speed measurements by installed anemometers are not equal to the speed of the wind flow hitting turbine blades [25]. These differences, which lead to a decrease in the correlation between the measured wind speed and the output power, and then scattering of the power curve, can be resolved by averaging the samples in a reasonable average period [25]. The SCADA data for this study was recorded with a 1 Hz frequency; as a result, it was possible to create multiple averaged sets for removing the

mentioned obstacle. According to a review conducted by Hanifi et al. [1], the maximum sampling rate used for wind speed and power forecasting in the previous research is 10 min. This is equivalent to an average time that the international standard for power performance measurements of electricity-producing wind turbines (IEC 61400-12-1) establishes for large wind turbines [37]. Based on the IEC 61400-12-1 and reviewed literature, the data presented here was averaged for each 10 min of data collection. Figure 10a,b show the wind power curves for the original and 10 min resampled data.

**Figure 10.** Wind power curves. (**a**) Original 1s data (**b**) and 10 min resampled data.

#### *2.10. Anomalies Detection and Treatment*

Outliers in a dataset are specific data points that are different or far from most other regular data points [38]. Undetected or improperly treated anomalies can adversely affect wind power forecasting applications. They may be biased with high prediction errors [38].

There are various reasons for having outliers among wind turbine and wind farm measurements, including wind turbine downtime [36], data transmission, processing or management failure [39], data acquisition failure [40], electromagnetic disturbance [36], wind turbine control system fault (such as the pitch control system fault) [41], damage of the blades or the existence of ice or dust [42], shading effect of neighbouring turbines, fluctuation of air density [43], etc.

Figure 11 shows four different types of anomalies in the current SCADA data. Category A points have negative, zero, or low values of generated power during speeds larger than the cut-in speed [25]. The leading causes of these outliers are wrong wind power measurements, wind turbine failure, and unexpected maintenance. Wind speed sensors and communication errors cause category B outliers. The mid-curve outliers (category C) represent power values lower than ideal—this is caused by the down-rating of the wind turbines and data acquisition. Outliers in category D are scattered irregular points due to faulty sensors exacerbated during harsh weather circumstances [36].

There are different methods for anomaly detection in machine learning, such as Density-Based Spatial Clustering of Applications with Noise (DBSCAN), IF, local outlier factor, and EE. In this study, three common methods for wind power forecasting are investigated. EE is used based on the assumptions described in [44]. IF, which is an unsupervised learning algorithm, recognises anomalies by isolating them in the data. This algorithm works based on two main features of anomalies, that they are few and different. The one-class support vector machine (OCSVM) is a common unsupervised learning algorithm for outlier detection, assuming rare anomalies create a boundary for most data, and considering data points out of the boundary as outliers [45]. This method of outlier detection and treatment chose the third method.

**Figure 11.** Observed anomalies coupled with the power curve of the 1 Hz original data. (A) Low power output in high wind speeds in turbine failure cases; (B) Outliers due to the wind speed sensor and communication errors; (C) Power outputs less than the rated power as a result of the turbine's down-rating; (D) Scattered outliers caused by sensor malfunctions or noise in signal processing.

#### **3. Experimental Results and Discussion**

This research employs packages and subroutines written in Python to implement the proposed algorithms. A PC with an Intel Core i5–7300 32.6 GHz CPU and 8 GB RAM (without any GPU processing) was used to run the experiments. Three outlier detection methods, which were described in Section 2.8, were used to detect and remove the outliers of the resampled dataset. The results of these treatments can be seen in Figures 12–14:

**Figure 12.** Elliptic envelope application for outlier detection and treatment. The blue points represent the normal data, and the red represents the detected anomalies.

This study considers six different preprocessing methods based on applying three different outlier detection methods and three approaches against the negative power values (Table 4). Different cases of preprocessed data are fed to the ARIMA and LSTM forecasting models. The grid search method is applied for the initial hyperparameter tuning; Table 4 shows the selected hyperparameters for the ARIMA and LSTM models. As expected, the values of the hyperparameters vary depending on the different employed preprocessing methods (Table 4).

**Figure 13.** Isolation forest application for outlier detection and treatment. The blue points represent the normal data, and the red represents the detected anomalies.

**Figure 14.** OCSVM application for outlier detection and treatment. The blue points represent the normal data, and the red represents the detected anomalies.

**Table 4.** Best ARIMA and LSTM hyperparameters resulting from the grid search.

\*: Mean value has been calculated after removing negative values.

After selecting the best ARIMA and LSTM prediction methods, both models were trained by the first 95% part of the dataset (as training data) to make predictions for the last 5% of the dataset. The predicted values were compared with the measured values to determine the RMSE of each forecasting process. Table 5 provides the RMSE values of the ARIMA, LSTM, and persistence methods.

**Table 5.** RMSE values of persistence, ARIMA, and LSTM models for six different treated case data.

1 Replaced by mean value calculated after removing negative values. 2 Replaced by nearest positive value.

Comparing the RMSE values of all three models (Table 5) for case data 1, 2, and 3 clarifies that the complete elimination of the negative values (without any replacement) will lead to worse forecasting. The highest RMSE value of case 3 means that removing the negative values will decrease the forecasting accuracy. One of the reasons for this performance drop can be the creation of discontinuity in the dataset.

Regarding the best specific value to be considered instead of negatives, a comparison of case data 1 and 2 proves that replacing the negative values with the average wind power values has a better impact than replacing them with the nearest (neighbour) positive value. Replacing the negative values with the average values can lead to about a 15% forecasting improvement for ARIMA and 11% for the LSTM models.

The results also highlight the importance of dealing with outliers in wind power forecasting. Cases 4, 5, and 6, representing the outlier removed data, show a significant enhancement of the accuracy rather than the other cases, without any action against the anomalies. Comparing the error levels of case data 3 with cases 4, 5, and 6 (for both ARIMA and LSTM models) shows a 30% to 38% forecasting improvement by the elimination of the outliers, either by isolation forest, elliptic envelope, or the one-class SVM outlier detection methods.

The assessment of the RMSE values of cases 4, 5, and 6 show that the IF and EE outlier detection methods overcome the OCSVM method. An elliptic envelope can improve forecasting performance up to 9.61% and 8.92% rather than OCSVM for ARIMA and LSTM methods. This performance enhancement can reach 9.96% and 9.64% for ARIMA and LSTM, respectively, by applying the isolation forest.

As shown in Table 5, the ARIMA and LSTM methods for all the treated case data have better performances than the persistence methods. This is understandable if one remembers that, in the persistence method, only one preceding step data is used for forecasting, whilst the ARIMA and LSTM models consider a more extensive range of prior data.

It is also clear that the LSTM performs better than the ARIMA almost for all approaches against the negative values and outliers. This is probably due to the fact that LSTMs are better equipped to learn long-term correlation. In addition, the LSTM can better capture the nonlinear dependencies between the features.

In this study, because of the better prediction performance of the LSTM model compared to the ARIMA model, the proposed optimisation algorithm is applied to the LSTM model to tune its hyperparameters even more. As discussed in Section 2.5, the hyperparameter ranges of the LSTM model are increased from what was examined in its grid search to the wider ranges shown in Table 2.

The six preprocessed case data are again divided into the first 95% as the training dataset and the rest 5% as the test data. These divisions were developed to establish the same conditions and logically compare the new and previous methods. The developed optimisation algorithm, with the two described strategies, including search and pruning, started the selection of different combinations to minimise the RMSE value. Table 6 shows the new hyperparameters found by the Optuna optimisation algorithm, and Figure 15 shows the measured power values of the turbine and prediction results of all the forecasting

methods, including ARIMA, LSTM–grid, and LSTM–Optuna, for one of the datasets (data 4—removed negative values and removed outliers with the EE method).

**Table 6.** Best LSTM hyperparameters resulted from Optuna optimisation.

\*: Mean value has been calculated after removing negative values.

**Figure 15.** Comparison of measured wind power and forecasted values by ARIMA, LSTM–grid, and LSTME–Optuna models for data 4 (removed negative values and removed outliers with EE method).

As can be seen in Figure 15, the LSTM model optimised by Optuna can predict more accurately by better learning the wind power's short-term and long-term dependencies. The diagram illustrated in Figure 16 is plotted to better compare the error levels of the different wind power forecasting methods. It can be recognised from this diagram that the LSTM–Optuna approach follows rules similar to the ARIMA and LSTM–grid models. To achieve a higher prediction accuracy, it is essential to eliminate the outliers and replace the negative power values with the average wind power value.

Building the LSTM models based on the new values of the hyperparameters, as shown in Table 6, improves the prediction accuracy of the LSTM model in a range from 1.22% to 2.65% for different cases of preprocessed data. These accuracy improvements can be seen in Table 7.

The results show that the highest accuracy improvement is related to case 5, a case in which negative values were replaced with the mean power value and the outliers were removed through the IF method. A comparison of the required search times to find the best combination of the hyperparameters in LSTM–grid and LSTM–Optuna proves the faster performance of the proposed method, as it spends from 13.79% to 20.59% less time adjusting the model for the most accurate prediction (Table 8).

**Figure 16.** Error comparison of persistence, ARIMA, LSTM, and LSTM optimised by Optuna forecasting methods.

**Table 7.** A comparison of RMSE, the LSTM–grid search, and LSTM–Optuna methods.

**Table 8.** A comparison of the required tuning time of the LSTM–grid search and LSTM–Optuna methods.

#### **4. Conclusions**

This study addresses issues regarding inaccurate wind power prediction using ML approaches. As discussed in the reviewed literature, most previous research applied ML without advanced model optimisation. At the same time, in this paper, a novel concept of Optuna–LSTM is reported to expedite the process of selecting the hyperparameters and tuning the wind power forecasting models. This model not only reduces the time complexity of creating reliable models, but also improves the accuracy of the predictions.

To accurately evaluate the proposed model, SCADA data of an offshore wind turbine was preprocessed by eliminating its negative values and outliers to help find the best preprocessing method. The performance of the proposed forecasting was demonstrated through comparisons with the persistence, ARIMA, and LSTM models, which were already tuned by grid search. This comparison proved the better performance of the proposed model, with a range up to 7.89, 5.9, and 2.65 percent compared to the persistence and conventional grid-search-tuned ARIMA and LSTM models.

This study also highlights the importance of eliminating negative values in the power recordings. The results of this study confirmed that replacing the negative values with the average power value has the most positive effect on the forecasting accuracy. In addition, comparisons between several data cases showed the significant impact of the outlier treatment methods on the forecasting performance. The results proved that removing the

outliers by the isolation forest method improves the forecast accuracy compared to the elliptic envelope and OCSVM methods. This novel forecasting method combining the capacity of the LSTM model in the prediction of nonlinearities and the optimisation tool for better tuning the hyperparameters can be used for different time-series-based predictions.

**Author Contributions:** Conceptualisation, S.H. and S.L.; methodology, S.H. and S.L.; investigation, S.H. and S.L.; writing—original draft preparation, S.H.; writing—review and editing, S.L., H.Z.-B., A.C. and S.H.; supervision, S.L.; data curation, S.H. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the EPSRC Doctoral Training Partnership (EP/R513222/1).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author. The data are not publicly available because it also forms part of an ongoing study.

**Acknowledgments:** The authors are grateful to Xiaolei Liu for his assistance and contribution to this research.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Nomenclature**

#### **References**

### *Article* **A Spiking Neural Network Based Wind Power Forecasting Model for Neuromorphic Devices**

**Juan Manuel González Sopeña 1,\*, Vikram Pakrashi 2,3,4 and Bidisha Ghosh 1,5**

**Abstract:** Many authors have reported the use of deep learning techniques to model wind power forecasts. For shorter-term prediction horizons, the training and deployment of such models is hindered by their computational cost. Neuromorphic computing provides a new paradigm to overcome this barrier through the development of devices suited for applications where latency and low-energy consumption play a key role, as is the case in real-time short-term wind power forecasting. The use of biologically inspired algorithms adapted to the architecture of neuromorphic devices, such as spiking neural networks, is essential to maximize their potential. In this paper, we propose a short-term wind power forecasting model based on spiking neural networks adapted to the computational abilities of Loihi, a neuromorphic device developed by Intel. A case study is presented with real wind power generation data from Ireland to evaluate the ability of the proposed approach, reaching a normalised mean absolute error of 2.84 percent for one-step-ahead wind power forecasts. The study illustrates the plausibility of the development of neuromorphic devices aligned with the specific demands of the wind energy sector.

**Keywords:** neuromorphic computing; spiking neural network; short-term wind power forecasting

#### **1. Introduction**

A large number of machine learning (ML) and deep learning (DL) models have been developed and applied to time series data of a varied nature for tasks such as forecasting [1], classification [2], and clustering [3]. This trend has been also been observed in the field of wind power forecasting (WPF) [4], particularly the use of artificial neural networks (ANNs) [5], which are usually trained with the backpropagation algorithm [6]. Recurrent neural networks, such as the gated recurrent unit (GRUs) [7] and long short-term memory (LSTM) neurons [8], can learn temporal features on wind data, whereas convolutional neural networks (CNNs) capture spatial ones [9]. Other ML algorithms that have been applied in the literature are support-vector machines [10], random forests [11], gradient boosting machines [12], and neuro-fuzzy models [13,14]. DL methods such as deep neural networks are built by stacking multiple layers between the input and output layers to extract higher-level features from the data [15]. Deep neural architectures such as deep belief networks [16], deep convolutional networks [17] and N-BEATS [18] have been applied in the WPF literature. Furthermore, the abilities of ML/DL as a modelling tool have proven valuable for solar power forecasting [19] and renewable energy systems [20].

Accurate WPFs can be estimated using ML/DL architectures considering different types of data collected at a wind farm [21,22]. However, such models may be associated

**Citation:** González Sopeña, J.M.; Pakrashi, V.; Ghosh, B. A Spiking Neural Network Based Wind Power Forecasting Model for Neuromorphic Devices. *Energies* **2022**, *15*, 7256. https://doi.org/10.3390/en15197256

Academic Editors: Paweł Piotrowski, Grzegorz Dudek and Dariusz Baczy ´nski

Received: 5 September 2022 Accepted: 28 September 2022 Published: 2 October 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

with a high computational cost, a critical factor for edge computing [23], including those applications for renewable energy [24]. For instance, the low latency inherent in neuromorphic devices can be critical for transmission system operators to manage the grid in real time and the decision-making process of traders participating in electricity markets, specifically to correct their positions in intraday markets. Thus, neuromorphic computing provides an alternative to the computational complexity of ML/DL models [25] with the development of devices inspired by the energy-efficient nature of biological systems such as the Intel's Loihi chip [26]. The architecture of spiking neural networks (SNNs) [27] resemble more closely biological neurons, and are thus adequate for implementation in neuromorphic devices to unleash their potential in terms of low latency and lower energy consumption. However, training SNNs remains a challenge, as the well-known backpropagation algorithm cannot be applied, due to the non-differentiable nature of spikes. The current approaches to train spiking DL algorithms can be broadly divided into *online* and *offline* approaches [28]. Online approaches first implement an SNN in neuromorphic hardware, leveraging on-chip plasticity to train the spiking network and evolve its parameters with the arrival of new data [29]. This approach includes online approximations of the backpropagation algorithm [30,31] and evolving SNNs [32]. On the other hand, the SNN is trained before deploying the model for offline approaches. These can be further divided into two categories, considering how the training stage is performed. One possibility is to train a conventional ANN using the backpropagation algorithm, and later map the parameters into an equivalent SNN model [33]. This approach is known as *ANN-to-SNN conversion*. Alternatively, a direct training approach uses a variation of error backpropagation to optimize directly the parameters of an SNN [34].

In addition, the research community has been developing specific software platforms to implement applications based on SNNs. For instance, Nengo [35] is a software based on the principles of the Neural Engineering Framework (NEF), a theoretical framework to implement large-scale neural models with cognitive abilities [36]. This same software was later extended with the sister library NengoDL [37], aiming to combine the principles of neuromorphic modelling with the well-known deep learning framework TensorFlow [38] to build deep spiking neural models by ANN-to-SNN conversion. Alternatively, other frameworks can directly train SNNs, such as the Spike Layer Error Reassignment (SLAYER) algorithm proposed by Shrestha and Orchard [39]. Recently, in October 2021, Intel's Neuromorphic Computing Lab released the first version of Lava [40], an open-source software framework, to implement neuromorphic applications for the Intel Loihi architecture [41].

The SNN features are maximized within the framework provided by neuromorphic computing. However, no study has realistically attempted to model short-term WPFs using SNNs while considering the current computational abilities of neuromorphic devices to date. We have to remember that neuromorphic computing is still in its infancy, so the goal is to reach an acceptable level of performance to build up our knowledge regarding the implementation of spiking-based models in WPF, and not to outperform the current well-established neural network models [28]. Therefore, we propose a SNN model for shortterm WPF, adapted to the hardware capacity of the current state-of-the-art neuromorphic devices, particularly the neuromorphic chip Loihi developed by Intel. The aim of this study is not solely constrained to achieving highly accurate WPFs, but also to efficiently design WPF models that leverage the neuromorphic processors' power efficiency. The proposed forecasting approach was designed by applying the modelling framework provided by NengoDL to build spiking neuron models, and NengoLoihi, a complementary library, to implement such models on Loihi hardware.

The rest of this paper is structured as follows. Section 2 describes the ANN-to-SNN conversion method used to train spiking neural networks, as well as the spiking model architecture tailored to Loihi hardware. Section 3 presents a case study using this methodology for short-term WPF, using real data from an Irish wind farm. Section 4 contains the concluding remarks and the scope for future research work.

#### **2. Methodology**

NengoDL [37] is a modelling framework that includes tools to design biological neuronal models and the optimization methods used to train ML/DL models. Such optimization methods are usually incompatible with SNNs, as spikes are not differentiable. NengoDL links SNNs and these optimization methods by performing the necessary transformations to apply the ANN-to-SNN conversion method proposed by Hunsberger and Eliasmith [42], which allows for the use of a rate-based version of the spiking model in the training stage and the SNN for inference. The design of this rate-based approximation is key to successfully mapping its parameters into a spiking network, so the parameters of the network must be carefully tuned to ensure a minimal loss of performance during the conversion, and the architecture of the model must be tailored to subsequently build the network on Loihi hardware. Typically, six steps were followed to build and evaluate the performance of our proposed SNN model within the framework provided by NengoDL, as follows:

In the remainder of this section, we introduce how the ANN-to-SNN conversion is performed and the model architecture chosen to forecast wind power.

**Figure 1.** Spiking ReLU activation profile (based on DeWolf et al. [43]).

#### *2.1. ANN-to-SNN Conversion*

The non-differentiable nature of spikes impedes the use of the backpropagation algorithm to train spiking neurons [47]. ANN-to-SNN conversion sorts this out by mapping the parameters of a trained ANN to an equivalent SNN. Thus, the main challenge is how to train the non-spiking model so that there is only a small loss of performance in the conversion process. The first point is choosing an adequate spiking activation function. Cao et al. [48] established an equivalence between the ReLU activation function [49] and the spiking neuron's firing rate. The ANN-to-SNN conversion method implemented in NengoDL was proposed by Hunsberger and Eliasmith [42]. This method is valid for both linear (such as ReLU) and non-linear activation functions such as leaky integrate-and-fire (LIF) by smoothing the equivalent rate equation employed to train the ANN. To understand this, let us look at the equation governing the dynamics of an LIF neuron:

$$
\pi\_{\rm RC} \frac{dv(t)}{dt} = -v(t) + I(t) \tag{1}
$$

where *τRC* is the membrane time constant, *v*(*t*) is the membrane voltage, and *I*(*t*) is the input current. The neuron will fire a spike if it reaches a certain threshold V and after the potential is reset during a certain period of time (known as refractory period *τref*). The dynamics of the neuron are recovered after the refractory period *τref* is ended. If a constant input current is given to the neuron, the steady-state firing rate (i.e., the time that it takes to the neuron to reach the threshold to fire a spike) can be determined as:

$$r(j) = \left[\tau\_{ref} + \tau\_{RC} \log(1 + \frac{\mathcal{V}}{\rho(j - \mathcal{V})})\right]^{-1} \tag{2}$$

where *ρ*(*x*) = *max*(*x*, 0). However, this function is not completely differentiable, so the LIF rate equation is softened to address this problem and allow for use of the backpropagation algorithm [42]. The hard maximum *ρ* is replaced by a soft maximum *ρ* defined as:

$$\rho'(\mathbf{x}) = \gamma \log(1 + e^{\mathbf{x}/\gamma}) \tag{3}$$

After training the conventional ANN, the parameters of the SNN are identical to its non-spiking counterpart, with only the neurons themselves changing. The performance of the spiking network can be further enhanced by tuning additional parameters. For instance, if using a linear activation function for the spiking forecasting model, the spiking firing rate can easily be increased after training by applying a scale to the input weights of the neurons to make them spike at a faster rate. The output of the network is divided for the

same scale to not affect the behavior of the trained network. This way of proceeding is not optimal for non-linear activation functions. Instead, the firing rates can be optimized during training with regularization, so the firing rates are encouraged to spike at a certain firing rate [43]. Furthermore, a synaptic filter can be applied to reduce any possible noise found in the output of the spiking network.

#### *2.2. Spiking Model Architecture*

The model architecture (Figure 2) is slightly different to conventional ANNs, as it has to be adapted to the requirements of the Loihi hardware. The first distinctive feature of this network is the *off-chip* layer. This layer is a prerequisite to transmitting any information with the hardware, as it only communicates with spikes. Thus, this initial layer is run off-chip and converts the input into spikes [44]. The rest of the network is run on the hardware. A convolutional (*conv-layer*) and a regular fully connected layer (*dense-layer*) are used to process the data and generate the forecast. The convolutional layer is constituted by filters in the form of convolutions:

$$(f \ast g)(t) = \int\_{-\infty}^{\infty} f(\tau)g(t - \tau)d\tau \tag{4}$$

where (*f* ∗ *g*) indicates the convolution between the functions *f* and *g*, in which the function *f* can be considered as a *filter* or *kernel* and *g* as the input data. On the right hand side, *g*(*t* − *τ*) indicates that the input data *g* are reversed and shifted to a certain time *t*. It is important to notice that not all types of neural networks are currently available in this ANNto-SNN conversion framework (e.g., LSTM neurons are not supported). The activation function of all these three layers will be a spiking ReLU activation for inference. The equivalent ANN used during training follows the same architecture, including the *off-chip* layer, although it only behaves as a regular convolutional layer in this case. ReLU activation functions are used instead during the training stage.

Following up on our previous work [46,50], this model architecture is applied to 10-min resolution wind power data. The data used in this paper were collected for a wind turbine of a site located in Ireland (the exact location cannot be disclosed due to confidentiality reasons) for a two-year-and-a-half period (from January 2017 to June 2019). As input, the model uses previous wind power observations to provide the one-step-ahead forecasts as the output.

The wind power data were preprocessed using the variational mode decomposition (VMD) algorithm [51]. In particular, the data were decomposed into 8 subseries (known as modes) with different levels of complexity [52], giving us the opportunity to examine and adapt the SNN architecture under varied conditions. The forecasts of each mode were later aggregated to subsequently estimate the WPF [53].

**Figure 2.** SNN model architecture.

#### **3. Results**

First, examples using a synthetic sine wave signal and load data are given to clarify some details of the steps that need to be followed to successfully convert an ANN model into a spiking one. Later, a case study is presented using data from an Irish wind farm.

#### *3.1. Synthetic Signal Forecasting*

Before applying the methodology to wind power data, let us present an example with a more simple signal (a synthetic sine wave) to clarify and further explain the details of tuning the parameters to achieve a good performance with the spiking network. For simplicity, the example using this signal was run within the NengoDL framework, so any additional parameters used to implement the model on Loihi hardware can be dismissed (such as the off-chip layer); therefore, a basic feedforward neural network (FFNN) model was used instead of the previously described model architecture, which suffices to accurately predict such a basic signal.

During the initial evaluation of the spiking network model (Steps 3 and 4), considering the discretization of the activation function required for Loihi hardware is of importance to posteriorly transfer our model without a significant drop in performance. Therefore, we must be particularly careful when scaling the firing rate of the spikes, as very high rates will not work on Loihi hardware. Let us examine the implications of disregarding this point with the example shown in Figure 3: we build the FFNN model (Step 1) and train it with a rate-based (i.e., non-spiking) ReLU activation (Step 2). Then, we replace the activation for its spiking counterpart, scaling the firing rate with a high enough value (Step 3). The neural activities of three neurons when presenting an input are shown in

Figure 3a,b, having replaced the ReLU activation function with the theoretical spiking ReLU and the discretized version with Loihi, respectively. Two of these neurons (shown in green and yellow) fire very fast in the first case, but their behavior is diminished in the second one due to the activation profile, impacting the performance of the model when all the input vectors conforming with the testing set are presented to the network (Step 4), as displayed in Figure 3c). Thus, the firing rate of this network should be lowered to satisfy the hardware specifications required in the following steps to implement the model on neuromorphic devices.

**Figure 3.** (**a**) Neural activities using a spiking ReLU activation for inference (one input vector is shown to the network during 50 timesteps), (**b**) neural activities using the discretized version of the spiking ReLU activation, and (**c**) predictions over the testing set.

The tuning of the firing rate scale, as well as the amplitude of the spikes, are essential to achieve a good forecasting accuracy while simultaneously trying to find a balance between the firing rates (enough spikes must be generated to transmit the information to the network) and the sparsity of spiking networks (leveraging the promise of low energy consumption by neuromorphic devices). Following the same example, let us a fix a certain spiking amplitude and experiment with different firing rate scales to find this trade-off, considering the Loihi-tailored spiking ReLU activation. The neural activities of the same three neurons are shown in Figure 4a for a scale of 1 (i.e., keeping the same input weights as the original SNN), in Figure 4b for a scale of 5 (a linear scale of 5 is applied to the inputs of the neurons), and in Figure 4c using a scale of 50. As expected, the spikes fire much faster when increasing this parameter, with the spikes being almost indistinguishable in the latter case, thus reducing the sparsity of this network. Between the neural activities shown in Figure 4a,b, the mean firing rates are low (6 and 30.9 Hz) and show a more sparse firing rate, meaning that both are, in principle, better-suited to this application. The preliminary results computed within NengoDL (Figure 4d) indicate that a scale of 5 provides a slightly better performance, being the most adequate value for this parameter. Naturally, tuning these parameters is a harder task when dealing with more complex data and more complex spiking architectures, as we will see in the following section.

**Figure 4.** (**a**) Neural activities setting an amplitude = 0.01 and a firing rate scale = 1, (**b**) Neural activities setting an amplitude = 0.01 and a firing rate scale = 5, (**c**) Neural activities setting an amplitude = 0.01 and a firing rate scale = 50, and (**d**) predictions over the testing set.

#### *3.2. Load Forecasting*

Let us set another example using real data to calculate one-step ahead forecasts. In particular, short-term load forecasting is of interest due to its close relation with WPF, as both are necessary to operate and maintain the stability of the electrical grid [54]. Furthermore, load demand data show regular daily and weekly patterns, which are not observed in wind power data [55], so a model architecture formed of CNNs is a good candidate to extract such features [56]. Records of aggregated hourly demand data from Ireland can be found on the European Network of Transmission System Operators for Electricity (ENTSO-E) website [57]. The available measurements were recorded between 2016 and 2018.

As usual, we built and trained the rate-based equivalent of the model, and subsequently the activation functions were replaced. Then, the spike parameters were tuned without specifying any hardware requirements, and we monitored the initial results to choose the best values for these parameters. Some of these initial forecasts are shown in Figure 5. The existing patterns in load data were captured by the model, and adjusting the spikes' parameters is fairly straightforward. The dashed red line (obtained using an amplitude of 0.05 and a firing rate scale of 50) more closely matches the test data than the rest, so these values were chosen for its implementation on Loihi's emulator (or the hardware itself, if available).

**Figure 5.** Preliminary one-step-ahead load forecasts, setting different spike amplitudes and firing rates.

As indicated in Step 5, the network must be further adjusted to be run on Loihi. In our particular case, we must indicate what layers are run on- and off-chip, but other adjustments might be needed for more complex networks, such as distributing the connections of the network over multiple cores on Loihi [44]. Figure 6a shows that neurons are effectively firing in each layer, whereas Figure 6b compares the initial forecasts thaat were obtained previously while tuning the spike parameters (the red dashed line) and the load forecasts emulating the Loihi chip (dash-dot green line). It can be observed that the model architecture translates well to the emulator after fine-tuning those hardware specifications, resulting in similar load forecasts with respect to the initial evaluation of Step 4.

**Figure 6.** *Cont.*

**Figure 6.** Results for one-step ahead load forecasts: (**a**) Neural activities of 5 neurons of each layer. One input vector is shown over 1000 timesteps. (**b**) Predictions over the testing set with the SNN architecture (dashed red line) and running the SNN on the Loihi emulator (dash-dot green line).

#### *3.3. Wind Power Forecasting*

Noting the computational abilities of current neuromorphic devices, let us apply the proposed model architecture to build the spiking forecasting models for each mode extracted from Irish wind-power data after preprocessing the data with the VMD algorithm. The exact location of the wind farm is not disclosed for confidentiality reasons. The library Nengo [35] was used to simulate neuromorphic algorithms, together with the extensions NengoDL [37] for deep learning and NengoLoihi to emulate the behavior of Loihi hardware. The original neural-network-based models were implemented using Keras with Tensorflow backend for the rest of the models [38,58].

Following the proposed methodology, we first built the model (Step 1), and trained the rate-based neural network model to set its network parameters (Step 2). Then, we transformed it into an SNN by switching the activation functions to spiking ones (Step 3). In Step 4, we set empiric values for the amplitude and firing rate of the spikes (Table 1) within the NengoDL framework until we obtained a reasonable performance from the spiking model. The spiking amplitude modulates the amount of information transmitted to the subsequent layers of the network, whereas the firing rate adjusts how fast the spikes are being fired. If the firing rate is high, the behavior will be closer to the non-spiking model, and thus the performance will increase, but at the cost of losing the characteristic temporal sparsity provided by the spikes [59]. In addition, a high firing rate will lead to detrimental results on Loihi because of the discrepancy resulting from discretizing the spiking activation function (as shown in Figure 1). The low mean firing rates of these preliminary results (Table 2) suggest that the selected parameters are potentially good for implementation on Loihi. Afterwards, we configured some additional parameters to run the model on the Loihi emulator (Step 5). In particular, we must indicate what part of the model is run off-chip (in this case, the *off-chip* layer we use to communicate with the chip) and how long each input vector is presented to the network (in our case, we show each one for 0.4 s).

**Table 1.** Main spiking network parameters.

**Table 2.** Mean firing rates (Hz) for each layer.

The information recorded in Steps 4 and 5 is shown in Figures 7–14 for modes 1–8, respectively. Part a of these figures shows the neural activities of each layer (limited to 5 neurons for illustrative purposes). These neural activities correspond to the first input vector fed to the model, and produce the first point forecast, shown on part b. This constitutes one of the main differences in comparison to ANNs. Neurons in regular ANNs are static entities, which are always activated every time a new input arrives to the model, whereas neurons of SNN models will only be activated if certain dynamic conditions are met. From an user perspective, the neural activities help us visualize the mean firing rates shown in Table 2: modes 1, 3, and 8 exhibit higher firing rates, which translates into a large number of spikes being generated during this timeframe, whereas the rest of the modes present a more sparse behavior, resulting in a lower generation of spikes. In some cases, such as mode 4 (Figure 10) and mode 5 (Figure 11), the neurons of the *off-chip* layer need a long time to settle and thus start to spike, delaying the neural response of subsequent layers. Even if temporal sparsity is a desirable feature in a spiking model, in the sense that a smaller number of spikes means a lower consumption of energy (as a non-activated neuron will consume no energy), it might occasionally be advisable to finetune the firing rate of the *off-chip* layer to propagate the information faster to the rest of the network, as delays in the neural responses could degrade model performance. A quote on exact power consumption on the chip is erroneous in the current context, since the entire hardware is active during implementation, while only a minute fraction is actually used for the proposed problem. Under such circumstances, power quotes become relevant with sector-customised chips and a better handling of SNN architecture for DL, a direction in which the industry and current research is quickly moving.

**Figure 7.** Results for mode 1: (**a**) Neural activities of 5 neurons of each layer. One input vector is shown over 1000 timesteps. (**b**) Predictions over the testing set with the SNN architecture (dashed red line) and running the SNN on the Loihi emulator (dash-dot green line).

**Figure 8.** Results for mode 2: (**a**) Neural activities of 5 neurons of each layer. One input vector is shown over 1000 timesteps. (**b**) Predictions over the testing set with the SNN architecture (dashed red line) and running the SNN on the Loihi emulator (dash-dot green line).

**Figure 9.** Results for mode 3: (**a**) Neural activities of 5 neurons of each layer. One input vector is shown over 1000 timesteps. (**b**) Predictions over the testing set with the SNN architecture (dashed red line) and running the SNN on the Loihi emulator (dash-dot green line).

**Figure 10.** Results for mode 4: (**a**) Neural activities of 5 neurons of each layer. One input vector is shown over 1000 timesteps. (**b**) Predictions over the testing set with the SNN architecture (dashed red line) and running the SNN on the Loihi emulator (dash-dot green line).

**Figure 11.** Results for mode 5: (**a**) Neural activities of 5 neurons of each layer. One input vector is shown over 1000 timesteps. (**b**) Predictions over the testing set with the SNN architecture (dashed red line) and running the SNN on the Loihi emulator (dash-dot green line).

**Figure 12.** Results for mode 6: (**a**) Neural activities of 5 neurons of each layer. One input vector is shown over 1000 timesteps. (**b**) Predictions over the testing set with the SNN architecture (dashed red line) and running the SNN on the Loihi emulator (dash-dot green line).

**Figure 13.** Results for mode 7: (**a**) Neural activities of 5 neurons of each layer. One input vector is shown over 1000 timesteps. (**b**) Predictions over the testing set with the SNN architecture (dashed red line) and running the SNN on the Loihi emulator (dash-dot green line).

**Figure 14.** Results for mode 8: (**a**) Neural activities of 5 neurons of each layer. One input vector is shown over 1000 timesteps. (**b**) Predictions over the testing set with the SNN architecture (dashed red line) and running the SNN on the Loihi emulator (dash-dot green line).

At this stage, the performance of our models can finally be examined (Step 6). The model is designed to provide one-step-ahead point forecasts (Figure 15). The dashed-red lines show the forecasts obtained while tuning the model using the NengoDL framework in Step 4. While this preliminary model is able to forecast increasing/decreasing trends of power generation, it is not as accurate for high or low power-generation scenarios. Nonetheless, this initial assessment allows for us to prepare our model for Loihi (dashdot green line), which demonstrates the same skill in detecting increasing/decreasing trends of power generation as the preliminary model, while showing a better ability to forecast high/low power-generation values. This difference in performance also arises from the model architecture itself. When the model is initially evaluated outside the Loihi framework, it cannot discern that the first layer is only set to start to generate spikes. Such nuance is captured when the model is configured for implementation on Loihi. Additionally, we observe that the forecasts are not as accurate compared to a non-spiking VMD-GRU model that we used in a previous study using the same data [50], and the outputs are generally noisier. However, this is an expected outcome due to the current limitations of neuromorphic hardware.

**Figure 15.** One-step ahead WPFs with the SNN architecture (dashed red line), running the SNN on the Loihi emulator (dash-dot green line) and a non-spiking VMD-GRU model (purple crosses) over the testing set.

In conclusion, we have successfully transformed a non-spiking neural model into a spiking one with a reasonably good performance, having achieved a 2.84% NMAE for onestep ahead forecasts with the model being adapted to neuromorphic hardware. This type of proof is not only helpful to prove that industrial applications such as WPF modelling can be transferred to non-von Neumann architectures such as neuromorphic computing, but to provide guidelines to the manufacturers of such hardware (e.g., Intel and the development of devices as Loihi) to cater to the industry's needs.

#### **4. Conclusions**

Neuromorphic computing provides a new paradigm to build energy-efficient, lowlatency algorithms in contrast to the current state-of-the-art ML/DL strategies, thus potentially reducing the computational cost of training and deploying artificial-intelligence-based models. In particular, SNNs aim to learn in a more biologically plausible manner [60] by more closely mimicking the spike-based transmission of information that occurs in the brain [61]. At present, the two major challenges for the use and implementation of SNNs are (1) the training of such models, as the well-established training strategies based on the backpropagation algorithm applied to ML/DL cannot be directly used, as spikes are not differentiable, and (2) the implementation of SNNs on neuromorphic hardware, as SNNs must be tailored to cater to the specific requirements of the hardware. The first challenge has been addressed with different approaches to date, such as ANN-to-SNN conversion, and using variations in error backpropagation to directly train SNNs. The second challenge is hardware-dependent, and should be addressed according to the requisites of the hardware used to implement the SNN. Additionally, there is currently a lack of studies applying neuromorphic computing for practical cases that are useful for both research and industrial practices, such as the design of WPF models [62].

In this paper, we adopt an ANN-to-SNN conversion approach to forecast wind power, and obtain these WPFs emulating or running the spiking model using the neuromorphic

hardware Loihi [41]. SNNs are designed using the framework provided by the software Nengo [35,37]. First, we build and train the non-spiking neural network. After training, we map the parameters and replace the activation functions for their spiking counterparts, which will be used during the prediction stage. Then, without considering hardware specific constraints, some preliminary results are evaluated to tune some spike-related parameters such as the firing rate or the amplitude of the spikes. Finally, the SNN is further adjusted to be run on the hardware emulator (or actually running the model on Loihi if available) to obtain the WPFs. Following all these steps, we managed to reach our goal of achieving a good level of performance with the proposed spiking architecture, obtaining a NMAE of 2.84% for one-step ahead forecasts when the model is emulated on Loihi.

As neuromorphic computing is not a well-established technology at present, there is room for future research. First, the proposed ANN-to-SNN conversion approach for short-term WPF can be further refined by tuning the firing rates of each layer individually and considering the use of synaptic filters to smooth the output. Second, the modelling of spiking neural networks can be improved by (1) directly training the network to increase the efficiency of neuromorphic devices from the energy point of view and (2) using online approximations of the backpropagation algorithm to adjust network parameters with the arrival of new data. Third, it will be possible to implement more complex biological inspired neural networks in the near future, as the computational capability of neuromorphic computing continues to increase with the development of new devices, such as Loihi 2 [63]. As reducing the computational cost is one of the main reasons to use neuromorphic devices, any future line of research must address model performance in terms of energy consumption to verify that it achieves a significant reduction in consumption compared to conventional computer architectures. To achieve that, the neural spiking models must not only be emulated but run on a real neuromorphic device to realistically measure this feature.

**Author Contributions:** Conceptualization, J.M.G.S., V.P. and B.G.; methodology, J.M.G.S. and B.G.; software, J.M.G.S.; validation, J.M.G.S., V.P. and B.G.; formal analysis, J.M.G.S.; investigation, J.M.G.S., V.P. and B.G.; resources, J.M.G.S., V.P. and B.G.; data curation, J.M.G.S., V.P. and B.G.; writing original draft preparation, J.M.G.S. and B.G.; writing—review and editing, J.M.G.S., V.P. and B.G.; visualization, J.M.G.S.; supervision, V.P. and B.G.; project administration, V.P. and B.G.; funding acquisition, V.P. and B.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** The authors acknowledge the funding of SEAI WindPearl Project 18/RDD/263.

**Data Availability Statement:** Relevant computed data are available from authors upon reasonable request to the corresponding author.

**Acknowledgments:** The authors like to thank George Vathakkattil Joseph and Aasifa Rounak for their technical support. Bidisha Ghosh would like to acknowledge the support of ENABLE (Grant number 16/SP/3804) and Connect Center (Grant number 13/RC/2077\_P2). Vikram Pakrashi would like to acknowledge the support of SFI MaREI centre (Grant number RC2302\_2), the resources and support of Intel Neuromorphic Research Community, Accenture NeuroSHM project with UCD, Science Foundation Ireland NexSys 21/SPP/3756, and UCD Energy Institute.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**

The following abbreviations are used in this manuscript:

#### **References**

### *Article* **A Comprehensive Study of Random Forest for Short-Term Load Forecasting**

**Grzegorz Dudek**

Department of Electrical Engineering, Cz ˛estochowa University of Technology, 42-200 Cz ˛estochowa, Poland<EMAIL_ADDRESS>
**Abstract:** Random forest (RF) is one of the most popular machine learning (ML) models used for both classification and regression problems. As an ensemble model, it demonstrates high predictive accuracy and low variance, while being easy to learn and optimize. In this study, we use RF for short-term load forecasting (STLF), focusing on data representation and training modes. We consider seven methods of defining input patterns and three training modes: local, global and extended global. We also investigate key RF hyperparameters to learn about their optimal settings. The experimental part of the work demonstrates on four STLF problems that our model, in its optimal variant, can outperform both statistical and ML models, providing the most accurate forecasts.

**Keywords:** random forest; regression tree; pattern representation of time series; short-term load forecasting

### **1. Introduction**

Electricity demand forecasting is extremely important for energy providers to ensure the secure, effective and economic operation of the power system. Short-term load forecasting (STLF) covers a forecast horizon of a few hours to a few days. STLF is necessary for generation resource planning to meet electricity demands and optimize the power flow on the transmission grid to avoid overloads. As electricity demand is a major driver of electricity prices, load forecasting plays a key role in competitive energy markets. The STLF accuracy directly affects the financial performance of energy market participants.

The importance of accurate electricity demand forecasts for the safe, reliable and effective operation of power systems is behind the great interest of researchers in this area. STLF problems are complex because electricity demand time series express a nonlinear trend, multiple seasonality, variable variance, significant random disruptions and changing daily profile. These challenging factors place high demands on STLF models.

### *1.1. Related Work*

Roughly, STLF methods can be divided into statistical and ML methods. The most popular representatives of the first group are: auto-regressive integrated moving average (ARIMA) [1], exponential smoothing (ETS) [2], linear regression [3], and Kalman filtering [4]. The main drawbacks of the statistical methods are their linear character, limited adaptability, limited ability to deal with complex seasonal patterns, and problems with capturing longterm dependencies in time series and introducing exogenous variables into the model [5].

ML models provide more flexibility in modeling nonlinear functions. Unlike statistical methods, they do not require strong assumptions about the mapping function, and they learn relationships between predictors and targets directly from historical data. Among ML methods for forecasting, neural networks (NNs) have gained the most popularity in recent years [6]. The multitude of architectural solutions and mechanisms to improve performance encourage the use of NNs to solve complex forecasting problems such as STLF. Classical NNs were investigated for suitability for STLF in [7]. To deal with triple seasonality in time

**Citation:** Dudek, G. A Comprehensive Study of Random Forest for Short-Term Load Forecasting. *Energies* **2022**, *15*, 7547. https://doi.org/10.3390/en15207547

Academic Editor: Andrzej Bielecki

Received: 15 September 2022 Accepted: 10 October 2022 Published: 13 October 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

series, patterns of the daily profiles were introduced, which filter out the trend and the weekly and yearly seasonality (a similar approach was used in this study, where different definitions of patterns are examined). Among the considered NN architectures, which included Multilayer Perceptron (MLP), Radial Basis Function NN, General Regression NN (GRNN), Fuzzy Counterpropagation NN, and Self-Organizing Maps, GRNN and MLP stand out as performing with the highest accuracy.

In recent years, a development in the NN field has been the move towards deep and, especially useful in forecasting, recurrent architectures [6]. Deep NNs (DNNs) are especially beneficial in learning the most useful data representation for modeling a given target function, while Recurrent NNs (RNNs) are beneficial in modeling complex, shortand long-term temporal relationships in data. New mechanisms and procedures introduced to RNNs such as delayed connections, attention, hybrid architecture, dynamic training sets, residual connections and flexible loss functions improve their learning capabilities and expressive power to solve forecasting problems [5]. Some examples of using DNNs and RNNs for STLF are [8], where Convolutional NNs (CNNs) are utilized to extract load and temperature features, which are fed as inputs into the bidirectional propagating RNN to perform hourly electrical load forecasting [9], where a recurrent inception CNN is proposed for STLF that combines RNN and 1-dimensional CNN [10], where RNN with attention significantly reduced forecasting errors as compared to the current state-of-the-art results; and [11], where deep residual networks integrate domain knowledge and researchers' understanding of the problem and enables probabilistic load forecasting using Monte Carlo dropout.

An effective way to increase the forecast accuracy and robustness of both statistical and ML models is ensembling. This combines multiple models for a common response to improve both the accuracy and stability of the final solution compared to a single model [12]. The theoretical properties of forecast combination investigated in [13] answer the question why a simple average of forecasts often outperforms forecasts from single models. They also prove that simple averages in many cases perform better than more complicated weighting schemes. The beneficial effects of the forecast aggregation on STLF accuracy are shown in many papers: in [14], several ML methods are aggregated in ensembles for one-day-ahead wind power forecasting; in [15], to forecast an interval-valued load, ensemble of RNNs is applied, which learns on the components of the bivariate empirical mode decomposition; in [5], ensembling of a hybrid model, which combines ETS and RNN, leads to a significant reduction in the forecast error; in [16], an ensemble of randomized DNN combined with a walk-forward decomposition is proposed; in [17], a stacking ensemble approach is used to combine DNNs, and in [18], several methods of aggregating base models (MLPs) are considered. Stacking, used in the last two papers, is a way of combining base models via meta-learning, i.e., a meta-model is trained on the predictions of the base models.

Alternatives to stacking are boosting and bagging. Popular representatives of these are gradient-boosted trees and random forest (RF), respectively. Both, used as forecasting models, are based on regression trees. It was shown in [19] that RF can compete with both classical models and NNs in STLF. It can deal with complex time series using appropriate data preprocessing, which produces normalized patterns of the daily profiles. Based on this research, in this study, to improve RF performance, we extend pattern definitions and introduce additional predictors. To enrich input information, in [20], the input patterns are extracted from electrical, meteorological and calendar data by temporal CNN. Fed with these patterns, a Light Gradient Boosting Machine, a type of gradient-boosted trees algorithm, was able to forecast very volatile industrial customer loads. A novel tree-based ensemble method called Warm-start Gradient Tree Boosting (WGTB) was proposed in [21]. It combines four different inference models and aggregates their outputs by a warm-start, bagging and boosting, which at the same time reduces bias and variance. The result proves the efficiency of the proposed strategy and shows an improvement in STLF accuracy over baseline models. Another type of tree-based ensemble, eXtreme Gradient Boosting (XGBoost), was used in [22] for forecasting electricity consumption by industrial customers. To deal with multiple seasonality, the time series were first decomposed using variational mode decomposition. Then, a linear regression model was applied for the trend series and a XGBoost regression model was applied for each fluctuation sub-series.

To close this section, we note that, according to some studies, tree-based ensembles are not inferior to NNs in terms of forecast accuracy. In [23], the authors have examined and reproduced a number of state-of-the-art DNNs for time series forecasting. DNNs were compared on different datasets to a Gradient Boosting Regression Tree (GBRT). The experimental results show that a conceptually simpler model such as GBRT can compete and sometimes outperform modern DNNs by efficiently feature-engineering the input and output structures of GBRT.

#### *1.2. Motivation and Contribution*

Tree-based methods are widely used as prediction models as they have very attractive properties such as a capacity for flexible nonlinear regression, which can capture complex interactions between variables and effectively handle multiple predictors (including exogenous ones) of various types (numeric, binary, and categorical). Moreover, they are robust against over-fitting of the training data, they are relatively simple to tune, and they are easy to implement with the available software. Their effectiveness has been confirmed in many forecasting competitions, for example those carried out on the Kaggle platform [24].

The excellent performance of tree-based approaches was demonstrated in the 2020 M5 forecasting competition. The top places in this competition, in terms of both accuracy and uncertainty, were dominated by entries that used tree-based ML methods such as gradientboosted trees [25]. Four out of the five winning models used a variant of the tree-based method and most of the other top 50 best-performing models adopted similar approaches to the winning submission by training recursive and non-recursive tree-based models [26]. Thus, tree-based forecasting models appear to be strong competitors to NNs, which in the form of deep learning-based models dominate the recent literature on forecasting methods [6].

In this study, motivated by the excellent results of tree-based models in forecasting competitions, we apply RF to the challenging problem of STLF. RF gives similar results to boosting, but is easier to train and tune [27]. The main contribution of this work is to examine RF models using a variety of time series preprocessing methods and training modes. In the local mode, RF learns on samples similar to the query sample, which enables the model to focus on the local features of the target function around the query pattern and improve accuracy in this region. In the global mode, to achieve the same goal, i.e., focusing on the proper region of the target function, we introduce additional calendar variables. By examining different methods of time series preprocessing, we find the most useful data representation for achieving the highest accuracy of the model. We empirically demonstrate that the proposed approach outperforms in terms of accuracy both standard statistical models as well as more sophisticated ML approaches.

The novelty of this work in relation to our previous work [19] is twofold. First, we extend pattern definition by introducing seven types of patterns based on the historical data. They incorporate daily and/or weekly seasonality, while in [19], the patterns captured only daily seasonality. Second, we introduce a global mode of training with additional predictors representing calendar data. In [19], only a local training was considered without calendar inputs.

The rest of the paper is organized as follows. In Section 2, we propose several data prereprocessing methods for electricity demand times series. Section 3 defines the forecasting problem and RF training modes. Section 4 describes the RF algorithm in application to STLF. The experimental framework used to evaluate the performance of the proposed model and compare it with baseline models is described in Section 5. Finally, Section 6 concludes the work.

#### **2. Data Preprocessing**

The power system load or electricity demand time series express a trend, triple seasonality (annual, weekly and daily) and random fluctuations. These components are dependent on the system size, size of the economy served, customer structure, as well as weather and climatic conditions. The daily load profiles that we focus on in STLF vary throughout the year and depend on the day of the week [5].

To forecast future demand with the least possible error, the forecasting model should be fed by the most relevant predictors. In our univariate STLF model, which produces forecasts for the next day, the predictors are selected from recent history and are preprocessed accordingly. The forecasting model based on RF is fed by input patterns **x***i* and produces encoded forecasts for hour *t* of day *i*, *yi*,*t* (MISO model).

Let {*zτ*}*M τ*=1 be an electricity demand time series with hourly resolution and vector **z***i* = [*zi*,1, ..., *zi*,24] represents its 24-hour-long sequence for day *i*. To capture the characteristic properties of the series, remove the trend and unify data, we define the input patterns as follows:

**r1** The input patterns are defined based on the weekly sequence which precedes forecasted day *i*:

$$\mathbf{x}\_{i} = \frac{\mathbf{s}\_{i} - \overline{\mathbf{s}}\_{i}}{||\mathbf{s}\_{i} - \overline{\mathbf{s}}\_{i}||} \tag{1}$$

where **x***i* ∈ R168 is the input pattern, **s***i* = [**z***i*−7, ..., **z***i*−1] is the demand sequence of the week preceding the forecasted day *i* and *si* is the mean of this sequence.

Input vectors (1), which for successive *i* represent overlapping weekly sequences shifted by one day, are normalized versions of centered vectors **s***i*. They all have zero mean, the same variance and the same unity length. However, they differ in shape. Thus, we assume that the weekly shape carries the information about the forecasted demand of the day following this week.

Figure 1 shows the sequences which are used for x-patterns construction and Figure 2 shows data used for construction patterns r6 and r7. Depending on the definition, x-patterns introduce different input information to the model. Pattern r1 introduces detailed information about the weekly sequence which precedes the forecasted day. Note that r1 expresses both daily and weekly seasonality unlike r2, which carries information only about the daily seasonality. To deal with weekly seasonality when r2-patterns are used, the model can be trained in the local mode, i.e., on the subset of x-patterns corresponding to the forecasting task (see Section 3).

**Figure 1.** Load time series points used for input patterns construction.

**Figure 2.** Cross-patterns: r6 (green + blue) and r7 (green + orange + blue).

Pattern r3 introduces information on the demand in the previous seven days at the same hour as the forecasted one. It expresses only weekly seasonality. Information about the daily seasonality is not included, so the local mode of training, i.e., training on the selected r3-patterns corresponding to the forecasting task (see Section 3), can help with dealing with daily seasonality. Similar information as in r3 is contained in pattern r4 but from a longer 3-week period. Pattern r5 shows neither daily nor weekly seasonality. It carries information about the demand at the same hours as forecasted in previous days of the same type as the forecasted day.

Cross-patterns r6 and r7 express both daily and weekly seasonalities as r1, but in a more sparing form, using respectively 30 or 44 instead of 168 components. In [28], we showed that STLF based on both daily and weekly patterns gives better results than forecasting based on separate daily or weekly patterns. In [28], we aggregated forecasts generated by two neural models: a daily pattern-based model and a weekly pattern-based model, while in this study, we combine daily and weekly patterns into one pattern and use only one model.

Examples of input patterns r1–r7 are depicted in Figure 3. Note different shapes of patterns, carrying different input information.

**Figure 3.** Examples of input patterns r1–r7.

The output data, i.e., the electricity demand at hour *t* of day *i*, is encoded as follows:

$$y\_{i,t} = \frac{z\_{i,t} - \overline{s}\_i}{||s\_i - \overline{s}\_i||} \tag{2}$$

where **s***i* is the demand sequence preceding forecasted day *i*, defined depending on the x-pattern type r1–r7 (this is the same sequence based on which pattern **x***i* was defined) and *si* is the mean of this sequence.

The output data are encoded similarly to the input data, using the same coding variables: *si* and **s***i* − *si*. Thus, Equation (2) like Equation (1) filters the data by removing the local trend (*si*) and unifying the variance (this is a function of the denominator of these equations, which can be thought of as a measure of diversity of the input sequence). Such filtered and unified data are predicted by the forecasting model (RF). Then, the real forecast is determined from transformed Equation (2):

$$\mathcal{Z}\_{i,t} = \mathcal{Y}\_{i,t} || \mathbf{s}\_i - \overline{\mathbf{s}}\_i || + \overline{\mathbf{s}}\_i \tag{3}$$

where *y*ˆ*i*,*t* is the model prediction and *z*ˆ*i*,*t* is the real forecast.

Note that (3) brings back the local current properties of the time series (level and dispersion), which were removed by (1) and (2) to simplify the relationships between input and output data. We have successfully used this kind of preprocessing of input and output data in our previous load forecasting models to deal with multiple seasonality, simplify the model and speed up training, see, e.g., [7,19,28–32].

#### **3. Forecasting Problem and Training Modes**

The forecasting task is defined as follows: predict electricity demand for hour *t* ∗ (1, ..., 24) of day *i* ∗ based on historical data. Day *i* ∗ represents day of the week *d*∗ (*Monday*, ..., *Sunday*). To maximize the forecasting performance and to make the most of all available training data up to day *i* ∗ (forecasted day), the forecasting model is trained individually for each forecasting task and it performs only one prediction: *y*ˆ*i*∗,*t*∗ . Note that the "global" generalization property of the model is not important because it is built to make only one prediction. What is important is the "local" performance in the neighborhood of pattern **x***i*∗ . To increase this property, we use two approaches. In the first one, we train the model in the local mode and in the second one, we extend input patterns with calendar variables when we use global learning. For comparison we also train the model in the standard global mode.

The full training set determined on the historical data is Ψ = {(**x***i*, *yi*,*t*)}, where *i* = 1, ..., *i* ∗ − 1, *t* = 1, ..., 24 and pair (**x***i*, *yi*,*t*) includes the input pattern and target defined according to r1–r7. The three training modes are as follows:

**Local** The model is trained on the subset of Ψ containing pairs (**x***i*, *yi*,*t*) which correspond to the forecasting task, i.e., the pairs which include targets representing the same day type as forecasted, *d*(*yi*,*t*) = *d*∗, and the same hour as forecasted, *t* = *t* ∗.

**Global** The model is trained on full training set Ψ.

**Global extended** The input data are extended with calendar information:

• season of the year encoded as follows [29]:

$$\mathbf{p}\_{i} = \left[ \sin \frac{2\pi \mathfrak{H}i}{366}, \cos \frac{2\pi \mathfrak{H}i}{366} \right] \tag{4}$$

where #*i* is the number of day *i* in the year,

The training set in the global extended mode is of the form: Ψ = {(**x***i*, **p***i*, *di*, *t* , *yi*,*t*)}, *i* = 1, ..., *i* ∗ − 1, *t* = 1, ..., 24.

In the local training mode, the model solves the forecasting task by learning on the samples expressing similar properties and relationships between input and output data as those expressed by input pattern **x***i*∗ and forecasted value *yi*∗,*t*∗ . That is, the training input patterns are limited to those that are similar in shape to pattern **x***i*∗ (further limitations in this regard can be made by selecting training data from the same period of the year as day *i* ∗ or by selecting training data based on similarity to pattern **x***i*∗ [7], but we have not employed these approaches in this study). The relationships between input and output data expressed in the local training set are limited to those corresponding to day type *d*∗ and hour *t* ∗, so we expect the model to more accurately approximate the relationship between **x***i*∗ and *yi*∗,*t*∗ than in case of global training on the full training set expressing the relationship for all day types and hours.

In the global training mode with extended inputs, the model has additional input information to more accurately solve the forecasting task, i.e., the calendar variables. Although the model is global, we expect that the calendar data will help to increase its local accuracy around pattern **x***i*∗ . A regression tree as a base forecasting model can more appropriately divide the input space using the calendar variables than without these variables, and therefore approximate locally the target function with greater accuracy. Although this will lead to a more complex model than in the other two training modes.

#### **4. Random Forest for STLF**

RF is a ensemble learning algorithm based on decision trees (CART [33]) as the base models [34]. It is suitable for either regression or classification problems. In this study, for forecasting problems, we focus on the regression RF based on regression trees.

RF is devoid of the well-known drawbacks of single trees such as unstable splits and a lack of smoothness [27]. It combines bagging [35] with a random subspace method [36]. The key idea in bagging is to average multiple noisy but approximately unbiased base models and thus reduce the variance. Trees as noisy and low biased models if they have grown sufficiently deep, are great candidates for bagging. The main goal of the random subspace method is to increase diversity between trees by restricting them to work on different random subsets of the full predictor space (more specifically, at each node of the tree, a random predictor subset is selected). Each tree in the forest is built from a bootstrap sample of the original dataset, which is an additional source of diversity. Random predictors selected in the nodes of bagged trees help to decorrelate the trees and improve prediction accuracy as well as reduce the model variance.

The RF algorithm draws a bootstrap sample Ψ*k* of size *N* from training set Ψ for each of *K* trees, *k* = 1, ..., *K*. For each bootstrapped sample, a tree *T* is grown by recursive partitioning the input space in each node until a minimum leaf size is reached. At each node, data splits based on *p* out of *n* predictors chosen at random are considered. The best split is determined by maximizing the reduction in mean squared error (MSE) over all splitting candidates and cutpoints. After all *K* trees are grown in this fashion, the RF predictor is [27]:

$$\hat{f}\_K(\mathbf{x}) = \frac{1}{K} \sum\_{k=1}^K T(\mathbf{x}; \Theta\_k) \tag{5}$$

where **x** is the input pattern and Θ*k* characterizes the *k*-th tree in terms of split predictors, cutpoints, and terminal-node values.

RF has the following hyperparameters:

In practice the best values for the hyperparameters are dependent on the problem, and they should be treated as tuning parameters.

The standard method of selecting split predictors [33] has two drawbacks. Firstly, it tends to miss important interactions between pairs of predictors and the response. Secondly, it tends to select continuous predictors that have many levels, which masks more important predictors that have fewer levels, such as categorical predictors. To mitigate selection bias and increase detection of important interactions, curvature or interaction tests can be applied [37,38]. Therefore, in this study we consider three methods of selecting the split predictors:

The algorithm of RF construction for STLF is shown in Algorithm 1. It produces a set of *K* trees, {*Tk*}*K k*=1. Based on them, to make a prediction for new point **x**, we use (5). Then, the real forecast is calculated from (3). Note that training set Ψ is prepared for the selected training mode and input pattern type.

#### **Algorithm 1** Random forest construction for STLF

**Input:** training set Ψ containing *N* samples, number of trees *K*, minimum leaf size *m*, number of predictors to select at random for each split *p*, split predictor selection method *s* **Output:** set of trees {*Tk*}*K k*=1

#### **Procedure:**

**for** *k* = 1 **to** *K* **do**

Draw a bootstrap sample Ψ*k* of size *N* from Ψ

Grow tree *Tk* to Ψ*k* by recursively repeating the following steps for each terminal node, until the minimum node size *m* is reached:

– Select *p* predictors from the *n* predictors

In the experimental study (Section 5), we use RF specified in Algorithm 1 in several variants depending on the data preprocessing method (r1–r7) and training modes, i.e., local, global and global extended. In the global extended mode, the predictor vector is composed of **x***i*, **p***i*, *di*, *t* . In the other training modes, the predictor vector is the same as input pattern **x** (1).

#### **5. Simulation Study**

In this section, we investigate RF variants with different data preprocessing methods and training modes. We compare RF performance with that of other models based on classical statistical methods and ML methods.

STLF for four countries is performed: Poland (PL), Great Britain (GB), France (FR) and Germany (DE). The real-world data was collected from ENTSO-E repository (www.entsoe. eu/data/power-stats; accessed on 6 April 2016). It details the hourly power system load in the period from 2012 to 2015. The last year of the data (2015) is treated as a test period. We predict the daily load profiles for each day of this period, excluding atypical days such as public holidays (between 10 and 20 days a year). RF models were optimized on the data from 2012 to 2014, with validation data composed of 100 patterns selected randomly from 2014 and training data preceding the validation pattern.

#### *5.1. Results for Different Preprocessing Methods and Training Modes*

Tables 1 and 2 show mean absolute percentage error (MAPE) and root mean square error (RMSE), respectively, for input patterns r1–r7 and different training modes. Figures 4–6 show the boxplots of MAPE. The results can be summarised as follows:

**Table 1.** Validation MAPE for different input patterns r1–r7 and training modes (lowest errors in bold, second lowest errors in italics).

**Table 2.** Validation RMSE for different input patterns r1–r7 and training modes (lowest errors in bold, second lowest errors in italics).

**Figure 4.** Local training mode: Boxplots of validation MAPE for different input patterns r1–r7.

**Figure 5.** Global training mode: Boxplots of validation MAPE for different input patterns r1–r7.

**Figure 6.** Global extended training mode: Boxplots of validation MAPE for different input patterns r1–r7.

Based on the results, the recommended training mode is global extended with r4 patterns for PL, FR and DE, and r6 patterns for GB. These variants of RF were used in the experiments described in the next sections.

#### *5.2. Tuning Hyperparameters*

In this experiment, we change the selected hyperparameter in the range shown in Figure 7, keeping the remaining hyperparameters at their constant values as follows: number of trees in the forest—*K* = 100, minimum number of leaf node observations *m* = 1, and number of predictors to select at random for each decision split—*p* = *n*/3.

Figure 7 shows the impact of hyperparameters on the forecasting error (MAPE). As expected, the error decreases with the number of trees in the forest. The reduction in MAPE when the RF size changes from 1 do 300 trees was from 38.7% for PL to 50.0% for DE. At the same time, a significant reduction in the forecast variance was also observed from 63.8% for PL to 82.5% for DE. It can be seen from Figure 7b that an increase in the minimum leaf size leads to a deterioration in the results. Small values of *m*, close to 1, are preferred. This means that trees as deep as possible are the most beneficial in RF. The optimal number of predictors selected in the nodes to perform a split varies from country to country (see Figure 7c). For PL and DE it is 15, for GB it is 20, and for FR it is 6. These values differ from the recommended *p* = *n*/3, which are 8 for PL, FR and DE, and 11 for GB.

**Figure 7.** Validation MAPE depending on hyperparameters: number of trees in the forest (**a**), minimum number of leaf node observations (**b**) and number of predictors to select at random for each decision split (**c**).

Using optimal values of the hyperparameters for each country, we investigate the methods of split predictor selection s1–s3. Table 3 shows the results, validation MAPE and RMSE. Both accuracy measures show similar results for all methods of split predictor selection. Therefore, s1 is recommended as a simple, standard method, which does not cause any additional computational burden.

**Table 3.** Validation MAPE and RMSE for different methods of split predictor selection (training mode: global extended, #trees: 300, minimum leaf size: 1; lowest errors in bold).

Figure 8 shows the "importance" or "predictive strength" of the predictors estimated on the out-of-bag data (this is discussed further in Section 5.4). As can be seen from this figure, when r4 extended pattern is used (PL, FR and DE), the most important predictor is the last component of the r-pattern, i.e., the predictor expressing electricity demand at forecasted hour *t* of the day preceding the forecasted day, *zi*−1,*t*. The importance of this predictor reaches 3.5 for FR and over 5 for PL and DE, while the importance of other demand predictors is usually below 2. Among the calendar predictors, the most important for r4 extended pattern are those coding the season of the year, especially for DE. For cross-pattern r6 (GB), the most important predictors are the calendar ones: day of the week and season of the year (*p*1). The next positions are occupied by predictors coding the demand for the last four hours of the day before the forecasted day (*zi*−1,24 is clearly the most important of these) and predictor coding demand at forecasted hour *t* week ago, *zi*−7,*t*. Note the low importance of the other predictors representing demand at hour *t* of the preceding days, *zi*−6,*t*–*zi*−2,*t*.

**Figure 8.** Importance of the predictors (*x*1 − *x*30—predictors expressing demand pattern r4 or r6, *p*1 and *p*2—components of vector **p** expressing season of the year, *d*—day of the week, and *t*—hour of the day).

Table 4 shows the forecasting results for the test set when using RF with the optimal values of hyperparameters. As performance metrics we use: MAPE, MdAPE (median of absolute percentage error), IqrAPE (interquartile range of APE), RMSE, MPE (mean percentage error), and StdPE (standard deviation of PE). MdAPE measures the mean error without the influence of outliers, while RMSE, as a square error, is especially sensitive to outliers.

**Table 4.** Results for test data (training mode: global extended, #trees: 300, minimum leaf size: 1, split predictor selection method: s1).

The MAPE and MdAPE values in Table 4 indicate that the most accurate forecasts were obtained for PL and DE, while the least accurate were for GB. MPE allows us to assess the forecast bias. Positive values of MPE indicate underprediction, while its negative values indicate overprediction. Note that for PL and DE the forecast bias was significantly smaller than for GB and FR. The same can be said about the forecasts dispersion measured by IqrAPE and StdPE.

#### *5.3. Results Comparison with Other Models*

We compare the performances of RF with other models including statistical models and ML models. The comparative models are outlined below (see [30,39] for further description). Their hyperparameters were selected on the data from 2012–2014 in grid search procedures using a variant of cross-validation or selected by experimentation (this applies to models with a large number of hyperparameters, which are difficult to optimize using standard methods due to the huge search space).

We also compare our model with competitive tree-based ensembles: XGBoost [40] and LightGBM [41]. Their predictors include both calendar data (hour of the day, day of the week, quarter, month, year, day of the year, day of the month and week of the year) and historical demands (demands at hour *t* of 21 consecutive days preceding the forecasted day).

Table 5 compares MAPE for RF and the baseline models. From this table, you can clearly see the better performance of RF compared to the other models. RF outperformed all other models in terms of accuracy for PL, GB and DE. For FR it took third place after RandNN and SVM. To confirm the results, a pairwise one-sided Giacomini-White test was performed (GM test) [42]. Its results, *p*-values, are shown in Figure 9 (we used GW test implementation from [43]). Small *p*-values, below 0.05 (green color), indicate that the model on the *X*-axis significantly outperforms in terms of accuracy the model on the *Y*-axis.

**Table 5.** MAPE comparison between RF and baseline models (lowest errors in bold).

**Figure 9.** Results of the Giacomini-White test for the proposed and baseline models (black color is for *p*-values larger than 0.10).

#### *5.4. Discussion*

In our previous work [19], we used a local training mode with input patterns r2, which express daily profiles. Our current research revealed that input patterns incorporating weekly seasonality (r4) or both daily and weekly seasonality (r6) combined with global training with extended inputs improve the results (note that in Tables 1 and 2, patterns r4 and r6 provide lower errors than r2 for all countries). Calendar data used as additional input in the global extended training helps the trees to properly partition the input space and thus approximate the target function with greater accuracy. It does not take place without costs: the complexity of the model increases due to learning on all data, not just the selected data as in local training.

Table 5 and Figure 9 show that the proposed RF outperforms classical statistical models (ARIMA and ETS), modern statistical model (Prophet), classical ML models (MLP, SVM, ANFIS, GRNN), modern ML models (LSTM, RandNN), similarity-based models (FNM, N-WE) as well as state-of-the-art ML models (MTGNN, ES-adRNNe) and boosted regression trees (XGBoost, LightGBM). The last two models, as well as the proposed RF model, also used calendar variables, even in larger numbers. However, in contrast to these models, our model uses specific time series preprocessing, which may be a decisive advantage. Our model also outperforms ES-adRNNe, which is a very sophisticated and complex model developed especially for STLF [39]. To increase its predictive power, it is equipped with a new type of RNN cell with delayed connections and inherent attention, it processes time series adaptively, learning their representation and it learns in the cross-learning mode (i.e., it learns from many time series in the same time). It reveals its strength with a large amount of data, numerous and long time series. In our case, this condition was not met—there were only four, relatively short series available for training the models. In this case, the proposed RF model, which learns from individual series, generated more accurate forecasts than ES-adRNNe.

It is worth noting that RF has few hyperparameters to tune, which makes it easy to optimize (compare with DNNs with many hyperparameters). The results of our experiments confirmed that the number of trees in the forest should be as large as possible and the mimimum leaf size can be set to one. Therefore, the key hyperparameter remains the number of predictors to sample in the nodes. Its optimal values significantly differed from

the recommended default values. Our attempt to increase the performance of the RF model through alternative methods of selecting predictors for split failed. Neither the curvature test nor the interaction test, which take into account the relationship between predictors and response when splitting data in nodes, improves the results significantly over the default CART method.

In our study, we used both continuous and categorical predictors. Such a mix causes many difficulties for other models such as ARIMA, ETS, NNs, SVM, LSTM and others. Categorical variables cannot be processed by these models directly. Such predictors must be converted into numerical data, so as to maintain the relationship between their values. The method of this conversion can be treated as an additional hyperparameter. RF has no problems with categorical variables, which is its big advantage. Moreover, RF can deal easily with any number of additional exogenous predictors and does not need to unify predictor ranges, which is often necessary for other models. RF can even deal with raw data because the predictors are not processed by the tree in any way, just selected in nodes, to construct a specific decision model (flowchart-like structure).

Regression tree provides fast one-pass training which does not need to repeatedly refer to the data. In contrast, NNs, which use a variant of the gradient descent optimization algorithm with multiple scanning of a dataset, are more time-consuming to train. Additionally due to the number of hyperparameters, they are also much more expensive to optimize in terms of time than RF. The training of a tree does not provide an optimal result because decisions about data split are made in nodes using a local rather than a global criterion, i.e., the split made may not be optimal from the point of view of the final result. However, the NN learning process also does not lead to optimal results due to sensitivity to the starting point and tendency to fall into the traps of the local minimum. Note that non-optimality of the trees is mitigated by their aggregation in the forest. Aggregation also smoothes out functions modeled by individual trees and reduces their variance. The learning process of RF can be easily paralleled because the individual trees learn independently.

One useful feature of RF is that it enables the generalization error to be estimated using out-of-bag (OOB) patterns, i.e., training patterns not selected for the bootstrap sample (approximately one third of the training patterns are left out in each bootstrap sample). Therefore, the time-consuming cross-validation that is widely used in other models for estimating the generalization error is not needed. Using OOB patterns, the generalization error can be estimated during one training session, along the way. Although for forecasting problems, where training patterns should precede validation to prevent data leakage, the OOB approach as well as standard cross-validation may be questionable. For this reason, we did not use the OOB approach in this study. Instead, we applied a different strategy. We chose a set of 100 validation patterns from 2014 and for each of them we trained RF on training patterns preceding the validation pattern.

A valuable feature of RF is its built-in mechanism for predictor selection. In each node, the predictor which improves the split-criterion the most is selected. The splitting criterion favors informative predictors over noisy ones, and can completely disregard irrelevant ones. Thus, in RFs an additional feature selection procedure is unnecessary. Based on the internal mechanism for selecting predictors, the predictor importance or strength can be estimated. The importance measure attributed to the splitting predictor is the accumulated improvement this predictor gives to the split-criterion at each split in each tree. RF also offers another method of estimating the predictor importance based on the OOB patterns [27]. When the tree is grown, the OOB patterns are passed down the tree, and the prediction error is recorded. Then the values for the given predictor are randomly permuted across the OOB samples, and the error is computed again. The importance measure is defined as the increase in error. This measure is computed for every tree, then averaged over the entire forest and divided by the standard deviation over the entire forest. Such a measure is presented in Figure 8. Note that information about the predictor importance is a key factor, which helps to improve the interpretability of the model and can be used for feature selection for other models.

Model interpretability is an emerging area in ML that aims to make the model more transparent and strengthen confidence in its results. This topic is also explored in electricity demand forecasting literature [44]. In [45], it was shown that the predictor importance is related to the model sensitivity to inputs and also to the method of importance estimation. An LSTM-based model, which is proposed in [45], is equipped with a built-in mechanism based on a mixture attention technique for temporal importance estimation of predictors. In the experimental study, this model demonstrated higher sensitivity to inputs than treebased models (RF and XGBoost) which showed very low sensitivity on the predictors except one, which strongly dominated (the authors used built-in functions of scikit-learn to calculate the predictor importance for tree-based models). In our study, the predictor importance is more diverse (see Figure 8), which may result from the fact that our trees are very deep and thus involve a great number of predictors. Note that tree-based models enhance interpretability not only through built-in mechanisms of predictor importance estimation, which show predictive power of individual predictors, but also through their flowchart-like tree structure. They can be interpreted simply by plotting a tree and observing how the splits are made and what is the arrangement of the leaves. It should be noted, however, that while following the path that a single tree takes to make a decision is trivial and self-explanatory, following the paths of hundreds of trees in the ensemble is much more difficult. To facilitate this, in [46], model compression methods were proposed that transform a tree ensemble into a single tree that approximates the same decision function.