Preventing wind turbine tower natural frequency excitation with a quasiLPV model predictive control scheme^{†}^{†}thanks: Short title: A qLPVMPC framework for preventing tower natural frequency excitation
Abstract
[Abstract]With the ever increasing size and costs of wind turbines, structural load minimization is crucial to extend the turbine lifetime, and to minimize operation and maintenance expenses. In practice, the turbine rotor has a mass imbalance, which induces a periodic and rotorspeed dependent tower sideside excitation during belowrated operation. Operating at this frequency generally degrades the expected structural life span, as the first tower natural frequency for larger turbines coincides with a belowrated operational rotor speed, and the sideward direction has negligible aerodynamic damping. To reduce this effect, earlier work has shown the effectiveness of active tower damping control strategies using collective pitch control. A more passive approach is frequency skipping by inclusion of speed exclusion zones, which avoids prolonged operation near the critical frequency. However, neither of the methods incorporate a convenient way of performing a tradeoff between energy maximization and fatigue load minimization. Therefore, this paper introduces a quasilinear parameter varying model predictive control (qLPVMPC) scheme, exploiting the beneficial (convex) properties of a qLPV system description. The qLPV model is obtained by a modulation transformation and augmented with a simple wind turbine model. Results show the effectiveness of the algorithm in synthetic and realistic simulations using the NREL 5MW reference wind turbine in highfidelity simulation code. Prolonged rotor speed operation at the tower natural frequency is prevented, whereas when the tradeoff is in favor of energy production, the algorithm decides to rapidly pass over the tower natural frequency to attain higher rotor speeds and power production rates.
Research Article
1]Sebastiaan Paul Mulders 2]Tobias Gybel Hovgaard 2]Jacob Deleuran Grunnet 1]JanWillem van Wingerden
S.P. MULDERS et al
Sebastiaan Paul Mulders, Delft Center for Systems and Control, Faculty of Mechanical Engineering, Mekelweg 2, 2628 CD Delft, The Netherlands.
1 Introduction
In practical scenarios, the center of mass of the wind turbine rotor assembly does not coincide with the actual rotor center as a result of, e.g., manufacturing imperfections, wear and tear, fouling and icing ^{1}. Moreover, vibrations are also induced by rotor aerodynamic imbalances caused by pitch errors and damage to the blade surface ^{2}. Consequently, during variablespeed belowrated operation, the bladepassing frequency may excite the structural sideside natural frequency. Because the turbine rotor provides negligible aerodynamic sideside damping, at an order of magnitude smaller than the foreaft damping ratio, small perturbations can lead to load fluctuations comparable to foreaft stresses ^{3}. As a result, excitation of the sideside mode possibly results in accelerated and accumulative fatigue damage.
Straightforward control implementations are available for reducing and mitigating the excitation of the tower foreaft/sideside modes. An active method for reducing tower motion is the use of an integrated nacelle acceleration signal in a proportional feedback structure. Depending on the measured acceleration direction, the resulting signals form an addition to the collective pitch^{4, 5} or generator torque^{6} control signal, for respective damping of foreaft and sideside vibrations. Another, more passive method, entails prevention of structural mode excitation by increasing the generator torque when the rotor speed approaches the excitation frequency^{7}. The method is often referred to as frequency skipping by inclusion of speed exclusion zones.
All of the abovedescribed active and passive methods complicate the controller design, by requiring extra proportionalintegralderivative (PID) feedback control loops with additional logic and speed setpoints. Also, the methods do not incorporate convenient and inherent tuning capabilities for a tradeoff between produced energy and fatigue loading. Therefore, more advanced control algorithms might form a solution by providing a more integrated way of controller synthesis, incorporating both power and load objectives. While an abundance of publications on advanced wind turbine control algorithms outline the possible benefits^{8}, to the authors’ knowledge, more sophisticated control methods do not yet see a widespread adoption in industrialgrade wind turbine control systems; PID control structures^{9} provide ease of implementation while resulting in a sufficient performance level.
An advanced control method, that has seen a substantial gain of interest from industry in the past decades, is model predictive control (MPC)^{10, 11}. The most evident benefits of MPC over PID control^{12} are (1) the ability of including constraints, (2) coping with the complexity of nonminimum phase systems, (3) robustness against deviations of the control model to the actual process, and (4) the convenient application to multivariable control problems. MPC has been considered for wind turbine load mitigations in the literature. A nonlinear MPC (NMPC) method is applied by assuming future wind speed knowledge using a lightdetection and ranging (LIDAR) system^{13}. Simulation results shows promising load reductions without affecting the energy production. Furthermore, a robust MPC (RMPC) implementation is compared to a nominal MPC control structure for the purpose of active tower foreaft damping^{14}. It is shown in numerical simulations that the former outperforms the latter mentioned, as in particular around rated operating conditions, physical actuations constraints form a limiting factor. In the same operating region, the benefits of NMPC using a future wind speed prediction^{15} are once again emphasized.
All of the above described MPC implementations focus on the active mitigation of structural loads. A more passive MPC implementation, providing frequency skipping capabilities, and thereby making an optimal tradeoff between turbine loads and energy production over the prediction horizon, does not seem to have been backed up by literature in the past. The complexity lies in the fact that fatigue loads are minimized by preventing operation at the tower natural frequency, while it is essential to cross the same frequency for attaining higher rotational speeds and power production rates. The conflicting objectives form a burden for describing the objective as a convex optimization problem. Moreover, NMPC for solving nonconvex problems is – because of its computational complexity – often considered unsuitable for realworld applications.
Imposing spectral constraints might form a possible solution path^{16}, by employing the shorttime Fourier transform (STFT) on the system output signal in a nonlinear MPC setting. A similar methodology^{17} uses the Selective Discrete Fourier Transform (SDFT) in an MPC approach to dampen oscillation modes in power system stabilizers (PSS). However, from an implementation and tuning perspective, a frequency domain approach seems to be unintuitive and nontrivial. Therefore, in this paper, another approach is considered, which consists of a model modulation transformation described for application and control in the field of Tapping ModeAtomic Force Microscope (TMAFM)^{18}. The model transformation is applied to the turbine tower model, and transfers frequencydependent magnitude and phase content to a quasi steadystate contribution. This is accomplished by converting a lineartime invariant (LTI) system into a linearparameter varying (LPV) model, scheduled on the excitation frequency. The technique shows similarities with the MultiBlade Coordinate (MBC) transformation, often used in individual pitch control (IPC) implementations for blade fatigue load reductions^{19, 20}.
An LPV system representation is frequently used for capturing nonlinear dynamics into a system description with a linear inputoutput mapping^{21}. An external scheduling variable varies the dynamics of the linear model. Now, consider the combination of an LPV model with MPC. The modelbased control method uses a mathematical system description to compute an optimal control signal over the prediction horizon. Unfortunately, for LPV systems, the considered model is subject to changes over time, described by the yet unknown scheduling trajectory. However, when the system is scheduled on state variables and/or input signals, the model is referred to as a quasiLPV (qLPV) system. Recently, an efficient MPC scheme for such qLPV systems is proposed^{22} by solving subsequent quadratic programs (QPs).
This paper subjects a tower model to the modulation transformation and augments it with a simplified wind turbine model, such that a qLPV model is obtained. The result is combined with an iterative MPC method, exploiting the beneficial properties of qLPV systems^{22}. Using the qLPVMPC framework, a methodology for performing an optimal tradeoff between produced energy and tower loads is presented, and thereby provides the following contributions:

•
Providing the derivation results of a model modulation transformation, for moving the magnitude and phase frequency content at the excitation frequency to a quasi steadystate contribution

•
Applying the transformation to a secondorder tower model and showcasing its working principles by an illustrative example;

•
Combining the transformed tower model with a simplified wind turbine model and linearizing at belowrated operating points, for obtaining a qLPV statespace system description

•
Discretizing and converting the qLPV model to its affine form

•
Formally deriving the efficient MPC approach for affine qLPV model structures

•
Showcasing the proposed approach in closedloop highfidelity simulations with different wind profiles, to clearly show its effectiveness and practical applicability.
The paper is organized as follows. Section 2 describes a methodology for transforming a nominal tower model into a modulated LPV system description. In Section 3, the obtained system is combined with a simplified wind turbine model, resulting in a qLPV system description after linearization. Next, in Section 4, the efficient MPC scheme is combined with the qLPV model to make an optimal and userdefined tradeoff between tower loads and generated energy for the prevailing environmental conditions. In Section 5, the qLPVMPC framework is evaluated with highfidelity simulations using the NREL 5MW reference wind turbine, subject to synthetic and realistic wind profiles. Finally, conclusions are drawn in Section 6.
2 Problem formalization and tower model modulation transformation
For performing a produced energy versus tower fatigue load tradeoff, a wind turbine model needs to be combined with a structural tower model. Section 2.1, describes the tower sideside dynamics by a secondorder massdamperspring system. Section 2.2 formalizes the problem statement, and shows that the combined system description results in nonconvexity. Therefore, in Section 2.3, this nominal turbine tower model is subject to a modulation transformation to facilitate its convexification. An analysis on the effects and implications of the transformation are clarified by an illustrative example in Section 2.4.
2.1 Modeling the tower dynamics as a secondorder system
In practical scenarios, the center of mass of a wind turbine rotor is likely to not coincide with the rotor center. In effect, as largescale stateoftheart wind turbines are operated with a variablespeed control strategy for belowrated conditions, the support structure is excited by a periodic and frequencyvarying centripetal force, as illustrated in Figure 1. The tower dynamics, excited by a rotorspeed dependent onceperrevolution (1P) periodic force, are modeled by a secondorder massdamperspring system
\displaystyle m\ddot{x}(t)+\zeta\dot{x}(t)+kx(t)  \displaystyle=a_{\mathrm{u}}\cos{(\psi(t))},  (1) 
in which {\left\{m,\,\zeta,\,k\right\}\in\mathbb{R}^{+}} are respectively the constant first mode modal mass, modal damping and modal stiffness, {\psi(t)\in[0,\,2\pi)} is the rotor azimuth angle, {a_{\mathrm{u}}\in\mathbb{R}^{+}} quantifies the periodic force amplitude, and {\left\{x,\,\dot{x},\,\ddot{x}\right\}\in\mathbb{R}} respectively represent the sideside displacement, velocity and acceleration in the hub coordinate system, illustrated in Figure 1. A secondorder system is taken, as this represents the first mode of the tower using the well known modaldecomposition model reduction technique^{23, 24}, and allows for a convenient assessment and derivation of the modulation transformation in the next section: the application to higherorder models is also possible and would result in a similar analysis. Furthermore, the force amplitude a_{\mathrm{u}} is assumed to be constant for all rotational speeds, however, as will be shown later, this assumption can be relaxed for mildly varying amplitude changes.
The system in Equation (1) is split in a set of firstorder differential equations by defining x_{1}=\dot{x}(t), x_{2}=x(t), such that it is rewritten in the standard statespace \dot{\boldsymbol{x}}=\mathbf{A_{\mathrm{g}}}\boldsymbol{x}+\mathbf{B_{\mathrm{% g}}}\boldsymbol{u} representation
\displaystyle\begin{bmatrix}\dot{x}_{1}\\ \dot{x}_{2}\end{bmatrix}=\begin{bmatrix}\mathrm{\zeta/m}&\mathrm{\omega}_{% \mathrm{n}}^{2}\\ 1&0\end{bmatrix}\begin{bmatrix}x_{1}\\ x_{2}\end{bmatrix}+\begin{bmatrix}\mathrm{a_{\mathrm{u}}}\\ 0\end{bmatrix}\cos{(\psi(t))}\qquad\text{and}\qquad G(s)\lx@stackrel{{% \scriptstyle s}}{{=}}\left[\begin{array}[]{cc}\mathbf{A_{\mathrm{g}}}&\mathbf% {B_{\mathrm{g}}}\\ \hline\mathbf{C_{\mathrm{g}}}&\mathbf{0}\end{array}\right],  (2) 
in which \mathbf{A_{\mathrm{g}}}\in\mathbb{R}^{n_{\mathrm{g}}\times n_{\mathrm{g}}}, \mathbf{B_{\mathrm{g}}}\in\mathbb{R}^{n_{\mathrm{g}}}, and \omega_{\mathrm{n}}=\sqrt{k/m} is the structural natural frequency. All states are assumed to be measured thus \mathbf{C_{\mathrm{g}}}=I_{\mathrm{n_{g}}}. The system G(s) has a statespace realization given by the quadruple \left(\mathbf{A_{\mathrm{g}}},\,\mathbf{B_{\mathrm{g}}},\,\mathbf{C_{\mathrm{g% }}},\,\mathbf{0}\right)^{25}, and s is the Laplace variable.
2.2 Problem formalization
This section formalizes the problem considered in this paper. The aim is to provide a tradeoff between energy generation efficiency and tower fatigue load reductions, by preventing rotor speed operation near the tower natural frequency. The considered nominal framework is graphically presented in Figure 2. The wind turbine model has a wind disturbance and a generator torque control input, the latter of which is subject to optimization. A cosine function acts on the azimuth position output from the wind turbine model, which results in a periodic input to the tower model. The load and energy outputs of the respective tower and wind turbine models, together with the torque input signal are included in the following cost function to optimize the energyload tradeoff
\displaystyle\underset{Torque}{\mathrm{arg\,min}}  \displaystyle\lambda_{\mathrm{1}}(Energy)+\lambda_{\mathrm{2}}(Loads)+\lambda% _{\mathrm{3}}(Torque),  (3) 
where \lambda_{i}~{}\text{and}~{}i=\{1,2,3\} are positive weighting constants determining the objective tradeoffs. The above given relation presents the optimization objective in an informal fashion for illustration purposes; later sections define the problem in a mathematical correct way. The load signal is a periodic and rotorspeed dependent measure for tower fatigue loading, caused by the presence of the trigonometric function. This forms a burden for describing the objective as a convex optimization problem.
The solution path employed in this paper is the aggregation of the nonlinear trigonometric function and the LTI tower model by a model modulation transformation. The transformed tower model results in an LPV system description. The subsequent combination with a wind turbine model, providing the rotor speed scheduling variable as an internal system state, results in a quasiLPV model. The next section provides theory for derivation of the model modulation transformation, whereas later sections elaborate on an efficient MPC method exploiting the beneficial properties of qLPV model structures
2.3 Theory on the tower model modulation transformation with periodic excitation towards an LPV representation
In signal analysis, modulation is the process of imposing a signal onto a carrier wave^{26}. However, in this section, modulation is not applied to a signal, but is performed on the tower model introduced in the previous section. The aim of the modulation transformation is to obtain a linear (but parameter varying) system description, which provides the frequency dependent dynamical behavior as a steadystate signal. The modulated signal is in a later section used to form a convex quadratic optimization problem in an MPC setting for computing the optimal control signal over the prediction horizon.
The model modulation transformation is inspired by the work of^{18} and only the main results are given in this section. The derivation is performed and validated with symbolic manipulation software for algebraic expressions^{27}, of which the code is made publicly available^{28}. The transformation relies on the assumption that the change in system response amplitude a_{\mathrm{y}}(\tau) and phase change \phi(\tau) is much slower than that of the driving excitation frequency \omega_{\mathrm{r}}. For this reason, the slower time scale is indicated by \tau as a substitute for the normal time scale t. Variables that are a function of the slowvarying time scale are assumed to be constant over a single period T_{\mathrm{r}}=2\pi/\omega_{\mathrm{r}} of the excitation, such that
\displaystyle\int^{T_{\mathrm{r}}}_{0}f(\tau)g(t)dt  \displaystyle=f(\tau)\int^{T_{\mathrm{r}}}_{0}g(t)dt.  (4) 
With respect to a driving periodic input, the steadystate response of the tower displacement is a magnitude and phase change at an equal frequency
\displaystyle x_{i}(t)  \displaystyle=a_{i}(\tau)\cos{(\omega_{\mathrm{r}}t+\phi(\tau))},  (5) 
where a_{i}\in\mathbb{R}^{+} is the amplitude and the subscript i\in\mathbb{Z}^{+} is a counter variable. By taking into account the introduced timescales, and using Euler’s formula e^{\mathrm{j}\phi}=\cos(\phi)+\mathrm{j}\sin{(\phi)}, the state variables are rewritten as
\displaystyle x_{i}(t)  \displaystyle=\Re{\left\{a_{i}(\tau)e^{\mathrm{j}(\omega_{\mathrm{r}}t+\phi(% \tau))}\right\}},  (6) 
with i=\{1,2\}, and \mathrm{j}=\sqrt{1} is the imaginary unit. The symbol \Re\{\cdot\} indicates the real part of a given expression, whereas \Im{\{\cdot\}} is used to represent the imaginary part. The slow varying term {X_{i}\in\mathbb{C}} is now written as a product with the fast harmonic function, with a fixed phase and amplitude, such that
\displaystyle x_{i}(t)  \displaystyle=\Re{\left\{a_{i}(\tau)e^{\mathrm{j}\phi(\tau)}e^{\mathrm{j}% \omega_{\mathrm{r}}t}\right\}}=\Re{\left\{X_{i}(\tau)e^{\mathrm{j}\omega_{% \mathrm{r}}t}\right\}}.  (7) 
and taking the first time derivative gives
\displaystyle\dot{x}_{i}(t)=\Re{\left\{\left(\dot{X}_{i}(\tau)+\mathrm{j}% \omega_{\mathrm{r}}X_{i}(\tau)\right)e^{\mathrm{j}\omega_{\mathrm{r}}t}\right% \}}.  (8) 
By substitution of Equations (7) and (8) in the nominal state space representation of Equation (2), the following expressions are obtained
\displaystyle\Re\left\{\left(\dot{X}_{1}(\tau)+j\omega_{\mathrm{r}}X_{1}(\tau)% +(\zeta/m)X_{1}(\tau)+\omega_{\mathrm{n}}^{2}X_{2}(\tau)\right)e^{j\omega_{% \mathrm{r}}t}a_{\mathrm{u}}e^{j\omega_{\mathrm{r}}t}\right\}=0,  (9)  
\displaystyle\Re{\left\{\left(\dot{X}_{2}(\tau)+j\omega_{\mathrm{r}}X_{2}(\tau% )X_{1}(\tau)\right)e^{j\omega_{\mathrm{r}}t}\right\}}=0.  (10) 
Furthermore, the following property of orthogonality is used, where
\displaystyle\int_{0}^{T_{\mathrm{r}}}\Re{\left\{Ce^{j\theta}\right\}}e^{j% \theta}d\theta=0,  (11) 
if and only if \{C\in\mathbb{C}\}=0. Thus, Equations (9) and (10) are multiplied with e^{\mathrm{j}\omega_{\mathrm{r}}t}, such that
\displaystyle\int_{0}^{T_{\mathrm{r}}}\Re\left\{\left(\dot{X}_{1}(\tau)+j% \omega_{\mathrm{r}}X_{1}(\tau)+(\zeta/m)X_{1}(\tau)+\omega_{\mathrm{n}}^{2}X_{% 2}(\tau)\right)e^{j\omega_{\mathrm{r}}t}a_{\mathrm{u}}e^{j\omega_{\mathrm{r}}% t}\right\}e^{\mathrm{j}\omega_{\mathrm{r}}t}dt=0,  (12)  
\displaystyle\int_{0}^{T_{\mathrm{r}}}\Re{\left\{\left(\dot{X}_{2}(\tau)+j% \omega_{\mathrm{r}}X_{2}(\tau)X_{1}(\tau)\right)e^{j\omega_{\mathrm{r}}t}% \right\}}e^{\mathrm{j}\omega_{\mathrm{r}}t}dt=0.  (13) 
Termbyterm integration of the integrals in Equation (13) using the mathematical property in Equation (4), gives the following result
\displaystyle\begin{bmatrix}\dot{X}_{1}\\ \dot{X}_{2}\end{bmatrix}  \displaystyle=\begin{bmatrix}\mathrm{j}\omega_{\mathrm{r}}X_{1}(\mathrm{\zeta% }/\mathrm{m})X_{1}\omega_{\mathrm{n}}^{2}X_{2}+a_{\mathrm{u}}\\ \mathrm{j}\omega_{\mathrm{r}}X_{2}+X_{1}\end{bmatrix}.  (14) 
Now, by defining a new state sequence \boldsymbol{q}=[q_{1},\,q_{2},\,q_{3},\,q_{4}]^{\mathrm{T}}=[\Re\{X_{1}\},\,% \Im\{X_{1}\},\,\Re\{X_{2}\},\,\Im\{X_{2}\}]^{\mathrm{T}}, the system is rewritten as {\boldsymbol{\dot{q}}(\omega_{\mathrm{r}})=\mathbf{A}_{\mathrm{h}}(\omega_{% \mathrm{r}})\boldsymbol{q}+\mathbf{B}_{\mathrm{h}}}, such that
\displaystyle\begin{bmatrix}\dot{q}_{1}\\ \dot{q}_{2}\\ \dot{q}_{3}\\ \dot{q}_{4}\end{bmatrix}=\begin{bmatrix}\mathrm{\zeta}/\mathrm{m}&\omega_{% \mathrm{r}}&\omega_{\mathrm{n}}^{2}&0\\ \omega_{\mathrm{r}}&\mathrm{\zeta}/\mathrm{m}&0&\omega_{\mathrm{n}}^{2}\\ 1&0&0&\omega_{\mathrm{r}}\\ 0&1&\omega_{\mathrm{r}}&0\end{bmatrix}\begin{bmatrix}{q}_{1}\\ {q}_{2}\\ {q}_{3}\\ {q}_{4}\end{bmatrix}+\begin{bmatrix}a_{\mathrm{u}}\\ 0\\ 0\\ 0\end{bmatrix}\qquad\text{and}\qquad H(s,\omega_{\mathrm{r}})\lx@stackrel{{% \scriptstyle s}}{{=}}\left[\begin{array}[]{cc}\mathbf{A_{\mathrm{h}}}&\mathbf% {B_{\mathrm{h}}}\\ \hline\mathbf{C_{\mathrm{h}}}&\mathbf{0}\end{array}\right],  (15) 
in which \mathbf{A_{\mathrm{h}}}\in\mathbb{R}^{n_{\mathrm{h}}\times n_{\mathrm{h}}}, \mathbf{B_{\mathrm{h}}}\in\mathbb{R}^{n_{\mathrm{h}}} and n_{\mathrm{h}}=2n_{\mathrm{g}}. Again, all states are measured thus \mathbf{C_{\mathrm{h}}}=I_{\mathrm{n_{h}}}, and the system H(s,\omega_{\mathrm{r}}) has a statespace realization given by the quadruple \left(\mathbf{A_{\mathrm{h}}},\,\mathbf{B_{\mathrm{h}}},\,\mathbf{C_{\mathrm{h% }}},\,\mathbf{0}\right)^{25}. The instantaneous amplitude and phase of the dynamic system response at frequency \omega_{\mathrm{r}} are given by
\displaystyle a_{\mathrm{y}}(\tau)  \displaystyle=\sqrt{q_{3}^{2}+q_{4}^{2}},  (16)  
\displaystyle\phi(\tau)  \displaystyle=\tan^{1}{(q_{4}/q_{3})}.  (17) 
It is also possible to write the result of the derivation using a summation of Kronecker products
\displaystyle H(s,\omega_{\mathrm{r}})  \displaystyle\equiv\boldsymbol{\dot{q}}=\left(\mathbf{A}_{\mathrm{g}}\otimes I% _{\mathrm{n_{g}}}+I_{\mathrm{n_{g}}}\otimes\begin{bmatrix}0&\omega_{\mathrm{r}% }\\ \omega_{\mathrm{r}}&0\end{bmatrix}\right)\boldsymbol{q}+\left(\mathbf{B_{% \mathrm{g}}}\otimes\begin{bmatrix}1\\ 0\end{bmatrix}\right).  (18) 
The nominal and transformed model representations G(s) and H(s,\omega_{\mathrm{r}}) are interchangeable: Figure 3 graphically summarizes the transformation of the nominal periodically excited secondorder tower model (Equation (2)) into an LPV model structure (Equation (15)). The amplitude a_{\mathrm{u}} of the periodic input is in the modulated model a direct input to the system. The outputs are the amplitude a_{\mathrm{y}} and phase shift \phi with respect to the input frequency. Note that the frequency \omega_{\mathrm{r}} is in the transformed case a scheduling variable to the LPV system, changing the system dynamics. The following section demonstrates and further explains the effects of the presented transformation by an illustrative example.
2.4 Illustrating the effects of the transformation
This section further analyses the properties of the model modulation transformation, and showcases its characteristics by an illustrative example. For setting up a simulation case, tower model parameters are taken sort of arbitrarily: the tower mass is m=1000 kg, the damping coefficient \zeta=100 kg s^{1} and the spring constant k=500 kg s^{2}.
Figure 4 shows Bode magnitude plots of the nominal plant and its modulated counterpart. As a result of the selected tower characteristics, the natural frequency is located at \omega_{\mathrm{n}}\approx 0.7 rad s^{1}, with a clearly present resonance peak at the same frequency. To obtain the amplitude output (Equation (16)) of the modulated model over {\omega\in\left\{\boldsymbol{\Omega}=[10^{2}\,\dots\,10^{1}]\right\}} rad s^{1}, the Euclidean norm of the frequency responses of q_{3} and q_{4} at each frequency point is taken.
For illustration purposes, a set of 4 evaluation frequencies is chosen as \left\{\omega_{\mathrm{r},1}\,\dots\,\omega_{\mathrm{r},4}\right\}=\left\{0,\,% 0.5,\,0.7,\,2.0\right\} rad s^{1} to showcase the effects of the transformation by indicative pointers in Figure 4. For \omega_{\mathrm{r},1}=0 rad s^{1} the transformed model reduces to the nominal case; for \left\{\omega_{\mathrm{r},i}\,\forall\,i>1\right\} the magnitude content is transferred to a DC contribution for each evaluation of H(\mathrm{j}\omega,\omega_{\mathrm{r},i}). Moreover, the nominal resonance peak at \omega_{\mathrm{n}} is for each frequency response split into two peaks with a 3 dB magnitude reduction. Thus, when the input amplitude a_{\mathrm{u}} of the transformed model is varied slowly, the magnitude from specific frequency points of the nominal model is mapped to a DC contribution in the transformed case; rapid variations will result in contributions from the resonances at higher frequencies. However, as in this paper a_{\mathrm{u}} is assumed constant, additional measures to reduce these effects, such as lowpass or notch filters, are dispensable.
The effect of the transformation in the timedomain is evaluated in Figure 5. Both the nominal and transformed model are fed with a linear increasing frequency, at a constant frequency acceleration rate of \dot{\omega}_{\mathrm{r}}=0.01 rad s^{2} for a period of 1200 s, starting from \omega_{\mathrm{r}}=0 to 1.2 rad s^{1}. This frequency range is chosen as modern largescale variablespeed wind turbines are controlled in this operating region. The transformed model shows a very close amplitude tracking of the nominal model dynamics. The earlier imposed restriction and assumption on the change in amplitude and phase by a slow time scale \tau, does not seem to limit the proposed method for applicability to the considered wind turbine control objective.
3 Wind turbine model augmentation and linearization
This section considers the derivation of a simple NREL 5MW (linear) model, for augmentation to the modulated tower model such that a quasiLPV model is obtained. Section 3.1 provides the simple firstorder wind turbine model. Next, in Section 3.2, the model is symbolically linearized and augmented with the transformed tower model in a qLPV representation. Section 3.3 provides linearization parameters over the complete below–rated operating region based on the properties of the NREL 5MW reference wind turbine^{29}. Finally, Section 3.4 validates the firstorder and affine linear models to simulation results of their nonlinear equivalent.
3.1 Simplified wind turbine system description
Because the dynamics of the transformed tower model H(s,\omega_{\mathrm{r}}) are scheduled by the input excitation frequency, which is in this case the (1P) rotor speed, it is a logical step to augment a wind turbine model adding this state to the overall system description. Aa system of which the scheduling variable is part of the state vector is known as a qLPV system description. The considered firstorder wind turbine model is
\displaystyle J_{\mathrm{r}}\dot{\omega}_{\mathrm{r}}  \displaystyle=\tau_{\mathrm{a}}\underbrace{N(\tau_{\mathrm{g}}+\Delta\tau_{% \mathrm{g}})}_{\tau_{\mathrm{s}}},  (19) 
in which J_{\mathrm{r}}\in\mathbb{R}^{+} is the total rotor inertia consisting out of the hub and 3 times the blade inertia, \left\{N\geq 1\right\}\subset\mathbb{R}^{+} is the gearbox ratio, and \tau_{\mathrm{a}} is the aerodynamic rotor torque defined as
\displaystyle\tau_{\mathrm{a}}  \displaystyle=\frac{1}{2}\rho_{\mathrm{a}}\pi R^{3}U^{2}C_{\mathrm{\tau}}(% \lambda,\beta),  (20) 
where \rho_{\mathrm{a}}\in\mathbb{R}^{+} is the air density, R\in\mathbb{R}^{+} the rotor radius, U\in\mathbb{R}^{+} the rotor effective wind speed and C_{\mathrm{\tau}}\in\mathbb{R} the torque coefficient as a function of the blade pitch angle \beta and the dimensionless tipspeed ratio \lambda=\omega_{\mathrm{r}}R/U. The system torque \tau_{\mathrm{s}}\in\mathbb{R}^{+} is a summation of the generator torque \tau_{\mathrm{g}}\in\mathbb{R}^{+} resulting from a standard Komegasquared torque control law, and \Delta\tau_{\mathrm{g}}\in\mathbb{R} is an additional torque contribution resulting from the MPC framework described later in this paper. The torque control law is taken as an integral part of the model, and is defined as
\displaystyle\tau_{\mathrm{g}}=K\omega_{\mathrm{r}}^{2}/N,  (21) 
in which K\in\mathbb{R}^{+} is the optimal mode gain
\displaystyle K=\frac{\pi\rho_{\mathrm{a}}R^{5}C_{\mathrm{p}}(\lambda,\beta)}{% 2\lambda^{3}},  (22) 
calculated for the lowspeed shaft (LSS) side.
3.2 Linearizing the augmented turbine and tower model
This section augments the wind turbine model from Section 3.1 to the modulated tower model H(s,\omega_{\mathrm{r}}), such that the following system is obtained
\displaystyle\begin{bmatrix}\dot{q}_{1}\\ \dot{q}_{2}\\ \dot{q}_{3}\\ \dot{q}_{4}\\ \dot{\omega}_{\mathrm{r}}\end{bmatrix}  \displaystyle=\begin{bmatrix}\mathrm{\zeta}/\mathrm{m}&\omega_{\mathrm{r}}&% \omega_{\mathrm{n}}^{2}&0&0\\ \omega_{\mathrm{r}}&\mathrm{\zeta}/\mathrm{m}&0&\omega_{\mathrm{n}}^{2}&0\\ 1&0&0&\omega_{\mathrm{r}}&0\\ 0&1&\omega_{\mathrm{r}}&0&0\\ 0&0&0&0&0\end{bmatrix}\begin{bmatrix}{q}_{1}\\ {q}_{2}\\ {q}_{3}\\ {q}_{4}\\ \omega_{\mathrm{r}}\end{bmatrix}+\begin{bmatrix}a_{\mathrm{u}}\\ 0\\ 0\\ 0\\ (\tau_{\mathrm{a}}N\left(\tau_{\mathrm{g}}+\Delta\tau_{\mathrm{g}}\right))/J_% {\mathrm{r}}\end{bmatrix},  (23)  
\displaystyle a_{\mathrm{y}}  \displaystyle=\sqrt{(q_{3}^{2}+q_{4}^{2})}.  (24) 
The abovegiven system description contains the nonlinear aerodynamic and generator torque input defined previously by Equations (20) and (21). Furthermore, the output a_{\mathrm{y}} is a nonlinear combination of state vector elements. The system is subject to linearization, where the desired linear state, input and outputs vectors are defined as
\displaystyle\boldsymbol{\hat{q}}(t)  \displaystyle=\left[\hat{q}_{1},\,\hat{q}_{2},\,\hat{q}_{3},\,\hat{q}_{4},\,% \hat{\omega}_{\mathrm{r}}\right]^{\rm T},  (25)  
\displaystyle\boldsymbol{\hat{u}}(t)  \displaystyle=\left[\hat{U},\,\Delta\hat{\tau}_{\mathrm{g}}\right]^{\rm T},  (26)  
\displaystyle\boldsymbol{\hat{y}}(t)  \displaystyle=\hat{A}_{\mathrm{y}},  (27) 
and the \hat{(\cdot)}notation indicates the deviation with respect to the considered linearization point. Now, the system is linearized by taking the partial derivatives of Equations (23) and (24) with respect to the state and inputs vectors, such that a linear statespace system is obtained
\displaystyle\boldsymbol{\dot{\hat{q}}}(t)  \displaystyle=\mathbf{A}(\boldsymbol{p})\boldsymbol{\hat{q}}(t)+\mathbf{B}(% \boldsymbol{p})\boldsymbol{\hat{u}}(t)  (28)  
\displaystyle\boldsymbol{\hat{y}}(t)  \displaystyle=\mathbf{C}(\boldsymbol{p})\boldsymbol{\hat{q}}(t), 
where the state, input and output matrices are defined as
\displaystyle\mathbf{{A}}(\boldsymbol{p})  \displaystyle=\begin{bmatrix}\mathrm{\zeta}/\mathrm{m}&\bar{\omega}_{\mathrm{% r}}&\omega_{\mathrm{n}}^{2}&0&\bar{q}_{2}\\ \bar{\omega}_{\mathrm{r}}&\mathrm{\zeta}/\mathrm{m}&0&\omega_{\mathrm{n}}^{% 2}&\bar{q}_{1}\\ 1&0&0&\bar{\omega}_{\mathrm{r}}&\bar{q}_{4}\\ 0&1&\bar{\omega}_{\mathrm{r}}&0&\bar{q}_{3}\\ 0&0&0&0&(\bar{k}_{\mathrm{\omega_{\mathrm{r}}}}N\bar{k}_{\tau_{\mathrm{g}}})/% J_{\mathrm{r}}\end{bmatrix},\quad\mathbf{B}(\boldsymbol{p})=\begin{bmatrix}0&0% \\ 0&0\\ 0&0\\ 0&0\\ \bar{k}_{\mathrm{U}}/J&N/J\end{bmatrix},\quad\mathbf{C}(\boldsymbol{p})=\frac% {1}{2}\begin{bmatrix}0\\ 0\\ \bar{q}_{3}(\bar{q}_{3}^{2}+\bar{q}_{4}^{2})^{1/2}\\ \bar{q}_{4}(\bar{q}_{3}^{2}+\bar{q}_{4}^{2})^{1/2}\\ 0\end{bmatrix}^{\mathrm{T}}.  (29) 
The aerodynamic rotor torque is linearized with respect to the rotor speed and wind speed
\displaystyle\hat{\tau}_{\mathrm{a}}  \displaystyle=\frac{\partial\tau_{\mathrm{a}}}{\partial\omega_{\mathrm{r}}}% \hat{\omega}_{\mathrm{r}}+\frac{\partial\tau_{\mathrm{a}}}{\partial U}\hat{U}=% \bar{k}_{\mathrm{\omega_{\mathrm{r}}}}({\omega_{\mathrm{r}}},~{}{\beta},~{}{U}% )\hat{\omega}_{\mathrm{r}}+\bar{k}_{\mathrm{U}}(\omega_{\mathrm{r}},\beta,U)% \hat{U},  (30) 
with
\displaystyle\bar{k}_{\mathrm{\omega_{\mathrm{r}}}}(\omega_{\mathrm{r}},\beta,U)  \displaystyle=\left.c_{\mathrm{r}}RU\frac{\partial C_{\mathrm{\tau}}(\lambda,% \beta)}{\partial\lambda}\right_{\omega_{\mathrm{r}}=\bar{\omega_{\mathrm{r}}}% ,\,\beta=\bar{\beta},\,U=\bar{U}},  (31)  
\displaystyle\bar{k}_{\mathrm{U}}(\omega_{\mathrm{r}},\beta,U)  \displaystyle=\left.2c_{\mathrm{r}}C_{\mathrm{\tau}}(\lambda,\beta)Uc_{% \mathrm{r}}\omega_{\mathrm{r}}R\frac{\partial C_{\mathrm{\tau}}(\lambda,\beta)% }{\partial\lambda}\right_{\omega_{\mathrm{r}}=\bar{\omega}_{\mathrm{r}},\,% \beta=\bar{\beta},\,U=\bar{U}},  (32) 
and c_{\mathrm{r}}=0.5\rho\pi R^{3} is a constant factor. Finally, the Komegasquared torque controller is linearized as
\displaystyle\hat{\tau}_{\mathrm{g}}(\omega_{\mathrm{r}})  \displaystyle=\left.\frac{\partial\tau_{\mathrm{g}}}{\partial\omega_{\mathrm{r% }}}\hat{\omega}_{\mathrm{r}}=\bar{k}_{\tau_{\mathrm{g}}}(\omega_{\mathrm{r}})% \hat{\omega}_{\mathrm{r}}=2K{\omega_{\mathrm{r}}}/N\right_{\omega_{\mathrm{r}% }=\bar{\omega}_{\mathrm{r}}}\hat{\omega}_{\mathrm{r}}.  (33) 
The \bar{\left(\cdot\right)}notation indicates the steadystate values of the corresponding operating points. The advantage of this approach is that for each operating point, corresponding steadystate values are substituted in the statespace matrices. This is done by a function \boldsymbol{p}=f(\omega_{\mathrm{r}}(t)):\mathbb{R}\rightarrow\mathbb{R}^{n_{% \mathrm{p}}}, which schedules the system \mathbf{A}(\boldsymbol{p}):\mathbb{R}^{n_{\mathrm{p}}}\rightarrow\mathbb{R}^{n% \times n}, input \mathbf{B}(\boldsymbol{p}):\mathbb{R}^{n_{\mathrm{p}}}\rightarrow\mathbb{R}^{n% \times m} and output \mathbf{C}(\boldsymbol{p}):\mathbb{R}^{n_{\mathrm{p}}}\rightarrow\mathbb{R}^{q% \times n} matrices. This leads to the description of nonlinear dynamics by a set of linear models, varying the system according to the operating point parameterized by \boldsymbol{p}\in\mathcal{P}. For the qLPV case, the scheduling variable is part of the state, which makes the system selfscheduling for each time step. In this paper, a finite number of linearizations is considered for operating conditions along the optimal power coefficient C_{\mathrm{p,max}}(\lambda^{*})=C_{\mathrm{\tau}}(\lambda^{*})\lambda^{*} corresponding to the set \boldsymbol{\mathcal{U}} of belowrated wind speeds.
The current form of the linear model in Equation (28) only describes deviations from the current operating point. To approach the actual states and outputs of the nonlinear model with a qLPV model structure, offsets for the state, input and output should be incorporated in the system description. The process of converting the LPV model to its affine form incorporating these operating point offsets is described in Appendix A. The same appendix also describes the employed fourth order RungeKutta statespace discretization method. When in the remainder of this paper is referred to the qLPV model, the system in its affine form is intended.
3.3 Completing the linearization for the NREL 5MW reference turbine
Description  Symbol  Value  Unit 

Blade inertia  J_{\mathrm{b}}  11.776\cdot 10^{6}  kg m^{2} 
Hub inertia  J_{\mathrm{h}}  115\,926  kg m^{2} 
Total rotor inertia  J_{\mathrm{h}}  35.444\cdot 10^{6}  kg m^{2} 
Torque coefficient fit (1/2)  \theta_{1,2}  \left.14.5924,\,42.7653\right.   
Torque coefficient fit (2/2)  \theta_{3,4}  \left.2.4604,\,0.0036\right.   
Gearbox ratio  N  97   
Air density  \rho_{\mathrm{a}}  1.225  kg m^{3} 
Fine pitch angle  \beta_{0}  1.9\cdot 10^{3}  rad 
Rotor radius  R  63  m 
Optimal mode gain (LSS)  K  2.1286\cdot 10^{6}  Nm (rad s^{1})^{2} 
Optimal tipspeed ratio  \lambda^{*}  7.7   
Input excitation amplitude  a_{\mathrm{u}}  1   
Tower mass  m  1000  kg 
Tower damping  \zeta  100  kg s^{1} 
Tower stiffness  k  500  kg s^{2} 
Tower natural frequency  \omega_{\mathrm{n}}  0.7071  rad s^{1} 
This section provides the data for linearization of the NREL 5MW turbine, and performs a validation of the resulting affine qLPV system to the nonlinear turbine model in highfidelity simulation code. All linearization parameters are summarized in Table (1).
First, an analytical fit is made to the NREL 5MW torque coefficient data as a function of the tipspeed ratio. This is needed as \bar{k}_{\mathrm{\omega_{r}}} and \bar{k}_{\mathrm{U}} are a function of the torque coefficient and its partial derivative with respect to the tipspeed ratio. The torque coefficient data is obtained using a graphical extension^{30} to NREL’s highfidelity wind turbine simulation software FAST v8.16^{31}, which includes blade element momentum (BEM) code^{3} for obtaining rotor characteristic data. As the framework being derived in this paper is focused on the belowrated region, and conventional wind turbine controllers keep the pitch angle fixed at finepitch angle \beta_{0} during partial load^{32}, the dependency of the torque coefficient on \beta is omitted. An often used parameterizable torque coefficient function is defined by
\displaystyle C_{\mathrm{\tau}}(\lambda)  \displaystyle={\rm e}^{\theta_{1}/\lambda}(\theta_{2}/\lambda\theta_{3})/% \lambda+\theta_{4},  (34) 
which is fitted by optimizing the values \theta_{i} using a nonlinear leastsquares routine, minimizing the sumofsquares between the fit and the datapoints. Figure 6 shows the torque coefficient trajectory as a function of the tipspeed ratio for \beta=\beta_{0}, and the fit to this data. Also, an evaluation of the analytically computed partial gradient with respect to the tipspeed ratio is given. Furthermore, the same figure shows the linearization parameters \bar{k}_{\mathrm{\omega_{r}}}, \bar{k}_{\mathrm{U}} and \bar{k}_{\mathrm{\tau_{g}}}. The evaluation is performed for all belowrated rotor speed conditions along the maximum power coefficient C_{\mathrm{p,max}} at an optimal tipspeed ratio of \lambda^{*}=7.7. The trajectories show smooth and linear behavior for all operating points.
In Figure 7, the steadystate values for \bar{q}_{1}, \bar{q}_{2}, \bar{q}_{3}, and \bar{q}_{4} are given as a function of rotor speed for the optimal power coefficient operating conditions. Compared to the previously presented linearization parameters, these trajectories show a more volatile behavior. At the tower natural frequency, two of the trajectories change signs, while the other two reach their extrema. This, as will be shown later, results in some erratic behavior when the qLPV model selfschedules itself around the natural frequency. Therefore, a fine grid of linear models should be taken in the LPV scheduling space, to minimize artifacts and to properly describe the nonlinear dynamics.
3.4 The qLPV model subject to a turbulent wind
The main advantage of a qLPV model, is that the scheduling parameter is part of the state vector. In this way, the scheduling signal is not exogenous, and the model is consequently selfscheduling according to its state evolution. To verify the validity of the derived affine qLPV model, a turbulent wind signal with a mean wind speed of \bar{U}=6.5\,m s^{1} is applied to:

1.
A nonlinear NREL 5MW aerodynamic model simulated in the highfidelity FAST code. A secondorder, first mode tower model G(s) is excited by a unity amplitude cosine as a function of the azimuth position;

2.
The firstorder linearized NREL 5MW wind turbine model with C_{\mathrm{\tau}}(\lambda) lookup table, driving the transformed tower model H(s,\omega_{\mathrm{r}}) by the rotor speed output;

3.
The qLPV model, incorporating both the linear wind turbine rotor and transfomed tower dynamics, selfscheduled by its rotor speed state.
The simulation results, based on the parameters in Table (1), are presented in Figure 8. The left plot shows the rotor speed simulation, which demonstrates that the firstorder and qLPV models accurately follow the FAST output. Subsequently, the right plot compares the tower sideside displacement responses, as a result of the rotor imbalance excitation. As concluded earlier by the frequency sweep in Figure 5, the nominal and transformed tower models show a good match. The additional qLPV response in Figure 8 shows a similar trajectory as the transformed model, apart from some minor artifacts between 700800 s, when the rotor speeds approaches the tower natural frequency. These anomalies are a result of the more erratic behavior of the steadystate gains \bar{q}_{14} in the region of \omega_{\mathrm{n}} (Figure 7), and the switching between the linear models by the scheduling parameter. However, as the response serves as a load indication, and the exact value is of less importance, the qLPV method is concluded being suitable for its intended purpose in the qLPVMPC framework, described in the next section.
4 QuasiLPV model predictive control
Nonlinear MPC is – because of its computational complexity – often considered as an unsuitable control method for application in fast realworld systems, such as wind turbines. Therefore, in this paper, an approach towards an efficient method for nonlinear MPC is employed, exploiting the inherent selfscheduling property of a qLPV system. This section describes the qLPVMPC framework, with the aim to provide a convex QP, defining a tradeoff between maximum power production efficiency, while minimizing tower excitation at its natural frequency. In practice this means that the rotor deviates from maximum power extraction trajectory when it approaches the rotor speed coinciding with the structural resonance frequency.
An economic MPC approach is used to directly optimize for the economic performance of the process ^{33}. For the considered case, that is: a predefined quadratic performance criterion specifies the tradeoff between energy production and load mitigations, and finds the optimal control signal in the prediction horizon resulting in minimization of the criterion. However, as for each time step the scheduling sequence over the prediction horizon is unknown, the nonlinear MPC control problem is solved by an iterative method^{22}. The method solves subsequent QPs minimizing the predefined cost, and uses the resulting predicted scheduling sequence as a warmstart for the next iteration. Each iteration uses a single QP solve. A norm of the consecutive output differences over the prediction horizon is used to determine whether the algorithm has converged.
By manipulation of the affine system representation equations in (43) and (44), an expression is derived for forward propagation of the qLPV model output, only requiring the initial state at time step k and the scheduling sequence over the prediction horizon
\displaystyle\boldsymbol{Y}_{k+1}  \displaystyle=\mathbf{H}(\mathbf{P}_{k})\left(x_{k}\breve{x}(p_{\mathrm{k}})% \right)+\mathbf{S}(\mathbf{P}_{k})\Delta\boldsymbol{U}_{k}(\mathbf{P}_{k})+% \left(\boldsymbol{\breve{Y}}_{k+1}(\mathbf{P}_{k})+\mathbf{L}(\mathbf{P}_{k})% \Delta\boldsymbol{\breve{X}}_{k}(\mathbf{P}_{k})+\mathbf{D}(\mathbf{P}_{k})% \Delta\boldsymbol{U}_{k+1}(\mathbf{P}_{k})\right),  (35) 
in which the matrices \left\{\boldsymbol{H},\,\boldsymbol{S},\,\Delta\boldsymbol{U}_{k},\,% \boldsymbol{\breve{Y}}_{k+1},\,\mathbf{L},\,\Delta\boldsymbol{\breve{X}}_{k},% \,\mathbf{D},\,\Delta\boldsymbol{U}_{k+1}\right\} are defined in Appendix B, and {\mathbf{P}_{\mathrm{k}}=\left[\boldsymbol{p}_{k},\,\boldsymbol{p}_{k+1}~{}% \cdots~{}\boldsymbol{p}_{k+\mathrm{N_{p}}}\right]\in\mathbb{R}^{n_{\mathrm{p}}% \times N_{\mathrm{p}}}} is the collection of scheduling variables at each time instant over the prediction horizon {N_{\mathrm{p}}\in\mathbb{Z}^{+}}. The \breve{(\cdot)}notation indicates steadystate offsets from the current operating point for the the states, in and outputs (Appendix A). The opportunity for defining a control horizon {N_{\mathrm{c}}\in\mathbb{Z}} is disregarded in this paper, and is chosen equal to N_{\mathrm{p}}. For sake of completeness, the above given propagation expression includes a direct feedthrough matrix \mathbf{D}, although it is not used for the considered problem.
At time instance k=0 only the initial state is assumed to be known, and the scheduling parameters are chosen constant over the prediction horizon, such that
\displaystyle\boldsymbol{P_{\mathrm{0}}}  \displaystyle=\boldsymbol{1}_{\mathrm{N_{p}}}\otimes\boldsymbol{p}_{\mathrm{0}},  (36) 
where {\boldsymbol{1}_{\mathrm{N_{p}}}\in\mathbb{R}^{N_{\mathrm{p}}}} is a onedimensional vector of ones. By assuming the initialization vector, the convex QP is solved with {\Delta\boldsymbol{\Theta}_{\mathrm{g},k+1}=\left[\Delta\tau_{\mathrm{g},k+1}% \cdots\Delta\tau_{\mathrm{g},k+N_{\mathrm{p}}}\right]\in\mathbb{R}^{N_{\mathrm% {p}}}} as the decision variable vector, minimizing the cost
\displaystyle\underset{\Delta\boldsymbol{\Theta}_{\mathrm{g},k+1}}{\mathrm{arg% \,min}}  \displaystyle\boldsymbol{Y}_{k+1}^{T}\mathbf{Q}\boldsymbol{Y}_{k+1}+\Delta% \boldsymbol{\Theta}_{\mathrm{g},k+1}^{T}\mathbf{R}\Delta\boldsymbol{\Theta}_{% \mathrm{g},k+1}  (37)  
subject to  \displaystyle\text{Dynamical system in Equation~{}}\eqref{eq:WT_MPC_Equation}, 
in which {\mathbf{Q}=\texttt{diag}(Q,\,Q~{}\cdots~{}Q)\in\mathbb{R}^{\mathrm{N_{p}}% \times\mathrm{N_{p}}}} and {\mathbf{R}=\texttt{diag}(R,\,R~{}\cdots~{}R)\in\mathbb{R}^{\mathrm{N_{p}}% \times\mathrm{N_{p}}}} are, respectively, weight matrices acting on the predicted tower amplitude and deviation from the optimal torque control signal. The latter term of the cost requires the assumption of optimal power production efficiency using the Komegasquared torque control strategy. Now, compare the above given minimization objective with the one introduced in the problem formalization by Equation (3). The first term of Equation (37) aims on fatigue load minimization, whereas the latter term is a combination of energy production maximization and penalization on the control input. Formulating the objective in this way, results in a convenient tradeoff between produced power and load reductions by varying the weight ratio between Q and R.
After the first solve with the initial scheduling sequence of Equation (36), the inherent qLPV property is exploited by using the predicted evolution of the state to form a warmstart initialization of \mathbf{P}_{k}^{j+1} in the next iteration. This iterative process is repeated until \left\left\boldsymbol{Y}_{k}^{j+1}\boldsymbol{Y}_{k}^{j}\right\right_{2}<\epsilon, or for a maximum number of iterations j_{\mathrm{n}}, with \epsilon being a predetermined error threshold. The algorithm is summarized using pseudocode in Algorithm 1.
An evaluation has shown that after convergence during initialization, warmstarting the scheduling sequence for the subsequent timesteps shows excellent results. That is, performing multiple iterations for k>0 shows no significant performance enhancements for the considered problem. Therefore, the described process is only performed in the initial time step k=0. The need for only a single QP for each step, makes the approach for solving the nonlinear MPC problem computationally efficient and tractable for realworld implementations.
5 Highfidelity simulation setup and results
This section implements the proposed qLPVMPC framework in conjunction with the nonlinear NREL 5MW model in the highfidelity FAST code, and the software implementation is made publicly available^{28}. As the sideside natural frequency of the NREL 5MW turbine is located outside the rotor speed operating region, the tower properties are modified. The tower thickness is scaled down by a factor 7.5, such that an effective turbine sideside resonance frequency of approximately 0.7 rad s^{1} is attained. Also, two of the three blades are configured to have an overall mass increase and decrease of 2% with respect to the reference blade. This mass imbalance induces a rotor eccentricity, exacerbating the excitation of the turbine sideside mode.
Furthermore, the simulation environment incorporates the modulated secondorder tower model from Equation (15). The transformed tower model is scheduled by the simulated rotor speed, and the resulting integrator states are, together with the rotor speed, used in each timestep to form the initial state. The initial state is, as shown in Algorithm 1, at each time instant used for forward propagation of the qLPV model by Equation (35).
The aim is now to showcase the framework capabilities of successfully preventing prolonged rotor speed operation near the tower resonance frequency. This is done by defining two separate simulation cases:

•
Case 1: initializing the wind turbine for operating conditions corresponding to a wind speed of U=5.5 m s^{1}, followed by a linearly increasing slope to a maximum wind speed of U=8.0 m s^{1} in approximately 250 s;

•
Case 2: operating the wind turbine in turbulent wind conditions with a mean wind speed \bar{U}=6.5 m s^{1} for 2000 s.
For both cases, the behavior of the qLPVMPC implementation is compared with standard Komegasquared torque control.
The employed wind signals are presented in Figure 9. Because the wind speed cannot assumed to be measurable in realworld scenarios, an effective immersion and invariance (I&I) rotor effective wind speed estimator^{34, 35} is used, which is also plotted in the same figure. Because the future wind speed is unknown at time instance k, the wind speed evolution is chosen to be constant and equal to the current estimated value over the prediction horizon. Also, the smoothened course of the estimated signal aids the qLPVMPC algorithm to prevent from overreacting to rapid variations. As the wind speed estimator takes the applied generator torque and measured rotor speed as inputs, and a rapid rotor speed and generator torque change occurs at this time instant, a discrepancy is seen after 100 s. Nonetheless, the estimator shows a quick recovery in consequent time steps.
Distinct sampling intervals are used for the simulation environment and MPC update actions. Simulation of the NREL 5MW in FAST requires a sampling time of t_{\mathrm{s,FAST}}=0.01 s to prevent numerical issues. The MPC sampling time is set to t_{\mathrm{s,MPC}}=1.0 s. Note that this rather low sampling interval is possible because the modulation transformation moves the load signal to a quasisteady state contribution. As a result of this transformation, the algorithm’s goal is to find the optimal operating trajectory, and not to actively mitigate a specific frequency. The low sampling interval is especially convenient for realworld applications, as this allows solving the QP less frequently, reducing the need for powerful control hardware.
The FAST simulation environment, implemented in MATLAB Simulink^{36}, simulates for t_{\mathrm{s,MPC}}, after which is simulation is paused and essential information is extracted. The data is provided to the MPC algorithm in MATLAB using CVX: a package for specifying and solving convex programs ^{37, 38}. After solving the optimization problem, the first input sample of the decision variable vector is updated in the simulation environment. The simulation is resumed and the input is held constant for the next t_{\mathrm{s,MPC}} seconds, after which it is paused again.
Figure 10, presents the results for simulation Case 1. For this case, the in and output weighing factors are chosen as Q=0.1, R=25, and the prediction horizon is set to N_{\mathrm{p}}=25. The simulation results show the ability of the algorithm to withhold the turbine from operating at a rotational speed exciting the tower natural frequency by increasing \Delta\tau_{\mathrm{g}}. Then, around 100 s, the wind speed is sufficient for the load and power tradeoff to be in favor of the latter mentioned. This is reflected by a swift reduction of the generator torque resulting in a rapid crossing of the critical rotor speed at \omega_{\mathrm{r}}=\omega_{\mathrm{n}}=6.75 rpm. The tower displacement shows a reduction in amplitude by excitation of the natural frequency for a shorter period of time. Obviously, this comes at the expense of generated power and energy.
The simulation results for Case 2 are given in Figure 11. By inspection of the rotor speed around the critical speed, it is being prevented from operating at this frequency for extended time periods. This results in a significant decrease of the tower displacement amplitudes. To further clarify this effect, Figure 12 shows histograms of the rotor speed and amplitude signals. Finally, Figure 13 shows the sidewards displacement spectra of both control strategies. A significant reduction of 18 dB is attained at the turbine sideside natural frequency.
The generator power and torque trajectories of both control strategies show a high degree of similarity, which indicates a minimal penalty on the produced energy production. This is confirmed by the evaluation of the generated energy over the total simulation time, resulting in 603.34 kWh and 601.19 kWh for the respective baseline and MPC cases, which turns out in a negligible produced energy reduction of 0.36 %. The tradeoff is conveniently tuned by the weight ratio between Q and R.
6 Conclusions
In practical scenarios, wind turbine rotors possess a mass and/or aerodynamic imbalance, which cause a periodic excitation possibly coinciding with the tower sideside natural frequency. To date, no efficient and intuitive MPC framework is available for preventing rotor speed operation at this frequency. In this paper, the dynamics of a wind turbine tower are subject to a modulation transformation, and thereby transformed into a quasiLPV system description. The resulting qLPV model, by aggregation with a wind turbine model, is reconciled with an MPC scheme. The combination exploits the inherent properties of the qLPV model, leading to an efficient method of solving a convex optimization problem. The qLPVMPC approach involves finding the qLPV scheduling sequence by performing multiple iterative QP solves for the first time step. Subsequent time steps only require a single QP solve using a scheduling sequence warm start originating from the previous time step. By imposing an additional torque contribution, the rotor speed is prevented from operating near the tower natural frequency at the expense of reduced aerodynamic efficiency. Simulation results with artificial sloped and realistic turbulent wind simulations show that the algorithm prevents excessive natural frequency excitation, by sacrificing an insignificant amount of produced energy. The current work only considers the exclusion of a single excitation frequency, however, the presented framework can be extended towards the exclusion of multiple resonances.
Conflict of interest
This project has been conducted in cooperation with Vestas Wind Systems A/S.
Supporting information
A software implementation is available as part of the online article^{28}.
Appendix A  The affine LPV model representation and discretization
This section presents the process of converting the LPV model derived in Section 3.2 to its affine form. For this, the steadystate offset values of the state, input and output values are saved for each linearization. The offsets are indicated by a \breve{(\cdot)}, and the following relations
\displaystyle\hat{\boldsymbol{q}}(t,\boldsymbol{p}^{*})  \displaystyle=\boldsymbol{q}(t)\breve{\boldsymbol{q}}(\boldsymbol{p}^{*}),  (38)  
\displaystyle\hat{\boldsymbol{u}}(t,\boldsymbol{p}^{*})  \displaystyle=\boldsymbol{u}(t)\breve{\boldsymbol{u}}(\boldsymbol{p}^{*}),  (39)  
\displaystyle\hat{\boldsymbol{y}}(t,\boldsymbol{p}^{*})  \displaystyle=\boldsymbol{y}(t)\breve{\boldsymbol{y}}(\boldsymbol{p}^{*}),  (40) 
are substituted in Equation (28), such that the affine form is obtained
\displaystyle\boldsymbol{\dot{q}}(t)  \displaystyle=\mathbf{A}(\boldsymbol{p})(\boldsymbol{q}(t)\breve{\boldsymbol{% q}}(\boldsymbol{p}^{*}))+\mathbf{B}(\boldsymbol{p})(\boldsymbol{u}(t)\breve{% \boldsymbol{u}}(\boldsymbol{p}^{*}))+\boldsymbol{\dot{\breve{q}}}(\boldsymbol{% p}^{*})  (41)  
\displaystyle\boldsymbol{y}(t)  \displaystyle=\mathbf{C}(\boldsymbol{p})(\boldsymbol{q}(t)\breve{\boldsymbol{% q}}(\boldsymbol{p}^{*}))+\breve{\boldsymbol{y}}(\boldsymbol{p}^{*}),  (42) 
in which \boldsymbol{p}(t)=\boldsymbol{p}^{*}(t) indicates the current linear model in the LPV scheduling space^{39}. Because a finite set of linear models is taken, the scheduling variable might fall between two model scheduling points. In this case, either the nearest offsets corresponding to the current scheduling value are taken, or a linear interpolation is performed. When the models are defined on a fine enough grid, the advantage of increased accuracy by interpolation diminishes, and therefore the nearest model approach is employed.
As this paper uses a samplebased and fixed timestep control setup, the continuoustime system is converted to its discretetime equivalent
\displaystyle\boldsymbol{q}(k+1)  \displaystyle=\mathbf{A}_{\mathrm{d}}(\boldsymbol{p}_{\mathrm{k}})(\boldsymbol% {q}(k)\breve{\boldsymbol{q}}(\boldsymbol{p}_{\mathrm{k}}^{*}))+\mathbf{B}_{% \mathrm{d}}(\boldsymbol{p}_{\mathrm{k}})(\boldsymbol{u}(k)\breve{\boldsymbol{% u}}(\boldsymbol{p}_{\mathrm{k}}^{*}))+\boldsymbol{\breve{q}}(\boldsymbol{p}_{% \mathrm{k}}^{*})  (43)  
\displaystyle\boldsymbol{y}(k)  \displaystyle=\mathbf{C}(\boldsymbol{p}_{\mathrm{k}})(\boldsymbol{q}(k)\breve% {\boldsymbol{q}}(\boldsymbol{p}_{\mathrm{k}}^{*}))+\breve{\boldsymbol{y}}(% \boldsymbol{p}_{\mathrm{k}}^{*}),  (44) 
where k is the discrete timestep variable, and the matrix subscripts (\cdot)_{\mathrm{d}} indicate the discrete time counterparts of the system and input matrices. Discretization of \mathbf{A} and \mathbf{B} is performed using a fourth order RungeKutta discretization method^{40}, of which the matrix transformation relations are given by
\displaystyle\mathbf{A}_{\mathrm{d}}  \displaystyle=\frac{1}{24}\mathbf{A}^{4}t_{\mathrm{s}}^{4}+\frac{1}{6}\mathbf{% A}^{3}t_{\mathrm{s}}^{3}+\frac{1}{2}\mathbf{A}^{2}t_{\mathrm{s}}^{2}+\mathbf{A% }t_{\mathrm{s}}+\mathrm{I_{n}},  (45)  
\displaystyle\mathbf{B}_{\mathrm{d}}  \displaystyle=\frac{1}{24}\mathbf{A}^{3}\mathbf{B}t_{\mathrm{s}}^{4}+\frac{1}{% 6}\mathbf{A}^{2}\mathbf{B}t_{\mathrm{s}}^{3}+\frac{1}{2}\mathbf{A}\mathbf{B}t_% {\mathrm{s}}^{2}+\mathbf{B}t_{\mathrm{s}}+\mathrm{I_{n}}.  (46) 
Note that the last term of Equation (43), originating from the lefthand side of the equation, is in the discretetime case taken at the current time instant, as the output for scheduling the next state offset \boldsymbol{\breve{q}} is unavailable at time step k.
Appendix B  LPV forward propagation matrices
This section defines the LPV forwardpropagation matrices of Equation (35):
\displaystyle\boldsymbol{Y}_{k+1}  \displaystyle=\begin{bmatrix}y_{k+1}\\ y_{k+2}\\ \vdots\\ y_{k+N_{p}}\end{bmatrix},~{}\boldsymbol{\breve{Y}}_{k+1}(\mathbf{P}_{k})=% \begin{bmatrix}\breve{y}(\boldsymbol{p}_{\mathrm{k+1}})\\ \breve{y}(\boldsymbol{p}_{\mathrm{k+2}})\\ \vdots\\ \breve{y}(\boldsymbol{p}_{\mathrm{k+N_{p}}})\end{bmatrix},~{}\Delta\boldsymbol% {U}_{k}(\mathbf{P}_{k})=\begin{bmatrix}u_{k}\breve{u}(\boldsymbol{p}_{\mathrm% {k}})\\ u_{k+1}\breve{u}(\boldsymbol{p}_{\mathrm{k+1}})\\ \vdots\\ u_{k+\mathrm{N_{p}}1}\breve{u}(\boldsymbol{p}_{\mathrm{k+\mathrm{N_{p}}1}})% \end{bmatrix}  
\displaystyle\mathbf{H}(\mathbf{P}_{k})  \displaystyle=\begin{bmatrix}C(\boldsymbol{p}_{k+1})A(\boldsymbol{p}_{k})\\ C(\boldsymbol{p}_{k+2})A(\boldsymbol{p}_{k+1})A(\boldsymbol{p}_{k})\\ \vdots\\ C(\boldsymbol{p}_{\mathrm{k+N_{p}}})A(\boldsymbol{p}_{k+\mathrm{N_{p}}1})% \dots A(\boldsymbol{p}_{k})\end{bmatrix},\quad\mathbf{D}(\mathbf{P}_{k})=% \begin{bmatrix}&D(\boldsymbol{p}_{k+1})&0&\cdots&0\\ &0&D(\boldsymbol{p}_{k+2})&\cdots&0\\ &\vdots&&\ddots&\vdots\\ &0&\cdots&&D(\boldsymbol{p}_{k+N_{p}})\end{bmatrix},  
\displaystyle\mathbf{S}(\mathbf{P}_{k})  \displaystyle=\begin{bmatrix}&C(\boldsymbol{p}_{k+1})B(\boldsymbol{p}_{k})&0&% \dots&0\\ &C(\boldsymbol{p}_{k+2})A(\boldsymbol{p}_{k+1})B(\boldsymbol{p}_{k})&C(% \boldsymbol{p}_{k+2})B(\boldsymbol{p}_{k+1})&\dots&0\\ &\vdots&&\ddots&\vdots\\ &C(\boldsymbol{p}_{\mathrm{k+N_{p}}})A(\boldsymbol{p}_{\mathrm{k+N_{p}}1})% \dots A(\boldsymbol{p}_{k+1})B(\boldsymbol{p}_{k})&\cdots&&C(\boldsymbol{p}_{% \mathrm{k+N_{p}}})B(\boldsymbol{p}_{\mathrm{k+N_{p}}1})\end{bmatrix},  
\displaystyle\mathbf{L}(\mathbf{P}_{k})  \displaystyle=\begin{bmatrix}&C(\boldsymbol{p}_{k+1})&0&\dots&0\\ &C(\boldsymbol{p}_{k+2})A(\boldsymbol{p}_{k+1})&C(\boldsymbol{p}_{k+2})&\dots&% 0\\ &\vdots&&\ddots&\vdots\\ &C(\boldsymbol{p}_{\mathrm{k+N_{p}}})A(\boldsymbol{p}_{k+\mathrm{N_{p}}1})% \dots A(\boldsymbol{p}_{k+1})&\cdots&&C(\boldsymbol{p}_{\mathrm{k+N_{p}}})\end% {bmatrix},\quad\Delta\boldsymbol{\breve{X}}_{k}(\mathbf{P}_{k})=\begin{bmatrix% }\breve{x}(\boldsymbol{p}_{\mathrm{k}})\breve{x}(\boldsymbol{p}_{\mathrm{k+1}% })\\ \breve{x}(\boldsymbol{p}_{\mathrm{k+1}})\breve{x}(\boldsymbol{p}_{\mathrm{k+2% }})\\ \vdots\\ \breve{x}(\boldsymbol{p}_{\mathrm{k+N_{p}1}})\breve{x}(\boldsymbol{p}_{% \mathrm{k+N_{p}}})\end{bmatrix}, 
References
 1 Hau E. Wind turbines: fundamentals, technologies, application, economics. Berlin, Germany: Springer Science & Business Media . 2013.
 2 GL Renewables Certification. Guideline for the Certification of Offshore Wind Turbines. 2012.
 3 Burton T, Jenkins N, Sharpe D, Bossanyi EA. Wind energy handbook. Chichester, United Kingdom: John Wiley & Sons . 2001.
 4 Bossanyi EA. Wind turbine control for load reduction. Wind Energy 2003; 6(3): 229–244.
 5 Bossanyi EA, Fleming P, Wright A. Field test results with individual pitch control on the NREL CART3 wind turbine: 1019–1028; American Institute of Aeronautics and Astronautics (AIAA) . 2012.
 6 Wright A, Fleming P, van Wingerden JW. Refinements and tests of an advanced controller to mitigate fatigue loads in the controls advanced research turbine: 815–827; American Institute of Aeronautics and Astronautics (AIAA) . 2011.
 7 Bossanyi EA. The design of closed loop controllers for wind turbines. Wind Energy 2000; 3(3): 149163.
 8 Bianchi FD, De Battista H, Mantz RJ. Wind turbine control systems: principles, modelling and gain scheduling design. Springer Science & Business Media . 2006.
 9 Mulders SP, van Wingerden JW. Delft Research Controller: an opensource and communitydriven wind turbine baseline controller. Journal of Physics: Conference Series 2018; 1037(3).
 10 Rawlings JB, Mayne DQ. Model predictive control: Theory and design. Nob Hill Pub. Madison, Wisconsin . 2009.
 11 Hovgaard TG, Boyd S, Jørgensen JB. Model predictive control for wind power gradients. Wind Energy 2015; 18(6).
 12 Holkar K, Waghmare L. An overview of model predictive control. International Journal of Control and Automation 2010; 3(4): 47–63.
 13 Schlipf D, Schlipf DJ, Kühn M. Nonlinear model predictive control of wind turbines using LIDAR. Wind energy 2013; 16(7): 1107–1129.
 14 Evans MA, Cannon M, Kouvaritakis B. Robust MPC tower damping for variable speed wind turbines. IEEE Transactions on Control Systems Technology 2015; 23(1): 290–296.
 15 Tofighi E, Schlipf D, Kellett CM. Nonlinear model predictive controller design for extreme load mitigation in transition operation region in wind turbines. In: 2015 Conference on Control Applications (CCA). IEEE. ; 2015: 1167–1172.
 16 Hours JH, Zeilinger MN, Gondhalekar R, Jones CN. Constrained spectrum control. IEEE Transactions on Automatic Control 2015; 60(7): 1969–1974.
 17 Jain A, Biyik E, Chakrabortty A. A model predictive control design for selective modal damping in power systems. In: 2015 American Control Conference (ACC). IEEE. ; 2015: 4314–4319.
 18 Keyvani A, Sadeghian H, van Wingerden JW, Tamer M, Goosen J, Keulen F. A Comprehensive Model for Transient Behavior of Tapping Mode Atomic Force Microscope. Nonlinear Dynamics 2019. Manuscript submitted for publication.
 19 Bossanyi EA. Individual blade pitch control for load reduction. Wind energy 2003; 6(2): 119–128.
 20 Geyler M, Caselitz P. Individual blade pitch control design for load reduction on large wind turbines. In: 2007 European Wind Energy Conference (EWEC). European Wind Energy Conference (EWEC). ; 2007.
 21 Hanema J. Anticipative model predictive control for linear parametervarying systems. PhD thesis. Technische Universiteit Eindhoven, Eindhoven, The Netherlands; 2018.
 22 Cisneros PS, Voss S, Werner H. Efficient nonlinear model predictive control via quasiLPV representation. In: 55th Conference on Decision and Control (CDC). IEEE. ; 2016: 3216–3221.
 23 Hansen MH. Aeroelastic stability analysis of wind turbines using an eigenvalue approach. Wind Energy 2004; 7(2): 133143.
 24 Dahleh M, Dahleh MA, Verghese G. Lectures on Dynamic Systems and Control  Chapter 12  Modal Decomposition of StateSpace Models. https://ocw.mit.edu/courses/electricalengineeringandcomputerscience/6241jdynamicsystemsandcontrolspring2011/readings/MIT6_241JS11_chap12.pdf; 2011. [Online; accessed 07May2019].
 25 Skogestad S, Postlethwaite I. Multivariable feedback control: analysis and design. New York, United States: John Wiley & Sons . 2007.
 26 Oppenheim A, Willsky A, Nawab S. Signals and Systems. London, United Kingdom: Pearson . 2013.
 27 MathWorks . Symbolic Math Toolbox. https://mathworks.com/products/symbolic.html; 2019. [Online; accessed 07May2019].
 28 Mulders SP, Hovgaard TG, Grunnet JD, van Wingerden JW. Software implementation: Preventing wind turbine tower natural frequency excitation with a quasiLPV model predictive control scheme (Version 1.0). Zenodo. http://doi.org/10.5281/zenodo.2672956; 2019. [Online; accessed 07May2019].
 29 Jonkman J, Butterfield S, Musial W, Scott G. Definition of a 5MW reference wind turbine for offshore system development. tech. rep., National Renewable Energy Lab (NREL); Golden, Colorado, United States: 2009.
 30 Bos R, Zaaijer M, Mulders SP, van Wingerden JW. FASTTool. https://github.com/TUDelftDataDrivenControl/FASTTool; 2019. [Online; accessed 07May2019].
 31 NREL  NWTC . FAST v8.16. https://nwtc.nrel.gov/FAST8; 2019. [Online; accessed 07May2019].
 32 Pao LY, Johnson KE. Control of wind turbines. IEEE Control Systems Magazine 2011; 31(2): 44–62.
 33 Rawlings JB, Angeli D, Bates CN. Fundamentals of economic model predictive control. In: 51st Conference on Decision and Control (CDC). IEEE. ; 2012: 3851–3861.
 34 Ortega R, MancillaDavid F, Jaramillo F. A globally convergent wind speed estimator for wind turbine systems. International Journal of Adaptive Control and Signal Processing 2013; 27(5): 413–425.
 35 Soltani MN, Knudsen T, Svenstrup M, et al. Estimation of rotor effective wind speed: A comparison. IEEE Transactions on Control Systems Technology 2013; 21(4): 1155–1167.
 36 MathWorks . Simulink: Simulation and ModelBased Design. https://nl.mathworks.com/products/simulink.html; 2019. [Online; accessed 07May2019].
 37 Grant M, Boyd S. CVX: Matlab Software for Disciplined Convex Programming, version 2.1. http://cvxr.com/cvx; 2014. [Online; accessed 07May2019].
 38 Grant M, Boyd S. Graph implementations for nonsmooth convex programs. In: Lecture Notes in Control and Information Sciences. SpringerVerlag Limited. 2008 (pp. 95–110). http://stanford.edu/~boyd/graph_dcp.html.
 39 MathWorks . Linear ParameterVarying Models. https://mathworks.com/help/control/ug/linearparametervaryingmodels.html; 2019. [Online; accessed 07May2019].
 40 Shampine LF, Allen RC, Pruess S. Fundamentals of numerical computing. New York, United States: John Wiley & Sons . 1997.