Improving wind farm flow models by learning from operational data

This paper describes a method to improve and correct an engineering wind farm flow model by using operational data. Wind farm models represent an approximation of reality and therefore often lack accuracy and suffer from unmodeled physical effects. It is shown here that, by surgically inserting error terms in the model equations and learning the associated parameters from operational data, the performance of a baseline model can be improved significantly. Compared to a purely data-driven approach, the resulting model encapsulates prior knowledge beyond that contained in the training data set, which has a number of advantages. To assure a wide applicability of the method – also including existing assets – learning here is purely driven by standard operational (SCADA) data. The proposed method is demonstrated first using a cluster of three scaled wind turbines operated in a boundary layer wind tunnel. Given that inflow, wakes, and operational conditions can be precisely measured in the repeatable and controllable environment of the wind tunnel, this first application serves the purpose of showing that the correct error terms can indeed be identified. Next, the method is applied to a real wind farm situated in a complex terrain environment. Here again learning from operational data is shown to improve the prediction capabilities of the baseline model.

1 Introduction

Knowledge of the flow at the rotor disk of each wind turbine in a wind power plant enables several applications, including wind farm control, the provision of grid services, predictive maintenance, the estimation of life consumption, the feed-in to digital twins, and power forecasting, among others.

This paper describes a new method to improve a wind farm flow model directly from standard operational data. The main idea pursued here is to use an existing wind farm flow model to provide a baseline predictive capability; however, as all models contain approximations and may lack the description of some physical phenomena, the baseline model is improved (or “augmented”, which is the term used in this work) by adding parametric correction terms. In turn, these extra elements of the model are learned by using operational data. The correction terms capture effects that are typically not present in standard flow models (such as, for example, secondary steering, Fleming et al., 2018; or wind farm blockage, Bleeg et al., 2018) or that are highly dependent on a specific site or difficult to model upfront (such as, for example, nonuniform inflow caused by local orography and vegetation).

Various wind farm flow models have been developed and are described in the literature. Whereas direct numerical simulation (DNS) is still out of reach for practical applications due to its overwhelming computational cost, large-eddy simulation (LES) methods are now routinely used for the modeling of wind farm flows (Fleming et al., 2014; Breton et al., 2017). Although invaluable for the understanding of the behavior of the atmospheric boundary layer and of wakes, LES is however still very expensive, so that its use outside of some specialized applications is limited. To reduce cost, one can resort to lower-fidelity computational fluid dynamics (CFD) models (Boersma et al., 2019), or to the extraction of reduced-order models (ROMs) from higher-fidelity ones (Bastine et al., 2014). Instead of deriving models from first principles, another widely adopted approach is to use engineering models, which are expressed in the form of parametric analytical formulas with a limited number of degrees of freedom and hence a much reduced numerical complexity (Frandsen et al., 2006; Gebraad et al., 2014; Bastankhah and Porté-Agel, 2016). The present paper uses this last family of methods, although ideas similar to the ones developed here could also be applicable to higher-fidelity models.

Even though engineering models are constantly improved and refined (Fleming et al., 2018), they will most likely always exhibit only a limited accuracy in many practical applications, for example whenever an important role is played by effects such as orography, (seasonal) vegetation, spatial variability of the wind, sea state roughness, the erection of other neighboring wind turbines, the presence of obstacles, and others. In addition, low-fidelity models often lack some physics, e.g. the flow acceleration caused by wake and rotor blockage, secondary steering, or others. The idea pursued in this paper is then to take a rather pragmatic approach: based on the realization that it will always be difficult – if not altogether impossible – to include all effects and all physics in a model of limited numerical complexity, a given model is corrected by unknown parametric terms, which are then learned by using operational data.

The idea of improving an existing model based on measurements is hardly new, and it is actually an important topic in the areas of controls and system identification. For example, in the field of wind farm flows, a Kalman filtering approach has been proposed by Doekemeijer et al. (2017) to update model predictions based on lidar measurements. Here again the present paper takes a more pragmatic approach, and model updating is based exclusively on data provided by the standard supervisory control and data acquisition (SCADA) systems that are typically available on contemporary wind turbines. On the one hand this has the advantage that the proposed method is applicable to existing assets, as it does not necessitate extra sensors. On the other hand, given that stored SCADA data typically represent 10 min averages, this also implies that the models obtained by this technique are of a steady-state nature. Although unsteady effects in wind farms are clearly important, steady-state models are still very valuable and can support many of the applications listed above. In addition, nothing prevents the generalization of the proposed approach to unsteady flow models, assuming that the relevant higher-frequency data sets are available, which is already the subject of ongoing work from these authors.

The contemporary literature – and not only in the field of wind energy – indicates an increasing interest in data-driven approaches. Just to give one single example related to wake modeling, a purely data-driven approach has been recently described by Göçmen and Giebel (2018). However, the current enthusiasm for data should not make one forget that physics-based and analytical models are also extremely valuable because they often encapsulate significant knowledge on a given problem, often corroborated by long experience. In fact, purely data-driven approaches suffer from a number of limitations that descend directly from a very simple and inevitable fact: a model that is exclusively based on data can only know what is contained in the data set that was used to build it. Typically, this means that a very significant amount of data is necessary to obtain a model that is sufficiently general and accurate. Furthermore, the data have to cover the entire spectrum of operation of the system. This also means that the model might have very poor knowledge (and hence poor performance) for rare situations or conditions that take place at the boundaries of the operating envelope, where few if any data points might be available.

An alternative to the purely data-driven approach is presented in this work, where a reference baseline model is augmented with parametric error terms, which are then identified using data. The baseline model already includes prior knowledge based on physics, empirical observations, and experience. Therefore, even prior to the use of data, a minimum performance can be guaranteed. The model is augmented with parametric error terms, whose choice is driven by physics and the knowledge of the limitations of the baseline model. Once the errors are identified using operational data, their inspection can clarify the causes of discrepancy between model and measurements. Eventually, this can be used to improve the underlying baseline model. Furthermore, by looking at the magnitude of the identified errors, significant deviations from the baseline model can be flagged to highlight issues with the model itself, the data, or the training process.

Finally, it should be noted that the identification of the error terms can be combined with the tuning of the parameters of the baseline model. This addresses yet another problem: tuning the parameters of a model that lacks some physics may lead to unreasonable values for the parameters, as the model is “stretched” to represent phenomena that it does not contain. By the proposed hybrid approach, the simultaneous identification of the parameters of the baseline model together with the ones of the error terms eases this problem, as unmodeled phenomena can be captured by the model-augmenting terms, thereby reducing the chances of nonphysical tuning of the baseline parameters.

The baseline model parameters and the extra correction terms have a different functional form in the augmented governing equations. Hence, they should be distinguishable from each other, as they imply different effects on the model. However, as for many identification problems, it is in general not possible to guarantee that all unknown parameters are observable and noncollinear given a set of measurements and, hence, given a certain informational content. To address this problem, the method proposed by Bottasso et al. (2014a) is used here, where the original unknown parameters are recast into a new set of statistically uncorrelated variables by using the singular value decomposition (SVD) of the inverse Fisher information matrix. Once the problem has been solved in the space of the orthogonal uncorrelated parameters, the solution is mapped back onto the original physical space. This approach not only avoids the ill-posedness of the original problem, but also allows one to clarify which physical parameters are visible given a certain data set.

The paper is organized as follows. First, the baseline model is introduced in Sect. 2.1, together with a detailed description of the proposed parametric corrections in Sect. 2.2. Next, the SVD-based parameter identification method is presented in Sect. 2.3. The approach is then applied in Sect. 3.1 to a cluster of scaled wind turbines operating in the atmospheric test section of the wind tunnel of the Politecnico di Milano (Bottasso et al., 2014b). The goal of this first application is to show that a correct identification of the error terms can be achieved. This is indeed possible in the controllable and repeatable conditions of a wind tunnel, where inflow and wake characteristics can be precisely measured, something that is hardly possible today in the field. Specifically, it is shown that the method can correctly learn the lack of uniformity of the wind tunnel inflow, which is akin to what happens in a real wind farm because of orographic effects. Similarly, it is shown that secondary steering, which is completely absent from the baseline model used here, can be learned by using turbine power measurements only. A more extended view on the wind tunnel results is reported in Appendix A. After having demonstrated the method in the known and controlled wind tunnel environment, a second application is developed in Sect. 3.2 that targets a real 43-turbine wind farm. Here results indicate that the augmented model has a markedly improved prediction capability when compared to the baseline one, thanks primarily to the identification of orographic effects on the inflow and the tuning of other model parameters. Finally, conclusions are drawn in Sect. 4.

2 Methods

2.1 Baseline wind farm flow model

The proposed method is applied here to the baseline wake model of Bastankhah and Porté-Agel (2016), implemented within the FLORIS framework (Doekemeijer and Storm, 2018). Given ambient wind conditions, steady-state velocities within a wind farm can be computed by this model, together with the corresponding operating states and power outputs of all its turbines. First, ambient conditions are estimated from un-waked machines operating in free stream, which are identified by the turbine yaw orientations and the wake model (Schreiber et al., 2018). Then, power and thrust of the upstream turbines are computed based on the turbine aerodynamic characteristics, regulation strategy, and alignment with the local wind direction. Next, the wakes shed by these turbines are calculated in terms of their trajectory and speed deficit. In turn, this yields the velocity at the rotor disks of the turbines immediately downstream. In the case of multiple wake impingements on a rotor, a combination model is used to superimpose multiple wake deficits. Similarly, an added turbulence model is used to estimate the turbulence intensity at a downstream turbine rotor disk, as this local ambient parameter affects the expansion of the wake. This process is repeated marching downstream throughout the wind farm until the last downstream turbine is reached.

In this work, the implementation uses the selfSimilar FLORIS velocity deficit model, the rans deflection model, the quadraticRotorVelocity wake combination model, and the crespoHernandez added turbulence model. The interested reader is referred to Bastankhah and Porté-Agel (2016), Crespo and Hernández (1996), and Doekemeijer et al. (2019) and references therein for detailed descriptions and derivations of these models.

Engineering wake models depend on a number of parameters, which should be tuned in order to obtain accurate predictions. For the specific model used in this work, these tunable factors are the wake parameters α, β, k_a, k_b, a_d, and b_d and the turbulence model parameters TI_a, TI_b, TI_c, and TI_d (Bastankhah and Porté-Agel, 2016).

In this work, the parameters are first set to an initial value, either taken from the literature or identified with ad hoc measurements; these initial values are held fixed throughout the analysis and not changed further. Corrections to the initial values are then expressed as

$$\begin{array}{}\text{(1)}& k=k^{\ast }+p_k,\end{array}$$

where k is a model parameter, k^∗ its initial value, and p_k the correction. Although this is not strictly necessary, this redundant notation helps highlight the changes to the nominal model parameters obtained by the proposed procedure.

2.2 Model augmentation

The engineering model described earlier is a rather simple approximation of a flow through a wind power plant and it is therefore bound to have only a limited fidelity to reality, with a consequently only limited predictive accuracy. Even for more sophisticated future models, it is difficult to imagine that all relevant physics will ever be precisely accounted for. But even if such a model existed, in practice one might simply not have all necessary detailed information on the relevant boundary and operating conditions that would be required. For example, one might not know with precision the conditions of the vegetation around and within a wind farm, with its effects on roughness and, hence, on the flow characteristics. In other words, it is safe to assume that all models are in error to some extent and probably always will be.

To address this problem, the model can be pragmatically augmented with correction terms. Here one could take two alternative approaches: either a generic all-encompassing error term is added to the model or “surgical” errors are introduced at ad hoc locations in the model to target specific presumed deficiencies. The first approach could be treated with a brute-force parametric modeling approach, for example by using a neural network. Here, the second approach was used, as it allows for more insight into the nature of the identified corrections. The specific parametric corrections used in the present paper are reviewed next. It is clear that these are only some of the many corrections that could be applied to the present baseline model, so that the following does not pretend to be a comprehensive treatment of the topic. Nonetheless, results indicate that some of these corrections are indeed significant and provide for a marked improvement of the baseline model.

• Nonuniform inflow. The inflow to a wind farm can exhibit spatial variability, mostly because of orographic and local effects, especially in complex terrain conditions. For example, commercial wind resource assessment tools include topographic speedup ratios customarily computed by CFD models (Jacobsen, 2019). In contrast to this established practice, no direct or equivalent modeling of orographic effects is at present available in engineering wake models. Another reason for inflow variability may be due to wind farm blockage effects (Bleeg et al., 2018). Indeed, current wake models such as the one used here assume that upstream turbines affect downstream ones through their wakes but do not model the effects of downstream machines on the upstream ones. In a wind farm, depending on the wind direction and cross-wind location considered, the number and operating state of downstream turbines vary, which may induce a cross-wind speed variability in the inflow.

  To capture some of these effects, the model ambient flow speed V_∞ is expressed here as a function of height above ground Z, cross-wind lateral position Y, and ambient wind direction Γ as

  $$\begin{array}{}\text{(2)}& \begin{array}{rl}& V_{\mathrm{\infty }}(Y,Z,\mathrm{\Gamma })=\\ & \left(\mathrm{1}+f_{\mathrm{augm},\mathrm{speed}}\left(Y,\mathrm{\Gamma },{\mathit{c}}_{\mathrm{speed}},{\mathit{p}}_{\mathrm{speed}}\right)\right)\\ & V_{\mathrm{\infty },\mathrm{0}}{\left({\displaystyle \frac{Z}{z_{\mathrm{h}}}}\right)}^{{\mathit{\alpha }}_{\mathrm{vs}}},\end{array}\end{array}$$

  where V_∞,0 is the reference (baseline uncorrected) ambient flow speed and z_h the reference height of the vertically sheared flow with exponent α_vs. Function $f_{\mathrm{augm},\mathrm{speed}}(Y,\mathrm{\Gamma },{\mathit{c}}_{\mathrm{speed}},{\mathit{p}}_{\mathrm{speed}})$ is the speed correction term. This function is defined in the 2D space $Y\in [Y_{min},Y_{max}]$, $\mathrm{\Gamma }\in [{\mathrm{\Gamma }}_{min},{\mathrm{\Gamma }}_{max}]$. For each value of the ambient wind direction Γ, Y is a lateral coordinate orthogonal to it that spans the width of the farm; hence, by selecting Γ_min and Γ_max a lateral inflow nonuniformity can be modeled for a given sector or the whole wind rose of directions. The (Y,Γ) space is discretized into rectangular cells with corner nodes ${\mathit{c}}_{\mathrm{speed}}=[\mathrm{\dots };(Y_i,{\mathrm{\Gamma }}_i);\mathrm{\dots }]$ (for an example, see Fig. 16). The corresponding unknown error nodal values are stored in vector p_speed, and bilinear shape functions interpolate the error in each cell based on the nodal values at its corners. Equation (2) could be extended to also include a longitudinal wind-aligned coordinate, similarly to the localized speedup ratios of Jacobsen (2019), to model wind farm blockage effects.

  Local orographic effects and blockage may also induce variability in the wind direction Γ. Similarly, the vertical shear exponent α_vs and turbulence intensity I may vary, for example on account of nonuniform roughness induced by vegetation or other obstacles. To include these effects in the farm flow model, the baseline quantities are augmented as

  $$\begin{array}{}\text{(3a)}& {\displaystyle }\mathrm{\Gamma }\left(Y\right)={\mathrm{\Gamma }}_{\mathrm{ref}}+Y{f}_{\mathrm{augm},\mathrm{dir}}\left({\mathrm{\Gamma }}_{\mathrm{ref}},{\mathit{c}}_{\mathrm{dir}},{\mathit{p}}_{\mathrm{dir}}\right),\text{(3b)}& {\displaystyle }{\mathit{\alpha }}_{\mathrm{vs}}\left(\mathrm{\Gamma }\right)={\mathit{\alpha }}_{\mathrm{vs},\mathrm{ref}}+f_{\mathrm{augm},\mathrm{shear}}\left(\mathrm{\Gamma },{\mathit{c}}_{\mathrm{shear}},{\mathit{p}}_{\mathrm{shear}}\right),\text{(3c)}& {\displaystyle }I\left(\mathrm{\Gamma }\right)=I_{\mathrm{ref}}+f_{\mathrm{augm},\mathrm{I}}\left(\mathrm{\Gamma },{\mathit{c}}_{\mathrm{I}},{\mathit{p}}_{\mathrm{I}}\right).\end{array}$$

  In all these expressions, (⋅)_ref indicates a baseline reference quantity, while function $f_{\mathrm{augm},(\cdot )}$ is a correction term. This function is defined on the 1D space $\mathrm{\Gamma }\in [{\mathrm{\Gamma }}_{min},{\mathrm{\Gamma }}_{max}]$, discretized with nodes ${\mathit{c}}_{(\cdot )}=[\mathrm{\dots };{\mathrm{\Gamma }}_i;\mathrm{\dots }{]}_{(\cdot )}$, using linear shape functions to interpolate the corresponding nodal values p_(⋅). Here again, by selecting Γ_min and Γ_max, corrections can be applied to the whole wind rose or just to a sector.

• Secondary steering. By misaligning a wind turbine rotor with respect to the incoming flow direction, the rotor thrust force is tilted, thereby generating a cross-flow force that laterally deflects the wake. As shown with the help of numerical simulations by Fleming et al. (2018), this cross-flow force induces two counter-rotating vortices that, combining with the wake swirl induced by the rotor torque, lead to a curled wake shape. As observed experimentally by Wang et al. (2018), the effects of these vortices result in additional lateral flow speed components, which are not limited to the wake itself but also extend outside of it. By this phenomenon, the flow direction within and around a deflected wake is tilted with respect to the upstream undisturbed direction. Therefore, when a turbine is operating within or close to a deflected wake, its own wake undergoes a change of trajectory – termed secondary steering – induced by the locally modified wind direction. Although models of this phenomenon are being developed (Martínez-Tossas et al., 2019), they significantly increase the computational cost and are not yet available in standard implementations of engineering wake models such as the one used here.

  The change of wind direction ΔΓ at a downstream turbine induced by secondary steering (indicated by the subscript ss) is modeled here as

  $$\begin{array}{}\text{(4)}& \mathrm{\Delta }\mathrm{\Gamma }\left(y\right)=f_{\mathrm{augm},\mathrm{ss}}\left(\stackrel{\mathrm{̃}}{y},{\mathrm{\Gamma }}_{\mathrm{init}},{\mathit{p}}_{\mathrm{ss}}\right),\end{array}$$

  where f_augm,ss is the correction term and $\stackrel{\mathrm{̃}}{y}=Y-y_{\mathrm{wc}}$ is the lateral distance to the wake centerline (see Fig. 1), defined in the baseline wind farm model as the locus of the points of minimum flow speed. According to the notation used in Eq. (6.12) of Bastankhah and Porté-Agel (2016), Γ_init indicates the initial wake direction of the closest upstream turbine. The correction term is expressed as the difference of two Gaussian functions and more precisely

  $$\begin{array}{}\text{(5)}& \begin{array}{rl}& f_{\mathrm{augm},\mathrm{ss}}\left(\stackrel{\mathrm{̃}}{y},{\mathrm{\Gamma }}_{\mathrm{init}},{\mathit{p}}_{\mathrm{ss}}\right)=\\ & {\mathrm{\Gamma }}_{\mathrm{init}}\left(p_{\mathrm{ss},\mathrm{1}}\exp \left(-\mathrm{0.5}{\left({\displaystyle \frac{\stackrel{\mathrm{̃}}{y}+\mathrm{sgn}\left({\mathrm{\Gamma }}_{\mathrm{init}}\right)p_{\mathrm{ss},\mathrm{3}}}{p_{\mathrm{ss},\mathrm{2}}}}\right)}^{\mathrm{2}}\right)\right.\\ & \left.-p_{\mathrm{ss},\mathrm{4}}\exp \left(-\mathrm{0.5}{\left({\displaystyle \frac{\stackrel{\mathrm{̃}}{y}+\mathrm{sgn}\left({\mathrm{\Gamma }}_{\mathrm{init}}\right)p_{\mathrm{ss},\mathrm{6}}}{p_{\mathrm{ss},\mathrm{5}}}}\right)}^{\mathrm{2}}\right)\right),\end{array}\end{array}$$

  where ${\mathit{p}}_{\mathrm{ss}}=(p_{\mathrm{ss},\mathrm{1}},p_{\mathrm{ss},\mathrm{2}},p_{\mathrm{ss},\mathrm{3}},p_{\mathrm{ss},\mathrm{4}},p_{\mathrm{ss},\mathrm{5}},p_{\mathrm{ss},\mathrm{6}})$ is the vector of free parameters, where parameters 1 and 4 are related to the amplitude, 3 and 6 to the standard deviation, and 2 and 5 to the location of the correction functions. Since the Gaussian functions are not centered at the wake centerline and the effect of secondary steering is assumed to be symmetric with respect to the misalignment angle, the correction term also depends on the direction of wake deflection sgn(Γ_init).

  This particular choice of the shape functions is motivated by the results shown in Fig. 8b of Wang et al. (2018). Indeed, LES simulations and measurements reveal the presence of a stronger lateral velocity component directed towards the wake on the leeward side of the wake itself, and of an opposite and weaker lateral component on the windward side. Such a distribution can be approximated by two Gaussian functions using Eq. (5).

  Note that the change in local wind direction also leads to a slight lateral deflection of the nonuniform wind farm inflow introduced previously. More precisely, for a turbine that is located ΔX behind an upstream turbine, the nonuniform inflow expressed by Eq. (2) is evaluated at Y+ΔXsin (ΔΓ) instead of Y.

  Figure 1a shows the hub height flow speed for two wind turbines modeled in FLORIS, with the turbine rotor disks being indicated with thick black lines. The wake centerlines and the undisturbed free-stream wind direction are indicated by black dotted and dashed lines, respectively. The upstream turbine is misaligned with respect to the incoming flow, and therefore its wake is deflected laterally. Using the baseline wake model, the downstream turbine wake develops along the free-stream wind direction. Panel (b) of the same figure shows the effects of the secondary steering correction term given by Eq. (5). The plot clearly shows that the downstream turbine wake path is affected by the locally changed wind direction.

• Non-Gaussian wake and flow acceleration. Engineering wake models are based, among other hypotheses, on assumed shapes of the speed deficit. For example, the present baseline model assumes a Gaussian distribution of the speed deficit within the wake. Another assumption is that the flow outside the wake is undisturbed and equal to the free stream. However, these assumptions can, at times, not be exactly satisfied, as already observed by Xie and Archer (2017) and Martínez-Tossas et al. (2019), among others. For example, aisle jets are local accelerations of the flow outside of the wake, produced by local blocking in the neighborhood of an operating turbine. It has been reported that aisle jets can induce local flow speedups in excess of 10 % of the undisturbed inflow (Dörenkämper et al., 2015).

  To account for such effects, the wake velocity V_wake of the baseline model is corrected as

  $$\begin{array}{}\text{(6)}& \begin{array}{rl}& V_{\mathrm{wake}}\left(d_{\mathrm{wc}}\right)=\\ & V_{\mathrm{wake},\mathrm{FLORIS}}\left(d_{\mathrm{wc}}\right)\left(\mathrm{1}+f_{\mathrm{augm},\mathrm{acc}}\left(d_{\mathrm{wc}},{\mathit{c}}_{\mathrm{acc}},{\mathit{p}}_{\mathrm{acc}}\right)\right),\end{array}\end{array}$$

  where V_wake,FLORIS is the baseline Gaussian wake speed profile, d_wc is the absolute distance to the wake center (which, at hub height, is equivalent to $\left|\stackrel{\mathrm{̃}}{y}\right|$), and f_augm,acc represents the correction term, which – similarly to the previous corrections – is modeled with linear shape functions characterized by node locations c_acc (in terms of d_wc) and nodal values p_acc.

• Reduced power extraction due to nonuniform wind turbine inflow. Numerical simulations conducted in FAST (Jonkman and Jonkman, 2018) using its blade element momentum (BEM) implementation yielded a slight reduction in the rotor power coefficient for horizontally sheared flow, when compared to unsheared conditions with the same hub wind speed. Even though BEM can only give a rough indication for such an effect, a correction of the power coefficient of the baseline model is introduced here in the form

  $$\begin{array}{}\text{(7)}& C_{\mathrm{P}}=C_{\mathrm{P},\mathit{\kappa }=\mathrm{0}}\left(\mathrm{1}+p_{\mathit{\kappa }}{\mathit{\kappa }}^{\mathrm{2}}\right),\end{array}$$

  where $C_{\mathrm{P},\mathit{\kappa }=\mathrm{0}}$ is the nominal power coefficient, κ the equivalent horizontal linear shear coefficient on the rotor disk, and p_κ the free correction parameter. The linear shear κ is either due to a lack of lateral uniformity of the inflow or due to the impingement of a wake, and it is evaluated accordingly within the farm model.

• Wind-speed-dependent power loss in yaw misalignment. The baseline formulation models the power extraction of a misaligned wind turbine using the cosine law $C_{\mathrm{P}}\left(\mathit{\gamma }\right)=C_{\mathrm{P}}\cos (\mathit{\gamma }{)}^{p_{\mathrm{P}}}$, where C_P is the power coefficient of the wind-aligned turbine, γ the misalignment angle with respect to the local flow direction, and p_P the power loss exponent. Different power loss exponents have been reported in the literature, ranging from the value of 1.4 found by Fleming et al. (2017) to 1.8 according to Schreiber et al. (2017), 1.9 for Gebraad et al. (2015), and all the way to the ideal value of 3 that is expected if only the rotor-orthogonal ambient flow component contributes to power extraction (Boersma et al., 2019). In addition, p_P might also depend on the regulation strategy used by the turbine controller. Here, the power coefficient in misaligned operation is augmented as

  $$\begin{array}{}\text{(8)}& C_{\mathrm{P}}=C_{\mathrm{P}}\cos {\left(\mathit{\gamma }+p_{\mathrm{P}\mathrm{0}}\right)}^{p_{\mathrm{P}}+p_{\mathrm{P},\mathrm{a}}\left(V-V_{\mathrm{rated}}\right)+p_{\mathrm{P},\mathrm{b}}},\end{array}$$

  where C_P is the power coefficient of the flow-aligned turbine (possibly reduced by shear effects, as argued above), p_P0 is the misalignment angle at which the turbine produces maximum power, and V and V_rated are, respectively, the rotor effective and rated wind speeds. Finally, p_P is the baseline exponent, while p_P,a and p_P,b are free parameters that model a linear wind speed dependency of the cosine law.

Figure 1 Effect of secondary steering on the trajectory of a downstream turbine. (a) Baseline wake model; (b) baseline model augmented with the empirical correction term of Eq. (5).

2.3 Parameter identification method

The parameters of the baseline model and of its correction terms are identified with the method developed by Bottasso et al. (2014a). The formulation of the parameter estimation problem is independent of whether the parameters belong to the baseline model or to its correction factors. In this sense, one can use the same method to just tune the baseline parameters without considering the correction terms, just identify the correction terms at the frozen baseline model, or concurrently identify both sets.

The formulation is based on the classical likelihood function, which describes the probability that a given set of noisy observations can be explained by a specific set of model parameters. By numerically maximizing this function, a set of parameters is identified that most probably explains the measurements. Bound constraints are used to guide the process and ensure convergence to meaningful results.

The accuracy with which the parameters can be estimated depends on how flat the likelihood function is with respect to changes in the parameters. For example, a flat maximum of the function implies that different nearby values of the model parameters are associated with similar values of the likelihood. These characteristics of the solution space are captured by the Fisher information matrix, which can be interpreted as a measure of the curvature of the likelihood function. Furthermore, it can be shown that the variance of the estimates is bound from below (Cramér–Rao bound) by the inverse of the Fisher matrix (Jategaonkar, 2015). Although the analysis of the Fisher information is useful for the understanding of the well-posedness of an estimation problem and of the quality of the identified model, it does not offer a constructive way of reformulating a given ill-posed problem. Indeed, a flat solution space and collinear parameters are to be expected in the present case, given the complex couplings and dependencies that may exist among the various parameters of a wind farm flow model and its correction terms.

To overcome this limitation of the classical maximum likelihood formulation, following Bottasso et al. (2014a), the original physical parameters of the model are transformed into an orthogonal parameter space, by diagonalizing the Fisher matrix using the SVD. This way, as the parameters are now statistically decoupled, one can set a lower observability threshold and in the analysis retain only the ones that are in fact observable given the available set of measurements. Once the problem is solved, the uncorrelated parameters are mapped back onto the original physical space.

As shown later on, this approach achieves multiple goals: it allows one to successfully solve a maximization problem with many free parameters, some of which might be interdependent on one another or not observable in a given data set; it reduces the problem size, retaining only the orthogonal parameters that are indeed observable; it highlights, through the singular vectors, the interdependencies that may exist among some parameters of the model, which provides for a useful interpretation tool that may guide the reformulation of parts of the model and its correction terms.

2.3.1 Maximum likelihood estimation of model parameters

A steady-state wind farm model can be mathematically expressed as

$$\begin{array}{}\text{(9)}& \mathit{y}=\mathit{f}(\mathit{p},\mathit{u}),\end{array}$$

where $\mathit{f}(\cdot ,\cdot ,\cdot )$ is the nonlinear static function describing the wind farm model, which depends on free parameters p∈ℝⁿ. These parameters can include both wake model parameters and/or model augmentation parameters. The model inputs $\mathit{u}\in {\mathbb{R}}^{n_{\mathrm{u}}}$ include ambient wind conditions (i.e. ambient wind speed, direction, air density, turbulence intensity) and control inputs (i.e. yaw misalignment, partialization factor, blade pitch, rotor speed of each turbine). The model outputs y∈ℝ^m represent quantities of interest for which measurements are available, in the present work these being the power outputs of each wind turbine in the farm. Experimental observations z of the simulated outputs y will in general result in a residual r∈ℝ^m, caused by measurement and process noise (e.g. plant–model mismatch), so that

$$\begin{array}{}\text{(10)}& \mathit{z}=\mathit{y}+\mathit{r}.\end{array}$$

Given a set $S=\mathit{\{}{\mathit{z}}_{\mathrm{1}},{\mathit{z}}_{\mathrm{2}},\mathrm{\dots },{\mathit{z}}_N\mathit{\}}$ of N independent observations, the likelihood function (Jategaonkar, 2015) can be defined as

$$\begin{array}{}\text{(11)}& \mathcal{L}\left(S{|}_{\mathit{p}}\right)=\prod _{i=\mathrm{1}}^{N}p\left(z_i{|}_{\mathit{p}}\right),\end{array}$$

where p(⋅) is the probability of S given p. Assuming the residuals r with covariance R to be statistically independent within the set of measurements (i.e. $E\left[{\mathit{r}}_i{\mathit{r}}_j^T\right]=\mathbf{R}{\mathit{\delta }}_{i,j}$, where δ_i,j is the Kronecker delta), the likelihood function can be written, following Jategaonkar (2015), as

$$\begin{array}{}\text{(12)}& \begin{array}{rl}& \mathcal{L}\left(S{|}_{\mathit{p}}\right)=\\ & {\left((\mathrm{2}\mathit{\pi }{)}^m\det \left({\mathbf{R}}^{-\mathrm{1}}\right)\right)}^{-N/\mathrm{2}}\exp \left(-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\sum _{i=\mathrm{1}}^N{\mathit{r}}_i^T\mathbf{R}{\mathit{r}}_i\right).\end{array}\end{array}$$

Maximizing ℒ (or minimizing its negative logarithm), a maximum likelihood estimate of the parameters can be obtained as

$$\begin{array}{}\text{(13)}& {\mathit{p}}_{\mathrm{MLE}}=\arg \underset{\mathit{p}}{min}J\left(\mathit{p}\right),\end{array}$$

where $J\left(\mathit{p}\right)=-\ln \left(\mathcal{L}\right(S{|}_{\mathit{p}})$. The measurement noise covariance matrix R can be estimated under mild hypotheses as $\mathbf{R}={\sum }_{i=\mathrm{1}}^N{\mathit{r}}_i^T{\mathit{r}}_i$, yielding J(p)=det(R), leading to an iteration between a solution at given covariance and a covariance update step (Jategaonkar, 2015). However, in this paper the measurement noise covariance matrix is estimated a priori and therefore assumed to be known. The cost function therefore becomes

$$\begin{array}{}\text{(14)}& J\left(\mathit{p}\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\sum _{i=\mathrm{1}}^N{\mathit{r}}_i^T{\mathbf{R}}^{-\mathrm{1}}{\mathit{r}}_i.\end{array}$$

To ensure reasonable and physically viable solutions, parameters can be forced to stay within predefined upper (subscript ub) and lower (subscript lb) bounds, by adding the corresponding inequality constraints ${\mathit{p}}_{\mathrm{lb}}\le \mathit{p}\le {\mathit{p}}_{\mathrm{ub}}$ to problem (13). As the parameter values and constraints can differ in magnitude, it is a good practice to scale all parameters such that a value of 1 corresponds to the upper bound p_ub and a value of −1 to the lower one p_lb. The optimization problem can finally be solved numerically by a suitable algorithm, such as sequential quadratic programming (SQP) (Nocedal and Wright, 2006).

2.3.2 Identifiability of parameters

The Fisher information matrix $\mathbf{F}\in {\mathbb{R}}^{n\times n}$ is defined as

$$\begin{array}{}\text{(15)}& \mathbf{F}=\sum _{i=\mathrm{1}}^N{\left[{\displaystyle \frac{\partial {\mathit{y}}_i}{\partial \mathit{p}}}\right]}^T{\mathbf{R}}^{-\mathrm{1}}\left[{\displaystyle \frac{\partial {\mathit{y}}_i}{\partial \mathit{p}}}\right]\end{array}$$

and describes the curvature of the likelihood function. It can be shown (Jategaonkar, 2015) that a lower bound (termed Cramér–Rao bound) of the covariance of the estimated parameter is given by

$$\begin{array}{}\text{(16)}& {\mathbf{F}}^{-\mathrm{1}}=\mathbf{P}\le \mathrm{Var}\left({\mathit{p}}_{\mathrm{MLE}}-{\mathit{p}}_{\mathrm{true}}\right),\end{array}$$

where p_true represents the true but unknown parameters. The kth diagonal element of P is a lower bound on the variance of the kth estimated parameter, while the correlation between different parameters is captured by the off-diagonal terms of that same matrix. The correlation coefficient between two parameters i and j is defined as

$$\begin{array}{}\text{(17)}& {\mathrm{\Psi }}_{p_i,p_j}={\displaystyle \frac{P_{i,j}}{\sqrt{P_{i,i}{P}_{j,j}}}},\end{array}$$

where P_i,j denotes the i,jth element (row, column) of P. By analyzing the estimated parameter variance, as well as the correlation between parameters, valuable insight into the well-posedness of the parameter identification problem can be readily obtained.

2.3.3 Problem transformation and untangling using the SVD

When some parameters are highly correlated or have large variance, the problem is ill-posed: it might exhibit sluggish convergence, or no convergence at all, and small changes in the inputs may lead to large changes in the estimates. Such situations are difficult to solve in physical space, because parameters are typically coupled together to some degree through the model.

To untangle the parameters, one may resort to the SVD (Golub and van Loan, 2013). By this approach (Hansen, 1987; Waiboer, 2007; Bottasso et al., 2014a), the original parameters are mapped into a new set of uncorrelated (orthogonal) parameters. Since the new unknowns are uncorrelated, one can set a threshold to their variance by using the Cramér–Rao bound and only retain those in the optimization that are observable within the given data set.

The Fisher matrix F is first factorized as F=M^TM, where $\mathbf{M}\in {\mathbb{R}}^{\mathrm{Nm}\times n}$ is defined as

$$\begin{array}{}\text{(18)}& \mathbf{M}=\left[\begin{array}{l}{\mathbf{R}}^{-\mathrm{1}/\mathrm{2}}\frac{\partial {\mathit{y}}_{\mathrm{1}}}{\partial \mathit{p}}\\ {\mathbf{R}}^{-\mathrm{1}/\mathrm{2}}\frac{\partial {\mathit{y}}_{\mathrm{2}}}{\partial \mathit{p}}\\ \mathrm{\dots }\\ {\mathbf{R}}^{-\mathrm{1}/\mathrm{2}}\frac{\partial {\mathit{y}}_N}{\partial \mathit{p}}\end{array}\right].\end{array}$$

Assuming a larger number of measurements than parameters (Nm>n), matrix M can be decomposed into

$$\begin{array}{}\text{(19)}& \mathbf{M}=\mathbf{U}\mathbf{\Sigma }{\mathbf{V}}^T,\end{array}$$

where $\mathbf{U}\in {\mathbb{R}}^{\mathrm{Nm}\times \mathrm{Nm}}$ and $\mathbf{V}\in {\mathbb{R}}^{n\times n}$ are the matrices of left and right, respectively, singular vectors, while

$$\begin{array}{}\text{(20)}& \mathbf{\Sigma }=\left[\begin{array}{l}\mathbf{S}\\ \mathbf{0}\end{array}\right],\end{array}$$

where $\mathbf{S}\in {\mathbb{R}}^{n\times n}$ is a diagonal matrix, whose entries s_i are the singular values sorted in descending order.

By using Eq. (19) and the factorization of F, the inverse of the Fisher information matrix can be written as

$$\begin{array}{}\text{(21)}& \mathbf{P}={\mathbf{VS}}^{-\mathrm{2}}{\mathbf{V}}^T.\end{array}$$

Note that the columns of the orthogonal matrix V are also the eigenvectors of P and $s_i^{-\mathrm{2}}$ the corresponding eigenvalues. Furthermore, P and F are symmetric and, based on the spectral theorem, diagonalizable.

The physical parameters p can now be transformed into a new set of orthogonal parameters Θ by a rotation performed with the right singular values:

$$\begin{array}{}\text{(22)}& \mathbf{\Theta }={\mathbf{V}}^T\mathit{p}.\end{array}$$

For the transformed set of parameters, the Cramér–Rao bound on the variance of the estimates is the diagonal matrix ${\mathbf{S}}^{-\mathrm{2}}\le \mathrm{Var}({\mathbf{\Theta }}_{\mathrm{MLE}}-{\mathbf{\Theta }}_{\mathrm{true}})$. Therefore, a small singular value s_i corresponds to a large uncertainty in the corresponding orthogonal parameter estimation.

To remove parameters that cannot be estimated with sufficient accuracy, matrix S can be partitioned as

$$\begin{array}{}\text{(23)}& \mathbf{S}=\left[\begin{array}{ll}{\mathbf{S}}_{\mathrm{ID}}& \mathbf{0}\\ \mathbf{0}& {\mathbf{S}}_{\mathrm{NID}}\end{array}\right],\end{array}$$

where S_ID contains the identifiable singular values, i.e. those such that $s_i^{-\mathrm{2}}<{\mathit{\sigma }}_{\mathrm{t}}^{\mathrm{2}}$, σ_t being a threshold on the highest acceptable standard deviation in the estimate. On the other hand, matrix S_NID contains singular values associated with parameters that cannot be identified with sufficient accuracy and are therefore discarded. Accordingly, V is also partitioned as $\mathbf{V}=[{\mathbf{V}}_{\mathrm{ID}},{\mathbf{V}}_{\mathrm{NID}}]$, while the orthogonal parameters are partitioned as $\mathbf{\Theta }=[{\mathbf{\Theta }}_{\mathrm{ID}}^T,{\mathbf{\Theta }}_{\mathrm{NID}}^T{]}^T$. Finally, the physical parameters are expressed in terms of the sole identifiable orthogonal parameters:

$$\begin{array}{}\text{(24)}& \mathit{p}\approx {\mathbf{V}}_{\mathrm{ID}}{\mathbf{\Theta }}_{\mathrm{ID}}.\end{array}$$

Given that the Fisher matrix depends on the values of the parameters p, an iterative procedure should be followed, where the diagonalization of the problem is repeated at each update of the parameter vector.

2.3.4 Identification method with variable measurement weights

In some cases, it may be useful to increase the importance of some measurements in the parameter estimation problem. This can be readily obtained by simply treating an observation with weight w as if it appeared w times in the observation data set (Karampatziakis and Langford, 2011). Cost function (14) then becomes

$$\begin{array}{}\text{(25)}& J\left(\mathit{p}\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\sum _{i=\mathrm{1}}^N{w}_i{\mathit{r}}_i^T{\mathbf{R}}^{-\mathrm{1}}{\mathit{r}}_i,\end{array}$$

where w_i is the relative weight of observation i and ${\sum }_{i=\mathrm{1}}^N{w}_i=N$. Similarly, the Fisher matrix becomes

$$\begin{array}{}\text{(26)}& \mathbf{F}=\sum _{i=\mathrm{1}}^N{w}_i{\left[{\displaystyle \frac{\partial {\mathit{y}}_i}{\partial \mathit{p}}}\right]}^T{\mathbf{R}}^{-\mathrm{1}}\left[{\displaystyle \frac{\partial {\mathit{y}}_i}{\partial \mathit{p}}}\right],\end{array}$$

and its factorization is

$$\begin{array}{}\text{(27)}& \mathbf{M}=\left[\begin{array}{l}\sqrt{w_{\mathrm{1}}}{\mathbf{R}}^{-\mathrm{1}/\mathrm{2}}\frac{\partial {\mathit{y}}_{\mathrm{1}}}{\partial \mathit{p}}\\ \sqrt{w_{\mathrm{2}}}{\mathbf{R}}^{-\mathrm{1}/\mathrm{2}}\frac{\partial {\mathit{y}}_{\mathrm{2}}}{\partial \mathit{p}}\\ \mathrm{\dots }\\ \sqrt{w_N}{\mathbf{R}}^{-\mathrm{1}/\mathrm{2}}\frac{\partial {\mathit{y}}_N}{\partial \mathit{p}}\end{array}\right].\end{array}$$

The remainder of the formulation is not affected by the introduction of weights.

3 Results

The proposed method is first applied in Sect. 3.1 to a wind tunnel experiment with a small cluster of three wind turbines and then in Sect. 3.2 to a real wind farm consisting of 43 wind turbines. The former example aims at a verification of the correctness of the identified augmentations, given the known and controllable conditions of the scaled experiments, whereas the latter is meant to offer a first glimpse of the practical applicability of the new method in the field.

3.1 Wind tunnel verification

Whether identified model corrections are indeed physical or only an artifact of the model–measurement mismatch is difficult to prove in general. From this point of view, wind tunnel experiments provide a unique opportunity to verify the concept proposed in this paper. Indeed, the overall flow within a cluster of turbines can be measured with good accuracy, and the experiments can be repeated in multiple desired operating conditions. The aim of this section is then to show that, even in the presence of multiple possibly overlapping model terms, the correct improvements to a baseline model can be learned from operational data only.

3.1.1 Experimental setup

The experimental setup is composed of a scaled cluster of three G1 wind turbines, each of them equipped with active yaw, pitch, and torque control. The turbines were operated in the boundary layer test section of the wind tunnel of the Politecnico di Milano. Details on the models and the wind tunnel are reported, among other publications, in Campagnolo et al. (2016a, b, c).

The turbines are labeled WT1, WT2, and WT3, starting from the most upstream one and moving downstream. The machines are mounted on a turntable, whose rotation is used to change the wind direction with respect to the wind farm layout. In the nominal configuration, i.e. for a turntable rotation γ_TT=0^∘, the three turbines are aligned with the wind tunnel main axis – and hence with the flow velocity vector. The turbines are installed with a longitudinal spacing of 5 diameters (D), as shown in Fig. 2 with a view looking down towards the wind tunnel floor. As indicated in the figure, positive turntable rotations are clockwise. For γ_TT≠0^∘, the longitudinal distance between the turbines decreases slightly. However, considering that in this work the largest investigated turntable angle was ±11.5^∘, the longitudinal distance varied only between 4.9D and 5D.

Figure 2 Wind farm layout for a null turntable rotation, looking down onto the wind tunnel floor.

A pitot probe was placed at hub height, 3D upstream of the first G1 in the nominal configuration. The probe was therefore not placed on the turntable, and its position remained fixed with respect to the wind tunnel test section. A wind-tunnel-fixed reference frame, used in the following to discuss the results, is also depicted in Fig. 2. Its origin is placed at the turntable center, while the frame x axis is aligned with the wind direction; the y axis points left, looking downstream; and hence Z points vertically up from the floor to complete a right-handed triad.

The yaw angle γ_WTi of the ith wind turbine is positive for a counterclockwise rotation looking down onto the floor, as shown for WT1 in Fig. 2, and null when the rotor disk is orthogonal to X and, therefore, to the nominal wind direction.

Figure 3 shows a photo of the cluster of turbines, looking downstream with WT1 in the foreground. The wind tunnel floor is blue, whereas the turntable is black.

Figure 3 View looking downstream of the cluster of three G1 turbines.

The ambient wind speed V_∞,0 measured by the pitot tube was, for all conducted experiments, between 5.20 and 5.75 m s⁻¹, which corresponds to slightly below-rated conditions. The ambient turbulence intensity was equal to 6.12 %, while the vertical shear was α_vs=0.144.

3.1.2 Model setup

The FLORIS model implementation used in this work is the one available online (Doekemeijer and Storm, 2018). All baseline model parameters are reported in Table 1 and taken from Campagnolo et al. (2019), where they were identified based on wake measurements of a single isolated G1 turbine.

Table 1 Initial FLORIS parameters for the G1 turbine.

Figure 4 shows the G1 power C_P and thrust C_T coefficients as functions of wind speed V. The curves were obtained from dynamic simulations conducted in turbulent inflow, using the same controllers implemented on the scaled models. The C_P and C_T vs. tip speed ratio (TSR) and blade pitch setting curves were obtained with a BEM formulation using experimentally tuned airfoil polars (Bottasso et al., 2014a). As the turbine controller does not consider variations in air density ρ, the coefficients shown in the figure exhibit a slight dependency on this ambient parameter. Within FLORIS, this effect is taken into account by interpolating within the coefficients based on the actual density measured in the wind tunnel during each experiment. For all reported test conditions, air density varied in the range $\mathit{\rho }\in [\mathrm{1.159},\mathrm{1.185}]$ kg m⁻³. The power loss exponent in misaligned conditions was evaluated experimentally to be p_P=2.1741, while for thrust the coefficient was found to be p_T=1.4248.

Figure 4 Power and thrust coefficients vs. wind speed for the G1 turbine.

The ambient wind speed was determined from the pitot tube. It was observed that, by using this value, the power of a free-stream turbine predicted by the FLORIS model was slightly underestimated, most probably due to the sheared flow. To correct for this effect, measurements provided by the pitot tube were scaled by the factor 1.0176, which was computed in order to match simulated and measured power. Furthermore, in the original FLORIS implementation the power of a turbine is computed as $P=\mathrm{1}/\mathrm{2}\mathit{\rho }A{V}_{\mathrm{avg}}^{\mathrm{3}}{C}_{\mathrm{P}}$, where V_avg is the average wind speed at the rotor disk and A the rotor disk area. Here, power was computed by integrating over the rotor disk area, i.e. $P=\mathrm{1}/\mathrm{2}\mathit{\rho }{\int }_A{V}^{\mathrm{3}}{C}_{\mathrm{P}}\mathrm{d}A$, which is probably slightly more accurate even though it involves a minor increase in computational effort.

3.1.3 Ranking of correction terms

To initially assess the role of the various parameters, a ranking analysis was conducted. The parameters were clustered in sets, depending on their role in the model. A first identification was performed using all parameter sets, yielding the presumed best value, denoted J_ref, of the cost function expressed by Eq. (14). The analysis was then repeated multiple times, each time removing one parameter set from the optimization. By looking at the resulting change in the value of the cost function, one may then rank the various parameter sets in order of importance. The analysis is based on a total of 190 experimental observations, as described in greater detail in the following.

Table 2 Definition of the parameters, together with their initial values, lower and upper bounds, and identified values.

All augmentation terms described in Sect. 2.2 were considered, except for the lateral variation in wind direction and the wind-direction-dependent vertical shear, as they are not applicable to the wind tunnel experiments. The nonuniform flow speed was modeled using five nodes located at ${\mathit{c}}_{\mathrm{speed}}\left(Y\right)=[-\mathrm{3},-\mathrm{2},-\mathrm{1},\mathrm{0},\mathrm{1}]$ m (which correspond to approximatively $[-\mathrm{2.7},-\mathrm{1.8},-\mathrm{0.9},\mathrm{0},\mathrm{0.9}]D$) and also indicated in Fig. 2 using × symbols. As only the turbine positions with respect to the flow are modified by rotating the turntable, a wind direction dependency was not included in this correction term. Table 2 reports the initial values and lower and upper bounds – chosen based on an educated guess – for the nonuniform inflow and secondary steering correction terms.

Figure 5 shows the relative increase in the cost function when eliminating one parameter set at a time. The figure clearly indicates that the most important parameters are the ones modeling laterally nonuniform speed and secondary steering. Indeed, this particular wind tunnel, due to its internal configuration and large width, does present a significant nonuniform flow speed, as already discussed by Campagnolo et al. (2019). Likewise, the effect of secondary steering is particularly important and should not be neglected for accurate predictions in misaligned conditions, as already reported in various publications. Based on these results, in the following only nonuniform inflow and secondary steering corrections are considered.

Figure 5 Relative increase in the optimization cost function when eliminating one parameter set at a time.

3.1.4 Results

A total of 451 observations were available, including 11 different turntable positions and thus wind farm layouts, with turbine yaw misalignments ranging from −40 to +40^∘. A total of 190 observations were used to identify the five parameters associated with nonuniform inflow speed and the six associated with secondary steering, whereas the remaining data points were used for model validation. The various tested configurations in terms of turbine misalignments and turntable positions are reported in the figures of Appendix A.

Among all the available measurements gathered at each operating condition, only the steady-state power of the wind turbines was utilized, mimicking what could be done at full scale in the field using SCADA data. The model outputs y (see Eq. 9) are defined as

$$\begin{array}{}\text{(28)}& \mathit{y}={\displaystyle \frac{\mathrm{1}}{P_{\mathrm{ref}}}}\left[\begin{array}{l}{P}_{\mathrm{WT}\mathrm{1}}\\ P_{\mathrm{WT}\mathrm{2}}\\ P_{\mathrm{WT}\mathrm{3}}\end{array}\right],\end{array}$$

where P_WTi is the power of the ith wind turbine and P_ref=37.6W is a reference value used as the scaling factor. Based on experience, a diagonal measurement noise covariance matrix R with all three terms equal to σ²=0.025² was specified.

The threshold of the highest acceptable standard variance ${\mathit{\sigma }}_{\mathrm{t}}^{\mathrm{2}}$ for the orthogonal parameters was set to 0.01. As the parameters are scaled within a range of $[-\mathrm{1},\mathrm{1}]$, the threshold corresponds to a relative variance of 2 %. Wind-aligned operating conditions (i.e. ${\mathit{\gamma }}_{\mathrm{WT}\mathrm{1}}={\mathit{\gamma }}_{\mathrm{WT}\mathrm{2}}={\mathit{\gamma }}_{\mathrm{WT}\mathrm{3}}=\mathrm{0}$^∘) were weighted with a factor of 2, to increase their importance in the parameter estimation process.

Figure 6 Variance of the orthogonal parameters before (a) and after (b) the first iteration. The identifiable orthogonal parameters are shown in red, whereas all others are shown in blue.

The constrained optimization problem (13) was solved in MATLAB using the fmincon function with the interior-point algorithm (Mathworks, 2019). As the baseline model with its initial nominal values (p=p_init) is far away from the optimal solution, a first optimization was performed including only the inflow correction. Afterwards, three iterations were conducted including all 11 parameters. At each iteration, a total of eight orthogonal parameters could be identified within the specified variance threshold. The method converged very quickly, as the identified parameters and the residual did not change significantly after the first iteration. Figure 6a shows the initial variance of all 11 orthogonal parameters, and panel (b) shows the variance computed after the first iteration. The horizontal black line indicates the threshold ${\mathit{\sigma }}_{\mathrm{t}}^{\mathrm{2}}$.

Table 3 Transformation matrix V^T after the first iteration. Each row corresponds to a different orthogonal parameter.

Interestingly, the 11th orthogonal parameter seems to have a very low observability. Table 3 shows the transformation matrix V^T that links the physical parameters to the orthogonal ones (Θ=V^Tp; see Eq. 22). The 11th orthogonal parameter is almost entirely associated with p_speed,5, which corresponds to the inflow speed augmentation node at position Y=1 m. Indeed, the location of this node is such that it has only a very marginal effect on the turbine outputs and, hence, a very low observability, as shown later in Fig. 7. The transformation matrix reported in Table 3 also shows that the other two orthogonal parameters with low observability (9 and 10) represent secondary steering modes, mainly associated with the second Gaussian function of the correction term.

Figure 7 Identified nonuniform inflow speed augmentation term (solid line) and associated standard deviation (whiskers). Hot-wire measurements at different heights above the floor are shown in thin solid lines. The upstream turbine (WT1) y position for all investigated turntable rotations is shown by × markers placed along the lower edge of the figure.

Table 4 Correlation coefficients Ψ after the first iteration.

Table 4 presents the correlation matrix Ψ (see Eq. 17) and shows a clear and to be expected dependency among neighboring inflow parameters. Among the secondary steering parameters, strong but less obvious correlations are present, which suggest that a simplification of the assumed correction term might be possible.

Figure 7 shows the identified inflow augmentation function. In the picture, whiskers indicate the parameter uncertainty σ_i, computed based on the Cramér–Rao lower error bound as $\mathit{\sigma }=\sqrt{\mathrm{diag}\left(\mathbf{P}\right)}$ (see Eq. 16). The same figure also reports measurements obtained with hot-wire probes in the empty wind tunnel at three different heights above the floor. These measurements, and especially the ones at hub height, are in good agreement with the estimates provided by the proposed method. The figure also reports (with × symbols) the lateral position of the upstream turbine for the investigated turntable rotations. Noting that all points are shifted to the left helps explain why the parameter associated with the inflow node at Y=1 m has a very low – but still finite – observability.

The identified secondary steering augmentation term is visualized in Fig. 8. The plot shows the wind direction change ΔΓ as a function of the distance $\stackrel{\mathrm{̃}}{y}$ to the wake centerline for a turbine misalignment of 20^∘. The gray shaded area shows the uncertainty band $p_{\mathrm{opt},i}\pm {\mathit{\sigma }}_i$. Consistently with the findings of Wang et al. (2018), the maximum change in wind direction is found at approximatively 0.3D on the leeward side of a deflected wake. The maximum magnitude of secondary steering in this operating condition is 1.9^∘, which is again comparable to the results of Wang et al. (2018).

Figure 8 Identified wind direction change ΔΓ due to secondary steering as a function of distance $\stackrel{\mathrm{̃}}{y}$ to the wake centerline for a turbine misalignment of 20^∘. The gray shaded area shows the uncertainty band.

Figure 9 Wake profiles 5D behind WT2 for various combinations of turbine yaw misalignment. Experimental values are indicated by the × symbols. Each subplot is accompanied by two flow visualizations based on the FLORIS model and its augmented version.

The validity of the augmentation terms, identified as explained, was assessed by comparing the results of the simulation model with experimental wake measurements from a different test campaign. The setup was identical to the one considered here, except for the fact that only the first two upstream wind turbines were installed in the wind tunnel. At the downstream distance where the third wind turbine should have been installed, flow velocity measurements were obtained at turbine hub height using hot-wire probes. Figure 9 shows wake profiles for the turntable position γ_TT=0^∘ for various combinations of turbine yaw misalignments, as indicated by the subplot titles. Each subplot is accompanied by two flow visualizations, one based on the baseline FLORIS model and the other on its augmented version. The figures also include the points at which the flow was measured with the probes.

In the left subplots, the improvements of the augmented model with respect to the baseline FLORIS are exclusively due to the inflow correction, as the upstream turbine is aligned with the flow and therefore there are no secondary steering effects. In the right subplots, the upstream turbine is misaligned (γ_WT1=30^∘) and secondary steering effects are present. Taking into account that model augmentation was obtained exclusively by turbine power measurements, the improved matching of the wake profiles is remarkable. Still, even with the extra correction terms some small model mismatches are present; these might be caused by the wake combination model, which was not augmented in this study.

The turbine power coefficients are computed as

$$\begin{array}{}\text{(29)}& C_{\mathrm{P},i}={\displaystyle \frac{P_{\mathrm{WT}i}}{\mathrm{0.5}\mathit{\rho }A{V}_{\mathrm{\infty }}{\left(Y_{\mathrm{WT}i},z_{\mathrm{h}},\mathrm{0}\right)}^{\mathrm{3}}}},\end{array}$$

where V_∞ is the augmented inflow function given by Eq. (2), evaluated at the respective turbine position Y_WTi and hub height z_h. A detailed overview of the results is offered by the figures of Appendix A, which report the power outputs and the model errors for all wind farm configurations. For readability, here a more synthetic overview of the results is presented, by condensing the information contained in Figs. A1, A2, and A3 in the probability density plots of Fig. 10. This figure shows the results for the baseline FLORIS model using a black dashed line, for the 11-parameter augmented model (i.e. including only nonuniform inflow speed and secondary steering corrections) using a red solid line, and for the 27-parameter augmented model (i.e. including all additional augmentation terms presented earlier) using a red dotted line. The root-mean-squared errors ϵ_RMS are shown in the respective legends.

Figure 10 Error distributions for each turbine for all tested configurations, for the baseline FLORIS model (black dashed line), the 11-parameter augmented model (red solid line), and the 27-parameter augmented model (red dotted line).

Note that the FLORIS error distribution shows two peaks for WT1 and WT3, indicating the presence of two uncorrelated errors. The 11-parameter model removes these peaks, even though a smaller pair of peaks remains for WT2 and WT3, indicating additional errors that only the 27-parameter augmented model is able to capture.

Here again the trend is clear: the addition of nonuniform speed and secondary steering substantially increases the accuracy of the baseline model, with additional small – but not insignificant – gains offered by the additional correction terms. Finally, there is still room for improvement, possibly through extra correction terms not yet explored.

Table 5 Turbine specifications.

3.2 Field application

In this section the model augmentation and identification method is applied to a full-scale wind farm, to test its applicability and usability in a realistic scenario. In such conditions, it is often difficult to assess whether the identified model corrections are indeed physical or not, due to a lack of knowledge of the actual ground truth. To deal with this problem, the classical approach of splitting the data set was used here: first, a relatively small subset of measurements is used for model and error identification; then, the rest of the data set is used for a verification of the generality of the identified model and of its improved performance with respect to the baseline one.

3.2.1 Wind farm and data preprocessing

The onshore wind farm is situated close to Sedini, on the Italian island of Sardinia, and it consists of 43 GE1.5s and GE1.5sle wind turbines, as specified in Table 5.

Figure 11 The 3D view of the Sedini wind farm with terrain elevation, as seen from Γ=260^∘.

The wind farm is located at a rather complex site, as shown in Fig. 11. Blue turbines are of the type GE1.5sle and black and red turbines are of the type GE1.5s, the latter being used as sensing turbines as explained later. Figure 12 shows a top view of the wind farm, including the turbine identifiers.

Figure 12 Top view of the Sedini wind farm with turbine identifiers. The gray arrows indicate the x and y axes for an ambient wind direction Γ=260^∘.

Historical 10 min SCADA data were made available for this research for a period of 24 months, throughout the years 2015 and 2016. The recorded turbine yaw orientations exhibit sudden jumps and long-term drifts. An ad hoc algorithm was developed for detecting and correcting these data issues. On average, for each turbine 45 % of the data points were missing, and 23 % were discarded because of low power output (<5 kW) or rotor speed (<1 rpm). As a result, about 33 700 data points were available for each turbine. Regarding the missing data points, it is unknown whether the turbines were operating or just not reporting. To avoid eliminating a large fraction of the data set, it was assumed that the turbines were indeed operational and thus shedding wakes. This way, even if recordings of one or more turbines were missing at a specific time instance, the data points of the other turbines could still be used.

Figure 13 Scaled number of measurement data points (10 min mean) within each speed and direction bin.

As no direct measurements of ambient conditions were available, the method described by Schreiber et al. (2018) was used to identify ambient wind speed and direction. The procedure works as follows. First, the ambient wind direction is estimated from turbine yaw orientations. Second, the ambient wind speed is estimated from the rotor effective wind speed of the free-stream turbines, computed from the turbine power curve below rated wind speed. For this purpose, the three sensing turbines A5-24, A5-25, and A5-26 indicated in red in Fig. 12 were used, checking that they were unwaked by using the flow model; the average of these speeds was attributed to the location of turbine A5-25. This way, 5667 ambient wind conditions could be processed for a range of wind directions $\mathrm{\Gamma }\in [\mathrm{184}{}^{\circ },\mathrm{320}{}^{\circ }]$. Based on the ambient wind conditions, the data of all turbines were aggregated in two-dimensional bins: ambient wind speed (bin width of 2 m s⁻¹) and ambient wind direction (bin width of 5^∘). Figure 13 shows the scaled number of measurements in each bin between 6 and 12 m s⁻¹.

3.2.2 Model setup

Here again the FLORIS implementation was based on the version available online (Doekemeijer and Storm, 2019). The initial values of both the wake and turbulence model parameters were set according to Bastankhah and Porté-Agel (2016) for $({\mathit{\alpha }}^{\ast },{\mathit{\beta }}^{\ast })$, Crespo and Hernández (1996) for $({\mathrm{TI}}_{\mathrm{a}}^{\ast },{\mathrm{TI}}_{\mathrm{b}}^{\ast },{\mathrm{TI}}_{\mathrm{c}}^{\ast },{\mathrm{TI}}_{\mathrm{d}}^{\ast })$, Niayifar and Porté-Agel (2015) for $(k_{\mathrm{a}}^{\ast },k_{\mathrm{b}}^{\ast })$, and Gebraad et al. (2014) for $(a_{\mathrm{d}}^{\ast },b_{\mathrm{d}}^{\ast })$, as reported in Table 6.

Table 6 Initial FLORIS parameters for the Sedini wind farm.

The required turbine power and thrust versus wind speed curves were provided by the turbine manufacturer. The vertical shear exponent of the inflow was set to α_vs=0.143 and the turbulence intensity to 14 %, which represent annual average values measured at 65 m of height by an on-site met mast. Air density was set to the constant value ρ=1.177 kg m⁻³.

Figure 14 Correlation between power output and hub height with respect to SL. (a) Power (× symbols and left y axis) and rotor height above SL (∘ symbols and right y axis) vs. lateral turbine position for a wind direction Γ=240^∘. (b) Power vs. rotor height above SL for $\mathrm{\Gamma }\in [\mathrm{220}{}^{\circ },\mathrm{275}{}^{\circ }]$ and $V_{\mathrm{\infty }}\in [\mathrm{8},\mathrm{10}]$ m s⁻¹. All conditions are free stream and all turbines of type GE1.5s.

The different turbine foundation heights were accounted for by accordingly increasing the tower heights, using the lowest foundation height as reference (turbine A1-02). Indeed, power measurements of the upstream turbines show a correlation with the actual turbine hub height with respect to sea level (SL), as shown in Fig. 14. As indicated by the only approximate correlation shown by the figure, it is clear that such simple correction might not provide satisfactory results for all wind directions and all turbines, because complex orthographic flow effects might also play a role. Nonetheless, this approximate correction seems to be a step in the right direction. In addition, some of these effects may be corrected by the lateral nonuniformity terms added to the augmented model. The reference height of the sheared inflow z_h (see Eq. 2) was set to the hub height of the sensing turbine A5-25.

3.2.3 Ranking of correction terms

As for the wind tunnel experiments, here again a first analysis was aimed at ranking the various correction terms. However, since the turbines were operated with a conventional wind-aligned strategy, secondary steering corrections were neglected. The ranking is based on data points in the range $V\in [\mathrm{8},\mathrm{10}]$ m s⁻¹, as described in greater detail in the following.

Figure 15 shows the relative increase in the cost function after optimization eliminating one set of parameters at a time. The results clearly indicate that the nonuniform wind farm inflow speed p_speed is the most important correction. In fact, this was to be expected, given that the Sedini wind farm is located at a rather complex site. Results also indicate a non-negligible effect of the wake deflection parameters for non-misaligned operation (a_d,b_d).

Figure 15 Relative increase in the optimization cost function for the Sedini wind farm when eliminating one parameter set at a time.

On the other hand, the additional model augmentation parameters $({\mathit{p}}_{\mathrm{TI}},{\mathit{p}}_{\mathrm{winddir}},{\mathit{p}}_{\mathrm{acc}},{\mathit{p}}_{\mathrm{shear}})$ do not seem to contribute to a significant extent. Note also the slight retuning of parameters $(\mathit{\alpha },\mathit{\beta },k_{\mathrm{a}},k_{\mathrm{b}})$ and $({\mathrm{TI}}_{\mathrm{a}},{\mathrm{TI}}_{\mathrm{b}},{\mathrm{TI}}_{\mathrm{c}},{\mathrm{TI}}_{\mathrm{d}})$, which can be explained with the fact that their initial values were taken from the literature and therefore apply to different turbine types and sites.

Given these results, the rest of the analysis is based only on the subset of parameters p_inflow, $(p_{a_{\mathrm{d}}},p_{b_{\mathrm{d}}})$, (p_α,p_β), $(p_{k_{\mathrm{a}}},p_{k_{\mathrm{b}}})$, and $(p_{{\mathrm{TI}}_{\mathrm{a}}},p_{{\mathrm{TI}}_{\mathrm{b}}},p_{{\mathrm{TI}}_{\mathrm{c}}},p_{{\mathrm{TI}}_{\mathrm{d}}})$. The augmentation term for nonuniform inflow speed is modeled using five nodes along the lateral position Y located at $[-\mathrm{2000};-\mathrm{1000};\mathrm{0};\mathrm{1000};\mathrm{2000}]$ m (which is approximatively $[-\mathrm{28};-\mathrm{14};\mathrm{0};\mathrm{14};\mathrm{28}]D_{\mathrm{GE}\mathrm{1.5}\mathrm{s}}$) and six nodes in wind direction Γ at $[\mathrm{180};\mathrm{210};\mathrm{140};\mathrm{270};\mathrm{300};\mathrm{330}]$^∘, resulting in 30 nodes. The Y-coordinate axis is orthogonal to the wind direction and its origin Y=0 m is located at the position of wind turbine A5-25, as shown in Fig. 12.

Table 7 Definition of the parameters, together with their lower and upper bounds, and initial and identified values. Bold italic numbers indicate vectors containing that number repeated as many times as the vector length.

The definitions of the correction parameter, together with their bounds and converged values, are reported in Table 7. Note that all parameters were set to zero at the beginning of the identification process.

3.2.4 Results

To identify the 40 parameters of Table 7, only aggregated mean power measurements for wind speeds $V\in [\mathrm{8},\mathrm{10}]$ m s⁻¹ were used. In addition, only one-third of all wind direction bins were employed,

Figure 16 Identified inflow augmentation parameters (a) and their uncertainties (b). Nodal points are indicated by the circle markers.

The model outputs y (see Eq. 9) were defined as

$$\begin{array}{}\text{(30)}& \mathit{y}={\displaystyle \frac{\mathrm{1}}{P_{\mathrm{ref}}}}\left[\begin{array}{l}{P}_{\mathrm{WT}\mathrm{1}}\\ \mathrm{\dots }\\ P_{\mathrm{WT}\mathrm{43}}\end{array}\right],\end{array}$$

where P_WTi is the power of wind turbine i and P_ref=1.11 MW a reference wind turbine value used as a scaling factor. A diagonal measurement noise covariance matrix R was used, with all diagonal terms equal to σ²=0.01². The threshold of the highest acceptable variance in the orthogonal parameter estimate was set to ${\mathit{\sigma }}_{\mathrm{t}}^{\mathrm{2}}=\mathrm{0.01}$, which corresponds to a relative variance of 2 %. The relative weight of each observation was set proportional to the number of measurement points within the respective bin. In a first iteration, 29 orthogonal parameters could be identified. In the second and third iterations only 23 and 25 orthogonal parameters fell below the threshold, although results changed only marginally after the first iteration.

The identified optimal parameter values p_opt,i are included in Table 7 and, for the inflow augmentation, are also reported in Fig. 16. The plot shows, according to the color map, the inflow augmentation function values $f_{\mathrm{augm},\mathrm{speed}}(Y,\mathrm{\Gamma },{\mathit{c}}_{\mathrm{speed}},{\mathit{p}}_{\mathrm{speed}})$ in panel (a). Each nodal point is indicated by a circle marker. The figure shows that significant variations in the inflow speed have been detected: for example, considering Γ=270^∘, the inflow speed at $Y=+\mathrm{1000}$ m (approximately at the location of wind turbines A3-19, A3-20, and A3-21) is 3.5 % smaller than the one measured at the reference turbines A5-24, A5-25, and A5-26. For the same wind direction, the speed at $Y=-\mathrm{1000}$ m (approximately located at the wind turbines A4-36, A4-37, and A4-38) is 4.8 % larger. These variations are expected to be mainly caused by terrain effects. Panel (b) of Fig. 16 shows the parameter uncertainty (Cramér–Rao bounds). The parameter at the nodal point ($Y=-\mathrm{2000}$ m; Γ=330^∘) is completely unobservable, because it lies far outside of the wind farm perimeter (see Fig. 12). Some of the outer nodal points at $Y=\pm \mathrm{2000}$ m do show significantly increased uncertainties. However, the corresponding augmentation parameters (panel a) are approximatively zero.

Figure 17 Power coefficient of each individual wind turbine, as indicated by the subplot title, as a function of wind direction Γ for wind speeds $V\in [\mathrm{8},\mathrm{10}]$ m s⁻¹. The gray shaded area indicates the standard deviation within the binned measurements. The number of measurements within each bin is reported in Fig. 13.

Figure 17 shows the power coefficient of each individual wind turbine, as indicated by the subplot title, as a function of wind direction. The power coefficient is computed as $C_{\mathrm{P}}=P/\left(\mathrm{0.5}\mathit{\rho }A{V}^{\mathrm{3}}\right)$, where ρ=1.177 kg m⁻³ is the constant air density, $A=\mathit{\pi }(\mathrm{70.5}/\mathrm{2}{)}^{\mathrm{2}}$ m² a reference rotor area, and V the corresponding estimated ambient wind speed. Blue crosses indicate SCADA data points, with the ones used for identification circled. The gray shaded area indicates the standard deviation within the binned measurements. The FLORIS (non-augmented) power estimates are shown by black dashed lines, whereas the augmented model results are shown using red solid lines.

Even though the baseline FLORIS power estimates already exhibit a reasonable correlation with the measurements for many turbines and wind directions, a significant improvement is achieved by the augmented model. Note that for Γ<210^∘ and Γ>300^∘ the number of measurement points within each bin is reduced (see Fig. 13), limiting the measurement quality and trustworthiness. More specifically, the augmented model shows improvements in the modeling of the free-stream turbine power, due to the effects of the wind farm inflow augmentation terms. Furthermore, the predictions of the wake-induced power deficits are corrected, improving in many cases the deficit depth as well as the deficit location in terms of wind direction.

Figure 18 Error probability density functions for different wind speed ranges.

The same results of Fig. 17 are also presented in a more synthetic form in terms of error probability densities in Fig. 18, where the error is defined as $\mathit{ϵ}=C_{\mathrm{P},\mathrm{Meas}.}-C_{\mathrm{P},\mathrm{FLORIS}/\mathrm{Augm}.}$. Each subplot shows the results for a different wind speed range. Note that the modeling error is also reduced for wind speed ranges that have not been used for model identification (i.e. $V\in [\mathrm{6},\mathrm{8}]$ m s⁻¹ and $V\in [\mathrm{10},\mathrm{12}]$ m s⁻¹). The overall root-mean-squared error is reported within the legend, showing error reductions of 14 %, 22 %, and 19 %, highlighting the generality of the identified model and augmentation parameters.

4 Conclusions

This paper has presented a new method to calibrate and augment parametric wind farm models. The proposed approach builds on the vast body of knowledge and experience embedded in available reduced wind farm flow models. However, recognizing that any such model will always have only a limited prediction accuracy, the present approach augments a baseline model with extra ad hoc terms designed to correct some of its presumed specific deficiencies. These additional elements of the model are then learned from operational data. Optionally, the baseline model parameters can also be tuned within a single integrated process. By design, the method has been exclusively based here on SCADA power measurements; therefore, it is readily applicable to most operational wind farms, whenever such data are available. However, the concept of model augmentation is very general and could clearly also be used with other measurements.

To limit the number of free parameters and to overcome the fact that the identification problem can be over-parameterized and hence ill-posed, a parameter transformation into an orthogonal space has been used. Thereby, only parameters that are sufficiently visible within a given data set enter into the identification process.

The method was first applied to a large data set obtained with scaled wind turbines operating in a boundary layer wind tunnel. Thereby, it was shown that a correct learning of the extra modeling terms is achieved. These conclusions are made possible by the fact that, in this case, the flow and wake characteristics are known with good accuracy. Next, the method was tested on a real wind farm, in a realistic and highly complex situation.

Based on the results shown here, the following conclusions can be drawn.

• Within the wind tunnel environment, a correct learning of nonuniform wind farm inflow speed and of secondary steering effects has been achieved. In particular, the latter shows a good match with detailed wake measurements in wind-misaligned conditions. It is remarkable, and very promising, that such detailed features of the solution could be inferred purely from operational power data, even when starting from a baseline model that does not at all consider secondary steering.

• The application to field data has shown that, as expected for the complex-terrain site analyzed here, orographic effects play a driving role. A marked model improvement could be observed, even in conditions where the model was used for extrapolating outside of the training conditions. It is worth noting that, in many practical onshore applications, orographic effects will be present, and the fact that one can learn them from simple and readily available operational data is very encouraging. Again, it should be explicitly pointed out that the baseline model did not include any orographic corrections.

• It has been shown that model tuning and the learning of extra correction terms can be achieved simultaneously. This reduces the risk of adapting the baseline parameters beyond their reasonable limits, driven by unmodeled physics.

• Although the augmented models show a much improved accuracy with respect to the baseline, some model mismatch still remains. Although these remaining errors may often be caused by issues in the data rather than in the model, additional improvements are thought to be possible.

Future work will apply the proposed method to other wind farms, to increase confidence in the obtained results. From longer and richer data sets, possibly in conjunction with meteorological reanalyses, it is presumed that yearly and seasonal variations could be observed. The integration of CFD analyses can be used to support and confirm the identification of orographic effects. Attention should also be paid to improved and additional forms of model corrections, including wake overlap models. Finally, it is worth pointing out again that an improved knowledge of the flow within a wind farm finds applicability in a potentially large range of digitally driven applications, including wind farm control, lifetime estimation, power forecasting, predictive maintenance, and others. Therefore, it is expected that methods for high-accuracy flow predictions in wind farms will be the subject of significant future research efforts.

Appendix A: Extended wind tunnel results

Figures A1, A2, and A3 report the power outputs of WT1, WT2, and WT3, respectively, for all tested configurations. In each figure, clusters of three subplots represent a unique turntable position, as indicated by the title and the wind farm layout sketch therein. The left part of each subplot shows the turbine power coefficient C_P,WTi as a function of γ_WT1 (x axis) and γ_WT2 (y axis). All measured configurations are indicated by a small cross symbol, whereas the measurements used for parameter identification are circled. The central part of each subplot shows the FLORIS model error ${\mathit{ϵ}}_{\mathrm{FLORIS}}=C_{\mathrm{P},\mathrm{Meas}.}-C_{\mathrm{P},\mathrm{FLORIS}}$, including an annotation of the root-mean-squared error ϵ_RMS. Similarly, the right part of each subplot shows the augmented model error ϵ_Augm..

For the first upstream wind turbine, WT1, the baseline FLORIS shows significant errors depending on the turntable position. For γ_TT<0^∘ the model underpredicts turbine power because of the lack of uniformity of the flow, as also shown in Fig. 7. The opposite behavior can be seen for γ_TT>0^∘. The augmented model however shows significant improvements, which are due to the inflow correction. Still, some underprediction for ${\mathit{\gamma }}_{\mathrm{TT}}=-\mathrm{11.5}$^∘ is present, which is probably caused by an excessively small number of parameters in the inflow augmentation function and/or by the third wind turbine power measurements, which are also strongly affected by lateral inflow variations.

The power of WT2, shown in Fig. A2, is only weakly affected or improved by the model corrections. In fact, in all investigated conditions, the second turbine lateral position remains almost constant, such that the inflow correction does not have a significant direct effect. However, secondary steering only slightly changes the inflow direction at WT2; for example, as shown in Fig. 8, a 20^∘ misalignment of WT1 changes the wind direction by about 1.9^∘. This leads to small misalignments and thus only very small changes in power output considering the cosine law. In addition, secondary steering also leads to a slight lateral deflection of the nonuniform inflow.

The power of WT3, reported in Fig. A3, shows significant improvements when using the augmentation terms. For example, for γ_TT>0^∘ the baseline model underpredicts the real flow velocities – and hence the power output – at WT3, an error that is corrected by the augmented model. In addition, for $\mid {\mathit{\gamma }}_{\mathrm{WT}\mathrm{1}}\mid >\mathrm{0}$, secondary steering augmentation affects the deflection of the second turbine wake (as shown in Fig. 8), leading to further improvements.

Figure A1 Wind turbine WT1. Each cluster of three subplots represents a unique turntable position, as indicated by the title and the wind farm layout sketch. Left subplot: turbine power coefficient C_P,WT1 as a function of γ_WT1 (x axis) and γ_WT2 (y axis). Middle subplot: FLORIS model error. Right subplot: augmented model error. Cross symbols: all measured configurations. Circles: conditions used for parameter identification.

Figure A2 Wind turbine WT2. Each cluster of three subplots represents a unique turntable position, as indicated by the title and the wind farm layout sketch. Left subplot: turbine power coefficient C_P,WT2 as a function of γ_WT1 (x axis) and γ_WT2 (y axis). Middle subplot: FLORIS model error. Right subplot: augmented model error. Cross symbols: all measured configurations. Circles: conditions used for parameter identification.

Figure A3 Wind turbine WT3. Each cluster of three subplots represents a unique turntable position, as indicated by the title and the wind farm layout sketch. Left subplot: turbine power coefficient C_P,WT3 as a function of γ_WT1 (x axis) and γ_WT2 (y axis). Middle subplot: FLORIS model error. Right subplot: augmented model error. Cross symbols: all measured configurations. Circles: conditions used for parameter identification.