Effective wind speed estimation based on a data-driven model of wind turbine tower deflection

Abstract

The deflection of the wind turbine tower can provide us with rich information about the effective wind speed. In this paper, a new method for effective wind speed estimation is proposed based on tower deflection. The tower vibration model is derived and a subspace identification method is used to identify the model parameters. An online estimator of aerodynamic thrust force based on the identified tower model is designed and then implemented using a Kalman filter together with a recursive least squares algorithm. The estimated aerodynamic thrust force is then used as an input to a neural network estimator, which is trained to invert the aerodynamic thrust force equation and estimate the effective wind speed. In order to show the performance of the proposed estimator, the estimated thrust force and wind speed are compared and verified with a third-party simulation data of a 1.5 MW wind turbine. The comparison shows close agreement between their values.

Introduction

The objectives of control algorithms in wind turbines include tracking maximum power points for a wind speed, which is less than the rated wind speed, and avoiding over-speeding of the rotor when wind speeds are higher than the rated values. The first goal can be accomplished by tuning the shaft speed through controlling the torque and the second goal can be achieved by adjusting the pitch angle through pitch control.

Since the wind speed is variable during the turbine operation, the wind speed is required to be measured or estimated for the control system. Because increasing the accuracy of wind speed measurement or estimation increases the performance of the control system, several techniques have been applied and evolved during the past years including the simplest approach of using anemometers. Anemometers are commonly used for wind speed measurement. However, since they are installed on top of the nacelle and out of the rotor plane, they do not provide the effective wind speed driving the turbine blades. Moreover, the wind speed, changed by the blades passing the front side of the anemometer, is not the wind speed on the turbine blades. Laser-based measurement systems include high cost and challenges with heavy rains, dense fog or dust conditions, and sunlight reflections.

The wind speed spatial average over the rotor plane gives the effective wind speed, and recent estimation methods have been trying to find that value. Rettenmeier et al. ¹ provided a comparison of single point anemometer measurements and estimations of the effective wind speed. Effective wind speed estimation has attracted a lot of attention, recently, and several methods other than anemometers have been proposed in order to provide more accurate estimates.^2,3

A number of wind speed estimation methods have been developed based on the dynamics of the driveline of a wind turbine. In those methods, several tools have been used in order to mathematically relate measurable dynamic parameters to the effective wind speed. Bhowmik et al. employed the polynomial that relates power coefficients to the tip speed ratio. ⁴ For any given turbine power and rotor speed, wind speed can be estimated by solving the polynomial equation. In another research study, ⁵ a multi-layered perceptron neural network was trained to also solve for wind speed from given power and rotor speed. In similar studies, an extreme learning machine was used.^6,7 Boukhezzar and Siguerdidjane ⁸ used a Kalman filter to estimate the aerodynamic torque and rotor speed, and then solved the aerodynamic torque equation by NR to estimate the effective wind speed. Hafidi and Chauvin ⁹ applied a linear state observer to estimate the aerodynamic torque, then a least squares method was used to estimate the effective wind speed by minimizing the aerodynamic torque estimation error. In summary, the existing methods used rotordynamics to estimate the wind speed. However, the high rotating inertia of the wind turbine means a potential time delay between the rotational response and the wind disturbance.^10,11 Additionally, the rotor dynamics system behaves like a low-pass filter which tends to filter out the rapidly varying components in the wind speed.

In order to estimate the effective wind speed, especially for higher-frequency disturbance content, we propose to include tower dynamics in addition to the rotor dynamics. The tower dynamics resembles an under-damped 2nd-order system because of its construction material. ¹² It has a quicker response time and also a wider passband. Using the tower dynamic model to estimate the wind speed makes estimation of the higher frequency wind speed fluctuations feasible. These can help improve the performance of control systems.

To construct the dynamics model for the deflections of wind turbine tower, both aerodynamics and structural dynamics should be included. For the first step, which is modeling the aerodynamics of the wind turbine, blade element momentum (BEM) Theory can be used.^13–15 The BEM method is fast and computationally cheap, which can provide reliable airfoil lift and drag coefficients, which leads to accurate data. ¹⁶ The structural dynamics of the wind turbine tower can also be modeled, and one can identify the model parameters using input-output data of the system. By formulating this model, an inverse problem can be solved to find the effective wind speed from tower deflections.

Solving an inverse problem to estimate the excitation forces based on the dynamic system response is frequently employed and several methods are proposed in the literature.^15,17–19 Nordstrom ²⁰ has used a dynamic programming algorithm to solve the least squares problem by minimizing the regularized squared error function and estimating the force history. The Sum of weighted accelerations (SWA) technique was developed by Carne and successfully applied to the force estimation problem. ²¹ SWA is based on the theory of a modal filter where a spatial filter constructed based on system mode shapes is used to isolate the rigid-body response and estimate the excitation force. The main drawback of these methods is that online load estimation is not possible i.e. the full sequence of the measurement signal is required before running the dynamic programming algorithm. ²² The Inverse structural filter technique^23,24 inverts the discrete-time form of the equations of motion, which results in a dynamic system that receives the response of the structure as an input and then gives back an estimation of the forces applied to the structure as an output. This method can be used for online force estimation. However, it is reported to be sensitive to measurement noise.^22,25 The above-mentioned methods are based on a deterministic approach assuming the error to be present only in the measurements. Some combined deterministic-stochastic methods are developed assuming noise on the measurements and the state variables. The unknown Input observer method is a combined deterministic-stochastic method that can estimate both unknown inputs and states simultaneously. ²⁶ The augmented “state observer” technique includes the unknown force in the state vector and uses a Kalman filter to estimate the unknown force in conjunction with states. ²⁷

Most of the recent research on wind speed estimation has focused on rotor speed measurements. There are many complications, such as brake engagement, generator load variations, and especially large inertia in the drivetrain of a wind turbine, which can affect the accuracy of the estimated wind speed. On the other hand, by recent advances in measurement tools,^28–30 such as wireless nodes and MEMS motion sensors, it is feasible to measure the deflection of the tower, which can be used for several purposes, such as tower fatigue prediction. In this work, we propose to use the wind turbine tower dynamics for wind speed estimation. Since the deflection of a wind turbine tower is caused by the wind, by creating a simple dynamics model of the wind turbine tower and then solving the inverse problem, the effective wind speed can be estimated. The thrust force of the wind turbine is calculated using the wind speed, rotor speed, blade pitch angle and BEM. Then, we use the aerodynamic thrust force as an input to the structural model of the tower. A Kalman filter along with a recursive least squares method, adopted from a different problem, ³¹ is then used to provide online estimation of the thrust force assuming a combined deterministic-stochastic system. In other words, the presented estimation method consists of two parts. In the first part we use a Kalman filter to estimate the innovation and residual sequences, ³² defined as the difference between the estimated and measured output, with and without the effect of the input. This results in an equation that can be solved in the second part using a recursive least squares method to estimate the unknown input. Having the estimated thrust force from the previous step, the second part of the model, which is governed by wind turbine aerodynamics, is inverted to estimate effective wind speed. The main challenge is that the aerodynamic thrust force, which we modeled using BEM theory, is nonlinear and must be solved iteratively. Artificial neural networks (ANN), which have been used successfully to solve nonlinear inverse problems, ³³ will be use to solve this inverse problem. Having a set of simulation data, we train an artificial neural network in order to obtain the relation between the inputs and outputs.

Theoretical background

Because of more complexity of the aerodynamics of the wind turbines, by splitting the system into two parts we can consider the dynamic model of the wind turbine as a so-called Hammerstein system. ³⁴ The Hammerstein structure describes systems that consist of a static nonlinearity in series with a linear dynamic part. As in Figure 1, the dynamic system of the wind turbine tower is represented as a Hammerstein system where the nonlinear mapping from wind speed $U$ , pitch angle $\beta$ , and tip speed ratio $\lambda$ , to aerodynamic thrust coefficients is considered static followed by a linear dynamic system of tower fore-aft motion. In that figure, $T_g$ is the generator torque. $F_T$ and $F_Q$ , which will be calculated later, are the blade aerodynamic thrust force and torque, respectively. $D_{SS},\mathrm{ }{D}_{FA},\mathrm{ }R,\mathrm{ }\mathrm{\Omega },$ and $\omega$ are the tower static and fore-aft deflections, blade radius, rotor speed, and the generator speed, respectively. The entire dynamic system of the wind turbine is shown under a closed-loop control. However, open-loop system identification methods still can be used for identification of the tower dynamics because there is no feedback from $D_{FA}$ to $F_T$ .

Figure 1.

The windturbine dynamic system.

By doing so, we can use the well-established LTI system identification techniques to identify a model for the wind turbine. The other advantage of this method is that we are not restricted to measure the input-output data around a fixed operating point, therefore, we can record longer sequences of data for system identification. Also, the identified system will be globally valid in the operating domain of the wind turbine, which is vital for wind speed estimation.

Tower fore-aft motion dynamics

In this section a dynamics model for wind turbine tower deflection is presented. It will be used in the following sections on system identification of the wind turbine tower deflection.

The fundamental equation describing the dynamic behaviour of the wind turbine tower deflection, as discretized by the Finite Element method, can be written as ³⁵

$$\mathbf{M}\ddot{\mathbf{q}}\left(t\right)+\mathbf{B}\dot{\mathbf{q}}\left(t\right)+\mathbf{Kq}\left(t\right)=\mathbf{g}\left(t\right)$$

(1)

where $\mathbf{M}$ , $\mathbf{B}$ and $\mathbf{K}$ are the mass, damping and stiffness matrices respectively, $t$ is the time, $\mathbf{q}\left(t\right)$ is the state vector that collects the structural displacement degrees of freedom, while $\mathbf{g}\left(t\right)$ is the vector that indicates the corresponding applied loads.

As it is derived in Appendix 2, equation (1) in its general form and subject to measurement and/or process noise can be described by a discrete-time LTI state-space model in an innovation form as follows

$${\mathbf{x}}_{k+1}=\mathbf{A}{\mathbf{x}}_k+\mathbf{B}{\mathbf{u}}_k+\mathbf{K}{\mathbf{e}}_k$$

(2)

$${\mathbf{y}}_k=\mathbf{C}{\mathbf{x}}_k+\mathbf{D}{\mathbf{u}}_k+{\mathbf{e}}_k$$

(3)

where $\mathbf{A}\in {\mathbf{R}}^{n\times n}$ is the system matrix given by equation (50), $\mathbf{B}\in {\mathbf{R}}^{n\times n_u}$ is the input matrix given by equation (51), $\mathbf{C}\in {\mathbf{R}}^{n_y\times n}$ is the output matrix given by equation (47), $\mathbf{D}\in {\mathbf{R}}^{n_y\times n_u}$ is the feedthrough matrix which is a zero matrix in this case, $\mathbf{K}\in {\mathbf{R}}^{n\times n_y}$ is the Kalman gain matrix, ${\mathbf{x}}_k\in {\mathbf{R}}^n$ is the state vector containing the wind turbine tower fore-aft displacements and velocities, ${\mathbf{u}}_k\in {\mathbf{R}}^{n_u}$ is the input vector containing the exerted aerodynamic thrust force, ${\mathbf{y}}_k\in {\mathbf{R}}^{n_y}$ is the output vector containing the measured tower tip fore-aft deflection and ${\mathbf{e}}_k\in {\mathbf{R}}^{n_y}$ is an unknown innovation with covariance matrix $E\left[{\mathbf{e}}_k{{\mathbf{e}}_k}^T\right]={\mathbf{R}}_e$ . In this representation index $k$ is the discrete time step, $n$ , $n_u$ and $n_y$ are the dimensions of the state, input and output vectors in the space of real numbers $\mathbf{R}$ .

The discretized state-space model is suitable for wind speed estimator design. However, it is hard to determine $\mathbf{A}$ , $\mathbf{B}$ and $\mathbf{C}$ matrices accurately using the presented method. The reason is that these matrices contain several physical parameters of the system which cannot be measured or estimated accurately. In the upcoming section, the system identification technique will be used to determine the model parameters precisely.

Blade element momentum theory

The BEM method calculates the aerodynamic torque and thrust force by splitting the wind turbine blades into several aerodynamically independent units, called blade elements, and using conservation of linear and angular momentum.

BEM consists of two different theories, namely, momentum theory and blade element theory. First, we derive two sets of equations based on the momentum and blade element theories for aerodynamic thrust force and torque. Then we equate aerodynamic torque and thrust force values from the two sets of equations, which results in an iterative algorithm for calculating induction factors and obtaining aerodynamic thrust force and torque. ³⁶ In this paper, the iterative algorithm, provided by Moriarty and Hansen, ³⁶ is implemented.

The first part of BEM, the momentum theory, considers an annular stream tube as a control volume. It uses conservation of linear and angular momentum to obtain the axial and angular induced velocities and to derive the forces on the wind turbine blade. The second part of BEM, the blade element theory, divides the blade into several annular elements operating independently under local aerodynamic conditions assuming that there is no interaction between blade elements. The blade elements are considered as two-dimensional airfoils whose aerodynamic forces are defined as a function of blade geometry using the drag and lift coefficients.

As a result of BEM, the wind turbine thrust force and torque can be defined in terms of non-dimensional thrust and torque coefficients, $C_T$ and $C_Q$ , which are only functions of pitch angle, $\beta ,$ and the tip speed ratio, $\lambda$ , which is defined as the ratio of blade tip speed, $\mathrm{\Omega }R,$ to wind speed, $U$

$$\lambda =\frac{\mathrm{\Omega }R}{U}$$

(4)

Here $R$ is the blade radius and $\mathrm{\Omega }$ is the wind turbine rotor speed. Hence, the thrust force and torque can be expressed as ³⁷

$$F_T=\frac{1}{2}\rho \pi R^2{C}_T(\lambda ,\beta )U^2$$

(5)

$$F_Q=\frac{1}{2}\rho \pi R^3{C}_Q(\lambda ,\beta )U^2$$

(6)

where $\rho$ is the air density. Also, the wind turbine power can be expressed similarly using the non-dimensional power coefficient $C_P$ ³⁷

$$P=\frac{1}{2}\rho \pi R^2{C}_P(\lambda ,\beta )U^3$$

(7)

Subspace system identification

To identify the resulting LTI system, there are two main system identification techniques, namely, prediction error methods and subspace identification methods. In prediction error methods, optimization methods are used to find the model parameters by minimizing the prediction error of the model. Alternatively, the subspace system identification methods use linear algebra techniques to determine the model parameters. The subspace system identification methods have some advantages over prediction error identification methods in practical cases. ³⁸ For example, they are fast and they avoid the difficulties of finding a model parametrization and performing a non-convex gradient based optimization, which may lead to numerical issues and local minima.^39,40 In this work the subspace system identification method is implemented based on Overschee and De Moor’s research ⁴¹ and will be applied to identify the dynamic model of the wind turbine tower fore-aft motion.

The notation and mathematical techniques required for the rest of this section are introduced in Appendix 3. The starting point to perform the state-space subspace system identification is to combine the recursive state-space innovation model into a single linear matrix equation by recursively substituting equation (2) into equation (3) (see Van Overschee and De Moor ⁴¹ and Qin ⁴² for details)

$${\mathbf{Y}}_f={\mathbf{\Gamma }}_h{\mathbf{X}}_f+{\mathbf{H}}_h^d{\mathbf{U}}_f+{\mathbf{H}}_h^s{\mathbf{E}}_f$$

(8)

where ${\mathbf{U}}_f$ , ${\mathbf{Y}}_f$ , ${\mathbf{X}}_f,{\mathbf{\Gamma }}_h,$ and ${\mathbf{E}}_f$ are the input, output, states sequence, extended observability, and noise block Hankel matrices, respectively, which are introduced in Appendix 3 and can be arranged in the same way as shown in equation (53). ${\mathbf{H}}_h^s$ and ${\mathbf{H}}_h^d$ are block-Toeplitz matrices defined in Appendix 3, as well. The first step of system identification is to determine ${\mathbf{\Gamma }}_h{\mathbf{X}}_f$ by eliminating ${\mathbf{H}}_h^d{\mathbf{U}}_f$ and ${\mathbf{H}}_h^s{\mathbf{E}}_f$ from the right hand side of equation (8), which will be done using oblique projection of subspaces. It can be proven that ${\mathbf{Y}}_f$ lies in the joint row space of ${\mathbf{U}}_f$ and ${\mathbf{W}}_p$ and also that ${\mathbf{E}}_f$ is perpendicular to it. ⁴¹ With these two assumptions and using equation (8), the oblique projection of ${\mathbf{Y}}_f$ along ${\mathbf{U}}_f$ on ${\mathbf{W}}_p$ will give ${\mathbf{\Gamma }}_h{\mathbf{X}}_f$

$${\mathbf{O}}_h={\mathbf{Y}}_f\underset{{\mathbf{U}}_f}{/}{\mathbf{W}}_p={\mathbf{\Gamma }}_h{\widehat{\mathbf{X}}}_f$$

(9)

$${\mathbf{O}}_{h+1}={\mathbf{Y}}_f^{-}\underset{{\mathbf{U}}_f^{-}}{/}{\mathbf{W}}_p^{+}={\mathbf{\Gamma }}_{h-1}{\widehat{\mathbf{X}}}_f^{+}$$

(10)

where hat ( $\widehat{}$ ) indicates the estimate. In the same way, ${\mathbf{\Gamma }}_{h-1}{\widehat{\mathbf{X}}}_f^{+}$ can be determined by oblique projection of ${\mathbf{Y}}_f$ along ${\mathbf{U}}_f$ on ${\mathbf{W}}_p$ as shown in equation (10). The second step is to use singular value decomposition (SVD) to split ${\mathbf{O}}_h$ into the observability matrix ${\mathbf{\Gamma }}_h$ and the state sequence matrix ${\mathbf{X}}_f$ . In this process, the order of the model $n$ can also be determined by separating the diagonal singular value matrix $\mathbf{S}$ into two parts in which ${\mathbf{S}}_1$ is a diagonal $n\times n$ matrix and contains the significant singular values and ${\mathbf{S}}_2$ contains the rest of the singular values

$${\mathbf{O}}_h=\mathbf{US}{\mathbf{V}}^T=\left(\begin{array}{ll}{\mathbf{U}}_1& {\mathbf{U}}_2\end{array}\right)\left(\begin{array}{ll}{\mathbf{S}}_1& 0\\ 0& {\mathbf{S}}_2\end{array}\right)\left(\begin{array}{l}{\mathbf{V}}_1^T\\ {\mathbf{V}}_2^T\end{array}\right)$$

(11)

$${\mathbf{\Gamma }}_h={\mathbf{U}}_1{\mathbf{S}}_1^{1/2}$$

(12)

$${\mathbf{X}}_f={\mathbf{S}}_1^{1/2}{\mathbf{V}}_1^T$$

(13)

Now, we obtain the ${\mathbf{\Gamma }}_{h-1}$ by removing the last block row from ${\mathbf{\Gamma }}_h$ . Then, from equation (10) using the least squares method, matrix ${\widehat{\mathbf{X}}}_f^{+}$ can be obtained as

$${\widehat{\mathbf{X}}}_f^{+}={\mathbf{\Gamma }}_{h-1}^{†}\mathrm{ }{\mathbf{O}}_{h+1}$$

(14)

where ${}^{†}$ indicates the pseudo-inverse. The last step is to estimate the state-space model parameters $\mathbf{A}$ , $\mathbf{B}$ , $\mathbf{C}$ and $\mathbf{D}$ by solving the following equation

$$\left(\begin{array}{l}{\widehat{\mathbf{X}}}_{i+1}\\ {\mathbf{Y}}_i\end{array}\right)=\left(\begin{array}{ll}\mathbf{A}& \mathbf{B}\\ \mathbf{C}& \mathbf{D}\end{array}\right)\left(\begin{array}{l}{\widehat{\mathbf{X}}}_i\\ {\mathbf{U}}_i\end{array}\right)+\left(\begin{array}{l}{\mathbf{\rho }}_w\\ {\mathbf{\rho }}_v\end{array}\right)$$

(15)

where ${\mathbf{\rho }}_w$ and ${\mathbf{\rho }}_v$ are state and output residual vectors, respectively. The state sequences ${\widehat{\mathbf{X}}}_i={\widehat{\mathbf{X}}}_f$ and ${\widehat{\mathbf{X}}}_{i+1}={\widehat{\mathbf{X}}}_f^{+}$ (shown in equations (59) and (60)) are estimated from the previous steps of the subspace system identification and the inputs ${\mathbf{U}}_i$ and outputs ${\mathbf{Y}}_i$ are available as given below

$${\mathbf{U}}_i=\left[\begin{array}{llll}{\mathbf{u}}_i& {\mathbf{u}}_{i+1}& \dots & {\mathbf{u}}_{i+j-1}\end{array}\right]$$

(16)

$${\mathbf{Y}}_i=\left[\begin{array}{llll}{\mathbf{y}}_i& {\mathbf{y}}_{i+1}& \dots & {\mathbf{y}}_{i+j-1}\end{array}\right]$$

(17)

The least squares solution of equation (15) will give us the state space model parameters

$$\left(\begin{array}{ll}\mathbf{A}& \mathbf{B}\\ \mathbf{C}& \mathbf{D}\end{array}\right)=\left(\begin{array}{l}{\widehat{\mathbf{X}}}_{i+1}\\ {\mathbf{Y}}_i\end{array}\right){\left(\begin{array}{l}{\widehat{\mathbf{X}}}_i\\ {\mathbf{U}}_i\end{array}\right)}^{†}$$

(18)

The input and output signals used for wind turbine tower system identification are then

$${\mathbf{y}}_k=D_{FA\mathrm{ }k}$$

(19)

$${\mathbf{u}}_k=F_{T\mathrm{ }k}$$

(20)

where $D_{FA\mathrm{ }k}$ and $F_{T\mathrm{ }k}$ are wind turbine tower tip fore-aft deflection and aerodynamic thrust force, respectively. It should be mentioned that output vector ${\mathbf{y}}_k$ could include tower deflections or accelerations at several points. Using the knowledge of wind speed $U$ , rotor speed $\mathrm{\Omega }$ , and pitch angle $\beta$ , the aerodynamic thrust force can be calculated using the BEM as given by equation (5).

In this paper, the $WP\mathrm{ }1.5MW$ wind turbine whose simulation data are provided by FAST ⁴³ is studied. For this turbine, the first four diagonal elements of the matrix $\mathbf{S}$ in equation (11) are significant and we can see a big drop in the value of the fifth diagonal element compared to the fourth element as shown in Figure 2. Therefore, the system has four significant singular values and the wind turbine tower system should be a second order model where ${\mathbf{S}}_1$ is a $4\times 4$ matrix.

Figure 2.

Diagonal elements of $S$ matrix SVD of $O_h$ matrix.

To collect the input-output data in the desired operating condition an adaptive wind turbine torque and pitch controller are designed and implemented.

It should be noted that using open-loop subspace system identification methods to identify systems that are operating under closed-looped control may lead to a biased estimation of model parameters. The reason for this is that the formulation used for open-loop system identification assumes that the input signal ${\mathbf{u}}_k$ to the system is uncorrelated with the past noise process ${\mathbf{e}}_k$ . Using closed-loop data for system identification violates the assumption, which means

$$\mathbb{E}\left[{\mathbf{u}}_k{\mathbf{e}}_j^T\right]\ne 0\mathrm{ }\mathrm{for}\mathrm{ }j<k$$

(21)

This in turn results in an input signal that is correlated with the process noise. Specifically, in subspace system identification algorithms, the past input-output data is used as the regressor to eliminate the noise in the matrix equation assuming that the past input-output data is uncorrelated with the noise.

However, despite the fact that wind turbines operate under closed-loop control, we still can use subspace system identification techniques to identify the model parameters of the wind turbine tower. The tower dynamics does not have any feedback loop. Therefore, considering thrust force $F_T$ as the input and tower fore-aft deflection $D_{FA}$ as the output, we can employ open loop system identification methods to identify the wind turbine tower dynamic model.

Thrust force estimation

In Kalman filter design, presented in Appendix 4, ⁴⁴ the innovation residual ${\tilde{\mathbf{y}}}_k$ is defined as the difference between the measured output and estimated output given by equation (65). We are going to use the estimated innovation residual sequence with and without the effect of the input force to estimate the unknown input ${\mathbf{F}}_k$ .

The forced innovation residual is defined as ${\widehat{\mathbf{z}}}_k={\tilde{\mathbf{y}}}_k$ . Using equations (65) and (63), we have

$${\widehat{\mathbf{z}}}_k={\mathbf{z}}_k-\mathbf{H}{\widehat{\mathbf{x}}}_{k|k-1}={\mathbf{z}}_k-\mathbf{HA}{\widehat{\mathbf{x}}}_{k-1|k-1}-\mathbf{HB}{\mathbf{F}}_{k-1}$$

(22)

where ${\widehat{\mathbf{x}}}_{k-1|k-1}$ is the Kalman filter estimation of the system states, which can be rewritten by combining equations (64), (65) and (63), as follows

$${\widehat{\mathbf{x}}}_{k|k}=(\mathbf{I}-{\mathbf{K}}_k{\mathbf{H}}_k)(\mathbf{A}{\widehat{\mathbf{x}}}_{k-1|k-1}+\mathbf{B}{\mathbf{F}}_{k-1})+{\mathbf{K}}_k{\mathbf{z}}_k$$

(23)

The innovation residual without the effect of the input ${\overline{\mathrm{z}}}_k$ can be defined by removing the effect of input ${\mathbf{F}}_{k-1}$ from equation (22)

$${\overline{\mathbf{z}}}_k={\mathbf{z}}_k-\mathbf{H}{\overline{\mathbf{x}}}_{k|k-1}={\mathbf{z}}_k-\mathbf{HA}{\overline{\mathbf{x}}}_{k-1|k-1}$$

(24)

where ${\overline{\mathbf{x}}}_{k|k}$ is the estimation of the system states without the effect of the input ${\mathbf{F}}_{k-1}$ , which can be defined by removing the input ${\mathbf{F}}_{k-1}$ from equation (23) as follows

$${\overline{\mathbf{x}}}_{k|k}=(\mathbf{I}-{\mathbf{K}}_k{\mathbf{H}}_k)\mathbf{A}{\overline{\mathbf{x}}}_{k-1|k-1}+{\mathbf{K}}_k{\mathbf{z}}_k$$

(25)

The unknown input can be estimated by defining $\mathrm{\Delta }{\mathbf{z}}_k={\overline{\mathbf{z}}}_k-{\widehat{\mathbf{z}}}_k$ which is the difference between the innovation residuals with and without the effect of the input. Therefore, by subtracting equation (22) from equation (24) and defining $\mathrm{\Delta }{\mathbf{x}}_k={\widehat{\mathbf{x}}}_{k|k}-{\overline{\mathbf{x}}}_{k|k}$ we have:

$${\overline{\mathbf{z}}}_k-{\widehat{\mathbf{z}}}_k=\mathbf{HA}\mathrm{\Delta }{\mathbf{x}}_k+\mathbf{HB}{\mathbf{F}}_{k-1}$$

(26)

By subtracting equation (23) from equation (25) and assuming the input to be zero before the time step $k=n$ ( ${\mathbf{F}}_k=0$ for $k<n$ ) and to be constant from $k=n$ to $k=n+l$ ( ${\mathbf{F}}_k=\mathit{Constant}$ over interval $k=n,n+1,\dots ,n+l$ ) we can express $\mathrm{\Delta }{\mathbf{x}}_k={\mathbf{M}}_k\mathbf{B}{\mathbf{F}}_{k-1}$ where ${\mathbf{M}}_k$ can be determined from the following recursive relation

$${\mathbf{M}}_k=\{\begin{array}{ll}(\mathbf{I}-{\mathbf{K}}_k{\mathbf{H}}_k)(\mathbf{A}{\mathbf{M}}_{k-1}+\mathbf{I})& n<k⩽n+l\\ 0& k<n\end{array}$$

(27)

Finally, substituting $\mathrm{\Delta }{\mathbf{x}}_k={\mathbf{M}}_k\mathbf{B}{\mathbf{F}}_{k-1}$ into equation (26) leads to the following equation that will be used to estimate unknown input ${\mathbf{F}}_k$

$${\overline{\mathbf{z}}}_k=\{\begin{array}{ll}{\widehat{\mathbf{z}}}_k& k⩽n\\ {\widehat{\mathbf{z}}}_k+{\mathbf{\Lambda }}_k{\mathbf{F}}_{k-1}& n<k⩽n+l\end{array}$$

(28)

where

$${\mathbf{\Lambda }}_k=\mathbf{H}(\mathbf{A}{\mathbf{M}}_{k-1}+\mathbf{I})\mathbf{B}$$

(29)

For input ${\mathbf{F}}_k$ , estimation equation (28) can be arranged in the following matrix form by assuming ${\mathbf{F}}_k$ to be constant over interval $n<k⩽n+l$

$${\mathbf{Y}}_N={\mathbf{\Phi }}_N\mathbf{F}+{\mathbf{\epsilon }}_N$$

(30)

where

$${\mathbf{Y}}_N={\left(\begin{array}{llll}{\overline{\mathbf{z}}}_{n+1}& {\overline{\mathbf{z}}}_{n+2}& \dots & {\overline{\mathbf{z}}}_{n+l}\end{array}\right)}^T$$

(31)

$${\mathbf{\Phi }}_N=\left(\begin{array}{c}\mathbf{H}\mathbf{B}\\ \mathbf{H}(\mathbf{A}{\mathbf{M}}_{n+1}+\mathbf{I})\mathbf{B}\\ ⋮\\ \mathbf{H}(\mathbf{A}{\mathbf{M}}_{n+l-1}+\mathbf{I})\mathbf{B}\end{array}\right)$$

(32)

$${\mathbf{\epsilon }}_N={\left(\begin{array}{llll}{\widehat{\mathbf{z}}}_{n+1}& {\widehat{\mathbf{z}}}_{n+2}& \dots & {\widehat{\mathbf{z}}}_{n+l}\end{array}\right)}^T$$

(33)

In equation (30), ${\mathbf{\epsilon }}_N$ can be treated as a disturbance with the following covariance matrix

$${\mathbf{\Sigma }}_N=\left(\begin{array}{llll}{\mathbf{s}}_{n+1}& 0& \dots & 0\\ 0& {\mathbf{s}}_{n+1}& \dots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & {\mathbf{s}}_{n+l}\end{array}\right)$$

(34)

where ${\mathbf{s}}_n$ is the covariance of ${\widehat{\mathbf{z}}}_k={\tilde{\mathbf{y}}}_k$ and can be determined from equation (68). It should be noted that ${\mathbf{\Sigma }}_N$ is a positive definite matrix because its diagonal elements defined by equation (68) are all positive.

The unknown thrust force $\mathbf{F}$ can be estimated by solving equation (30). By applying the least squares method ⁴⁵ to equation (30) and using ${{\mathbf{\Sigma }}_N}^{-1}$ as the weighting matrix – inspired by the fact that ${\mathbf{\Sigma }}_N$ is a positive definite matrix – we can estimate the unknown thrust force $\mathbf{F}$ using $N$ measurements as follows

$$\widehat{\mathbf{F}}=[{{\mathbf{\Phi }}_N}^T{{\mathbf{\Sigma }}_N}^{-1}{\mathbf{\Phi }}_N{]}^{-1}{{\mathbf{\Phi }}_N}^T{{\mathbf{\Sigma }}_N}^{-1}{\mathbf{Y}}_N$$

(35)

where ${{\mathbf{\Sigma }}_N}^{-1}$ is defined as equation (36) and $0<\gamma ⩽1$ is the fading factor, which causes the recent measurements to be weighted more heavily than the less recent ones.

$${{\mathbf{\Sigma }}_N}^{-1}=\left(\begin{array}{cccc}{{\mathbf{s}}_{n+1}}^{-1}{\gamma }^{l-1}& 0& \dots & 0\\ 0& {{\mathbf{s}}_{n+2}}^{-1}{\gamma }^{l-2}& \dots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & {{\mathbf{s}}_{n+l}}^{-1}\end{array}\right)$$

(36)

Using the least squares method to include an additional measurement ${\overline{\mathbf{z}}}_{N+1}$ and using $N+1$ measurements to estimate a new value of the thrust force $F$ , we need to repeat the procedure, which will be very expensive in terms of computational cost, especially if we want to use it for real time estimation of $\mathbf{F}$ . The recursive least squares method ⁴⁵ can alleviate this problem by utilizing the previous estimate of ${\widehat{\mathbf{F}}}_N$ to estimate the new ${\widehat{\mathbf{F}}}_{N+1}$ as follows

$${\widehat{\mathbf{F}}}_{N+1}={\widehat{\mathbf{F}}}_N+{{\mathbf{K}}_{\mathit{RLS}}}_{N+1}[{\overline{\mathbf{z}}}_{N+1}-{\mathbf{\Lambda }}_{N+1}{\widehat{\mathbf{F}}}_N]$$

(37)

where ${{\mathbf{K}}_{\mathit{RLS}}}_{N+1}$ and ${{\mathbf{P}}_{\mathit{RLS}}}_{N+1}$ are defined^52,42 as

$${{\mathbf{K}}_{\mathit{RLS}}}_{N+1}={\gamma }^{-1}{{\mathbf{P}}_{\mathit{RLS}}}_N{{\mathbf{\Lambda }}_{N+1}}^T[{\mathbf{\Lambda }}_{N+1}{\gamma }^{-1}{{\mathbf{P}}_{\mathit{RLS}}}_N{{\mathbf{\Lambda }}_{N+1}}^T+{\mathbf{s}}_{N+1}{]}^{-1}$$

(38)

$${{\mathbf{P}}_{\mathit{RLS}}}_{N+1}=[\mathbf{I}-{{\mathbf{K}}_{\mathit{RLS}}}_{N+1}{\mathbf{\Lambda }}_{N+1}]{\gamma }^{-1}{{\mathbf{P}}_{\mathit{RLS}}}_N$$

(39)

Then, we can use equation (37) for recursive estimation of the thrust force using the wind turbine tower deflection measurement. The tower deflection measurement can be practically performed by using an array of low-cost MEMS accelerometers and IMUs installed inside and along the tower of the wind turbine.^28–30

Effective wind speed estimation

From equation (5), given the estimated aerodynamic thrust force, the blade pitch angle $\beta$ and the rotor speed $\mathrm{\Omega }$ , we can solve the inverse thrust force function to estimate the wind speed.

If we train an ANN with the input vector composed of thrust force, rotor speed and pitch angle $[F_T,\mathrm{\Omega },\beta ]$ and the output composed of wind speed $U$ , we can invert equation (5) to estimate the wind speed from these three parameters. However, some manipulations in equation (5) can improve the performance of wind speed estimation. By substituting wind speed $U$ from equation (4), we can rewrite equation (5) as

$$\frac{2F_T}{\rho \pi R^4{\mathrm{\Omega }}^2}=\frac{C_T(\beta ,\lambda )}{{\lambda }^2}$$

(40)

The left hand side of equation (40) can be calculated knowing the thrust force $F_T$ and rotor speed $\mathrm{\Omega }$ . We can use the vector $[\frac{2F_T}{\rho \varphi R^4\mathrm{\Omega }},\beta ]$ as the input and tip-speed ratio $\lambda$ as the output to train the ANN. Using this approach, we can calculate the estimated tip-speed ratio $\widehat{\lambda }$ from the ANN, which can be used to estimate the effective wind speed $\widehat{U}=\frac{R\mathrm{\Omega }}{\widehat{\lambda }}$ .

Solving equation (40) for the wind speed estimate has some advantages over solving equation (5).⁷ Using equation (40), we use two inputs for training the ANN while using equation (5) we use three inputs for training the ANN. Having a smaller number of inputs to the ANN requires fewer number of hidden layers in training the ANN. Also, having the same number of hidden layers, the ANN with fewer inputs has better computational performance. Moreover, in equation (5), we have the air density inside the mapping function. Therefore, if the air density changes from the situation in which we collected the training data, we should train the neural network again. Using equation (40), the air density is inside the input vector and the neural network is independent of the air density. It should be noted that the air density varies due to changes in height and temperature which can cause significant change in wind turbine power generation. ⁴⁶

In this work, the MATLAB neural network toolbox ⁴⁷ is used to train a MLP neural network using the backpropagation technique. The trained neural network has ten neurons in its only hidden layer, two inputs $[\frac{2F_T}{\rho \pi R^4\mathrm{\Omega }},\beta ]$ and one output $\lambda$ where the effective wind speed is determined as $\widehat{U}=\frac{R\mathrm{\Omega }}{\widehat{\lambda }}$ . The structure of MLP is shown Figure 3.

Figure 3.

Neural network configuration. ⁴⁸

Results and discussion

In this section, we will use the developed wind speed estimation algorithm to estimate the effective wind speed in different wind conditions. Some numerical results for effective wind speed estimation will be presented and the performance of the estimator will be tested and evaluated.

Wind speed estimation is performed for the $WP\mathrm{ }1.5MW$ wind turbine using the simulation data provided by FAST, ⁴³ which is a certified wind turbine simulation package developed by National Renewable Energy Laboratory.

First, we use a uniform time-varying wind profile, which is spatially constant in the rotor plane and only varies by time, to evaluate the estimator. Then the estimator will be used to estimate the effective wind speed in fully turbulent wind fields where the wind profile varies temporally and spatially over the rotor plane. In addition to turbulence, wind speed mean value increases as we go further from the ground as a result of wind shear effect.

For a uniform wind profile, since the wind speed varies uniformly over the rotor plane, we can compare the estimated wind speed with the input wind speed to verify the estimation result. Figure 4 shows the performance of the proposed estimator by comparing the estimated wind speed and aerodynamic thrust force to the input wind speed profile and the thrust force from the FAST environment. The wind profile used for evaluation in this case is simulated to be spatially constant in the rotor plane and only varies in time.

Figure 4.

(a) Thrust force estimator validation; (b) wind speed estimator validation.

For a turbulent wind profile the effective wind speed is expected to be a weighted average of the wind speeds at different points of the wind field where the weight of each point is determined based on the aerodynamic performance of the wind turbine blade section swiping that point. Hence, there is no reference wind speed with which to compare the estimated effective wind speed and evaluate the results.

Therefore, a verification method based on a correlation coefficient is proposed to evaluate the estimated effective wind speed for turbulent wind profiles. ⁴⁸ Considering ${\sigma }_A$ and $\mathit{Cov}\left(A,B\right)$ to be the standard deviation of $A$ and the covariance of $A$ and $B$ respectively, the correlation coefficient is defined by

$$\rho \left(A,B\right)=\frac{\mathit{Cov}\left(A,B\right)}{{\sigma }_A{\sigma }_B}$$

(41)

Its value varies from $-1$ to $1$ and higher absolute values of the correlation coefficient show higher dependency between the signals.

To evaluate the wind speed estimation method, we determine the correlation coefficient of the estimated effective wind speed with wind speeds at different points in the rotor plane obtained from the FAST simulation environment. We expect to see higher values of the correlation coefficient inside the rotor plane.

In the simulations, the strength of turbulence in a wind field is measured by Turbulence Intensity (TI) defined as follows ⁴⁹

$$TI=\frac{{\sigma }_V}{\overline{V}}\times 100\mathrm{\%}$$

(42)

where ${\sigma }_V$ is the standard deviation of the wind field and $\overline{V}$ is the mean wind speed.

In the rest of this section, we use our wind speed estimation method to estimate effective wind speed for two different turbulent wind fields. Simulation conditions for the wind fields are presented in Table 1. To better visualize the applied wind speed distribution, in Figure 5, two snapshots of a turbulent wind field are shown illustrating that wind speed is time varying and has a nonuniform distribution over a field. The circles in Figures 5 and 6 show the rotor plane of the wind turbine. Furthermore, the variation of the mean value of the wind speed over a wind field is presented in Figure 6, which shows that the mean value of the wind speed increases as we go higher from the ground. This variation in wind speed mean value is a result of a wind shear effect that causes the wind speed to change with distance from the ground. For the wind field presented in Figures 5 and 6, $TI=16\mathrm{\%}$ , which is a relatively high turbulence intensity.

Table 1.

Wind field and simulation conditions. ⁴⁸

Test #	1	2
$TI$ (%)	$16\mathrm{\%}$	$12\mathrm{\%}$
Mean wind speed (m/s)	$18$	$25$
Wind field size (m)	$120\times 120$	$120\times 120$
Rotor plane diameter (m)	$70\mathrm{ }$	$70\mathrm{ }$
Wind sampling time step (s)	$0.01\mathrm{ }$	$5$
Simulation period (s)	$1000\mathrm{ }$	$1000\mathrm{ }$
Simulation time step (s)	$0.01\mathrm{ }$	$0.01\mathrm{ }$

Figure 5.

Spatial distribution of wind speed in a turbulent wind field: (a) $t=50\mathrm{ }\mathrm{s};$ (b) $t=65\mathrm{ }\mathrm{s}$ .

Figure 6.

Spatial distribution of the mean value of the wind speed in the simulated turbulent wind field.

Next, in order to verify the effective thrust force estimation results, we have compared the estimated values for thrust force with thrust force output of FAST. The results of comparisons for test #1 and test #2 are presented in Figure 7. We can see that test #1 shows more fluctuations in both estimated thrust force and FAST output. The reason is that for test #1 we have a higher value of $TI$ and higher wind speed sampling frequency, which together cause more fluctuations in thrust force. The mean absolute error in the estimated effective thrust force is presented in Table 2 showing that estimation error for test #1 is higher than test #2. The reason is that the wind speed profile in test #1 has a greater amount of high frequency content falling beyond the bandwidth of tower dynamics. This shows that the thrust force estimator has good performance even for highly turbulent wind conditions and its performance gets better for wind fields with lower turbulence intensity.

Figure 7.

Comparison of estimated values of effective thrust force and FAST thrust force output for: (a) test #1; (b) test #2.

Table 2.

Mean absolute value of effective thrust force estimation error.

Test #	1	2
Thrust force estimation error (mean)	$13.8\mathrm{\%}$	$7.3\mathrm{\%}$

The estimated effective wind speed is compared to hub-height wind speed and the results are presented in Figure 8. The mean value of the effective wind speed is expected to be close to the mean value of the wind speeds over the rotor plane. For the simulated cases, these two values are also close to the mean value of hub-height wind speeds as shown in Table 3.

Figure 8.

Comparison of estimated values of effective wind speed and hub-height wind speed for: (a) test #1; (b) test #2.

Table 3.

Comparison of wind speeds mean values.

Test #	1	2
Effective wind speed (m/s)	$18.54$	$24.7\mathrm{ }$
Wind speed mean value over rotor plane (m/s)	$18.7$	$24.4$
Hub-height wind speed m/s)	$18.36$	$25.31$

Moreover, in Figure 9, the correlation of effective wind speeds and simulated wind speeds over the wind field is presented. The color bar shows the dimensionless values of the correlation coefficient. It is observed that correlation coefficients are higher for test #2 where we have lower turbulence intensity. Also, from Figure 9, we have a higher value of correlation coefficient for wind speeds inside the rotor plane.

Figure 9.

Correlation coefficient of the turbulent wind field and the effective wind speed for: (a) test $\#1;$ (b) test $\#2$ .

Based on the definition of effective wind speed, which is the spatial average of wind speeds over the rotor plane, we can conclude that the presented results meet the main purpose of this research, which is the development of a tower deflection-based effective wind speed estimation method.

As mentioned before, the wind speed estimation in this study has two stages. In the second stage, an MLP was used. In order to investigate the performance of other types of machine learning tools, the estimation errors obtained using a radial basis function (RBF) and a support vector machine (SVM) are compared to MLP results. In order to limit the errors to the error of the second stage, the real thrust force and pitch angle values were given as inputs to these three machine learning tools, and the estimated wind speed is then compared to the hub-height wind speed. In Table 4, the mean squared error (MSE) and the mean absolute error (MAE) are reported.

Table 4.

Mean squared and mean absolute value of effective wind speed estimation errors of three machine learning algorithms.

Machine learning algorithm	MSE ${\left(\mathrm{m}/\mathrm{s}\right)}^2$	MAE $\left(\mathrm{m}/\mathrm{s}\right)$
MLP	0.0370	0.1312
RBF	0.0375	0.1351
SVM	0.2090	0.1674

MLP: multi-layer perceptron; RBF: radial basis function; SVM: support vector machine.

The results in Table 4 show that the estimation errors of these methods are close to each other and low compared to the speed values. Therefore, the final speed estimation accuracy is not significantly affected by selection of the machine learning algorithm.

Finally, it should be mentioned that effective wind speed has some components in addition to the original upstream wind speed, which are created due to the contribution of the wind turbine itself. For example, wind turbine blades experience different wind speed values while sweeping the rotor plane, causing fluctuations in effective wind speed with the rotational frequency of the wind turbine $P$ and more significantly its multiples such as $3P$ (for three-bladed wind turbines). Also, wind turbine blade vibration contributes to effective wind speed by changing the relative wind speed. To investigate these effects, we can study the frequency domain representation of the effective wind speed.

The frequency domain comparison presented in Figure 10 shows that the effective wind speed and hub-height wind speed have similar frequency content at low frequency range. However, the effective wind speed for both test #1 and test #2 has some higher frequency components that are not present in the hub–height wind speed. Furthermore, the additional components have the same frequency of $f_1\simeq 0.4\hspace{0.17em}\mathrm{Hz}$ , $f_2\simeq 0.6\hspace{0.17em}\mathrm{Hz}$ and $f_3\simeq 1.2\hspace{0.17em}\mathrm{Hz}$ for test $\#1$ and test $\#2$ .

Figure 10.

Frequency domain representation of: (a) hub-height wind speed, test #1; (b) effective wind speed, test #1; (c) hub-height wind speed, test #2; (d) effective wind speed, test #2.

Figure 11.

Rotor speed for test $\#2$ which is maintained around $\mathrm{\Omega }=12\mathrm{ }\mathrm{r}/\min$ using pitch control system.

Modal analysis of the wind turbine structural dynamics shows that these additional high frequency components correspond to the natural frequencies of the wind turbine blades $f\simeq 1.2\hspace{0.17em}\mathrm{Hz}$ and tower $f\simeq 0.4\hspace{0.17em}\mathrm{Hz}$ . Furthermore, the rotor speed is maintained at $\mathrm{\Omega }=12\mathrm{rpm}$ in our simulations using the designed pitch control system (see Figure 11). Consequently, the rotational frequency is $P=0.2\hspace{0.17em}\mathrm{Hz}$ showing that the frequency peak of $f=0.6\hspace{0.17em}\mathrm{Hz}$ is associated with $3P$ fluctuations caused by turbulence, wind shear and tower shadow.

In summary, the presented results in this section show that the proposed tower deflection-based wind speed estimation method can predict the effective thrust force and wind speed successfully for turbulent wind fields.

Conclusions

This paper has presented a method for estimation of the effective wind speed based on wind turbine tower deflection measurement. The wind speed estimator consists of two main components: A Kalman filter along with a recursive least squares estimator to estimate the thrust force from deflection measurement, and an MLP neural network for solving the inverse problem of effective wind speed estimation from thrust force. The Kalman filter was designed based on a data-driven state space model of wind turbine tower dynamics that is identified using the subspace system identification technique. It was used to estimate the aerodynamic thrust force, which was then used to estimate the effective wind speed for turbulent wind fields. The results of the proposed wind speed estimator showed good agreement between the estimated thrust force and the thrust force obtained from the FAST simulation environment. The results also showed a high correlation between the wind speed close to the hub position and the estimated wind speed. The speed estimation error caused by the applied MLP was compared with the errors of an RBF and an SVM estimators, which showed no significant differece. Comparisons in the frequency domain were also promising for low-frequency components of the wind speeds. The results showed that this proposed method for estimating the effective wind speed can be an effective and successful alternative to other available methods.

An optimal control system might need to be designed to further exploit the faster estimation of effective wind speed. It is expected from this effective wind speed estimation that pitch control systems will benefit from using this fast estimation algorithm. Testing this will require designing a pitch control system that uses wind speed estimation in a feed-forward setting. Further investigation seems to be necessary though in order to study the performace of the algorithm in the case of rapidly changing wind speed mean value.