Preventing wind turbine tower natural frequency excitation with a quasi-LPV model predictive control scheme

TL;DR

This paper presents a quasi-LPV model predictive control scheme to prevent wind turbine tower excitation at natural frequencies, balancing energy output and structural fatigue in flexible, tall towers.

Contribution

Introduction of a quasi-LPV-MPC method that enables effective trade-offs between energy maximization and fatigue load reduction for wind turbines.

Findings

Prevents rotor speed from operating at tower natural frequency.

Enables rapid passage over natural frequency for higher energy production.

Effective in both synthetic and realistic high-fidelity simulations.

Abstract

With the ever increasing power rates of wind turbines, more advanced control techniques are needed to facilitate tall towers that are low-weight and cost effective, but in effect more flexible. Such soft-soft tower configurations generally have their fundamental side-side frequency in the below-rated operational domain. Because the turbine rotor practically has or develops a mass imbalance over time, a periodic and rotor-speed dependent side-side excitation is present during below-rated operation. Persistent operation at the coinciding tower and rotational frequency degrades the expected structural life span. 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 demodulation transformation, and is subsequently 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 side-side natural frequency is prevented, whereas when the trade-off is in favor of energy production, the algorithm decides to rapidly pass over the natural frequency to attain higher rotor speeds and power productions.

Tables

Table 1. Table 1: Parameters for linearization and simulation of the qLPV model in the below-rated operating region.

Description	Symbol	Value	Unit
Blade inertia	$J_b$	$11.776\cdot {10}^6$	kg m2
Hub inertia	$J_h$	$\mathrm{115\hspace{0.17em}926}$	kg m2
Total rotor inertia	$J_h$	$35.444\cdot {10}^6$	kg m2
Torque coefficient fit (1/2)	${\theta }_{1,2}$	$14.5924,\mathrm{\hspace{0.17em}42.7653}$	-
Torque coefficient fit (2/2)	${\theta }_{3,4}$	$2.4604,\mathrm{\hspace{0.17em}0.0036}$	-
Gearbox ratio	$N$	$97$	-
Air density	${\rho }_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 }^{\ast }$	$7.7$	-
Input excitation amplitude	$a_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 }_n$	$0.7071$	rad s-1

Equations

$$m\ddot{x}(t)+\zeta \dot{x}(t)+kx(t)$$

$$\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)\stackrel{{\scriptstyle s}}{{=}}\left[\begin{array}[]{c|c}\mathbf{A_{\mathrm{g}}}&\mathbf{B_{\mathrm{g}}}\\ \hline\cr\mathbf{C_{\mathrm{g}}}&\mathbf{0}\end{array}\right],$$

$$\underset{Torque}{\mathrm{argmin}}$$

$${\int }_0^{T_{\mathrm{r}}}f(\tau )g(t)dt$$

$$x_i(t)$$

$${\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\}.$$

$$\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,$$

$$\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.$$

$${\int }_0^{T_{\mathrm{r}}}\Re \left\{C{e}^{j\theta }\right\}{e}^{j\theta }d\theta =0,$$

$${\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,$$

$${\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.$$

$$\left[\begin{array}{c}{\dot{X}}_1\\ {\dot{X}}_2\end{array}\right]$$

$$\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}})\stackrel{{\scriptstyle s}}{{=}}\left[\begin{array}[]{c|c}\mathbf{A_{\mathrm{h}}}&\mathbf{B_{\mathrm{h}}}\\ \hline\cr\mathbf{C_{\mathrm{h}}}&\mathbf{0}\end{array}\right],$$

$$a_{\mathrm{y}}(\tau )$$

$$\varphi (\tau )$$

$$H(s,{\omega }_{\mathrm{r}})$$

$$J_{\mathrm{r}}{\dot{\omega }}_{\mathrm{r}}$$

$${\tau }_{\mathrm{a}}$$

$${\tau }_{\mathrm{g}}=K{\omega }_{\mathrm{r}}^2/N,$$

$$K=\frac{\pi {\rho }_{\mathrm{a}}{R}^5{C}_{\mathrm{p}}(\lambda ,\beta )}{2{\lambda }^3},$$

$$\left[\begin{array}{c}{\dot{q}}_1\\ {\dot{q}}_2\\ {\dot{q}}_3\\ {\dot{q}}_4\\ {\dot{\omega }}_{\mathrm{r}}\end{array}\right]$$

$$a_{\mathrm{y}}$$

$$\hat{\mathbf{q}}(t)$$

$$\hat{\mathbf{u}}(t)$$

$$\hat{\mathbf{y}}(t)$$

$$\dot{\hat{\mathbf{q}}}(t)$$

$$\hat{\mathbf{y}}(t)$$

$$\mathbf{A}(\mathbf{p})$$

$${\hat{\tau }}_{\mathrm{a}}$$

$${\bar{k}}_{{\omega }_{\mathrm{r}}}({\omega }_{\mathrm{r}},\beta ,U)$$

$${\bar{k}}_{\mathrm{U}}({\omega }_{\mathrm{r}},\beta ,U)$$

$${\hat{\tau }}_{\mathrm{g}}({\omega }_{\mathrm{r}})$$

Full text

††articletype: Research Article

1 Introduction

The tower makes up a substantial part of the total turbine capital costs, and therefore finding an optimum between its mass and manufacturing expenses is a critical trade-off.1 For conventional towers, diameters are limited because of land-based transportation constraints. This aspect dictates the increase of wall thickness for the production of taller towers, and consequently leads to increased weight and costs. Conventional tower designs are soft-stiff to locate the tower fundamental frequency outside the turbine variable-speed operational range, and thereby eliminate the possibility of exciting a tower resonance by the rotor rotational or blade-passing frequency. However, with the ever increasing wind turbine power rates, a combination of technical solutions should enable future, low-cost, tall towers, by relaxing this frequency constraint. Soft-soft tower configurations form an opportunity for tall towers by their smaller tower diameters and reduced wall thickness. As a result, soft-soft towers are less stiff, and have their natural frequency in the turbine operational range. Therefore more advanced control solutions will be key in avoiding the excitation of these frequencies for extended periods of time.

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.2 Moreover, vibrations are also induced by rotor aerodynamic imbalances caused by pitch errors and damage to the blade surface.3 Consequently, during variable-speed below-rated operation, the rotor rotational or blade-passing frequency may excite the structural side-side natural frequency. Small perturbations can lead to load fluctuations comparable to fore-aft stresses, because the turbine rotor provides negligible aerodynamic side-side damping, at an order of magnitude smaller than the fore-aft damping ratio.4 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 and 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 pitch5, 6 or generator torque7 control signal, for respective damping of fore-aft and side-side vibrations. Another, more passive method, entails the prevention of structural mode excitation by manipulating the generator torque when the rotor speed approaches the excitation frequency.8 This 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 power and load objectives. While an abundance of publications on advanced wind turbine control algorithms outline the possible benefits9, 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 structures10 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).11, 12 The most evident benefits of MPC over PID control13 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 in the literature for wind turbine load mitigations. A nonlinear MPC (NMPC) method is applied by assuming future wind speed knowledge using a light-detection and ranging (LIDAR) system.14 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 damping.15 In numerical simulations, the former outperforms the latter mentioned, as particularly around rated operating conditions, physical actuations constraints form a limiting factor. The benefits of NMPC using a future wind speed prediction are once again emphasized for similar operating conditions.16

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 loads and energy production over the prediction horizon, does not seem to have been backed up by literature in the past. For tall soft-soft tower configurations, the complexity lies in the fact that fatigue loads are minimized by preventing operation at the natural frequency, while it is essential to cross the same frequency for attaining higher rotational speeds and power productions. 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 ineligible for real-time applications.

Imposing spectral constraints might form a possible solution path, by employing the short-time Fourier transform (STFT) on the system output signal in a nonlinear MPC setting.17 A similar methodology18 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. The method involves a model demodulation transformation described for application and control in the field of tapping mode-atomic force microscopy (TM-AFM).19 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 description into a linear-parameter varying (LPV) model, scheduled on the excitation frequency. The technique shows similarities with the multiblade coordinate (MBC) transformation20, often used in individual pitch control (IPC) implementations for blade fatigue load reductions.21, 22, 23

An LPV system representation is frequently used for capturing nonlinear dynamics into a system description with a linear input-output mapping.24 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 proposed by solving subsequent quadratic programs (QPs).25

This paper subjects a tower model to the earlier introduced demodulation transformation, and augments it with a simplified wind turbine model, such that a qLPV model is obtained. The result is combined with the efficient MPC method, exploiting the beneficial properties of qLPV systems.25 The proposed qLPV-MPC framework provides a methodology for performing an optimal trade-off between produced energy and tower loads is presented, and thereby presents the following contributions:

• •

Providing the derivation results of a model demodulation transformation, for moving the magnitude and phase 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 demodulated 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 energy production 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 demodulation 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 explains why straightforward combination of wind turbine and tower models results in nonconvexity. Therefore, in Section 2.3, the nominal tower model is subject to a demodulation transformation to facilitate convexification. The effects and implications of the transformation are analyzed and clarified by an illustrative example in Section 2.4.

2.1 Modeling the tower dynamics as a second-order system

In practical scenarios, the center of mass of a wind turbine rotor is likely to be unaligned 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 ($1$P) periodic force, are modeled by a second-order mass-damper-spring system

[TABLE]

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 tower-top displacement, velocity and acceleration in the hub coordinate system, illustrated in Figure 1. A second-order system is taken to represent the tower first mode using the well known modal-decomposition model reduction technique26, 27, and allows for a convenient assessment and derivation of the demodulation transformation in the next section. Application of the transformation 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 Eq. (1) is split in a set of first-order differential equations by defining $x_{1}=\dot{x}(t)$ and $x_{2}=x(t)$, respectively representing the tower-top velocity and displacement, 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

[TABLE]

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}}}$. Using the operator $\stackrel{{\scriptstyle s}}{{=}}$ to equate the state-space system description to the transfer function with Laplace variable $s$, is a notation taken from Skogestad et. al.28 The notation means that the transfer function $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)$. An explicit definition of the transfer function is omitted in this work, since the notation is only used as a convenient way of parameterizing and referring to the state-space system.

2.2 Problem formalization

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

[TABLE]

in which $\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 to subject the combined nonlinear trigonometric function and LTI tower model by a demodulation transformation. The transformed tower model results in an LPV system description. The subsequent aggregation with a wind turbine model results in a quasi-LPV model, as the internal rotor speed state serves as the scheduling variable. The next section provides theory for derivation of the demodulation transformation, whereas later sections elaborate on an efficient MPC method exploiting the beneficial properties of qLPV model structures. The process of obtaining a demodulated version of a periodically excited LTI model, is referred to as the demodulation transformation in the remainder of this paper.

2.3 Theory on the tower model demodulation transformation with periodic excitation towards an LPV representation

In signal analysis, modulation is the process of imposing an information-bearing signal onto a second signal, often referred to as the carrier, carrier signal or carrier wave.29 Subsequently, the procedure of recovering the signal of interest from the modulated signal is called demodulation. The modulation-demodulation scheme is in communication systems mostly performed using dedicated hardware. A second method for performing the operation is the derivation of a mathematical framework for obtaining a demodulated model, which can be applied when the dominant dynamics of the nominal process are known. The latter mentioned approach is employed in this work, as it will appear useful for analysis, controller design, and forward propagation in MPC.

The aim of the demodulation transformation is to obtain a linear (but parameter varying) system description, which provides the frequency dependent dynamical behavior as a steady-state signal. The demodulated signal is used in a later section to form a convex quadratic optimization problem in an MPC setting for computing the optimal control signal over the prediction horizon. The demodulation transformation is applied to the assumed tower model, introduced in Section 2.1.

The model demodulation transformation is inspired by the work of19, and only the main results are given in this section. The derivation is performed and validated with symbolic manipulation software for algebraic expressions30, and the code is made publicly available.31 The transformation relies on the assumption that the changes in system response amplitude $a_{\mathrm{y}}(\tau)$ and phase change $\phi(\tau)$ are 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:

[TABLE]

Driving a linear system with a periodic input, results in a periodic response with the same frequency, however, with a certain phase shift and magnitude relative to that of the input, which is characterized by:

[TABLE]

in which $\phi\in\mathbb{R}$ is the phase shift, $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

[TABLE]

with $\mathrm{j}=\sqrt{-1}$ being the imaginary unit, and $i=\{1,2\}$, where $i=1$ relates to velocity and $i=2$ to displacement. 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:

[TABLE]

and taking the first time derivative gives

[TABLE]

By substitution of Eqs. (9) and (10) in the nominal state space representation of Eq. (4), the following expressions are obtained

[TABLE]

Furthermore, the following property of orthogonality is used, where

[TABLE]

if and only if $\{C\in\mathbb{C}\}=0$. Thus, Eqs. (11) and (12) are multiplied with $e^{\mathrm{j}\omega_{\mathrm{r}}t}$ as follows:

[TABLE]

Term-by-term integration of the integrals in Eq. (15) using the mathematical property in Eq. (6), gives the following result

[TABLE]

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 the following expression is obtained:

[TABLE]

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}}}$. 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)$.28 The instantaneous amplitude and phase of the dynamic system response at frequency $\omega_{\mathrm{r}}$ are given by

[TABLE]

It is also possible to write the result of the derivation using a summation of Kronecker products

[TABLE]

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 (Eq. (4)) into an LPV model structure (Eq. (19)). The amplitude $a_{\mathrm{u}}$ of the periodic input is in the demodulated 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 analysis.

2.4 Illustrating the effects of the transformation

This section adds context to the rather abstract derivation of the demodulation transformation in Section 2.3. Therefore, the first part of this section presents a frequency domain analysis of the transformation properties. Then, in an illustrative time domain simulation case, a frequency sweep is applied to the nominal and transformed models. The analysis and exemplary simulation clarify the characteristics and applicability of the transformation for the considered objective in this paper.

The nominal and transformed tower models from Eqs. (4) and (20) are parameterized by the following quantities: A modal mass of $m=1000$ kg, a modal damping coefficient of $\zeta=100$ kg s-1, and a modal spring constant of $k=500$ kg s-2. Resulting from the somewhat arbitrarily selected parameters, the first tower mode is located at $\omega_{\mathrm{n}}\approx 0.71$ rad s-1, with a clearly present resonance peak at the same frequency. Section 5 modifies the NREL 5-MW reference turbine tower to move its side-side fundamental frequency to the same location. Figure 4 shows Bode magnitude plots of the nominal plant and its demodulated counterpart in the frequency range ${\boldsymbol{\Omega}=\left\{\omega\,|\,\omega\subset\mathbb{R},~{}10^{-2}\leq\omega\leq 10^{1}~{}rad\,s\textsuperscript{-1}\right\}}$. To obtain the amplitude output $a_{\mathrm{y}}$ of the demodulated model, the Euclidean norm of the frequency responses of $q_{3}$ and $q_{4}$ at each frequency point is taken.

Figure 4 showcases the frequency domain effects of the transformation. The frequency responses are evaluated for 4 rotor speeds $\omega_{\mathrm{r},i}\in\boldsymbol{\Omega_{\mathrm{r}}}$, $i=\left\{1,\,2,\,3,\,4\right\}$, defined by the set $\boldsymbol{\Omega_{\mathrm{r}}}=\left\{0,\,0.5,\,0.7,\,2.0\right\}~{}\text{rad\,s\textsuperscript{-1}}$. The rotor speed elements parameterize the transformed model $H(s,\omega_{\mathrm{r}})$. The plots in Figure 4 show arrows, indicating that frequency dependent magnitude information (left) is transferred to a steady-state contribution (right). Note that for $\omega_{\mathrm{r},1}=0$ rad s-1 the transformed model reduces to the nominal case. Moreover, the right plot shows that the nominal resonance peak at $\omega_{\mathrm{n}}$ is for each frequency response split into two peaks with a $3$ dB magnitude reduction. In effect, when the input amplitude $a_{\mathrm{u}}$ of the transformed model is constant or varied slowly, the magnitude at specific nominal model frequency points is mapped to a DC contribution in the transformed case; rapid variations will result in contributions from the resonances at higher frequencies. However, in this paper, additional measures to reduce these effects, such as low-pass or notch filters, are dispensable, because $a_{\mathrm{u}}$ is assumed to be constant.

Figure 5 shows the time-domain characteristics of the transformation. For this, a frequency sweep is applied to the nominal and transformed models. For the total simulation time of $1200$ seconds, the signal has a linearly increasing frequency, with a constant increase rate of ${\dot{\omega}_{\mathrm{r}}=10^{-3}}$ rad s-2, 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 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 (linear) NREL 5-MW model, for augmentation to the demodulated 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 turbine.32 Finally, Section 3.4 validates the first-order 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 ($1$P) rotor speed, it is a logical step to augment a wind turbine model adding this state to the overall system description. A 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

[TABLE]

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

[TABLE]

in which $\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 strategy8, and $\Delta\tau_{\mathrm{g}}\in\mathbb{R}$ is an additional torque contribution resulting from the MPC framework described later in this paper. The K-omega-squared torque control law is taken as an integral part of the model, and is defined as

[TABLE]

in which $K\in\mathbb{R}^{+}$ is the optimal mode gain

[TABLE]

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 demodulated tower model $H(s,\omega_{\mathrm{r}})$, such that the following system is obtained:

[TABLE]

The above-given system description contains the nonlinear aerodynamic and generator torque input defined previously by Eqs. (24) and (25). 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 output vectors are defined as

[TABLE]

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 Eqs. (27) and (28) with respect to the state and inputs vectors, such that a linear state-space system is obtained

[TABLE]

in which the state, input and output matrices are defined as

[TABLE]

The aerodynamic rotor torque is linearized with respect to the rotor speed and wind speed

[TABLE]

with

[TABLE]

and $c_{\mathrm{r}}=0.5\rho\pi R^{3}$ is a constant factor. Finally, the K-omega-squared torque controller is linearized as

[TABLE]

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 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 self-scheduling at 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 Eq. (32) 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 incorporating these operating point offsets, converting the LPV model to its affine form, 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

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 nonlinear 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 operational rotor and wind speed. The torque coefficient data is obtained using a graphical extension33 to NREL’s high-fidelity wind turbine simulation software FAST v8.1634, which includes blade element momentum (BEM) code4 for obtaining rotor characteristic data. As the framework being derived in this paper focuses on the below-rated region, and conventional wind turbine controllers keep the pitch angle fixed at fine-pitch angle $\beta_{0}$ during partial load35, the dependency of the torque coefficient on $\beta$ is omitted. An often used parameterizable torque coefficient function is defined by

[TABLE]

which is fitted by optimizing the values $\theta_{i}$ using a nonlinear 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 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 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:

A nonlinear 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 the linear wind turbine rotor and transformed 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-top 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 steady-state gains $\bar{q}_{1-4}$ being more sensitive 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.

4 Quasi-LPV model predictive control

Nonlinear MPC is – because of its computational complexity – often considered as an unsuitable control method for application in fast real-time systems. Therefore, in this paper, an approach towards an efficient method for nonlinear 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 maximizing power production efficiency, and minimizing tower natural frequency excitation. Practically, this means that the rotor deviates from the 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.36 For the considered case, a predefined quadratic performance criterion specifies the trade-off between energy production maximization and resonance excitation minimization, and the optimizer finds the optimal corresponding control signal in the prediction horizon. 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.25 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 on the consecutive output differences is used to determine whether the algorithm has converged.

By manipulation of the affine system representation defined by Eqs. (47) and (48), an expression is derived for forward propagation of the qLPV model output, only requiring the initial state at time instant $k$ and the scheduling sequence over the prediction horizon:

[TABLE]

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 to equal $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 instant $k=0$, only the initial state is assumed to be known, and the scheduling parameters are chosen constant over the prediction horizon, such that:

[TABLE]

in which ${\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

[TABLE]

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-top displacement 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 Eq. (5). The first term of Eq. (41) 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 power production and load reductions by varying the weight ratio of $Q$ and $R$.

After the first solve with the initial scheduling sequence of Eq. (40), 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 predefined error threshold. The algorithm is summarized using pseudocode in Algorithm 1.

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 time instants $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 in each step, makes the approach for solving the nonlinear MPC problem computationally efficient and tractable for real-world implementations.

5 High-fidelity simulation setup and results

This section implements the proposed qLPV-MPC framework in conjunction with the NREL 5-MW reference wind turbine model in the high-fidelity FAST code. The software implementation is made publicly available.31 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 wall thickness is scaled down by a factor $7.5$ to mimic the characteristics of a tall, more flexible, soft-soft tower configuration. As a result of the reduced thickness, an effective turbine side-side resonance frequency of approximately $0.71$ rad s-1 is attained, equal to $\omega_{\mathrm{n}}$ defined in Section 2.4. 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 demodulated second-order tower model from Eq. (19). 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 Eq. (39).

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.

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 estimator is used37, 38, which is also plotted in the same figure. Because the future wind speed is unknown at time instant $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 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 after $100$ s, a discrepancy is seen at this time instant. 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 reference turbine 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 demodulation 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 Simulink39, simulates for $t_{\mathrm{s,MPC}}$, after which is simulation is paused, and essential information is extracted. The simulation data is provided to the MPC algorithm in MATLAB using CVX: A package for specifying and solving convex programs.40, 41 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-top 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 produced energy.

The simulation results for Case 2 are given in Figure 11. By inspection of the rotor speed around the resonance frequency, it shows that the qLPV-MPC implementation prevents operation at this speed for extended time periods. This operational strategy results in a significant decrease of tower-top displacement amplitudes. To further clarify this effect, Figure 12 shows histograms of the rotor speed and displacement amplitude signals. Furthermore, Figures 13a and 13b show the sidewards displacement spectra of the two control strategies. A significant reduction of $18$ dB is attained at the turbine side-side natural frequency. Finally, a fatigue assessment is performed on the tower base side-side moment by evaluating the damage equivalent loads (DEL) from the corresponding time-domain signals using MLife.42 The DEL measure quantifies the amplitude of a certain harmonic load variation that would cause the same damage level when repeated for a given amount of cycles.43, 44 In this fatigue analysis, a Wöhler-exponent of $4$ is chosen, which is a typical value for steel.3 Figure 13c shows the normalized DEL values with respect to the baseline case: A significant $52$ % DEL reduction is attained. Since the analysis is based on the single load case performed in this section, the short-term DELs are calculated using the $1$ Hz equivalent load.4 This implies that the number of cycles is equal to the simulation time, and consequently makes the analysis independent of the simulation runtime.

The generator power and torque trajectories of the control strategies show a high degree of similarity, which indicates a minimal penalty on the overall energy production. The observation is confirmed by the evaluation of the produced energy over the total simulation time, resulting in $603.34$ kWh and $601.19$ kWh for the respective baseline and qLPV-MPC cases, which turns out in a negligible produced energy reduction of $0.36$ %. The trade-off is conveniently tuned by varying 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 side-side excitation. For future turbines with higher power ratings and taller towers to be weight and cost effective, soft-soft tower configurations are being considered. Such towers are more flexible, and have their side-side natural frequency in the variable-speed domain, possibly coinciding with the rotor rotational or blade passing frequencies.

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 demodulation 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 profiles show that the algorithm prevents persistent excitation of the tower fundamental frequency, 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.31

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

[TABLE]

are substituted in Eq. (32), such that the affine form is obtained:

[TABLE]

in which $\boldsymbol{p}(t)=\boldsymbol{p}^{*}(t)$ indicates the current linear model in the LPV scheduling space.45 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 setup, the continuous-time system is converted to its discrete-time equivalent

[TABLE]

in which $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 method46, of which the matrix transformation relations are given by

[TABLE]

Note that the last term of Eq. (47), 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 Eq. (39):

[TABLE]