Predictive control of wind turbine for load reduction during ramping events

Electrical Power and Energy Systems 93 (2017) 135–145

Contents lists available at ScienceDirect

Electrical Power and Energy Systems journal homepage: www.elsevier.com/locate/ijepes

Predictive control of wind turbine for load reduction during ramping events Weipeng Liu a, Changgang Li a, Yutian Liu a,⇑, Qiuwei Wu a,b a b

Key Laboratory of Power System Intelligent Dispatch and Control of Ministry of Education, Shandong University, 17923 Jingshi Road, Jinan 250061, China Centre for Electric Power and Energy, Department of Electrical Engineering, Technical University of Denmark, Kgs. Lyngby 2800, Denmark

a r t i c l e

i n f o

Article history: Received 6 January 2017 Received in revised form 26 April 2017 Accepted 14 May 2017

Keywords: Predictive control Wind turbine load Wind power ramping events Multi-objective optimization Neural network prediction

a b s t r a c t With increasing penetration of wind power, the impact of its intermittence and volatility on power systems becomes more severe. A predictive control strategy for wind turbines (WTs) is proposed to deal with wind power ramping events and reduce WT load on the blades. The blade load model is based on the Blade Element Momentum (BEM) theory. The generator speed and pitch angle are simultaneously regulated to realize the control objectives. A two-stage optimization is designed in order to reduce the computational complexity. The objectives of the first stage are minimizing the ramping rate and maximizing the power generation. A trade-off is made between the two contradictory objectives by setting weight coefficients. The second stage reduces the WT load and meanwhile guarantees the power reference from the first stage is tracked. Feedback is designed based on neural network prediction to compensate the error of the prediction model. Case studies with a 1.5 MW WT were conducted to demonstrate the efficacy of the proposed predictive control strategy. Simulation results show that the proposed control can reduce the WT load during ramping events and the risk of ramping events. Ó 2017 Published by Elsevier Ltd.

1. Introduction Wind power has been developing very fast over the world. With high wind power penetration, the impact of wind power fluctuation on power systems becomes more severe [1–3]. Under some serious circumstances, a ramping event will lead to frequency deviation and influence security of power systems [4]. Apart from regulating the output of conventional generators, it is sometimes necessary to use demand side response services, activate reserves and curtail wind power output [5–7]. These operation actions inevitably increase system operation cost. With high penetration of wind power and weak grid interconnection, the negative effect is more severe and the risk of tripping incidents and load shedding is significantly higher. On February 26, 2008, a large ramp-down event of wind power generation occurred in Texas, together with a quick load ramp-up. The Electric Reliability Council of Texas (ERCOT) had to curtail 1108 MW loads [8] to recover the system frequency. The ramping event in Texas shows that it is crucial to deal with wind power ramping events. Research on wind power ramping events mainly focuses on two areas: wind power ramping events ⇑ Corresponding author. E-mail addresses: [email protected] (W. Liu), [email protected] (C. Li), [email protected] (Y. Liu), [email protected] (Q. Wu). http://dx.doi.org/10.1016/j.ijepes.2017.05.025 0142-0615/Ó 2017 Published by Elsevier Ltd.

prediction and wind power ramping control. Wind power output control can be divided at three levels: wind turbine (WT) control, wind farm (WF) control and wind farm cluster (WFC) control [9]. The WT control includes operation stability improvement [10,11], damping control [12], etc. The existing WT power ramping control strategies can be divided into two categories. The first category is associated with an energy storage system to smooth output power [13]. The construction and maintenance cost of such a system is usually high. The second category utilizes a power rate limit to reduce the ramping rate [14]. In these methods, WT load reduction has not been considered. As the WT size increases, its structure becomes lighter and more flexible to reduce the cost. Consequently, load reduction has become an important objective of WT control and design [15,16]. WTs will be shut down at the cut-out speed due to heavy load, resulting in a large-scale ramp-down event. Hence, it is necessary to combine WT ramping control and WT load reduction. This paper intends to decrease wind power intermittence by WT ramping control with load reduction. Studies on WT load cover several aspects such as extreme load reduction [17], asymmetric load compensation [16,18], tower load reduction [19] and vibration reduction [20]. Pitch control is the most effective way to reduce WT load although it decreases output power of WTs. In other words, minimizing load and maximizing power generation are two contradictory objectives. To strike a balance between the

136

W. Liu et al. / Electrical Power and Energy Systems 93 (2017) 135–145

two objectives, control of the generator speed and pitch angle needs to be optimized together. A number of methods have been developed to control the generator speed and pitch angle. The gain selection of proportional-integral-derivative (PID) control is based on a linearized model at one working point [21]. Its control performance will deteriorate in real operation due to the highly nonlinear characteristic of WTs. The sliding mode control (SMC) [18] and fuzzy logic control (FLC) [22] are suitable to handle nonlinear systems because they require little model information. The main drawback of the SMC is the chattering. The stability of a fuzzy system is difficult to be guaranteed due to lack of mathematical description [23]. Another effective method to control WTs is predictive control. With the receding horizon principle, it can control the generator speed and pitch angle in a wide operation region with input and output constraints. The predictive control of WTs was studied in [24] and [25]. An offset-free predictive control was proposed in [24] to guarantee the maximum power tracking under rated wind speed. Ref. [25] proposed a multiple model predictive control to improve power quality and increase lifetime of WTs by mitigating transient load of the drivetrain. To reduce the computational complexity, the prediction models were linearized in both [24] and [25]. To compensate the model-plant mismatch caused by linearization, an observer was applied to the state space model in [24] and a prediction model bank with scheduling signal was designed in [25]. However, the blade load reduction has not been studied in predictive control. In this paper, a predictive control strategy of WTs with blade load reduction is proposed to deal with ramping events. There are three objectives: minimizing the ramping rate, maximizing power generation and minimizing the WT load. The blade load model is based on the blade element moment (BEM) theory [26]. The main contribution of this paper is summarized as follows: (A) develop a predictive control strategy for WTs to strike a balance among the three objectives; (B) design a two-stage optimization for the proposed predictive control which largely reduces the computational complexity since the WT model and WT load model remain nonlinear in this paper. This paper is organized as follows. Sections 2 and 3 present the wind energy conversion system model and WT load model, respectively. Section 4 introduces the overall control framework. Section 5 describes the predictive control strategy based on the two-stage optimization. Case studies are presented and discussed in Section 6, followed by conclusions. 2. Wind energy conversion system model In this paper, the WT is considered as a system which transforms wind energy into electrical energy. The high frequency part of wind speed is damped by the WT rotor due to its large inertia. Thus a lowpass filter is included in the model, whose transfer function is,

Gr ðsÞ ¼

1 1 þ sr s

ð1Þ

where sr is the time constant of the rotor. The WT output power can be represented as,

P¼

1 q  pR2  v 3  C p ðk; bÞ 2

xg ¼ k  xr

ð2Þ

xopt ¼ k

kopt v R

ð5Þ

The pitch actuator system is shown as,

b¼

1 b 1 þ ssb ref

bmin 6 b 6 bmax ; b_ min 6 b_ 6 b_ max

ð6Þ ð7Þ

where sb is the time constant of the pitch actuator system, b is the pitch angle, bref is the reference pitch angle, and ⁄min or ⁄max is the minimum or maximum limit of ⁄. 3. Wind turbine load model The WT load studied in this paper is the aerodynamic load on blades. There are three methods to model it: the BEM method, computational fluids dynamics method and vortex wake method [28]. With respect to the calculation time, the BEM is the most suitable method to be implemented online. Wind speed and aerodynamic characteristics along the rotor radius are various. Thus modeling blade load directly through aerodynamics is complicated. BEM imposes the idea of finite element analysis. The blade is divided into many blade elements with thickness dr, as Fig. 1 shows. The momentum theorem is applied to one blade element to obtain its load. The load of whole rotor is calculated by numerical integration. In this section, blade load modeling is briefly demonstrated and the characteristics of the blade load concerning generator speed and pitch angle is analyzed to illustrate how to reduce it. 3.1. Load modeling based on BEM The inducing factors are the most important parameters in the BEM method. The axial inducing factor denoted as a, describes the speed loss of air flow passing through the rotor. This part of kinetic energy loss is absorbed by WT. The tangential inducing factor denoted as b, describes the influences of air flow on the rotational rotor. For every blade element, they have to be calculated iteratively before calculating the load. The inducing factors are represented as,

8 a < 1a ¼ 2Bcpr  b 1þb

Cn 4F sin2 u

¼ 2Bcpr  4F sinCut cos u

ð8Þ

where B is the blade number, c is the chord of the blade element, r is the distance between the blade element and the rotor hub, Cn is the normal force coefficient, Ct is the tangential force coefficient, and F is the Prandtl loss factor. As Fig. 1 shows, u is the inflow angle represented as,

u ¼ arctan ð3Þ

ð4Þ

where xg is the generator speed, k is the gear box ratio, and xr is the rotor speed. The rotor speed can be controlled by controlling the generator speed. There is an optimal tip-speed ratio kopt with the maximal power coefficient Cpmax. Accordingly, the optimal generator speed is shown as,

:

where q is the air density, R is the rotor radius, v is the wind speed, k is the tip-speed ratio, and Cp is the power coefficient given as the following empirical formulas [27],


 116 12:5 C p ¼ 0:22  0:4b  5 e ki ki 1 1 0:035  ¼ ki k þ 0:08b b3 þ 1

The rotor shaft and generator shaft are both considered as a rigid body. The generator speed and the rotor speed are coupled by the drivetrain, whose model is simplified as,

Vð1  aÞ

xr rð1 þ bÞ

ð9Þ

where V is the upstream wind speed. If a > ac = 0.2, (8) needs to be replaced by the Glauert’s empirical formula:

137

W. Liu et al. / Electrical Power and Energy Systems 93 (2017) 135–145

V

V 1 a

Chord plane rr 1 b

dr

β

φ

V 1 a

V0

r

Blade element

Blade Fig. 1. BEM theory.

a¼

 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 2 þ Kð1  2ac Þ  ðKð1  2ac Þ þ 2Þ2 þ 4ðKa2c  1Þ 2

ð10Þ

K¼

8prF sin u BcC n

ð11Þ

2

Usually, set a = 0.3 and b = 0.2 as the initial value. The convergence condition is,

ja j  aj1 j < e and jb  b j

j1

j
ð12Þ

where j is the iteration number. Usually set e = 0.0001. The iteration steps are as follows: (i) Compute u using (9); (ii) Compute a and b using (8); (iii) Check the error between two iterations, if it satisfies the convergence condition, stop; otherwise, go to i). The relative inflow speed V0 is given by,

V0 ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 V 2 ð1  aÞ2 þ x2r r 2 ð1 þ bÞ

ð13Þ

The blade root bending moment is computed by integration,

Mx ¼

My ¼

1 2 1 2

Z

R

Rhub

Z

R

Rhub

qV 20 cC n rdr

ð14Þ

qV 20 cC t rdr

ð15Þ

where Mx is the in-the-plane bending moment, My is the out-ofplane bending moment, and Rhub is the hub radius. The blade load can be calculated from (8) to (14) and (15), with the knowledge of the look-up table for the chord c, the normal force coefficient Cn and the tangential force coefficient Ct. 3.2. Influence analysis The out-of-plane moment is much smaller than the in-theplane moment, and these two kinds of load share the same changing trend [26]. Therefore, only the in-the-plane moment is studied in this paper. The blade root bending moment is related to the upstream wind speed V, rotor speed xr and pitch angle b:

M ¼ f ðV; xr ; bÞ

ð16Þ

A WT with 35 m blades [29] is used to illustrate the influence of rotor speed xr and pitch angle b in Figs. 2–4. The whole blade is composed of NACA63-421, NACA63-418 and NACA63-415 airfoils. In these three figures, the WT load is shown as a surface of two control variables xr and b. As the color bar shows, red means heavy load and blue means light load. As the three figures show, with wind speed increasing, the red area expands and the blue area shrinks. Thus there is a positive correlation between the upstream wind speed V and blade root bending moment M. The black line in

the figures represents the rated power constraint. With wind speed above the rated wind speed, WTs deliver rated power, which means its working points have to be on the line. For each certain wind speed, blade root bending moment can be considered as a function defined on this line concerning variables rotor speed xr and pitch angle b. The optimal working point of WT should be on the line with the smallest blade load. The number of extremums and its location of the blade load are decided by the aerodynamics of the blades. For blade airfoils used in this paper, under a certain wind speed, the WT load has two extremums: X1, Y1 in Fig. 2, X2, Y2 in Fig. 3 and X3, Y3 in Fig. 4. Y1, Y2 and Y3 are around the working points of maximum power point tracking (MPPT). In Fig. 2, X1 is the global minimum and Y1 is the local minimum. However, in Figs. 3 and 4, Y2, Y3 is the global minimum and X2, X3 is the local minimum respectively. This indicates: (1) the optimization needs to avoid the local minimum and find the global minimum; (2) the optimal reference working point of WT probably suddenly changes with wind speed fluctuation. However, under real operation, WT’s working point cannot immediately jump to the new optimal working point. Thus the WT has to operate in a suboptimal region under transition. For different airfoils, there are similar questions of blade load. Thus before the online application, the number of extremums and which one is the global minimum should be calculated offline. 4. Control framework The main characteristic that distinguishes predictive control from other control methods is the receding horizon principle. In this section, the basic concept of predictive is presented, and the control framework is illustrated based on this concept. 4.1. Basic concept of predictive control Predictive control is a methodology based on certain basic principles [30]. All predictive control methods consist of three parts: prediction, online optimization and feedback correction. The future state or output of the system is acquired by the prediction model, which can describe the states and responses of the system. As Fig. 5 shows, given different optimization strategies, different outputs of the system will be obtained. In predictive control, the optimization problem only covers a finite horizon. The optimization yields an optimal control sequence, in which the first control is applied to the system [31]. This is the main difference between predictive control and traditional optimal control [32]. At each sampling instant, the optimization module solves an open-loop optimal problem with a finite horizon. As such, the computational complexity is significantly reduced. It is much easier for predictive control to handle complicate dynamic industrial processes with multiple variables and multiple constraints. As the optimization proceeds, the finite horizon optimal control problem is equivalent to the original problem

138

W. Liu et al. / Electrical Power and Energy Systems 93 (2017) 135–145 x10e5

1.2

Y1 Mx=7.3561×105 Nm

9 8

Rotor speed (p.u.)

1.1

7

1.0

6 5

0.9 X1 Mx=5.9055×105 Nm

4 3

0.8

2

0.7

0

4

8

16

12

20

24

28

32

36

Pitch angle (deg) Fig. 2. Blade root bending moment with wind speed 13 m/s.

x10e5

1.2

Y2 Mx=4.7456×105 Nm

9 8

Rotor speed (p.u.)

1.1

7

1.0

6 5

0.9

4

X2 Mx=5.5265×105 Nm

3

0.8

2

0.7

0

4

8

16

12

20

24

28

32

36

Pitch angle (deg) Fig. 3. Blade root bending moment with wind speed 14 m/s.

x105

1.2

9

Rotor speed (p.u.)

8

Y3 Mx=3.2493×105 Nm

1.1

7

1.0

6 5

0.9

4

X3 Mx=5.0049×105 Nm

0.8

3 2

0.7

0

4

8

12

16

20

24

28

32

36

Pitch angle (deg) Fig. 4. Blade root bending moment with wind speed 15 m/s.

or a modified one with the infinite horizon [32]. At each sampling instant, the feedback updates or corrects the prediction model before optimization to make it more accurate.

4.2. Predictive control framework Based on the concept explained previously, a predictive control strategy is proposed to deal with ramping events with load reduction. Its control framework is shown in Fig. 6. The input and output of every module and control signals are depicted in this figure. The

prediction model predicts the wind speed, which determines the optimization boundary. The online optimization has three objectives: minimizing the ramping rate, maximizing power generation and minimizing WT load. Since it is a multi-objective optimization problem, and the second objective is conflicting with the other two, it is handled as a two-stage optimization to make a trade-off and reduce the computational complexity. The objectives of the first stage are minimizing the ramping rate and maximizing power generation because they are both directly related to the WT output. The maximum available power is used as a constraint. The second stage is

139

W. Liu et al. / Electrical Power and Energy Systems 93 (2017) 135–145

past

plant is sufficient for overall system stability. This property confirms that for an open-loop stable plant, the closed-loop stability is obtained automatically for any stable controller. This is also true even with the nonlinear controller as long as it is input-output stable. In this paper, the model-plant mismatch is small enough to be ignored. With constraints of control variables, the online optimization is a bounded-input bounded-output controller. In this paper, only the transfer relation between pitch angle, rotor speed and power output is concerned. Therefore the WT can be considered as an open-loop stable plant. The proposed predictive control meets all the conditions of the DSC and its stability is guaranteed.

future Corresponding output of A Corresponding output of B

Output y Prediction horizon

Control strategy A

Input u

Control strategy B

k k+1

...

5. Predictive control strategy

t

k+P

rolling Fig. 5. Predictive control concept.

for load reduction and guarantees that the result from the first stage is satisfied. As such, decreasing WT load will not make WT output deviate from the power reference. However, if the power references are given by the WF, the online optimization will skip the first stage, as shown in Fig. 7. After each optimization, the WT output power reference Pref is obtained along with the generator speed and pitch angle references (xg_ref, bref). xg_ref can be tracked by torque control. The widely used vector control of DFIG comprises of two decoupled loops: the first is the rotor flux control loop and the second is the torque control loop which imposes electromagnetic torque to give PWM signals to rotor side converter in order to control generator speed. The vector control ensures a fast and accurate response. bref is tracked by the pitch actuator shown in Fig. 6. All the three objectives are influenced by wind speed. Apparently, the wind speed prediction cannot be absolutely accurate. Feedback is designed to compensate the wind speed error in order to avoid the control performance significantly deviating from the reference. The proposed predictive control has the internal model control (IMC) structure. The properties and stability of IMC was analyzed by Garcia and Morari [33]. The dual stability criterion (DSC) points out that when the model is exact, stability of both controller and

wind farm

The prediction model not only forecasts wind speed but also calculates WT available output power according to (2). Wind speed is an essential part of the prediction model since both the output power and the blade root bending moment are directly related to wind speed. With the development of wind speed prediction technologies and improved accuracy, single turbine prediction has become possible. The artificial neural network (ANN) and Markov chain (MC) were used to develop a new ANN-MC model for wind speed prediction in [34]. The short-term patterns of wind speed are captured by the ANN and the long-term patterns are obtained by the MC. Therefore, ultra-short-term wind speed can be predicted in different prediction horizons. The mean absolute

reference power

βref

online optimization

predicted wind speed

corrected wind speed

+

-

feedback

ωg_ref

βmax dβ/d /dttmmax max & d ax 1

5.1. Wind speed prediction

wind speed prediction

wind speed

pitch actuator

According to the control framework mentioned in the previous section, the detailed control strategy consisting of wind speed prediction, online optimization and feedback is developed in this section. The accuracy of wind speed prediction can affect the selection of the prediction horizon and the optimization boundary. The twostage optimization design is clarified with specific objectives, constraints and control variables. The optimization process including working point transition and hysteresis characteristic is illustrated in Section 5.3 with input and output of every optimization stage. Time series prediction is imposed in feedback to fit the nonlinear mapping between output power error and wind speed error.

1 s

power reference

Pref

·/·

+

_

Tg_ref

βmin & ddβ/d /dttmmin in

torque control

β

PWM signal

measured output power

converter ωr

gear box

ωg DFIG

rotor Fig. 6. Control framework of proposed control strategy.

grid

140

W. Liu et al. / Electrical Power and Energy Systems 93 (2017) 135–145

Corrected wind speed

Wind farm

Stage 1

Pref

Pref

(ωg_ref, βref)

Stage 2

Wind turbine

Pref Online optimization Fig. 7. Power reference in the online optimization.

percentage error (MAPE) of this method is less than 7.5% in 7.5 s. In this paper, the prediction horizon Tp = 5 s. 5.2. Online optimization The objective functions comprise three important terms, which are the ramping rate, power generation and WT load. The specific objectives are given as follows.

The key step of converting multiple objectives into a single objective is the determination of weight coefficients. f1 and f2 are converted into per-unit values with the base value Prate2 and Prate respectively to determine weight coefficients. Weight coefficients show the preference of a decision maker. The transformed optimization model is,

(

min f ¼

(a) Ramping rate

min f 1 ¼

 Np
X DP i 2 i¼1

Dt

ð17Þ

ð18Þ

(b) Power generation

min f 2 ¼

Np X ðPi Þ  Dt

ð19Þ

i¼1

where Pi is the WT power output of the ith time interval. N 1 X q V 20 cC n ðRhub þ nDrÞ  Dr 2 n¼1

ð20Þ

where f3 is the discrete form of (14), N is shown in (21) and Dr is the step of numerical integration, usually Dr = 0.1 m.


 R  Rhub N ¼ int Dr

Np1 X

P2rate

ð21Þ

The online optimization is a multi-objective optimization problem. The generator speed reference xg_ref and pitch angle reference bref are control variables. Usually, the three objectives are combined into one objective. However, the optimization problem is too complicated to be calculated online. In order to simplify the calculation, the whole optimization process is divided into two stages: the power output optimization stage and the load optimization stage. The power output optimization stage is to minimize the ramping rate and maximize power generation. Pi is the control variables and the optimal output power can be obtained at this stage. The load optimization stage is to minimize load by setting xg_ref and bref. The power reference is an equality constraint so that the load optimization will not affect the WT power output. Stage1: power output optimization The optimization problem of the first stage is a multi-objective programming problem. Multiple objectives are merged and then solved by a single objective programming algorithm.

ðPi  Pi1 Þ2 þ P22 þ P2Np1  2P1 P 2

i¼3



) Np k e  Dt X ðPi Þ þ Prate i¼2

ð22Þ

s.t.

8 > < P1 ¼ P0 0 < Pi < Pi pred > : PNp ¼ PNp pred

i ¼ 2; 3;    ; Np  1

ð23Þ

where P0 is obtained from the previous optimization and Pi_pred is the predicted maximum power of the ith time interval. Due to the receding horizon principle of predictive control, PNp will not affect the control performance. P1 and PNp are given fixed values and removed from the objective function. As such, the objective function can be transformed into a standard quadratic programming form as,

min

(c) WT load

min f 3 ¼

Dt  2

 2PNp PNp1

where DPi is the power change in time interval Dt, and Np is shown in (18). int() is the bracket function.


 Tp Np ¼ int Dt

"

ks


 1 T x Hx þ CT x 2

ð24Þ

where x is a vector consisting of Pi. Matrix H, known as the Hessian matrix, and C are given as follows. H is a (Np-2)-dimensional square matrix and C is a (Np-2)-dimensional vector.

0

2

B 1 B B H¼B B B @

1

1 2

1 .. . 1

2 1

T

C ¼

C C C C  2ks C Dt 2  P 2 C rate 1 A

ð25Þ

2

ke  Dt ke  Dt ke  Dt þ ; ;; N  Prate P2rate  Dt 2 N  Prate N  Prate 2P 1  ks

! ð26Þ

This quadratic programming problem is solved by the function quadprog in the MATLAB optimization toolbox. The Trust-regionreflective algorithm is chosen to compute the quadratic programming due to its shortest convergence time among three available algorithms of quadprog. Stage2: WT load optimization The optimization objective is to minimize load by setting xg_ref and pitch angle bref. This optimization model can be represented as,

min M ¼

N 1 X q V 20 cC n ðRhub þ nDrÞ  Dr 2 n¼1

ð27Þ

W. Liu et al. / Electrical Power and Energy Systems 93 (2017) 135–145

erence, keep xref,i = xopt and change the pitch angle reference according to (2).

s.t.

8 xmin < xg ref ;i < xmax > > > > > > < 0 < jxg ref ;i  xg ref ;i1 j < Dxlim bmin < bref ;i < bmax > > > 0 < jbref ;i  bref ;i1 j < Dblim > > > : pi ðv i ; xg ref ;i ; bref ;i Þ  Pref ;i ¼ 0

5.4. Feedback

i ¼ 2; 3;    ; Np

ð28Þ

where Pref,i is the power reference from the first stage optimization, and the equality constraint guarantees that the power reference is tracked during the second stage optimization. Limiting the generator speed and the pitch rate can prevent load from fluctuating too much, which ensures the fatigue of WTs does not increase. This nonlinear programming problem is solved by the function fmincon in the MATLAB optimization toolbox. There are two suitable algorithms of fmincon for this problem: interior point (IP) and sequential quadratic programming (SQP). The calculation time of the IP is around 0.33 s per set-point on average, while the SQP’s is 0.24 s. Therefore, the SQP is chosen to compute (27).

It is necessary to design feedback to compensate the prediction error because the prediction model cannot be absolutely accurate. The WT output power error at the sampling instant t + 1 is shown as,

eðt þ 1Þ ¼ pðt þ 1Þ  p1 ðt þ 1jtÞ

The Stage2 optimization problem is a non-convex optimization. According to the analysis in Section 3.2, the optimal operating point of WTs may suddenly change with various wind speed. In this paper, the critical point of wind speed is 13.75 m/s. When V  13.75 m/s, the optimal working point is X (for example, X1 in Fig. 2); when V > 13.75 m/s, the optimal working point is Y (for example, Y2, Y3 in Figs. 3 and 4). In order to ensure the optimization converge to the global minimum, the feasible region should be constrained in the neighborhood of X or Y respectively by changing the limits of variables in (28). To avoid frequent change between the distant two working points, the hysteresis characteristic is high introduced, as Fig. 8 shows. vlow thre and vthre is the low and high threshold value of hysteresis characteristics, respectively. vc is the critical wind speed point. When the optimal working point is changed, WTs will reach the new optimal working point after a brief transition. During the transition, the set-points given by the proposed control strategy should satisfy the equality constraint in (28), which means the power reference is still tracked. Each set-point (xg_ref,i, bref,i) is closer to the optimal operation point (xg_opt, bopt) than the previous set-point (xref,i-1, bref,i-1). For example, in Fig. 3, the power reference constraint is described as the black line. If wind speed varies from 13 m/s to 14 m/s, the working point of WT will move on the line from X towards Y. The entire optimization process is shown in Fig. 9. Since Stage2 is a quite complicate optimization problem with multiple constraints, sometimes there may be no feasible solutions found by fmincon. This situation will lead to output power deviating from the power reference and control performance deteriorating. In this situation, the power reference from Stage1 or WF is used to determine the generator speed and pitch angle. The generator speed is firstly calculated while the pitch angle reference remains unchanged. If the generator speed reaches the optimal value xopt and the output power is still smaller than the power ref-

low vthre

Y transition

X high vthre

Fig. 8. Hysteresis characteristic.

ð29Þ

where p(t + 1) is the measured WT output power, and p1(t + 1|t) is the predicted WT output power of sampling instant t. e(t + 1) contains the uncertainty not considered in the prediction model. It can be utilized to predict the future output error as a supplement. It is hard to identify the exact reason causing the error. Therefore, a heuristic approach is usually used. The corrected wind speed prediction can be represented as,

vcor ðt þ 2Þ ¼ v1 ðt þ 2jt þ 1Þ þ e1 ðt þ 1Þ 5.3. Working point transition

vc

141

ð30Þ

where v1(t + 2|t + 1) is the wind speed prediction of the sampling instant t + 1. e1(t + 1) is usually derived from e1(t + 1) = he(t + 1), where h is the correction vector comprised of weight coefficients. The input of feedback is power output error and the output of feedback is wind speed error. There is a highly nonlinear relation between the input and output. In this paper, feedback is managed as a time series prediction. The neural network is chosen as the time series model because it can fit highly nonlinear mapping. The training is equivalent to the determination of h. The neural network model is established by the neural network time series tool of MATLAB. Feedback is deemed as a nonlinear autoregressive problem with an external input (NARX). The architecture of this neural network is shown in Fig. 10. The training is done offline by the Levenberg-Marquardt algorithm. The sensitivity analysis with different numbers of delays and neurons shows that the NARX neural network with 6 delays and 20 neurons has the least mean squared error (MSE). The detailed results of the sensitivity analysis are shown in Section 6.

6. Case study In this section, a 1.5 MW WT [29] is used to demonstrate the performance of the proposed control. The wind speed data are extracted from 14:21 to 15:00 on December 24, 2011 in the Xiangfeng wind farm of the State Grid Jibei Electric Power Company Limited of China. The data resolution is 1 s. The case study is divided in three parts. Firstly, the performances of the feedback module are illustrated. The last two parts are two scenarios to test the efficacy of the proposed control corresponding to the power reference selection in Fig. 7.

6.1. Feedback performance Wind speed data from 14:21 to 14:50 are used for training and data from 14:50 to 15:00 are used for simulation. In this paper, the training dataset is divided into three subsets: 70% for training, 15% for validation and 15% for test. Table 1 shows the mean MSE with 50 training vectors for different architectures of the NARX neural network. The least MSE is achieved with 6 delays and 20 neurons and the architecture of proposed neural network is determined. Predicted wind speed, measured wind speed and corrected wind speed are shown in Fig. 11. Corrected wind speed follows measured wind speed more closely. The MAPE of the corrected

142

W. Liu et al. / Electrical Power and Energy Systems 93 (2017) 135–145

Start

Input wind speed and blade airfoil aerodynamic data

yes

no

Wind farm power command ?

Power output optimization

Use wind farm power command as reference power

Use optimal power output as reference power

Stage1 Stage2 yes

no

high ? thre

v

Wind speed

yes

Change initial value to ωg0,i=ωmppt,i β0,i=βmppt,i

Wind speed

low ? vthre

no Change feasible region to a neighborhood of Y

Set initial value ωg0,i=ωg0,i-1, β0,i=β0,i-1

Change initial value to ωg0,i=0.7 p.u., β0,i=0°

Set feasible region as previous optimization

Change feasible region to a neighborhood of X

Solve optimization problem (27)

yes

no Under transition? no

No feasible solution found?

Determine set-point (ωg_ref,i, βref,i) that satisfies equality constraint in (28)

yes ωg_ref,i=ωopt Calculate βref,i according to (2) Output results

End Fig. 9. Optimization flowchart.

Table 1 Sensitivity analysis for different numbers of delays and neurons.

neurons

time delays e(k+1)

1:N

activation function

Delays

activation function e1(k+1)

1:N

Neurons

5

6

7

16 20 24

0.0408 0.0637 0.0724

0.0359 0.0270 0.0348

0.0707 0.0886 0.1146

6.2. Control performance compared with MPPT b-bias vector

b-bias vector

Fig. 10. NARX neural network.

wind speed is 0.9% compared with 7.5% of the predicted wind speed, which is significantly reduced. The feedback module makes wind speed prediction more accurate.

Weight coefficients in (22) are set as ke = 1 and ks = 1000. The huge difference between these two coefficients is because the value of objective (b) is much bigger than that of objective (a) with the per unit value form. The average calculation time of one setpoint is around 0.24 s, which is fast enough for the online application.

143

W. Liu et al. / Electrical Power and Energy Systems 93 (2017) 135–145

wind speed (m/s)

predicted wind speed measured wind speed corrected wind speed 16 14 12 10 0

50

100

150

200

250

300

350

400

450

500

time (s) Fig. 11. Wind speed.

WT output power (kW)

power reference MPPT predictive control 1600 1400 1200 1000 800 600 0

50

100

150

200

250

300

350

400

450

500

time (s) Fig. 12. Power generation.

x 105

MPPT predictive control

load Mx (N*M)

9 8 7 6 5 4 3 0

50

100

150

200

250

300

350

400

450

500

time (s) Fig. 13. WT load optimization.

The power reference shown in Fig. 12 is from the first optimization stage. The ramping rate of the predictive control is reduced and the energy production is 190.7549 kWh. The energy production of the MPPT control is 197.6927 kWh. In the derated section, the ramping rate is reduced at the cost of approximately 7.56% wind energy reduction. The predictive control (the solid line in Fig. 12) tracks the optimal power curve closely, which also validates the necessity of the feedback module. The WT load optimization result is shown in Fig. 13. The hysteresis characteristic can effectively prevent the WT working point from changing frequently. The WT working point is not changed until simulation time 350 s. During simulation time 1–349 s and 433–500 s, WT operates at working point X. During simulation time 362–427 s, WT operates at working point Y, which is close to the MPPT working point. During simulation time 350–361 s and 427–432 s, the WT working point is under transition. As Fig. 12 shows, the power reference is still tracked during the transition process. When the WT operates at the optimal working point

Y, the WT load is 0.6% smaller than the MPPT’s. When the WT operates at the optimal working point X, the WT load is 12.62% smaller than MPPT’s. Fig. 14 shows the power spectral density (PSD) of the blade root bending moment Mx. The PSD of the WT load can show its fatigue damage to the mechanical structure. The blade root bending moment is slightly deceased between 0.1 Hz and 0.5 Hz after the optimization. Even though the fatigue reduction is not included as an objective in the online optimization module, the rate change constraints of the generator speed and pitch angle prevent load vibration from increasing. Both generator speed and pitch angle references are optimized to achieve the control objectives. The control performances of generator speed and pitch angle are shown in Fig. 15 and 16. During simulation time 1–349 s, the optimized generator speed deviates from the upper limit, and the pitch angle of the predictive control is much smaller than the MPPT control’s, which correspond to the optimal operating point X of the WT. With the MPPT control, the generator speed reaches the upper limit and the pitch angle is rel-

144

W. Liu et al. / Electrical Power and Energy Systems 93 (2017) 135–145

PSD of Mx (Nm 2 /Hz)

1014

MPPT predictive control

1010

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

frequency (Hz) Fig. 14. Power spectral density of in-the-plane moment.

generator speed (p.u.)

MPPT predictive control 1.2 1.1 1 0.9 0.8 0.7 0

50

100

150

200

250

300

350

400

450

500

time (s) Fig. 15. Generator speed control performance.

pitch angle (deg)

MPPT predictive control 16 12 8 4

0

50

100

150

200

250

300

350

400

450

500

time (s) Fig. 16. Pitch angle control performance.

wind turbine output power (kW)

atively higher, corresponding to the local minimum working point Y. During 362–427 s, the WT operates at the optimal working point Y, the WT load is still slightly smaller than MPPT’s. Heavy load is one of the important reasons to cut out WTs to avoid mechanical damages. With significant load reduction, large-scale ramp-down events caused by WTs being cut out simultaneously are less likely to happen.

6.3. Control performance compared with PRL When a ramp-up event happens, the WF often needs to limit power ramping rate by giving power references to WTs. The PRL is a commonly used ramping control method. The power reference is shown in Fig. 17. The ramping rate is limited to 6 kW/s. The proposed predictive control can regulate WTs to track power reference

power reference predictive control

1600 1400 1200 1000 800 600 0

50

100

150

200

250

300

time (s) Fig. 17. Power reference tracking.

350

400

450

500

145

W. Liu et al. / Electrical Power and Energy Systems 93 (2017) 135–145

9

x 105 PRL predictive control

load Mx (N*M)

8 7 6 5 4 0

50

100

150

200

250

300

350

400

450

500

time (s) Fig. 18. WT load optimization with power reference.

accurately. Fig. 18 shows the load optimization results with the power reference from the WF. Compared with the PRL, the WT load of the proposed control is reduced by 24.5%. 7. Conclusions This paper proposes a predictive ramping control strategy for WTs with load reduction. The proposed predictive control minimizes the ramping rate, maximizes power generation and reduces WT load simultaneously. The ramping rate and power generation are both taken into account to determine the power reference by setting proper weight coefficients. Load reduction is realized on the premise of tracking the power reference. With the two-stage optimization, not only is the computational complexity of the online optimization reduced, but also it can find a balance among the three objectives. When the optimal working point changes, the proposed control can find the new one and control the WT to operate at the set-point. Due to the significant load reduction, the risk of large-scale ramp-down events due to WTs being cut out can be reduced with the proposed predictive control. Acknowledgements The authors would like to thank the financial support by National Basic Research Program of China (973 Program) under grant #2012CB215101. References [1] Senjyu T, Sakamoto R, Urasaki N, Funabashi T, Fujita H, Sekine H. Output power leveling of wind turbine generator for all operating regions by pitch angle control. IEEE Trans Energy Convers 2006;21(2):467–75. [2] Wang H, Liu X, Wang C. Probabilistic production simulation of a power system with wind power penetration based on improved UGF techniques. J Mod Power Syst Clean Energy 2013;1(2):186–94. [3] Liu Y, Fan R, Terzija V. Power system restoration: a literature review from 2006 to 2016. J Mod Power Syst Clean Energy 2016;4(3):332–41. [4] Zhang H, Li C, Liu Y. Quantitative frequency security assessment method considering cumulative effect and its applications in frequency control. Int J Electr Power Energy Syst 2015;65:12–20. [5] Cui M, Ke D, Sun Y, Gan D, Zhang J, Hodge B. Wind power ramp event forecasting using a stochastic scenario generation method. IEEE Trans Sustain Energy 2015;6(2):422–33. [6] Moiseeva E, Hesamzadeh M, Biggar DR. Exercise of market power on ramp rate in wind-integrated power systems. IEEE Trans Power Syst 2015;30 (3):1614–23. [7] Wang S, Yu D, Yu J. A coordinated dispatching strategy for wind power rapid ramp events in power systems with high wind power penetration. Int J Electr Power Energy Syst 2015;64:986–95. [8] Ela E, Kirby B. ERCOT event on February 26, 2008: lessons learned. National Renewable Energy Laboratory; 2008 [NREL/TP-500-43373]. [9] Qi Y, Liu Y. Wind power ramping control using competitive game. IEEE Trans Sustain Energy 2016;7(4):1516–24.

[10] Xu L, Wang Y. Dynamic modeling and control of DFIG-based wind turbines under unbalanced network conditions. IEEE Trans Power Syst 2007;22 (1):314–23. [11] Zhi D, Xu L. Direct power control of DFIG with constant switching frequency and improved transient performance. IEEE Trans Energy Convers 2007;22 (1):110–8. [12] Huang H, Chung CY. Coordinated damping control design for DFIG-based wind generation considering power output variation. IEEE Trans Power Syst 2012;27(4):1916–25. [13] Kasem AH, El-Saadany EF, El-Tamaly HH, Wahab MAA. Power ramp rate control and flicker mitigation for directly grid connected wind turbines. IET Renew Power Gener 2010;4(3):261–71. [14] Ackermann T. Wind power in power systems. New York (USA): John Wiley & Sons, Inc.; 2005. [15] Bossanyi EA. Un-freezing the turbulence: application to LiDAR-assisted wind turbine control. IET Renew Power Gener 2013;7(4):321–9. [16] Bossanyi EA. Wind turbine control for load reduction. Wind Energy 2003;6 (3):229–44. [17] Bossanyi EA, Kumar A, Hugues-Salas O. Wind turbine control applications of turbine-mounted LIDAR. J Phys Conf Ser 2014;555(1):329–70. [18] Saravanakumar R, Jena D. Validation of an integral sliding mode control for optimal control of a three blade variable speed variable pitch wind turbine. Int J Electr Power Energy Syst 2015;69:421–9. [19] Xiao S, Geng H, Yang G. Non-linear pitch control of wind turbines for tower load reduction. IET Renew Power Gener 2014;8(7):786–94. [20] Schuler S, Schlipf D, Cheng PW, Allgower F. ‘1-optimal control of large wind turbines. IEEE Trans Control Syst Technol 2013;21(4):1079–89. [21] Hand MM. Variable-speed wind turbine controller systematic design methodology: a comparison of non-linear and linear model-based designs. National Renewable Energy Laboratory; 1999 [NREL/TP-500-25540]. [22] Chekkal S, Lahaçani NA, Aouzellag D, Ghedamsi K. Fuzzy logic control strategy of wind generator based on the dual-stator induction generator. Int J Electr Power Energy Syst 2014;59:166–75. [23] Zhao H, Wu Q, Rasmussen CA, Blanke M. Adaptive speed control of a small wind energy conversion system for maximum power point tracking. IEEE Trans Energy Convers 2014;29(3):576–84. [24] Dang DQ, Wang Y, Cai W. Offset-free predictive control for variable speed wind turbines. IEEE Trans Sustain Energy 2013;4(1):2–10. [25] Soliman M, Malik OP, Westwick DT. Multiple model predictive control for wind turbines with doubly fed induction generators. IEEE Trans Sustain Energy 2011;2(3):215–25. [26] Liu W, Cui Y, Liu Y. Predictive control of wind power ramping based on wind turbine load optimization. In: IET renewable power generation conference. p. 1–4. [27] Han B, Zhou L, Yang F, Xiang Z. Individual pitch controller based on fuzzy logic control for wind turbine load mitigation. IET Renew Power Gener 2016;10 (5):687–93. [28] Leishman JG. Challenges in modelling the unsteady aerodynamics of wind turbines. Wind Energy 2000;5(2):85–132. [29] Parameters of a 1.5 MW wind turbine, Baidu Wenku; 2016. . [30] Camacho EF, Alba CB. Model predictive control. 2nd ed. Berlin (German): Springer Science and Business Media; 2013. [31] Ye H, Liu Y. Design of model predictive controllers for adaptive damping of inter-area oscillations. Int J Electr Power Energy Syst 2013;45(1):509–18. [32] Mayne DQ, Rawlings JB, Rao CV, Scokaert POM. Constrained model predictive control: stability and optimality. Automatica 2000;36(6):789–814. [33] Garcia CE, Morari M. Internal model control.1. A unifying review and some new results. I&EC Process Des Dev 1982;21(2):308–23. [34] Kani SAP, Ardehali MM. Very short-term wind speed prediction: a new artificial neural network–Markov chain model. Energy Convers Manage 2011;52(1):738–45.