Stochastic predictive control of battery energy storage for wind farm dispatching: Using probabilistic wind power forecasts

Renewable Energy 80 (2015) 286e300

Contents lists available at ScienceDirect

Renewable Energy journal homepage: www.elsevier.com/locate/renene

Stochastic predictive control of battery energy storage for wind farm dispatching: Using probabilistic wind power forecasts Peng Kou a, *, Feng Gao b, Xiaohong Guan c a

State Key Laboratory of Electrical Insulation and Power Equipment, School of Electrical Engineering, Xi'an Jiaotong University, Xi'an 710049, China Systems Engineering Institute, SKLMS, Xi'an Jiaotong University, Xi'an 710049, China c Systems Engineering Institute, MOE KLINNS, Xi'an Jiaotong University, Xi'an 710049, China b

a r t i c l e i n f o

a b s t r a c t

Article history: Received 20 December 2013 Accepted 3 February 2015 Available online

The limited dispatchability of wind energy poses a challenge to its increased penetration. One technically feasible solution to this challenge is to integrate a battery energy storage system (BESS) with a wind farm. This highlights the importance of a BESS control strategy. In view of this, a stochastic model predictive control scheme is proposed in this paper. Based on the forecasted wind power distributions, the proposed scheme ensures the optimal operation of BESS in the presence of practical system constraints, thus bringing the wind-battery combined power output to the desired dispatch levels. The salient feature of the proposed scheme is that it takes into account the non-Gaussian wind power uncertainties. In this scheme, a probabilistic wind power forecasting model is employed as the prediction model, which quantifies the non-Gaussian uncertainties in wind power forecasts. Using chance constraints, the quantified uncertainties are incorporated into the controller design, thus forming a chance constrained stochastic programming problem. Using warping function, this problem is recast as a convex quadratic optimization problem, which is tractable both theoretically and practically. This way, the proposed control scheme handles the non-Gaussian uncertainties in wind power forecasts. The simulation results on actual data demonstrate the effectiveness of the proposed scheme. The data used in the simulation are obtained in the real operation of a wind farm in China. © 2015 Elsevier Ltd. All rights reserved.

Keywords: Wind energy Power dispatch Stochastic model predictive control Battery energy storage Probabilistic wind power forecast Stochastic optimization

1. Introduction In the new century, the electric power industry is facing unprecedented challenges arising from the continuous depletion of fossil fuels and the worsening global environment. These challenges have spurred a great interest in integrating renewable energies into existing power systems. Among a wide variety of the renewable energy technologies, wind power is considered as the most promising one, because of its mature technology and low operation costs. The growth of wind power has surpassed the most optimistic expectations, the worldwide installed capacity of wind generation reached 273 GW in December 2012 [1]. Nevertheless, although wind power has significant environmental and economic advantages, the stochastic nature of wind impedes the large scale integration of wind power into the grid. Due to the erratic nature of the atmosphere, the wind power output is always

* Corresponding author. Tel.: þ86 29 82667679; fax: þ86 29 82667856. E-mail address: [email protected] (P. Kou). http://dx.doi.org/10.1016/j.renene.2015.02.001 0960-1481/© 2015 Elsevier Ltd. All rights reserved.

fluctuating. In Fig. 1(a), the square marked line shows the typical power output profile of a large wind farm. It is seen that the generated power has steep rises and sudden drops in a short time. Interconnecting such a highly fluctuating power sources to the grid would pose a series of problems, such as stability, protection, flow control, power quality, and power dispatching [2,3]. Therefore, it is highly desired to reduce the fluctuation of wind power output and improve its dispatchability. Recent advances in battery and power electronic technologies provide a technically feasible solution to this issue, that is, to integrate a battery energy storage system (BESS) with a wind farm [4e6]. This way, the power fluctuation of a wind farm can be compensated by the charging and discharging of BESS. Based on this, we can drive the wind-battery combined power output to a desired dispatch curve. In this sense, BESS makes a wind farm more dispatchable like other conventional generation plants. This approach has been implemented in practical systems. For instance, in 2011, the State Grid of China put the Zhangbei national energy storage and transmission demonstration project into operation. This project combines 100 MW of wind energy generation, 36 MWh of BESS, and a smart power

P. Kou et al. / Renewable Energy 80 (2015) 286e300

287

Fig. 1. Illustration of (a) control performance degradation caused by wind power forecast uncertainties, (b) quantification of wind power forecast uncertainties using predictive distributions.

transmission system [7,8]. It serves as a demonstration of a stable solution for transferring large amounts of dispatchable, predictable, and dependable renewable electricity to the grid on an unprecedented scale. The integration of wind farm and BESS requires an optimal control strategy, which determines the optimal charge and discharge commands for the BESS. Otherwise, the BESS will lose its function as shock absorbers once it is fully charged or discharged [9]. Due to the presence of the BESS power and energy capacity constraints, this control problem is a constrained one. Furthermore, because wind power is perpetually varying, the power fluctuation of a wind farm has to be mitigated continuously. These issues give rise to an online constrained control problem. Model predictive control (MPC) is an ideal framework to tackle this problem, since it explicitly takes into account the system constraints and repeatedly optimizes performance over a receding horizon. For this reason, there is a significant body of literature dealing with this subject. Teleke et al. [10] proposed an open-loop MPC scheme for the windbattery system. This scheme based on the 1-h ahead wind power forecasts and a linear battery model, it makes the power injected to the grid as smooth as possible, and so the wind farm can be dispatched on an hourly basis. Khalid and Savkin [11] developed an autoregressive model to forecast the wind power output, and proposed an MPC scheme that makes the output of the wind farm as close as possible to the forecasted values. A later paper by the same authors [12] used a similar scheme to solve the frequency control problem of a wind farm. Khatamianfar et al. [13] further developed this scheme. In the electricity market environment, they designed a fuzzy controller to generate the economic dispatch curve, which is a function of the electricity price. Then the MPC scheme in Ref. [11] is adopted to track this curve, i.e., to sell more electricity at peak price periods and store it at off-peak periods. Qi et al. [14] proposed a supervisory MPC system for the optimal operation of a hybrid windesolareBESS system. In this system, MPC determines the optimal reference power signals for the wind and solar subsystems, and then the local controllers drive the wind and solar subsystems to the desired reference signals. In Ref. [15], the same author extended this scheme to a distributed setting. In order to coordinate the actions of wind and solar subsystems accordingly, they developed two supervisory MPC architectures, one based on sequential distributed control and the other based on iterative distributed control, respectively. Venayagamoorthy et al. [9] utilized a neural network to forecast the active and reactive power of a wind farm. Using the forecasted values, an MPC scheme based on adaptive critic design was proposed, with the goal of smoothing wind power fluctuations.

All the aforementioned works are based on the deterministic MPC (DMPC) framework. In this framework, a prediction model is exploited to forecast the expected value of wind power output, i.e., to provide the point forecasts of wind power output. Based on these point forecasts, at each time stamp a predictive control performance index is optimized under operating constraints w.r.t. a sequence of future control actions. In this sense, the control performances of the preceding schemes strongly depend on the wind power point forecast accuracies. However, due to the variable and stochastic nature of the wind, the accuracy of such forecasts cannot be guaranteed and fairly low on average [16]. That is to say, the performances of the aforementioned schemes depend on the forecasting results that have inherent uncertainties. Fig. 1(a) gives a graphic illustration of this fact. In this figure, one can see that the actual wind power measurements at time 7.5 h and 40.5 h are quite far from the point forecasts. In this case, the control actions computed by DMPC may result in the failure of tracking the dispatched target curve and the violation of operating constraints. To address this, a stochastic MPC (SMPC) scheme is needed [17e19]. Rather than merely point forecasts, the prediction model in SMPC can give probabilistic forecasts, i.e., the wind power predictive distributions in Fig. 1(b) [16,20,21]. In this figure, one can see that, the actual wind power measurements at time 7.5 h and 40.5 h all fall within the regions predictive from the predictive distributions. This way, the wind power forecast uncertainties are quantified by these predictive distributions. At each time stamp, these quantified uncertainties are explicitly considered in the SMPC problem formulation, thus forming a stochastic optimization problem. Then the optimal control actions are computed by solving this problem. Since the resulting control law takes into account the wind power forecast uncertainties, SMPC can achieve more robust control performance than DMPC. From Fig. 1(b), it is seen that there are two key factors for applying SMPC to wind-battery system, i.e., the quantification and handling of wind power forecast uncertainties. These uncertainties are non-Gaussian. The reason lies in the fact that, the wind power is mainly a function of wind speed, and the turbine power curve leads a nonlinear transformation from the uncertainties in wind speed to the uncertainties in wind power. Due to the sigmoid shape of the turbine power curve, the transformed uncertainties (i.e., the uncertainties in wind power) are no longer Gaussian [22]. Such nonGaussian uncertainty imposes two requirements on the SMPC scheme: (i) its prediction model should have the ability to provide the non-Gaussian predictive distributions; (ii) the design of the controller should consider the non-Gaussian disturbances. These two requirements pose a challenge to the existing SMPC schemes,

288

P. Kou et al. / Renewable Energy 80 (2015) 286e300

since few of them considers the non-Gaussian predictive distributions and non-Gaussian disturbances. In these schemes, a common assumption when facing uncertainties is to model them as Gaussian predictive distributions and Gaussian disturbances [17,23e26]. Warped Gaussian process (WGP) offers a possible solution to these two requirements [27]. In our previous work [28], we established a sparse online warped Gaussian process (SOWGP) model, which can provide reliable non-Gaussian predictive distributions of wind power output, thus accurately quantifying the nonGaussian uncertainties in wind power. Furthermore, SOWGP introduced a warping function, with which the non-Gaussian disturbances can be converted to Gaussian ones in the latent space, thus easy to handle. For these reasons, the integration of SOWGP and SMPC has the potential to overcome the above challenge. Therefore, it is the objective of this work to study the integration of SOWGP and SMPC, and to apply the integrated control scheme to the wind-battery system. Consider above issues, this paper presents a new SMPC scheme for the wind-battery system, with the goal of improving the wind power dispatchability and reducing its fluctuation. With the proposed scheme, a wind farm can be dispatched based on the forecasted wind conditions. The main advantage of this scheme lies in its ability to handle the non-Gaussian wind power uncertainties. In the proposed scheme, the wind power forecast uncertainties are modeled as non-Gaussian disturbances. Employing SOWGP as the prediction model, the statistics of these disturbances can be quantified by the output of the prediction model, i.e., its nonGaussian predictive distributions. Using chance constraints, the quantified uncertainty is explicitly incorporated into the controller design. This way, we obtain a set of chance constraints with nonGaussian coefficients. Using the warping function in SOWGP, these constraints can be further transformed to deterministic convex constraints, and thus can be readily incorporated into the formulation of the MPC optimization problem. The resulting online optimization problem can be solved using standard quadratic programming routines. By these means, the proposed SMPC scheme handles the non-Gaussian uncertainties in wind power forecasts. The remainder of this paper is organized as follows. Section 2 presents the system model and the precise problem formulation. In Section 3, we give an overview of SOWGP, which plays the role of the prediction model. The detailed description of the proposed SMPC scheme is given in Section 4. Section 5 investigates the advantage of the proposed scheme. Section 6 contains the description of the simulations and the discussion of the results. The data used in the simulations are obtained from a wind farm in Inner Mongolia, China. Finally, we conclude the paper in Section 7. 2. Wind-battery system description and control problem statement 2.1. Description of the wind-battery system The overall structure of the wind-battery system considered in this work is depicted in Fig. 2. This system is based on the aforementioned Zhangbei national energy storage and transmission demonstration project [7,8]. It consists of a wind farm, a battery bank, an AC/DC converter, two transformers, an SMPC controller, and a probabilistic wind power forecasting model. The system is connected to the main grid through transmission lines. The battery bank (i.e., the BESS) is connected to the system at the point of common coupling and is charged/discharged through an AC/DC converter to smooth the net power injected to the grid. As mentioned before, our goal is to make a wind farm more dispatchable like the other conventional power plants. In other

Fig. 2. Schematic diagram of a grid-connected wind-battery system, together with the proposed SMPC controller, and the probabilistic wind power forecasting model.

words, we want to bring the power output of a wind farm as close as possible to a desired dispatch curve PD. To achieve this goal, the operating mechanism of the system in Fig. 2 is as follows. At each time stamp k, the probabilistic wind power forecasting model makes an H-steps ahead forecast, and thus providing the predictive distributions p(PW(k þ 1)), p(PW(k þ 2)), …, p(PW(k þ H  1)). Based on these predictive distributions, the SMPC controller computes the appropriate amount of BESS charging/discharging power PB(k), which compensates the differences between the dispatch level PD(k þ 1) and the actual wind power output PW(k þ 1). Meanwhile, the BESS power and energy capacity constraints are considered and satisfied. At the next time stamp k þ 1, this procedure is repeated. This way, the net power injected to the grid (i.e., the wind-battery combined power output PG(k)) will follow the desired dispatch curve. For a wind farm, the average of its forecasted power output for the next few hours is a good choice for its desired power dispatch level PD [10,29]. Therefore, in this work, the next 4 h average of the wind power point forecasts is selected as PD. Correspondingly, the dispatch interval is set to be 4 h. Notice that, devising optimal dispatch curves for wind farms is also a very important topic of active research, but it is not the focus of this paper. 2.2. System modeling Based on the structure shown in Fig. 2, a linear dynamic model for the wind-battery system can be formulated as follows

PG ðk þ 1Þ ¼ PB ðkÞ þ PW ðkÞ; EB ðk þ 1Þ ¼ EB ðkÞ  hDTB PB ðkÞ;

(1)

where PG(k), net power injected to the grid, i.e., the combined power output; PB(k), power control signal. When PB(k) > 0, it represents the amount of the BESS discharging power. Conversely, when PB(k) < 0, it represents the amount of the BESS charging power; PW(k), original wind power output; EB(k), state of charge (SOC) of the BESS. The SOC of a BESS is its available capacity expressed as a percentage of its rated capacity; h, coefficient of BESS charge/ discharge efficiency; DTB, MW-to-MWh conversion factor, e.g., when the sampling time is 5 min, then DTB ¼ 5 min/60 min ¼ 1/12. To facilitate the subsequent analysis and design, all of PG, PB, PW, and EB are normalized with respect to the nominal capacity of the wind farm, i.e., PW ¼ 100% means nominal output power, PW ¼ 0 means zero output power, etc. This model closely follows the models presented in Refs. [11,13]. The primary difference between them lies in the treatment of the wind farm power output. In Refs. [11,13], the wind farm power

P. Kou et al. / Renewable Energy 80 (2015) 286e300

output PW is considered as a deterministic variable, hence the uncertainties in wind power are neglected. However, as discussed in Section 1, such uncertainties may significantly affect the control performance of the wind-battery system. Therefore, unlike in Refs. [11,13], we treat PW as a stochastic variable in this paper. Subsequently, the state space form of this linear dynamic model can be written as follows

xðk þ 1Þ ¼ AxðkÞ þ b1 uðkÞ þ b2 wðkÞ; yðkÞ ¼ cxðkÞ;




x1 ðkÞ



x2 ðkÞ control input uðkÞ ¼ PB ðkÞ;


¼

PG ðkÞ EB ðkÞ

2.3.1. SOC constraint Overcharge or overdischarge would adversely impact the lifetime of the battery units, and make the BESS lose its function as a shock absorber. Therefore, the SOC of the BESS should be kept within a suitable range. In this work, this range is empirically set as

 ; (3)

additive disturbance wðkÞ ¼ PW ðkÞ; system output yðkÞ ¼ PG ðkÞ;




0 0



  0 1 1 ; b2 ¼ ; b1 ¼ ; c ¼ ½1 1 hDTB 0

0 :

(4)

In model (2)e(4), the original power output of the wind farm is modeled as an additive stochastic disturbance w(k). The distribution of this disturbance represents the uncertainties associated with the wind power. Once the state space model is obtained, we need to check its controllability. Since the controllability matrix Q ¼ [b1, Ab1] is full rank, this system is controllable. 2.3. MPC formulation MPC is a powerful control technique that has been successfully applied in many fields [30]. In MPC, at each time stamp, a finite horizon optimal control problem is formulated and solved over a prediction horizon, subject to some state and input constraints. The result is a sequence of control actions, which satisfies the current and future constraints of the plant while optimizing some given performance index. Only the first control action from this sequence is implemented. At the next time stamp, the procedure is repeated using the new measurements of the states. In terms of the windbattery system control, this means that at every time stamp, a BESS charging/discharging plan is formulated for the next several time stamps, based on the forecasts of the future wind power outputs. The first step of this plan is applied to the BESS. At the next time stamp, this procedure is repeated using the new actual wind power measurement. In the MPC framework, a performance index has to be specified to find the optimal control actions. Since our goal is to improve the wind power dispatchability and reduce its fluctuation, it is desired that the net power injected to the grid follows the dispatch curve. Furthermore, in order to extend the lifetime of the BESS, it is also desired that the amount of charging/discharging power is minimized. In this case, the performance index can be represented by a quadratic cost function as follows

J¼b

H X h¼1

ðPG ðk þ hjkÞ  PD ðk þ hjkÞÞ2 þ ð1  bÞ

socmin  EB ðk þ hjkÞ  socmax ;

h ¼ 1; 2; …; H;

(6)

where socmin ¼ batcap  20% and socmax ¼ batcap  80%, respectively. batcap is the rated capacity of the BESS, which is normalized with respect to the nominal capacity of the wind farm.

together with the system matrices

A¼

trade-off parameter between the above two objectives. For simplicity, both the prediction horizon and the control horizon are set to H. This cost function penalizes the deviations from the dispatch power curve while minimizing the control effort. The minimization of the cost function (5) is subject to a set of system constraints. For a wind-battery system, the major constraints having a great influence on its behavior are the BESS energy and power limits. These limits lead to the following constraint.

(2)

with

state vector xðkÞ ¼

289

H1 X

uðk þ hjkÞ2 ;

2.3.2. Charge and discharge power constraint To protect the BESS from being damaged in charging and discharging process, a limitation needs to be put on its maximum change and discharge power, that is

PB;max  PB ðk þ hjkÞ  PB;max ;

h ¼ 1; 2; …; H;

(7)

where PB,max is the maximum charging/discharging power of the BESS. It is mainly decided by the DC voltage level and the maximum charge/discharge current of the BESS, which is usually taken as batcap/2, e.g., a current of 10 A with a capacity of 20 Ah. Hence, PB,max can be calculated as PB,max ¼ VBbatcap/2, where VB is the DC voltage level of the BESS [31]. Combining (5)e(7), together with the notation in (2), the MPC problem at time step k can be expressed as follows:

min J ¼ b

uðkþhjkÞ

H  X

2 yðk þ hjkÞ  yref ðk þ hjkÞ

h¼1

þ ð1  bÞ

H1 X

(8a) 2

uðk þ hjkÞ ; subject to

h¼0

socmin  x2 ðk þ hjkÞ  socmax ;

h ¼ 1; 2; …; H;

(8b)

PB;max  uðk þ hjkÞ  PB;max ;

h ¼ 1; 2; …; H;

(8c)

where yref is the reference trajectory, i.e., the dispatched target curve. In problem (8), both the state x and the output y are depending on the disturbance w, as indicated in (2). Therefore, to solve the problem (8) there are three key issues: (i) quantify the uncertainties associated with w; (ii) explicitly incorporate the quantified uncertainties into (8); and (iii) solve the resulting stochastic optimization problem. We will discuss these three issues separately in the two subsequent sections. 3. Probabilistic wind power forecasting using SOWGP

h¼0

(5) where (k þ hjk) denotes the value predicted for time stamp k þ h based on the information available at time stamp k, b 2 (0,1) is a

In the proposed SMPC scheme, the uncertainties in wind power forecasts (i.e., the uncertainties associated with w) are quantified by the prediction model. As discussed in Section 1, such uncertainties cannot be quantified by Gaussian, so the prediction model should

290

P. Kou et al. / Renewable Energy 80 (2015) 286e300

have the ability to provide the non-Gaussian predictive distributions. To address this, in the proposed SMPC scheme, an SOWGP model is utilized as the prediction model of the plant. SOWGP is an online probabilistic wind power forecasting model. A key feature of this model is that it can provide the non-Gaussian predictive distributions, with which the non-Gaussian uncertainties in wind power forecasts can be quantified [28]. In this section, we give a brief description of this model. From the perspective of machine learning, wind power forecast is a regression problem. The input explanatory variables form the D-dimension input vector v, while the actual wind power generation is the corresponding real valued target w. Since we use SOWGP as a prediction model in MPC, we focus on the short-term forecasting. For very short forecast horizons that from a few minutes up to 3 h, it has been stated that the dominant explanatory variables for the wind power forecasting are historical wind speed observations rather than meteorological data [32e34]. For forecast horizons greater than or equal to 3 h, numerical weather prediction (NWP) data should be used as features [32]. Therefore, in SOWGP, the historical local wind speed observations, the historical wind speed observations from multiple nearby sites, as well as the NWP data, are used as the explanatory variables. SOWGP is an extension of the standard Gaussian process (GP). Given a dataset D ¼ ðV; wÞ consisting of N input vectors N V ¼ fvn gN n¼1 and corresponding targets w ¼ fwn gn¼1 . In standard GP, we assume that the relation between the input vector and the target is given by

wn ¼ f ðvn Þ þ 3 n ;

(9)

where f is the underlying function, 3 n is Gaussian i.i.d. noise that follows N ð0; s2w Þ. Furthermore, we assume that f ¼ [f(v1), f(v2), …, f(vN)]T behaves according to a Gaussian process, that is pðfjv1 ; v2 ; …; vN ¼ N ð0; KÞÞ. Here K is the kernel matrix with elements Kij ¼ k(vi, vj), k being the kernel function parameterized by the hyperparameters q. Then the distribution of w conditioned on the values of f is given by an isotropic Gaussian

pðwjf; VÞ ¼ N

  f; s2w I :

(10)

From the property of the Gaussian distribution, we can see that the marginal distribution of w is

Z pðwjV; qÞ ¼

pðwjf; VÞpðfjVÞdf ¼ N

  0; K þ s2w I :

(11)

For a testing point v*, the predictive distribution given by standard GP is also Gaussian, that is

 pðw* jv* ; q; D Þ ¼ N

 1 k* T K þ s2w I w; s2w

 1  þ k**  k* T K þ s2w I k* :

However, in some applications such as the wind power forecasts, the targets w and the noise 3 cannot be assumed to be Gaussian. In this case, standard GP becomes unsuitable. WGP handles this problem by transforming the targets w from the original observation space to a latent space, thus forming a latent targets z ¼ fzn gN n¼1 . This transformation is implemented by a warping function g(w), that is

zn ¼ gðwn Þ ¼ wn þ a tanhðbðwn þ cÞÞ;

(12)

where J ¼ {a, b, c} are the warping parameters to be estimated, with a, b  0. In the latent space, z can be well-modeled by a standard GP, whose log likelihood is given by

1  N 1  1   ln pðzjV; qÞ ¼  zT K þ s2z I z  ln K þ s2z I   ln 2p: 2 2 2 (13) Substituting (12) into (13), the log likelihood becomes

 1  1 1   ln pðwjV; q; JÞ ¼  gðwÞT K þ s2z I gðwÞ  ln K þ s2z I  2 2 N X vgðwn Þ N þ log  ln 2p: vwn 2 n¼1 (14) Training of WGP is simply achieved by maximizing (14) with respect to q and J. This way, both the hyperparameters q and the warping parameters J can be learned simultaneously. When predicting for a testing point w*, in latent space, the predictive distribution is a standard Gaussian distribution, which is given by

pðz* jv* ; D Þ ¼ N

    1 b z*; b s 2z* ¼ N k* T K þ s2 I u; s2z þ k**  1   k* T K þ s2z I k* ;

(15)

where k* ¼ [k(x1, x*), …, k(xN, x*)]T and k** ¼ k(x*, x*). Untransforming (15) through the warping function, we can get the predictive distribution for the original observation space, that is

g 0 ðw* Þ pðw* jv* ; D Þ ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffi exp 2pb s 2z*



 2 ! 1 gðw* Þ  b z* : b z* 2 s

(16)

The probability density function (pdf) in (16) is the direct output of the WGP. In general, this predictive pdf is asymmetric, and therefore quite far from Gaussian. The WGP model described above handles the non-Gaussian process and non-Gaussian noise, which are two important features of the wind power series. Therefore, it is suitable for the probabilistic wind power forecasting. The application of WGP to probabilistic wind power forecasting is diagrammatically illustrated in Fig. 3. Roughly, it involves three stages. Firstly, WGP transforms the non-Gaussian target value (i.e., the wind power series w) to a latent series z, which is well-described by a Gaussian process. Secondly, in the latent space, WGP models the relation between latent target series z and explanatory variable v using a latent GP model, and makes predictions with this latent GP. Finally, the forecasting results are untransformed back to the original observation space by the inverse warping function. In this way, we get the non-Gaussian predictive distributions for the wind power series. Wind generation is a process whose characteristics change over time, so both the underlying function f in (9) and the warping function g in (12) cannot be assumed constant over time. This timevarying characteristic makes any offline prediction model to be a valid approximation only in a certain span of time. To address this, we introduce a sparse online algorithm to WGP. This algorithm

P. Kou et al. / Renewable Energy 80 (2015) 286e300

291

Fig. 3. Illustration of the application of WGP to probabilistic wind power forecasting.

permits the WGP to track to the time-varying characteristic of the wind generation, and reduce its computational costs. For further details on this online algorithm, the reader is referred to Ref. [28]. This way, we obtain a sparse online WGP forecasting model, referred to as SOWGP. Fig. 4 shows the examples of the probabilistic wind power forecasts obtained using SOWGP. From Fig. 4(a), one can see that the actual wind power values all fall within the regions predicted from the forecast distributions. It can also be seen that the predictive pdf is asymmetric in general, thus quite far from Gaussian. A

comparison between the predictive distribution and corresponding empirical distribution of the wind power is given in Fig. 4(b). Here the empirical distributions are obtained from the historical wind power series, e.g., in the middle subfigure, the actual wind speed is 10.02 m/s, then the empirical distribution of wind power in this figure is approximated by the historical wind power data points whose corresponding wind speed belong to [9.02, 11.02] m/s. One can see that, the predictive distributions are in general close to the empirical distributions. These results demonstrate that SOWGP yields reliable probabilistic forecasting performance. Therefore, it

Fig. 4. Probabilistic wind power forecasts given by SOWGP: (a) Non-Gaussian predictive pdf (solid line), 80%-level prediction interval (region between two vertical dotted lines) obtained from SOWGP, together with the actual wind generation measurements (*); (b) predictive pdf (solid line) and empirical distribution (histogram) of the wind power, when the actual wind power generation is relatively low (top), medium (middle), and high (bottom). The empirical distribution of wind power is represented by the historical wind power observations whose corresponding wind speed is in the region between two vertical dotted lines, “” denotes the empirical relation between wind speed and the power output. The data used for forecasting are obtained in the real operation of a wind farm in China.

292

P. Kou et al. / Renewable Energy 80 (2015) 286e300

can effectively quantify the uncertainties in wind power forecasts, i.e., the uncertainties associated with w.

a quantile of the cost, i.e., the value-at-risk in Ref. [37]. This new SMPC optimization problem has the following form

min

4. SMPC algorithm

uðkþhjkÞ;eðkþhjkÞ

J ¼b

4.1. Handling wind power forecast uncertainties using chance constraints

uðkþhjkÞ;eðkþhjkÞ

b

H X

eðk þ hjkÞ2 ; subject to

(17a)

h ¼ 1; 2; …; H: (17b)

Since the system output y(k þ hjk) is subjected to the random disturbances w, y(k þ hjk) is also a stochastic variable. In this case, the constraint (17b) can be replaced by a chance constraint

 i h   Pr yðk þ hjkÞ  yref ðk þ hjkÞ  eðk þ hjkÞ  a;

k¼1 H1 X

(19a) 2

uðk þ hjkÞ ; subject to

h¼0

k¼1

    yðk þ hjkÞ  yref ðk þ hjkÞ  eðk þ hjkÞ;

eðk þ hjkÞ2

þ ð1  bÞ

Once the uncertainties in wind power forecasts (i.e., the uncertainties associated with w) have been quantified, the next issue is to incorporate the quantified uncertainties into the controller design. The Chance constrained approach is very well suited for addressing this issue [35]. In optimization problems, the chance constraint is an efficient way to handle the uncertainties. This type of constraints explicitly contains stochastic parameters, which represent the uncertainties caused by disturbances. For an MPC controller in the presence of uncertain disturbance, if we enforce a chance constraint on its output, then this output can be, with a certain probability, kept within a reasonable bound around the reference trajectory [17]. This characteristic is similar in spirit to our control objective, that is, to maintain the combined power output y within a feasible range around the dispatch curve yref, in the presence of the uncertain wind power forecasts w. Furthermore, from (2), we see that y is linear in w, so the uncertainties in w can be completely reflected by the uncertainties in y. Therefore, the uncertainties in w can be incorporated into (8) by enforcing chance constraints on y. To achieve this, we first rewrite the first penalty term of (8a) in a probabilistic form. To do so, we introduce a set of auxiliary variables e(k þ 1jk), e(k þ 2jk), …,e(k þ Hjk), so that the penalty term P 2 b H h¼1 ðyðk þ hjkÞ  yref ðk þ hjkÞÞ can be recast as a squared summation of these auxiliary variables, together with a new constraint on y(k þ hjk) [36]. That is

min

H X

h ¼ 1; 2; …; H;

 i h   Pr yðt þ kjtÞ  yref ðt þ kjtÞ  eðt þ kjtÞ  a;

h ¼ 1; 2; …; H; (19b)

socmin  x2 ðk þ hjkÞ  socmax ;

h ¼ 1; 2; …; H;

(19c)

PB;max  uðk þ hjkÞ  PB;max ;

h ¼ 1; 2; …; H:

(19d)

Problem (19) is a chance-constrained optimization problem. Since e(k þ hjk) is a linear function of y(k þ hjk), minimizing the new cost function (19a) is equivalent to minimizing the original cost. As stated earlier, the uncertainty in w can be represented by the uncertainty in y. Therefore, as indicated by (19b), the SMPC problem formulation (19) takes into account the uncertainties in w, thus incorporating the wind power forecast uncertainties into the controller design. By this means, at every time stamp, the proposed SMPC controller computes the control actions based on the predictive distribution of wind power output, i.e., it optimizes the performance index based on a distribution of future wind power output. This way, the proposed SMPC controller can achieve better performance than the previous DMPC controllers, particularly when the wind power point forecast accuracy is poor. This is confirmed by the simulation results in Section 5. To facilitate the following derivation, we recast (19) in a matrix

 ~ ¼ uT eT T . To do so, form using an augmented decision variable u we introduce the following notations

2 3 2 3 3 b þ 0jkÞ uðk þ 0jkÞ xðk þ 1jkÞ wðk 6 wðk 7 6 6 xðk þ 2jkÞ 7 7 6 b þ 1jkÞ 7; 7; u ¼ 6 uðk þ 1jkÞ 7; w X¼6 5 b ¼4 4 5 4 5 « « « b þ H  1jkÞ uðk þ H  1jkÞ xðk þ HjkÞ wðk 2

(20) 3 2 3 3 2 eðk þ 1jkÞ yðk þ 1jkÞ yref ðk þ 1jkÞ 7 6 7 6 yðk þ 2jkÞ 7 6 7; yref ¼ 6 yref ðk þ 2jkÞ 7; e ¼ 6 eðk þ 2jkÞ 7; y¼6 5 4 5 5 4 4 « « « eðk þ HjkÞ yðk þ HjkÞ yref ðk þ HjkÞ 2

(18) where Pr[E] is the probability of event E occurring, and a2ð0; 1Þ is the confidence parameter. The chance constraint (18), which explicitly takes into account the uncertainties in y, is guaranteed to be satisfied with a prescribed probability a. This means that, we formulate the chance constraints on y as individual chance constraints, i.e., at each time stamp, (18) has to be individually fulfilled with the probability a. Using (17a) and (18), the nominal cost formulation in (8a) can easily be replaced by a probabilistic formulation, which explicitly takes the uncertainties into account. To do so, we simply replace P 2 the penalty term b H h¼1 ðyðk þ hjkÞ  yref ðk þ hjkÞÞ in (8) with (17a), and add the chance constraint (18). Consequently, the resulting probabilistic cost function can be regarded as minimizing

(21) 2 6 6 C¼6 4

6 7 6 7 7; J ¼ 6 4 5

c 1 c

cb ¼ ½ 0

2

3

c

2

6 6 1 ; Cb ¼ 6 4

:::

0

3

b1 «

1

« «

7 7 7; 5

AM2 b1

:::

b1

b1

0

Ab1 « AM1 b1

cb

3 7 7 7; 5

cb 1 cb

(22)

P. Kou et al. / Renewable Energy 80 (2015) 286e300

2 3 3 A b2 0 ::: 0 2 6 Ab2 6A 7 b2 « 7 ~ ~ w 6 7 7 w; 6 b G ¼ AxðkÞ þB 2 b ¼ 4 « 5xðkÞ þ 4 « « 1 « 5 H H1 H2 A b2 A b2 ::: b2 A 2

(23) Dy;U ¼ ½ CJ IH ; Dy;L ¼ ½ CJ IH ; Db ¼ ½ Cb J 0H ; DPb ¼ ½ IH 0H ;
G¼

bIH 0

 0 ; ð1  bÞIH

(24)

(25)

where IH denotes a H  H unit matrix, 0 is a matrix whose entries are all zero. Substituting (20)e(25) into (19), we then obtain T

~ Gu ~ ; subject to min J ¼ u

(26a)

~ ðkþhjkÞ u


  1    1 þ a e ~ ~ b r  CB e  CJu  C AxðkÞ þ y ;  Pr ½ w 2 ref 2 r r ¼ 1; 2; …; H; (26b)
  1    1 þ a e ~ ~ b r  CB  e  CJu  C AxðkÞ þ y ;  Pr ½ w 2 ref 2 r r ¼ 1; 2; …; H; (26c) h i ~ ~ r  socmax  Cb AxðkÞ ; ½Db u r

~ r  ½Db u

h

r ¼ 1; 2; …; H;

i ~  socmin þ Cb AxðkÞ ;

~ r  PB;max ; ½DPb u

r

Therefore, we need to recast this problem as a tractable deterministic one, i.e., to convert the chance constraints (26b) and (26c) to their deterministic equivalents. It is well known that if the disturbance is Gaussian distributed, then individual chance constraints can be analytically converted to the equivalent deterministic ones [38,39]. However, this Gaussian assumption is not valid for the wind power disturbance, which is clearly non-Gaussian [22]. Fortunately, this difficulty can be overcome by using warping function. In the proposed scheme, the non-Gaussian predictive distribution of wind power (i.e., the non-Gaussian predictive distribution of w) is provided by an SOWGP model. From Section 3, we see that in the SOWGP model, for each non-Gaussian predictive distribution pðw* jv* ; D Þ in original observation space, there exists a unique Gaussian counterpart pðz* jv* ; D Þ in latent space. That is, for each non-Gaussian predictive distribution (16), there is a Gaussian counterpart (15). In addition, there is a monotonic mapping between pðw* jv* ; D Þ and pðz* jv* ; D Þ, i.e., the warping function (12). From these two facts, if we map the chance constraints (26b) and (26c) to the latent space via the warping function, then the resulting mapped chance constraints in latent space will be equivalent to (26b) and (26c). Furthermore, in the mapped chance constraints, the stochastic parameter z is Gaussian distributed, so these new chance constraints can be analytically converted to deterministic ones. This way, the chance constrained optimization problem (26) is recast as a deterministic optimization problem, which is tractable both theoretically and practically. We now give the details of this conversion. First, we consider the chance constraint (26b). SOWGP is an online prediction model, so the warping function is also time-varying. For this reason, we denote the warping function at time stamp k as gk(w). Mapping the chance constraint (26b) to the latent space via gk(w), we obtain


  1   

~ ~ b r  gkþr e  CJu  CAxðkÞ þ yref CB Pr gkþr ½ w 2 r

r ¼ 1; 2; …; H;

r ¼ 1; 2; …; H;

~ r  PB;max ; ½DPb u

(26d)

293

r ¼ 1; 2; …; H:

1þa ; r ¼ 1; 2; …; H:  2

(26e)

(27)

(26f) (26g)

Here [m]r denotes the rth row of the vector m. As indicated by the chance constraints (26b) and (26c), this formulation explicitly incorporates the uncertainties associated with w, i.e., the uncertainties in wind power forecasts. In addition, compared to the problem formulation in (19), this matrix formulation is more convenient for the following derivations. So far, we have obtained the formulation of the SMPC optimization problem, which explicitly incorporates the wind power forecast uncertainties. This problem includes chance constraints on w, i.e., the chance constraints (26b) and (26c). Since an optimization problem can only be solved if all variables are deterministic, these chance constraints need to be converted into the equivalent deterministic constraints. However, as discussed in Section 3, the distribution of w is non-Gaussian. Therefore, the chance constraints (26b) and (26c) cannot be directly converted into deterministic ones. In the next subsection, we will show how to overcome this difficulty by using warping function.

4.2. Transformation of the non-Gaussian chance constraints The SMPC optimization problem (26) is a chance constrained stochastic programming problem, which is intractable in general.

Since gk(w) is monotonic, the mapped chance constraint (27) is equivalent to the original one in (26b). In the latent space, b r Þ is Gaussian distributed, so (27) can be recast as an zkþr ¼ gkþr ð½ w equivalent deterministic constraint, which is given by

gkþr




1 

~ e  CJu  CAxðkÞ þ yref   1þa ; r ¼ 1; 2; …; H;  F1 kþr 2 ~ CB 2

  r

(28)

where Fkþr denotes the cumulative density function (cdf) of zkþr, i.e., the r-step ahead probabilistic forecast in latent space. Since zkþr is Gaussian distributed, Fkþr can be easily obtained from (15). Unmapping (28) via the inverse warping function gk1 ðzÞ, we obtain




 ~ e  CJu  CAxðkÞ þ yref r    1 1 1 þ a ; r ¼ 1; 2; …; H:  gkþr Fkþr 2

~ CB 2

1 

1 ðF1 ðð1 þ aÞ=2ÞÞ; gU;a ðkÞ ¼ ½gkþ1 kþ1

Introducing

a

1 ðF1 ðð1 gkþ2 kþ2

1 ðF1 aÞ=2ÞÞ; …; gkþH kþH1 ðð1

þ

notation

be rewritten in a compact form

(29)

þ aÞ=2ÞÞT , then (29) can

294

P. Kou et al. / Renewable Energy 80 (2015) 286e300

h i ~ ~ g ðkÞ ; ½CJur  yref þ e  CAxðkÞ  CB 2 U;a r

r ¼ 1; 2; …; H: (30)

From (30), we have

h i ~ ~ g ðkÞ ; ½CJu  er  yref  CAxðkÞ  CB 2 U;a r

r ¼ 1; 2; …; H: (31)

Making use of the notation (24), we can rewrite (31) as a deter~ , in the ministic constraint on the augmented decision variable u form



h i  ~ ~ g ðkÞ ; ~  yref  CAxðkÞ  CB Dy;U u 2 U;a r

r ¼ 1; 2; …; H:

r

(32) In this way, the chance constraint (26b) is analytically converted to an equivalent deterministic constraint (33). This deterministic ~. constraint is a linear inequality, so it is a convex constraint on u Using the same approach, the chance constraint (26c) can also be converted to a deterministic counterpart, that is



h i  ~ ~ g ðkÞ ; ~   yref þ CAxðkÞ þ CB Dy;L u 2 L;a r r

r ¼ 1; 2; …; H: (33)

~. This is also a linear convex constraint on u So far, all the chance constraints in problem (26) have been converted to deterministic ones. Making use of these converted deterministic constraints, the stochastic optimization problem (26) can be recast as a tractable deterministic one. To achieve this, we replace the chance constraint (26b) and (26c) by the deterministic constraints (32) and (33), respectively. This leads to the following deterministic optimization problem T

~ Gu ~ ; subject to min J ¼ u

(34a)

~ ðkþhjkÞ u



h i  ~ ~ g ðkÞ ; ~  yref  CAxðkÞ  C B Dy;U u 2 U;a r h i  ~ ~ g ðkÞ ; ~   yref þ CAxðkÞ þ CB Dy;L u 2 L;a r r

In this section, we discuss the advantage of the proposed SMPC scheme over the conventional DMPC scheme. We first investigate the DMPC scheme. For a wind-battery system, a DMPC controller only considers the wind power point forecasts, which is the b i.e., Eð wÞ b ¼w b ¼ ½ wðk b þ 0jkÞ; wðk b þ 1jkÞ; …; expectation of w, b þ H  1jkÞT . Neglecting the operation constraints, the optiwðk mization objective of the DMPC scheme can be written as

min J ¼

uðkþhjkÞ

H  X

r ¼ 1; 2; …; H;

2 yðk þ hjkÞ  yref ðk þ hjkÞ

(35)

h¼1

Using the auxiliary variables e(k þ 1jk), e(k þ 2jk), …, e(k þ Hjk), (35) can be rewritten as

min

uðkþhjkÞ;eðkþhjkÞ

H X

eðk þ hjkÞ2 ; subject to

(36a)

k¼1

    yðk þ hjkÞ  yref ðk þ hjkÞ  eðk þ hjkÞ;

h ¼ 1; 2; …; H; (36b)

Notice that, although (36) has a similar form to (17), in (36) the wind power output is regarded as a deterministic variable. Substituting (20)e(25) into (36), this problem can then be formulated as

uðkþhjkÞ;eðkþhjkÞ

(34b)

5. Discussion

min

r ¼ 1; 2; …; H;

r

applied to the wind-battery system. At the next time stamp, this procedure is repeated using the new actual wind power measurement. That is, computing PB(k þ 1) and EB(k þ 1) using the actual wind power measurement PW(k). By this receding horizon approach, feedback is introduced into the controller design, since the new optimization problem for the next time stamp is a function of the new measured states at that point in time.

H X

eðk þ hjkÞ2 ; subject to

(37a)

k¼1

i h ~ ~ w  eðk þ rjkÞ; CJu þ CAxðkÞ  yref þ CB 2b r

r ¼ 1; 2; …; H; (37b)

(34c) i h ~ ~ r  socmax  Cb AxðkÞ ½Db u ; r

~ r  ½Db u

h

r ¼ 1; 2; …; H;

~ r  PB;max ; ½DPb u

i h ~ ~ w  eðk þ rjkÞ; CJu þ CAxðkÞ  yref þ CB 2b r

r ¼ 1; 2; …; H: (37c)

i

~  socmin þ Cb AxðkÞ ;

~ r  PB;max ; ½DPb u

(34d)

r

r ¼ 1; 2; …; H; r ¼ 1; 2; …; H:

r ¼ 1; 2; …; H;

(34e) (34f) (34g)

~ , in In problem (34), the cost function is a quadratic function of u ~. addition, all constraints are linear inequalities in terms of u Therefore, problem (34) is a convex quadratic optimization problem, and thus can be solved using standard quadratic programming algorithms. At every time stamp, solving problem (34) yields an optimal control sequence. The first control action from this sequence is

Solving problem (37), we obtain a control sequence u* ¼ ½u* ðk þ 0jkÞ; u* ðk þ 0jkÞ; …; u* ðk þ H  1jkÞT and a set of tracking error bounds e* ¼ ½e* ðk þ 0jkÞ; e* ðk þ 0jkÞ; …; e* ðk þ H  1jkÞT . The first control action u* ðk þ 0jkÞ is implemented. However, in actual applications, the accuracy of the wind power point forecast cannot be guaranteed and fairly low on average. Therefore, there may be significant differences between the actual b þ hjkÞ. Assuming wind power output w(k þ hjk) and its forecast wðk that the actual wind power output at time stamp k þ 1 is wr(k þ 1jk), which is significantly higher than its forecast b þ 1jkÞ. In this circumstance, the control action u* ðk þ 0jkÞ may wðk lead to the violation of the following constraint

P. Kou et al. / Renewable Energy 80 (2015) 286e300

i h i h ~ ~ wr ðk þ 1jkÞ  e* ðk þ 1jkÞ: ½CJk u* ðk þ 0jkÞ þ CAxðkÞ þ CB 2 r

r

(38) This violation means that, the combined power output y(k þ 1jk) will greatly deviate from the desired trajectory yref(k þ 1jk), thus leading to a poor tracking performance. Conversely, in the proposed SMPC scheme, the uncertainties b are taken into account in the controller design. associated with w This way, the constraints (37) become

h

~ ~ g CJu þ CAxðkÞ  yref þ CB 2 U;a

i r

 eðk þ rjkÞ;

r ¼ 1; 2; …; H; (39a)

h

~ ~ g CJu þ CAxðkÞ  yref þ CB 2 L;a

i r

 eðk þ rjkÞ;

r ¼ 1; 2; …; H: (39b)

Here gU,a can be regarded as the upper (1 þ a)/2 fractile of the wind power predictive distribution. Since the point forecast value b þ hjkÞ is the mean of this distribution, gU,a is larger than wðk b þ hjkÞ. Similarly, gL,a, being the lower (1 þ a)/2 fractile of the wðk b þ hjkÞ. As a result, for a predictive distribution, is smaller than wðk same tracking error bound e* ðk þ hjkÞ, the constraint (39) is stricter than the constraints (38). In other words, with the proposed SMPC scheme, as long as the actual wind power output lies between gL,a and gU,a, the constraint (39) will be satisfied. Thus, a small tracking error is achieved. In contrast, with the conventional DMPC, the constraint (38) will only be satisfied when the actual wind power output is very close to its forecast. In this sense, the proposed SMPC scheme is more likely to achieve small tracking errors than the DMPC scheme. Of course, this improvement in tracking performance comes at the cost of a higher control effort. In terms of a wind-battery system, this means a deeper depth of charge/discharge. In spite of this, the proposed SMPC scheme sill yields a good compromise between tracking performance and control efforts, as is supported by the simulation results in Section 6. 6. Simulation results In this section, in order to verify the effectiveness of the proposed SMPC scheme, several simulations are carried out using an actual wind farm data. For comparison purposes, a conventional DMPC scheme is also simulated under the same conditions. All simulations are implemented in Matlab 7.8, run on a personal computer with 2.52 GHz CPU and 4 GB memory. 6.1. Simulation setup An actual wind farm, geographically located in Inner Mongolia, China, is considered in the following simulations. The nominal capacity of this wind farm is 120 MW. It comprises 80 Goldwind 1.5 MW wind turbines, with a hub height of 65 m. In our simulations, the 5-min wind power data from this wind farm is used as PW, i.e., the original wind power outputs. The data from September 14th, 2011 to September 21st, 2011 are used, i.e., the whole simulation period is one week, which covers 2016 data points. Among various types of BESS, the sodiumesulfur (NaS) battery technology appears to be promising for integration with wind farms, due to its high efficiency, high energy capacity, and long life span [13,40,41]. For these reasons, a 48 MWh NaS battery bank is assumed for this simulation. Based on the information from the

295

battery manufacturer [42] and the discussion in Section 2, the BESS charge/discharge efficiency h is set to 90%, the battery's power rating constraint is set to ±24 MW, the lower and upper SOC limits are selected as 20% and 80%, respectively. Regarding the SMPC controller, the prediction horizon is chosen to be 1e12 h, and the sampling time is set to 5 min. For the chance constraints (18), the probability of constraint satisfaction is set to 99%. For the prediction model SOWGP, its explanatory variables contain the historical local wind speed observations, the historical wind speed observations from multiple nearby sites, as well as the NWP data. The NWP data is obtained from the NCEP/NCAR reanalysis dataset [43]. As suggested in Refs. [10,29], the reference trajectory yref, i.e., the dispatched target curve PD, is selected as the next 4 h average of the wind power point forecasts, which is provided by SOWGP. Correspondingly, each dispatch interval is assumed to be 4 h. The key feature of the proposed SMPC controller is that it incorporates the wind power forecast uncertainties, so we compare it with a conventional DMPC controller, which does not consider such uncertainties. This way, we examine whether incorporating these uncertainties would improve the control performance. The primary goal of our controller is to improve the dispatchability of a wind farm, i.e., to minimize the deviations between the combined power output PG and the dispatched curve PD. To assess these deviations, we select a tracking performance criterion as follows

PI ¼

N X

xn jPG ðnÞ  PD ðnÞj:

(39)

n¼1

This criterion is based on the performance index in Ref. [10], it sums the unacceptable power deviations that are larger than a certain threshold. Here xn represents the occurrence of the unacceptable deviations. Assuming that the deviations up to ±6 MW (i.e., ±5% of Pn) are acceptable, then xn is set to 1 if jPG(n)  PD(n)j < 6 MW and zero otherwise. Moreover, we denote the number of unacceptable deviations (i.e., the number of data points that satisfy jPG(n)  PD(n)j  6 MW) by ND. 6.2. Probabilistic wind power forecasts The probabilistic wind power forecast is the stepping stone for our SMPC scheme. Therefore, before proceeding with the control system simulation, we first evaluate the forecasting performance of SOWGP, which provides the probabilistic wind power forecasts in the proposed scheme. Fig. 5 shows the 4-h-ahead forecasting results given by SOWGP, in the form of a fan chart, i.e., a set of prediction intervals with increasing confidence level and fading color. Notice that, to visually illustrate the probabilistic forecasting results, we plot the prediction intervals rather than the predictive pdf. Furthermore, in order to avoid visual clutter, this figure does not show the results for the whole simulation period (i.e., one week), but only the results for the first two days. These prediction intervals depict the ranges of possible conditional range metric values for future wind power data. The color of each interval varies with the respective probability of the wind power lying within the interval. From this figure, we can see that almost all actual wind power measurements fall within the 60%-level prediction intervals, and the prediction intervals are generally not symmetric around the point predictions. Furthermore, the width of the prediction intervals is clearly influenced by the volatility of the wind power. When the generation fluctuates sharply, the prediction intervals tend to become wider. Conversely, when the generation is fairly flat, the corresponding prediction intervals are tighter.

296

P. Kou et al. / Renewable Energy 80 (2015) 286e300 80

90% 80% 70% 60% 50% 40% 30% 20% 10% mean actual

Power (% of Pn)

60

40

20

0 0

10

20

30

40

50 Time (hour)

60

70

80

90

Fig. 5. Four-hours-ahead wind power forecasting results given by SOWGP, together with the actual measured values.

Table 1 Probabilistic and point forecasting performance of SOWGP. Horizon (h)

1

2

3

4

6

12

NLPD MAPE (%)

0.85 10.18

0.88 11.85

0.91 12.98

0.93 13.40

0.95 14.74

1.08 16.65

Table 1 reports both the probabilistic and the point forecasting results. Here we use the average negative log predictive density (NLPD) [44] to evaluate the quality of the probabilistic forecasting results. The NLPD is a means of evaluating the amount of probability that the model assigns the target. It penalizes predictions that are either under or over confident. Moreover, since SOWGP also provides point forecasts, we adopt the mean absolute percentage error (MAPE) to evaluate its accuracy. From this table, one can see that SOWGP provides reliable probabilistic forecasts, as well as satisfactory point predictive accuracy. These results demonstrate that, the probabilistic forecasts provided by SOWGP effectively quantify the uncertainties in wind power forecasts. The main reason for this is that, SOWGP captures the non-Gaussian characteristics of the wind power uncertainties, as well as the time varying nature of wind generation.

6.3. Trade-off between two control objectives In this subsection, we study the influence of the trade-off parameter b. This parameter is introduced to balance the aforementioned two control objectives, i.e., the minimization of tracking error and the minimization of the amount of BESS charging/discharging power. To investigate the influence of b, we vary its value

Fig. 7. PI and ND of the proposed SMPC controller, as a function of b.

from 0.05 to 0.95 at a step of 0.05, and set the prediction horizon to 4 h. The results of this investigation are shown in Figs. 6e8. As stated earlier, our primary control objective is to make the combined power output follows the dispatch curve as closely as possible. Fig. 6 presents the influence of b on the achievement of this objective. In this figure, the combined power output PG with different values of b is plotted, together with the dispatch target curve PD. Notice that, for sake of clarity, we only plot the case of b ¼ 0.05, 0.5, 0.8, and 0.95. For the same reason, this figure does not show the results for the whole simulation period (i.e., one week), but only the results for the first two days. It can be shown that, when b > 0.8, the combined power output closely follows the

80 Dispatched level P D P G, β=0.05 P G, β=0.5

Power (% of Pn)

60

P G, β=0.8 P G, β=0.95

40

20

0

0

10

20

30

40

50 Time (hour)

60

70

80

90

Fig. 6. Combined power output PG with different values of b, together with the dispatch target curve PD.

P. Kou et al. / Renewable Energy 80 (2015) 286e300

297

Fig. 8. BESS performance as a function of b: (a) BESS charging/discharging power, (b) sum of the squares of the BESS charging/discharging power, together with the standard deviation of the SOC.

that the charge/discharge cycle is approximately every 3e4 h, but most of the time, both charging and discharging are shallow. From the above results, it is found that b ¼ 0.8 represents a good compromise between the aforementioned two objectives. Therefore, we set b to 0.8 in the following simulations.

desired dispatch curve in most of the time. Fig. 7 depicts the evolution of PI and ND, as a function of b. It is obvious that both of them decrease rapidly when b<0.8 and then levels off. When b  0:8, the value of ND is about 30, which is relatively very small compared to the length of the simulation period (i.e., 2016 data points). Our secondary control objective is to minimize the control effort PB, i.e., to minimize the amount of energy charged in or discharged from the BESS. Fig. 8 illustrates the influence of b on the achievement of this objective. In Fig. 8(a), the evolutions of PB are plotted when b ¼ 0.05, 0.5, 0.7, and 0.95. We see that, as b increases, PB tends to fluctuate more vigorously. This is confirmed by Fig. 8(b), P which shows the PB ðkÞ2 and std(EB) as a function of b. In addition, Fig. 8(a) shows that the BESS charging/discharging power is kept between 24 MW and 24 MW, as desired. The figure also indicates

6.4. Comparison with DMPC In this subsection, we compare the proposed SMPC scheme to the conventional DMPC scheme. We vary the prediction and control horizon H2f4; 8; 12; 16; 24; 48g, i.e., 1 h, 2 h, 3 h, 4 h, 6 h, and 12 h. In DMPC, the control actions are computed based on the wind power point forecast, which is also provided by SOWGP. Specifically, we use the mean of the SOWGP predictive distributions as

100 90% prediction interval Wind power output PW

Power (% of Pn)

80

Dispatched level P D Net injected power PG

60 40 20 0

0

10

20

30

40

50 Time (hour)

60

70

80

90

(a) 100 Point wind power forecast Wind power output PW

Power (% of Pn)

80

Dispatched level PD Net injected power P G

60 40 20 0

0

10

20

30

40

50 Time (hour)

60

70

80

90

(b) Fig. 9. Dispatching of wind farm power with BESS, using (a) SMPC and (b) DMPC, with prediction horizon H ¼ 4 h.

298

P. Kou et al. / Renewable Energy 80 (2015) 286e300

50 DMPC SMPC

40 30

6MW, upper bound for acceptable deviation

Power (% of Pn)

20 10 0 -10 -20

-6MW, lower bound for acceptable deviation

-30 -40 -50

0

20

40

60

80 100 Time (hour)

120

140

160

Fig. 10. Power tracking error of SMPC and DMPC, with prediction horizon H ¼ 4 h.

the point forecasts. The results of this comparison are reported in Figs. 9e13 and Table 2. Taking as an example that H ¼ 4 h, Fig. 9 shows the comparison between SMPC and DMPC in terms of tracking performance, thus comparing their effectiveness in improving the dispatchability of a wind farm. In order to avoid cluttering, this figure does not show the results for the whole simulation period (i.e., one week), but only the results for the first two days. Fig. 9(a) shows the results of the proposed SMPC scheme. Here, the shaded region represents the 90%-level prediction interval of wind power generation, i.e., the probabilistic wind power forecasts. It is seen that, almost all actual wind power measurements fall within this prediction interval. Based on these probabilistic forecasts, the control actions of the proposed SMPC controller are computed. At first glance, one can see that, by implementing the computed control actions, the combined power output PG becomes smoother than the original wind power output PW. More significantly, after an initial transient period, the combined power output PG follows the desired dispatch curve PD perfectly most of the time, except from hour 16 to hour 25 when the tracking error is relatively large. The reason for these relatively large tracking errors is twofold. First, during hour 16 to hour 25, the original wind power output fluctuates sharply, as indicated by PW (dotted line in Fig. 9(a)). Second, as an example,

consider the maximum tracking error occurs at hour 20. One can see that, although the actual wind power measurement at this time stamp falls within the 90%-level prediction interval, it lies at the edge of the interval, so the confidence for this prediction is low. For comparison, the result of DMPC is shown in Fig. 9(b). Here, the dashed line denotes the wind power point forecasts, based on which the control law is computed. It is obviously that, the tracking performance of DMPC highly depends on the point forecast accuracy. When the wind power point forecast is accurate, the combined power output closely follows the desired dispatch curve, e.g., from hour 33 to hour 50, and hour 54 to hour 73. Conversely, when the point forecast accuracy is poor, the tracking performance of DMPC is significantly degraded, e.g., at hour 20, hour 53, and hour 57. In Figs. 10 and 11, we compare the power tracking error of SMPC and DMPC in detail. Fig. 10 shows the tracking error of SMPC and DMPC, i.e., the deviations between PG and PD. One can see that, the overall tracking performance of SMPC is quite satisfactory. It is observed that, after an initial transient, almost all tracking errors are within the acceptable region (i.e., ±6 MW). The initial transient is due to the initial value of PG, i.e., we set PG(0) ¼ 0. Correspondingly, Fig. 11 presents the histogram of the tracking errors. As Fig. 11(a) illustrates, for the proposed SMPC controller, about 70% of the tracking errors are within ±2 MW, and larger errors rarely occur. The lower outliers in this histogram are caused by the initial transient, as discussed earlier. Fig. 11(b) shows the histogram of the tracking errors of DMPC. At first sight, this histogram is much flatter than that of the SMPC, the central part of the distribution is lower and its tails longer. This indicates that, for the conventional DMPC controller, 89% tracking errors are within the acceptable region (i.e., ±6 MW), only 50% of the tracking errors are within ±2 MW. Moreover, its maximum power tracking error reaches up to 11.97% of Pn after the initial transient. The result in Fig. 9 shows that, the combined power output PG becomes smoother than the original wind power output after mitigation by the BESS. Fig. 12 illustrates this smoothing effect in detail. In this figure, we use the maximum fluctuation of PG in each dispatch interval (i.e., 4 h) as the smoothness measure, and then compare the smoothness properties of the proposed SMPC, the conventional DMPC, and the original wind power output PW. It is obvious to see that, compared to the original wind power output PW, both the SMPC and the DMPC significantly reduce the power fluctuations, after an initial transient period. In addition, compared to the DMPC, the power output of SMPC is smoother. Specifically, for the proposed SMPC controller, the maximum fluctuation is as low as 3% of Pn, except for some outliers. In contrast, for the

Fig. 11. Histogram of the power tracking error, (a) SMPC and (b) DMPC, with prediction horizon H ¼ 4 h.

P. Kou et al. / Renewable Energy 80 (2015) 286e300

299

Table 2 Control performance for different prediction horizons.

Fig. 12. Maximum output power fluctuations in each dispatch interval, with prediction horizon H ¼ 4 h.

conventional DMPC controller, the maximum fluctuations exceed 5% of Pn in many cases. Notice that, for both SMPC and DMPC, the maximum power fluctuations in the fifth dispatch interval are significantly larger than those in other intervals. The reason for this lies in the poor wind power forecast at hour 20, as indicated in Fig. 5. Fig. 13 shows the evolution of BESS SOC. One can see that, due to the SOC constraint (8b), after an initial transient the SOC of BESS is always maintained between 20% and 80%, as desired. As a result, the overcharge and overdischarge are avoided, and thus the operational lifetime of BESS is prolonged. Notice that, compared to the conventional DMPC, the proposed SMPC controller makes the SOC fluctuation greater. However, it should be noted that there is a trade-off between power tracking performance and SOC fluctuation, so improving the tracking performance will generally increase the SOC fluctuation, and vice versa. In this sense, one can say that the proposed SMPC scheme offers a good compromise between these two objectives, due to the fact that the primary objective of a wind-battery system is to improve the dispatchability of a wind farm.

Prediction horizon (h)

1

2

3

4

6

12

SMPC

PI (% of Pn) ND

4.55 46

4.53 42

4.43 39

4.34 36

4.33 36

4.40 40

DMPC

PI (% of Pn) ND

6.19 80

6.05 78

5.93 76

5.93 76

5.92 76

5.96 79

Finally, we assess the performance of SMPC with different prediction horizons. The results are reported in Table 2. From this table, it is obvious that the performance of SMPC is better than that of DMPC. With the proposed SMPC controller, the number of unacceptable tracking errors does not exceed 50 in all scenarios. In contrast, with the conventional DMPC controller, this number is greater than 70. This demonstrates the high performance of the proposed SMPC scheme in improving wind power dispatchability. Furthermore, when the prediction horizon is less than 6 h, for both SMPC and DMPC, it appears that the longer the prediction horizon, the better the performances are. The reason for this lies in the fact that, by considering a longer prediction horizon, the controller can better oversee the consequences of its actions [45]. However, it is also noted that, when the prediction horizon is greater than 6 h, the control performance is deteriorated. The reason for this may be that, the wind power forecasting performance decreases significantly as the prediction horizon exceeds 6 h, as indicated in Table 1. 7. Conclusion In this work, a novel stochastic MPC controller is established to improve the wind power dispatchability and reduce its fluctuation. With this controller, the wind farm power output can be brought to a desired dispatch curve, so that a wind farm can be dispatched like a conventional generation plant. This controller employs a probabilistic wind power forecasting model as the prediction model, which quantifies the uncertainties in wind power forecasts. Using chance constraints, the quantified uncertainties are incorporated into the controller design. This way, the proposed controller handles the non-Gaussian uncertainties in wind power forecasts. Simulation results on actual wind farms show that the performance of the proposed SMPC controller is satisfactory, it follows the desired dispatch curve quite closely. Meanwhile, it keeps the SOC and the power output of the BESS within its operational limits, thus avoiding overcharging and depletion. Comparison with the conventional deterministic MPC controller indicates that, the proposed controller can achieve better performance in terms of power tracking errors, thus demonstrating the benefits of incorporating wind power forecast uncertainties. Future work includes developing a more precise model for BESS, investigating the relation between control performance and probabilistic forecast quality, and adding the switching constraints to reduce the number of charge/discharge cycles. Furthermore, the principle of the proposed scheme can also be applied to some other energy storage systems. In summary, the proposed scheme is an alternative with some competitive advantages regarding the existing techniques. Although it is not yet a sophisticated control system, the described methodology can be integrated into currently existing systems, and serve as a stepping stone for further progress. Acknowledgments

Fig. 13. SOC curves with the SMPC and DMPC controllers, with prediction horizon H ¼ 4 h.

The authors would like to thank Dr. Edward Snelson for providing the code for WGP. The authors also gratefully acknowledge the financial support received from the National Natural

300

P. Kou et al. / Renewable Energy 80 (2015) 286e300

Science Foundation of China (Project Nos. 61403303, 51177125, 60974101, 60736027, 61221063, 51477130), the China Post-Doctoral Science Foundation (Projects No. 2014M560776), and the State Key Laboratory of Electrical Insulation and Power Equipment (Project No. EIPE15306).

References [1] World wind energy association, half-year report. 2012. Available online at: http://www.wwindea.org/webimages/Half-year_report_2012.pdf. [2] Wang J, Botterud A, Bessa R, Keko H, Carvalho L, Issicaba D, et al. Wind power forecasting uncertainty and unit commitment. Appl Energy 2011;88(11): 4014e23. [3] Parsons B, Milligan M, Zavadil B, Brooks D, Kirby B, Dragoon K, et al. Grid impacts of wind power: a summary of recent studies in the United States. Wind Energy 2004;7(2):87e108. [4] Divya K, Østergaard J. Battery energy storage technology for power systems e an overview. Electr Power Syst Res 2009;79(4):511e20. [5] Hadjipaschalis I, Poullikkas A, Efthimiou V. Overview of current and future energy storage technologies for electric power applications. Renew Sustain Energy Rev 2009;13(6):1513e22. [6] Jayasinghe S, Vilathgamuwa D, Madawala U. Direct integration of battery energy storage systems in distributed power generation. IEEE Trans Energy Convers 2011;26(2):677e85. [7] Gao M, Hui D, Gao Z, Lei W, Li J, Wang Y. Presentation of national wind/ photovoltaic/energy storage and transmission demonstration project and analysis of typical operation modes. Autom Electr Power Syst 2013;37(1):59e64. [8] Wang H, Tang Y, Chen J, Hou J. Modeling and equivalence of integrated power generation system of wind, photovoltaic and energy storage. In: Proc. of AsiaPacific power and energy engineering conference, Shanghai, China; March 2012. [9] Venayagamoorthy G, Rohrig K, Erlich I. One step ahead: short-term wind power forecasting and intelligent predictive control based on data analytics. IEEE Power Energy Mag 2012;10(5):70e8. [10] Teleke S, Baran M, Bhattacharya S, Huang A. Optimal control of battery energy storage for wind farm dispatching. IEEE Trans Energy Convers 2010;25(3): 787e94. [11] Khalid M, Savkin V. A model predictive control approach to the problem of wind power smoothing with controlled battery storage. Renew Energy 2010;35(7):1520e6. [12] Khalid M, Savkin V. An optimal operation of wind energy storage system for frequency control based on model predictive control. Renew Energy 2012;48(12):127e32. [13] Khatamianfar A, Khalid M, Savkin V, Agelidis G. Improving wind farm dispatch in the Australian electricity market with battery energy storage using model predictive control. IEEE Trans Sustain Energy 2013;4(3):745e55. [14] Qi W, Liu J, Chen X, Christofides D. Supervisory predictive control of standalone wind/solar energy generation systems. IEEE Trans Control Syst Technol 2011;19(1):199e207. [15] Qi W, Liu J, Christofides D. Distributed supervisory predictive control of distributed wind and solar energy systems. IEEE Trans Control Syst Technol 2013;21(2):504e12. [16] Pinson P, Kariniotakis G. Conditional prediction intervals of wind power generation. IEEE Trans Power Syst 2010;25(4):1845e56. [17] Schwarm T, Nikolaou M. Chance-constrained model predictive control. AIChE J 1999;45(8):1743e52. [18] Couchman D, Cannon M, Kouvaritakis B. Stochastic MPC with inequality stability constraints. Automatica 2006;42(12):2169e74. [19] Kouvaritakis B, Cannon M, Rakovic V, Cheng Q. Explicit use of probabilistic distributions in linear predictive control. Automatica 2010;46(10):1719e24. [20] Nikhil C, Larry H. Forecasting the reliability of wind-energy systems: a new approach using the RL technique. Appl Energy 2012;96(8):422e30.

[21] Bessa R, Miranda V, Botterud A, Zhou Z, Wang J. Time-adaptive quantilecopula for wind power probabilistic forecasting. Renew Energy 2012;40(1): 29e39. [22] Lange M. On the uncertainty of wind power predictionsdanalysis of the forecast accuracy and statistical distribution of errors. J Sol Energy Eng 2005;127(2):177e84. [23] Ma Y, Borrelli F. Fast stochastic predictive control for building temperature regulation. In: Proc. of American control conference, Montreal, Canada; June 2012. [24] Ono M. Joint chance-constrained model predictive control with probabilistic resolvability. In: Proc. of American control conference, Montreal, Canada; June 2012. [25] Oldewurtel F, Jones N, Morari M. A tractable approximation of chance constrained stochastic MPC based on affine disturbance feedback. In: Proc. of IEEE conference on decision and control, Cancun, Mexico; December 2008. [26] Li P, Wendt M, Wozny G. A probabilistically constrained model predictive controller. Automatica 2002;38(7):1171e6. [27] Snelson E, Rasmussen C, Ghahramani Z. Warped Gaussian processes. In: Advances in neural information processing systems, vol. 14. Cambridge, MA: MIT Press; 2003. p. 337e44. [28] Kou P, Gao F, Guan X. Sparse online warped Gaussian process for wind power probabilistic forecasting. Appl Energy 2013;108(8):410e28. [29] Ahlstrom M, Jones L, Zavadil R, Grant S. The future of wind forecasting and utility operations. IEEE Power Energy Mag 2005;3(6):57e64. [30] Maciejowski M. Predictive control with constraints. London: Prentice Hall; 2002. [31] Jiang Q, Gong Y, Wang H. A battery energy storage system dual-layer control strategy for mitigating wind farm fluctuations. IEEE Trans Power Syst 2012;28(3):3263e73. [32] Giebel G, Brownsword R, Kariniotakis G, Denhard M, Draxl C. The state of the art in short-term prediction of wind power: a literature overview. 2nd ed. Roskilde, DK: Riso National Lab.; 2011 Available online at: http://anemos.cma. fr [accessed Jan. 2011]. [33] Costa A, Crespo A, Navarro J, Lizcano G, Madsen H, Feitosa E. A review on the young history of the wind power short-term prediction. Renew Sustain Energy Rev 2008;12(6):1725e44. [34] Wu J, Guan X, Zhou X, Zhou Y. Estimation and characteristic analysis of aggregated generation of geographically distributed wind farms. In: Proc. of 2011 IEEE power and energy society general meeting, Detroit, MI, USA; July 2011. [35] Charnes A, Cooper W. Deterministic equivalents for optimizing and satisficing under chance constraints. Oper Res 1963;11(1):18e39. [36] Schildbach G, Calafiore GC, Fagiano L, Morari M. Randomized model predictive control for stochastic linear systems. In: Proc. of American control conference, Montreal, Canada; June 2012. [37] Shapiro A, Dentcheva D, Ruszczynski A. Lectures on stochastic programming: modeling and theory. Philadelphia: SIAM; 2009. [38] Prekopa A. Stochastic programming. 3rd ed. Dordrecht, Netherlands: Kluwer Academic Publishers; 1995. [39] Calafiore C, Laurent G. Linear programming with probability constraints e part 1. In: Proc. of American control conference, New York, NY, USA; July 2007. [40] Yao D, Choi S, Tseng K, Lie T. A statistical approach to the design of a dispatchable wind power-battery energy storage system. IEEE Trans Energy Convers 2009;24(4):916e25. [41] Shakib D, Spahic E, Balzer G. Computation dimensioning of a sodium-sulphur battery as a backup of large wind farm. In: Proc. of the 16th power system computation conference, Glasgow, Scotland; July 2008. [42] Sodium sulfur battery catalog. Mizuho, Nagoya, Japan: NGK Insulators Ltd.; 2009. [43] Kalnay E, Kanamitsu M, Kistler R, Collins W, Deaven D, Gandin L, et al. The NCEP/NCAR 40-year reanalysis project. Bull Am Meteorol Soc 1996;77(3): 437e70. [44] Good I. Rational decisions. J R Stat Soc Ser B 1952;14(1):107e14. [45] Negenborn R. Multi-agent model predictive control with applications to power networks [Ph.D. thesis]. Delft University of Technology; 2007.