A Novel Hybrid Method for Multi-Step Short-Term 70 m Wind Speed Prediction Based on Modal Reconstruction and STL-VMD-BiLSTM

Abstract

In the context of achieving the goals of carbon peaking and carbon neutrality, the development of clean resources has become an essential strategic support for the low-carbon energy transition. This paper presents a method for the modal decomposition and reconstruction of time series to enhance the prediction accuracy and performance regarding the 70 m wind speed. The experimental results indicate that the STL-VMD-BiLSTM hybrid algorithm proposed in this paper outperforms the STL-BiLSTM and VMD-BiLSTM models in forecasting accuracy, particularly in extracting nonlinearity characteristics and effectively capturing wind speed extremes. Compared with other machine learning algorithms, including the STL-VMD-LGBM, STL-VMD-SVR and STL-VMD-RF models, the STL-VMD-BiLSTM model demonstrates superior performance. The average evaluation criteria, including the RMSE, MAE and R², for the proposed model, from t + 15 to t + 120 show improvements to 0.582–0.753 m/s, 0.437–0.573 m/s and 0.915–0.951, respectively.

1. Introduction

In the context of addressing global climate change as the greatest non-traditional security challenge, carbon reduction has become one of the most important principles and policies. Electricity generation accounts for more than 70% of greenhouse gas emissions and more than 90% of carbon gas emissions [1,2]. In the context of achieving the goals of carbon peaking and carbon neutrality, the comprehensive development and utilization of low-carbon clean energy to meet the increasing demand for carbon-emitting energy, and the gradual and comprehensive replacement of the existing stock, have become important strategic choices for the realization of the sustainable development of society [3].

The 2024 report of the Global Wind Energy Council (GWEC) indicates that the total installations of the global wind industry reached 117 GW in 2023, representing a 50% year-on-year increase from 2022. The future global demand for renewable energy will gradually increase, and the wind power installations in the 2024–2030 period will grow upwards by 10%, especially increasing in Europe, Asia and America [4]. Therefore, wind energy, as an important zero-carbon equivalent source of electricity, will play a key role in the energy transition. However, the wind speed is characterized by randomness, volatility and intermittency [5], which pose many challenges to the safe and stable operation and effective dispatch of power systems [6]. Therefore, the accurate prediction of the wind speed is crucial to improve the utilization of wind energy resources and enhance the ability of meteorological services to ensure the energy supply. In this regard, it is of great practical importance to improve the accuracy of wind speed prediction at different time scales, to ensure the safe and stable operation of the power grid and promote the development of the green economy by vigorously developing energy storage technology.

The forecasting methods for wind speed series are mainly based on physical or statistical methods. Physical methods predict the atmospheric wind field based on the equations of motion and the interaction laws among the physical parameters influencing the wind field. There are several global and regional models based on physical methods, including the European Center for Medium-Range Weather Forecasts (ECMWF) model and the China Global/Regional Assimilation and Prediction System (CGRAPS) model, which possess a temporal and spatial resolution of 3 h and 9 km, respectively. With the ongoing advancements in numerical weather prediction toward refinement [7], regional models such as GRAPES-MESO and WRF, which utilize background field data along with initial and boundary conditions provided by global models, have enhanced their temporal and spatial resolutions to 1 h and 3 km or even finer. The application of the four-dimensional assimilation forecast technique utilizing a numerically dense observation network can significantly improve both the spatial and temporal resolutions, as well as the forecast accuracy, of the wind speed [8,9]. However, in regions characterized by pronounced local climate features and limited monitoring resources, such as Northwest China, physical methods employing numerical weather prediction still exhibit a disparity between the forecast values and actual observations [10,11]. In addition, numerical models that rely on discrete computational systems possess high complexity and require time-consuming computation for short-range wind speed forecasts in complex terrain due to the physical parameterization schemes that depict the atmospheric dynamics and thermal processes as a complex nonlinear system. This complexity is further exacerbated by the high sensitivity of the atmospheric equations of motion to the boundary conditions and topographic parameters.

In contrast, statistical methods possess relatively lower time consumption and costs and are lightweight. Utilizing historical wind speed data, forecasting models can reconstruct the wind speed through traditional statistical methods or artificial intelligence methods. For smooth wind speed prediction, models such as the Autoregressive Moving Average (ARMA) are used, which can effectively fit and predict time series data that exhibit smoothness and linear trends [12]. Traditional statistical methods like the Autoregressive Integrated Moving Average (ARIMA) model can effectively capture the nonsmooth and nonlinear characteristics of the wind speed [13]. However, forecasting models that rely on traditional statistical methods tend to struggle with accuracy when dealing with wind speed series characterized by high complexity and volatility. This is primarily evident through the lag in the forecast results, unstable forecasting performance, and challenges in further improving the prediction accuracy [14]. Recently, artificial intelligence methods have emerged as innovative tools in the meteorological domain [15]. In contrast to physical and traditional statistical methods, artificial intelligence techniques (such as recurrent neural networks (RNNs) [16], convolutional neural networks (CNNs) [17] and the Self-Attention Mechanism Transformer [18]) demonstrate significant potential for the analysis of extensive datasets. In tpursuit of improved accuracy in wind power forecasting, Viet et al. (2023) propose a double-optimization approach aimed at developing a short-term wind power generation forecasting tool, which has shown heightened accuracy and can be widely implemented [19].

Combining high-temporal-resolution and high-spatial-resolution meteorological observation data with multi-source live fusion analysis technology enhances the performance of artificial intelligence methods for data mining through self-learning capabilities. This improves the ability to detect nonlinear wind speed variation features. Techniques such as integrated machine learning and deep learning methods based on RNN networks have demonstrated strong performance. Bentsen et al. (2023) focused on enhancing the forecasting precision of the multistep spatiotemporal wind speed by using a graph neural network (GNN) architecture to capture spatial dependencies, utilizing different update functions to learn temporal correlations [20]. Their study indicated that the GNN-FF Transformer could achieve superior outcomes for the 10 min and 1 h ahead forecasts. To mitigate the impact of noise on the characterization of nonlinear variations, the nonlinear single-time-series decomposition methods of Empirical Mode Decomposition (EMD) [21] and Extended Empirical Mode Decomposition (EEMD) [22] were implemented to establish a forecasting model. In [23], the EMD-CC Transformer connected an encoder and decoder through an attention mechanism. In this method, the EMD algorithm is used to decompose the time series to reveal the temporal changes across different time scales. Although the studies mentioned above concentrate on the decomposition and reconstruction of the wind speed series, there are several limitations associated with using a single-time-series decomposition method for wind speed forecasting: (1) the wavelet transform is prone to producing spurious harmonics that lack original physical significance when breaking down the original time series into multiple superpositions of wavelet functions [24]; (2) while EMD and EEMD can adaptively reflect the localized eigenfunctions of the original time series based on its nonlinear and smooth characteristics, they are also susceptible to issues such as modal aliasing phenomenon and endpoint effects. Thus, it is essential to develop prediction models that integrate various time series decomposition algorithms to enhance the short-term forecasting accuracy. To address the limitations of the aforementioned decomposition techniques, Dragomiretskiy et al. (2014) proposed the variational modal decomposition (VMD) method [25]. VMD can realize the effective separation of the intrinsic modal components (IMFs) and frequency domain segmentation of the signal, obtain the practical decomposition components of a given signal, and finally obtain the optimal solution of the variational problem [26].

From the above discussion, the main purpose of this work is to achieve credible wind speed forecasts at a 70 m height of the fan wheel hub, and a range of time series forecasting techniques in conjunction with modal reconstruction are constructed and assessed. This study introduces the numerical forecast as the background field factor for the nonlinear time series of the 70 m wind speed with poor regularity and high complexity, observed by a wind tower at a wind farm in Gansu Province. We adopt the forecast idea of “decomposition–prediction–reconstruction” to construct the 0–2 h short-term forecast model and carry out corresponding validation experiments. In summary, the main research content is as follows.

(1)

The STL method decomposes the original wind speed series and obtains the trend component, seasonal component, and residual component.

(2)

The VMD decomposition of the wind speed is utilized to obtain the high-frequency and low-frequency intrinsic modal components; the BiLSTM model of short-term wind speed prediction is constructed based on the results of the STL-VMD dual-time-series decomposition.

(3)

The feasibility and superiority of the STL-VMD dual-time-series decomposition method are verified through prediction experiments using the proposed method and the STL-BiLSTM, VMD-BiLSTM and BiLSTM methods with different datasets and prediction step sizes.

(4)

Through comparison tests between the proposed method, support vector regression (SVR), light gradient boosting machine (LGBM) and the random forest algorithm (RF), the prediction performance and superiority of the proposed model are validated across four datasets with different time step sizes.

2. Materials and Methods

2.1. Materials

The following two datasets are selected as research materials to carry out the analysis: (1) the continuous observation of the wind speed at a 70 m height acquired from the anemometer tower of a wind farm in Gansu Province during 2022, with a time resolution of 15 min, where the wind speed time series contained a total of 34,721 sample points, and (2) the ECWMF forecast products, 10-meter U (V) wind component and 100-meter U (V) wind component obtained at 12:00 UTC each day during 2022, denoted as 10 m U (V) and 100 m U (V), with a spatial resolution of 0.125° × 0.125°. Equation (1) presents the formula for the absolute wind speed at a 10 m and 100 m height, labeled as WS10m and WS100m, respectively. The WS10m and WS100m are used as input features instead of the 10-meter U (V) wind component and 100-meter U (V) wind component. The sample data are mapped to the interval [1] using the Min–Max normalization method, as shown in Equation (2).

$$WS=\sqrt{U^2+V^2}$$

(1)

$$x^{\prime }={\displaystyle \frac{x-x_{\min }}{x_{\max }-x_{\min }}}$$

(2)

where the normalized data is denoted by $x^{\prime }$; the maximum value and minimum value of the sample data x are denoted by x_max and x_min, respectively.

2.2. Methods

A novel hybrid model (STL-VMD-BiLSTM) is proposed in this study, with the flowchart shown in Figure 1. The technology roadmap for the STL-VMD-BiLSTM model includes four modules in the overall architecture, which are as follows.

(1)

In Section 1 (top left), the original wind speed series observed at the height of 70 m is decomposed by STL, and the trend component, seasonal and residual component are obtained.

(2)

In Section 2 (middle left), the VMD method is used to decompose the actual wind speed again, and the high-frequency and low-frequency intrinsic modal components of the actual wind speed data are obtained.

(3)

In Section 3 (bottom left), the ECMWF grid data are interpolated to the site via a linear interpolation scheme; the temporal resolution of the site data is 15 min.

(4)

In Section 4 (right), the training set and validation set of the time series model are constructed. In order to verify the feasibility and superiority of the proposed method, experiments and a contrastive analysis are carried out to train and optimize the models mentioned in this paper.

2.2.1. STL

Seasonal and trend decomposition (STL) is a time series decomposition algorithm based on the locally weighted regression smoothing method (LOESS) for the estimation of nonlinear relationships [27]. Assuming additive decomposition written as in Equation (3), the STL method consists of an inner loop and an outer loop. Using the LOESS algorithm for smoothing and low-pass filtering, the original wind speed time series data $Y_{\upsilon }$ are decomposed into the trend component $T_{\upsilon }$, the seasonal component $S_{\upsilon }$ and the residual component $R_{\upsilon }$ at the period $\upsilon$.

$$Y_{\upsilon }=T_{\upsilon }+S_{\upsilon }+R_{\upsilon }$$

(3)

The STL inner loop consists of a seasonal smoothing process and a trend smoothing process, which are employed to update the seasonal component and the trend component, respectively. There are 6 steps in the inner loop.

Step 1: Detrending involves subtracting the trending component of the results of the previous round $Y_{\upsilon }-{T_{\upsilon }}^{\left(k\right)}$ and the initial value ${T_{\upsilon }}^{\left(0\right)}$ is 0 in the first round.

Step 2: For cycle sub-series smoothing, LOESS is used to conduct regression smoothing for each cycle sub-series and extend them one cycle forward and one cycle backward, denoted as $C_{\upsilon }^{(k+1)}$, $\upsilon =-n(p)+1$ to $N+n(p)$. The LOESS parameters are chosen as $q=n_{(s)}$ and $d=1$.

Step 3: The low-pass filtering of the smoothed cycle sub-series is performed, which consists of running averages and a LOESS smoothing process with length parameters of $n_{(i)}$ and 3, respectively. The result is written as $L_{\upsilon }^{(k+1)}$, where $\upsilon$ is between 1 and $N.$ This step is equivalent to extracting the low pass of the cycle sub-series.

Step 4: The detrending of the smoothed cycle sub-series is performed to obtain the final seasonal component, described as $S_{\upsilon }^{(k+1)}=C_{\upsilon }^{(k+1)}-L_{\upsilon }^{(k+1)}$.

Step 5: Deseasonalization represents subtracting the seasonal component, described as $Y_{\upsilon }-S_{\upsilon }^{(k+1)}$.

Step 6: Trend smoothing involves applying the LOESS method with the parameters $q=n_{(t)}$ and $d=1$ to the series to obtain the trend component, labeled as $T_{\upsilon }^{(k+1)}$.

The outer loop can regulate the robustness weight by calculating the robustness weight for noise reduction based on the results of the inner loop, thus reducing the impact of noise on the next inner loop. The robustness weight (see Equation (4)) prevents the data series from having outliers, which results in a significant residual term and a large bias on the regression. Its bisquare function B and h are defined as in the following Equations (5) and (6):

$${\rho }_{\upsilon }=B(|R_{\upsilon }|/h)$$

(4)

$$h=6\ast median(|R_{\upsilon }|)$$

(5)

$$B(u)=\{\begin{array}{ll}{(1-u^2)}^2& 0\le u<1\\ 0& u>1\end{array}$$

(6)

where $n_{(i)}$ inner loops and $n_{(\mathrm{o})}$ outer loops are performed in the STL decomposition process. During Steps 2 and 6 in the inner loop of each iteration of LOESS regression, the neighborhood weight needs to be multiplied by the robustness weight to minimize the effect of outliers on LOESS.

The process of STL decomposition for the original wind speed time series is completed when the estimates of the trend and seasonal terms converge after multiple iterations of looping inside and outside.

2.2.2. VMD

The VMD process for wind speed time series consists of the construction and solution of constrained variational problems. Given that each intrinsic modal function (IMF) has a finite bandwidth with distinct center frequencies, while ensuring that the sum of all intrinsic modal components equals the original signal, the critical step in the VMD technique process is to determine the number k of modes in such a way that minimizes the sum of the estimated bandwidths for each IMF (Equation (7)) [28]. The wind speed time series is decomposed into k different time-varying and nonlinear frequency components.

$$\begin{array}{l}\underset{\{u_k,w_k\}}{\min }\{{\displaystyle \sum _k\Vert {{\partial }_t\left[\left(\sigma (t)+\frac{j}{\pi t}\right)\ast u_k(t)\right]e^{-j{w}_{k}t}\Vert }_2^2}\}\\ \begin{array}{ll}s.t.& {\displaystyle \sum _k{u}_k=f(t)}\end{array}\end{array}$$

(7)

where $u_k=\{u_1,\dots ,u_k\}$ and $w_k=\{w_1,\dots ,w_k\}$ are the k mode components and the corresponding center frequencies obtained from the VMD of the wind speed time series. $\sigma (t)$ is the Dirichlet function, and the coefficient $j$ is equal to $\sqrt{-1}$.

In order to determine the optimal solution of the constrained variational model, a penalty factor $\alpha$ and a Lagrange multiplier $\lambda (t)$ are introduced to transform this constrained variational problem into an unconstrained variational problem, as represented by Equation (8). The multiplicative operator alternating direction method is used to solve the variational problem. Decomposition is employed to extract the respective intrinsic modal components and center frequencies.

$$\begin{array}{ll}\hfill L(\{u_k\},\{w_k\},\lambda )=& \alpha {\displaystyle \sum _k\Vert {{\partial }_t\left[\left(\sigma (t)+\frac{j}{\pi t}\right)\ast u_k(t)\right]e^{-j{w}_{k}t}\Vert }_2^2}+\Vert {f(t)-{\displaystyle \sum _k{u}_k(t)}\Vert }_2^2\hfill \\ & +\langle \lambda (t),f(t)-{\displaystyle \sum _k{u}_k(t)}\rangle \hfill \end{array}$$

(8)

In this work, in order to find the optimal modal number k for VMD, the decomposition accuracy evaluation index MAPE (which is the average value of the ratio of the residual to the original value) is adopted, represented by Equation (9).

$$MAPE={\displaystyle \frac{1}{N}}{\displaystyle \sum _{t=1}^N\left|\frac{f-{\displaystyle \sum _K{u}_k(t)}}{f}\right|}$$

(9)

where $f$ is the raw wind speed time series, N is the length of the wind speed series and $f-{\displaystyle \sum _K{u}_k(t)}$ denotes the residual term.

2.2.3. BiLSTM

The bidirectional long short-term memory (BiLSTM) neural network, proposed by Graves et al. (2005), comprises forward and inverse LSTM layers [29]. Unlike the unidirectional LSTM neural network, BiLSTM can explore the intrinsic connections of wind time series data to the maximum extent, rather than simply reversing the input sequence [30]. In order to enhance the predictive abilities of traditional neural networks and overcome problems such as gradient disappearance with the growth of the learning period, BiLSTM realizes the storage and control of historical information by introducing the concepts of memory units and gating mechanisms to regulate the flow of information in the network layers, as with LSTM [31].

In the forward and inverse LSTM layers, the BiLSTM unit structure consists of 3 types of gating (denoted as in Equations (10)–(12)). The input gate controls the information added to the current unit, using the hyper tangent function to generate a new candidate quantity, the memory unit ${\tilde{c}}_t$, and thus control the added information (see Equation (13)). The forgetting gate controls the degree to which the previous unit has been forgotten, and its inputs are the outputs of the previous unit as well as the inputs of the current unit. Finally, the output gate controls the degree to which the current unit filters the information.

$$i_t=\sigma (W_i·\left[h_{t-1},x_t\right]+b_i)$$

(10)

$$f_t=\sigma (W_f·\left[h_{t-1},x_t\right]+b_f)$$

(11)

$$o_t=\sigma (W_o·\left[h_{t-1},x_t\right]+b_o)$$

(12)

$${\tilde{c}}_t=\tanh (W_c·\left[h_{t-1},x_t\right]+b_c)$$

(13)

$$c_t=f_t\odot c_{t-1}+i_i\odot {\tilde{c}}_t$$

(14)

$$\sigma (x)={\displaystyle \frac{1}{1+e^{-x}}}$$

(15)

For the input gate $i_t$, forgetting gate $f_t$, output gate $o_t$ and memory unit ${\tilde{c}}_t$, the corresponding weight matrixes are $W_i$, $W_f$, $W_c$, $W_o$ and the bias vectors are $b_i$, $b_f$, $b_c$, $b_o$; $c_t$ is the metacellular state and $\sigma$ is the activation function with value domain [0, 1], and these are described as in Equations (14) and (15). The total output of the BiLSTM computational unit at moment t is calculated as follows in Equations (16)–(18).

$$\overrightarrow{h_t}=\sigma (W_{F1}{x}_t+W_{F2}{\overrightarrow{h}}_{t-1}+\overrightarrow{b})$$

(16)

$${\overleftarrow{h}}_t=\sigma (W_{B1}{x}_t+W_{B2}{\overleftarrow{h}}_{t-1}+\overleftarrow{b})$$

(17)

$$y_t=W_{O1}{\overrightarrow{h}}_t+W_{O2}{\overleftarrow{h}}_{t-1}+b_y$$

(18)

where $x_t$, $\overrightarrow{h_t}$, ${\overleftarrow{h}}_t$, $y_t$ are the input vectors and output vectors of the implicit layer of forward and inverse propagation at time step $t$, respectively; $W_{F1}$, $W_{F2}$, $W_{B1}$, $W_{B2}$, $W_{O1}$, $W_{O2}$ are the matrix coefficients; $\overrightarrow{b}$, $\overleftarrow{b}$, $b_y$ are the bias vectors.

2.2.4. Performance Evaluation Criteria

In this study, three evaluation indexes are selected to comprehensively evaluate the prediction performance of the models mentioned. The root mean square error (RMSE) and the mean absolute error (MAE) are used to quantify the number of errors. The goodness of fit R-squared (R²) is used to measure the performance of the proposed wind speed prediction model and other involved prediction models. The R² score reflects the extent to which the variance of the prediction variable explains the variance of the observed variable. The closer the RMSE and MAE are to 0, the closer the R² is to 1 and the better the performance of the prediction model. The three evaluation criteria are calculated as follows:

$$MAE={\displaystyle \frac{1}{N}}{\displaystyle \sum _{i=1}^n\left|y_i-{\widehat{y}}_i\right|}$$

(19)

$$RMSE=\sqrt{{\displaystyle \frac{1}{N}}{{\displaystyle \sum _{i=1}^n(y_i-{\widehat{y}}_i)}}^2}$$

(20)

$$\begin{array}{l}{R}^2=1-{\displaystyle \frac{SSE}{SST}}\\ SSE={{\displaystyle \sum _{i=1}^N(y_i-{\widehat{y}}_i)}}^2\\ SST={{\displaystyle \sum _{i=1}^N(y_i-\overline{y})}}^2\end{array}$$

(21)

where $y_i$, ${\widehat{y}}_i$, $\overline{y}$ are the actual value, the prediction value and the expectancy value of the wind speed, respectively. N represents the length of the wind speed sequence. The sum of squares due to error (SSE) is the sum of squares of the residual difference between the observed and predicted values, representing the unexplained variability in the variable. The total sum of squares (SST) represents the sum of squares of the mean deviation between the observed and predicted values, signifying the total variation degree of the variable.

3. Results

3.1. Experimental Data Description

Given the seasonal variability, the 70 m wind speed dataset is segmented into four seasonal subsets. Figure 2 illustrates the Weibull distributions of these wind speed subsets, with K serving as the shape parameter, indicating how closely the data align with the mean (or peak). A larger (or smaller) K signifies a more concentrated (or dispersed) distribution. The scale parameter dictates the absolute magnitude of the distribution, reflecting the data’s range or scale. The location parameter (Loc) is utilized to adjust the standard distribution. The probabilty density function of Weibull minimum continuous random variable are denoted by the red lines. In spring, the maximum wind speed recorded is 24.80 m/s, while the mean wind speed is 7.116 m/s. As listed in

Table 1. Overview of some statistical values of wind speed data in 2022.

Case	Maximum (m/s)	Std (m/s)	Mean (m/s)	K	Loc	Scale
Spring	24.80	4.383	7.116	1.731	−0.012	8.033
Summer	20.99	3.579	6.921	2.055	−0.008	7.844
Fall	19.67	3.299	5.688	1.830	−0.028	6.452
Winter	20.280	2.84	4.529	1.732	−0.140	5.245

, the shape parameter K is the lowest (1.731) in spring, with the highest scale parameter (8.033), suggesting that the Weibull distribution for spring is the most dispersed, as illustrated in Figure 2. In contrast, the mean wind speed during winter stands at 2.84 m/s, which is the lowest value among all seasons. The results indicate that the wind speed reaches the largest in spring and has the most vital characteristic of gusty fluctuations, while it is the smallest in winter but still has the most robust volatility characteristics.

3.2. Parameter Selection

To validate the proposed model, the dataset is segmented into groups representing three consecutive months. The first subset, spanning from January to March 2022, is designated as dataset 1. Similarly, the final subset, covering October to December 2022, is labeled dataset 4. Each case is further divided into the training set (85%) and the test set (15%).

Table 2. Optimal input parameters of prediction models based on BiLSTM.

Model Parameter	Test Parameters	Optimal Parameter
Window_size L	(8,12,16,20,24,28)	20
Batch_size	(64,128,256,512)	256
Epoch	(10,15,20)	15
Neurons of fully connected layer	8
Rate	0.2
Activation function	ReLu
Objective function	MAE
Optimizer	Adam

presents the optimal parameters of the prediction models based on BiLSTM. Each model consists of three LSTM layers, two dropout layers with a dropout rate of 0.2, one fully connected layer, and an output layer. The optimal values of window_size (20), batch_size (256) and epoch (15) are determined through the grid parameter optimization algorithm. The input features of the BiLSTM model include wind speed forecast data at 10 m and 100 m (WS10m and WS100m), along with observed wind speed data at 70 m (WS70m). For instance, in the fourth case, the input and output data shapes are (7379, 20, 3) and (7379, 8), respectively, intended for the BiLSTM, SVR, LGBM and RF models. The optimal input parameters for the prediction models utilizing the machine learning techniques are detailed in

Table 3. Optimal input parameters of prediction models based on machine learning methods.

Model	Parameter
SVR	kernel	rbf
	epsilon	0.2
	shrink	True
	tol	0.001
LGBM	boosting_type	gbdt
	num_leaves	31
	learning_rate	0.02
	feature_fraction	0.9
RF	bagging_fraction	0.8
	bagging_freq	5
	estimators	10
	max_depth	5
	random_state	0

. The experiments were conducted on a 13th-Gen Intel (R) Core (TM) i7-1360P.

Table 4. Optimal input features of prediction models.

Type	Model	Optimal Features ($1<\mathit{k}\le \mathit{L},\mathit{L}=20$)	Label ($1\le \mathit{k}\le \mathit{N},\mathit{N}=8$)
1	RF SVR LGBM BiLSTM	WS10m(t − k), …, WS10m(t) WS100m(t − k), …, WS100m(t) WS70m(t − k), …, WS70m(t)	WS70m(t + k), …, WS70m(t + N)
2	STL-RF STL-SVR STL-LGBM STL-BiLSTM	WS10m(t − k), …, WS10m(t) WS100m(t − k), …, WS100m(t) WS70m(t − k), …, WS70m(t) STL-Trend(t − k), …, STL-Trend(t)	STL-Trend(t + k), …, STL-Trend(t + N)
		WS10m(t − k), …, WS10m(t) WS100m(t − k), …, WS100m(t) WS70m(t − k), …, WS70m(t) STL-Seasonal(t − k), …, STL-Seasonal(t)	STL-Seasonal(t + k), …, STL-Seasonal(t + N)
3	VMD-RF VMD-SVR VMD-LGBM VMD-BiLSTM	WS10m(t − k), …, WS10m(t) WS100m(t − k), …, WS100m(t) WS70m(t − k), …, WS70m(t) VMD-IMF(t − k), …, VMD-IMF(t)	VMD-IMF(t + k), …, VMD-IMF(t)
4	STL-VMD-RF STL-VMD-SVR STL-VMD-LGBM STL-VMD-BiLSTM	WS10m(t − k), …, WS10m(t) WS100m(t − k), …, WS100m(t) WS70m(t − k), …, WS70m(t) STL-Trend(t − k), …, STL-Trend(t) STL-Seasonal(t − k), …, STL-Seasonal(t) VMD-IMF(t − k), …, VMD-IMF(t)	VMD-IMF(t + k), …, VMD-IMF(t + N)

lists the optimal features for the prediction models. Specifically, for the first type of model, such as BiLSTM, the input features include the WS10m, WS10m(t) and WS70m(t − k) vectors spanning from t − k to t, while the output vector is the WS70m from t + 15 to t + 120. The second (and third) types of models are fed by the WS10m, WS10m(t) and WS70m(t − k) vectors from t − k to t, along with the decomposition components generated by STL and VMD. The outputs of these two model types consist of STL or VMD subcomponents from t + 15 to t + 120, which are then used to reconstruct the wind speed sequences. In contrast, the label of the fourth type of model is based on the VMD component vectors, such as IMF1.

3.3. Analysis of Proposed Models

3.3.1. Decomposed Results for STL

The wind speed time series have the characteristics of trends, seasonality, periodicity and random noise. The STL algorithm is used to decompose the original wind speed series (top) into the smoothed trend component (middle top), the seasonal component (middle lower) and the residual component (bottom), which are shown in Figure 3. The trend component curve changes more gradually and smoothly.

Table 5. Overview of statistical values of wind speed components decomposed by STL.

Case	STL-Trend (m/s)			STL-Seasonal (m/s)		STL-Remainder (m/s)		FT	FS
	Max	Min	Mean	Max	Min	Max	Min
Spring	23.674	0.570	7.116	1.790	−1.841	7.004	−5.474	0.971	0.292
Summer	20.164	−1.841	6.921	2.264	−2.199	6.281	−3.709	0.956	0.308
Fall	18.387	0.912	5.688	1.943	−1.952	3.554	−3.566	0.966	0.324
Winter	19.198	−0.022	4.529	1.251	−1.429	3.276	−2.830	0.968	0.337

lists some statistical values of the wind speed components decomposed by the STL algorithm; the strength of the trend (F_T) and seasonality (F_S) is defined as in Equations (22) and (23), respectively. The trend component, with a gentle feature, explains the low-frequency variation of the wind speed series and represents the long-term variation trend. For the strongly trended data, in four cases, the maximum ranges from 18.387 m/s to 23.674 m/s, and the strength of the trend FT is between 0.956 and 0.971. The seasonal component ranges from −2.199 m/s to 2.264 m/s and denotes the periodic variation of the time series. The strength of the seasonal component FS is between 0.292 and 0.337. The residual component represents changes not attributed to the periodic component and the trend component and is more random and approximately normally distributed with 0. The residual component is regarded as unpredictable random noise.

$$F_T=\max (0,1-{\displaystyle \frac{Var(R_t)}{Var(T_t+R_t)}})$$

(22)

$$F_S=\max (0,1-{\displaystyle \frac{Var(R_t)}{Var(S_t+R_t)}})$$

(23)

3.3.2. Decomposed Results for VMD

The trend components of the wind speed series are categorized into 10 groups using a three-month running period. The VMD algorithm is applied to decompose the observed wind speed data into several intrinsic modal components (IMFs) with the optimal decomposition modal number, denoted as K. In this study, the modal number is deemed optimal when two conditions are met: (1) the evaluation accuracy measured by the MAPE is less than 1% and (2) the MAPEs are calculated for various modal numbers K, ranging from 2 to 20. As K increases, the center frequency of each decomposed component tends to decrease. Furthermore, the change in the MAPEs when transitioning from K to K + 1 is minimal, indicating that K is considered to be optimal choice. Therefore, in the specific VMD process, we began by establishing the penalty coefficients, $\alpha =2000$, $\tau =0$ and $N\le 20$, followed by determining the modal number K of the VMD test. Finally, we calculated and compared the changes in the MAPE index for different values of the modal number K (see

Table 6. The variation characteristics of the MAPE with different modal numbers K.

Groups	Modal Number K
	2	3	4	5	6	7	8	10	20
1–3	4.546	2.490	1.614	1.136	0.784	0.688	0.461	0.292	0.144
2–4	4.844	2.594	1.793	1.177	0.867	0.67	0.457	0.300	0.079
3–5	4.575	2.409	1.716	1.203	0.757	0.612	0.442	0.292	0.134
4–6	4.260	2.285	1.593	1.145	0.914	0.574	0.420	0.276	0.137
5–7	4.162	2.305	1.552	1.107	0.722	0.606	0.412	0.282	0.095
6–7	4.156	2.317	1.489	1.164	0.699	0.588	0.407	0.283	0.074
7–9	4.339	2.369	1.683	1.262	0.744	0.625	0.411	0.300	0.145
8–10	4.195	2.079	1.659	1.148	0.717	0.598	0.392	0.269	0.136
9–11	4.534	2.395	1.429	1.170	0.889	0.572	0.413	0.280	0.076
10–12	4.658	2.672	1.846	1.241	0.880	0.796	0.421	0.278	0.154

). The MAPE index exhibits a rapid decline as the modal number K increases. When the modal number K is 6, the MAPE values for each dataset range from 0.600 to 0.914, remaining below 1%. As the modal number (K) continues to rise from 7 to 20, the reduction in the MAPE index becomes less pronounced, suggesting that the optimal modal number for VMD in wind speed analysis is 6. For groups 4–6, a modal number of (K = 7) can be used without significant time cost considerations.

The observed wind speed time series alongside the VMD results are shown in Figure 4. IMF1 represents the low-frequency variability in the wind speed data, reflecting their long-term trend. IMF2–4 are the detail components, which indicate the changes in the wind speed across various time scales. IMF5–6 are the stochastic components, which characterize the wind speed time series with randomness features.

3.3.3. Prediction Results with Different Time Series Decomposition Schemes

We next assess the feasibility and reliability of the time series dual hybrid decomposition algorithms used in this paper to enhance the performance of the wind speed short-range proximity model. The comparative analyses of the prediction performance of the models, including the proposed STL-VMD-BiLSTM model and the VMD-BiLSTM, STL-BiLSTM and BiLSTM models, are carried out across the four datasets from t + 15 to t + 120.

Table 7. Prediction results of models at specified steps based on BiLSTM in four cases.

Case	Model	t + 15			t + 60			t + 120
		RMSE	MAE	R2	RMSE	MAE	R2	RMSE	MAE	R2
(a) Spring	STL-VMD-BiLSTM	0.668	0.483	0.958	0.754	0.550	0.946	0.920	0.696	0.918
	VMD-BiLSTM	0.735	0.544	0.954	0.830	0.623	0.943	1.048	0.800	0.906
	STL-BiLSTM	0.915	0.668	0.919	1.184	0.861	0.868	2.108	1.516	0.548
	BiLSTM	1.320	0.982	0.865	2.073	1.521	0.643	2.725	2.008	0.290
(b) Summer	STL-VMD-BiLSTM	0.665	0.504	0.946	0.749	0.567	0.929	0.835	0.631	0.910
	VMD-BiLSTM	0.742	0.543	0.931	0.837	0.623	0.909	0.991	0.751	0.860
	STL-BiLSTM	1.058	0.785	0.845	1.136	0.850	0.814	1.899	0.423	0.392
	BiLSTM	1.401	1.070	0.758	2.017	1.519	0.399	2.271	1.740	0.118
(c) Fall	STL-VMD-BiLSTM	0.561	0.439	0.961	0.627	0.481	0.951	0.715	0.558	0.939
	VMD-BiLSTM	0.568	0.428	0.964	0.656	0.499	0.952	0.803	0.661	0.914
	STL-BiLSTM	0.783	0.587	0.932	0.881	0.679	0.909	1.362	1.166	0.669
	BiLSTM	1.026	0.786	0.878	1.542	1.138	0.667	1.919	1.555	0.236
(d) Winter	STL-VMD-BiLSTM	0.435	0.321	0.937	0.470	0.347	0.922	0.542	0.405	0.893
	VMD-BiLSTM	0.439	0.323	0.934	0.508	0.380	0.909	0.610	0.460	0.859
	STL-BiLSTM	0.621	0.444	0.859	0.696	0.522	0.809	1.189	0.941	0.298
	BiLSTM	0.836	0.589	0.750	1.223	0.953	0.317	1.458	1.180	−0.339

lists prediction error scores of the BiLSTM models integrating time series decomposition methods at specified steps for these four sets. The performance of the short-range prediction model using different datasets is evaluated in terms of three main aspects: the goodness-of-fit R², which evaluates how well the models fit the data; the RMSE, which measures the extent of deviation between the observed and predicted values; and the MAE, which reflects the absolute error magnitude between the predicted and observed values, as well as the degree of deviation between the predicted values.

Figure 5, Figure 6 and Figure 7 illustrate the forecasted and observed values (represented as curved lines), with the scatter plots demonstrating the four models developed by integrating different time series decomposition methods. In Figure 5, the observed results and the predictions from the four methods at time t + 15 are closely aligned. The proposed STL-VMD-BiLSTM model exhibits the highest R² values, ranging from 0.937 to 0.958. In contrast, the BiLSTM model shows the lowest R² values, which range from 0.750 to 0.878.

As illustrated in Figure 6, none of the four models effectively and thoroughly replicates the gusty high winds or calm winds at t + 60, a discrepancy that is particularly evident in spring, when gusty high winds are prevalent, and in winter, characterized by calm winds. The R² for the STL-VMD-BiLSTM model, noted for its robust descriptive capabilities, ranges from 0.922 to 0.951. The R² scores for the VMD-BiLSTM and STL-BiLSTM models range from 0.909 to 0.952 and from 0.809 to 0.909 at t + 60, respectively. The BiLSTM model exhibits an R² value of less than 0.667.

Figure 7 presents a comparison of the forecast and observed curves at the t + 120 timeframe, indicating a gradual decline in the predictive abilities of all models with the increase in the forecast length. Following the application of the decomposition and reconstruction methods, the forecasting capability of the STL-VMD-BiLSTM model shows a significant enhancement, with an R² range of 0.893 to 0.939. The VMD-BiLSTM model follows with an R² range of 0.859 to 0.914, while the STL-BiLSTM model slightly outperforms the BiLSTM model. The challenges in capturing the characteristics of wind speed changes and accurately forecasting wind speed extremes are principally reflected in the overestimation of extreme minima and the hysteresis evident in the forecast curves, which represent typical limitations of RNN forecasting models.

Figure 8 displays the forecast errors of the four models based on the BiLSTM method integrated with various decomposition methods across the four datasets and forecast steps. Throughout the forecasting period from t + 15 to t + 120, the four forecasting models demonstrate more pronounced errors in spring and summer than in fall and winter. As the forecasting step increases, the predictive capabilities of all models exhibit a declining trend. From spring to winter, the RMSE (MAE) evaluation indices for the BiLSTM model at t + 15 are between 0. 836 and 1.401 m/s (0.589 and 1.07 m/s). Following the implementation of the STL and VMD time series decomposition algorithms, the RMSEs (MAEs) of the STL-BiLSTM and VMD-BiLSTM models are reduced to 0.621–1.058 (0.621–1.05) m/s and 0.439–0.735 (0.323–544) m/s at t + 15. Furthermore, through the synthetic application of the STL and VMD time series decomposition algorithms, the STL-VMD-BiLSTM model possesses a better overall predictive effect than the models mentioned above, with the RMSE and MAE reduced to 0.435–0.668 m/s and 0.321–0.483 m/s, respectively.

At t + 60 and t + 120 (Figure 6 and Figure 7), the STL-VMD-BiLSTM model achieves the best prediction performance, exhibiting the lowest RMSE and MAE scores. At t + 60, the STL-VMD-BiLSTM model records RMSEs from 0.470 to 0.753 m/s and MAEs from 0.347 to 0.567 m/s. The RMSEs and MAEs of the VMD-BiLSTM and STL-BiLSTM models are lower than those of the BiLSTM model. As the prediction step increases, the RMSEs and MAEs of the STL-VMD-BiLSTM model rise to 0.542–0.92 m/s and 0.405–0.696 m/s at t + 120. The BiLSTM model exhibits the worst prediction performance, with the RMSE and MAE values reaching 1.18–1.2.725 m/s and 1.18–2.008 m/s at t + 120.

In summary, the prediction performance of the three models during the periods of t + 15~t + 120 indicates that as the prediction step size increases, the predictive advantage of the STL-VMD-BiLSTM model becomes more significant, effectively capturing the wind speed variability. The synergistic benefits of the dual decomposition technique significantly enhance the short-range wind speed prediction performance of the BiLSTM model.

3.3.4. Prediction Results with Different Models

To further verify the forecasting advantages of the proposed method, machine learning algorithms, including SVR, LGBM, and RF, are combined with the STL, VMD, and STL-VMD time series decomposition techniques. Corresponding comparative tests and analyses are conducted. Figure 8 illustrates the distribution of the accuracy evaluation criteria for the RMSEs and MAEs from t + 15 to t + 120 across the four data subsets.

Table 8. Prediction performance of machine learning models at specified steps in spring (a) and summer (b) sets.

Case	Model	t + 15			t + 60			t + 120
		RMSE	MAE	R2	RMSE	MAE	R2	RMSE	MAE	R2
(a) Spring	STL-VMD-SVR	0.678	0.516	0.960	0.860	0.658	0.931	1.221	0.943	0.840
	VMD-SVR	0.688	0.523	0.959	0.858	0.653	0.932	1.177	0.908	0.855
	STL-SVR	0.942	0.681	0.919	1.727	1.252	0.700	2.576	1.893	0.235
	SVR	1.228	0.894	0.862	2.070	1.510	0.547	2.738	2.030	0.119
	STL-VMD-LGBM	0.909	0.654	0.902	1.080	0.784	0.853	1.382	1.036	0.749
	VMD-LGBM	0.924	0.662	0.897	1.083	0.788	0.852	1.368	1.028	0.753
	STL-LGBM	1.019	0.741	0.873	1.389	0.988	0.754	2.211	1.585	0.270
	LGBM	1.151	0.853	0.850	1.979	1.446	0.497	2.658	1.974	0.042
	STL-VMD-RF	0.743	0.544	0.950	1.013	0.748	0.900	1.428	1.089	0.797
	VMD-RF	0.737	0.537	0.950	1.021	0.756	0.899	1.381	1.052	0.802
	STL-RF	0.877	0.633	0.927	1.358	0.979	0.820	2.258	1.616	0.448
	RF	1.071	0.786	0.900	1.999	1.465	0.613	2.716	2.006	0.255
(b) Summer	STL-VMD-SVR	0.758	0.561	0.927	0.903	0.682	0.889	1.168	0.890	0.796
	VMD-SVR	0.803	0.593	0.916	0.950	0.711	0.874	1.144	0.869	0.804
	STL-SVR	0.983	0.733	0.869	1.597	1.207	0.611	2.058	1.575	0.217
	SVR	1.300	0.963	0.776	1.925	1.453	0.409	2.188	1.678	0.095
	STL-VMD-LGBM	0.915	0.705	0.859	1.005	0.777	0.820	1.175	0.913	0.735
	VMD-LGBM	0.916	0.707	0.859	1.014	0.785	0.815	1.189	0.928	0.722
	STL-LGBM	1.049	0.800	0.800	1.222	0.937	0.713	1.839	1.445	0.242
	LGBM	1.265	0.957	0.736	1.974	1.514	0.225	2.352	1.855	−0.300
	STL-VMD-RF	0.826	0.614	0.911	0.964	0.732	0.869	1.181	0.890	0.787
	VMD-RF	0.828	0.615	0.911	0.965	0.732	0.868	1.180	0.889	0.784
	STL-RF	0.963	0.707	0.874	1.217	0.910	0.777	1.860	1.437	0.387
	RF	1.212	0.889	0.817	1.973	1.461	0.416	2.281	1.750	0.061

and

Table 9. Prediction performance of machine learning models at specified steps in fall (a) and winter (b) sets.

Case	Model	t + 15			t + 60			t + 120
		RMSE	MAE	R2	RMSE	MAE	R2	RMSE	MAE	R2
(a) Fall	STL-VMD-SVR	0.585	0.446	0.962	0.719	0.551	0.940	0.983	0.763	0.877
	VMD-SVR	0.615	0.466	0.958	0.747	0.570	0.935	0.963	0.746	0.883
	STL-SVR	0.738	0.558	0.936	1.304	0.968	0.775	1.963	1.497	0.364
	SVR	1.005	0.746	0.880	1.612	1.198	0.629	2.118	1.632	0.196
	STL-VMD-LGBM	0.781	0.604	0.910	0.874	0.676	0.883	1.078	0.831	0.811
	VMD-LGBM	0.781	0.604	0.910	0.877	0.678	0.882	1.087	0.838	0.804
	STL-LGBM	0.863	0.661	0.885	1.060	0.812	0.817	1.712	1.319	0.429
	LGBM	1.066	0.813	0.826	1.648	1.254	0.492	2.160	1.684	−0.148
	STL-VMD-RF	0.649	0.488	0.952	0.826	0.634	0.916	1.058	0.821	0.853
	VMD-RF	0.650	0.488	0.952	0.823	0.632	0.917	1.067	0.830	0.847
	STL-RF	0.722	0.542	0.940	1.036	0.772	0.865	1.767	1.315	0.528
	RF	0.972	0.719	0.891	1.587	1.183	0.644	2.086	1.598	0.175
(b) Winter	STL-VMD-SVR	0.446	0.327	0.927	0.557	0.424	0.872	0.751	0.589	0.711
	VMD-SVR	0.466	0.341	0.919	0.577	0.438	0.863	0.739	0.577	0.742
	STL-SVR	0.598	0.453	0.860	0.986	0.781	0.552	1.340	1.089	−0.074
	SVR	0.797	0.591	0.757	1.206	0.959	0.302	1.433	1.162	−0.303
	STL-VMD-LGBM	0.604	0.493	0.828	0.674	0.548	0.769	0.806	0.648	0.625
	VMD-LGBM	0.604	0.493	0.828	0.683	0.556	0.762	0.828	0.664	0.601
	STL-LGBM	0.692	0.553	0.761	0.841	0.682	0.614	1.253	1.028	−0.166
	LGBM	0.851	0.660	0.657	1.289	1.051	−0.078	1.545	1.283	−1.155
	STL-VMD-RF	0.471	0.349	0.921	0.594	0.454	0.861	0.758	0.589	0.741
	VMD-RF	0.471	0.350	0.921	0.594	0.454	0.861	0.763	0.597	0.731
	STL-RF	0.583	0.436	0.872	0.758	0.591	0.767	1.186	0.950	0.225
	RF	0.756	0.537	0.794	1.221	0.962	0.286	1.474	1.201	−0.422

show the prediction outcomes of the models at specified steps for these four sets. The RMSEs (MAEs) from t + 15 to t + 120 for the methods based on SVR in conjunction with the time series decomposition algorithms are as follows: the RMSEs (MAEs) of the SVR model range from 0.797 to 2.738 (0.591 to 2.118) m/s, and the RMSEs (MAEs) of the STL-SVR and VMD-SVR models range from 0.598 to 2.576 (0.453 to 1.632) m/s and 0.604 to 1.144 (0.493 to 0.869) m/s, while the STL-VMD-SVR exhibits RMSE (MAE) scores of 0.446–0.751 (0.327–0.589) m/s. The methods combined with time series decomposition outperform the individual approach (SVR).

The RMSEs (MAEs) for the prediction results generated by the hybrid methods based on LGMB and RF from t + 15 to t + 120 are within the ranges of 0.604–2.658 (0.493–1.974) m/s and 0.604–2.716 (0.493–2.006) m/s, respectively. The STL-VMD-LGBM model yields RMSE scores ranging from 0.604 to 1.382 m/s and MAE scores from 0.493 to 1.036 m/s, while the STL-VMD-RF model achieves RMSE scores from 0.471 to 1.428 m/s and MAE scores from 0.393 to 1.089 m/s.

The average RMSE and MAE values for the STL-VMD-SVR, STL-VMD-LGBM and STL-VMD-RF models from t + 15 to t + 120 are 0.617–1.031 (0.463–0.796) m/s, 0.804–1.110 (0.614–0.857) m/s and 0.672–1.105 (0.498–0.847) m/s, respectively. Additionally, the average R² values of the STL-VMD-SVR model range from 0.806 to 0.944, indicating superior performance as well. The contrastive analysis indicates that the STL-VMD-SVR model delivers better prediction performance than the methods based on other machine learning methods.

For all methods, the RMSE and MAE values are higher during spring and summer than in fall and winter (Figure 8a–d). The STL-VMD-BiLSTM model exhibits the best performance, with average RMSE, MAE and R² values ranging from 0.582 to 0.753 m/s, 0.437 to 0.573 m/s and 0.915 to 0.951, respectively, from t + 15 to t + 120. The time cost of the STL-BiLSTM (STL-VMD-BiLSTM) model is two to three (approximately six to seven) times higher than that of BiLSTM, where the time taken to train and predict is 12–18 (42–49) minutes and 0.34 s, respectively, when computed sequentially. Overall, the STL-VMD-BiLSTM model outperforms the hybrid algorithms based on machine learning methods, albeit at a higher computational cost.

4. Conclusions

This paper introduces a hybrid method, STL-VMD-BiLSTM, for the modal decomposition and reconstruction of wind speed time series, taking into account their nonlinear and nonstationary characteristics. This approach aims to enhance the accuracy of short-term forecasts for time intervals of 0 to 2 h by analyzing both temporal and spatial variability features. The comparative analysis of the validation experiment’s results regarding the models’ prediction performance across four data subsets leads to the following conclusions.

A novel model, STL-VMD-BiLSTM, is established by integrating the STL-VMD hybrid decomposition method with BiLSTM, effectively breaking down the wind speed series into high-frequency and low-frequency intrinsic modes, thereby improving the prediction performance for nonlinear variations.

The STL-VMD-BiLSTM model excels in extracting the characteristics of wind speed changes and effectively capturing extreme wind speed forecasts. The STL-VMD-BiLSTM model achieves average RMSE, MAE and R² values ranging from 0.582 to 0.753 m/s, 0.437 to 0.573 m/s and 0.915 to 0.951 from t + 15 to t + 120. On the other hand, the forecasting performance of the BiLSTM model, when not utilizing any time series model reconstruction and decomposition methods, indicates the lowest performance.

To further assess the prediction capabilities and superiority of the proposed methods, SVR, LGBM and RF models were developed for wind speed forecasting in conjunction with the STL-VMD decomposition method. The STL-VMD-BiLSTM model demonstrates superior performance, with RMSEs below 0.92 m/s, although it requires a comparatively longer training time. The average RMSEs (MAEs) for the STL-VMD-SVR, STL-VMD-LGBM and STL-VMD-RF models over the same forecasting periods of t + 15 to t + 120 are 0.617 to 1.031 (0.463 to 0.796) m/s, 0.804 to 1.110 (0.614 to 0.857) m/s and 0.672 to 1.105 (0.498 to 0.847) m/s, respectively, which are significantly higher than those of the proposed method.

5. Limitations and Future Research Directions

The proposed method has shown the superiority of short-term wind speed prediction results with a 2 h forecast horizon. However, there remain a few limitations that are to be addressed in further investigation, which are summarized as follows.

(1)

In this study, the historical period of the wind speed data used to build the model is only one year, and the optimization algorithms are not applied to adjust the parameters and enhance the accuracy of the short-term wind speed forecasting results. In the future, an optimization strategy and longer training and testing sets should be used for prediction and verification.

(2)

This study does not predict the wind speed at a forecast horizon longer than 2 h. In the future, the longer-sequence time series forecasting of the wind speed should be conducted to achieve prediction at longer forecast horizons.