Control-oriented modelling of wind turbines using a Takagi-Sugeno model structure

For a horizontal-axis wind turbine (HAWT), a dynamic nonlinear model with four degrees of freedom is derived and transformed into a Takagi-Sugeno (TS) model structure using the sector nonlinearity approach. Thereby, an exact transformation of the nonlinear model is obtained as a weighted combination of linear models. This structure allows for a convenient design of controller and observer structures. The maps of the rotor thrust and torque coefficients can be implemented in the model as look-up tables or, alternatively, as analytical nonlinear functions. Open-loop simulation results of the derived TS model for a reference model turbine are compared to those obtained with the aero-elastic code FAST. The small deviations obtained demonstrate the high model quality of the control-oriented TS model. In future work, the derived TS model shall be used as a basis for the design of fault detection and isolation (FDI) concepts.


I. INTRODUCTION
As wind turbines are gaining a growing significance in the global energy supply, the demand for advanced control strategies aiming at higher efficiency and availability as well as load reduction also rises. Therefore, it is essential to develop model-based control algorithms that exceed the possibilities of classical frequency-domain control design concepts. Based on dynamic wind turbine models with few degrees of freedom, various control design concepts have already been employed to address different issues arising in wind turbine control. For example, LQG controller concepts were used in [1] and [2] to reduce loads. Classical gain-scheduling techniques have been replaced by more rigorous linear parametervariable (LPV) approaches to design the wind turbine speed control for the different operating regions in [3].
The purpose of this paper is to present a wind turbine model that shall be used as a basis for fault detection and isolation (FDI) concepts in wind turbine control in future work. By using the Takagi-Sugeno model structure with sector nonlinearities, an exact representation of the nonlinear model can be obtained as a weighted combination of linear statespace models. For this structure, both controllers and observers can be designed by solving linear matrix inequalities (LMIs) [4]. Fault detection schemes can be developed by extending the standard TS observer structures to sliding mode TS observers [5]. In [6], this is used for the detection of sensor faults in the pitch system of wind turbines. The same methodology shall be applied in order to detect and isolate more sensor and actuator faults occurring in the operation of wind turbines, for which the model presented here serves as a starting point.
The derived wind turbine model is validated against the aero-elastic simulation code FAST by NREL [7] using a 5 MW reference wind turbine described in [8].
There are a lot of papers that make use of a TS model structure for wind turbine models and controller design. For example, in [9] maximum wind power tracking is tackled by designing a TS controller along with an observer for a turbine model including an explicit generator/converter model. In [10] a robust power tracking controller is designed using linear matrix equalities to account for disturbances and parametric uncertainties. In [11] a Fuzzy clustering approach is applied to obtain a TS model in order to achieve maximum energy extraction. In most cases, however, the controllers are not validated against more detailed structural simulation models (like FAST) normally used in the wind energy industry. This, however, is important in order to test both the quality of the design models and the controller performance with more realistic models. This paper is organised as follows: In Section II, the wind turbine model is derived and written in state-space form. In Section III, the state-space model is transformed into a TS model structure. In Section IV, model parameters for the 5 MW reference wind turbine are given including the nonlinear aero maps for the rotor thrust and torque coefficients. In addition, these maps are approximated by nonlinear functions. In section V, simulation results are presented which are compared to the results obtained with the FAST aero-elastic simulation. Conclusion and outlook are given in section VI.

II. WIND TURBINE MODEL
The aero-elastic codes normally used for wind turbine load simulation contain models with more than twenty degrees of freedom which are derived by employing modal analysis techniques. For control design purposes, these models are too complicated and also capture dynamic effects that are not directly influenced by the control action. For these reasons, the models used for control design have to be as simple as possible but must capture the dominant system dynamics [3]. Fig. 1. Schematics of the wind turbine model with the submodels for aerodynamics, mechanics and pitch and their respective inputs and outputs. v: wind speed; F T , Tr: Rotor thrust and torque; ωr: Rotor angular velocity; Tg: Generator torque; β: Pitch angle; β d : Demanded pitch angle; x: System state vector. The model in [12] or variations of it are useful for this aim. In this paper, the model as presented in [3] is used, with slight modifications. The complete model consists of three submodels: the mechanical submodel, which includes drivetrain and structure, the aerodynamics submodel, and the pitch actuator submodel. These submodels are coupled, which is illustrated in Fig. 1. An explicit generator/converter submodel is omitted here, i.e. an ideal frequency converter is considered, where the demanded generator torque is equal to the actual generator torque. However, it would be straightforward to include also converter dynamics, for example as a first order delay model.

A. Mechanical Submodel
For the mechanical submodel, four degrees of freedom are considered: rotor and generator angles, as well as fore-aft tower top deflection and flapwise blade tip deflection.
The vector of the generalised coordinates is therefore given by q := (y T y B θ r θ g ) T , where y T and y B denote the tower top and blade tip deflections, while θ r and θ g denote the rotor and generator rotational angles. The vector of external forces is given by where F T denotes the aerodynamic rotor thrust force, T r the aerodynamic rotor torque and T g the applied generator torque.
For the tower and blade degrees of freedom, the mechanical models are illustrated in Figs   coefficients for the tower top and blade tip motion. d T and d B denote the damping coefficients for tower and blade, which can be estimated from the aerodynamic rotor damping (see section IV).
The drive-train model is illustrated in Fig. 4. The rotor and generator are modelled as inertias J r and J g , respectively. ω r and ω g denote the rotor and generator angular velocities, which are the time-derivatives of the rotational angles θ r and θ g . k S and d S denote the torsional stiffness and damping coefficients of the shaft. As gearbox dynamics is not considered here, a gearbox ratio could easily be included into the model by appropriately modifying the generator inertia J g and the generator torque T g .
The derivation of the differential equations governing the behaviour of the mechanical model is done by employing the formalism of Lagrangian mechanics. The derivation is omitted here, as it can be found, for example, in [3].
The equations of motion, which describe the dynamics of the mechanical model, are obtained in matrix form as where the mass matrix M, the damping matrix D and the stiffness matrix K are given by

B. Pitch actuator submodel
The pitch actuator dynamics is modelled as a first-order delay system [3], [12]: where β and β d denote the actual and demanded pitch angle and τ denotes the delay time constant. This model is a simplified, but reasonable assumption for the overall pitch actuator dynamics, which subsumes the dynamics of all internal control structures normally used in hydraulic or electromechanical pitch systems.

C. Aerodynamics submodel
The aerodynamics submodel comprises the expressions for the rotor thrust force F T and the rotor torque T r . These depend on the aero maps C T and C Q for the thrust and torque coefficients [3], [13]: where R denotes the rotor radius, ρ the air density, v the wind speed and λ = R ωr v the tip speed ratio. 2

D. State-space model
In order to design controllers or observers for the wind turbine model, it is necessary to write the model equations in state-space form. i.e.,ẋ = f (x, u), where f is a function of the state vector x and the input vector u. The output equation y = C x describes which states of the system are actually measurable. The measurable states form the output vector y. The choice of the output matrix C will become important for an observer design. For the purpose of this paper, it can be chosen arbitrarily, as only the dynamics of the state vector is considered.
A state-space form for the complete wind turbine model can be obtained by inserting the aerodynamic equations (6) and (7) into the mechanical model equations (1) and adding the pitch model (5).
The resulting state-space model iṡ with the state-vector and the output vector y = x, i.e. the output matrix is given by C = I 9×9 , where I denotes the unit matrix. The wind speed v is considered as a disturbance input.
The linear part of the state-space model is given by the system matrix and the input matrix In addition, there is a vector that nonlinearly depends on the system states and the wind speed: III. TAKAGI-SUGENO MODEL STRUCTURE As the model (8) is nonlinear, standard state-space techniques for controller or observer design cannot be directly applied. The main advantage of a TS model structure is its flexibility. In this paper, a TS model is derived by means of the sector nonlinearity approach which provides a way to obtain an exact representation of the full nonlinear model as a weighted combination of linear submodels, where the nonlinearities of the system are shifted into the membership functions. Another approach would be to derive an approximated TS model by linearising the nonlinear model equations around different operating points and using predefined fuzzy membership functions to combine the linear submodels to an overall nonlinear model. For a system in TS model structure, stable controllers can be designed by solving LMIs where conditions for the overall decay rate or for optimal controller design can be included. Because of these advantages, TS model structures are an ideal basis for the design of controllers and observers for FDI concepts.
The general TS model structure for systems without direct feed-through is [4] x with N r = 2 N l , where N l denotes the number of differing nonlinearities in the model. z is the vector of known premise variables, which may be functions of the state variables, external disturbances and/or time [4]. The matrices A i , B i and C i have constant coefficients and the membership functions h i fulfill the condition The nonlinear wind turbine model (8) will now be transformed into a TS model structure, thus the nonlinear dependency of g(x, v) in (12) has to be transformed into the membership functions h i , such that the subsystems A i x + B i u are linear state-space models. In order to achieve this, g(x, v) can be written as (15) such that (8) can be rearranged aṡ In principle, the nonlinear function entries can be placed in any column of the matrix A N L , they only have to be divided by the respective coordinate of the state-vector, such that the vector g(x, v) is correctly reproduced. Here, the nonlinear functions were placed in the seventh column, because the state x 7 =θ r = ω r , i.e. the rotor angular velocity, can always safely be assumed to fulfill ω r > 0. Even if starting procedures of the wind turbine shall be explicitly examined, the initial rotor angular velocity can be chosen as a small positive value, so that ω r > 0 is fulfilled for the complete simulation time. However, other states would not be equally well-suited as a dividing factor for the nonlinear functions in the matrix A N L . For example, the states x 1 = y T , x 2 = y B , x 5 =ẏ T , x 6 =ẏ B can all oscillate around 0 in normal operation and therefore cannot be used.
In the combined matrixÃ = (A + A N L ), at positions (6,7) and (7,7) there are now two scalar, nonlinear functions: If these two functions are bounded, i.e. f i ∈ [f i , f i ], they can be written as sector functions in order to arrive at a TS model structure: where the weighting functions w jk are defined by This is the sector nonlinearity approach [4], [14]. Equations (19) and (20) are exact expressions for the functions f 1 and f 2 if these are bounded, which is the case for f 1 and f 2 as defined in (17) and (18). From equations (17) and (18) Here, a maximum wind speed of v max = 60 m s and a minimum rotor angular velocity of ω r,min = 0.01 rad s was assumed. From (19), (20) and (21) it follows that the weighting functions fulfill the convexity condition It also follows that (14) for the TS membership functions is fulfilled. Using the expressions (23) and (24) instead of (17) and (18) in the matrixÃ and multiplying all other matrix entries by which is an exact representation of (8). The entries of the submatrices A i are the coefficients f j , f j of (23) and (24) for the matrix positions (6,7) and (7,7). For all other positions, the entries of A i are equal to those ofÃ.
As the nonlinearities were included in the new matrixÃ, the term B u, as well as the output equation y = C x, remain unaltered. This special TS structure with only one input matrix B is advantageous for controller design, as the resulting number of submodels in the closed-loop system is equal to N r , compared to the more general TS structure (13), where N 2 r submodels are obtained for the closed-loop system.

IV. MODEL PARAMETERS
For simulation studies, the TS wind turbine model (25) is validated against an aero-elastic simulation with FAST of the NREL 5 MW reference turbine defined in [8]. To this end, the model parameters for the 5 MW reference turbine are also chosen for the TS wind turbine model. Some parameters can be directly taken from [8] and from NREL example input and log files for the 5 MW turbine. These are: N = 3, R = 63 m, ρ = 1.225 kg m 3 , m B = 17740 kg, m T = 644240 kg, 3 , J r = 38759227 kg m 2 , J g = 5025347 kg m 2 , k s = 867637000 Nm, d s = 6215000 Nm s.
As a gearbox ratio is included in the FAST simulation model, the generator inertia given in the FAST model (J g = 534.1 kg m 2 ) was multiplied by the square of the gearbox ratio (n g = 97) to obtain the correct value for the generator inertia J g in the TS model, where no gearbox ratio is included.
The damping factors d T and d B for tower top and blade tip motion are mainly determined by the aerodynamic damping of the rotor, which can be estimated as [13], [15]: where the coefficient d * depends on the tip speed ratio. For a rotor with a similar design tip speed ratio (λ d ≈ 6) as the 5 MW reference turbine (λ d ≈ 7.6), the values of d * lie between -0.5 (in the stall region around λ ≈ 4) and 2 (from [16], chapter 25, lit. 10). For the simulation studies, d * was set to 1 (independent of the tip speed ratio). This is a reasonable approximation, as it is hard to determine exact values for the damping and the dynamic behaviour will be well approximated if the orders of magnitude for the damping coefficients are correct. Following the same argument, it is reasonable to calculate the aerodynamic damping factor in (26) using a constant wind speed. Here, a wind speed of 8 m/s was used to obtain d aero = ρ 2 R 2 v ≈ 2 · 10 4 Ns m and the damping factors d T and d B for tower and blade were set equal to this value.
The stiffness coefficients k T and k B for the tower top and blade tip motion cannot be directly taken from [8], as the tower and blade bending stiffness coefficients for the FAST model are given separately for each tower and blade section. It would be possible to obtain effective bending stiffness coefficients for the whole tower and the whole blade by applying a transition matrix technique used in structural mechanics calculations. This was omitted here, as it would be beyond the scope of this work. Instead, k T and k B were chosen such that the mean tower top and blade tip deflections y T and y B for an idling wind turbine at a stationary operating point λ ≈ 17.85 assume approximately the same values as the respective degrees of freedom in the FAST simulation. This procedure was done for a constant wind speed of v = 8 m s at a constant pitch angle of β = 0 • and generator torque T g = 0 Nm. The values obtained are: k T = 1962000 N m and k B = 39333 N m . The delay time constant for the pitch actuator model (5) was estimated considering a pitch velocity limit of 10 deg/s (as in [2]) and a reference step jump of 1 deg and can thus be set to τ = 0.1 s.

Aero Maps for Rotor Thrust and Torque Coefficient
The aero maps C T (λ, β) and C Q (λ, β) for the rotor thrust and torque coefficients were extracted from FAST simulations for different pitch angles. In each simulation, the generator torque was set to 0 Nm, and the initial rotor angular velocity to 0 rpm, such that the wind turbine assumes all values of the tip speed ratio from λ = 0 up to the stationary idling tip speed ratio (λ ≈ 17.85), where an equilibrium of accelerating and decelerating aerodynamic forces has built up. In this way, the aero maps can be obtained for the whole range of tip speed ratios and pitch angles as look-up tables 4 . For pitch angles between 0 • and 20 • , the obtained curves for C Q are depicted in Fig. 5a, those for C T in Fig. 6a.

Analytical Approximation of the Aero maps
For the simulation of the TS wind turbine model (25), the thrust and torque coefficients can be calculated in each step by simply interpolating the values in the look-up tables of the aero maps, because C Q and C T appear in the membership functions h i .
However, it is useful to find good analytical approximations for C Q and C T , especially if the controller design shall be based on linearised submodels. In that case, derivatives of the aero maps are needed and differentiating analytical functions is advantageous to calculating numerical derivatives of tabulated look-up tables. For the rotor torque map C Q , the function given in [18] ( [19], [20]) for approximating the power coefficient was modified to be used for C Q : with (λ > 0) and To exclude unphysical negative terms, this expression was limited to values ≥ 0: Using the MATLAB R curve fitting toolbox and manual corrections, the following coefficients were obtained :   5b shows the plots of function (29) with these coefficients for pitch angles between 0 • and 20 • compared to the plots of the curves from the look-up tables in Fig. 5a. For the rotor thrust map C T , the following nonlinear function was used to approximate the look-up table.
As for C Q , this expression was limited to values ≥ 0: The following coefficients were obtained:  Fig. 6b shows the plots of function (31) with these coefficients compared to the plots of the curves from the C T look-up tables in Fig. 6a.

V. SIMULATION RESULTS
The TS wind turbine model (25) was implemented in MATLAB/Simulink R and simulation runs of the open-loop dynamics were done around different operating points, i.e., fixed values of pitch angle and generator torque, using IEC wind gusts as input signals. 5 . For the C Q and C T curves, both the look-up tables and the analytical approximations (29) and (31) were used. The results were compared to FAST simulations at the same operating points. 6 For wind gusts around 8 m/s and 18 m/s, results are depicted in Fig. 7, with time-series of tower top and blade tip displacement (y T and y B ), rotor angular velocity (ω r ) and the torsion angle (θ s = θ r − θ g ).
The results show a good agreement between the FAST simulations and the simulations of the TS model using the C Q and C T look-up tables. However, there are some deviations: 1) The amplitudes of the oscillations in the tower and blade degree of freedom are higher in the TS model than in FAST. This is partly due to the fact that no dynamic correction of the wind speed was included. Including such a correction in turn leads to higher deviations in the rotor rotational velocity. 2) For the blade tip deflection, the deviation between the TS model and FAST is very pronounced for high tip speed ratios (not shown here). In FAST, a centrifugal stiffening of the blades is taken into account by including terms proportional to ω 2 r in the blade stiffness coefficients. In first tests, it could be verified that including centrifugal terms in the state-space model (8) yields better agreements with the FAST simulation results. In Fig. 7 it can also be seen that the blade tip deflection in the FAST simulation, where each blade is treated separately, oscillates with the rotor angular frequency, which is probably due to gravitational effects that have also not been included in the TS model.   Using the analytical approximations for the C Q and C T maps, the agreements to the FAST simulation results are still acceptable, however, the TS model simulations show higher deviations, which is simply due to the difficulty of finding analytical functions that yield very precise fits to the complicated nonlinear aero maps. Still, the functions given in (29) and (31) should be good candidates in cases where analytical expressions for the aero maps are needed.

VI. CONCLUSION AND OUTLOOK
In this paper, a four-degree of freedom wind turbine model was derived and transformed into a Takagi-Sugeno (TS) model structure using the sector nonlinearity approach. The TS model exactly represents the nonlinear model as a weighted combination of linear models.
The open-loop dynamics of the TS model were compared to results obtained using the simulation code FAST for the NREL 5 MW reference turbine. The overall agreement for the different degrees of freedom is good, which demonstrates that the reduced 4 DOF model is appropriate for controller and observer design. To further improve the model quality, centrifugal and gravitational effects could be included. However, each model improvement must be carefully traded off against the aim of using the model for controller and observer design, as additional nonlinearities increase the number of TS submodels and can thereby lead to more conservative LMI design results.
The TS model structure can be extended to design TS sliding mode observers and therefore serves as an ideal startingpoint to develop FDI concepts for wind turbines in future work.