Comparison Study of AI-based Methods in Wind Energy

Abstract

Wind energy forecasting is particularly important for wind farms due to cost-related issues, dispatch planning, and energy market operations. Thus, improving forecasting accuracy becomes an urgent task for researchers in the field of wind energy. However, there is limited research to discuss an overall comparison among various forecasting types, which is a foundation for future works with respect to wind energy prediction because this comparison may reveal whether there is a best model for a specific forecasting type or specific data in this field. For the purpose of laying a strong foundation for wind energy research, this chapter introduces five basic forecasting models, which are Autoregressive Moving Average Model (ARMA), Back-Propagation Neuron Network (BPNN), Support Vector Regression (SVR), Extreme Learning Machine (ELM), and Adaptive Network-Based Fuzzy Inference System (ANFIS) with implement codes before comparing the forecasting effectiveness of five different models in three wind farms based on five forecasting types. Comparison results indicate that each model has great divergent forecasting results in different wind farms and every forecasting type has its own “best model.”

1. Introduction

The important environmental advantages of renewable energy sources have been significantly noticed, which results in most industrialized countries committed to developing the installation of wind power plants [1]. The share of total installed capacity of wind power in China has reached approximately 27% in the global capacity, which is 96.37 GW, due to the historically high installation of new wind power capacity.¹

N.E. Administration, http://www.nea.gov.cn/2015-02/12/c_133989991.htm (2015).

Moreover, a renewable-energy-oriented power system was proposed as the fundamental aim of China’s energy transformation during the Asia-Pacific Economic Cooperation (APEC) Conferences.²

N.E. Administration, http://www.nea.gov.cn/2015-06/08/c_134305870.htm (2015).

Advanced prediction techniques are urgently needed to integrate wind energy into the electrical power grid in a manner that benefits both Transmission System Operators (TSOs) and Independent Power Producers (IPPs) [2,3].

Since wind energy has the inherently intermittent nature and stochastic nonstationarity, it brings significant levels of uncertainty to system operators [4]. Thus, accurate wind forecasts are of primary importance to solve operational, planning, and economic problems in the growing wind power scenario [4,5]. Current wind power forecasting research has been divided into point forecasts (also called deterministic predictions) [6–8] and uncertainty forecasts [9,10]. Deterministic forecasts deliver specific amounts of wind power and focus on reducing the forecasting error [11]. By contrast, it is essential to decision-making processes and electricity market trading strategies [12] that uncertainty forecasts provide uncertainty information for system operators to manage the wind power generation of wind farms [11,13].

The existing approaches published in the literature with respect to wind energy prediction can be divided into three categories—artificial intelligent model, physical model, and statistical model—and sometimes a hybrid one, which integrates advantages of different categories, is involved. Researchers often utilize them to forecast various types of wind speed and wind power by various types including the multi-step forecasting, long-term wind speed (power) forecasting, and so on. However, few literature exist to discuss an overall comparison among various forecasting types, which is a foundation for future works of wind energy prediction because this comparison may reveal whether there is a best model for a specific forecasting type or specific data in the field of wind energy.

The following parts in this paper are demonstrated by four sections: methodology briefly introduces forecasting approaches that this paper adopts; data collection specifically illustrates three types of wind data; simulation and result displays and evaluates the final results of forecasting effectiveness of each model; and conclusion draws the main results that this paper investigates.

2. Methodology

In this section, we will introduce five basic models including Autoregressive Moving Average Model (ARMA), Back-Propagation Neuron Network (BPNN), Support Vector Regression (SVR), Extreme Learning Machine (ELM), and Adaptive Network-Based Fuzzy Inference System (ANFIS) as forecasting methods in this paper to predict wind speed. For the purpose of clearly demonstrating procedures of forecasting wind speed using these methods, this chapter will attach the Matlab code after introduction of each model and the wind speed series is defined as follows:

$Ws=\left(\begin{array}{cccc}w{s}_1,& w{s}_2,& \cdots & ,w{s}_N\end{array}\right)$E1

where N is a positive integer and $w{s}_n\in \left(0,+\infty \right)\subset ℝ$ represents the wind speed of time n.

2.1. Autoregressive moving average model

2.1.1. Brief introduction of ARMA

This section displays ARMA model and the ARMA(p, q) model can be expressed as follows [14]:

$y_t=\delta +{\displaystyle \sum _{i=1}^p{\phi }_i{y}_{t-i}}+{\displaystyle \sum _{j=1}^q{\varphi }_j{e}_{t-j}}+e_t$E2

Where $\delta$ is the constant term of the ARMA model, ${\phi }_i$ is the ith autoregressive coefficient, ${\varphi }_j$ is the jth moving average coefficient, e_t is the error term at time period t, and represents the value of wind speed observed or forecasted at time period t. Thus, the first step in applying ARMA model to forecast wind speed is to identify p and q, which is related to the stationarity of the time series, which means that the stationarity assumption of observed series should be checked first. For this purpose, inspection of the run plots and Auto Correlation Function (ACF) plots can be used for deciding on the order of differencing.

2.1.2. Implement for wind speed forecasting

Since the Matlab has the ARMA package, we will implement this model to forecast wind speed using five steps as follows:

1. Identify the domain of p and q. In this chapter, we set $p\in \left\{1,2,3,4,5\right\}$ and $q\in \left\{1,2,3,4,5\right\}$.

2. For each pair (p, q), calculate corresponding model’s AIC value.

3. Select the pair (p, q), which makes corresponding model’s AIC value to be minimum one among all of pairs of (p, q) as the best (p, q) and denote this pair as (p_test, q_test).

4. Apply (p_test, q_test) to establish ARMA model.

5. Forecast m steps using the model established in step (4).

The detailed Matlab codes are listed in Code 1.

Code 1. Matlab Codes of ARMA.

function [y_output,p_test,q_test] = arma_dynmc(Ws,m)

z = Ws; step = m;

test = [];

for p = 1:5

    for q = 1:5

      m = armax(z,[p,q]);

      AIC = aic(m);

      test = [test;p q AIC];

    end

end

[min_aic,min_i] = min(test(:,3));

p_test = test(min_i,1);q_test = test(min_i,2);

m = armax(z,[p_test,q_test]);

P = predict(m,z,step);

y_output = P(end-step+1:end,1)';

2.2. Back-propagation neuron network

2.2.1. Brief introduction of BPNN

Mccelland and Rumelhart developed the BP neural network model in 1985. There are three layers in a particular network: the input, hidden, and output layers. Each layer has designed nodes, whose functions aim to calculate the inner product of the input vector and weight vector by the transfer function [15]. The process of BPNN is demonstrated by the following steps:

1. Initialize. Assume that the input layer has n nodes, and hidden layer has l nodes, and output layer has m nodes. Let the network distribute values randomly for each threshold value θ_j, γ_t and the connection weight w_ij, v_jt, i = 1, 2, …, n, j = 1, 2, …, l, and t = 1, 2, …, m.

2. Calculate the output of hidden layer. Assume X$\subset {\left[-1,1\right]}^n$ is the input space and Y$\subset {\left[-1,1\right]}^m$ is the output space, which means X is an n-dimensional space and Y is an m-dimensional space. We denote each element in as x_i and each element in Y is y_t. Thus, the output of hidden layer can be calculated using the following function:

   $H_j=f\left({\displaystyle \sum _{i=1}^n{w}_{ij}{x}_i}-{\theta }_j\right),j=1,2,\mathrm{...},l$E3
   $f(•)$ is the transfer function and is expressed as follows:
   $f\left(z\right)=\frac{1}{1+e^{-z}}$

3. Calculate the output of the network. Using H_j, v_jk, and γ_t, the output of the network can be calculated as follows:

   $O_t={\displaystyle \sum _{j=1}^l{v}_{jk}{H}_j}-{\gamma }_t,t=1,2,\mathrm{...},m$E4

4. Error calculation. The error between y_t and O_t can be expressed as e_t = O_t – y_t.

5. Weights updating. According to back-propagation algorithm, the weights w_ij, v_jt can be updated by using the following functions:

   $w_{ij}=w_{ij}+\eta H_j\left(1-H_j\right)x_i{\displaystyle \sum _{t=1}^m{w}_{jt}{e}_t}$E5
   $v_{jt}=v_{jt}+\eta H_j{e}_t$E6
   where $\eta$ is the learning rate and i = 1, 2, …, n, j = 1, 2, …, l, and t = 1, 2, …, m.

6. Biases updating. Similarly, the biases θ_j, γ_t can be updated by using the following functions:

   ${\theta }_j={\theta }_j+\eta H_j\left(1-H_j\right){\displaystyle \sum _{t=1}^m{w}_{jt}{e}_t}$E7
   ${\gamma }_t={\gamma }_t+e_t$E8

7. Iteration. If the error is more than the expected value, then execute steps (2)–(6). Otherwise, denote θ_j, γ_t, w_ij, and v_jt as the optimized parameter of this network with respect to X and Y.

2.2.2. Implement for wind speed forecasting

In this part, how to apply BPNN to forecast wind speed will be introduced in detail and the MATLAB code for predicting wind speed will be attached. There are six steps to forecast wind speed applying BPNN as follows:

1. Design forecasting mode. In this chapter, AI-based model will apply the following mode to forecast wind speed of different wind databases

   $\left(w{s}_{t-1},w{s}_{t-2}\right)\stackrel{forecast}{\to }w{s}_t$E9
   Eq. (8) indicates that the input space X for BPNN model is a two-dimensional space and the output space Y is a one-dimensional space, meaning that n = 2 and m = 1 as described in Section 2.2.1.

2. Assign training samples. Based on step (1), the input and output of training samples are expressed as follows:

   $Train\_input=\left(\begin{array}{cccc}w{s}_1& w{s}_2& \cdots & w{s}_{N-2}\\ w{s}_2& w{s}_3& \cdots & w{s}_{N-1}\end{array}\right)$E10
   $Train\_output=\left(\begin{array}{cccc}w{s}_3& w{s}_4& \cdots & ws{}_N\end{array}\right)$E11

   The Matlab codes of this step are listed in Code 2.

   Code 2. Matlab codes to construct training samples.

   n = length(Ws);

   for i = 1:2

       Train_input(i,:) = Ws(i:n+i-3);

   End

   Train_output = Ws(3:n)

3. Normalize Train_input and Train_output. Since $X\subset {\left[-1,1\right]}^n$ and $Y\subset {\left[-1,1\right]}^m$, Train_input and Train_output need to be adjusted to satisfy this condition. First, the maximum and minimum values of each row in Train_input need to be calculated using the following functions:

   $Mi{n}_{input}^{\left(i\right)}=\min \left\{i\text{th}\text{row}\text{in}Train\_input\right\}$E12
   $Ma{x}_{input}^{\left(i\right)}=\max \left\{i\text{th}\text{row}\text{in}Train\_input\right\}$E13

   Let $D_{input}^{\left(i\right)}=Ma{x}_{input}^{\left(i\right)}-Mi{n}_{input}^{\left(i\right)}$ and normalized Train_input can be obtained by the following equation (i = 1, 2):

   $Train\_inputm=\left(\begin{array}{cccc}\frac{w{s}_1-Mi{n}_{input}^{\left(1\right)}}{D_{input}^{\left(1\right)}}& \frac{w{s}_2-Mi{n}_{input}^{\left(1\right)}}{D_{input}^{\left(1\right)}}& \cdots & \frac{w{s}_{N-2}-Mi{n}_{input}^{\left(1\right)}}{D_{input}^{\left(1\right)}}\\ \frac{w{s}_2-Mi{n}_{input}^{\left(2\right)}}{D_{input}^{\left(2\right)}}& \frac{w{s}_3-Mi{n}_{input}^{\left(2\right)}}{D_{input}^{\left(2\right)}}& \cdots & \frac{w{s}_{N-1}-Mi{n}_{input}^{\left(2\right)}}{D_{input}^{\left(2\right)}}\end{array}\right)$E14

   Similarly, normalized Train_output can be obtained by the following equations:

   $Mi{n}_{output}^{\left(i\right)}=\min \left\{i^{th}\text{row}\text{in}Train\_output\right\}$
   $Ma{x}_{output}^{\left(i\right)}=\max \left\{i^{th}\text{row}\text{in}Train\_output\right\}$
   $D_{output}^{\left(i\right)}=Ma{x}_{output}^{\left(i\right)}-Mi{n}_{output}^{\left(i\right)}$
   $Train\_outputm=\left(\begin{array}{cccc}\frac{w{s}_3-Mi{n}_{output}^{\left(1\right)}}{D_{output}^{\left(1\right)}}& \frac{w{s}_4-Mi{n}_{output}^{\left(1\right)}}{D_{output}^{\left(1\right)}}& \cdots & \frac{w{s}_N-Mi{n}_{output}^{\left(1\right)}}{D_{output}^{\left(1\right)}}\end{array}\right)$E15

   The Matlab codes of this step are listed in Code 3.

   Code 3. Matlab codes to normalize training samples.

   [Train_inputm, inputs] = mapminmax(Train_input);

   [Train_outputm, outputs] = mapminmax(Train_output);

4. Assign parameters and train the network. We assign that the network has three layers and the input layer has two nodes, the hidden layer has two nodes and the output layer has one node. Thus, according to Section 2.2.1, the optimized weights and biases of this network can be calculated and denoted by θ_j, γ_t, w_ij, v_jt, i = 1, 2, j = 1, 2, and t = 1. The Matlab codes of this step are listed in Code 4.

   Code 4. Matlab codes to build and train BPNN.

   Net = newff(Train_input, Train_output, 2);

   Net = train(Net, Train_input, Train_output);

   Net is the trained network and can be applied to forecast wind speed in the following steps:

5. Forecast $w{s}_{N+1}$. Based on the forecasting mode in step (1), we need to apply $w{s}_{N-1}$ and $w{s}_N$ to forecast $w{s}_{N+1}$, which means that the Test_input is ${\left(w{s}_{N-1},w{s}_N\right)}^T$. Then, $D_{input}^{\left(i\right)},Ma{x}_{input}^{\left(i\right)}\text{and}Mi{n}_{input}^{\left(i\right)}$ obtained in step (3) will be applied to normalize Test_input and the normalized Test_input can be expressed using Test_inputm as follows:

   $Test\_inputm=\left(\begin{array}{c}\frac{w{s}_{N-1}-Mi{n}_{input}^{\left(1\right)}}{D_{input}^{\left(1\right)}}\\ \frac{w{s}_N-Mi{n}_{input}^{\left(2\right)}}{D_{input}^{\left(2\right)}}\end{array}\right)$E16

   Applying Eqs. (3) and (4) in Section 2.2.1, we can obtain the output of the network and denote it by $\mathit{O}\in ℝ$. Thus, the forecasting value of $w{s}_{N+1}$ is Ws_forecast = $\mathit{O}\times {\mathit{D}}_{output}^{\left(1\right)}+\mathit{M}i{n}_{output}^{\left(1\right)}$.

   The Matlab codes of this step are listed in Code 5.

   Code 5. Matlab codes to forecast using established BPNN (single step ahead)

   Test_inputm = mapminmax(‘apply’, Test_input, inputs);

   O=sim(Net, Test_inputm);

   Ws_forecast= mapminmax(‘reverse’, O, outputs);

6. Forecast $w{s}_{N+2}$. Based on the forecasting mode in step (1), we need to apply $w{s}_N$ and forecasting value of $w{s}_{N+1}$, Ws_forecast, to forecast $w{s}_{N+2}$, which means that the Test_input2 is ${\left(w{s}_N,Ws\_forecast\right)}^T$. Then, $D_{input}^{\left(i\right)},Ma{x}_{input}^{\left(i\right)}\text{and}Mi{n}_{input}^{\left(i\right)}$ obtained in step (3) will be applied to normalize Test_input2 and the normalized Test_input2 can be expressed using Test_inputm2 as follows:

   $Test\_inputm2=\left(\begin{array}{c}\frac{w{s}_N-Mi{n}_{input}^{\left(1\right)}}{D_{input}^{\left(1\right)}}\\ \frac{Ws\_forecast-Mi{n}_{input}^{\left(2\right)}}{D_{input}^{\left(2\right)}}\end{array}\right)$E17

   Applying Eqs. (3) and (4) in Section 2.2.1, we can obtain the output of the network and denote it by $O_2\in ℝ$. Therefore, the forecasting value of $w{s}_{N+2}$ is Ws_forecast2 = $O_2\times D_{output}^{\left(1\right)}+Mi{n}_{output}^{\left(1\right)}$.

   The Matlab codes of this step are listed in Code 6.

   Code 6. Matlab codes to forecast using established BPNN (two-step ahead).

   Test_inputm2 = mapminmax(‘apply’, Test_input2, inputs);

   O2=sim(Net, Test_inputm2);

   Ws_forecast2= mapminmax(‘reverse’, O, outputs);

2.3. Support vector regression

2.3.1. Brief introduction of SVR

SVR is the most common application of support vector machines (SVMs) and constructs a hyperplane that separates examples with maximum margin, to categorize or forecast series [16]. Given the training data $\left\{\left(x_1,y_1\right),\dots ,\left(x_l,y_l\right)\right\}\subset \mathbf{{\rm X}}\times ℝ$, where $\mathbf{{\rm X}}={ℝ}^d$ denotes the space of input patterns, then the goal of SVR is to find a function $f\left(x\right)$ that is as flat as possible with at most $\epsilon$ deviation from the actually obtained targets $y_i$ for all the training data. The simplest, linear formula for the output of a linear SVR is defined as [17]

$f\left(x\right)=⟨w,x⟩+b,w\in X,b\in ℝ$E18

where $⟨w,x⟩$ denotes the dot product in ${\rm X}$. One way to get the largest flatness is to minimize the $w$ in Eq. (18) and this problem can be written as Eq. (19) by introducing slack variables ${\xi }_i$, ${\xi }_i^{*}$ into a convex optimization problem function:

$\begin{array}{l}minimize\frac{1}{2}{\Vert w\Vert }^2+C{\displaystyle {\sum }_{i=1}^l\left({\xi }_i+{\xi }_i^{*}\right)}\\ subjectto\{\begin{array}{l}{y}_i-⟨w,x_i⟩-b\le \epsilon \\ ⟨w,x_i⟩+b-y_i\le \epsilon \\ {\xi }_i,{\xi }_i^{*}\ge 0\end{array}\end{array}$E19

The constant $C>0$ determines the trade-off between the flatness of f and the amount up to which deviations larger than $\epsilon$ are tolerated.

To figure out the support vector regression function, the dual problem of Eq. (19) is defined in Eq. (20) by constructing a Lagrange function from the objective function:

$\begin{array}{l}maximize\{\begin{array}{l}\frac{1}{2}{\displaystyle {\sum }_{i,j=1}^l\left({\alpha }_i-{\alpha }_i^{*}\right)\left({\alpha }_j-{\alpha }_j^{*}\right)⟨x_i,x_j⟩}\\ -\epsilon {\displaystyle {\sum }_{i=1}^l\left({\alpha }_i+{\alpha }_i^{*}\right)+{\displaystyle {\sum }_{i=1}^l{y}_i\left({\alpha }_i-{\alpha }_i^{*}\right)}}\end{array}\\ subjectto{\displaystyle {\sum }_{i=1}^l\left({\alpha }_i-{\alpha }_i^{*}\right)=0and{\alpha }_i,{\alpha }_i^{*}\in \left[0,C\right]}\end{array}$E20

With the solution of $\left(\alpha ,{\alpha }^{*}\right)$ from Eq. (20), the support vector regression function can be written as follows:

$f\left(x\right)={\displaystyle {\sum }_{i=1}^l\left({\alpha }_i-{\alpha }^{*}\right)⟨x_i,x⟩+b}$E21

When it comes to the nonlinear regression problem, the training patterns $x_i$ can be mapped into a high-dimensional space, where the nonlinear regression problem is transformed into a linear one. The expansion of Eq. (21) is defined as Eq. (22):

$f\left(x\right)={\displaystyle {\sum }_{i=1}^N\left({\alpha }_i^{*}-{\alpha }_i\right)k\left(x_i,x;g\right)+b}$E22

where ${\alpha }_i^{*}$ and ${\alpha }_i$ are Lagrange multipliers and $k\left(x_i,x;g\right)$ is the kernel function, in which g is a parameter and generally set as 1/d.

2.3.2. Implement for wind speed forecasting

We use the libsvm package (Version 3.17) to implement forecasting task of wind speed using SVR. There are six steps and steps (1)–(3) are same as in steps (1)–(3) in Section 2.2.2. Thus, we start to introduce from step (4).

(4) Assign parameters and train the SVR model. We assign that C is 4 and g is 0.5. Thus, according to Section 2.3.1, the optimized solution $\left(\alpha ,{\alpha }^{*}\right)$ can be calculated. The Matlab codes of this step are listed in Code 7.

Code 7. Matlab codes to establish and train SVR.

Model = svmtrain(Train_output’, Train_input’, ‘-C 4 –g 0.5’);

Model is the trained SVR model and can be applied to forecast wind speed by the following steps:

(5) Forecast $w{s}_{N+1}$. Based on the forecasting mode in step (1), we need to apply $w{s}_{N-1}$ and $w{s}_N$ to forecast $w{s}_{N+1}$, which means that the Test_input is ${\left(w{s}_{N-1},w{s}_N\right)}^T$. Then, $D_{input}^{\left(i\right)},Ma{x}_{input}^{\left(i\right)}\text{and}Mi{n}_{input}^{\left(i\right)}$ obtained in step (3) will be applied to normalize Test_input and the normalized Test_input can be expressed using Test_inputm as follows:

$Test\_inputm=\left(\begin{array}{c}\frac{w{s}_{N-1}-Mi{n}_{input}^{\left(1\right)}}{D_{input}^{\left(1\right)}}\\ \frac{w{s}_N-Mi{n}_{input}^{\left(2\right)}}{D_{input}^{\left(2\right)}}\end{array}\right)$E23

Applying Eqs. (3) and (4) in Section 2.2.1, we can obtain the output of the network and denote it by $O\in ℝ$. Thus, the forecasting value of $w{s}_{N+1}$ is Ws_forecast = $O\times D_{output}^{\left(1\right)}+Mi{n}_{output}^{\left(1\right)}$. The Matlab codes of this step are listed in Code 8.

Code 8. Matlab codes to forecast using established SVR (single step ahead).

Test_inputm = mapminmax(‘apply’, Test_input, inputs);

O=svmpredict(0, Test_inputm’,Model);

Ws_forecast= mapminmax(‘reverse’, O, outputs);

(6) Forecast $w{s}_{N+2}$. Based on the forecasting mode in step (1), we need to apply $w{s}_N$ and the forecasting value of $w{s}_{N+1}$, Ws_forecast, to forecast $w{s}_{N+2}$, which means that the Test_input2 is ${\left(w{s}_N,Ws\_forecast\right)}^T$. Then, $D_{input}^{\left(i\right)},Ma{x}_{input}^{\left(i\right)}\text{and}Mi{n}_{input}^{\left(i\right)}$ obtained in step (3) will be applied to normalize Test_input2 and the normalized Test_input2 can be expressed using Test_inputm2 as follows:

$Test\_inputm2=\left(\begin{array}{c}\frac{w{s}_N-Mi{n}_{input}^{\left(1\right)}}{D_{input}^{\left(1\right)}}\\ \frac{Ws\_forecast-Mi{n}_{input}^{\left(2\right)}}{D_{input}^{\left(2\right)}}\end{array}\right)$E24

Applying Eqs. (3) and (4) in Section 2.2.1, we can obtain the output of the network and denote it by $O_2\in ℝ$. Therefore, the forecasting value of $w{s}_{N+2}$ is Ws_forecast2 = $O_2\times D_{output}^{\left(1\right)}+Mi{n}_{output}^{\left(1\right)}$. The Matlab codes of this step are listed in Code 9.

Code 9. Matlab codes to forecast using established SVR (two-step ahead).

Test_inputm2 = mapminmax(‘apply’, Test_input2, inputs);

O2=svmpredict(0,Test_inputm2’,Model);

Ws_forecast2= mapminmax(‘reverse’, O, outputs);

2.4. Extreme learning machine

2.4.1. Brief introduction of ELM

Huang et al. [18] investigated the ELM, which is an effective and efficient learning algorithm. The ELM aims to randomly initialize the weights and biases of SLFN and then to explicitly calculate the hidden layer output matrix and hence the output weights. Because of the nature of non-adjusted weights and bias, the network can be established using a very low computational cost [19]. Then, the ELM with l hidden neurons and transfer function $\phi (·)$ can approximate the N samples with zero error as

${\displaystyle \sum _{j=1}^l{\beta }_j\phi \left(w_j{x}_k+b_j\right)=y_k,k=1,2,\mathrm{...},N}$E25

which can be written as HB=Y, with

$H\left(w_1,\mathrm{...},w_l,b_1,\mathrm{...},b_l,x_1,\mathrm{...},x_N\right)=\left(\begin{array}{ccc}\phi \left(w_1{x}_1+b_1\right)& \cdots & \phi \left(w_l{x}_1+b_l\right)\\ ⋮& \ddots & ⋮\\ \phi \left(w_1{x}_N+b_1\right)& \cdots & \phi \left(w_l{x}_N+b_l\right)\end{array}\right)$E26

where w_i is the weight vector connecting the ith hidden neuron and the input nodes, ${\beta }_i$ is the weight vector connecting the ith hidden neuron and the output neurons, and $b_i$ is the threshold of the ith hidden neuron, and $\mathit{B}={\left({\beta }_1^T\dots {\beta }_l^T\right)}^T,\mathit{Y}={\left(y_1^T\dots y_N^T\right)}^T$

Then, the output weights B can be calculated from the hidden layer output matrix H and the target values Y as B^ = H†Y, where H† is the Moore-Penrose generalized inverse of the matrix H.

2.4.2. Implement for wind speed forecasting

We use Matlab to compile ELM model and provide two *.m files, which are Elmtrain.m and Elmpredict.m, in Appendix. Elmtrain function is similar to train function for BPNN and svmtrain function for SVR and Elmpredict function is similar to sim function for BONN and svmpredict function for SVR. In detail, there are six steps for forecasting wind speed using ELM model and steps (1)–(3) are same as in steps (1)(3) in Section 2.2.2. Thus, we start to introduce them from step (4).

(4) Assign parameters and train the ELM model. We design that the input layer of ELM model has two nodes, and the hidden layer of ELM has two nodes, and the output layer of ELM model has one node. Thus, according to Section 2.4.1, the optimized B can be calculated. The Matlab codes of this step are listed in Code 10.

Code 10. Matlab codes to establish and train ELM.

[W,b,B,TF,TYPE] = elmtrain(Train_input, Train_output, 2);

W is the weight matrix of the input layer and hidden layer, and b is the bias vector of input layer and hidden layer, and B is B^ in Section 2.4.1, and TF is the transfer function $\phi (·)$, and TYPE is model’s functions (classification or regression). We select sigmoid function as the transfer function $\phi (·)$ of ELM, which is same as that of BPNN.

(5) Forecast $w{s}_{N+1}$. Based on the forecasting mode in step (1), we need to apply $w{s}_{N-1}$ and $w{s}_N$ to forecast $w{s}_{N+1}$, which means that the Test_input is ${\left(w{s}_{N-1},w{s}_N\right)}^T$. Then, $D_{input}^{\left(i\right)},Ma{x}_{input}^{\left(i\right)}\text{and}Mi{n}_{input}^{\left(i\right)}$ obtained in step (3) will be applied to normalize Test_input and the normalized Test_input can be expressed using Test_inputm as follows:

$Test\_inputm=\left(\begin{array}{c}\frac{w{s}_{N-1}-Mi{n}_{input}^{\left(1\right)}}{D_{input}^{\left(1\right)}}\\ \frac{w{s}_N-Mi{n}_{input}^{\left(2\right)}}{D_{input}^{\left(2\right)}}\end{array}\right)$E27

Applying Eqs. (3) and (4) in Section 2.2.1, we can obtain the output of the network and denote it by $O\in ℝ$. Thus, the forecasting value of $w{s}_{N+1}$ is Ws_forecast = $O\times D_{output}^{\left(1\right)}+Mi{n}_{output}^{\left(1\right)}$. The Matlab codes of this step are listed in Code 11.

Code 11. Matlab codes to forecast using established ELM (single step ahead).

Test_inputm = mapminmax(‘apply’, Test_input, inputs);

O=elmpredict(Test_inputm, W,b,B,TF,TYPE);

Ws_forecast= mapminmax(‘reverse’, O, outputs);

(6) Forecast $w{s}_{N+2}$. Based on the forecasting mode in step (1), we need to apply $w{s}_N$ and forecasting value of $w{s}_{N+1}$, Ws_forecast, to forecast $w{s}_{N+2}$, which means that the Test_input2 is ${\left(w{s}_N,Ws\_forecast\right)}^T$. Then, $D_{input}^{\left(i\right)},Ma{x}_{input}^{\left(i\right)}\text{and}Mi{n}_{input}^{\left(i\right)}$ obtained in step (3) will be applied to normalize Test_input2 and the normalized Test_input2 can be expressed using Test_inputm2 as follows:

$Test\_inputm2=\left(\begin{array}{c}\frac{w{s}_N-Mi{n}_{input}^{\left(1\right)}}{D_{input}^{\left(1\right)}}\\ \frac{Ws\_forecast-Mi{n}_{input}^{\left(2\right)}}{D_{input}^{\left(2\right)}}\end{array}\right)$E28

Applying Eqs. (3) and (4) in Section 2.2.1, we can obtain the output of the network and denote it by $O_2\in ℝ$. Therefore, the forecasting value of $w{s}_{N+2}$ is Ws_forecast2 = $O_2\times D_{output}^{\left(1\right)}+Mi{n}_{output}^{\left(1\right)}$. The Matlab codes of this step are listed in Code 12.

Code 12. Matlab codes to forecast using established ELM (two-step ahead).

Test_inputm2 = mapminmax(‘apply’, Test_input2, inputs);

O2=svmpredict(Test_inputm2, W,b,B,TF,TYPE);

Ws_forecast2= mapminmax(‘reverse’, O, outputs);

2.5. Adaptive network-based fuzzy inference system

2.5.1. Brief introduction of ANFIS

ANFIS is introduced to compensate for the disability of conventional mathematical tools to address uncertain systems, such as human knowledge and reasoning processes. There are two contributions of FIFs restructured: proposing a standard method for transforming ill-defined factors into identifiable rules of FIS and using an adaptive network to tune the membership functions. We assume that there are two fuzzy if-then rules contained in the system, two inputs (x and y) and one output (z), and the processes of ANFIS are described in Figure 1 [20–21].

Figure 1.

Five processes included in this figure. In an ANFIS architecture, circles are fixed nodes without parameters and squares represent adaptive nodes whose parameters are determined by training data and a gradient-based learning procedure.

Layer I: Mapping a certain input x to a fuzzy set $O_i^{\left(1\right)}$ for every node i by the member functions ${\mu }_A$, which is usually bell-shaped with a parameter set $\left\{a_i,b_i,c_i\right\}$, as is y

$O_i^{\left(1\right)}={\mu }_A\left(x\right),\text{where}{\mu }_A\left(x\right)=\frac{1}{1+{\left[\left(\frac{x-c_i}{a_i}\right)\right]}^{b_i}}\text{or}{\mu }_A\left(x\right)=e^{-{\left(\frac{x-c_i}{a_i}\right)}^2}$E29

Layer II: In this layer, each circle node performs the connection “AND” and multiplies inputs, as well as sends the product out:

$O_i^{\left(2\right)}={\omega }_i={\mu }_A\left(x\right)\times {\mu }_B\left(x\right)$E30

Layer III: Every circle node in this layer calculates a normalized firing strength, namely a ratio of the rule’s firing strength to the sum of all rules’ firing strengths:

$O_i^{\left(3\right)}={\overline{\omega }}_i=\frac{{\omega }_i}{{\displaystyle \sum {\omega }_i}}$E31

Layer IV: Assume the rules of this system are as follows:

Rule 1: If x is A₁ and y is B₁, then $f_1=p_{1}x+q_{1}y+r_1$

Rule 2: If x is A₂ and y is B₂, then $f_2=p_{2}x+q_{2}y+r_2$

Then, the outputs of the adaptive nodes in this layer are computed by

$O_i^{\left(4\right)}={\overline{\omega }}_i{f}_i={\overline{\omega }}_i\left(p_{1}x+q_{1}y+r_1\right)$E32

Layer V: The overall output is the weighted average of all incoming signals:

$O_i^{\left(5\right)}={\displaystyle \sum _i{\overline{\omega }}_i{f}_i}=\frac{{\displaystyle {\sum }_i{\omega }_i{f}_i}}{{\displaystyle {\sum }_i{\omega }_i}}$E33

Particularly, in this case,

$O_i^{\left(5\right)}={\displaystyle \sum _{i=1}^2{\overline{\omega }}_i{f}_i}={\displaystyle \sum _{i=1}^2\left({\overline{\omega }}_{i}x\right)p_i+\left({\overline{\omega }}_{i}y\right)q_i+{\overline{\omega }}_i{r}_i}$E34

2.5.2. Implement for wind speed forecasting

We use the genfis3 function and ANFIS function in Matlab to implement forecasting task of wind speed using ANFIS. There are six steps and steps (1)–(3) are same as in steps (1)–(3) in Section 2.2.2. Thus, we start to introduce from step (4).

(4) Assign parameters and train the ANFIS model. We set the number of iteration of ANFIS as 100. Thus, the optimized parameters can be calculated. The Matlab codes of this step are listed in Code 13.

Code 13. Matlab codes to establish and train ANFIS.

fismat = genfis3(Train_input’, Train_output’);

out_fis1 = anfis([Train_input’ Train_output’],fismat,100);

The out_fis1 is the trained ANFIS and can be applied to forecast wind speed.

(5) Forecast $w{s}_{N+1}$. Based on the forecasting mode in step (1), we need to apply $w{s}_{N-1}$ and $w{s}_N$ to forecast $w{s}_{N+1}$, which means that the Test_input is ${\left(w{s}_{N-1},w{s}_N\right)}^T$. Then, $D_{input}^{\left(i\right)},Ma{x}_{input}^{\left(i\right)}\text{and}Mi{n}_{input}^{\left(i\right)}$ obtained in step (3) will be applied to normalize Test_input and the normalized Test_input can be expressed using Test_inputm as follows:

$Test\_inputm=\left(\begin{array}{c}\frac{w{s}_{N-1}-Mi{n}_{input}^{\left(1\right)}}{D_{input}^{\left(1\right)}}\\ \frac{w{s}_N-Mi{n}_{input}^{\left(2\right)}}{D_{input}^{\left(2\right)}}\end{array}\right)$E35

Applying Eqs. (3) and (4) in Section 2.2.1, we can obtain the output of the network and denote it by $O\in ℝ$. Thus, the forecasting value of $w{s}_{N+1}$ is Ws_forecast = $O\times D_{output}^{\left(1\right)}+Mi{n}_{output}^{\left(1\right)}$. The Matlab codes of this step are listed in Code 14.

Code 14. Matlab codes to forecast using established ANFIS (single step ahead).

Test_inputm = mapminmax(‘apply’, Test_input, inputs);

O=evalfis(Test_inputm’, out_fis1);

Ws_forecast= mapminmax(‘reverse’, O, outputs);

(6) Forecast $w{s}_{N+2}$. Based on the forecasting mode in step (1), we need to apply $w{s}_N$ and forecasting value of $w{s}_{N+1}$, Ws_forecast, to forecast $w{s}_{N+2}$, which means that the Test_input2 is ${\left(w{s}_N,Ws\_forecast\right)}^T$. Then, $D_{input}^{\left(i\right)},Ma{x}_{input}^{\left(i\right)}\text{and}Mi{n}_{input}^{\left(i\right)}$ obtained in step (3) will be applied to normalize Test_input2 and the normalized Test_input2 can be expressed using Test_inputm2 as follows:

$Test\_inputm2=\left(\begin{array}{c}\frac{w{s}_N-Mi{n}_{input}^{\left(1\right)}}{D_{input}^{\left(1\right)}}\\ \frac{Ws\_forecast-Mi{n}_{input}^{\left(2\right)}}{D_{input}^{\left(2\right)}}\end{array}\right)$E36

Applying Eqs. (3) and (4) in Section 2.2.1, we can obtain the output of the network and denote it by $O_2\in ℝ$. Therefore, the forecasting value of $w{s}_{N+2}$ is Ws_forecast2 = $O_2\times D_{output}^{\left(1\right)}+Mi{n}_{output}^{\left(1\right)}$. The Matlab codes of this step are listed in Code 15.

Code 15. Matlab codes to forecast using established ANFIS (two-step ahead).

Test_inputm2 = mapminmax(‘apply’, Test_input2, inputs);

O2=evalfis(Test_inputm2’, out_fis1);

Ws_forecast2= mapminmax(‘reverse’, O, outputs);

3. Data collection

We select three types of wind data to train and test forecasting models, which are 10-minute wind speed data in wind farm A (this database is denoted by WFD1), 15-minute wind speed data in wind farm B (this database is denoted by WFD2), and 15-minute wind speed data in wind farm C (this database is denoted by WFD3), and wind farms B and C are in the same city in China.

Figure 2.

Descriptive information and histograms of three wind farm databases.

In this section, the descriptive information of each database will be provided and we will test the similarity between WFD2 and WFD3 using Friedman test. In the beginning, Figure 2 shows the three databases and their histograms.

From Figure 2, this is apparent that WFD1 and WFD2 have similar maximum and minimum values because they are in the same city. However, WFD1 and WFD2 have different mean and variance values. Thus, we use Friedman test to figure out whether both of two time series are significantly different, which is shown in Table 1. In Table 1, the p-value of Friedman test is very small (2.46e–11), which indicates that WFD1 and WFD2 have a significant difference.

Source	SS	df	MS	Chi-squared	Prob>Chi-squared
Columns	73.3	1	73.2825	44.57	2.46e–11
Interaction	61568.7	17,085	3.6037
Error	22,649	34,172	0.6628
Total	84,291	68,343

Table 1.

The result of Friedman test between WFD1 and WFD2.

For WFD3, its mean value is significantly different from others and the histogram of WFD3 is also different from those of other databases. From the description of these three databases, it can be seen that they are very different, which lays a strong foundation for the comparison study.

4. Simulation and results

For these different databases, we will test 1-h-ahead, 4-h-ahead, 1-h average, 4-h average, and 1-day average wind speed forecasting effectiveness by using five methods separately, which is to say that we conduct five experiments of different time scales or time horizons in this comparison study. To test the overall forecasting effectiveness of each model, the testing period is 30 weeks in 2014 and we apply three criteria to evaluate forecasting error, which are the mean absolute percent error (MAPE), root-mean-squared error (RMSE), and mean absolute error (MAE). The formulas of three criteria are expressed as follows:

$e_t={\widehat{y}}_t-y_t$E37

$MAPE=\frac{1}{N}{\displaystyle \sum _{t=1}^N\frac{\left|e_t\right|}{y_t}}\times 100\%$E38

$RMSE=\frac{1}{N}\sqrt{{\displaystyle \sum _{t=1}^N{e}_t^2}}$E39

$MAE=\frac{1}{N}{\displaystyle \sum _{t=1}^N\left|e_t\right|}$E40

where t represents the time t, $y_t$ is the actual value in time t, ${\widehat{y}}_t$ is the forecasting value in time t, and N is the length of testing period.

4.1. Experiment I: 1-h-ahead wind speed forecasting

In this experiment, we will test the 1-h-ahead forecasting effectiveness of five models. It should be pointed that 1-h-ahead prediction means a four-step forecast for WFD1 and WFD2 and a means six-step forecast for WFD3 because WFD3 is constructed by 10-min wind speeds. The main results are shown in Table 2, which contains MAPE, RMSE, and MAE of each model when forecasting 1-h-ahead wind speed.

Wind farms	Criteria	ANFIS	BPNN	ELM	SVM	ARMA
WFD1	MAE	1.4193	0.8066	0.8335	0.9471	1.4293
	MAPE	43.67%	19.25%	19.51%	30.74%	30.96%
	RMSE	1.5092	0.9083	0.9375	1.0521	1.5574
WFD2	MAE	1.8110	0.7651	0.7649	0.8590	2.9846
	MAPE	49.82%	20.21%	19.76%	28.11%	83.01%
	RMSE	1.8837	0.8611	0.8607	0.9550	3.0937
WFD3	MAE	1.2618	0.9781	0.8571	0.9055	3.0764
	MAPE	25.22%	20.26%	16.29%	19.05%	64.92%
	RMSE	1.3735	1.1019	0.9797	1.0290	3.1724

Table 2.

The forecasting results of Experiment I.

From Table 2, it is apparent that BPNN has the best performance among these five models in WFD1 and ELM outperforms others in WFD2 and WFD3 by three criteria. By contrast, the forecasting accuracy of ANFSI and ARMA is poor and even cannot be accepted in some wind farms, such as the performance of ARMA in WFD1 and WFD2 and the performance of ANFIS in WFD1 and WFD2. These results show that each wind farm has the corresponding best forecasting model and there is no single model that can perform best in every wind farm by three criteria.

4.2. Experiment II: 4-h-ahead wind speed forecasting

Experiment II shows the 4-h-ahead forecasting effectiveness of wind speed, and the 4-h-ahead forecasting is equal to a 16-step forecast for WFD1 and WFD2 and is a 24-step forecast for WFD3. Table 3 descripts forecasting results of five models in three wind farms.

In Table 3, BPNN is the best model in WFD1 and WFD2 by three criteria and outperforms other models in three databases by MAE and RMSE. ANFIS has the lowest MAPE in WFD3 among these models but has high MAPE, MAE, and RMSE in WFD1 and WFD2. It is similar to Experiment I that ARMA does not have a decent forecasting results in three wind farms. These results indicate that models forecasting effectiveness varies a lot when the wind farm changes (such as the effectiveness of ANFIS).

Wind farms	Criteria	ANFIS	BPNN	ELM	SVM	ARMA
WFD1	MAE	1.9441	1.4385	1.4725	1.6147	2.4799
	MAPE	54.41%	35.36%	36.67%	52.28%	55.59%
	RMSE	2.1531	1.6664	1.7062	1.8376	2.7478
WFD2	MAE	1.9937	1.3251	1.3492	1.4246	2.8887
	MAPE	54.96%	34.71%	36.36%	45.85%	78.64%
	RMSE	2.1770	1.5440	1.5647	1.6301	3.1138
WFD3	MAE	1.3250	1.3176	1.3575	1.3951	2.6248
	MAPE	26.04%	26.80%	27.20%	30.12%	53.30%
	RMSE	1.5494	1.5347	1.5924	1.6122	2.8495

Table 3.

The forecasting results of Experiment II.

4.3. Experiment III: 1-h average wind speed forecasting

In Experiment III, we will forecast 1-h average wind speed in three wind farms using ANFIS, BPNN, ELM, SVM, and ARMA, and 1-h average wind speed forecasting is a one-step forecast for three wind farms.³

It should be pointed that every day has 24 1-h average wind speed.

Table 4 shows the forecasting effectiveness of each model in this experiment.

Wind farms	Criteria	ANFIS	BPNN	ELM	SVM	ARMA
WFD1	MAE	1.5664	0.8586	0.8315	0.9069	1.2883
	MAPE	44.82%	20.15%	19.33%	24.92%	27.58%
	RMSE	2.0981	1.2376	1.1882	1.2805	1.8759
WFD2	MAE	1.7780	0.7924	0.7607	0.7919	2.5731
	MAPE	47.44%	20.30%	19.35%	21.92%	80.67%
	RMSE	2.3788	1.1838	1.1298	1.1712	3.3596
WFD3	MAE	0.8873	0.7964	0.7943	0.7914	1.1274
	MAPE	16.90%	15.07%	14.95%	14.91%	20.75%
	RMSE	1.1568	1.0568	1.0548	1.0539	1.4833

Table 4.

The forecasting results of Experiment III.

From Table 4, ELM and SVM have decent forecasting accuracy in WFD3, but ELM is better than SVM and other models in WFD1 and WFD2. ANFIS has a good performance in WFD3 and cannot forecast 1-h average wind speed with a reasonable accuracy. Similarly for ARMA, it works well in WFD1 and WFD3 but has a high MAPE in WFD2.

4.4. Experiment IV: 4-h average wind speed forecasting

It is similar to Experiment III, but this experiment shows the one-step-ahead wind speed forecasting effectiveness of each model. Table 5 demonstrates the results of five models obtained in three wind farms.

Wind farms	Criteria	ANFIS	BPNN	ELM	SVM	ARMA
WFD1	MAE	2.0516	1.6435	1.6209	1.5747	2.1725
	MAPE	55.75%	41.64%	38.76%	36.76%	48.81%
	RMSE	2.6411	2.1411	2.1478	2.1241	2.8977
WFD2	MAE	2.1538	1.4282	1.3489	1.3422	2.4222
	MAPE	55.53%	36.02%	34.27%	30.86%	74.71%
	RMSE	2.8318	1.9498	1.8597	1.8933	3.1035
WFD3	MAE	1.2564	1.2015	1.1483	1.1689	1.6348
	MAPE	23.37%	22.33%	21.13%	21.19%	29.41%
	RMSE	1.5951	1.5249	1.4856	1.5074	2.0953

Table 5.

The forecasting results of Experiment IV.

The analyzing results in Table 5 are quite different from that in Table 4. Specifically, SVM model has better performance in WFD1 and WFD2 than other models, and ELM is the best forecasting model in WFD3. For each model, wind speed in WFD3 can be forecasted more accurately than that in WFD1 and WFD2, and AI-based models (ANFIS, BPNN, ELM, and SVM) outperform ARMA in WFD3.

4.5. Experiment V: 1-day average wind speed forecasting

In this experiment, we will test the one-step-ahead forecasting effectiveness based on 1-day average wind speed database. Similar to other experiments, Table 6 demonstrates each model's forecasting effectiveness.

Table 6 shows that SVM has the best forecasting performance on three criteria in three wind farms. By contrast, ARMA and ANFIS have the higher MAPE, MAE, and RMSE than other models.

Wind farms	Criteria	ANFIS	BPNN	ELM	SVM	ARMA
WFD1	MAE	2.1001	1.9069	1.8377	1.7675	2.5753
	MAPE	52.78%	42.99%	45.70%	36.81%	60.59%
	RMSE	2.5729	2.3784	2.3053	2.2692	3.1655
WFD2	MAE	1.9384	1.9710	1.9038	1.7653	1.9271
	MAPE	51.81%	47.82%	48.19%	38.04%	52.61%
	RMSE	2.4685	2.4758	2.3848	2.3573	2.4679
WFD3	MAE	1.9288	1.8385	1.5886	1.5487	1.8763
	MAPE	32.30%	30.40%	26.39%	24.79%	30.47%
	RMSE	2.4229	2.3061	1.9152	1.9159	2.2568

Table 6.

The forecasting results of Experiment V.

5. Conclusion

The inherently intermittent nature and stochastic nonstationarity of wind sources bring great levels of uncertainty to system operators and are an urgent problem in the wind-forecasting field. This chapter introduces some basic forecasting approaches and corresponding procedures to forecast wind speed (detailed codes are attached). After simulating these methods (the testing period is 30 weeks) in MATLAB based on three databases, conclusions can be drawn as follows:

1. None of ANFIS, BPNN, ELM, SVM, and ARMA has the best forecasting performance in all experiments.

2. Based on the experimental results of different forecasting types (1-h-ahead forecasting, 4-h-ahead forecasting, 1-h-average-ahead forecasting, 4-h-average-ahead forecasting, and 1-day-average-ahead forecasting), the best model varies in time scales and time horizons.

3. The forecasting effectiveness differs a lot from database to database (in Experiment II, ANFIS is the best model of which MAPE is 26.04 in WFD3 but has a high MAPE, 54.41%, in WFD1).

4. AI-based models are more suitable for wind speed forecast than ARMA.

According to these conclusions, it is obvious that wind speed forecasting models or systems need to be built up under specific conditions in wind farms and model selection in wind speed forecasting plays a significant role to improve the accuracy.

Acknowledgments

The work was supported by the National Natural Science Foundation of China (grant no. 71573034).