On the Potential of Reduced Order Models for Wind Farm Control: A Koopman Dynamic Mode Decomposition Approach

Abstract

The high dimensions and governing non-linear dynamics in wind farm systems make the design of numerical optimal controllers computationally expensive. A possible pathway to circumvent this challenge lies in finding reduced order models which can accurately embed the existing non-linearities. The work presented here applies the ideas motivated by non-linear dynamical systems theory—the Koopman Operator—to an innovative algorithm in the context of wind farm systems—Input Output Dynamic Mode Decomposition (IODMD)—to improve on the ability to model the aerodynamic interaction between wind turbines in a wind farm and uncover insights into the existing dynamics. It is shown that a reduced order linear state space model can reproduce the downstream turbine generator power dynamics and reconstruct the upstream turbine wake. It is further shown that the fit can be improved by judiciously choosing the Koopman observables used in the IODMD algorithm without jeopardizing the models ability to rebuild the turbine wake. The extensions to the IODMD algorithm provide an important step towards the design of linear reduced order models which can accurately reproduce the dynamics in a wind farm.

1. Introduction

Renewable energy production has increased substantially over recent years, with a 5 to 6 fold increase since the 1960s [1]. Wind energy is one of the fast growing technologies, representing 19% of the total renewable energy production in 2018, compared to 8% in 2010 [1]. The total installed wind power capacity in Europe has been steadily increasing since 2009 [2]. Statistics indicate an increase from 77 GW (75 onshore and 2 offshore) of cumulative wind power capacity in 2009 to 205 GW (183 onshore and 22 offshore) in 2019.

The trend for going offshore is explained by its advantages, namely the high availability of resources and potential for reduced costs. Nevertheless, as wind turbines are sited together in wind farms, the aerodynamic interaction between turbines may have a negative effect on the total electrical power production. These effects are caused by the wakes, which are flow structures that form behind each turbine and are characterized by reduced flow velocity and increased turbulence intensity. Wind plant control techniques emerge as a possible solution to mitigate the aerodynamic negative effects by maximizing the power output of the plant instead of each individual turbine. The concept of Wind Farm Control (WFC) was firstly introduced in 1988 in [3], showing that slight increases in the overall energy capture are possible. Since then, different control techniques have emerged, aiming to increase power production and decrease turbine loads.

Wind Farm Control can be categorized into three distinct technologies [4], namely wake steering, Axial Induction Control (AIC) and Dynamic Induction Control (DIC). The first focuses on redirecting the wake interaction, the second on modifying interactions and impacts by derating upstream or uprating downstream turbines and the third on inducing additional wake recovery by dynamically uprating and downrating upstream turbines on short time scales.

The focus of this work is on AIC and DIC. Despite initial promising results, recent studies with high-fidelity simulation models and scaled wind tunnel experiments indicate that the power gain that can be achieved with steady state induction control is limited to non existing. In [5], it is shown that the control settings that yield gain at the wind plant level are dependent on particular atmospheric conditions, the wind plant configuration and turbine characteristics. Simulations in SOWFA, a simulator of high fidelity, further suggest the existence of certain circumstances where this strategy becomes infeasible and that the power lost on the first turbine is not regained at the second turbine. Wind tunnel tests point again in the same direction: in [6], only very limited power gains are measured, showing a modest increase in total power of circa 1% when two turbines are considered.

Alternative approaches that focus on time-varying control inputs which purposely influence the inherently dynamic nature of the wake have emerged. Recent studies [7,8], have shown that DIC, where the induction factor is varied over time, improves power production in small to medium-sized wind farms by generating a turbulent wind flow which enables enhanced wake recovery. In the former, power gains for wind farms with optimal coordinated control in the order of 8–21% compared with a greedy controlled case are found. In the latter, these results are explored in more depth, where the simpler approach of varying the induction sinusoidally is taken. A best parameter for the sinusoidal variations is found, corresponding to an amplitude $A^{*}$ of $1.5$ and a non-dimensional frequency $S{t}^{*}$ of 0.25, proving to be robust to both wind-turbine spacing and turbulence intensity. Wind tunnels experiments with DIC detailed in [9] prove the potential of such control strategy in real-world experiments, showing a potential increase in power production up to 4% in a three-scaled wind turbine setup.

Finding the optimal control settings in wind turbines that account for the aerodynamic interaction at the wind farm level in real time has thus become a challenging and active topic of research. However, developing optimal control frameworks becomes computationally costly and unsuitable for real time control. Reduced Order Models (ROMs) provide a possible pathway to design and implement in real time optimal controllers by using model order reduction techniques.

Proper Orthogonal Decomposition (POD) and Dynamic Mode Decomposition (DMD) are examples of such approaches. In [10], POD is used to decompose and analyze the growth and development of structures in the wake. The POD is then applied to sets of modes of common rank from each separate POD basis to show the existing sub-modal structures. This technique, termed Double Proper Orthogonal Decomposition (DPOD), offers a meaning of characterizing the spatial developments of POD modes, presenting insights beyond the snapshot POD in addition to filtering noise and errors deriving from the initial eigenvalue decompositions. In [11], a data-driven ROM is derived through a polynomial expansion which relates the POD coefficients to their respective time derivatives through a least-squares minimization. Polynomial coefficients up to cubic order are then used to couple the dynamic evolution of the coefficients with Ordinary Differential Equations (ODE). Fluctuating velocity fields in the wake are then predicted by solving the coupled set of ODEs.

In [12], both POD and DMD are used to compute a ROM to detect and predict relevant wake dynamics. Wake vorticity structures are captured and compared using data from two simulations, one using a classical actuator line model and the other including turbine tower and nacelle. It is found that both modal decomposition techniques produce similar results, with the advantage that DMD modes of interest can be selected to formulate a time-marching algorithm that can be used to predict wake behaviour. In [13], a DMD model is embedded in a Kalman filter and it is shown that the entire wake flow behind a turbine can be predicted with a global accuracy of 4%. This application is only for a single turbine whilst a two turbine setting arises in [14], where the findings indicate that a ROM can be obtained and used as a predictive model where a forced input is used.

The purpose of this work is to explore the ideas of non-linear dynamical system theory in the context of data-driven dynamical system to improve the capabilities of reduced order models in reproducing the nonlinear dynamics inherent in WFC strategies, more specifically in induction control. This task represents a pivotal step in the design of accurate predictive controllers which can be implemented with low computational cost.

Two data sets containing information of the power production and the wake characteristics of a two turbine wind farm where collective pitch control is tested are used. One of the data sets is used to compute the ROMs (identification data set) and the other is used to validate the ROMs (validation data set). Several variables that characterize the wake behaviour are used to compute the ROMs, such as different velocity components and non-linear functions of the latter. The ROMs are evaluated by quantifying the accuracy in predicting the power of both turbines and the wake velocity.

This article is divided into four sections: Section 2 introduces the mathematical reduced order modelling techniques used for control as well as the simulation setup used to gather data. Section 3 presents the ROMs derived using different techniques and in Section 4 results and future research directions are discussed.

2. Methods

2.1. Input Output Dynamic Mode Decomposition

Dynamic Mode Decomposition (DMD) is a data-driven technique that approximates the nonlinear dynamics of a system in an equation-free way. It originated in the fluid dynamics community as a means to decompose flows into a simple representation [15]. DMD extracts meaningful spatio-temporal patterns which dominate dynamic activity in high dimensional non-linear systems and it has been applied to a variety of flow geometries and proven its value as a new quantitative flow analysis tool when compared to other methods, such as POD [16]. The standard DMD algorithm has been applied in different fields, such as neuroscience [17] and epidemiology [18], providing a deeper understanding of the systems being studied. Furthermore, different variations have surged since its first publication [19], generalizing DMD to a larger class of datasets, including non sequential time series.

In the context of control, DMD has been developed to account for exogenous inputs [20] and input–output information [14,21]. The latter approach, referred to as Input Output Dynamic Mode Decomposition (IODMD), combines DMD with standard subspace identification often used in the controls literature to obtain a characterization of the system. The discrete time state space representation of a Linear Time Invariant (LTI) system is used:

$$\begin{array}{c}\hfill {\mathit{x}}_{k+1}=\mathbf{A}{\mathit{x}}_k+\mathbf{B}{\mathit{u}}_k\\ \hfill {\mathit{y}}_k=\mathbf{C}{\mathit{x}}_k+\mathbf{D}{\mathit{u}}_k\end{array}$$

(1)

where ${\mathit{x}}_k\in {\mathbb{R}}^n$, ${\mathit{u}}_k\in {\mathbb{R}}^q$ and ${\mathit{y}}_k\in {\mathbb{R}}^l$ are the state, input and output vectors respectively. Matrices $\mathbf{A}\in {\mathbb{R}}^{n\times n}$, $\mathbf{B}\in {\mathbb{R}}^{n\times q}$, $\mathbf{C}\in {\mathbb{R}}^{l\times n}$ and $\mathbf{D}\in {\mathbb{R}}^{l\times q}$ are the state matrices that together form the state space representation. Data regarding state, input and output information is collected and used to reverse engineer the governing equations of the system. For that purpose, the IODMD algorithm can be divided into three distinct steps:

(1)

Organize snapshot matrices: the series of data containing relevant information of the three dimensional vector flow field at each time instant k is taken and flatten into a single column vector ${\mathit{x}}_k$, referred to as a snapshot of data. All time instants are combined into one single matrix. The latter is divided into the data matrix X and its time shifted version X’:

$$\mathbf{X}=\left[\begin{array}{cccc}\mid & \mid & & \mid \\ {\mathit{x}}_1& {\mathit{x}}_2& \cdots & {\mathit{x}}_{m-1}\\ \mid & \mid & & \mid \end{array}\right]{\mathbf{X}}^{\prime }=\left[\begin{array}{cccc}\mid & \mid & & \mid \\ {\mathit{x}}_2& {\mathit{x}}_3& \cdots & {\mathit{x}}_m\\ \mid & \mid & & \mid \end{array}\right]$$

(2)

where the number of rows n is equal to the number of spatial points saved per time snapshot and the number of columns m is equal to the number of snapshots taken. The input $\mathsf{\Psi }$ and output $\mathbf{Y}$ snapshot matrices are similarly gathered:

$$\mathsf{\Psi }=\left[\begin{array}{cccc}\mid & \mid & & \mid \\ {\mathit{u}}_1& {\mathit{u}}_2& \cdots & {\mathit{u}}_{m-1}\\ \mid & \mid & & \mid \end{array}\right]\mathbf{Y}=\left[\begin{array}{cccc}\mid & \mid & & \mid \\ {\mathit{y}}_1& {\mathit{y}}_2& \cdots & {\mathit{y}}_{m-1}\\ \mid & \mid & & \mid \end{array}\right]$$

(3)

(2)

Take the singular value decomposition of X: by appropriately choosing a truncation value r of the number of singular values to retain, X can be factorized as $\mathbf{U}\Sigma {\mathbf{V}}^{*}$, where $\mathbf{U}\in {\mathbb{R}}^{n\times r}$, $\Sigma \in {\mathbb{R}}^{r\times r}$, ${\mathbf{V}}^{*}\in {\mathbb{R}}^{r\times m-1}$ and ${}^{*}$ is the conjugate transpose.

(3)

Formulate least squares problem by projecting system to reduced dimension space: Equation (1) is then rewritten in terms of the data matrices in (2) and the input and output snapshot matrices in (3) in the compact format:

$$\left[\begin{array}{c}{\mathbf{X}}^{\prime }\\ \mathbf{Y}\end{array}\right]=\left[\begin{array}{cc}\mathbf{A}& \mathbf{B}\\ \mathbf{C}& \mathbf{D}\end{array}\right]\left[\begin{array}{c}\mathbf{X}\\ \mathsf{\Psi }\end{array}\right]$$

(4)

The state is then projected onto the subspace defined by an orthonormal basis given by the POD modes of X. The columns of U specify the orthonormal basis to project the high order state. The low order state matrices are then $\tilde{\mathbf{A}}={\mathbf{U}}^{*}\mathbf{A}\mathbf{U}$, $\tilde{\mathbf{B}}={\mathbf{U}}^{*}\mathbf{B}$, $\tilde{\mathbf{C}}=\mathbf{C}\mathbf{U}$ and $\tilde{\mathbf{D}}=\mathbf{D}$. Substituting in (4) and rearranging, the reduced order state matrices are obtained by minimizing the error of the Frobenius norm [22]:

$$\underset{\left[\begin{array}{cc}\tilde{\mathbf{A}}& \tilde{\mathbf{B}}\\ \tilde{\mathbf{C}}& \tilde{\mathbf{D}}\end{array}\right]}{\min }{∥\left[\begin{array}{c}{\mathbf{X}}^{\prime }\\ \mathbf{Y}\end{array}\right]-\left[\begin{array}{cc}\mathbf{U}& \mathbf{0}\\ \mathbf{0}& \mathbf{I}\end{array}\right]\left[\begin{array}{cc}\tilde{\mathbf{A}}& \tilde{\mathbf{B}}\\ \tilde{\mathbf{C}}& \tilde{\mathbf{D}}\end{array}\right]\left[\begin{array}{cc}{\mathbf{U}}^{*}& \mathbf{0}\\ \mathbf{0}& \mathbf{I}\end{array}\right]\left[\begin{array}{c}\mathbf{X}\\ \mathsf{\Psi }\end{array}\right]∥}_F^2$$

(5)

where $\tilde{\mathbf{A}}\in {\mathbb{R}}^{r\times r}$, $\tilde{\mathbf{B}}\in {\mathbb{R}}^{r\times q}$, $\tilde{\mathbf{C}}\in {\mathbb{R}}^{l\times r}$, and $\tilde{\mathbf{D}}\in {\mathbb{R}}^{l\times q}$ form the state space representation of the ROM. The reduced order state matrices are then obtained by least squares, taking into account that the reduced order representation of X and X’ are, respectively, ${\mathbf{U}}^{*}{\mathbf{X}}^{\prime }$ and $\Sigma \mathbf{V}$:

$${\Theta }_{\mathrm{IODMD}}=\left[\begin{array}{cc}\tilde{\mathbf{A}}& \tilde{\mathbf{B}}\\ \tilde{\mathbf{C}}& \tilde{\mathbf{D}}\end{array}\right]=\left[\begin{array}{c}{\mathbf{U}}^{*}{\mathbf{X}}^{\prime }\\ \mathbf{Y}\end{array}\right]{\left[\begin{array}{c}\Sigma \mathbf{V}\\ \mathsf{\Psi }\end{array}\right]}^{†}$$

(6)

where † is the pseudoinverse. Using this methodology, an accurate mapping for input–output information can be obtained, as well as an approximation of the full order state, by ${\mathit{x}}_k=\mathbf{U}{\tilde{\mathit{x}}}_k$.

In the standard DMD algorithm, described in [23], an eigendecomposition is performed to the reduced order state space matrix $\tilde{\mathbf{A}}$, $\tilde{\mathbf{A}}\mathbf{W}=\mathbf{W}\mathsf{\Lambda }$, where columns of W are eigenvectors and $\mathsf{\Lambda }$ is a diagonal matrix containing the corresponding eigenvalues ${\mathbf{\lambda }}_k$. These can be investigated for fundamental properties of the underlying system such as growth modes and resonance frequencies. To approximate the eigenvalues and eigenvectors of A without its explicit computation, an approximation is calculated based on the so called DMD modes, which correspond to the columns of $\mathsf{\Phi }$, $\mathsf{\Phi }$ = UW, as originally proposed in [24]. These representations of the high order eigenvectors based on the low order eigenvectors describe the dynamics observed in the time series in terms of oscillatory components.

2.2. Koopman Operator Theory

The Koopman operator, first introduced in 1931 [25] and recently making a resurgence in the context of modern data-driven dynamical systems, is defined as an infinite-dimensional linear operator that evolves fields of physical observables defined on a state space of a dynamical system. More explicitly, let the generic dynamical system in Equation (7) be considered:

$$\frac{d\mathit{x}}{dt}=\mathbf{f}(\mathit{x};t;\mu )$$

(7)

where $\mathit{x}\left(t\right)\in {\mathbb{R}}^n$ is a vector representing the state of the dynamical system at time t, $\mu$ contains parameters of the system and f (·) represents the dynamics, defined on a state space M. This induces a discrete-time dynamical system given by the flow map ${\mathbf{F}}_t:\mathbf{M}\to \mathbf{M}$, ${\mathit{x}}_{k+1}={\mathbf{F}}_t\left({\mathit{x}}_k\right)$. A function defined as $g:\mathbf{M}\to \mathbb{R}$ is called an observable of the system. The Koopman operator ${\mathcal{K}}_t$ acts on observable functions g as ${\mathcal{K}}_t=g\circ {\mathbf{F}}_t$ where ∘ is the composition operator, so that:

$${\mathcal{K}}_{t}g\left({\mathit{x}}_k\right)=g\left({\mathbf{F}}_t\left({\mathit{x}}_k\right)\right)=g\left({\mathit{x}}_{k+1}\right)$$

(8)

$$g\left({\mathit{x}}_{k+1}\right)={\mathcal{K}}_{t}g\left({\mathit{x}}_k\right)$$

(9)

It is now evident in (9) that the observation of the state, or, in other words, a function of the measured state variable, is advanced to the next time step by the Koopman operator. This is true for any observable function and for any point ${\mathit{x}}_k\in \mathbf{M}$, hence the infinite-dimension. The spectral decomposition of the Koopman operator is instrumental in representing solutions to a dynamical system of interest. Considering the eigenvalue problem $\mathcal{K}{\phi }_k={\lambda }_k{\phi }_k$, the functions ${\phi }_k\left(\mathit{x}\right)$ are called Koopman eigenfunctions and they define a set of intrinsic measurement coordinates on which it is possible to advance these measurements with a linear dynamical system. The new measurement space formed by a vector of observables $\mathit{g}\left(\mathit{x}\right)$ can be expressed in terms of the Koopman eigenfunctions:

$$\mathit{g}\left(\mathit{x}\right)=\left[\begin{array}{c}{g}_1\left(\mathit{x}\right)\\ ⋮\\ g_p\left(\mathit{x}\right)\end{array}\right]=\sum _{k=1}^{\infty }{\phi }_k\left(\mathit{x}\right){\mathit{v}}_k$$

(10)

where ${\mathit{v}}_k$ is the kth Koopman mode of the map F corresponding to the observable $\mathit{g}$, associated with the kth Koopman eigenfunction, i.e., the weighting of each observable on the eigenfunction. The expression in (10) can be interpreted as a linear combination of the eigenfunctions ${\phi }_k$ of $\mathcal{K}$ where ${\mathit{v}}_k$ is the coefficient in the expansion. The nonlinear system can either be evolved in the original state space or in the measurement space as in (10).

$$\mathcal{K}\sum _{k=1}^{\infty }{\phi }_k\left(\mathit{x}\right){\mathit{v}}_k=\sum _{k=1}^{\infty }\mathcal{K}{\phi }_k\left(\mathit{x}\right){\mathit{v}}_k=\sum _{k=1}^{\infty }{\lambda }_k{\phi }_k\left(\mathit{x}\right){\mathit{v}}_k$$

(11)

Future solutions can be computed by multiplication with the Koopman eigenvalues, ${\lambda }_k\in \mathbb{C}$, which characterize the temporal behaviour of the corresponding Koopman mode ${\mathit{v}}_k$.

The connections of DMD to the underlying non-linear dynamics of a system through the Koopman operator theory were firstly established in [26], the that DMD produces approximations to eigenvalues of the Koopman operator and their corresponding modes. DMD is thus a manifestation of Koopman theory for the specific case when the observable functions are the identity or a linear transformation of the underlying state space. DMD provides a tractable computational alternative, as the finite-dimensional matrices it yields are an approximation to the infinite dimensional operator.

Approaches combining DMD with the ideas of the Koopman operator by augmenting DMD to include a richer set of non-linear observables and find the best set of coordinates to describe the non-linear dynamics have since been developed [27,28]. The latter approach is designated by Extended Dynamic Mode Decomposition (EDMD). The major challenges it faces concern the choosing of the non-linear observables, as shown in [29,30], and the loss of interpretability. An alternative approach is to select a parsimonious set of observables based upon knowledge of the physics of the system, as described in [30] and proved to yield better results when compared to POD for prediction purposes [31]. The linear representation provided by the Koopman Operator is not only useful for prediction but also for control, as explored in the works of [32,33], where a structured approach is given to derive linear models suited for predictive control by doing a nonlinear transformation of the data, termed the lifting-and a linear least squares problem in the lifted space, based on the ideas of the Koopman Operator.

In this work, the ideas of the Koopman Operator are explored in the context of WFC. Three modifications are devised to the standard IODMD algorithm:

(1)

IODMD using different observables: the IODMD algorithm is tested using different variables of the velocity flow field or combination of variables. The common variable used as a state, the streamwise velocity component, u, is replaced by others, such as the spanwise velocity component v, the vertical velocity component w or a non-linear combination of the three. This procedure attempts to find the coordinates best suited to describe the dynamics of AIC. The results of this formulation are referred to as IODMD${}_{xx}$, where $xx$ is the observable used. In terms of the IODMD described in Section 2.1, only the snapshot matrices organisation step differs, with the data matrices being generalized.

$$\mathbf{X}=\left[\begin{array}{cccc}\mid & \mid & & \mid \\ g\left({\mathit{x}}_1\right)& g\left({\mathit{x}}_2\right)& \cdots & g\left({\mathit{x}}_m\right)\\ \mid & \mid & & \mid \end{array}\right]{\mathbf{X}}^{\prime }=\left[\begin{array}{cccc}\mid & \mid & & \mid \\ g\left({\mathit{x}}_1^{\prime }\right)& g\left({\mathit{x}}_2^{\prime }\right)& \cdots & g\left(x_m^{\prime }\right)\\ \mid & \mid & & \mid \end{array}\right]$$

(12)

where g is a function of the velocity vector variable. This approach has been applied to Wake Redirection Control, another WFC strategy, in [34].

(2)

Augment state vector by another flow variable: the IODMD algorithm is tested by augmenting the data matrices with different observables simultaneously, as in [30]. This formulation is termed Extended Input Output Dynamic Mode Decomposition (EIODMD${}_{xx,yy}$) where $xx$ is the first observable and $yy$ the second observable used to augment the original snapshot matrix. In this approach, the data matrices are generalized to:

$$\mathbf{X}=\left[\begin{array}{cccc}\mid & \mid & & \mid \\ \mathit{g}\left({\mathit{x}}_1\right)& \mathit{g}\left({\mathit{x}}_2\right)& \cdots & \mathit{g}\left({\mathit{x}}_m\right)\\ \mid & \mid & & \mid \end{array}\right]{\mathbf{X}}^{\prime }=\left[\begin{array}{cccc}\mid & \mid & & \mid \\ \mathit{g}\left({\mathit{x}}_1^{\prime }\right)& \mathit{g}\left({\mathit{x}}_2^{\prime }\right)& \cdots & \mathit{g}\left(x_m^{\prime }\right)\\ \mid & \mid & & \mid \end{array}\right]$$

(13)

where $\mathit{g}={\left[\begin{array}{cccc}{g}_1\left(x\right)& g_2\left(x\right)& \cdots & g_p\left(x\right)\end{array}\right]}^T$ is a vector of functions of velocity field variables.

(3)

Augment the state vector by any variable: the third hypothesis focuses on augmenting the flow field information with turbine information, following the ideas laid out in [33]. The additional states are referred to as deterministic states ${\mathbf{X}}_d$ and the corresponding data matrices are organized as in (2). The least squares problem to solve boils down to:

$${\Theta }_{\mathrm{EIODMD}}=\left[\begin{array}{c}{{\mathbf{X}}^{\prime }}_d\\ {\mathbf{U}}^{*}{\mathbf{X}}^{\prime }\\ \mathbf{Y}\end{array}\right]{\left[\begin{array}{c}{\mathbf{X}}_d\\ \Sigma \mathbf{V}\\ \mathsf{\Psi }\end{array}\right]}^{†}$$

(14)

2.3. Simulation

In light of the pure data-driven approach employed, a simulation is performed to obtain measurements to be used for IODMD. The goal is to obtain a ROM which can accurately describe the generator power and wake dynamics of two wind turbines interacting in an offshore wind farm.

The input–output data set to be used for IODMD is obtained by performing a Computational Fluid Dynamics (CFD) simulation in Simulator for Offshore Wind Farm Applications (SOWFA) [35]. SOWFA implements a Large Eddy Simulation (LES) of the three-dimensional wind flow around one or more turbine rotors in the Atmospheric Boundary Layer (ABL) developed by the United States National Renewable Energy Laboratory (NREL). In the LES method the larger scales of the flow field are resolved by solving the temporally and spatially unsteady Navier–Stokes equations in order to limit the computational cost of the simulation. The rotating rotor blades are modelled through an actuator line approach. The actuator lines are coupled with the FAST turbine aeroelastic simulator tool which calculates the loads, power and rotor speed of each turbine. SOWFA is a widely used tool in the literature. More information is available in [35].

In this work, wind with a uniform inflow profile is used to demonstrate the working principles of IODMD. The case of a wind farm with two turbines is studied, to assess the gain in energy in the second one. The downstream turbine is situated five rotor diameters (5D) behind the upstream turbine. The full domain as well as inflow characteristics used for the CFD simulation are the same as the ones used in [36], where the complete wind farm domain size corresponds to 2.5 km × 1 km × 0.6 km and the uniform inflow wind speed is 9 m/s. The simulation is of a neutral ABL in which the inflow was generated through a so-called precursor simulation. The surface roughness and turbulence intensities are set to zero for the purposes of this simulation. The static uniform inflow profile and absence of turbulence are not representative of true wind conditions present in wind farms. Nevertheless, the results gathered using this simplified scenario provide a starting point for more sophisticated configurations.

The turbine used for the simulation is the DTU 10MW turbine [37], with a rotor diameter of 178.3 m. The domain size used for the purposes of IODMD has dimensions of approximate 1.71 km in the x direction, 343 m in the y direction and 300 m in the z direction. A total of 550 points in x, 111 points in y and 99 points in z are sampled. A total of 6 million points, each one containing information on the $(u,v,w)$ velocity components is taken.

A pre processing step of the data is performed in a cluster, aiming to retrieve only information at each fourth point, decreasing the total number of points from 6 million to roughly 100 thousand samples, decreasing computational effort. Flow field information is sampled at each 2 s. In addition to flow field data, information of the blades, nacelle, rotor, tower and generator power is available. Randomized Binary Signals (RBS) provide a setpoint for the collective blades pitch angle of the first turbine ${\theta }_{t1}$. These serve as inputs and are changed step-wise from 0 to 4 degrees. The blades pitch angle of the downstream turbine remains the same throughout the simulation.

Three pre processing steps are done before using simulation data to perform systems identification: (1) only flow field information between the first and second turbine is used (2) turbine information is resampled to be consistent with flow field data (3) trends and offsets are manually removed from data, so that only variations of power with relation to the operating point where the blades are not pitched are evaluated. Consequently, the generator power of both turbines corresponding to this point is subtracted from the data, which coincides to a value of 5.485 MW for the upstream turbine and of 0.773 MW for the downstream turbine.

The input–output information, used for identification and validation of the models where a collective pitch strategy is used, is depicted in Figure 1. In the latter, the influence of the collective blade pitch step-wise changes on the power produced in each turbine is visible. Taking the identification data set as an example, it is seen that a variation from 0 to 4 degrees, near the beginning of the simulation, leads to an almost immediate decrease of approximately 1 MW in power produced by the upstream turbine. As a consequence of less energy being extracted from the wind, the power production of the downstream turbine increases. This effect is only felt 2 min later, and the power production increases 0.8 MW approximately 2 min and 24 s later. However, the power lost on the first turbine cannot be fully recovered by the downstream turbine, as the last graph shows.

Momentarily gains are observable when the blades pitch back to 0 degrees and a rapid increase in power captured in the first turbine is registered. This represent a gain of 0.9 MW in total power produced, when compared to the initial turbine operating point. As the wake suffers alterations as a consequence of the blades pitching back to zero degrees and the wake travels and hits the downstream turbine rotor, a peak in power is registered, which almost reaches 2.4 MW. As the wake continues to propagate, the power productions of both turbines return to the baseline operating levels.

3. Results

In this section, the techniques introduced in Section 2.1 and Section 2.2 are used to obtain the desired input–output mapping illustrated in Figure 1.

3.1. Input Output Dynamic Mode Decomposition with Different Observables

The streamwise velocity component u is firstly used to derive the IODMD model, referred to as IODMD${}_u$. Singular values are increasingly added to compute models. A total of r models, corresponding to the truncation level used for the SVD are computed and validated.

Figure 2 depicts the fit of each one, measured by the Variance Accounted For (VAF) of both outputs. A certain number of singular values (or modes) is required to properly recover the dynamics of generator power, which is especially true for the downstream turbine. Furthermore, adding more modes does not directly lead to more accurate models. This patterns is recognizable in Figure 2 where the gap between VAF for identification and validation widens as more modes are used, possibly explained by overfitting, meaning that the models correspond too closely to the data set used for identification and are not able to reproduce the dynamics of the new data set (validation data set) with the same accuracy. This is also supported by the observation that the accuracy of the models tends to saturate as more singular values are added, for both the identification data set (accuracy saturates around 80%) and validation data set (accuracy saturates around 65%).

The criteria used to choose the best performing model is the fitness of the outputs, quantified by the VAF. Nevertheless, the power dynamics of the downstream turbine are more difficult to model when compared to the upstream turbine, hence the final decision criteria is the fitness of the downstream turbine power (second model output). Based on the information in Figure 2, the best performing model has 38 states and presents a fit of 73.59% for the downstream turbine generator power, depicted in Figure 3.

The IODMD${}_u$ model can also be evaluated in its ability to reproduce the upstream turbine wake characteristics. The state space trajectory of the ROM is written in the high order coordinates by ${\mathit{x}}_k=\mathbf{U}{\tilde{\mathit{x}}}_k$ and the streamwise velocity component u can be entirely reconstructed and compared to SOWFA validation measurements, as represented in Figure 4 and Figure 5.

Figure 4 and Figure 5 illustrate the prediction capabilities of the ROM. It is seen in Figure 4 that the ROM can provide accurate estimates of the normalised streamwise velocity at hub height by visually comparing the reconstructed wake and SOWFA data. Nevertheless, it is also noticeable how the ROM fails to clearly define the borders of the wake and how the different velocity regions are not precisely distinguishable.

Figure 5 quantifies the deviations between the ROM predictions and the validation data set by computing the absolute relative deviation at different grid points located at the downstream turbine rotor plane. The first row shows a highly accurate prediction, with deviations in the order of 10% appearing only at certain locations. However, the blades collective pitch angle of the upstream turbine change, inducing alterations in the wake. The middle snapshot shows a critical moment when the altered wake hits the downstream turbine rotor plane, and deviations of 30% are registered at the outer grid points.

It is hypothesized that the transient regime characterized by high non-linear wake behaviour is more difficult to predict by using a simple linear reduced order model. As the wake is again fully developed, now with lower velocity deficit as a consequence of the adjusted pitch angles of the first turbine, the predictions are again more accurate, although sporadic large deviations are still found.

The deviations between IODMD${}_u$ streamwise velocity reconstruction and SOWFA measurements are quantified for each time instant by computing the Normalized Root Mean Squared Error (NRMSE). Taking all time instants into account, the deviations can be summarized resorting to an average NRMSE, (which is evaluated only for the best performing model, as measured by the fitness of the outputs). The best performing model results are summarized in

Table 1. Models validation results using IODMD${}_u$ algorithm, where u is the streamwise velocity component. First column has the data chosen for the states. Second column has the fit criteria VAF (the second turbine generator power is the most relevant, hence acts as a primary decision criteria). Third column has the number of states of the model. Fourth column the fit of the upstream turbine generator power for the chosen model. Fifth and sixth columns present the average NRMSE for identification and validation data.

Model Properties	Model Obtained by IODMD
Observable	$\mathrm{VAF}{\left(\mathrm{WT}2\right)}_{\max }$	Model size	${\mathrm{VAF}\left(\mathrm{WT}1\right)|}_r$	$\overline{\mathrm{NRMSE}}$ (id)	$\overline{\mathrm{NRMSE}}$ (val)
u	73.59%	38	98.08%	3.66%	3.75%

.

The dynamical properties of the ROM are summarized in

Table 2. Dynamical properties of IODMD${}_u$ model. The first column entries corresponds to the numerology chosen for the modes, which are ordered by increasing natural frequency. The second column entries correspond to their order in terms of eigenvectors of state matrix $\tilde{\mathbf{A}}$. The third and fourth columns have information of the modes natural frequencies and fifth column the damping ratio.

Modes		Dynamical Properties
#	DMD	${\omega }_n$ [Hz]	${\omega }_n$ St [-]	$\xi$
1	27	0	0	1
2	34	0.0014	0.0283	1
3	37–38	0.0030	0.0585	0.7608
4	35–36	0.0032	0.0637	0.2041
5	32–33	0.0064	0.1268	0.3072
6	30–31	0.0088	0.1735	0.0275
7	28–29	0.0109	0.2162	0.1192
8	25–26	0.0129	0.2559	0.0497
9	23–24	0.0158	0.3124	0.0712
10	21–22	0.0194	0.3830	0.0140
11	19–20	0.0207	0.4087	0.0278
12	17–18	0.0247	0.4884	0.0358
13	15–16	0.0263	0.5204	0.0280
14	13–14	0.0301	0.5959	0.0066
15	11–12	0.0343	0.6787	0.0025
16	9–10	0.0404	0.7996	0.0026
17	7–8	0.0472	0.9331	0.0100
18	5–6	0.0560	1.1071	0.0040
19	3–4	0.1479	2.9251	0.0091
20	1–2	0.2249	4.4473	0.0482

. The location of the poles dictating the dynamical properties of each of the modes can be visualised in Figure 6, where they have been labelled according to the numerology used in Table 2. Frequencies are provided in both continuous time (Hz) and as a dimensionless number (St) to investigate the possible existence of certain wake phenomena appeaering at certain frequencies.

Each of the eigenvalues in Figure 6 corresponds to an eigenvector of the high dimensional state matrix A, which, in turn, corresponds to a mode or coherent structure in space which inherits the dynamical properties indicated in Table 2. Figure 7, Figure 8 and Figure 9 allow to visualise these structures in different frequency ranges. Isosurfaces of various values have been chosen so that the structures can be visualised and interpreted qualitatively.

The modes, excited by certain inputs, will contribute to the overall evolution of the flow field and thus model the wake. It should be noted that these modes are merely representations of the high order eigenvectors and are not representative of the control action influence. Nevertheless, their study provides a deeper understanding of the wake dynamical properties.

The purpose of this work is to leverage the ideas behind the Koopman Operator in order to find coordinates where the existing non-linear dynamics can be better described. The three pathways mentioned in Section 2.2 are attempted. Having the hypothesis in mind that other variables may exist, several candidates are tested. The results are compiled in

Table 3. Models validation results using the IODMD algorithm with different observables. The information in the first 4 columns is the same as in Table 1. The last two columns show the improved performance in terms of fit and the improvements in terms of the size of model. Reduced models are less computationally expensive for the purposes of model predictive control, hence they are written in green. $\Delta$ VAF ${\left(\mathrm{WT}2\right)}_{\max }$ is the difference of fit with relation to the best model obtained using IODMD${}_u$. $\Delta$ size is the relative difference with relation to the best model obtained using IODMD${}_u$.

Model Properties	Model Obtained by IODMD_u			$\Delta$ IODMD_u
Observable	$\mathrm{VAF}{\left(\mathrm{WT}2\right)}_{\max }$	Model size	${\mathrm{VAF}\left(\mathrm{WT}1\right)|}_r$	$\Delta$ VAF ${\left(\mathrm{WT}2\right)}_{\max }$	$\Delta$ size
v	74.78%	39	97.64%	+1.19%	+2.63%
w	81.71%	39	96.99%	+8.12%	+2.63%
$u^{\prime }$	43.16%	29	93.99%	−30.43%	−23.68%
$v^{\prime }$	43.09%	49	96.71%	−30.50%	+28.95%
$w^{\prime }$	55.57%	49	96.69%	−18.02%	+28.95%
$u^{\prime 2}$	51.83%	57	95.58%	−21.76%	+50.00%
$w^{\prime 2}$	35.18%	104	94.78%	−38.41%	+41.32%
$u^2$	74.71%	24	96.76%	+1.12%	−36.84%
$v^2$	72.38%	112	99.33%	−1.21%	+194.74%
$w^2$	77.67%	106	99.29%	+4.08%	+178.95%
$u·v$	75.56%	59	99.03%	+1.97%	+55.26%
$u·w$	76.54%	35	97.42%	+2.95%	−7.89%
$v·w$	75%	120	99.25%	+1.41%	+63.07%
$u·v·w$	72.83%	72	99.25%	−0.76%	+89.47%
$u^2+w^2$	75.67%	24	96.79%	+2.08%	−36.84%
$u^2+v^2+w^2$	78.65%	24	96.89%	+5.09%	−36.84%

. The choosing of the observables is motivated by the Navier–Stokes non-linear Partial Differential Equations (NPDE) which govern the behaviour of a fluid and consequently of the travelling wake. The triplet $(u,v,w)$ corresponds to the velocity field, with u being the streamwise, v the spanwise and w the vertical velocity component. The velocity components can be separated into their mean and fluctuating parts $u=\overline{u}+u^{\prime }$, $v=\overline{v}+v^{\prime }$ and $w=\overline{w}+w^{\prime }$. The triplet $(\overline{u},\overline{v},\overline{w})$ refers to the base flow and the triplet $(u^{\prime },v^{\prime },w^{\prime })$ is computed by taking the difference between the velocity vector at each time instant and the base flow. The fluctuating terms may be interpreted as unsteady flow, similar to the terms appearing in the famous Reynold Average Navier–Stokes (RANS) equations. The choosing of the base flow is subjective and in the context of this work it corresponds to the fully developed flow for a blade pitch angle of 0 degrees. Subtracting a mean of the data reduces DMD to the temporal discrete Fourier transform, which may be restrictive and undesirable, as indicated in [38].The alternative approach of subtracting an equilibrium point of the dynamics appears to have advantages over using the mean and is attempted by removing the base flow.

Using the observable w shows gains of 8.12%, indicating that variations in the vertical velocity component may be more suited to describe the captured power dynamics at the downstream turbine. Furthermore, choosing an observable that includes measurements of all velocity components, such as $u^2+w^2$ or $u^2+v^2+w^2$ appears to be fruitful as well, showing increase in fit of 2.08% and 5.98% respectively, and relative decreases in the number of states of the best performing model of −36.85% for both cases. The results found show similarities with the ones in [34], where using the observables $u^2+v^2$ and $u^2+v^2+w^2$ yielded increases in fit of the downstream generator power output-albeit smaller for the case of Wake Redirection Control (WRC) by yaw misalignment, another wind farm control strategy-in the order of 1% when compared to the best performing IODMD${}_u$ model.

3.2. Extended Input Output Dynamic Mode Decomposition with Flow Field Information

In light of the second modification proposed in Section 2.2, the snapshot matrix containing data of the streamwise velocity component u is augmented to include other observables. The observables used to augment the original state vector are summarized in

Table 4. Models validation results using the EIODMD algorithm with different observables. The second observables column has the state chosen to augment the initial data matrix with u information. The last two columns represent the average NRMSE of the first observable, using both identification and validation data.

Model Properties		Model Obtained by EIODMD			$\Delta$ IODMD${}_{\mathit{u}}$
Observable	Second observables	$\mathrm{VAF}{\left(\mathrm{WT}2\right)}_{\max }$	Model size	${\mathrm{VAF}\left(\mathrm{WT}1\right)|}_r$	$\Delta$ VAF ${\left(\mathrm{WT}2\right)}_{\max }$	$\Delta$ size	$\overline{\mathrm{NRMSE}}$ (id)	$\overline{\mathrm{NRMSE}}$ (val)
u	v	79.00%	32	96.70%	+5.41%	−15.79%	-0.05%	−0.07%
u	w	81.74%	30	97.27%	+8.15%	−21.05%	−0.12%	−0.33%
u	v, w	79.76%	28	97.22%	+6.17%	−2.19%	−0.08%	−0.13%

.

The results indicate that augmenting the data matrix with additional information can have benefits. Augmenting u with w data translates into increases in fit of 8.15% when compared to the best performing model using IODMD${}_u$ and a reduction in the number of necessary states to accurately describe the downstream turbine generator power dynamics to 30. Furthermore, the EIODMD${}_{u,w}$ is capable of rebuilding the upstream turbine wake wit an average NRMSE of 3.42% for the validation data set. It should be noted that this process can become computationally quite expensive due to the systems dimensions, leading to data matrices containing approximately 80 million entries. A comparison of the IODMD${}_u$ and the EIODMD${}_{u,w}$ is portrayed in Figure 10, which shows the ability of the EIODMD${}_{u,w}$ to better capture the spikes in the downstream generator power as a consequence of the blades pitching back to zero degrees and again inducing modifications in the wake’s characteristics.

Both models can also be compared based on their response in the frequency domain. The EIODMD${}_{u,w}$ shows a resonance frequency ${\omega }_n$ at approximate 0.2 Strouhal, which can be visualised in Figure 11. The system’s frequency response shows at this frequency for the downstream generator power shows a magnitude of 10.43 dB, higher than any other frequency.

3.3. Extended Input Output Dynamic Mode Decomposition with Turbine Information

In light of the third modification proposed in Section 2.2, the snapshot matrix containing relevant information of the velocity components is augmented with turbine measurements, referred to as deterministic states. The rotor speeds of the upstream and downstream turbines are chosen as deterministic states, being hypothesized that these contain relevant information for the generator power dynamics and are usually available as measurements in wind turbines. The results are summarized in

Table 5. Models validation results using the EIODMD algorithm with different base observables and deterministic states.

Model Properties		Model Obtained by EIODMD			$\Delta$ IODMD${}_{\mathit{u}}$
Observable	Deterministic	$\mathrm{VAF}{\left(\mathrm{WT}2\right)}_{\max }$	Model size	${\mathrm{VAF}\left(\mathrm{WT}1\right)|}_r$	$\Delta$ VAF ${\left(\mathrm{WT}2\right)}_{\max }$	$\Delta$ size	$\overline{\mathrm{NRMSE}}$ (id)	$\overline{\mathrm{NRMSE}}$ (val)
u	${\mathsf{\Omega }}_1,{\mathsf{\Omega }}_2$	76.71%	12	98.97%	+3.12%	−68.42%	+0.17%	+0.45%
w	${\mathsf{\Omega }}_1,{\mathsf{\Omega }}_2$	76.79%	30	98.93%	+3.20%	−21.05%	-	-

.

The results indicate that augmenting data matrices with turbine information is a trade-off. On the one hand, the power predictions are more accurate and the ROMs have fewer states. On the other hand, the prediction of the wake is less accurate.

This is beneficial for the purposes of outputs prediction and control, as this procedure provides models which are better able to capture the generator power dynamics (when compared to the IODMD${}_u$ approach). Additionally, it does not incur on an extra computational cost, as only a few lines of data are added to the data matrices. Nevertheless, small increases in average NRMSE were registered to 3.83% and 4.20% for identification and validation data sets, respectively, hence the disadvantage for the purposes of wake behaviour prediction.

4. Conclusions and Future Work

The present work explored ideas motivated by non-linear dynamical systems theory in the context of reduced order modelling of wind farm systems. Data from high fidelity simulations where Axial Induction Control is tested are used to compute a linear state space Reduced Order Model (ROM).

It is firstly shown that a ROM with 38 states is able to reproduce the downstream turbine generator power with a fit of 74% and reconstruct the upstream turbine wake with an average NRMSE of 3.75% for validation data. The approach taken is quite straightforward, where the streamwise velocity component u is used and a simple mapping from collective pitch angle to generator power of both turbines is derived. Furthermore, the ROM can be investigated so that the relevant frequencies can be analysed and the physical mechanisms behind the observed gains in energy extraction studied.

Secondly, it is shown that taking other velocity components as states or functions of these can improve the model’s fit, indicating that for a given Wind Farm Control strategy there might exist some coordinate which is better suited to describe the observed dynamics. Using w led to a ROM able to predict the downstream turbine generator power with a VAF of 81.71%, 8.12% higher when compared to the baseline case where u is used. Non-linear observables, meaning non-linear functions of the velocity components provided improvements of 5.09% and 4.08%, using $w^2$ and $u^2+v^2+w^2$, respectively.

Thirdly, it is shown that by judiciously choosing a parsimonious set of observables to augment the data matrices improves the model’s fit. Initial tests showed benefits of using the observable w instead of u. Therefore, the adding of w to the data matrix containing u, as opposed to its replacement, led to an increase in 8.15% in fit as measured by VAF of the downstream generator power dynamics. Furthermore, this is obtained by a ROM with 30 states—eight states fewer than the IODMD${}_u$ model—despite using data matrices with increased size. The ability of the model to reconstruct the wind turbine wake is not affected, as the first state u is preserved. In fact, decreases in average NRMSE for the validation data set in the order of −0.33% were registered, indicating that this procedure improves both the model’s prediction capabilities of the output and state trajectory. A comparison of the IODMD${}_u$ and EIODMD${}_{u,w}$ response in the time and frequency domains further shows that the EIODMD${}_{u,w}$ model can better capture the spikes in the downstream turbine generator power and a new resonance peak at approximately 0.20 $St$ in the EIODMD${}_{u,w}$ frequency response appears. At this point it can only be hypothesized that wake dynamics of relevance appearing at this frequency are encoded in the w velocity component, which in turn are of use to describe the peaks registered at the downstream turbines power production as the upstream turbine blades pitch back and induce dynamic modifications in the wake as it hits the downstream turbine rotor plane.

Lastly, the procedure of augmenting the data matrices with turbine information led to increases in the model’s fit of the output in the order of approximately 3%. However, the prediction capabilities of the state trajectory worsened, showing an increase in average NRMSE of 0.45% for validation data.

It should be noted that the results found are valid for the operating regime and inflow conditions studied in the context of this work. Nevertheless, based on the results found, it is believed that a careful choosing of the observables to use for the systems identification task can have a tremendous impact on the model’s ability to recover the existing non-linear dynamics. The choosing of the observables should ideally be led by expert knowledge and a deeper understanding of the dynamics inherent to the Wind Farm Control being studied.

Additionally, it should be emphasized that different wind conditions and wind farm configurations may lead to more sophisticated scenarios and hence to a more complex modelling task. Different wind conditions, such as increased turbulence intensity and varying incoming wind speed direction may add to the non-linear behaviour of the wake and lead to oscillations in the produced power of both turbines. It can naturally be hypothesized that the prospects of using a linear reduced order model to capture the increasing non-linear behaviour inherent in these more sophisticated scenarios decrease. Still, different combinations of observables may be used, in an attempt to maintain the same level of accuracy as the simplified operating regime, analysed within the context of this work, although the trade-off between computational effort and increases in the accuracy of the model should be studied. The same principle applies to scenarios with multiple wind turbines and wake interaction phenomena. The large amount of information from a wind farm with four to five wind turbines may require the usage of a computer cluster to perform the data-driven techniques discussed, as the Singular Value Decomposition of the data matrices may become too costly for a normal computer. Nevertheless, techniques extending DMD to exploit low-rank tensor decompositions of high dimensional systems have been developed [39]. Such approach is of interest to explore within the context of IODMD in order to reduce the computational complexity and the amount of memory required to store the data, as the authors of the former article mention. As a final remark, it is worth mentioning that linear reduced order models will most likely see their wake prediction capabilities diminished in the more evolved scenarios as a consequence of increased non-linearities, although a clear quantification of these deviations has to be further investigated.

Future use of reduced order models may fall into three categories: pure study of the system dynamics, prediction and control. The first concerns increasing ones understanding of the dynamical phenomena caused by specific Wind Farm Control strategies, which can be done by analysing the frequencies dominating dynamical activity and the spatio temporal patterns associated with it. These patterns can then be compared to known wake behaviours, such as wake meandering, for example, which is associated with certain frequencies. This is of special interest as current research is increasingly interested in dynamic wind farm control. The second concerns predicting relevant information of wind farm systems, such as power, loads and the flow field. Accurate DMD-based ROMs may be used for real time applications, as evidence gathered in [31] suggests, where DMD outperformed POD in terms of efficiency due to its iteration-free feature. Based on the results of this work, the prospects of using ROMs for Wind Farm Control in real time are considerable. By designing accurate estimators, the ROM states can be estimated using wind turbine power measurements, which means that the knowledge of the wake at each time is not necessary to predict generator power in future instants. The third application concerns control. Many different control problems can be formulated, and the user may optimise one variable, such as power yield, or several simultaneously, including loads or even actively manipulating the wake. The development of optimal controllers is of course naturally dependent on the accuracy of the models. As an example, in [34], a Model Predictive Control (MPC) is developed based on a ROM and implemented in SOWFA in order to track a certain collective power reference by using yaw control to redirect the wake. The MPC developed, with a time horizon of 350 time samples and with minimum, maximum and rate constraints is able to solve the quadratic programming optimisation problem in a matter of seconds.

Reduced Order Models, such as the one presented here, should ultimately be used for real-time control strategies and optimization of wind power plants, aiming at reducing the cost of electricity production. Therefore, evaluating the accuracy of the presented ROM under more realistic wind farm conditions is of the essence. Other data sets including non uniform inflow profiles, terrain effects, non aligned wind turbines and wake overlap should be used to further assess the accuracy of the ROM using the described Koopman motivated IODMD algorithms. One question to further investigate is what observables will be beneficial to compute the ROM under more complex scenarios and evaluate if the addition of w to the data matrix continues to increase the fitness of the model as it did for the simpler case. Furthermore, it would also be of interest to investigate, in the context of more complex scenarios, the possibility of increasing the accuracy of the wake prediction (and possible reduction in the number of states needed) by including a simplified model for small-scale turbulence, as in [40], where evidence that the small-scale dynamics can be captured by adding a homogenous turbulent field to the developed POD model is presented. Other future research directions might focus on applying the data-driven modelling tools and techniques presented here to other control strategies, such as Individual Pitch Control (IPC), where the generator power dynamics may be revealed to be more complex. Moreover, an additional step is the design of predictive controllers based on the computed ROMs for axial induction control to maximize power production in a wind farm and implementation in high fidelity simulators.