Augmented Kalman filter with a reduced mechanical model to estimate tower loads on an onshore wind turbine: a digital twin concept

The paper presents an application of the Kalman filtering technique to estimate loads on a wind turbine. The approach combines a mechanical model and a set of measurements to estimate signals that are not available in the measurements, such as: wind speed, thrust, tower position, and tower loads. The model is several fold faster than real-time and is intended to be run online, for instance, to evaluate real-time fatigue life consumption of a field turbine using a digital twin. The mechanical model is built using a Rayleigh-Ritz approach and a set of joint coordinates. The paper presents a general method and illustrates it 5 using a 2 degrees of freedom model of a wind turbine, and, using rotor speed, generator torque, pitch, and tower-top acceleration as measurement signals. The different components of the model are tested individually. The overall method is evaluated by computing the errors in estimated tower bottom equivalent moment from a set of simulations. From this preliminary study, it appears that the tower bottom equivalent moment is obtained with about 10% accuracy. The limitation of the model and the required steps forwards are discussed. 10 DOF degrees of freedom


Introduction
Wind turbines are designed and optimized for a given site using numerical tools, and, a statistical assessments of the environmental conditions the turbine will experience. The uncertainty on the tools and data are accounted for using multiplicative safety factors, which are determined from a combination of experience and specifications by the standards. Overconservative The first part presents the different components required for this work: the augmented KF, the numerical model of the turbine, and the estimators for the wind speed, thrust, tower load, and fatigue. Simple illustrations and validation results for the different components of the model are provided in a second part. The third part presents full applications but limited to simulations. Discussions and conclusions follow.
2 Description of the models 60 2.1 Example for a 2DOF wind turbine model We start this section by an illustrative example, before describing the different parts of the model in their general form. A wind turbine is here modelled using 2 DOF: 1) the generalized coordinate associated with the fore-aft (FA) bending of the tower, q t ; 2) the shaft rotation, ψ. The tower bending is associated with a shape function Φ t (z), such that the FA displacement of a point, at height z, and at time t, is given by u(z, t) = q t (t)Φ t (z). The shape function is normalized to unity at the tower-top, and q t is 65 then equal to the FA displacement at the tower-top (see Figure 1). The equations of motion of the system are:  where: M , C, K are the generalized mass, damping and stiffness associated with the FA DOF; J is the drivetrain inertia; T * a and Q a are the aerodynamic thrust and torque; and Q g is the generator torque. A star is used as upper-script of the thrust to indicate that using the thrust directly is a rough approximation. A more elaborate expression of the generalized force acting 70 on q t is given in subsection 3.3. The determination of M , C and K is discussed in Branlard (2019). For the NREL-5MW turbine, the values are: M = 4.4e 5 kg, D = 2.5e 4 kg/s, K = 2.7e 6 kg/s 2 and J = 4.3e 7 kg.m 2 . In this example, the system of equations is only coupled via the aerodynamic loads.
The following measurements are usually readily available on any commercial wind turbine: the generator power, P g ; the blade pitch angle, θ p ; the rotor rotational speed, Ω =ψ. and the tower-top acceleration in the FA direction,q t ; The knowledge 75 of the generator power, speed and losses allows to estimate the generator torque Q g . In this study, the generator torque is assumed known. We will use an augmented KF concept to combine these measurements with the mechanical model to estimate the state of the system. The KF algorithm requires linear state and output equations. The state vector is assumed to be x = [q t , ψ,q t ,ψ, Q a ]. The fact that some of the loads were included into the state vector is referred to as state augmentation. The choice of loads to include in the state vector is not unique and will lead to different state equations. Using this choice for x, 80 Equation 1 is written into the following state equation:  where for simplicity the time derivatives of the aerodynamic torque is assumed to be zero, an assumption referred to as "random walk" force model. This assumption accounts of saying that the estimate of the torque at the next time step is likely to be close to the one at the current time step. Improvements on this will be discussed in section 5. The thrust is determined based on 85 the rotor speed, the aerodynamic torque, and the pitch angle, using tabulated data, as described in subsection 2.4. The output equation relates the measurements to the states and inputs as follows:  Equation 2 and Equation 3 are used within a KF algorithm to estimate the states vector based on the measurements. The estimated time series of q t , together with its associated shape function Φ t , are used to determine the bending moments within 90 the tower and estimate tower fatigue loads, based on the method presented in subsection 2.5. Results from this simple model will be provided in section 3. The remaining paragraphs of this section generalize the approach presented.

Mechanical model of the wind turbine
The wind turbine is described using a set of degrees of freedom (DOF) that consist of joints coordinates and shape functions coordinates. The method was described in previous work (Branlard, 2019), and the source code made available online, via a li-95 brary called YAMS (Branlard). Similar approaches are for instance used in the elastic codes Flex and OpenFAST (OpenFAST).
The advantage of the method is that the system can be described with few DOF. The number of DOF is between 2 and 30 DOF whereas traditional finite element methods require in the order of one thousand DOF.
The only joint coordinate retained in the current model is the shaft azimuthal position, noted ψ. The shaft torsion, and nacelle yaw and tilt joints can be added without difficulty. The tower and blades are represented using a set of shape functions taken 100 as the first mode shapes of these components. The shape functions of the tower are assumed to be the same in the FA and side side (SS) directions, which are respectively aligned with the x and y directions (see Figure 1). The number of shape functions are noted n xt , n yt and n b for the tower FA, tower SS, and blade respectively. Writing B the number of blades, n s the number of DOF representing the shaft, the total number of DOF is: n q = n s + Bn b + n xt + n yt . The tower DOF are written q xt,i with i ∈ [1..n xt ] and q yt,i with i ∈ [1..n yt ]. Similar notations are used for the blade DOF.

105
The equation of motions are established using Lagrange's equation. The example presented in subsection 2.1 corresponds to n s = 1, n b = 0, n t,SS = 0 and n t,FA = 1. An example for a 5DOF system, with n s = 1, B = 2, n b = 1, n t,SS = 0 and n t,FA = 1, is given in Branlard (2019). In the general case, the equation of motions are described as: where M , C, K are the mass, damping, and stiffness matrices; q is the vector of DOF; and f is the vector of generalized 110 loads acting on the DOF. An inconvenience of the method is that the mass matrix is a non-linear function of the DOF. The main assumption of this work is that the non-linearities can be discarded as a first approximation. This assumption is further discussed in section 5.

Augmented Kalman filter applied to a mechanical system
A description of the standard KF can be found e.g. in the textbook of Grewal and Andrews (2014) or Zarchan and Musoff 115 (2015). The algorithm will not be detailed in this paper. The method expect a state and output equations of the following form: where x, u and y are the state, input, and measurement vectors ; X x , X y , Y x and Y u are Jacobian matrices describing the expected relationships between measurements, states and inputs ; and w x and w y are Gaussian uncorrelated noises as-120 sociated with the state-space model and measurements respectively, of which the associated covariance matrices are noted , with E[w x w t y ] = 0 and E the expected value operator. We will develop these equations in the case of a mechanical system that follows the general form of Equation 4. Specific applications will be given in section 3 and section 4. Different approaches can be used to write Equation 4 in the form of Equation 5, depending how the force vector is to be 125 treated. In a first approach, the forces can be considered to be inputs f = u, in which case Equation 4 is directly in the form of Equation 5, with x = [q,q]. This implies that we have full knowledge of the forces acting on the system at every time step, which is unlikely. In a second approach, the forces can be assumed to be part of the system noise, w x , which would lead to x = [q,q], and B = 0. This is obviously a crude approximation since the forces acting on the system are non stochastic, and, we likely have some knowledge on them. In the intermediate approach introduced by Lourens et al. (Lourens et al., 2012), 130 some of the forces are included in the system noise, and others as part of the states. The reduced set of loads that are part of the states is written p, of length n p , and the full force vector is assumed to be approximated by: f ≈ S p p, where S p is a matrix of dimension n q × n p . The reduced set of forces, p, is integrated into the state vector as: x = [q,q, p]. This process is referred to as state augmentation.
We introduce a generalized approach and assume that the forces are a combination of states, inputs and unknown noise: where the F • matrix represent the Jacobian of the force vector with respect to vector •, and w f are unknown forces that are assumed to be part of system disturbance w x . The terms F q and Fq are linearized stiffness and damping terms. These terms are zero if their contributions are already included in the definitions of K and C. In practice, the linearization of the force vector may not be possible, and assumed relationships or engineering models are used. As an example, if p contain the thrust 140 force and f the moment at the tower base, the appropriate element of F p could be set with the lever arm between tower-top and tower base.
This approach allows us to use the knowledge we have of some of the main loads acting on the system and express their dynamics into the state-space equation. The forces may for instance be assumed to follow a first order system as follows: where the P • matrices are obtained from a knowledge of the force evolution. Second order system could also be introduced, in which case the state needs to be augmented with both p andṗ ("random walk" force model). For simplicity, the applications used in this work will assumeṗ = 0, but future work will investigate the benefit of using first order systems for the evolution of the forces.

Inserting Equation 7 into Equation 4
, introducing x = [q,q, p], and using Equation 8, a state equation of the form of Equa-150 tion 5 is obtained: The measurements are assumed to be a combination of the acceleration, velocity, displacements, loads and inputs: Inserting the accelerationq from Equation 4 into Equation 10, an output equation of the form of Equation 6 is obtained, with: Equation 9 and Equation 11 form the bridge between the definition of the mechanical model and the state and output equations needed by the KF algorithm.
Equation 5 and Equation 6 are in continuous form, whereas the KF algorithm uses discrete forms. The discrete form of the matrices perform the time integration of the states from one time step to the next, namely: the subscript d indicates the discrete form of the matrices and k is the time step index. The matrix X x,d is referred to as the fundamental matrix. For time-invariant systems, this matrix may be obtained using Laplace tranform or by Taylor-series expansion (Zarchan and Musoff, 2015). For a given time step ∆t, the discrete matrices corresponding to X x and X u are: The approximation in Equation 12 is effectively a first order forward Euler time integration. The matrix Y x and Y u remain unchanged by the discretization since the ouput equation is an algebraic equation involving quantities at the same time step.
Many choices are possible as to how the model may be formulated: which forces should be accounted for in the reduced set p, which forces should be assumed to be obtained from the inputs, which models to use for the P matrices, etc. Since the study is limited to onshore wind turbines, the main loads are the aerodynamic thrust and torque. A subtlety to account for, is that 170 some of the forces of the model presented in Equation 4 are generalized forces, which are projection of loads onto the shape functions (Branlard, 2019). An example will be given in subsection 3.3.
The Jacobian matrices introduced should be determined by linearization about an operating point. The mass matrix should also be linearized about such point. In the current work, the non-linearities are either neglected, or directly inserted into the expression presented without performing a linearization. This crude simplification will be discussed in section 5, in light of the 175 results presented in section 3 and section 4.

Wind speed and thrust estimation
In this paragraph, Q a , θ p and Ω are assumed to be given. The aerodynamic power and thrust coefficients, C P and C T , are also assumed to be known as function of the pitch angle and tip speed ratio, λ = ΩR/U 0 , where R is the rotor radius and U 0 the wind speed. The functions C P (λ, θ p ) and C T (λ, θ p ) are estimated by running a parametric set of simulations at constant 180 operating conditions. Some uncertainty is here present as to whether the real turbine does performs as predicted by these functions. This question will be considered in section 5. The aerodynamic torque is computed from the tabulated data as: The wind speed is obtained by solving the following non-linear constraint equation for u est : where ρ is the air density, which is another potential source of uncertainty to be considered when dealing with measurements.
The wind speed determined by this method is assumed to be the effective wind speed acting over the rotor area. A correction for nacelle displacements is discussed in section 5. The aerodynamic thrust is estimated from this wind speed as:

Tower loads and fatigue estimation 190
The deflection of the tower, U , in the x or y directions, at a given height z, and a given time t, is given by the sum of the tower shape functions scaled by the tower degrees of freedom: The curvature, κ, is obtained by differentiating the deflection twice, giving: The bending moments along the tower heights are then obtained from the curvatures using Euler beam theory: where EI is the bending stiffness of a given tower cross section. The time series of bending moment are processed using a rain flow counting algorithm to estimate the equivalent loads and damage (WG3, 2005). devised for this case. The mean relative error for the entire time series is = 2.5%. The estimated wind speed is seen to follow the challenging trends of this time series, matching both the low and high frequencies. In the top zoom, it is seen that no phase lag is observed in the estimated wind speed, but the estimated value is seen to overshoot.

215
There are several potential sources of errors in the current methodology. One concern is whether the unsteady aerodynamic torque, can be determined using a look-up table that uses instantaneous values. The relative error between the unsteady torque and is uncorrelated to the error in torque, yet a stronger correlation is seen when the turbine is producing power. Looking at the probability density functions given in the right of Figure 2, it is seen that the errors in torque and wind speed are centered on zero. The fact that both errors are centered on zero, indicate that when the unsteady torque can indeed be obtained using instantaneous values and tabulated data, the (estimated) effective wind speed is close to the average wind speed at the rotor.

230
Other sources of errors are discussed in section 5.
A more thorough study on the questions raised above are left open for future work. Overall, the results from the test case are encouraging. Wind speed estimation is a standard feature of most wind turbine controllers, and it is likely that more advanced features are implemented by manufacturers. Any improvement on the methodology used in the current paper will be beneficial for the procedure of loads estimation presented.

Thrust estimation
The wind speed estimated in subsection 3.1 is used to estimate the thrust T a,est with Equation 15. In Figure 3  operating conditions considered. Using the estimated wind speed is seen to produce thrust values closer to the reference thrust 240 than if u ref is used. In line with the discussions of subsection 3.1, this could support the fact that the estimated wind speed provides an effective velocity, consistent with the instantaneous state of the rotor, but different from the rotor averaged wind speed. Yet, it is also possible that compensating errors are at play, or, that the thrust is less sensitive to changes of wind speed or drive-train dynamics than the torque. Despite these open questions, we continue by assuming the method provide thrust estimates with sufficient accuracy. dz (L t ). The shape functions are assumed to be normalized at their extremity, i.e. Φ t,i (L t ) = 1, so that the generalized force is:

255
The main forces acting at the tower-top are assumed to the be aerodynamic thrust and the gravitational force from the rotor nacelle assembly (RNA) mass, M RNA . The loads are then obtained as: where: θ tilt is the tilt angle of the nacelle; N R is the vector from the tower-top to the rotor center, where the thrust is assumed to act; N G is the vector from the tower-top to the RNA center of mass; g is the acceleration of gravity; and α y is the y-rotation 260 of the tower-top due to the bending of the tower (see Figure 1). For a single tower mode α y (t) = q t (t)ν 1 . The linearization of

Equation 19 and Equation 20
for small values of q t leads to: where: the term in parenthesis is the main contribution, which justifies the use of T a in Equation 1; the term in curly brackets is seen to act as a stiffness term. The presence of T a in this term introduce an undesired coupling and this term is kept on 265 the right-hand-side of Equation 1. It is noted that the vertical force F z,N contributes to the softening of the tower. The main softening effect attributed to the RNA mass is included in the stiffness matrix, as described in Branlard (2019).  Figure 4. The rotational speed is well captured, indicating that 275 the rotational inertia is properly set, but also indicating that the drive-train torsion does not have a strong impact for this simulation. The overall trend of the tower-top displacements is also well captured, though more differences are present due to missing contributions from additional blade and tower DOFs, missing non-linearities and quadratic velocity forces.
The method from subsection 2.5 is used to estimate the bending moments along the tower from the tower-top displacement.
The results shown on the right of Figure 4 indicate that the overall trends and load levels are well estimated, but some offsets  are observed, which are function of height. A contribution to the moment may be missing in the current model. This will be taken into consideration when analysing the results from the KF analysis.

Application to wind turbine tower loads estimation
Some of the individual models presented in section 2 were briefly validated in section 3. The augmented KF described in subsection 2.3 is now used, combining the different models together with the measurements. The state and output equations 285 given in Equation 2 and Equation 3 are implemented. The state equation is discretized according to Equation 12 and provided to the KF algorithm. Results from the KF simulation, which combines a set of measurements with a model, will be referred to as "KF estimation". The values used for the covariance matrices, P and Q, are discussed in section 5.

Ideal cases without noise
The same simulation as the one presented in subsection 3.1 is used, which extends from region 0 to region 3. The measurements 290 sampled at 20 Hz are here directly taken from the OpenFAST simulation and not from a field experiment. This is obviously an ideal situation since no noise is present in the measurements. Further, the OpenFAST and KF models are based on the same parameters such as the mass and stiffness distribution. States and tower loads estimated using the KF model are compared with the simulation results in Figure 5. The signals are seen to be well estimated by the KF model over the entire range of operating regions. The error observed for the tower bottom moment is in the range of errors observed for the isolated mechanical system A turbulent simulation is run, at an average wind speed of 14 m/s with turbulence intensity of 0.14, to illustrate the differences in the power spectral density of the signals. The results are given in Figure 6 and commented further.  are not in the mechanical system (e.g. the second FA mode and the DT torsion) are still "captured" by the estimator via the measurements. The rotational speed is directly observable by the KF, so the signal is obviously well estimated. The thrust, is 300 estimated based on the rotational speed and thus exhibits similar frequencies as the rotational speed, which is not the case for the reference thrust signal. The integration of the acceleration into the TT position (q t ) shows a higher frequency content than the reference signal. The second FA frequency has a strong energy content in the estimated q t signal. This frequency content 13 https://doi.org/10.5194/wes-2020-55 Preprint. Discussion started: 13 March 2020 c Author(s) 2020. CC BY 4.0 License. comes from the acceleration signal, but it is not sufficiently captured and damped by the model which does not represent the 2nd mode. A moving average filter of period 1 s was introduced to reduce the high-frequency content of the acceleration. The 305 results are labelled "Estimation w/filt." on the figure. The analysis of the moment spectrum given on the right of Figure 6 indicates that the frequencies are well captured but the overall content at frequencies beyond the 1st FA mode is too high. This is indicated by the values of the equivalent loads which are respectively 20 MNm and 30 MNm for the reference and estimated signal, using a Wöhler slope of m = 5. The low-pass filter on the acceleration signal greatly improves the spectrum of M y . The error in equivalent loads is further quantifies in the next paragraph.

Simulations with noise
The simulations presented in subsection 4.1 used as measurements the simulated values from OpenFAST. In this section, a Gaussian noise is added to each of the OpenFAST signals in order to account for measurement uncertainties. The noise level is taken a 10% of the standard deviation of the signal simulated by OpenFAST. A noise level of 20% will be referred to as "Large noise". Simulations were performed with OpenFAST for 10 wind speeds, with six different turbulent seeds for each wind speed.

315
A noise level was applied to these simulation results, prior to feeding them to the KF estimator. Cases with or without applying the low-pass filter on the (noisy) acceleration input were tried. Results for the error in equivalent load and standard deviation of the TB moment are shown in Figure 7. The equivalent loads are estimated using a Wöhler slope of m = 5. As expected,  the errors in standard deviation and equivalent loads follow similar trends. Errors without filtering are several fold larger than when the acceleration is filtered. Without noise, the equivalent loads are estimated with ±8% error. The error increases with the noise level and the equivalent loads appear to be mostly overestimated. Further tuning of the filter and of the covariance matrices involved in the KF may reduce the error. Further discussions are provided in section 5.

Computational time
The framework is written in the noncompiled Python language. The code was run on a single CPU. The average computational time for a 10 min period of measurements at 20 Hz was 37 s. Doubling the frequencies and the number of DOF would still 325 keep the computational time several fold smaller than real time. The expensive part of the algorithm is the non-linear solve needed to find the optimal wind speed (Equation 14).

Discussions
Measurements The results presented in the current study remained within the simulation realm. The accuracy of the method under uncertain conditions was partly quantified using various noise models. Yet, future work will evaluate the model using 330 field measurement data.
Model choices As mentioned in subsection 2.3, a certain level of choice is present as to whether the loads are placed as an input or within the state vector. A consequence is that different load models may also be implemented, for instance, models of higher order that the one used in Equation 8. In the current study, a "random walk" force model was used for the torque, and the thrust was set as a dependent variable of the torque. Yet, these loads are functions of the axial inductions, which typically 335 are assumed to follow a second order model referred to as dynamic wake. A linearization of this model could be applied to the aerodynamic thrust and torque and potentially improve the performance prediction of the Kalman filter.
Non-linearities and time-invariance This study assumed a linear form of the equation of motion and that the system matrices were time-invariant. Despite this crude assumption, reasonable results were obtained. Yet, further improvement are likely to be obtained if these assumptions are lifted. A simple approach would consist in updating the system matrices at some given interval 340 based on a slow moving average on the wind speed or the tower-top position. A more advanced method would use filtering methods that are adapted to non-linear systems, such as extended Kalman filters or particles filters. This approach would yet greatly increase the computational time. A shortcoming of the current approach is that the linear form of the equation was established "by hand". A more systematic approach will be considered in the future, using the linearized form of the state matrices returned by OpenFAST, which would include aerodynamic damping directly.

345
Degrees of freedom and offshore application The general formalism presented in section 2 can be applied to more degrees of freedom than the 2DOF model used: adding more shape function for the tower, and including side-side motion, yaw, tilt, shaft torsion, and blade motions. The results from the 2DOF model appeared encouraging enough to limit ourselves to this set, but future work will consider the inclusion of additional DOF. The extension of the method to offshore application could be done by adding extra degrees of freedom for the substructure, or, by using shape functions that represent the entire support structure. The generalized force due to the wave loading would need to be included. This force may be modelled based on the wind speed, or assumed of the model noise (see subsection 2.3).
Model tuning Apart from the choices of degrees of freedom and model formulation, there remains a part of model tuning, through: the choice of covariance matrices, and, the potential filtering done on the measurements. As shown in subsection 4.2, the filtering of the acceleration was seen to greatly improve the performance of the model. A time constant of 1s was chosen 355 empirically for the filter, but this value may need to be adapted for other applications. The choice of values used for the covariance matrices is usually the main source of criticism for KF based models. Indeed, these values have a strong influence on the results, and they are usually tuned empirically. For the current method to be successfully applied on various wind plants, an automatic tuning procedure is required. In the current study, the covariance matrices of the process were set automatically based on the value of the standard deviation of the simulated signal at rated conditions. For the measurements, these values 360 were divided by two. It was found that this procedure lead to satisfactory results. A sensitivity study should be considered in future work to give further insight on the procedure, in particular if more states and measurements are used.
Wind speed estimation and standstill/idling condition The wind speed estimation model presented in subsection 2.4 is limited to cases where the turbine is operating. Also, the accuracy of this model is crucial for the determination of the thrust, which in turns determine the tower-top position and the tower loads. The nacelle velocity was for instance omitted in the current 365 study and could be considered in future studies. Wind speed estimation is a field in which the industry has a great expertise.
Improvements on the algorithm would benefit the model presented in this paper.
Airfoil performance The performance of the airfoils is a large source of uncertainty which was not addressed. The thrust was determined using tabulated C T data, which may be significantly affected by the airfoil performance, which in turn are affected by blade erosion or other roughness sources, and, additional uncertainty on the aerodynamic modelling. Further improvement 370 of the model is thus required to provide an accurate determination of the thrust that would account for such unknowns. Air density should also be considered for a correct account of the loading if a tabulated approach is used.

Conclusions
The paper presented a general approach using Kalman filtering to estimate loads on a wind turbine, combining a mechanical model and a set of readily available measurements. An open source framework was established in hope to be further applied for 375 real-time fatigue estimation of wind turbine loads, providing inspiration for a digital twin concept. As an example, the equations for a 2DOF system of a wind turbine were presented, and this system was used throughout the article. The study focused on the estimation of tower bending moment and in particular the associated damage equivalent load. Based on simulation results, the estimator was seen to be able to capture the damage equivalent loads with an accuracy in the order of 10%. Future work will address the following points: use of field measurements, offshore application of the method, increased number of 380 DOF, automatic covariance tuning, improved wind speed estimation in standstill, improved thrust determination in off-design conditions, and use of a linearized model obtained from an aero-servo-elastic tool.