This paper addresses linear quadratic regulation (LQR) for variable speed variable pitch wind turbines. Because of the inherent nonlinearity of wind turbines, a set of operating conditions is identified and then a LQR controller is designed for each of the operating points. The feedback controller gains are then interpolated linearly to get a control law for the entire operating region. In addition, the aerodynamic torque and effective wind speed are estimated online to get the gainscheduling variable for implementing the controller. The potential of this method is verified through simulation with the help of MATLAB/Simulink and GH Bladed. The performance and mechanical load when using LQR are also compared with those obtained when using a PI controller.
I. INTRODUCTION
Recently, variable speed variable pitch wind turbines (VSVP WT) have been continuously increasing their market share. This configuration is the best for multimegawatt machines because it can maximize the energy captured over a wide range of wind speeds and reduce the mechanical fatigue by using aerodynamics control systems. In this configuration, the turbine operates with a variable speed and a fixed pitch if the wind speed is in below the rated region (region 2) to achieve maximum aerodynamic efficiency. When the wind speed is above the rated region (region 3) the rotor is regulated at its rated speed by varying the pitch angle to ensure that the mechanical limitations are not exceeded. The controller should be designed intelligently to transit smoothly between the two regions and to ensure the other requirements during transition such as limiting bladedtip noise, minimizing output power fluctuations, etc.
The thrust force acting on the rotor and torque developed by the wind turbine are nonlinear functions of the wind speed, rotor speed and pitch angle. In addition, modern, largesize wind turbines are usually equipped with individual pitch actuators at each blade and force/moment sensors or accelerometers on the tower, nacelle and blades. These inputs and outputs combine with the structural modes to make the machine an inherently nonlinear multiinput multioutput (MIMO) system. However, most of the control methods applied to these multicontrolobjective MIMO systems are implemented by using multiple singleinput singleoutput (SISO) loops
[1]
. The basic multipleSISOloop structure of a wind turbine system is illustrated in
Fig. 1
. This configuration has two loops which operate independently of each other. The top part of this figure is the pitch PI control loop which has a major role for regulating the rotor speed in the above rated region. Below the rated wind speed, the blade pitch angle is pitched off at an optimal value so as to maintain maxCp operation. In addition, in this region, varying the rotor speed proportionally to the wind speed is a function of the torque control loop presented in the bottom part of the figure.
Control loop of a MW Wind turbine.
The linear quadratic regulation (LQR) technique has been applied to wind turbines in
[2]

[4]
. Most of these studies focus on the high wind speed region, using only the drivetrain model or considering only the tower vibration mode in the foreaft direction. The purpose of this work is to design a MIMO LQR controller for a multimegawatt wind turbine.
The proposed controller has the same strategy for the entire operating region. The controller is synthesized with the objective of reaching a tradeoff between maximizing the energy captured from the wind and mitigating the mechanical load of tower in both the foreaft and sidetoside directions. Because of the nonlinearity of wind turbines, the controller is designed for specific operating points. The feedback gains are then interpolated or extrapolated for the whole operating region. The potential of the LQR controller is verified by the commercial wind turbine simulation package GH Bladed, and compared with a PI controller.
This paper is presented as follows. Section II presents the nonlinear dynamics of the wind turbine. Section III derives a linear model for control design and to develop the linear quadratic optimal control which aims for a tradeoff between the control objectives. Finally, in section IV, the proposed controller is illustrated in a highfidelity simulation environment on a representative 2MW wind turbine.
II. WIND TURBINE MODEL
A block diagram for VSVP WT.
Fig. 2
shows a basic block diagram for the entire variable speed variable pitch wind energy conversion system which can be structured as several interconnected subsystems. In that figure,
F_{a}
is the thrust force,
β
is the blade pitch angle, and
T_{a}
and
T_{g}
are the aerodynamic torque and the reaction generator torque, respectively.
Ω_{r}
and
Ω_{g}
are the rotor speed and the generator speed,
v
is the wind speed which is described by the mean wind speed
v_{m}
, and turbulent wind
v_{t}
, is given by (1).
The wind speed
v_{r}
seen by the rotor plane is defined by (2), where
n_{fa}
is the nacelle displacement in the foreaft direction.
 A. Aerodynamic Conversion
The thrust force acting on the entire rotor, the useful torque developed by the wind turbine and the aerodynamic power captured from the wind are expressed by the highly nonlinear equations as follows
[5]
:
where
p
is the air density,
R
is the rotor radius,
C_{T}
is the force coefficient, and
C_{P}
is called the power coefficient which represents the wind turbine power conversion efficiency.
C_{T}
and
C_{P}
are functions of the blade pitch angle and the tipspeed ratio,
λ
, defined by:
The coefficients
C_{T}
and
C_{P}
are very important in terms of the turbine control system design. The characteristics for different values of the tipspeed ratio and the pitch angle are illustrated in
Fig. 3
and
Fig. 4
.
Fig. 3
indicates that there is one specific set of
λ
and
β
where the turbine power coefficient is at its maximum, C
_{P}
max. This means that if the pitch angle is fixed at the optimal value
β_{0}
and the rotor speed is varied proportionally to wind speed to keep tipspeed ratio at
λ_{0}
, where
C_{P}
coefficient is at its maximum, the wind turbine will extract the maximum power from the wind. However, this operation condition occurs only in the low wind speed region. When the wind speed is high, the tipspeed ratio and the pitch angle are not at their optimal values to ensure that the power and rotor speed are being regulated around their rated values. This is done in order to avoid excessive electrical and mechanical stresses.
Power coefficient CP.
Force coefficient CT.
 B. Drivetrain Dynamics
The drivetrain system can be modeled as two inertias interconnected by a springdamper and a gearbox, as schematically represented in
Fig. 5
[7]
. The rotor of the wind turbine is mounted at the low speed side of the drivetrain system. On the high speed side, the generator is mounted giving the opportunity to control the reaction torque from the generator. The dynamic model of the drivetrain system is modeled in the form of equations (7) and (8). In these equations,
Ω_{r}
and
Ω_{g}
are the rotor speed and generator speed, respectively,
J_{r}
and
J_{g}
are the rotor inertia and generator inertia, respectively,
B_{r}
and
B_{g}
are the damping of the low speed shaft and the high speed shaft, respectively, and
k_{s}
and
c_{s}
are the torsional stiffness and torsional damping of the drivetrain axis.
A drivetrain model.
The gearbox ratio
N
is defined as
N
=
Ω_{g}
/
Ω_{r}
.
 C. Tower Dynamics
In order to simply get a sufficient dynamic model of the tower for designing the controller in this paper, the first bending mode in both the foreaft direction and sidetoside direction are considered. When neglecting the tower top rotation, torsion deformation, yawing effects, higher bending modes and tilt influences, the tower top displacements of the first bending mode in both directions can be modeled by two common massspringdamper systems as in (9) and (10)
[5]
. It is also noted that the effects of the sideward aerodynamic fore on the sideward vibration of the tower can be neglected because it is much smaller than the reaction torque of the nacelle.
M_{T}
is the tower top effective mass, which is determined by the sum of the nacelle mass including the rotor and the top equivalent tower mass.
C_{T}
and
K_{T}
are the structural damping and bending stiffness, respectively.
are the tower top displacement, velocity and acceleration in the foreaft direction, respectively.
represent these values in the sidetoside direction.
h
is the hub height. The reaction torque of the nacelle,
T_{nac}
, is given by:
 D. Generator Dynamics
The generator torque is controlled by a power converter which also manages the active and reactive powers of the generator. Because the high speed switching power electronics are able to set the electric generator torque almost instantaneously with respect to the mechanical dynamics, the generator dynamics are sufficiently modeled by the first order transfer function presented in (12):
In the above equation,
T_{g}(s)
and
T_{gCMD}(s)
are the generator torque and generator torque commands, respectively.
τ_{g}
is the time constant of the generator dynamics.
 E. Pitch Actuator
When wind turbine systems operates in the high wind speed region, a high rotor speed that may cause mechanical damage which can no longer be managed by increasing the generated power because this would lead to overloading the generator and converter
[4]
. As seen in
Fig. 3
,
C_{P}
can be reduced by varying the pitch angle
β
to maintain the output power at its rated power and to regulate the rotor speed at its rated value, instead of changing the rotor speed proportionally to the wind speed as in the low wind speed region.
A second order pitch actuator model.
The pitch actuator is a nonlinear servo that generally rotates all of the blades or part of them in unison. In a closed loop the pitch actuator can be modeled as a linear loworder dynamic system with saturation in the amplitude and derivative of the output signal.
Fig. 6
shows a block diagram of a secondorder actuator model. The dynamic behavior of a pitch actuator operating in its linear region is described by the following transfer function:
where
β
and
β_{CMD}
are the actual and desired pitch angles, respectively.
ω_{β}
and
ξ_{β}
are the natural frequency and damping ratio of the pitch actuator dynamics, respectively. For a multiMW wind turbine, the pitch angle ranges from 3
^{o}
to 90
^{o}
and varies at a maximum rate of ±10
^{o}
/s.
III. CONTROLLER DESIGN
 A. Model for the Controller Design
Equations (7) through (13), describing the dynamics of a wind turbine, can be written in compact form as:
By defining the state, input and output vectors:
Equation (14) is a highly nonlinear model, due to the expression of the extracted thrust force
F_{a}
and the aerodynamic torque
T_{a}
as in (3) and (4). In order to design the controller, the global model can be linearized around the operating points by linearizing the aerodynamic torque and thrust force. The deviations of the aerodynamic torque and thrust force from the steady state values are expressed in (14).
where the operator δ corresponds to the deviation of the values from the linearization point OP(x0, u0, v0), and the coefficients in the above equations are defined by:
Thus, the wind parameterized linear model of the wind turbine around an operating point can be set on the state space representation as in (18):
where
δx
=
x

x
_{0}
;
δu
=
u

u
_{0}
;
δy
=
y

y
_{0}
.
The state space matrices A(v0), B(v0) and B
_{v}
(v0) depend on the operating points.
 B. Target Trajectory
As mention above, when the wind speed is below the rated region, the blade pitch angle is constantly maintained at its optimal value β
_{0}
. The rotor speed is changed proportional to the wind speed so as to maintain the tip speed ratio at a constant value λ
_{0}
. On the other hand, above the rated region, the pitch angle is varied in order to regulate the rotor speed and to generator torque/power at the rated values.
Fig. 7
illustrates the goal trajectories for rotor speed (the first window), blade pitch angle (the second window), generator torque (the third window) and electrical power (the fourth window).
Schedule for regulation set point.
 C. LQR Design
With the operating points determined, a set of controllers can be synthesized by applying a LQR for the linear model presented in (18), with a quadratic cost function for the regulation problem at an operating point defined as:
where Q(
v_{0}
) and R(
v_{0}
) are nonnegative and symmetric matrices of the weights. The control law optimizing the above criterion
J
, is a state feedback law with an optimal gain matrix determined by solving an LQR problem. In the implementation, the intermediate controllers are then interpolated linearly from the discrete number of controllers, which have been designed for specific operating points.
Fig. 8
illustrates the structure of the wind turbine system. The LQR feedback gains are scheduled by the effective wind speed, and estimated by a wind speed estimator. This estimation will be presented in the following subsection. The drivetrain vibration is damped by a drivetrain damper.
LQR control system.
IV. WIND SPEED ESTIMATION
As explained above, the target trajectories are scheduled by the wind speed. However, there is a very rough measurement of the wind speed available. Therefore, it is necessary to have a good wind speed estimator to get the controllers implementable.
Fig. 9
shows a schematic diagram of a wind speed estimator, which consists of three consecutive processing modules. The first module is the aerodynamic torque estimation, following by a 3D lookup table to calculate the effective wind speed, and the next is a low pass filter
[6]
.
The governing equations of the drivetrain model (7) and (8) can be combined into one equation as:
Wind speed estimator.
where
J_{t}
=
J_{r}
+
N^{2}J_{g}
, and
T_{L}
represent all of the mechanical losses.
Then, augmenting (21) with the unknown aerodynamic torque Ta, the above equation can be expressed in the state space form of (22)
[6]
,
[7]
:
where
w_{g}
and
w_{v}
represent the input process and the output sensor noise. A Kalman filter is applied to estimate the aerodynamic torque, which has the structure of (23) where the variables with a hat are to be estimated and
L
is the Kalman filter gain.
The second stage of the wind speed estimation is to calculate wind speed through a 3D lookup table, which has three inputs. These inputs are the estimated aerodynamic torque, measured rotor speed and measured pitch angle. The 3D lookup table is build by using equation (4), which is reformulated in the form of:
where
is the estimated tipspeed ratio.
Fig. 10
illustrates the estimated aerodynamic torque and the estimated wind speed with a mean value of 8m/s and a turbulence intensity (TI) of 23.2%. The real signals and the estimated values are presented in separate windows. As shown in the fourth window of
Fig. 10
, the estimated wind speed looks like a spatially average wind speed, which is different than the wind speed measured at the hub (dashdot line). The high frequency components in the turbulent wind are filtered out to get the estimated wind speed as a good scheduling input variable.
Estimated aerodynamic torque and wind speed.
V. NUMERICAL SIMULATION
The proposed controller is validated through simulations with stochastic wind input in the entire operating region. The performance and mechanical load of proposed method are compared with those of a PI controller. The control law was developed in MATLAB/Simulink and compiled into a Dynamic Link Library (DLL). The external controller in the DLL format will be used with the GH Bladed package to perform a full system simulation of the wind turbine.
A typical example for below the rated region is performed for a wind speed of 8m/s, perturbed by a turbulence intensity of 23.2%. The trace of hub height wind speed is shown in the top window of
Fig. 11
.
Fig. 11
also shows the rotor speed behavior, pitch response, generator torque, electric power, tower sidetoside moment (Mx
_{T}
) and foreaft moment (My
_{T}
) versus time. The solid lines in
Fig. 11
represent the responses of the wind turbine when using the LQR controller, while the dashdot lines are the wind turbine responses when using a PI controller (as designed in
[6]
). The final two plots are the tower foreaft bending moments (My
_{T}
) for the LQR and PI controls, respectively. By applying the LQR control, the peaktopeak magnitude of (My
_{T}
) decreases a lot. The performance differences of these two responses between the two controllers are summarized in
Table I
. As shown, the rotor rpm of the LQR control is a little more oscillatory than that of the PI control, but the output powers of two cases are almost the same.
PERFORMANCE DATA AT WIND SPEED OF 8M/S TI 23.2%
PERFORMANCE DATA AT WIND SPEED OF 8M/S TI 23.2%
It is difficult to differentiate between the structural responses such as the blade or tower bending moment in the time domain. The damage equivalent loads (DEL),
M_{eq}
, can be a quantitative measure for this ambiguous situation.
M_{eq}
is given by (25):
where
n_{k}
is the number of cycles in the mechanical load range
M_{k}
and
n_{tot}
is the total number of cycles in the mechanical load signal.
m
is a material specific number, for example
m
= 3.5 for a steel tower structure and
m
= 10 for a fiber glass blade [10]. The larger the DEL value, the more prone it is to fatigue failure. The statistics for the DEL data of the structural load responses at a mean wind speed of 8m/s for both the LQR and PI controllers are summarized and compared in
Table II
. The DEL of the blade bending moment in the inplane (Mx
_{B}
) direction as well as the tower root bending moment in the foreaft (My
_{T}
) and sidetoside direction (Mx
_{T}
) are significantly improved.
Simulation result for wind speed of 8m/s 23.2% TI.
MECHANICAL LOAD AT WIND SPEED OF 8M/S TI 23.2%
MECHANICAL LOAD AT WIND SPEED OF 8M/S TI 23.2%
Fig. 12
shows the time domain response simulation for the above rated wind speed region. The wind speed in this simulation has a mean value of 18m/s with a TI of 16.97%. The first window of
Fig. 12
shows the wind speed measured at the hub. The second window shows the rotor speed. The next windows gives the pitch, generator torque, and electric power. Again, the solid lines are for the LQR, and the dashdot lines are for the PI control.
Fig. 13
illustrates the inplane blade root bending moment (Mx
_{B}
, the first window), the outofplane blade root bending moment (My
_{B}
, the second window), the tower root bending moment in the sideward direction (Mx
_{T}
, the third) and the tower root bending foreaft moment (My
_{T}
, the last one). The final two windows are enlarged plots for My
_{T}
during the 100 second time period. These are prepared to compare the effectiveness of the LQR control with the PI control. As can be seen, the foreaft tower root bending moment (My
_{T}
) for the LQR control becomes much smaller than that of the PI control. The differences in performance and mechanical stresses between the two controllers are summarized in
Table III
and
IV
. It can be seen in
Table III
, that the LQR has a slight decrease in electrical power. However, the power fluctuations decrease remarkably. Furthermore, there is more than a 14% DEL decrease in the tower bending moment in the foreaft direction and nearly 9% in the sideward direction, as shown in
Table IV
.
PERFORMANCE AT WIND SPEED OF 18M/S TI 16.97%
PERFORMANCE AT WIND SPEED OF 18M/S TI 16.97%
MECHANICAL LOAD AT WIND SPEED OF 18M/S TI 16.97%
MECHANICAL LOAD AT WIND SPEED OF 18M/S TI 16.97%
Simulation result for wind speed of 18m/s 16.97% TI.
Mechanical loads at wind speed of 18m/s 16.97% TI.
The statistical data for the electrical power and tower mechanical loads in the entire operating region are presented in
Fig. 14
,
15
and
16
. It is noted that the simulation at a wind speed of 7m/s has been done with a very high turbulence intensity so that it covers all of the low wind speed region in the operating region of the wind turbine. Although the amount of power captured slightly decreases, the mechanical loads of the tower in both the foreaft direction and sidetoside direction are improved in the whole region. How much of a decrease in the mechanical loads or how much power loss, can be managed by appropriate values for the Q and R matrix in the quadratic cost function. However, there is a no direct relashionship between the performance and mechanical stresses. The proposed LQR controller is designed to get a similar performance response with a PI controller in addition to improving the tower load reduction. This design process requires some trial and errors in selecting the Q and R matrix.
Electrical power.
Tower root bending moment in sideward direction.
Tower root bending moment in foreaft direction.
VI. CONCLUSIONS
This paper investigated the designing of a windscheduling linear quadratic controller with online wind speed estimation for multiMW size wind turbines. The controller was designed by linearising the nonlinear wind turbine along the operating point trajectory. The potential of the controller was checked by simulations with GH Bladed software. The responses of the proposed method were compared with those of the conventional PI controller. The statistical data show that the proposed controller has tradeoffs between capturing energy and alleviating mechanical loads.
Acknowledgements
This work was supported by the New & Renewable Energy program (No. 2012T100201670) and the Human Resources Development program(No. 20134030200240) of the Korea Institute of Energy Technology Evaluation and Planning(KETEP) grant funded by the Korea government Ministry of Trade, Industry and Energy. And, this study was supported by 2013 Research Grant from Kangwon National University (No. 120131305).
BIO
Yoonsu Nam was born in Suwon, Korea, in 1958. He received his B.S. in Nuclear Engineering from Seoul National University, Seoul, Korea, in 1981. He received his M.S. from the Department of Machine Design, Seoul National University, in 1983. He received his Ph.D. in Mechanical Engineering from Georgia Tech, Atlanta, GA, USA, in 1991. From 1983 to 1986, he was with the Central Research Institute of GoldStar (LG Electronics) Inc., Korea. From 1992 to 1997, he was with the Flight Dynamics and Control Laboratory of the Agency for Defense Development, Daejeon, Korea. In 1997 he joined the School of Mechanical and Mechatronics Engineering, Kangwon National University, Chuncheon, Korea. His current research interests include mechanical load reduction control of MW class wind turbines.
Pham Trung Kien was born in Vietnam, in 1979. He received his B.S. and M.S. from the School of Electrical Engineering, Hanoi University of Science and Technology, Hanoi, Vietnam, in 2002 and 2006, respectively. His current research interests include control system theory and design, wind turbine control system design and wind power systems.
YoHan La was born in Gyeonggido, Korea, in 1988. He received his B.S. in Mechanical and Mechatronics Engineering from Kangwon National University, Chuncheon, Korea, in 2011. He is currently working toward his M.S. at Kangwon National University. His current research interests include mechanical load reduction control of MW class wind turbines.
Bossanyi E. A.
2000
“The design of closed loop controllers for wind turbines,”
Wind energy
3
(3)
148 
163
DOI : 10.1002/we.34
Ostergaard K. Z.
,
Brath P.
,
Stoustrup J.
2007
“Gainscheduling linear quadratic control of wind turbines operating at high wind speed,”
16th IEEE International Conference on Control Applications
Singapore
276 
281
Stol K.
,
Balas M.
2001
“Fullstate feedback control of a variableSpeed wind turbine: A comparison of periodic and constant gains,”
ASME J. Solar energy engineering
123
(4)
319 
326
DOI : 10.1115/1.1412237
Muhando E. B.
,
Senjyu T.
,
Kinjo H.
,
Funabashi T.
2008
“Augmented LQG controller for enhancement of online dynamic performance for WTG system,”
Renewable energy
33
(8)
1942 
1952
DOI : 10.1016/j.renene.2007.12.001
van der Hooft E. L.
,
Shaak P.
,
van Engelen T. G.
2003
“Wind turbine control algorithms,”
ECNC03111
Nam Y.
,
Kim J.
,
Paek I.
,
Moon Y. H.
,
Kim S. J.
,
Kim D. J.
2011
“Feedforward Pitch control using wind speed estimation,”
Journal of Power Electronics
11
(2)
211 
217
DOI : 10.6113/JPE.2011.11.2.211
Bossanyi E. A.
2009
“GH Bladed theory manual(version 3.81),” Garrad Hassan and Partners
282BR009