A hybrid model of modal decomposition and gated recurrent units for short-term load forecasting

Variational mode decomposition

VMD decomposes the original series into several intrinsic mode functions (IMFs), each of which is a sub-sequence of the frequency modulation and amplitude modulation (Dragomiretskiy & Zosso, 2014; Zhang & Guo, 2020). VMD demodulates the IMF to its own fundamental frequency bandwidth and aims to minimize total modal bandwidth to find the optimal IMF. The VMD includes both variational construction and a variational solution.

To obtain the analytical signals of each IMF and the corresponding unilateral spectrum, The modal function obtained by decomposition, represented as: u_k(t), k=1 , 2, …, K is processed with Hilbert transform (Huang et al., 1998) as follows: (1)${\displaystyle \left(\delta \left(t\right)+\frac{\mathrm{j}}{\pi t}\right)\ast u_k\left(t\right)}$

where δ(t) is the unit impulse signal.

Center frequency, represented as ω_k, is then multiplied by the exponential term e^−jω_kt to modulate the spectrum of the mode to its fundamental frequency: (2)${\displaystyle \left[\left(\delta \left(t\right)+\frac{j}{\pi t}\right)\ast u_k\left(t\right)\right]e^{-j{\omega }_{k}t}}$

The bandwidth of the mode is then calculated, using the solution of the 2-Norm of the modulated signal gradient to solve the variational constraint problem: (3)${\displaystyle \begin{array}{c}\underset{u_k,{\omega }_k}{min}\left\{\sum _{k=1}^K{∥{\partial }_t\left[\left(\delta \left(t\right)+\frac{j}{\pi t}\right)u_k\left(t\right)\right]e^{-j{\omega }_{k}t}∥}_2^2\right\}\hfill \\ s.t.\sum _{k=1}^K{u}_k\left(t\right)=f\left(t\right)\hfill \end{array}}$

where f (t) is the original signal.

The variational constraint problem is then transformed into an unconstrained problem using the Lagrange multiplier method and quadratic multiplication operator alternation algorithm. By introducing the Lagrange multiplication operator θ(t) and penalty factor C into the constraint problem, the unconstrained problem is as follows: (4)${\displaystyle L\left(u_k,{\omega }_k,\theta \right)=C\sum _{k=1}^K{∥{\partial }_t\left[\left(\delta \left(t\right)+\frac{j}{\pi t}\right)u_k\left(t\right)\right]e^{-j{\omega }_{k}t}∥}_2^2+{∥f\left(t\right)-\sum _{k=1}^K{u}_k\left(t\right)∥}_2^2+〈\theta \left(t\right),f\left(t\right)-\sum _{k=1}^K{u}_k\left(t\right)〉.}$

Then, the optimization problem of u_k can be obtained by using the multiplication operator alternating algorithm: (5)${\displaystyle u_k^{n+1}=\underset{u_k\in X}{\mathrm{argmin}}\left\{\theta {∥{\partial }_t\left[\left(\delta \left(t\right)+\frac{j}{\pi t}\right)u_k\left(t\right)\right]e^{-j{\omega }_{k}t}∥}_2^2+{∥f\left(t\right)-\sum _{k=1}^K{u}_k\left(t\right)∥}_2^2\right.\left.+{∥f\left(t\right)-\sum _{i=1}^K{u}_i\left(t\right)+\frac{\theta \left(t\right)}{2}∥}_2^2\right\}}$

where i is the iteration control parameter. By using the Parseval/Plancherel Fourier equidistant method under 2-norm, it can be obtained as follows: (6)${\displaystyle {\hat{u}}_k^{n+1}=\underset{{\hat{u}}_k,u_k\in X}{argmin}\left\{{\int }_0^{\infty }4\theta {\left(\omega -{\omega }_k\right)}^2{\left|{\hat{u}}_k\left(\omega \right)\right|}^2+2{\left|\hat{f}\left(\omega \right)-\sum _{i=1}^K{\hat{u}}_i\left(\omega \right)+\frac{\hat{\theta }\left(\omega \right)}{2}\right|}^{2}d\omega \right\}.}$

Therefore, the optimized solution for this quadratic problem is: (7)${\displaystyle {\hat{u}}_k^{n+1}\left(\omega \right)=\frac{\hat{f}\left(\omega \right)-\sum _{i\ne k,i=1}^K{\hat{u}}_i\left(\omega \right)+\frac{\hat{\theta }\left(\omega \right)}{2}}{1+2C{\left(\omega -{\omega }_k\right)}^2}.}$

Finally, the center frequency can calculated using the following quadratic formula: (8)${\displaystyle {\omega }_k^{n+1}=\frac{{\int }_0^{\infty }\omega {\left|{\hat{u}}_k\left(\omega \right)\right|}^2\mathrm{d}\omega }{{\int }_0^{\infty }{\left|{\hat{u}}_k\left(\omega \right)\right|}^2\mathrm{d}\omega }.}$

Adaptive VMD

As a non-recursive decomposition method, the number of modes in VMD needs to be set in advance. When the number of modes is too small, there is insufficient decomposition; also, some modal functions can result in false spectrum and spectrum breakage (Jiang, Shen & Shi, 2018). In our previous work, the number of modes is set according to the cross-correlation coefficient (Shang, Li & Wu, 2023). In the proposed adaptive VMD model, to improve reliability the number of modes is jointly determined by the cross-correlation coefficient and the average sample entropy.

GRU network

The GRU network is a type of recurrent neural network (RNN; Cho et al., 2014) that has been proposed to solve the problems of long-term memory in back propagation. The GRU network performs in a similar way to the LSTM network but is computationally cheaper. The structure of a GRU network is shown in Fig. 1 where x_t is the input at the current node t, h_t-1 is the hidden state of transmission at the previous node t-1, y_t is the output, and h_t is the hidden state at t.

The GRU network has two gate states, as shown in Fig. 2. Here, r and z are the reset gate and the update gate, respectively, and can be written as: (16)${\displaystyle r=\sigma \left(W_r\left[h_{t-1},x_t\right]+b_r\right)}$ (17)${\displaystyle z=\sigma \left(W_z\left[h_{t-1},x_t\right]+b_z\right)}$

where σ is the sigmoid function which constrains the data in the interval [0,1]; W and b are the weight and threshold of the networks.

According to the reset gate r and the hidden state of h_t-1, the reset signal can be obtained, as follows: (18)${\displaystyle h_{t-1}^{{}^{\prime }}=h_{t-1}\otimes r}$

where ⊗ is the Hamiltonian operator.

The transmission state $h_t^{{}^{\prime }}$ is written as: (19)${\displaystyle h_t^{{}^{\prime }}=tanh\left(W_g\left[h_{t-1}^{{}^{\prime }},x_t\right]+b_g\right).}$

Lastly, the transmission state of h_t can be written as: (20)${\displaystyle h_t=\left(1-z\right)\otimes h_{t-1}+z\otimes h^{\prime }.}$

The improved optimization algorithm

The Adam algorithm updates weights by calculating the first and second moments of the gradient, improving the slow convergence problem caused by the fixed learning rate in the gradient descent method. The Adam algorithm was used on the random adjustment parameters in this study to effectively improve the convergence rate.

Firstly, initializing the learning rate µ, and using optimization parameters W_t, the gradient g_t, the first moment m_t of the gradient and the second moment v_t of the gradient can be calculated, iteratively, as follows: (21)${\displaystyle g_t=\nabla W_{t}f\left(W_t\right)}$ (22)${\displaystyle m_t={\beta }_1{m}_t+\left(1-{\beta }_1\right)g_t}$ (23)${\displaystyle v_t={\beta }_2{v}_{t-1}+\left(1-{\beta }_2\right)g_t^2}$

where β₁ and β₂ are the decay rates of the first and the second moments, respectively.

Then, the deviations of the first and second moments of the gradient can be calculated as: (24)${\displaystyle m_t^{{}^{\prime }}=m_t/\left(1-{\beta }_1\right)-{\eta }_1{g}_t}$ (25)${\displaystyle v_t^{{}^{\prime }}=v_t/\left(1-{\beta }_{\mathit{2}}\right)-{\eta }_2{g}_t^2}$

where η₁ and η₂ are the adjustment parameters of the first and second moments and the random numbers on interval [0,1], respectively.

Lastly, the parameters can be updated using the following formula: (26)${\displaystyle w_t=w_{t-1}-\mu m_t^{{}^{\prime }}/\left(\sqrt{v_t^{{}^{\prime }}}+\varepsilon \right)}$

where ɛ is the allowable error to prevent a zero value in the iterative process.

Forecasting process

In this work, a hybrid model of adaptive VMD and the GRU network is proposed to reduce frequency aliasing and eliminate the randomness of the load series. The electrical load series is represented as {x(1), x(2), …, x(N)}, where N is the number of samples. When C_i is the decomposed mode, the prototype mode is calculated, as follows: (27)${\displaystyle C_i=\left\{c_i\left(1\right),c_i\left(2\right),\dots ,c_i\left(N\right)\right\},i=1,2,\dots ,M.}$

The forecasting value of the prototype mode at time (N +1) is ${\hat{c}}_i\left(N+1\right)$ using the GRU networks. Reconstructing other components after removing the residual sequence, the forecasting result at time N +1 is calculated, as follows: (28)${\displaystyle \hat{x}\left(N+1\right)=\sum _{i=2}^M{\hat{c}}_i\left(N+1\right).}$

Figure 3 illustrates the load forecasting process of the hybrid model. First, the load series is decomposed using the adaptive VMD model. Then, each mode is forecasted in the next time using the improved GRU model. Finally, the sum of the IMF is the forecasting result of the load series. The forecasting process is as follows:

(a) The IMF of each modal component is obtained by initializing the number of modes and decomposing the load data using the adaptive VMD model.

(b) The cross-correlation coefficient ρ_c between the noise residue and IMF is calculated, using Eq. (11).

(c) The sample entropy of each mode and the modal residue, and the average sample entropy (ASE) are calculated using Eqs. (9) to (15) under different numbers of modes.

(d) The decomposition number K is selected, corresponding to the minimum of ASE and ρ_c.

(e) The modal components at the next time point are then forecasted using the improved GRU network.

(f) The final result is obtained by reconstructing the modal forecasting components.

Results of the adaptive VMD model

The ASE of each component, except the residual and cross-correlation coefficient ρ_c, were calculated, and then the minimum values of the mode number K were identified corresponding to the minimum ASE and ρ_c. The minimum ASE indicates that the similarity of each IMF was high, and that the sequence was more “orderly.” The minimum ρ_c means that the correlation between the residual sequence and the modal sequence was the smallest, meaning the reconstructed mode, except the residual, is closer to the real load sequence. Figure 4 shows the ASE and the cross-correlation coefficient ρ_c at different time points; ASE reached the minimum value when the number of decompositions was 5, and ρ_c was the smallest when K = 5.

The proposed adaptive VMD model decomposes the original electrical load series. Figure 5 shows the decomposition results of the load series of a day, consisting of 96 samples. In this figure, mode IMF4 has a small amplitude and violent random noise, giving it a residual sequence. The frequencies of the other modes decrease from top to bottom, revealing the short-term and long-term characteristics of the load series.

Results of the improved GRU network

A training set was created from 100 groups of load series, with the other load series used as the test set. The default values of related parameters are given in Table 1.

The forecasting performance of the improved GRU network was compared with the original GRU model, with both methods using the same training and test data. Table 2 shows the comparison results of both forecasting error and time. Each error index was based on the forecasting value of 672 datapoints. There were 100 groups of error index values in total, and the average value was the error result shown in Table 2. The forecasting error value of the improved GRU model proposed in this article was less than the forecasting error value of the original GRU model. The forecasting time of the improved GRU model was also significantly reduced compared to the original GRU model because the Adam algorithm uses a random adjustment of parameters.

Table 3 shows the forecasting error for different numbers of modes. When the number of modes was four, the forecasting error was the smallest, verifying the reliability of the adaptive VMD model.

Results of the hybrid model

The proposed forecasting model was then verified using different types of measured data. Figure 6A and 7A show the actual valuables and forecasting results (original GRU model and our hybrid model) on a working day and a non-working day, respectively. In Fig. 6A, the RMSE of the original GRU model was 335.7 MW and the RMSE of the proposed model was 334.5 MW. In Fig. 7A, the RMSE of the original GRU model and the proposed model were 335.9 MW and 334.6 MW, respectively. Since these differences were small, Figs. 6B and 7B further illustrate the comparisons with smaller units on the load (y-axis) and less forecasting points (x-axis). These figures show that the proposed hybrid model had better forecasting performance than the original GRU model.

The number of hidden neurons is an important parameter that can affect forecasting performance. Table 4 shows the forecasting results of different numbers of neurons of the hidden layer, with all other parameters being optimal. The forecasting error was the smallest when the number of neurons was 40, but the optimal number of neurons was different with different datasets.

The proposed hybrid forecasting model was also compared with the following classical statistical models: the ARIMA model, support vector regression (SVR), machine learning (Elman neural network), and the combined model. The parameter settings of the ARIMA model were set based on Lee & Ko (2021), and the selection of model order was based on the AIC criterion. The kernel function of the SVR model was the Gaussian radial basis function and the kernel parameters were optimized based on Sina & Kaur (2020). The parameter settings of the Elman neural network were set according to Xie et al. (2020), the optimization algorithm adopted the traditional gradient descent method, and the number of neurons in both hidden layers was 40. The single method selection and parameter settings of the combined model were based on Li & Chang (2018).

Table 5 illustrates the comparison of forecasting error of the above methods. Compared with the traditional statistical model and machine learning, the hybrid model proposed in this work had a higher forecasting accuracy. The traditional statistical learning method only obtained the evolution characteristics of the time series based on a limited number of samples, making it difficult to forecast long-term evolution and reversal characteristics of a time series. The machine learning methods, such as the radial basis function (RBF) and back-propagation (BP) neural network, have poor forecasting ability of a time series, and the Elman network is unable to forecast the long-term dependence of a time series. The proposed hybrid model, based on the GRU network, is a deep learning method, which obtains the evolution characteristics of sequences based on a large number of data, so the forecasting accuracy is higher. Because it is a deep learning method, the training time of the hybrid model based on the GRU network is much longer than that of the traditional statistical model and machine learning method.

Finally, the forecasting performance of the hybrid model in this work was compared with three state of the art methods: the LSTM network, the GRU network and the QRNN model. The LSTM algorithm and parameter settings were based on Rafi et al. (2021). The selection and parameter value of the GRU network were based on Shen et al. (2021). The setting and parameter values of the QRNN network were based on Cannon (2011). Table 6 compares the forecasting performance of the above methods. The average forecasting error value shows that, compared with the LSTM forecasting method combined with EMD and VMD, the AVMD model proposed in this article decomposed the sequence more accurately and improved the accuracy of the GRU network forecasting results. The maximum and minimum errors also verified the stability and reliability of the hybrid model. Compared with other deep learning networks, the parallel structure of the hybrid model based on the GRU network significantly shortened the training time, making this model suitable for short-term power load forecasting.

Because it is a deep learning method, the hybrid model requires a large number of training samples, and therefore the training time is longer than that of machine learning methods. Training time could be further reduced by reducing the number of training samples and batch size through a consideration of the periodicity of the load series. Furthermore, based on actual data, a reasonable network structure could be optimized, such as the number of network layers and nodes, without significantly reducing forecasting accuracy.