Preventing wind turbine tower natural frequency excitation with a quasi-LPV model predictive control schemeShort title: A qLPV-MPC framework for preventing tower natural frequency excitation

Preventing wind turbine tower natural frequency excitation with a quasi-LPV model predictive control schemethanks: Short title: A qLPV-MPC framework for preventing tower natural frequency excitation

[    [    [    [ \orgdivDelft Center for Systems and Control, Faculty of Mechanical Engineering, \orgnameDelft University of Technology, \orgaddress\countryThe Netherlands \orgdivVestas Power Solutions, \orgnameVestas Wind Systems A/S, \orgaddress\countryDenmark

[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 rotor-speed dependent tower side-side excitation during below-rated operation. Operating at this frequency generally degrades the expected structural life span, as the first tower natural frequency for larger turbines coincides with a below-rated 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 trade-off between energy maximization and fatigue load minimization. Therefore, this paper introduces a quasi-linear parameter varying model predictive control (qLPV-MPC) 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 5-MW reference wind turbine in high-fidelity simulation code. Prolonged rotor speed operation at the tower natural frequency is prevented, whereas when the trade-off 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.

tower natural frequency skipping, model predictive control, quasi-linear parameter varying, model modulation transformation

Research Article

1]Sebastiaan Paul Mulders 2]Tobias Gybel Hovgaard 2]Jacob Deleuran Grunnet 1]Jan-Willem 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 variable-speed below-rated operation, the blade-passing frequency may excite the structural side-side natural frequency. Because the turbine rotor provides negligible aerodynamic side-side damping, at an order of magnitude smaller than the fore-aft damping ratio, small perturbations can lead to load fluctuations comparable to fore-aft stresses 3. As a result, excitation of the side-side mode possibly results in accelerated and accumulative fatigue damage.

Straightforward control implementations are available for reducing and mitigating the excitation of the tower fore-aft/side-side 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 pitch4, 5 or generator torque6 control signal, for respective damping of fore-aft and side-side vibrations. Another, more passive method, entails prevention of structural mode excitation by increasing the generator torque when the rotor speed approaches the excitation frequency7. The method is often referred to as frequency skipping by inclusion of speed exclusion zones.

All of the above-described active and passive methods complicate the controller design, by requiring extra proportional-integral-derivative (PID) feedback control loops with additional logic and speed set-points. Also, the methods do not incorporate convenient and inherent tuning capabilities for a trade-off 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 benefits8, to the authors’ knowledge, more sophisticated control methods do not yet see a wide-spread adoption in industrial-grade wind turbine control systems; PID control structures9 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 control12 are (1) the ability of including constraints, (2) coping with the complexity of non-minimum phase systems, (3) robustness against deviations of the control model to the actual process, and (4) the convenient application to multi-variable control problems. MPC has been considered for wind turbine load mitigations in the literature. A non-linear MPC (NMPC) method is applied by assuming future wind speed knowledge using a light-detection and ranging (LIDAR) system13. 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 fore-aft damping14. 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 prediction15 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 trade-off 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 non-convex problems is – because of its computational complexity – often considered unsuitable for real-world applications.

Imposing spectral constraints might form a possible solution path16, by employing the short-time Fourier transform (STFT) on the system output signal in a non-linear MPC setting. A similar methodology17 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 non-trivial. 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 Mode-Atomic Force Microscope (TM-AFM)18. The model transformation is applied to the turbine tower model, and transfers frequency-dependent magnitude and phase content to a quasi steady-state contribution. This is accomplished by converting a linear-time invariant (LTI) system into a linear-parameter varying (LPV) model, scheduled on the excitation frequency. The technique shows similarities with the Multi-Blade Coordinate (MBC) transformation, often used in individual pitch control (IPC) implementations for blade fatigue load reductions19, 20.

An LPV system representation is frequently used for capturing non-linear dynamics into a system description with a linear input-output mapping21. An external scheduling variable varies the dynamics of the linear model. Now, consider the combination of an LPV model with MPC. The model-based 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 quasi-LPV (qLPV) system. Recently, an efficient MPC scheme for such qLPV systems is proposed22 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 systems22. Using the qLPV-MPC framework, a methodology for performing an optimal trade-off 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 steady-state contribution

  • Applying the transformation to a second-order 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 below-rated operating points, for obtaining a qLPV state-space 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 closed-loop high-fidelity 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 user-defined trade-off between tower loads and generated energy for the prevailing environmental conditions. In Section 5, the qLPV-MPC framework is evaluated with high-fidelity simulations using the NREL 5-MW 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 trade-off, a wind turbine model needs to be combined with a structural tower model. Section 2.1, describes the tower side-side dynamics by a second-order mass-damper-spring system. Section 2.2 formalizes the problem statement, and shows that the combined system description results in non-convexity. 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 second-order system

Figure 1: A rotor imbalance excites the turbine support structure due to the centripetal force F=a_{\mathrm{u}}\cos{(\psi(t))} at the once-per-revolution and rotor-speed dependent (1P) frequency. The tangential speed of the imbalance is denoted as v_{\mathrm{t}}, and the side-side nacelle displacement x is given in the hub coordinate 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 large-scale state-of-the-art wind turbines are operated with a variable-speed control strategy for below-rated conditions, the support structure is excited by a periodic and frequency-varying centripetal force, as illustrated in Figure 1. The tower dynamics, excited by a rotor-speed dependent once-per-revolution (1P) periodic force, are modeled by a second-order mass-damper-spring 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 side-side displacement, velocity and acceleration in the hub coordinate system, illustrated in Figure 1. A second-order system is taken, as this represents the first mode of the tower using the well known modal-decomposition model reduction technique23, 24, and allows for a convenient assessment and derivation of the modulation transformation in the next section: the application to higher-order 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 first-order differential equations by defining x_{1}=\dot{x}(t), x_{2}=x(t), such that it is rewritten in the standard state-space \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}[]{c|c}\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 state-space 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


Figure 2: In the upper diagram, a wind turbine model is driven by wind disturbance and torque control inputs, and has energy production and azimuth position as outputs. The latter mentioned output is taken by a cosine function, and serves as an input to the tower model, resulting in a fatigue load signal. The presence of the trigonometric function forms a barrier for describing the energy-load trade-off as a convex optimization problem. Therefore, in the lower diagram, the trigonometric function is combined with the tower model by a model modulation transformation, resulting in an LPV system description. Joining the wind turbine model with the transformed tower model provides a quasi-LPV system description, as the rotor speed state serves as the scheduling variable.

This section formalizes the problem considered in this paper. The aim is to provide a trade-off 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 energy-load trade-off

\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 trade-offs. 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 rotor-speed 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 non-linear 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 quasi-LPV 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 wave26. 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 steady-state 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 of18 and only the main results are given in this section. The derivation is performed and validated with symbolic manipulation software for algebraic expressions27, of which the code is made publicly available28. 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 slow-varying 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 steady-state 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 time-scales, 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)
Figure 3: Left: the nominal tower model is periodically excited at a certain frequency and amplitude. For the linear case, the response is scaled and phase-shifted with respect to the driving input signal. Right: after the modulation transformation, the input amplitude is a direct input to the system, whereas its frequency changes the system dynamics. The resulting outputs give the response amplitude and phase shift as a quasi-steady state signal.

Term-by-term 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}[]{c|c}\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 state-space 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 second-order 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: Left: frequency response of the nominal tower model G(\mathrm{j}\omega). A clear resonance peak is observed at \omega_{\mathrm{n}}\approx 0.7 rad s-1, and a -40 dB/decade roll-off at higher frequencies. A set of 4 comparison points \boldsymbol{\omega_{\mathrm{r}}}=\left\{0,\,0.5,\,0.7,\,2.0\right\} rad s-1 is chosen for evaluation of the nominal and modulated model. Right: frequency responses of the transformed model H(\mathrm{j}\omega,\omega_{\mathrm{r},i}) for the set of comparison points: the magnitude content at the indicated frequencies in the left plot is transferred to a steady-state contribution in the transformed case. When the input signal a_{\mathrm{u}} to the transformed model is considered constant or slowly varying, the additional resonances at higher frequencies do not contribute to the output.
Figure 5: Frequency sweep applied to the nominal and transformed model, from \omega_{\mathrm{r}}=0 to 1.2 rad s-1 with a constant acceleration in 1200 s. The transformed model shows a very close amplitude tracking of the nominal model magnitude response.

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 low-pass or notch filters, are dispensable.

The effect of the transformation in the time-domain 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 large-scale variable-speed 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 5-MW (linear) model, for augmentation to the modulated tower model such that a quasi-LPV model is obtained. Section 3.1 provides the simple first-order 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 5-MW reference wind turbine29. Finally, Section 3.4 validates the first-order and affine linear models to simulation results of their non-linear 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 first-order 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 tip-speed 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 K-omega-squared 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 low-speed 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 above-given system description contains the non-linear aerodynamic and generator torque input defined previously by Equations (20) and (21). Furthermore, the output a_{\mathrm{y}} is a non-linear 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 state-space 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)


\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)U-c_{% \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 K-omega-squared 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 steady-state values of the corresponding operating points. The advantage of this approach is that for each operating point, corresponding steady-state values are substituted in the state-space 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 non-linear 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 self-scheduling 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 below-rated 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 non-linear 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 Runge-Kutta state-space 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 5-MW reference turbine

Table 1: Parameters for linearization and simulation of the qLPV model in the below-rated operating region.
Description Symbol Value Unit
Blade inertia J_{\mathrm{b}} 11.776\cdot 10^{6} kg m2
Hub inertia J_{\mathrm{h}} 115\,926 kg m2
Total rotor inertia J_{\mathrm{h}} 35.444\cdot 10^{6} kg m2
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 tip-speed 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
Figure 6: Left: torque coefficient curve of the NREL 5-MW wind turbine as a function of the dimensionless tip-speed ratio. The fit according to the model structure proposed by Equation (34) shows a close fit to the data points. The fit allows the derivation and evaluation of the partial gradient. Right: the linearization parameters defining the LPV model at each scheduling instant.

This section provides the data for linearization of the NREL 5-MW turbine, and performs a validation of the resulting affine qLPV system to the non-linear turbine model in high-fidelity simulation code. All linearization parameters are summarized in Table (1).

First, an analytical fit is made to the NREL 5-MW torque coefficient data as a function of the tip-speed 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 tip-speed ratio. The torque coefficient data is obtained using a graphical extension30 to NREL’s high-fidelity wind turbine simulation software FAST v8.1631, which includes blade element momentum (BEM) code3 for obtaining rotor characteristic data. As the framework being derived in this paper is focused on the below-rated region, and conventional wind turbine controllers keep the pitch angle fixed at fine-pitch angle \beta_{0} during partial load32, 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 non-linear least-squares routine, minimizing the sum-of-squares between the fit and the data-points. Figure 6 shows the torque coefficient trajectory as a function of the tip-speed ratio for \beta=\beta_{0}, and the fit to this data. Also, an evaluation of the analytically computed partial gradient with respect to the tip-speed 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 below-rated rotor speed conditions along the maximum power coefficient C_{\mathrm{p,max}} at an optimal tip-speed ratio of \lambda^{*}=7.7. The trajectories show smooth and linear behavior for all operating points.

In Figure 7, the steady-state 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 self-schedules 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 non-linear dynamics.

Figure 7: State-state gains \bar{q}_{1} to \bar{q}_{4} as a function of the rotor speed scheduling variable. Around the natural tower frequency, the gains show a more erratic behavior raising the need for an LPV model set on a fine scheduling grid.

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 self-scheduling 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. 1.

    A non-linear NREL 5-MW aerodynamic model simulated in the high-fidelity FAST code. A second-order, first mode tower model G(s) is excited by a unity amplitude cosine as a function of the azimuth position;

  2. 2.

    The first-order linearized NREL 5-MW wind turbine model with C_{\mathrm{\tau}}(\lambda) look-up table, driving the transformed tower model H(s,\omega_{\mathrm{r}}) by the rotor speed output;

  3. 3.

    The qLPV model, incorporating both the linear wind turbine rotor and transfomed tower dynamics, self-scheduled 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 first-order and qLPV models accurately follow the FAST output. Subsequently, the right plot compares the tower side-side 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 700-800 s, when the rotor speeds approaches the tower natural frequency. These anomalies are a result of the more erratic behavior of the steady-state gains \bar{q}_{1-4} 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 qLPV-MPC framework, described in the next section.

Figure 8: Simulation results showing the rotor speed and tower displacement ampltide for different models driven by a turbulent wind disturbance input signal. The FAST model excites the nominal tower model and serves as a baseline simulation case. The first-order wind turbine model is combined with the non-linearized transformed tower model. The qLPV system is a single scheduled system description of both the rotor and tower dynamics, and shows – apart from minor artifacts around the tower natural frequency – a close resemblance to its non-linear companion.

4 Quasi-LPV model predictive control

Non-linear MPC is – because of its computational complexity – often considered as an unsuitable control method for application in fast real-world systems, such as wind turbines. Therefore, in this paper, an approach towards an efficient method for non-linear MPC is employed, exploiting the inherent self-scheduling property of a qLPV system. This section describes the qLPV-MPC framework, with the aim to provide a convex QP, defining a trade-off 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 trade-off 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 non-linear MPC control problem is solved by an iterative method22. The method solves subsequent QPs minimizing the predefined cost, and uses the resulting predicted scheduling sequence as a warm-start 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 steady-state 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 one-dimensional 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 K-omega-squared 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 trade-off 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 warm-start 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.

Algorithm 1 - Pseudocode for iteratively finding the scheduling vector \mathbf{P}_{\mathrm{k}} in the first time step, and warm-starting for subsequent time instances.
k\leftarrow 0, j\leftarrow 1, j_{\mathrm{n}}\leftarrow 5
Define Q, R, N_{\mathrm{p}}
Initialize matrices \mathbf{Q}, \mathbf{R}
Initialize state \mathbf{X}, output \mathbf{Y}, and scheduling \boldsymbol{P} matrices as empty \mathbf{0}-matrices
\mathbf{P}^{j}_{k}\leftarrow 1_{\mathrm{N_{p}}}\otimes f(\boldsymbol{x}_{0}), \mathbf{X}(:,k)=\boldsymbol{x}_{0}
for time instance k do
     while j\leq j_{\mathrm{n}} do
         Construct matrices \mathbf{H}(\mathbf{P}^{j}_{k}), \mathbf{S}(\mathbf{P}^{j}_{k}), \mathbf{L}(\mathbf{P}^{j}_{k}), \Delta\boldsymbol{U}(\mathbf{P}^{j}_{k}), \boldsymbol{\breve{Y}}_{k+1}(\mathbf{P}^{j}_{k}), \Delta\boldsymbol{\breve{X}}_{k+1}(\mathbf{P}^{j}_{k})
         Solve for \Delta\boldsymbol{\Theta}_{\mathrm{g}} as in Equation (37) with \mathbf{X}(:,k) as initial state
         Simulate the qLPV model with \Delta\boldsymbol{\Theta}_{\mathrm{g}} for N_{\mathrm{p}} samples to find the state evolution \boldsymbol{\mathcal{X}}^{j}
         Define \mathbf{P}^{j+1}_{k}=f(\boldsymbol{\mathcal{X}}^{j})
         j\leftarrow j+1
     end while
     j\leftarrow 1, j_{\mathrm{n}}\leftarrow 1
     Take the first sample of \Delta\boldsymbol{\Theta}_{\mathrm{g}} to apply in high-fidelity code and simulate for t_{\mathrm{s,FAST}}
     Save resulting state and output data: \mathbf{X}(:,k+1)\leftarrow\boldsymbol{x}_{k}, \mathbf{Y}(:,k)\leftarrow\boldsymbol{y}_{k}
     Define \mathbf{P}^{j}_{k+1}=f(\boldsymbol{\mathcal{X}}^{\mathrm{end}}) as a warm start for the next time instance
     k\leftarrow k+1
end for

An evaluation has shown that after convergence during initialization, warm-starting the scheduling sequence for the subsequent time-steps 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 non-linear MPC problem computationally efficient and tractable for real-world implementations.

Figure 9: A linearly increasing slope and turbulent wind profile employed for the two simulation cases. The rotor effective wind speed is estimated by a wind speed estimator. A discrepancy of the estimated sloped wind speed is observed after 100 s, which is a result of sudden changes in applied generator torque and measured rotor speed.

5 High-fidelity simulation set-up and results

This section implements the proposed qLPV-MPC framework in conjunction with the non-linear NREL 5-MW model in the high-fidelity FAST code, and the software implementation is made publicly available28. As the side-side natural frequency of the NREL 5-MW 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 side-side 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 side-side mode.

Furthermore, the simulation environment incorporates the modulated second-order 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 time-step 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 qLPV-MPC implementation is compared with standard K-omega-squared torque control.

Figure 10: Simulation Case 1 shows a comparison to conventional torque control with a linearly increasing wind speed. The proposed algorithm prevents the rotor speed from prolonged operation at the tower’s natural frequency by imposing an additional generator torque demand. Then, when the wind speed is sufficient for operation at a higher rotor speed, the additional generator torque is rapidly reduced to facilitate a swift crossing of the critical frequency. The strategy is beneficial for reducing periodic tower loading, at the expense of generated power.
Figure 11: Simulation Case 2 shows a comparison to conventional torque control with a realistic turbulent wind profile. It is shown that the tower loading extremes are significantly reduced by preventing prolonged operation at the critical rotor speed. The algorithm shows to have minimal impact on the generated power.
Figure 12: Histograms of the rotor speed and displacement amplitude occurrence for simulation Case 2. The rotor speed histogram (left) clearly shows that the qLPV-MPC algorithm prevents operation at the critical speed. Consequently, the amplitude histogram (right) shows a reduced maximum occurrence, whereas smaller amplitudes happen more frequently, which is beneficial from a fatigue loading viewpoint.

The employed wind signals are presented in Figure 9. Because the wind speed cannot assumed to be measurable in real-world scenarios, an effective immersion and invariance (I&I) rotor effective wind speed estimator34, 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 qLPV-MPC algorithm to prevent from over-reacting 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 5-MW 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 quasi-steady 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 real-world applications, as this allows solving the QP less frequently, reducing the need for powerful control hardware.

The FAST simulation environment, implemented in MATLAB Simulink36, 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 trade-off 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 side-side 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 trade-off is conveniently tuned by the weight ratio between Q and R.

Figure 13: Side-side displacement spectra of the NREL 5-MW turbine subject to high-fidelity turbulent wind simulations. The tower steel thickness is modified, such that a turbine side-side natural frequency of approximately 0.7 rad s-1 is obtained. The blades are given a dissimilar mass to induce rotor eccentricity. Spectra for the K-omega-squared and qLPV-MPC implementations show a significant reduction of the dominant resonance. The content in the dashed box is enlarged in the right plot.

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 side-side 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 quasi-LPV 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 qLPV-MPC 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 article28.

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 steady-state 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 space39. 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 sample-based and fixed time-step control set-up, the continuous-time system is converted to its discrete-time 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 time-step 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 Runge-Kutta discretization method40, 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 left-hand side of the equation, is in the discrete-time 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 forward-propagation 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},


  • 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): 149-163.
  • 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 open-source and community-driven 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 parameter-varying systems. PhD thesis. Technische Universiteit Eindhoven, Eindhoven, The Netherlands; 2018.
  • 22 Cisneros PS, Voss S, Werner H. Efficient nonlinear model predictive control via quasi-LPV 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): 133-143.
  • 24 Dahleh M, Dahleh MA, Verghese G. Lectures on Dynamic Systems and Control - Chapter 12 - Modal Decomposition of State-Space Models.; 2011. [Online; accessed 07-May-2019].
  • 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.; 2019. [Online; accessed 07-May-2019].
  • 28 Mulders SP, Hovgaard TG, Grunnet JD, van Wingerden JW. Software implementation: Preventing wind turbine tower natural frequency excitation with a quasi-LPV model predictive control scheme (Version 1.0). Zenodo.; 2019. [Online; accessed 07-May-2019].
  • 29 Jonkman J, Butterfield S, Musial W, Scott G. Definition of a 5-MW 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.; 2019. [Online; accessed 07-May-2019].
  • 31 NREL - NWTC . FAST v8.16.; 2019. [Online; accessed 07-May-2019].
  • 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, Mancilla-David 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 Model-Based Design.; 2019. [Online; accessed 07-May-2019].
  • 37 Grant M, Boyd S. CVX: Matlab Software for Disciplined Convex Programming, version 2.1.; 2014. [Online; accessed 07-May-2019].
  • 38 Grant M, Boyd S. Graph implementations for nonsmooth convex programs. In: Lecture Notes in Control and Information Sciences. Springer-Verlag Limited. 2008 (pp. 95–110).
  • 39 MathWorks . Linear Parameter-Varying Models.; 2019. [Online; accessed 07-May-2019].
  • 40 Shampine LF, Allen RC, Pruess S. Fundamentals of numerical computing. New York, United States: John Wiley & Sons . 1997.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description