On the potential of reduced order models for wind farm control: a Koopman dynamic mode decomposition approach

: 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 ﬁnding reduced order models which can accurately embed the existing non linearities. The work here presented 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 - 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 ﬁt 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.


Introduction
Renewable energy production has increased substantially in the last 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 later, 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 St * 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 3 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] 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 and a two turbine setting arises in [11], 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 reduced order model's capabilities 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.
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.

Input Output Dynamic Mode Decomposition
Dynamic Mode Decomposition (DMD) is a data driven technique which 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 [12]. 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 [13]. The standard DMD algorithm has been applied in different fields, such as neuroscience [14] and epidemiology [15], providing a deeper understanding of the systems being studied. Furthermore, different variations have surged since its first publication [16], 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 [17] and input-output information [18,19]. The later 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: where x k ∈ R n , u k ∈ R q and y k ∈ R l are the state, input and output vectors respectively. Matrices A ∈ R n×n , B ∈ R n×q , C ∈ R l×n and D ∈ R l×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: and rearranging, the reduced order state matrices are obtained by minimizing the error of the Frobenius norm [20]: whereÃ ∈ R r×r ,B ∈ R r×q ,C ∈ R l×r , andD ∈ R l×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 U * X' and ΣV: 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 x k = Ux k .
In the standard DMD algorithm, described in [21], an eigendecomposition is performed to the reduced order state space matrixÃ,ÃW = WΛ, where columns of W are eigenvectors and Λ is a diagonal matrix containing the corresponding eigenvalues λ 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 is explicit computation, an approximation is calculated based on the so called DMD modes, which correspond to the columns of Φ, Φ=UW, as originally proposed in [22]. 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.

Koopman Operator Theory
The Koopman operator, firstly introduced in 1931 [23] and recently resurged 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: where x(t) ∈ R n is a vector representing the state of the dynamical system at time t, µ 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 F t : M → M, x k+1 = F t (x k ). A function defined as g : M → R is called an observable of the system. The Koopman operator K t acts on observable functions g as K t = g • F t where • is the composition operator, so that: 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 x k ∈ 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 Kϕ k = λ k ϕ k , the functions ϕ k (x) are called Koopman eigenfunctions and they define a set of intrinsic measurement coordinates on which it is Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 16 October 2020 doi:10.20944/preprints202010.0352.v1 possible to advance these measurements with a linear dynamical system. The new measurement space formed by a vector of observables g(x) can be expressed in terms of the Koopman eigenfunctions: where v k is the kth Koopman mode of the map F corresponding to the observable 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 ϕ k of K where 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).
Future solutions can be computed by multiplication with the Koopman eigenvalues, λ k ∈ C, which characterize the temporal behaviour of the corresponding Koopman mode v k . The connections of DMD to the underlying non linear dynamics of a system through the Koopman operator theory were firstly established in [24], where it is shown 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 [25,26]. The later 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 [27,28], 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 [28] and proved to yield better results when compared to POD for prediction purposes [29]. The linear representation provided by the Koopman Operator is not only useful for prediction but also for control, as explored in the works of [30,31], 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 Subsection 2.1, only the snapshot matrices organisation step differs, with the data matrices being generalized.

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) [33]. 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 [33].
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 5 rotor diameters (5D) behind the upstream turbine. The full domain as well as inflow characteristics and turbulence intensity used for the CFD simulation are specified in [34]. The turbine used for the simulation is the DTU 10MW turbine [35], 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 seconds. 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 θ t1 . These serve as inputs Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 16 October 2020 doi:10.20944/preprints202010.0352.v1 varying 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 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, is depicted in Figure 1.

Results
In this section, the techniques introduced in Subsection 2.1 and Subsection 2.2 are used to obtain the desired input-output mapping illustrated in Figure 1. 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. 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 x k = Ux k and the streamwise velocity component u can be entirely reconstructed and compared to SOWFA validation measurements, as represented in Figure 4 and   Taking all time instants into account, the deviations can be summarized resorting to an average NRMSE. The best performing model results are summarized in Table 1. The dynamical properties of the ROM are summarized in Table 2. 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.  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. Figures 7, 8 and 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 Subsection 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. 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 The spatial patterns present here appear to also explain some of the dynamics more related with the downstream turbine, as opposite to background and low frequency modes. velocity components can be separated into their mean and fluctuating parts u =ū + u , v =v + v and w =w + w . The triplet (ū,v,w) refers to the base flow and the triplet (u , v , w ) 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 [36].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 [32], 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. Table 3. Models validation results using 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. ∆ VAF (WT2) max is the difference of fit with relation to the best model obtained using IODMD u . ∆ size is the relative difference with relation to the best model obtained using IODMD u . In light of the second modification proposed in Subsection 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. Table 4. Models validation results using 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. 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, where it is shown 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 ω 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. In light of the third modification proposed in Subsection 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. The results indicate that augmenting data matrices with turbine information is a trade-off. For the purposes of outputs prediction and control, the procedure provides models which are better able to capture the generator power dynamics when compared to the IODMD u approach. Furthermore, this procedure does not incur on an extra computational cost, as only a few lines of data are added to the snapshot matrices. Nevertheless, small increases in average NRMSE were registered to 3.83% and 4.20% for identification and validation data sets, respectively.

Conclusions and Future Work
The present work explored ideas motivated by non linear dynamical systems theory in the context of reduced order data driven modelling of wind farm systems. Data from high fidelity simulations where Axial Induction Control is tested is 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 it solemn 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, 8 states less when compared to 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 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 further domains 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.
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.
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 reveal to be more complex. Moreover, an additional step is the design of predictive controllers based on the computed ROMs to maximize power production in a wind farm and implementation in high fidelity simulators.