Random vector functional link neural network based ensemble deep learning for short-term load forecasting

The EWT is an automatic signal decomposition algorithm with solid theoretical foundations and remarkable effectiveness in decomposing non-stationary time series data [24]. Unlike discrete wavelet transform (DWT) and EMD [25], EWT precisely investigates the time series in the Fourier domain after fast Fourier transform (FFT). It realizes the spectrum separation using band-pass filtering with the data-driven filter banks.

Figure 1 shows the EWT’s regular procedures. In the EWT, limited freedom is provided for selecting wavelets. The algorithm employs Littlewood-Paley and Meyer’s wavelets because of the analytic accessibility of the Fourier domain’s closed-form expression [26]. In [24], the formulations of these band-pass filters are denoted by Equations 1 and 2

with a transitional band width parameter $\gamma$ satisfying $\gamma\leq\min_{n}\frac{\omega_{n+1}-\omega_{n}}{\omega_{n+1}+\omega_{n}}$. The most common function $\zeta(x)$ in Equation 1 and 2 is presented in Equation 3. This empowers the formulated empirical scaling and wavelet function $\{\hat{\phi}_{1}(\omega),\{\hat{\psi}_{n}(\omega)\}_{n=1}^{N}\}$ to be a tight frame of $L^{2}(\mathbb{R})$ [27].

It can be observed that $\{\hat{\phi}_{1}(\omega),\{\hat{\psi}_{n}(\omega)\}_{n=1}^{N}\}$ are used as band-pass filters centered at assorted center frequencies.

Plentiful works utilize signal decomposition techniques as a feature engineering block for the forecasting algorithms [9, 10, 28, 7, 29, 8, 30], however, most do not implement the decomposition in a proper way [8, 30]. As mentioned in [7, 8, 30], direct application of signal decomposition algorithm to the whole time series causes the data leakage problem in terms of forecasting. The decomposed data are actually the output from convolution operations and the future data definitely are involved during the convolution. Therefore, decomposition of the whole time series is incorrect and improper, especially for establishing forecasting models.

Some solutions are proposed to avoid the future data leakage problem for decomposition-based forecasting models, such as the data-driven padding [7], moving window strategy [30] and walk-forward decomposition [8]. The data-driven padding approach is to train a simple learning algorithm which aims at padding its forecast to the end of the time series [7]. The moving window strategy only decomposes the data located in the window (order) and then the decomposed series are fed into forecasting models [30]. Different from the moving window strategy, only part of the decomposed sub-series are used as input in the walk-forward decomposition. The moving window strategy is a subset of the walk-forward decomposition. When the order is equal to the window, the moving window strategy and the walk-forward decomposition are the same.

This paper adopts the walk-forward decomposition for the EWT. The walk-forward EWT decomposes the data in a rolling window $w$, which consists of $x(t-1),x(t-2),......x(t-w)$, into $k$ scales with the aim to predict $x(t)$. Then only the last order data points are used as input for the forecasting model. Therefore, only historical observations are involved both in the decomposition process and the model’ training.

Inspired by the deep representation learning, the deep RVFL is an extension of the RVFL with a shallow structure [23]. The deep RVFL is established by stacking multiple enhancement layers to achieve deep representation learning. The clean data are fed into each enhancement layer to guide the random features’ generation. In this fashion, the enhancement features of hidden layers are generated based on the information from the clean data and the features from the previous layer. A diverse set of features is generated with the help of hierarchical structures. Ensemble learning is introduced into the deep RVFL architecture to formulate the ensemble deep RVFL (edRVFL). Different from the popular deep learning models with a single output layer, the edRVFL trains multiple output layers based on all the hidden features. Finally, the forecasts from all output layers are combined for forecasting.

For the sake of presentation simplicity, we only present the edRVFL with a structure of $L$ enhancement layers and there are $N$ enhancement nodes in each layer. Figure 2 shows the architecture of the edRVFL network.

Suppose that the input data is $\mathbf{X}\in\mathbb{R}^{n\times d}$, where $n$ and $d$ represent the number of samples and feature dimension, respectively. $d$ is the time lag (order) for the time series forecasting model. The features generated by the first enhancement layer are defined as

where $\mathbf{W_{1}}\in\mathbb{R}^{d\times N}$ represents the weight vector of the first enhancement layer, $\mathbf{H^{1}}\in\mathbb{R}^{n\times N}$ denotes the enhancement features and $g()$ is a non-linear activation function. The readers can refer to [31] for a comprehensive evaluation of different activation functions. Then, for the deeper enhancement layer $l$, the enhancement features can be computed as

where $\mathbf{W_{\textit{l}}}\in\mathbb{R}^{(d+N)\times N}$ and $\mathbf{H^{\textit{l}}}\in\mathbb{R}^{n\times N}$. The enhancement weight vectors $\mathbf{W_{1}}$ and $\mathbf{W_{\textit{l}}}$ are randomly initialized and remained fixed during training.

The edRVFL computes the output weights by splitting the task into $l$ small tasks. The output weights are calculated separately for each layer. There are several differences from using the last layer’s features and all layers’ features for decisions. Most deep learning models only use the last layer’s features for decisions, however, the information from the intermediate features is lost. Using all layers’ features requires a computation on the feature matrix with a huge dimension. Moreover, both of the above architectures only train one network, but our method benefits from the ensemble approach, which reduces the uncertainty of a single model.

The loss function of $l^{th}$ enhancement layer is defined as

where $\beta_{\textit{l}}$ denotes the output vector of $l^{th}$ layer and $\lambda$ is the regularization parameter. The minimization of $Loss_{\textit{l}}$ can be solved via a closed-form solution based on ridge regression [32].

where $\mathbf{D}=\mathbf{[H^{\textit{l}},X]}$. After computing all $\beta_{\textit{l}}$, the deep network can output $L$ forecasts. The final forecast is an ensemble of all outputs. Any forecast combination approach can be applied to this procedure [33]. According to the suggestions in [33], the mean or median operation is always likely to improve the forecast combination’s performance. Therefore, we use the mean and median as the combination operator. Correspondingly, two different edRVFLs are proposed, the Mea-edRVFL and Med-edRVFL.

The model EWT-edRVFL consists of two blocks, the walk-forward EWT decomposition and the edRVFL. The walk-forward EWT is first applied to the load data to extract some features in a causal fashion. Then the raw data concatenated with the sub-series are fed into the edRVFL with $L$ enhancement layers for learning purposes. The output weights $\beta_{\textit{l}}$ of the $l^{th}$ enhancement layer are computed according to Equation 7. Finally, we ensemble the $L$ forecasts with mean or median operation to obtain the output $\hat{y}$. Correspondingly, two different EWT-edRVFLs are proposed, the EWTMea-edRVFL and EWTMed-edRVFL.

Since the higher enhancement layer’s performance depends on the lower ones’, the hyper-parameters of the whole model are tuned in a layer-wise fashion. Once the shallow layer’s hyper-parameters are determined, then they are fixed and the cross-validation approach is applied to the next layer. Layer-wise cross-validation offers a different set of hyper-parameters for each layer. Therefore each enhancement layer has its own regularization parameter, which helps the overall edRVFL learns a diverse set of output layers.

Table I summarizes the descriptive statistics of the twenty load time series. These load data are collected from the states of South Australia (SA), Queensland (QLD), New South Wales (NSW), Victoria (VIC), and Tasmania (TAS) of the year 2020, which is significantly affected by Covid-19. Four months, January, April, July, and October are selected to reflect the four seasons’ characteristics as in [8, 34, 9]. The data are recorded every half an hour. Therefore, there are 48 data points per day.

A suitable and correct data pre-processing approach helps the machine learning model generate accurate outputs. We utilize the max-min normalization to pre-process the raw data. We assume that the maximum and minimum of the training set are $x_{max}$ and $x_{min}$, respectively. The data are transformed into the range [0,1] using the following equation:

where $x_{normalized}$ and $x$ represent the normalized and original time series, respectively.

All datasets are split into three sets, the training, validation and test set, to adopt the cross-validation [35]. The validation and test set account for 10% and 20% of the dataset, respectively. The remaining data are used as the training set.

Three forecasting error metrics are employed to appraise the accuracy of these models. The first error metric is the regular root mean square error (RMSE) whose definition is

where $L_{test}$ is the size of the test set, $x_{j}$ and $\hat{x_{j}}$ are the raw data and predictions. The second error metric implemented in the paper is the mean absolute scaled error (MASE) [36]. The definition of MASE is

where $L_{train}$ represents the size of training set. The denominator of MASE is the mean absolute error of the in-sample naive forecast. The third error metric is the Mean Absolute Percentage Error (MAPE) whose definition is

We compare the proposed model with many classical and state-of-the-art models. These models are Persistence model [2], ARIMA [3], SVR [5], MLP [13], LSTM [37], Temporal CNN (TCN) [38], hybrid EWT fuzzy cognitive map (FCM) learned with SVR (EWTFCMSVR) [7], Wavelet High-order FCM (WHFCM) [29], Laplacian ESN (LapESN) [39], EWTRVFL [8] and RVFL [6]. The previous one day, 48 data points are used as input for all the models as in [8]. To achieve a fair comparison, all models’ hyper-parameters are optimized by cross-validation. The hyper-parameter search space is presented in Table II. The decomposition level for the walk-forward EWT is set to 2 according to the conclusion and suggestions in [8]. Some parameters are not involved in the optimization process and they are set to the same values for all the relevant models, which include the batch size equals to 32, learning rate equals to 0.001 and epochs equal to 200.

Tables III, IV and V summarize the performance on the test sets. The numbers in bold represent the corresponding model’s performance is the best on the specific time series. Figures 3, 4, 5, 6 and 7 present the comparison of raw data and the forecasts generated by the proposed model. It is clear to find that the proposed model anticipates future trends, cycles, and fluctuations accurately. Statistical tests are implemented to investigate the difference among all the models further. We first implement the Friedman test, and the $p$-value is smaller than 0.05, which represents that these forecasting models are significantly different on these twenty datasets. Therefore, a post-hoc Nemenyi test is utilized to distinguish them [40]. The critical distance of the Nemenyi test is calculated by:

where $q_{\alpha}$ is the critical value coming from the studentized range statistic divided by $\sqrt{2}$, $k$ represents the number of models and $N_{d}$ is the number of datasets [40]. Figure 8 represents the Nemenyi test results. The figures show that the models that achieve excellent performance are at the top, whereas the model with the worst performance is at the bottom. Some consistent conclusions can be drawn from the Nemenyi test results of three error metrics. The Persistence method is the tailender because it learns nothing about the patterns. ARIMA is a penultimate because of its simple linear structure. The LSTM model outperforms many benchmark models except the EWTRVFL and the model proposed in this paper. Figure 8 demonstrates the superiority of the proposed models because they are always at the top. Another finding is that the edRVFL with mean ensemble operator is better than the median operator. A pair-wise Nemenyi post-hoc statistical comparison is further conducted and the $p$ values are shown in Tables VII, VIII and IX. The $p$ values smaller than 0.05 indicate that the two corresponding models are significantly different. The negative one in the diagonal positions represent that it is meaningless to compare the model with itself. The proposed EWTMea-edRVFL is significantly different from Persistence, ARIMA, SVR, MLP, LSTM, TCN, EWT-FCM-SVR, WHFCM, LapESN, and RVFL. The Mea-edRVFL and Med-edRVFL do not show significant superiority over LSTM, WHFCM, Lap ESN, RVFL, EWT-RVFL, and the EWT-based edRVFL models.

Table X records the simulation time for optimization and training time. It is worth noting that the optimization time is the time of the cross-validation using grid-search. The training time represents the time that the model is trained using the hyper-parameters selected by the cross-validation. The time for RVFL-related models is the summation of twenty runs. Several phenomenons are concluded according to Table X. The most time-consuming model is the LSTM because of its recurrent structure which processes the data in a sequential order. The hybrid RVFL model with EWT is more time-consuming than RVFL-related models. For example, the EWT-RVFL and EWT-edRVFL are more time-consuming than RVFL and edRVFL, respectively. Therefore, the main computation is in the walk-forward EWT decomposition block because it happens at each step.