An improved wind power forecasting model based on eliminating outliers

Introduction

For the past few years, the technology of wind power has been developed and is relatively mature, and the proportion of wind energy in the power grid shows a rising trend. However, the characteristics of wind such as intermittent, volatile and low energy density (Wang et al., 2011; Yao et al., 2011) always lead to intermittent and volatile power generation, which increases the difficulty of dispatching and the requirement of network preparation dosage (Wang et al., 2009). Therefore, developing wind-power research actively to improve the accuracy of forecasting has a realistic significance on dispatching and reducing the operation cost of the system (Yang et al., 2005).

There are lots of wind power forecasting models such as the time series analysis method (Liu et al., 2013; Qing and Yang, 2014), the neural network method (Sideratos and Hatziargyriou, 2007), the Kalman filtering method (Cassola and Burlando, 2012; Xiu and Guo, 2013), wavelet analysis (Grinstead et al., 2004; Tai et al., 2003), etc. In 2008, Pan Yifu developed a mixture forecasting model by integrating a time series with the Kalman filter, which could solve the problem of time delay and mine the regulation among dynamic data greatly. However, the defect is that the forecasting results are non-ideal when there are large leaps among the data.

What’s more, fewer models have considered the effects of abnormal observation data on the accuracy of a model, for example, it would lead to a large deviation from realistic results. In order to improve the accuracy of the forecast, it is necessary to eliminate the abnormal data first. In 1997, Wang Zhengming developed the one-by-one method of outlier rejection in a linear regression model. Nevertheless, the fact is that the observation data usually show the characteristic of a nonlinear relationship.

In this paper, we research a pointwise elimination method for the abnormal data of wind power based on the nonlinear model, and we obtain a criterion of stopping-rule. Furthermore, we propose the exponential smoothing method which considers the factor of trend to amend the large leaps among the data, and we develop an improved short-term forecasting model, which has a high accuracy. Finally, a numerical experiment shows the effectiveness of our model based on the data from some power supply companies.

This paper is divided into three sections. In the next section, we develop the pointwise elimination method for the abnormal data under the nonlinear model and show an improved stopping-rule. In the third section, we obtain an improved short-term forecasting model which integrates the exponential smoothing method with a time series and the Kalman filter. In the last section, a numerical experiment shows simulation results and illustrates the effectiveness of our method.

Outlier pointwise elimination method

Generally, there will inevitably exist some abnormal observation data among statistical data and it would have a negative impact on wind-power forecasting results. For the purpose of reducing the errors of forecast, it is very important to detect and eliminate the abnormal observation data.

Considering the nonlinear model

Y i = X i β + e i

(1)

where X i = ( 1, x i , x i 2 , ⋯ , x i n ) , X = X m × 1 = ( X 1 , X 2 , ⋯ , X m ) T , x i is the time at the point i , Y = ( Y 1 , Y 2 , ⋯ , Y m ) T , Y i is the observation data at the time x i , β = ( β 0 , β 1 , ⋯ , β n ) T , β is the n + 1 dimensional parameter vector, e = ( e 1 , e 2 , ⋯ , e m ) T , e i is the random error, e i ∼ N ( 0 , δ 2 ) and ( i = 1 , 2 , ⋯ , m ) .

Let

{ H = ( h i j ) m × m = X ( X T X ) − 1 X T σ = ( σ 1 , σ 2 , ⋯ , σ m ) T = ( I m − H ) Y ζ = ( ζ 1 , ζ 2 , ⋯ , ζ m ) T , ζ i = σ i 2 1 − h i h i = ∑ j = 1 m h i j 2 , 1 ⩽ i , j ⩽ m

(2)

Remark 2.1 Based on the data in M and N , we have the models

Y i = X i β + e i , e i ∼ N ( 0 , δ 2 ) Y j = X j β + η j , η i ∼ N ( 0 , k δ 2 )

Lemma 2.1 (Tong and Zhou, 2001) If h < 1 / 4 , k > 1 , we have

E ζ i < E ζ j , i ∈ M , j ∈ N

where M and N are two index sets and the total number of data in these two sets is m . If ζ j ≫ ζ i , h j ≪ 1 / 4 .

Lemma 2.1 shows that the greater the value of ζ i , the greater the probability of i ∈ N ; similarly, the greater the value of ζ j , the greater the probability of j ∈ M . We can identify the abnormal data by comparing the value of ζ i . The aim of the pointwise elimination method is to find out the maximum ζ i and eliminate the abnormal data one by one.

Lemma 2.2 (Wang and Wang, 1997) Under the assumption of model (1)

P { max 1 ⩽ i ⩽ m ζ i ⩾ w ∥ σ ∥ 2 } ⩾ α

(3)

For a given stopping-parameter w , from equation (3), the probability of normal data being misjudged as abnormal is less than α .

According to Lemma 2.2, the stopping-parameter w is the key to reserve normal data. In this paper, we obtain a method to determine the stopping-parameter and show an improved stopping-rule combining with the T test under the nonlinear model.

Theorem 2.1 For nonlinear model (1), let α be the confidence level and σ = ( I − H ) e , ζ i = σ i 2 / ( 1 − h i ) and 1 ⩽ i ⩽ m . If max 1 ⩽ i ⩽ m ζ i ⩾ w ∥ σ ∥ 2 and the stopping-parameter w = ( t 1 − α 2 ( m − 1 ) ) / ( m − 1 ) , then Y i is abnormal data.

Proof. According to polynomial fitting, we have

y ^ i = β ^ 0 + β ^ 1 x + β ^ 2 x 2 + … + β ^ n x n

where y ^ i is the fitted value of y i and β ^ = ( β ^ 0 , β ^ 1 , ⋯ , β ^ n ) is the fitted value of β .

Let

σ i = y i − y ^ i = ( I − H ) e i ∼ N ( 0 , δ 2 )

(4)

we have

y i − y ^ i δ ∼ N ( 0 , 1 )

(5)

From equations (4) and (5), there exists the test statistic T i

T i = σ i ∑ j ≠ i m σ j 2 / ( m − 1 ) ∼ t ( m − 1 )

(6)

where ∑ j ≠ i m σ j 2 = ∥ σ ∥ 2 and the critical of T i is t ( m − 1 )

From equation (2), rewriting ζ i = σ i 2 / ( 1 − h i ) as

σ i 2 = ( 1 − h i ) ζ i , 1 ⩽ i ⩽ m ,

(7)

Combined with equations (6) and (7) and Lemma 2.1, we obtain

T i = ζ i ∥ σ ∥ 2 m − 1 ∼ t ( m − 1 )

(8)

By equation (8) and Lemma 2.2, it follows that Y i is the abnormal data when

T i = ζ i ∥ σ ∥ 2 m − 1 ⩾ t 1 − α ( m − 1 )

(9)

Thus, it is obvious that equation (9) satisfies

ζ i ⩾ t 1 − α 2 ( m − 1 ) m − 1 ∥ σ ∥ 2

then we have

w ∥ σ ∥ 2 = t 1 − α 2 ( m − 1 ) m − 1 ∥ σ ∥ 2

which implies that

w = t 1 − α 2 ( m − 1 ) m − 1

This completes the proof of Theorem 2.1.

Improved short-term forecasting model

In this section, we research a low order Auto-Regressive and Moving Average Model (ARMA) which is a time series model and can reflect the signal variation. Under the prediction equation from the ARMA model, the Kalman filtering state equation and measurement equation can be deduced. Then, we take advantage of the Kalman filtering prediction iterative equation to forecast. In order to reduce the errors, we utilize the exponential smoothing method to smooth the large leaps among the data before forecasting.

The ARMA time series forecasting

Let { x t } be a zero mean stationary series, which satisfies the following model

x t − φ 1 x t − 1 − … − φ p x t − p = ε t − θ 1 ε t − 1 − … − θ q ε t − q

(10)

where ε t is the stationary white noise. The mean and variance of ε t are zero and σ ε 2 , ( t = 0 , ± 1, ⋯ ) . Model (10) can be denoted by ARMA( p , q ) model, p and q are the orders of the model.

Applying operators φ ( L ) , θ ( L ) the model can be simplified as

φ ( L ) y t = θ ( L ) ε t

(11)

where the solutions of φ ( L ) and θ ( L ) must be out of the unit circle. Therefore, it is necessary to test the stationarity of raw data at the beginning of forecasting. For the unsteady data, we can utilize the difference method to get the steady series.

Considering the Augmentation Dickey–Fuller(ADF) testing model

Δ x t = α + β t + δ x t − 1 + ∑ t = 1 p β i Δ x t − i + ε t

(12)

where ∆x_t is the first-order difference series, α is the constant term, β t is the time trend term, δ x t − 1 is the first-order lag item and ∑ t = 1 p β i Δ x t − i is the difference lag item. In practice, the value of p has to ensure the minimum ε t .

After testing, we estimate the parameters of the prediction equation of the ARMA model through the maximum likelihood estimation.

The fundamental of Kalman filtering

The Kalman filtering method describes the filter through the state space model of the linear stochastic system which is composed of the state equation and observation equation. Besides, it utilizes the recurrence properties of the state equation and the linear unbiased minimum variance estimation criterion to estimate the state variable of filter optimally.

The general linear discrete system is denoted as

X ( k + 1 ) = Φ ( k + 1, k ) X ( k ) + Γ ( k + 1, k ) ω ( k )

(13)

Z ( k + 1 ) = H ( k + 1 ) X ( k + 1 ) + ν ( k + 1 )

(14)

where X ( k ) is the n-dimensional state vector and X ( k ) contains the matrix P ( k ) which means the forecasting error covariance matrix for the moment k (Pan et al., 2008), Z ( k ) is the m-dimensional observation vector, ω ( k ) is the p-dimensional system noise vector, ν ( k ) is the m-dimensional measurement noise vector, Φ ( k + 1, k ) is the excitation transfer matrix from the moment k to the moment k + 1 , Γ ( k + 1, k ) is the incentive transfer matrix from the moment k to the moment k + 1 and H ( k + 1 ) is the predicted output transfer matrix at the moment k + 1 . Equation (13) is the state equation, and equation (14) is the measurement equation.

Mixture forecasting model based on ARMA with Kalman filter

Through converting the model (10) into state space, we would obtain the state equation and the measurement equation of the Kalman filter.

First, we rewrite equation (10) as

x t = φ 1 x t − 1 + … + φ p x t − p + ε t − θ 1 ε t − 1 − … − θ q ε t − q − 1

(15)

where { x t } is the zero mean time series.

Then, we convert equation (15) into the state space model

{ X t = A X t − 1 + B E t Y t = C X t

(16)

where

X t = ( x t , x t − 1 , ⋯ , x t − m − 1 )

E t = ( ε t , ε t − 1 , ⋯ , ε t − m − 1 )

C = ( 1 , 0 , ⋯ , 0 )

A m × m = ( a * a m I m − 1 0 )

a * = ( a 1 , a 2 , ⋯ , a m − 1 )

B = ( b * 0 )

b * = ( 1 , b 1 , ⋯ , b m − 1 )

and I m − 1 is the autoregressive coefficient of the state space model.

Improvement of the model

In general, the forecasting results of the mixture model have been satisfactory. However, facing the large leaps among data, the results would be non-ideal. Hence we take advantage of the three times exponential smoothing method, which considers the factor of trend to amend the large leaps among the data, and integrate it with the mixture model to develop an improved short-term forecasting model.

Considering the three times exponential smoothing moving model

Y ^ t + m = a t + b t m + c t m 2 , m = 1 , 2 , ⋯

(17)

{ a t = 3 S t ( 1 ) − 3 S t ( 2 ) + S t ( 3 ) b t = ρ 2 ( 1 − ρ ) 2 [ ( 6 − 5 ρ ) S t ( 1 ) − 2 ( 5 − 4 ρ ) S t ( 2 ) + ( 4 − 3 ρ ) S t ( 3 ) ] c t = ρ 2 2 ( 1 − ρ ) 2 [ S t ( 1 ) − 2 S t ( 2 ) + S t ( 3 ) ]

where S t ( 1 ) , S t ( 2 ) , S t ( 3 ) are the first, second and third exponential smoothing values, Y ^ t + m is the forecasting value at the moment t + m , a t , b t , c t are the parameters and ρ is the weight parameter which means the proportion of new data of the forecast one.

In fact, the rising or falling trend of some parts of the data will lead to the forecast behind the actual trend. Thus, we introduce the trend factor into the model equation (17).

The adjusted model is

Y ′ T = Y ^ t + m + T t + m

T t + m = ( 1 − ρ ) T t + m − 1 + ρ ( Y ^ t + m − Y ^ t + m − 1 )

where Y ′ T is the forecasting value include trend at the moment t + m , T t + m is the time trend after smoothing at the moment t + m and T t + m − 1 is the time trend after smoothing at the moment t + m − 1 .

A numerical example

In this experiment, we take the wind-power data per five minutes (216 in all) from some wind farms in Zhejiang province on January 1, 2009 as the research object.

First of all, we conduct the three order polynomial fitting on the data per hour, and adopt Theorem 2.1 to eliminate the abnormal data. Figure 1 shows the fitting curve before and after elimination.

OPEN IN VIEWER

Figure 1.

The fitting curves before and after elimination.

In accordance with the new fitting curve, we update the raw data (refer to Table 1). The trend curve of wind power (Figure 2) shows clearly that there exist some large leaps among the data. Thus, we apply the three times exponential smoothing method to optimize the data as a series { x t } and let ρ = 0 . 5 according to the fluctuation of wind power.

OPEN IN VIEWER

Table 1.

Updated data of wind power.

Moment	Wind power
0.00–0.55	56	56	58	57	54	50	50	51	51	54	55	56
0.60–1.15	60	59	58	63	70	69	71	75	73	75	78	75
0.20–1.75	77	90	96	93	92	92	90	83	83	83	83	82
0.80–2.35	79	78	77	77	75	77	80	86	89	89	89	95
0.40–2.95	94	93	93	93	94	95	96	97	97	97	95	92
0.00–3.55	94	95	92	89	85	83	78	77	73	71	69	68
0.60–4.15	67	66	67	68	70	71	74	77	79	82	83	83
0.20–4.75	83	88	87	88	87	84	82	81	77	76	80	80
0.80–5.35	77	80	81	84	85	85	86	84	83	83	83	80
0.40–5.95	81	80	76	74	76	75	73	73	72	73	76	79
0.00–6.55	78	76	78	78	80	82	84	85	89	88	88	87
0.60–7.15	90	91	90	92	93	93	95	91	93	94	94	96
0.20–7.75	96	93	90	96	97	97	93	100	95	98	97	99
0.80–8.35	97	94	97	99	98	97	99	93	95	90	93	94
0.40–8.95	97	89	91	93	86	87	76	75	78	79	76	73
0.00–9.55	78	80	76	74	70	67	62	61	60	62	62	61
0.60–10.15	62	60	64	75	74	66	58	53	48	46	45	48
0.20–10.75	54	60	63	62	63	62	58	56	52	51	53	52

OPEN IN VIEWER

Figure 2.

The trend curve of wind power.

Then, combining the first 20 autocorrelation coefficients of the sequence { x t } (Figure 3) and the ADF test method, we can estimate that the series { x t } is unsteady. Utilizing the difference method on the series { x t } , we have a stationary series { y t } . Figure 4 shows the first 20 autocorrelation coefficients of the series { y t } .

OPEN IN VIEWER

Figure 3.

The first 20 autocorrelation coefficients of { x t } .

OPEN IN VIEWER

Figure 4.

The first 20 autocorrelation coefficients of { y t } .

Furthermore, the orders of model (10) could be identified as ARMA(3,2) by the Akaike information criterion (Akaike, 1974) and the parameters would be estimated by maximum likelihood method as follows

y t = − 0 . 1972 + 0 . 8682 y t − 1 − 1 . 0675 y t − 2 + 0 . 2387 y t − 3 + ε t + 0 . 7060 ε t − 1 + 0 . 9872 ε t − 2

Let { y t } be the zero mean series { z t } satisfying

z t = 0 . 8682 y t − 1 − 1 . 0675 y t − 2 + 0 . 2387 y t − 3 + ε t + 0 . 7060 ε t − 1 + 0 . 9872 ε t − 2

(18)

Converting equation (18) into state space, we can get the state equation (19) and the measurement equation (20)

( z t z t − 1 z t − 2 ) = ( 0 . 8682 − 1 . 0675 0 . 2387 1 0 0 0 1 0 ) ( z t − 1 z t − 2 z t − 3 ) + ( 1 0 . 7060 − 0 . 9872 0 0 0 0 0 0 ) ( ε t ε t − 1 ε t − 2 )

(19)

z t = ( 1 , 0 , 0 ) ( z t , z t − 1 , z t − 2 ) T + ν ( t )

(20)

where ν ( t ) is the measurement of additional noise which also named white noise. Ordinarily, the initial values of z ( 0 ) , ν ( 0 ) are z ( 0 ) = [ 0 ] , p ( 0 ) = 10 I .

Finally, we realize the improved short-term forecasting model through MATLAB to forecast the data. Figure 5 shows the curves of the comparison between actual data and forecasting data for the next hour.

OPEN IN VIEWER

Figure 5.

The wind-power forecasting result under the hybrid model.

The experiment shows the following.

1. The mixture model which integrates the time series with the Kalman filter can catch the dynamic regulation of wind-power data effectively. Based on the updated data, we can forecast the data for the next hour and get an error precision of 1.69%.

2. The mixture model does not fit the large variational data such as big leap points. It is effective to utilize the three times exponential smoothing method to improve the accuracy of the forecasting results. Under the improved short-term forecasting model, the error precision of the results can decrease to 0.94% 0 . 94 % , which shows a high significance of practical application.

3. Although we can obtain high-precision results through the improved short-term forecasting model, the length of forecast is short. With the increase in the length of forecasting, the accuracy of results will decrease.