Towards real-time optimal control of wind farms using large-eddy simulations

Large-eddy simulations (LESs) are commonly considered too slow to serve as a practical wind farm control model. Using coarser grid resolutions, this study examines the feasibility of LES for real-time, receding-horizon control to optimize the overall energy extraction in wind farms. By varying the receding-horizon parameters (i.e. the optimization horizon and control update time) and spatiotemporal resolution of the LES control models, we investigate the trade-off between computational speed and controller performance. The methodology is validated on the TotalControl Reference Wind Power Plant using a fine-grid LES model as a wind farm emulator. Analysis of the resulting power gains reveals that the performance of the controllers is primarily determined by the receding-horizon parameters, whereas the grid resolution has minor impact on the overall power extraction. By leveraging these insights, we achieve near-parity between our LES-based controller and real-time computational speed, while still maintaining competitive power gains up to 40 %.

1 Introduction

Turbine–wake interactions can significantly impact the efficiency of energy extraction when many turbines are clustered together in large-scale wind farms. Standard control strategies do not account for the wake interactions but maximize performance at the turbine level, resulting in significant power deficits and increased loading in downstream regions. In the last decade, much research has been done into dynamic receding-horizon optimal control strategies to mitigate these effects, both through axial induction control and wake redirection (see Meyers et al., 2022 for a recent review). More recently, Goit and Meyers (2015), Goit et al. (2016) and Munters and Meyers (2017, 2018b) have developed an optimal control framework for wind farm power maximization based on high-fidelity large-eddy simulations (LESs) of the wind farm boundary layer. In their latest work combining overinduction and wake redirection control, Munters and Meyers (2018b) report energy gains of up to 34 % for an aligned 4×4 wind farm. Despite its accuracy, LES is usually deemed impractical for real-time applications due to its high computational cost. In that sense, the aforementioned optimal control studies were intended as benchmarking studies, aimed at exploring the potential of LES for power optimization in wind farms.

To achieve real-time optimal control, one option is to resort to models that are less computationally intensive (compared to 3D LES). For instance, in some dynamic flow models, the vertical dimension of the flow is either disregarded or approximated to cope with the computational cost of 3D wake dynamics (see, e.g., Soleimanzadeh et al., 2014; Rott et al., 2017; Boersma et al., 2018). This results in a 2D LES-like model suitable for online wind farm control. Instead of LES, more simplified formulations of the governing equations, such as the 2D dynamic wake meandering model (Jonkman et al., 2017) or the Reynolds-averaged Navier–Stokes equations (Iungo et al., 2015), can be employed to accelerate the computations. In Shapiro et al. (2017), an even simpler 1D wake model is proposed for closed-loop receding-horizon control. However, these expeditious engineering models potentially lack the necessary physical intricacies inherent to 3D LES and may not capture the actual turbulent wake dynamics.

The present paper is a first investigation on the feasibility of using LES as a real-time plant model for receding-horizon wind farm control. To overcome the challenge of computational speed, this study aims to leverage the insights from the earlier work of Bauweraerts and Meyers (2019). In the context of turbulent flow forecasting in the atmospheric boundary layer, they demonstrated that prediction errors only slowly increase when coarsening the grid. By resorting to coarser-grid formulations and incorporating an efficient spatial parallelization, they were able to reduce LES wall times up to a factor of 300 compared to simulated time. We envisage a similar approach, but focusing on LES-based receding-horizon control. By varying the spatiotemporal grid resolution of the LES plant model, we investigate the trade-off between computational speed and performance of the controller. In view of the latter, we also study the influence of the parameters of the receding-horizon framework, i.e. the optimization horizon, the control update time and number of optimization iterations. To take into account the computational times for the optimization and allow for a practical, real-time control action, the framework from Munters and Meyers (2017) is applied in a time-decoupled fashion where the control signals are computed ahead of time (based on a prediction of the future state; see, e.g., Grüne and Pannek, 2017). The proposed methodology is demonstrated on the TotalControl Reference Wind Power Plant (TCRWP) (see, e.g., Andersen et al., 2018), combining yaw and induction control strategies.

In the context of wind farm modeling, the coarse LES models envisioned in this work may not fully capture all the relevant dynamics in the turbulent wakes. In general, finer grids are required to accurately represent secondary flow features such as (the breakdown of) tip vortices and helical structures. For example, in the dynamic induction control (DIC) strategy proposed by Munters and Meyers (2018a), one of the main mechanism to enhance wake mixing is the periodic shedding of vortex rings from front-row turbines synchronized to the turbulent inflow. These vortex rings cannot be accurately represented on a coarse resolution. More novel approaches, such as the helix approach, also require finer grids to represent the helical structures in the wakes, especially near the turbine rotor (Frederik et al., 2020). However, for other phenomena that are mostly triggered by large-scale motions in the flow, coarser grids may suffice. Examples of the latter include the overall deflection and gross behavior of the wakes and potentially also wake meandering triggered by the time-varying inflow conditions or incited by dynamically yawing the turbines (Meyers et al., 2022).

The paper is organized as follows. In Sect. 2, we first summarize some important aspect of the LES modeling and introduce the time-decoupled receding-horizon framework, as well as how this framework is adapted to incorporate the coarsening strategy from Bauweraerts and Meyers (2019). Next, the TCRWP test case and simulation setup are discussed in Sect. 3. In Sect. 4, we present the results and discuss the influence of the grid resolution and receding-horizon parameters on the performance and computational time of the LES-based controllers. Section 5 then leverages these insights to design a competitive controller (in terms of power gains) as close to real time as possible. Section 6 concludes the paper and summarizes the main contributions.

2 Methodology

We first discuss the time-decoupled receding-horizon optimal control framework in Sect. 2.1. Next, we describe the wind farm optimization problem in Sect. 2.2 and highlight some aspects of turbine modeling in Sect. 2.3. The optimization method and gradient computation are discussed in Sect. 2.4. Finally, the grid-coarsening strategy is elaborated in Sect. 2.5.

2.1 Receding-horizon optimal control

Previous LES-based wind farm control studies (see, e.g., Goit and Meyers, 2015; Goit et al., 2016; Munters and Meyers, 2017, 2018b) have adopted a simplified model-predictive control (MPC) framework where the optimization is performed on an accurate control model assuming perfect knowledge of the system state; see Fig. 1a. This control loop was applied in a receding-horizon fashion (see Fig. 1b, where time is divided into overlapping windows of length T_A. In each consecutive window, the controls φ were optimized over a prediction horizon T through a series of fine-grid forward and adjoint LES simulations and then applied to the wind farm over a control update time T_A (with T_A<T). These controllers, however, are not realizable in practice, because full state information is not available and the fine-grid LES control model is unfeasible due to its excessive computational cost. Moreover, in their framework, the computational time to solve the optimization problems is ignored, since controls are first computed and then applied to the same time interval, which is not possible in real time.

Figure 1 MPC approach from earlier studies (Goit and Meyers, 2015; Goit et al., 2016; Munters and Meyers, 2017, 2018b). (a) Control loop. (b) Receding-horizon approach. In every window i, the optimization stage is represented by a series of forward and adjoint simulations, resulting in controls φ_i that are applied to the farm over a control update time T_A. Figures adapted from Munters and Meyers (2018b).

To account for computational time and allow real-time control, we propose a time-decoupled approach where controls are computed with a predefined offset corresponding to the control update time T_A from previous studies. The time-decoupled MPC loop, including state estimation, is shown in Fig. 2a (Grüne and Pannek, 2017). For $t\in [t^k,t^k+T_A]$, the estimator yields an estimate $\widehat{\mathit{q}}\left(t^k\right)$ of the instantaneous flow at t^k=kT_A. The estimate is computed using lidar or SCADA measurements h(q(t)) that were collected over the past estimation window $[t^k-T_{\mathrm{SE}},t^k]$ and stored in the estimation buffer (with h(⋅) the measurement function), assuming moving horizon estimation with horizon T_SE (note that the MPC loop may differ on the estimation side when using other approaches such as a Kalman filter). Next, the predictor uses a flow model to propagate the estimate, yielding a prediction $\widehat{\mathit{q}}(t^k+T_A)$ of the future wind farm state. This allows the optimizer to compute the optimal controls ahead of time, resulting in φ_k+1(t) for $t\in [\mathrm{0},T]$ (with T the optimization horizon). This set of controls is then stored in the actuator buffer. In the next window ($t\in [t^k+T_A,t^k+\mathrm{2}{T}_A]$), the subset $\mathit{\varphi }\left(t\right)={\mathit{\phi }}_{k+\mathrm{1}}(t-(k+\mathrm{1}\left)T_A\right)$ of controls corresponding to that window is released by the buffer and applied to the farm. In a real-time setting, the update time T_A should be long enough to compensate for the computational time of estimation, prediction and optimization. Similar schemes have been proposed in, e.g., Chen et al. (2000), Findeisen and Allgöwer (2004) and Su et al. (2013) for low-dimensional generic demonstration problems.

Figure 2 (a) Time-decoupled MPC assuming moving horizon estimation. In the window $t\in [t^k,t^k+T_A]$, a state estimate $\widehat{\mathit{q}}\left(t^k\right)$ is computed based on past measurements h(q(t)) stored in the estimation buffer; controls are computed with offset T_A based on a prediction $\widehat{\mathit{q}}(t^k+T_A)$ of the future state and buffered until they are applied in the next window $t\in [t^k+T_A,t^k+\mathrm{2}{T}_A]$. (b) Time-decoupled MPC approach considered in this work, assuming perfect state information. Controls and prediction are computed using coarse-grid LES; fine-grid LES serves as an emulator for the actual farm.

The time-decoupled control loop can be applied in a receding-horizon framework. The sequence of real-time computations versus the corresponding receding-horizon computations (in simulated time) is outlined in Fig. 3. As described above, the state estimate is typically generated based on past flow measurements over some estimation horizon T_SE, whereas the optimization is performed over a time horizon T with offset T_A.

Figure 3 Sequence of real-time computations (estimation, prediction and optimization) versus corresponding receding-horizon computations in simulated time, corresponding to the control loop from Fig. 2a.

To reduce the computational time, as in Bauweraerts and Meyers (2019), the prediction and optimization are performed using a coarse-grid wind farm model. The actual wind farm in Fig. 2a is represented by an emulator in the form of fine-grid LES, to which the control action ϕ is applied. In the present study, we assume perfect state knowledge and hence omit the state estimator. Instead, we just introduce a restriction operator to filter the exact wind farm state onto the coarse prediction and optimization grids. Under this simplification, the exact state (but restricted to the coarse grid) is then propagated forward in time by T_A in the predictor to generate the prediction $\widehat{\mathit{q}}(t^k+T_A)$ of the future wind farm state. The restriction operator, used to map system feedback from the fine LES grid to the coarser prediction and optimization grid, is discussed in more detail in Sect. 2.5. The time-decoupled control loop, omitting state estimation, is shown in Fig. 2b. Note that this approach introduces two sources of model mismatch between the coarse models and the fine-grid emulator. Firstly, a restriction error arises from filtering the LES data from the emulator to the coarser resolution of the predictor and optimizer. Secondly, the predictor introduces a prediction error that is inherent to time-decoupled MPC and depends on the update time T_A.

2.2 Wind farm optimization problem

In the receding-horizon approach, in every optimization window, the overall wind farm power is optimized over the optimization horizon T. The optimization problem is formulated as in Munters and Meyers (2018b):

$$\begin{array}{}\text{(1)}& {\displaystyle }{\min }_{\begin{array}{c}\mathit{\phi }\left(t\right),\mathit{q}\left(t\right)\\ \mathrm{0}<t\le T\end{array}}\mathcal{J}(\mathit{\phi },\mathit{q})=-\underset{\mathrm{0}}{\overset{T}{\int }}\sum _{m=\mathrm{1}}^{N_t}{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{C}_{P,m}^{\prime }{\stackrel{\mathrm{‾}}{V}}_m^{\mathrm{3}}{A}_m\mathrm{d}t,\text{(2)}& {\displaystyle }\begin{array}{rl}\mathrm{s}.\mathrm{t}.{\displaystyle \frac{\partial \mathit{u}}{\partial t}}& +(\mathit{u}\cdot \mathrm{\nabla })\mathit{u}=-{\displaystyle \frac{\mathrm{\nabla }p}{\mathit{\rho }}}-\mathrm{\nabla }\cdot {\mathit{\tau }}_{\mathrm{sgs}}\\ & -\sum _{m=\mathrm{1}}^{N_t}{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\widehat{C}}_{\mathrm{T},m}^{\prime }{\stackrel{\mathrm{‾}}{V}}_m^{\mathrm{2}}{\mathcal{R}}_m\left(\mathit{x}\right){\mathit{e}}_{⟂,m}\mathrm{in}\mathrm{\Omega }\times (\mathrm{0},T],\end{array}\text{(3)}& {\displaystyle }\mathrm{\nabla }\cdot \mathit{u}=\mathrm{0}\mathrm{in}\mathrm{\Omega }\times (\mathrm{0},T],\text{(4)}& {\displaystyle }\mathit{\tau }{\displaystyle \frac{\mathrm{d}{\widehat{C}}_{\mathrm{T},m}^{\prime }}{\mathrm{d}t}}=C_{\mathrm{T},m}^{\prime }-{\widehat{C}}_{\mathrm{T},m}^{\prime }m=\mathrm{1}\mathrm{\dots }{N}_t\mathrm{in}(\mathrm{0},T],\text{(5)}& {\displaystyle }{\displaystyle \frac{\mathrm{d}{\mathit{\theta }}_m}{\mathrm{d}t}}={\mathit{\omega }}_{m}m=\mathrm{1}\mathrm{\dots }{N}_t\mathrm{in}(\mathrm{0},T],\text{(6)}& {\displaystyle }{C}_{\mathrm{T},\min }^{\prime }\le C_{\mathrm{T},m}^{\prime }\le C_{\mathrm{T},\max }^{\prime }m=\mathrm{1}\mathrm{\dots }{N}_t\mathrm{in}(\mathrm{0},T],\text{(7)}& {\displaystyle }{\mathit{\omega }}_{\min }\le {\mathit{\omega }}_m\le {\mathit{\omega }}_{\max }m=\mathrm{1}\mathrm{\dots }{N}_t\mathrm{in}(\mathrm{0},T],\text{(8)}& {\displaystyle }\mathit{u}(\mathit{x},\mathrm{0})={\widehat{\mathit{u}}}_{\mathrm{0}}^{\mathrm{pred}}\mathrm{in}\mathrm{\Omega }.\end{array}$$

The flow through the wind farm is governed by the LES equations (Eqs. 2–3), with u and p the velocity and pressure. Subgrid scales are represented by the stress tensor τ_sgs using a standard Smagorinsky model with constant coefficient C_s=0.14 including wall damping. For the discretization, we employ a pseudo-spectral method in the streamwise and spanwise direction and a fourth-order energy-conservative finite-difference scheme in the vertical direction. For the time stepping, we use an explicit fourth-order Runge–Kutta scheme. Here, all four Runge–Kutta stages are stored on disc (as opposed to the aforementioned LES studies that only store the first stage). This results in a more accurate representation of the gradients using the adjoint method (see Sect. 2.4). In the vertical direction, a high Reynolds number wall stress boundary condition and symmetry boundary condition are imposed on the bottom and top surface of the domain respectively. All simulations are performed using the in-house simulation code SP-Wind (for more information, see, e.g., Goit and Meyers, 2015; Munters and Meyers, 2017, 2018b).

Each turbine m (=1 … N_t) is controlled by a time-dependent thrust coefficient set point $C_{\mathrm{T},m}^{\prime }$ and yaw rate ω_m, both subject to box constraints (Eqs. 6–7). Put together, they constitute the vector of optimization variables,

$$\begin{array}{ll}{\displaystyle }\mathit{\phi }\left(t\right)& {\displaystyle }=\left[C_{\mathrm{T},\mathrm{1}}^{\prime }\left(t\right),\mathrm{\dots },C_{\mathrm{T},N_t}^{\prime }\left(t\right),{\mathit{\omega }}_{\mathrm{1}}\left(t\right),\mathrm{\dots },{\mathit{\omega }}_{N_t}\left(t\right)\right]\\ {\displaystyle }& {\displaystyle }=\left[{\mathit{C}}_{\mathrm{T}}^{\prime }\left(t\right),\mathit{\omega }\left(t\right)\right],\end{array}$$

for $\mathrm{0}<t\le T$. To model the turbine response time, an exponential filter (Eq. 4) with time constant τ is applied to the thrust coefficient set points, resulting in the time-filtered disc-based thrust coefficients ${\widehat{\mathit{C}}}_{\mathrm{T}}^{\prime }\left(t\right)$ (Munters and Meyers, 2017). The yaw rates control the turbine yaw angles θ(t) through the yaw equation (Eq. 5). All together, this results in the vector of state variables $\mathit{q}=\left[\mathit{u}\right(\mathit{x},t),p(\mathit{x},t),{\widehat{\mathit{C}}}_{\mathrm{T}}^{\prime }(t),\mathit{\theta }(t\left)\right]$. The disc-based thrust coefficients and yaw angles determine the thrust force exerted on the flow and the power extracted by the turbine (see Sect. 2.3). For a more detailed explanation on all the terms and equations, see Goit and Meyers (2015) and Munters and Meyers (2017, 2018b).

Note that controls are computed ahead of time in line with the time-decoupled MPC loop from Sect. 2.1. In receding-horizon interval k (for $t\in [t^k,t^k+T_A]$, see Fig. 3), the exact state from the fine-grid emulator is first restricted to the coarse grid, resulting in $\widehat{\mathit{q}}\left(t^k\right)$. Next, the prediction $\widehat{\mathit{q}}(t^k+T_A)$ is computed by propagating $\widehat{\mathit{q}}\left(t^k\right)$ over T_A using the same coarse-grid LES model described by Eqs. (2)–(7). This estimate is used as initial condition ${\widehat{\mathit{u}}}_{\mathrm{0}}^{\mathrm{pred}}$ for the optimization; see Eq. (8). Also note that the controls ϕ(t), which are actually applied to the wind farm emulator in the next time interval (for $t\in [t^k+T_A,t^k+\mathrm{2}{T}_A]$), only comprise the first part of length T_A from the optimal controls φ_k+1(t) (with T_A≤T).

2.3 Wind turbine modeling

For the turbine modeling, a standard non-rotating actuator disc model is used. Based on actuator disc theory, the turbines exert a force on the flow: ${\mathit{f}}_m=-\frac{\mathrm{1}}{\mathrm{2}}{\widehat{C}}_{\mathrm{T},m}^{\prime }{\stackrel{\mathrm{‾}}{V}}_m^{\mathrm{2}}{\mathcal{R}}_m\left(\mathit{x}\right){\mathit{e}}_{⟂,m}$, where ℛ_m is a smoothed footprint of the rotor on the LES grid, ${\mathit{e}}_{⟂,m}$ the unit vector perpendicular to the rotor plane and ${\stackrel{\mathrm{‾}}{V}}_m=\frac{M}{A_m}{\int }_{\mathrm{\Omega }}{\mathcal{R}}_m\mathit{u}\cdot {\mathit{e}}_{⟂,m}\mathrm{d}\mathit{x}$ the (corrected) disc-averaged velocity (with A_m the rotor disc area and M a correction factor defined below). The power extracted from turbine m is then given by $P_m=\frac{\mathrm{1}}{\mathrm{2}}{C}_{\mathrm{p},m}^{\prime }{\stackrel{\mathrm{‾}}{V}}_m^{\mathrm{3}}{A}_m$, where $C_{\mathrm{p},m}^{\prime }$ denotes the disc-based power coefficient.

Figure 4 Empirical power coefficient versus disc-based thrust coefficient for different correction strategies at different yaw angles θ: linear scaling $C_{\mathrm{P}}^{\prime }=a{\widehat{C}}_{\mathrm{T}}^{\prime }$ and factor from Shapiro et al. (2019) in Eq. (9) for θ=0^∘ (a) and θ=30^∘ (b), and lookup table approach for θ=0, 12.5 and 27.5^∘ (c). For every simulation, uniform inflow velocity U_∞=8 m s⁻¹ is prescribed using a fringe region spanning the final 20 % of the domain, with domain size L_x=26.92 D, L_y=13.46 D and L_z=8.41 D on the coarse grid $\mathrm{\Delta }x=\mathrm{\Delta }y=\mathrm{1.6}\mathrm{\Delta }z=\mathrm{80}$ m. In each plot, we also show the empirical power coefficient (including Shapiro correction) at the reference resolution $\mathrm{\Delta }x=\mathrm{\Delta }y=\mathrm{2}\mathrm{\Delta }z=\mathrm{13.33}$ m.

On present-day grid resolutions, power is typically overestimated due to the diffuse smearing of the rotor disc on the simulation grid by the rotor footprint (Martínez-Tossas et al., 2015; Shapiro et al., 2019). To account for this, Munters and Meyers (2017, 2018b) have proposed to set $C_{\mathrm{p}}^{\prime }=a{\widehat{C}}_{\mathrm{T}}^{\prime }$, where a is selected based on fitting LES data to 1D momentum theory. While linear scaling is effective on intermediate grid resolutions for unyawed turbines, its performance deteriorates on coarser grids and for yawed turbines. This is illustrated in Fig. 4a and b, where we show the empirical power coefficient $C_{\mathrm{p}}=P/\left(\frac{\mathrm{1}}{\mathrm{2}}{U}_{\mathrm{\infty }}^{\mathrm{3}}A\right)$ for the DTU 10 MW turbine for different values of a at yaw angles θ=0^∘ and θ=30^∘, computed using LES with uniform inflow (U_∞=8 m s⁻¹) on a grid resolution $\mathrm{\Delta }x=\mathrm{\Delta }y=\mathrm{1.6}\mathrm{\Delta }z=\mathrm{80}$ m (the coarsest resolution in this work). We also show the effect of the correction from Shapiro et al. (2019), which is expressed as a correction factor on the disc-averaged velocity $V_m=\frac{\mathrm{1}}{A_m}\underset{\mathrm{\Omega }}{\int }{\mathcal{R}}_m\mathit{u}\cdot {\mathit{e}}_{⟂,m}\mathrm{d}\mathit{x}$. In particular, the corrected disc-averaged velocity is given by ${\stackrel{\mathrm{‾}}{V}}_m=M^{\mathrm{sh}}{V}_m$, with

$$\begin{array}{}\text{(9)}& M^{\mathrm{sh}}={\left(\mathrm{1}+{\displaystyle \frac{{\widehat{C}}_{\mathrm{T},m}^{\prime }}{\mathrm{4}}}{\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{3}\mathit{\pi }}}}{\displaystyle \frac{\mathrm{\Delta }}{R}}\right)}^{-\mathrm{1}}\end{array}$$

the Shapiro factor, R the rotor radius and Δ the filter width of the Gaussian filtering kernel. For comparison, the power coefficient (for a=1.0) on the reference resolution $\mathrm{\Delta }x=\mathrm{\Delta }y=\mathrm{2}\mathrm{\Delta }z=\mathrm{13.33}$ m (the resolution of the emulator) is also included, where we also applied the Shapiro factor to better replicate 1D momentum theory.

In Fig. 4a, the correction a=0.775 performs relatively well compared to the reference. However, under yaw misalignment of 30^∘ in Fig. 4b, the same factor overestimates power for higher thrust coefficients. The Shapiro correction consistently underestimates power on the coarse grid. In an optimal control setting, any mismatch between model and reference will result in suboptimal thrust coefficients and yaw angles. Therefore, we take a=1 and propose a lookup table approach to account for discrepancies in power prediction on different grid resolutions. In particular, the disc-averaged velocity on the coarse grids is corrected based on a lookup value M that depends on the thrust coefficient and the yaw angle:

$$\begin{array}{}\text{(10)}& {\stackrel{\mathrm{‾}}{V}}_m=M\left({\widehat{C}}_{\mathrm{T},m}^{\prime },{\mathit{\theta }}_m\right)V_m.\end{array}$$

The lookup table is constructed in such a way that, for uniform inflow, the corrected disc-averaged velocity (and consequently the power output) on the coarse grids matches the disc-averaged velocity from the reference for given (${\widehat{C}}_{\mathrm{T}}^{\prime },\mathit{\theta }$). This approach allows us to match power calculation on the coarse grids arbitrarily well with a given reference turbine (depending on the lookup table resolution in terms of ${\widehat{C}}_{\mathrm{T}}^{\prime }$ and θ). In the current work, we propose a lookup range ${\widehat{\mathit{C}}}_{\mathrm{T}}^{\mathrm{table}}\times {\mathit{\theta }}^{\mathrm{table}}=\mathit{\{}\mathrm{0.0},\mathrm{0.5},\mathrm{\dots },\mathrm{2.5},\mathrm{3.0}\mathit{\}}\times \mathit{\{}\mathrm{0},\mathrm{5},\mathrm{\dots },\mathrm{35},\mathrm{40}{}^{\circ }\mathit{\}}$, and we use bilinear interpolation on the lookup table to compute the correction for given ${\widehat{C}}_{\mathrm{T}}^{\prime }$ and θ. The resulting lookup tables are tabulated in Appendix A, where we also briefly examine the performance of the lookup table approach for turbulent inflow. The corrected power coefficient is shown in Fig. 4c.

2.4 Optimization method and gradient computation

As in Munters and Meyers (2018b), the optimization problem is solved in a reduced fashion by explicitly substituting state equations (Eqs. 2–5) into the cost function, i.e. by minimizing $\stackrel{\mathrm{̃}}{\mathcal{J}}\left(\mathit{\phi }\right)=\mathcal{J}(\mathit{\phi },\mathit{q}(\mathit{\phi }\left)\right)$ subject to box constraints (Eqs. 6–7). To solve the optimization problem, we use the limited-memory Broyden–Fletcher–Goldfarb–Shanno method with box constraints (L-BFGS-B). In each iteration, this quasi-Newton method constructs a quadratic approximation of the objective function, where the inverse Hessian is approximated by the BFGS formula. A search direction is then generated based on the minimum of the quadratic model (Nocedal and Wright, 2006). We use a line search method to determine a step length in the search direction that satisfies the strong Wolfe conditions. In every optimization window, this procedure results in a sequence of function and gradient evaluations (see Fig. 3). To compensate for the large-scale nature of the problems at hand, we use the limited-memory version of BFGS that only stores a limited number of correction pairs. We use the L-BFGS-B Fortran library from Zhu et al. (1997) and Morales and Nocedal (2011); for more information on the algorithm, the reader is referred to Byrd et al. (1995).

In contrast to previous LES-based wind farm control studies (Goit and Meyers, 2015; Goit et al., 2016; Munters and Meyers, 2017, 2018b) that relied on a continuous adjoint approach to compute gradients, we follow up on the work in Yilmaz and Meyers (2019) by using a temporally discrete adjoint method. In the continuous adjoint method, the adjoint equations are first derived based on the description of the optimization problem in Eqs. (1)–(8) and subsequently discretized and solved using LES similar to a forward simulation. Conversely, the discrete adjoint method first discretizes and linearizes the state equations and then formulates the discrete adjoint of the linearized equations. Consequently, the discrete adjoint method obtains the gradient of the discretized cost functional, whereas the continuous adjoint method yields a discrete approximation of the gradient of the continuous cost functional. In the limit of infinite grid resolution, both methods are equivalent (Giles and Pierce, 2000).

Below we directly formulate the temporally discrete adjoint Runge–Kutta 4 scheme, derived from a fourth-order Runge–Kutta discretization of the state equations (Eqs. 2–5) and cost function. More details on the discretization are provided in Appendix B1. For the Navier–Stokes equations, the thrust coefficient filter equation and the yaw equation, this results in the following (where i and j subscripts denote the Runge–Kutta stages, m is the turbine number, and n the discrete time instant, i.e. tⁿ=nΔt):

$$\begin{array}{}\text{(11)}& {\displaystyle }\begin{array}{rl}{\mathit{\xi }}_i^n/\mathrm{\Delta }t& =\left(-{\left(\mathrm{\nabla }{\mathit{u}}_i^n\right)}^T+\left.({\mathit{u}}_i^n\cdot \mathrm{\nabla }\right)\right){\widehat{\mathit{\xi }}}_i^n-\mathrm{\nabla }{\mathit{\pi }}_i^n/\mathit{\rho }\\ & -\mathrm{\nabla }\cdot {\mathit{\tau }}_{\mathrm{sgs}}^{*}\left({\mathit{u}}_i^n,{\widehat{\mathit{\xi }}}_i^n\right)+\sum _{m=\mathrm{1}}^{N_t}{\mathit{f}}_{m,i}^{{*}^n}i=\mathrm{1}\mathrm{\dots }\mathrm{4},\end{array}\text{(12)}& {\displaystyle }\begin{array}{rl}{\mathrm{\nabla }}^{\mathrm{2}}{\mathit{\pi }}_i^n/\mathit{\rho }& =\mathrm{\nabla }\cdot \left[\left(-{\left(\mathrm{\nabla }{\mathit{u}}_i^n\right)}^T+\left({\mathit{u}}_i^n\cdot \mathrm{\nabla }\right)\right){\widehat{\mathit{\xi }}}_i^n\right.\\ & \left.-\mathrm{\nabla }\cdot {\mathit{\tau }}_{\mathrm{sgs}}^{*}\left({\mathit{u}}_i^n,{\widehat{\mathit{\xi }}}_i^n\right)+\sum _{m=\mathrm{1}}^{N_t}{\mathit{f}}_{m,i}^{{*}^n}\right]i=\mathrm{1}\mathrm{\dots }\mathrm{4},\end{array}\text{(13)}& {\displaystyle }{\mathit{\xi }}^n={\mathit{\xi }}^{n+\mathrm{1}}+\sum _{i=\mathrm{1}}^{\mathrm{4}}{\mathit{\xi }}_i^n,\text{(14)}& {\displaystyle }\begin{array}{rl}{\mathit{\sigma }}_{m,i}^n/\mathrm{\Delta }t=& {\displaystyle \frac{\mathrm{1}}{\mathit{\tau }}}\left({\widehat{\mathit{\sigma }}}_{m,i}^n-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\stackrel{\mathrm{‾}}{V}}_{m,i}^{n^{\mathrm{2}}}\left(b_i{\stackrel{\mathrm{‾}}{V}}_{m,i}^n-{\widehat{X}}_{m,i}^n\right)A_m\right.\\ & -{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\widehat{C}}_{{\mathrm{T}}_m}^{\prime n}{\stackrel{\mathrm{‾}}{V}}_{m,i}^n{V}_{m,i}^n{\displaystyle \frac{\partial M\left({\widehat{C}}_{{\mathrm{T}}_m}^{\prime n},{\mathit{\theta }}_m^n\right)}{\partial {\widehat{C}}_{\mathrm{T}}^{\prime }}}\\ & \left.\left(\mathrm{3}{b}_i{\stackrel{\mathrm{‾}}{V}}_{m,i}^n-\mathrm{2}{\widehat{X}}_{m,i}^n\right)A_m\right)i=\mathrm{1}\mathrm{\dots }\mathrm{4},\end{array}\text{(15)}& {\displaystyle }{\mathit{\sigma }}_m^n={\mathit{\sigma }}_m^{n+\mathrm{1}}+\sum _{i=\mathrm{1}}^{\mathrm{4}}{\mathit{\sigma }}_{m,i}^n,\text{(16)}& {\displaystyle }\begin{array}{rl}& {\mathit{\eta }}_{m,i}^n/\mathrm{\Delta }t=-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\widehat{C}}_{{\mathrm{T}}_m}^{\prime n}{\stackrel{\mathrm{‾}}{V}}_{m,i}^n\\ & \left[\underset{\mathrm{\Omega }}{\int }\left(\left(\mathrm{3}{b}_i{\stackrel{\mathrm{‾}}{V}}_{m,i}^n-\mathrm{2}{\widehat{X}}_{m,i}^n\right)M{\mathit{u}}_i^n-{\stackrel{\mathrm{‾}}{V}}_{m,i}^n\mathit{\xi }\right)\right.\\ & \left.\cdot \left({\mathit{e}}_{\left|\right|,m}{\mathcal{R}}_m+{\mathit{e}}_{⟂,m}{\mathcal{D}}_m\right)d\mathit{x}\right]-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\widehat{C}}_{{\mathrm{T}}_m}^{\prime n}{\stackrel{\mathrm{‾}}{V}}_{m,i}^n\\ & V_{m,i}^n{\displaystyle \frac{\partial M\left({\widehat{C}}_{{\mathrm{T}}_m}^{\prime n},{\mathit{\theta }}_m^n\right)}{\partial \mathit{\theta }}}\left(\mathrm{3}{b}_i{\stackrel{\mathrm{‾}}{V}}_{m,i}^n-\mathrm{2}{\widehat{X}}_{m,i}^n\right)A_m\\ & i=\mathrm{1}\mathrm{\dots }\mathrm{4},\end{array}\text{(17)}& {\displaystyle }{\mathit{\eta }}_m^n={\mathit{\eta }}_m^{n+\mathrm{1}}+\sum _{i=\mathrm{1}}^{\mathrm{4}}{\mathit{\eta }}_{m,i}^n,\end{array}$$

with a_ij and β_i the Runge–Kutta coefficients. In these equations, ξⁿ, πⁿ, ${\mathit{\sigma }}_m^n$ and ${\mathit{\eta }}_m^n$ are the adjoint state variables associated with uⁿ, pⁿ, ${\widehat{C}}_{\mathrm{T},m}^{{}^{\prime }n}$, and ${\mathit{\theta }}_m^n$, and ${\mathit{\xi }}_i^n$, ${\mathit{\pi }}_i^n$, ${\mathit{\sigma }}_{i,m}^n$, and ${\mathit{\eta }}_{i,m}^n$ for i=1 … 4 are the corresponding adjoint Runge–Kutta stages. Furthermore, ${\widehat{\mathit{\xi }}}_i^n$, ${\widehat{\mathit{\sigma }}}_{i,m}^n$ and ${\widehat{\mathit{\eta }}}_{i,m}^n$ are auxiliary variables (see Appendix B). As in Munters and Meyers (2018b), ${\mathit{f}}_m^{*}$ denotes the adjoint turbine force, X_m is the disc-averaged adjoint velocity, ${\mathit{e}}_{\parallel ,m}$ is the rotor-parallel unit vector and 𝒟_m is the rotational rotor footprint. The temporally discrete adjoint equations are derived in detail in Appendix B. For a detailed definition of all terms related to the turbine modeling, see Goit and Meyers (2015) and Munters and Meyers (2018b). However, note that these studies only used the first forward Runge–Kutta stage, whereas here all four (forward) Runge–Kutta stages are used in the adjoint scheme, resulting in a more accurate method.

Finally, the adjoint variables are used to compute the gradients with respect to the thrust coefficient set points and yaw rates:

$$\begin{array}{}\text{(18)}& \begin{array}{rl}{\mathrm{\nabla }}_{{\mathit{\phi }}_m^n}{\stackrel{\mathrm{̃}}{\mathcal{J}}}^N& =\left(\begin{array}{c}\frac{\partial {\stackrel{\mathrm{̃}}{\mathcal{J}}}^N}{\partial {C_{\mathrm{T}}}_{,m}^{\prime n}}\\ \frac{\partial {\stackrel{\mathrm{̃}}{\mathcal{J}}}^N}{\partial {\mathit{\omega }}_m^n}\end{array}\right)\\ & =\left(\begin{array}{c}-\mathrm{\Delta }t{\sum }_{i=\mathrm{1}}^{\mathrm{4}}\left(b_i{\mathit{\sigma }}_m^{n+\mathrm{1}}+{\sum }_{j=\mathrm{1}}^{i-\mathrm{1}}{a}_{ij}{\mathit{\sigma }}_{m,i}^n\right)\\ -\mathrm{\Delta }t{\sum }_{i=\mathrm{1}}^{\mathrm{4}}\left(b_i{\mathit{\eta }}_m^{n+\mathrm{1}}+{\sum }_{j=\mathrm{1}}^{i-\mathrm{1}}{a}_{ij}{\mathit{\eta }}_{m,i}^n\right)\end{array}\right).\end{array}\end{array}$$

Note that, as mentioned above, Eq. (18) is an exact expression for the gradient of the discretized objective function. The discrete adjoint approach and gradient computation are validated in Appendix B5.

2.5 Coarse-grid optimization and coarsening strategy

As discussed in Sect. 1, we investigate the influence of the spatiotemporal grid resolution of the LES wind farm control model on the overall power gain and computational speed. To that end, as in Bauweraerts and Meyers (2019), we define three coarse grid resolutions for prediction and optimization, as well as a fine reference grid for the wind farm emulator. For the grid specifications, the reader is referred to Sect. 3.1, where we discuss the case setup. Below, we discuss the coarsening strategy in view of the time-decoupled MPC loop from Fig. 2b.

2.5.1 Coarse-grid prediction and control methodology

At the start of every window, feedback (i.e. the 3D flow field) from the fine-grid wind farm emulator (the reference) is provided to the coarse predictor through a restriction operator (see Fig. 2b). In view of computational time, we allow different domain sizes for the LES models in the coarse-grid predictor and optimizer. In particular, for the predictor, we propose an upstream domain length proportional to the total prediction horizon T+T_A, i.e. $L_{x,\mathrm{upstream}}=\mathit{\alpha }(T+T_A)U_{\mathrm{\infty }}$ (where α≥1 is a safety factor). The predictor uses this domain to propagate the restricted reference field, which in turn yields the initial condition for the optimization. For the optimization, we take $L_{x,\mathrm{upstream}}=\mathit{\alpha }T{U}_{\mathrm{\infty }}$. Thus, the approach is characterized by cropping and restricting the reference field to the prediction domain, and subsequently another cropping from the prediction to the optimization domain, as graphically illustrated in Fig. 5. The horizon-dependent upstream domain lengths ensure that inflow never reaches the front-row turbines within the prediction and optimization windows, rendering fringe regions and turbulent inflow generation superfluous in the predictor and optimizer. The fringe region from the fine-grid reference simulation is therefore excluded upon restriction, which, along with the smaller domain sizes (for a given grid resolution), entails significant computational speedups. We note that the lack of proper inflow may influence the optimized controls towards the end of the optimization window. However, choosing an update time T_A<T should suffice to counter these effects. Also note that, besides the horizon-dependent streamwise cropping, we also allow for a spanwise and vertical cropping to further accelerate prediction and optimization.

Figure 5 Coarsening strategy: the fine-grid reference field is cropped (removing the fringe region and retaining an upstream domain length of αTU_∞) and restricted to the coarser resolution for the optimization. The resulting field is used as initial condition for the optimization.

2.5.2 Restriction from reference to optimization resolution

Given the pseudo-spectral discretization in SP-Wind, the reference flow field from the emulator must be transformed from Fourier space into real space before applying the cropping. The cropped reference velocity in real space, ${\widehat{\mathit{u}}}^{\frac{\mathrm{3}}{\mathrm{2}},\mathrm{ref}}$, is then restricted to the coarser resolution used in the predictor and optimizer using linear interpolation (superscript “ref” denotes the reference):

$$\begin{array}{}\text{(19)}& \begin{array}{rl}{\widehat{\mathit{u}}}^{\frac{\mathrm{3}}{\mathrm{2}}}\left(x^{\frac{\mathrm{3}}{\mathrm{2}}},y^{\frac{\mathrm{3}}{\mathrm{2}}},z\right)=& \sum _{i=\mathrm{1}}^{\frac{\mathrm{3}}{\mathrm{2}}{N}_x^{\mathrm{ref}}}\sum _{j=\mathrm{1}}^{\frac{\mathrm{3}}{\mathrm{2}}{N}_y^{\mathrm{ref}}}\sum _{k=\mathrm{1}}^{N_z^{\mathrm{ref}}}{\widehat{\mathit{u}}}^{\frac{\mathrm{3}}{\mathrm{2}},\mathrm{ref}}\left(x_i^{\frac{\mathrm{3}}{\mathrm{2}},\mathrm{ref}},y_j^{\frac{\mathrm{3}}{\mathrm{2}},\mathrm{ref}},z_k^{\mathrm{ref}}\right)\\ & \max \left(\mathrm{1}-\left|{\displaystyle \frac{z-z_k^{\mathrm{ref}}}{\mathrm{\Delta }{z}^{\mathrm{ref}}}}\right|,\mathrm{0}\right)\\ & \max \left(\mathrm{1}-\left|{\displaystyle \frac{y^{\frac{\mathrm{3}}{\mathrm{2}}}-y_j^{\frac{\mathrm{3}}{\mathrm{2}},\mathrm{ref}}}{\frac{\mathrm{2}}{\mathrm{3}}\mathrm{\Delta }{y}^{\mathrm{ref}}}}\right|,\mathrm{0}\right)\\ & \max \left(\mathrm{1}-\left|{\displaystyle \frac{x^{\frac{\mathrm{3}}{\mathrm{2}}}-x_i^{\frac{\mathrm{3}}{\mathrm{2}},\mathrm{ref}}}{\frac{\mathrm{2}}{\mathrm{3}}\mathrm{\Delta }{x}^{\mathrm{ref}}}}\right|,\mathrm{0}\right),\end{array}\end{array}$$

where ${\widehat{\mathit{u}}}^{\frac{\mathrm{3}}{\mathrm{2}}}$ denotes the restricted velocity (in real space) and Δx^ref, Δy^ref and Δz^ref denote the grid resolutions of the reference grid. Note that the restriction takes place in real space, where the grid is a factor of $\mathrm{3}/\mathrm{2}$ finer in the horizontal directions for de-aliasing using Orszag's rule (hence the factors of $\mathrm{3}/\mathrm{2}$ in the equation above).

3 Case setup

This section provides a detailed description of the simulation cases and numerical setup that are used to evaluate the proposed wind farm controller. Throughout this work, we consider the TotalControl Reference Wind Power Plant (TCRWP) (Andersen et al., 2018) consisting of 32 DTU 10 MW turbines arranged in an 8×4 aligned pattern, as illustrated in Fig. 6. Turbines are separated by an intermediate spacing of 5 D in streamwise and spanwise directions, where D=178.3 m is the rotor diameter and turbines are placed at hub height z_h=119 m (based on the DTU 10 MW reference turbine, as reported in Bak et al., 2013).

Figure 6 Layout of the TotalControl Reference Wind Power Plant. Axes in rotor diameter units, with D=178.3 m. Figure adapted from Andersen et al. (2018).

3.1 Simulation setup

The grid resolutions for the numerical discretization on the three coarseness levels and the reference are summarized in Table 1. Since SP-Wind requires a constant integration time step, $T/\mathrm{\Delta }t$ should be an integer. Furthermore, we use the same Δtⁱ for all cases on grid level i. Therefore, the time steps in Table 1 are selected as the largest possible ones meeting these requirements and adhering to a Courant–Friedrichs–Lewy (CFL) condition of 0.8. Note that based on the CFL condition only, the time step on the coarsest grid levels could in principle still be higher, hence speeding up the computation. The latter case is investigated in Sect. 5, where we design a controller that is as close to real time as possible.

Table 1 Grid resolutions for the different coarseness levels.

For the fine-grid reference simulation in the wind farm emulator, we take the simulation setup from Andersen et al. (2018) and Sood and Meyers (2020), consisting of $\mathrm{1200}\times \mathrm{1200}\times \mathrm{225}$ grid cells and a simulation domain with dimensions $L_x\times L_y\times L_z=\mathrm{16}\times \mathrm{16}\times \mathrm{1.5}$ km³, where the final 6.25 % of the domain in the streamwise direction is used as a fringe region to impose the inflow conditions. The spatial extent of the domain suffices to keep blockage effects negligible in the reference simulation. The simulations are performed using a standard offshore roughness length $z_{\mathrm{0}}=\mathrm{2}\times {\mathrm{10}}^{-\mathrm{4}}$ m, and the flow is driven by a pressure gradient ${\partial }_x{p}_{\mathrm{\infty }}/\mathit{\rho }=\mathrm{5.2267}\times {\mathrm{10}}^{-\mathrm{5}}$ m s⁻² resulting in a friction velocity u_τ=0.28 m s⁻¹, which is a typical value in offshore boundary layers.

Table 2 Grid specifications and domain sizes for the different coarseness levels as a function of the optimization horizon T and control update time T_A for the optimizer and predictor. Domain sizes in spanwise and vertical directions are equal for all cases: L_y=6.24 km and L_z=1.5 km. Number of grid cells in spanwise and vertical direction ($N_y^i$ and $N_z^i$) are equal for the predictor and optimizer.

Figure 7 Snapshot of the initial condition for optimal control (on the reference resolution). Colors represent the instantaneous velocity magnitude [m s⁻¹]. The black dots represent wind turbine locations.

Prior to the wind farm simulations, a pressure-driven precursor simulation with periodic boundary conditions was run on the same (reference) domain to generate turbulent inflow. With the roughness length and pressure gradient specified above, this results in a freestream wind speed roughly equal to 9.4 m s⁻¹ at the turbine hub height. The precursor data (with a detailed overview of the precursor simulation setup) is publicly available in Munters et al. (2019). Using this precursor, the flow is then advanced through the wind farm for a spin-up period of 60 min to account for startup transients. Note that here turbines are modeled by non-rotating actuator discs, and the turbine locations are shifted backwards in the streamwise direction in comparison to Sood and Meyers (2020). The latter is required to accommodate the entire flow field encompassed in the longest prediction windows (see Sect. 2.5). The resulting flow field, depicted in Fig. 7, is used as the starting point for the optimization. Figure 8 shows the initial flow field after restriction to the resolutions of the coarse control models from Table 1.

Figure 8 Snapshot of the initial condition for optimal control on the reference resolution and coarser optimization resolution from Table 1. z slices are taken at hub height, and y slices are taken at the location of the first column of turbines. Colors represent the instantaneous velocity magnitude [m s⁻¹]. The black dots represent wind turbine locations. Grid resolutions are also indicated.

As explained in Sect. 2.5, the upstream domain lengths for the coarse models in the predictor and optimizer are chosen proportional to the optimization horizon T and control update time T_A. In this work, we consider four different horizons (see also Sect. 3.1.1): $T,T_A\in \mathit{\{}\mathrm{50},\mathrm{150},\mathrm{250},\mathrm{350}\mathit{\}}$. The corresponding domain sizes and grid specifications for each combination of T and T_A are summarized in Table 2 for both the predictor and optimizer. In all simulations, a high-Reynolds number wall model is used at the bottom of the domain with roughness length $z_{\mathrm{0}}=\mathrm{2}\times {\mathrm{10}}^{-\mathrm{4}}$ m, and the top boundary is treated by a stress-free condition. For the prediction and optimization, we use the same coarse grid resolutions from Table 1 and omit the fringe region (see Sect. 2.5).

3.1.1 Receding-horizon optimal control setup

Wind farm operation is optimized over a total horizon of T_tot=30 min, which comprises just under three through flows for the given wind farm (given a free-stream velocity U_∞≈9.4 m s⁻¹). In the receding-horizon framework, turbine controls are optimized over windows of horizon T with offset T_A (the control update time) corresponding to the prediction horizon of the predictor. Table 3 summarizes the different combinations of T and T_A considered here. The longest optimization window (T=350 s) allows the optimizer to account for wake interactions over three consecutive rows. With every horizon reduction of 100 s in Table 3, the optimizer looses control authority over one row of turbines in the wakes. For T=50 s, wakes cannot propagate to the next row within the optimization window; case 10 may therefore be considered as an “uncoordinated” control case. Note that optimization and prediction horizons are limited due to the natural divergence of trajectories in chaotic flows. In practice, for our setup, we find that gradients are still accurately represented for T=350 s (the longest optimization horizon considered in this work; see also Appendix B5 for the gradient verification).

Table 3 Optimization horizons T and control update times T_A.

All cases are initialized from the same flow field depicted in Fig. 7, which was generated by advancing the flow on the reference grid using constant thrust coefficients for all turbines ($C_{\mathrm{T}}^{\prime }=\mathrm{2}$, no yawing ω=0^∘ s⁻¹), until a statistically stationary state is achieved. As explained in Sect. 2.5, before each optimization run, the current flow field is taken from the reference simulation and restricted to the coarser prediction grid and propagated over the control update time T_A in the predictor (Fig. 5). Turbine controls are then optimized over the optimization horizon T, starting from initial guess $C_{\mathrm{T}}^{\prime }=\mathrm{2}$ and ω=0^∘ s⁻¹ for all turbines, until a stopping criterion is met. The convergence criterion used here is based on the relative improvement of the objective function over the L-BFGS-B iterations, i.e. $({\mathcal{J}}^{k-\mathrm{1}}-{\mathcal{J}}^k)/{\mathcal{J}}^{k-\mathrm{1}}\le \mathrm{5}\times {\mathrm{10}}^{-\mathrm{6}}$. The optimized controls are then applied to the reference in the next time window. Since the time-dependent controls are optimized on the temporal grid of the optimizer, the optimized controls are first interpolated onto the finer temporal reference grid using a simple zero-order hold rule.

3.1.2 Turbine control cases

Three control scenarios are examined; see Table 4. First, we define a steady reference case (R) where turbines operate at Betz-optimal thrust coefficients $C_{\mathrm{T}}^{\prime }=\mathrm{2}$, aligned with the mean-flow direction (θ=0^∘ and ω=0^∘ s⁻¹). Next, we consider a combined induction and yaw control case (IY) with a maximum yaw rate ω_max=0.4^∘ s⁻¹. The induction control part is restricted to the underinduction regime (i.e. $C_{\mathrm{T},\max }^{\prime }=\mathrm{2}$) to avoid bias in the results due to the inherent inaccuracy of the ADM in the overinduction regime. A response time τ=15 s is adopted for the time filtering of the thrust coefficient set points. Finally, we also consider the steady yaw control case from Sood and Meyers (2022), who used the recursive wake merging methodology from Lanzilao and Meyers (2022) on the Bastankhah wake model (Bastankhah and Porté-Agel, 2016) in a basic optimization framework to determine the yawing set points for the TCRWP, subject to a maximum yawing angle of 30^∘. For the steady yaw case here, we simply take their set points and initialize the turbines using these set points at the start of the simulation.

Table 4 Specifications for the turbine control cases.

4 Results and discussion

This section presents and discusses the results of the optimal control cases. All simulations are conducted on the wICE supercomputing platform of the VSC (Vlaams Supercomputer Centrum), using Ice Lake nodes containing two Intel Xeon Platinum 8360Y CPUs (36 cores each).

4.1 Allocation of resources

The focus of our study is to evaluate the performance of the control models in relation to their computational cost. Computational cost is measured in wall times, which depend on the number of cores used and the spatial parallelization. For the spatial parallelization, we employ a 2D domain decomposition, similar to the method used in earlier studies involving SP-Wind (see, e.g., Goit and Meyers, 2015; Goit et al., 2016; Munters and Meyers, 2018b, 2017).

For each grid resolution, since the time integration in the forward and adjoint simulations is the predominant contribution in the overall wall time, we select the number of compute cores that minimizes the wall time per Runge–Kutta step. Using a maximum of one whole compute node, these scaling tests reveal an optimum of 30, 54 and 72 cores for grid level 0, 1 and 2 respectively. Note that, for every simulation, we reserve the full compute node and then allocate the optimal number of cores using a “bunch” processor mapping to distribute the cores evenly over both sockets of the Ice Lake node.

4.2 Convergence behavior

Table 5 reports the average number of PDE evaluations (i.e. sum of forward and adjoint simulations) per optimization window required for formal convergence as specified by the convergence criterion from Sect. 3.1.1. Figure 9 shows the L-BFGS-B iterations versus the number of PDE evaluations for the optimization window starting at $t+T_A=\mathrm{750}$ s. To maintain clarity, only cases 4, 7, 9 and 10 (for which T_A=50 s) are depicted in the figure. As expected, the number of PDE evaluations increases with the optimization horizon, as this increases the number of optimization variables. In general, higher resolutions also require more function evaluations.

Table 5 Number of PDE evaluations and L-BFGS-B iterations per optimization run, averaged over the windows, for the different grid resolutions. The standard deviation on the number of PDE evaluations is also shown. The first window is excluded to account for the startup of the controller.

Figure 9 Cost function versus number of PDE evaluations for cases 4, 7, 9 and 10 (with T_A=50 s) on each grid level for the optimization window starting at $t+T_A=\mathrm{750}$ s. Circles mark L-BFGS-B iterations. For every grid level i, the cost function ${\widehat{\mathcal{J}}}^i$ is scaled by ${\widehat{\mathcal{J}}}_{\mathrm{0}}^i$ of the first iteration.

4.3 Power gains versus computational time

Figure 10 shows the performance of the proposed controller versus the computational cost for the three grid resolutions from Table 1 and the different combinations of T and T_A. Error bars indicating the 95 % confidence intervals are also depicted.

The performance of the controllers is measured in terms of the power gains η_P and farm efficiency η_farm:

$$\begin{array}{}\text{(20)}& {\mathit{\eta }}_{\mathrm{P}}={\displaystyle \frac{{\stackrel{\mathrm{‾}}{P}}_{\mathrm{farm}}}{{\stackrel{\mathrm{‾}}{P}}^{\mathrm{ref}}}},{\mathit{\eta }}_{\mathrm{farm}}={\displaystyle \frac{{\stackrel{\mathrm{‾}}{P}}_{\mathrm{farm}}}{N_t{\stackrel{\mathrm{‾}}{P}}_{\mathrm{R}\mathrm{1}}^{\mathrm{ref}}}}.\end{array}$$

The power gain compares the overall power extraction ${\stackrel{\mathrm{‾}}{P}}_{\mathrm{farm}}$ against the power ${\stackrel{\mathrm{‾}}{P}}^{\mathrm{ref}}$ extracted by the Betz-optimal reference case. The farm efficiency evaluates performance against a fictional farm where all turbines operate in the free-stream flow. In that case, the overall power extraction is compared against the average power ${\stackrel{\mathrm{‾}}{P}}_{\mathrm{R}\mathrm{1}}^{\mathrm{ref}}$ extracted by a Betz-optimal, front-row reference turbine. For the power computations, we only consider turbine operation after t_startup=300 s to take into account the startup time of the controllers. Error bars on η_P and η_farm are computed starting from t_startup=300 s using block bootstrapping with window length 600 s.

Computational times are measured in terms of the real-time factor

$$\begin{array}{}\text{(21)}& \mathrm{RT}={\displaystyle \frac{{\stackrel{\mathrm{‾}}{t}}_{\mathrm{wall}}}{T_A}},\end{array}$$

where ${\stackrel{\mathrm{‾}}{t}}_{\mathrm{wall}}$ is the wall time per optimization run, averaged over all optimization windows of the corresponding case. Error bars on the real-time factor are based on the deviations of the wall times over the different optimization windows. In the time-delayed MPC loop from Figs. 2b and 3, for real-time operation, all computations for a given receding-horizon window should be performed within a time interval of length T_A (the control update time). In the present study, we only consider the computational time for the optimization and omit details of the state estimation. In practice, depending on the method, estimation may take as long as the optimization of the controls, such that RT<0.5 is expected to be sufficiently fast for real-time operation. Note that the flow prediction subsequent to the state estimation (see Fig. 3) corresponds to one forward simulation on the coarse grid; computational time for the prediction is hence negligible compared to that of optimization (and possibly estimation).

4.3.1 Analysis of computational cost

First of all, as can be appreciated in Fig. 10, all real-time factors are bigger than one, ranging from 1.79 to 270. Three observations can be made: the real-time factor and hence the computational cost of the optimization increases if (a) the optimization horizon T increases, (b) the update time T_A is reduced or (c) the grid is refined. Case 10 on grid level 0 – with T=50 s, T_A=50 s and the lowest number of grid points of all simulations – is therefore the fastest and only a factor of 1.79 (on average) slower than real time. Conversely, case 0 on grid level 2 is the most challenging in terms of wall time, with a real-time factor of 270. In terms of RT, Fig. 10 reveals that, in general, the relative order of the cases remains unchanged when refining the grid.

Figure 10 Gain factor (left axis) and farm efficiency (right axis) versus real-time factor for the optimal control cases, including error bars. The Betz-optimal reference case and steady yaw control case are also indicated.

4.3.2 Analysis of power gains and farm efficiency

The power gains from Fig. 10 are more clearly summarized in Fig. 11. As can be expected for the fully aligned farm layout, the farm efficiency for the uncontrolled Betz-optimal reference case is relatively low at approximately 45 %. From Figs. 10 and 11, it can be seen that all optimal control cases improve on the uncontrolled reference, except for case 10 due to its short optimization horizon (T=50 s, T_A=50 s). The highest power gain is observed for case 4 (T=350 s, T_A=50 s) on grid level 1: η_P=1.51. Two main trends can be observed: (a) increasing the optimization horizon T increases the power extraction and (b) decreasing the control update time T_A increases the power extraction.

Figure 11 Gain factor (left axis) and farm efficiency (right axis) for the optimal control cases, including error bars. Betz-optimal reference case and steady yaw control case are also indicated (dotted lines).

Interestingly, the grid resolution has no clear impact on the power gains. In some cases, refining the grid increases the power extraction, and sometimes the power extraction decreases. However, the performance gains and losses that come from refining or coarsening the grid (for a given T and T_A) are only marginal compared to the effects of changing T and T_A and mostly fall within the bounds of the confidence interval. Overall, the influence of the receding-horizon parameters is hence much bigger than that of the grid resolution.

Figure 12 Time evolution of filtered thrust coefficients ${\widehat{C}}_{\mathrm{T}}^{\prime }$ for turbine column C1 for different optimal control cases for grid level 0. For cases 8–10 (not shown), the thrust coefficients do not deviate from the Betz-optimal value ${\widehat{C}}_{\mathrm{T}}^{\prime }=\mathrm{2}$.

4.4 Yaw and induction characteristics

Figures 12 and 13 illustrate the time evolution of respectively the filtered thrust coefficients ${\widehat{C}}_{\mathrm{T}}^{\prime }$ and yaw angles θ for the eight turbines in column C1 (see Fig. 6) for optimal control cases 1, 4, 5, 7, 8, 9 and 10 after optimization on grid level 0. We only show results for column C1, since the observations for the other turbine columns are similar. The thrust coefficients and yaw angles for the other grid levels are shown in Appendix C, but the trends observed there are similar to the ones for grid level 0. For comparison, the steady yaw angles from Sood and Meyers (2022) are shown in Appendix D.

Figure 13 Time evolution of yaw angles θ for turbine column C1 for different optimal control cases for grid level 0.

For cases 8–10, characterized by the shortest time horizons T, wakes cannot propagate from one row of turbines to the next row within the optimization window. In those cases, power is therefore maximized at the level of individual turbines: all turbines operate at the Betz limit ${\widehat{C}}_{\mathrm{T}}^{\prime }=\mathrm{2}$ while oscillating around the flow-aligned yaw angle of 0^∘. The magnitude of the oscillations increases for downstream turbines to account for the higher flow angles that exist in downstream regions of the wind farm due to the unsteadiness in the local flow. For these cases, the control update time T_A has no impact on the controls.

As the time horizon increases, significant yaw angles emerge for upstream turbines, with even some quasi-static yawing behavior for the cases with short update times. This is particularly evident for case 4 (T=350 s, T_A=50 s), where front-row turbines are immediately redirected to a yaw angle of ±30^∘. Downstream turbine rows 2–4 exhibit similar behavior but at lower misalignment angles and with more complex oscillations in response to the local unsteadiness of the flow. In contrast, the last few rows again operate around the unyawed position. Note that the misalignment angles of ±30^∘ for the front-row turbines match the value obtained using the static yaw controller from Sood and Meyers (2022). However, the main difference compared to the steady yaw case is the distinct yawing of downstream turbines (starting in row 2 already), resulting in significant gains. Case 7 behaves similar to case 4 in terms of yawing, but the quasi-static yawing is less pronounced due to the shorter optimization horizon. Also note that transitions of the yaw angle between −30 and 30^∘, as observed for front-row turbines for case 4 and case 7, are propagated in the mean flow to downstream turbines, as indicated by the red line in Fig. 13.

Increasing the control update time T_A (see, e.g., case 1 and 5) eliminates the quasi-static yawing behavior observed in the upstream turbines (as discussed above for case 4). This is because the controller redirects yawed turbines to the unyawed position at the end of the optimization window, as the corresponding wakes cannot propagate to the next turbine row anymore, rendering yaw control disadvantageous. This end-of-time effect is detrimental to the long-term power extraction, since it is merely an artifact of the finite optimization horizon. For case 4 (T_A=50 s, and analogously for case 7), the end-of-time effect is mitigated by the shorter control update time (T_A≪T). In contrast, for case 1 (T_A=350 s), upstream turbines are steered towards a yawed position at the start of every optimization window but then redirected at the end of the window due to the end-of-time effect. This process repeats for the consecutive windows. Case 5 exhibits similar behavior, but the time spent in the yawed position is shorter due to the shorter optimization window (a larger portion of the windows is affected by the end-of-time effect). Cases 2 and 3 and case 6 (not shown in Figs. 12 and 13) display behavior similar to cases 1 and 5 respectively, but the end-of-time effect is less pronounced due to the shorter control update time.

Overall, yaw control emerges as the dominant control mechanism for the TotalControl wind farm. Induction control is only used for longer optimization horizons (T≥250 s), resulting in minor deviations from the Betz-optimal thrust coefficient. Interestingly, for cases 1, 4, 5 and 7 and for the front-row turbines, it seems that the dips in C_T for the upstream turbines roughly coincide with the zero-crossings of the yaw angle of those turbines. This shows the connection between yaw and induction control: for the given setup, yaw control is the dominant control mechanism, and induction control is used as an additional control mechanism when there is no (instantaneous) yawing. Simulations (not shown here) suggest that yaw control only (i.e. disabling induction control) does not entail a significant performance reduction.

Note that the discussion in this section was limited to the effect of the receding-horizon parameters on the thrust coefficients and yaw angles. More general aspects of wind farm control – the curtailing of power of front-row turbines in favor of downstream turbines, the typical absence of yawing for downstream turbines, etc. – are similar to previous wind farm control studies such as the ones in Munters and Meyers (2017, 2018b); the reader is referred there for a more fundamental view on the physics of wind farm control.

4.5 Discussion

Figure 11 clearly shows that all optimal control cases outperform the Betz-optimal reference case in terms of power gain, except case 10 (due to its short optimization horizon T=50 s). If the time horizon is long enough, i.e. T>150 s (cases 1–7), we also obtain significant improvements over the steady yaw controller from Sood and Meyers (2022). This is a non-trivial observation, given the coarseness of the control models and the accompanying model mismatch compared to the fine-grid reference, especially in the context of the time-decoupled MPC loop including the prediction step. This suggests that large-scale structures (in space and time) in the wind farm boundary layer may already suffice to extract an adequate control signal. Short-term evolutions of the boundary layer are not accurately captured by the coarser control models, resulting in only limited gains (for cases 8 and 9) or even losses (case 10) compared to the Betz-optimal reference. Extending the time horizon allows the optimizer to tailor the controls based on the large-scale spatiotemporal structures, which are described sufficiently accurately by the control model, hence resulting in significant improvements. It must also be noted that increasing the horizon (either through T or T_A) increases the variability on the results, resulting in larger error bars on the power gains (compare, e.g., the error bars of cases 1–4 to those of cases 8–10). Furthermore, increasing T_A entails performance losses not only due to the end-of-time effect, but also due to the increased prediction error in the predictor.

The substantial power gains compared to the steady yaw controller from Sood and Meyers (2022) originate from the dynamic yaw steering throughout the whole farm in response to the turbulent inflow. As shown in Appendix D for the framework of Sood and Meyers (2022), only front-row turbines are effectively yawed to ±30^∘, whereas we also observe significant yawing in the downstream regions of the farm. On the one hand, for case 4 in particular, there is the quasi-static wake steering – with occasional turnovers from 30 to −30^∘ and vice versa depending on the inflow – that persists downstream in turbine rows R2–R4/R5. On the other hand, towards the end of the farm, we also observe dynamic yawing around a mean angle of approx. 0^∘ to optimally align the turbines to the turbulent inflow that has become increasingly unsteady due to the superimposed wakes from the upstream turbines. This also entails a significant power gain that cannot be captured with a steady yaw controller. Finally, we note that the proposed LES-based controller is able to account for secondary steering effects that are not included in the Bastankhah wake model that was used to optimize the yaw set points in Sood and Meyers (2022).

Interestingly, refining the optimization resolution does not significantly improve the results in terms of power gain and can even be disadvantageous in some cases (at least in the range of resolutions considered in this paper). This effect may be attributed to the model mismatch: upon refining, the controls are tuned to account for the additional small(er)-scale variations, but these (incremental) adjustments are not necessary optimal on the fine-grid reference. In Bauweraerts and Meyers (2019) (in the context of turbulent forecasting), it was shown that modeling errors even slightly decrease with grid coarsening, due to the decreasing subgrid-scale errors in the LES and the decreased effect of chaotic divergence of solution trajectories on coarser grids. In other words, it may be better to only optimize for the large scales than to also take into account smaller scales that may be inaccurately modeled. However, it must be noted that even the finest optimization resolution (level 2) is still more than 3 times coarser than the reference resolution. It is expected that further grid refinements would eventually improve on the coarse-grid results, when the actual small-scale variations are sufficiently accurately described by the control models. However, these kinds of simulations would be prohibitive due to excessive computational costs.

Finally, it must be noted that the coarse control models are incapable of capturing all the turbulent dynamics governing the optimal wind farm control problem. On the one hand, it seems that they can accurately model large-scale motions in the flow, such as the general deflection of the wake under yawed conditions and the gross behavior of the wakes. On the other hand, the shedding of vortex rings, which play a crucial role in enhancing wake mixing when dynamically controlling turbine thrust (Munters and Meyers, 2018a), cannot be represented on coarse grids. As such, the proposed methodology is less eligible for dynamic induction control, and yaw control emerges as the dominant mechanism for coarse-grid LES-based control, where in that case the gains mainly originate from dynamically steering away the wakes from downstream turbines, synchronized to the turbulent inflow.

5 Towards real-time optimal wind farm control

The discussion in Sect. 4.3 has revealed that, on the one hand, the optimization horizon T needs to be long enough to take into account wake interactions over subsequent turbine rows. On the other hand, the control update time needs to be short enough to discard end-of-time effects. Furthermore, it was shown that refining the grid resolution significantly increases the computational cost without a significant improvement in the performance of the controller in terms of power gain. By leveraging these insights, we now design a competitive (in terms of power gains) controller that is as close to real time as possible.

5.1 Case setup

For the receding-horizon parameters, we propose T=300 s and T_A=120 s. An optimization horizon of T=300 s allows the controller to account for wake interactions during the optimization. With T_A=120 s, the final 180 s of the optimization window are discarded, which is roughly the portion of the window affected by end-of-time effects. Based on the analysis in Sect. 4.3, longer update times are expected to reduce performance due to end-of-time effects. Although further decreasing T_A could potentially improve the power gains by mitigating the model mismatch, the benefits would be relatively insignificant compared to those that come from addressing the end-of-time effect and would hence result in an undue increase in the real-time factor. Regarding the convergence criteria, we consider three cases: case (a) that uses the same convergence criterion as in Sect. 3.1.1 (i.e. $({\mathcal{J}}^{k-\mathrm{1}}-{\mathcal{J}}^k)/{\mathcal{J}}^{k-\mathrm{1}}\le \mathrm{5}\times {\mathrm{10}}^{-\mathrm{6}}$), and cases (b) and (c) that additionally impose a maximum number of optimization iterations, respectively $N_{\mathrm{opt},\max }=\mathrm{10}$ and $N_{\mathrm{opt},\max }=\mathrm{5}$. Furthermore, given the limited contribution of induction control in the overall power gains, all thrust coefficients are kept constant to the Betz-optimal value $C_{\mathrm{T}}^{\prime }=\mathrm{2}$. By disabling induction control, the number of optimization variables decreases by a factor of 2, which may potentially speed up the convergence of the optimization problems.

To further minimize the real-time factor, simulations are performed on the coarsest grid level (level 0). We use the streamwise domain length from the cases with horizon T=250 s (see Table 2): L_x=11.52 km. Earlier results (not shown) suggest that this domain is still big enough to prevent recycling of the wakes into the front-row turbines. Furthermore, in comparison to Table 2, we apply an additional cropping of the optimization domain in the spanwise and vertical directions to further reduce computational costs, resulting in L_y=4.16 km and L_z=0.75 km. To keep a friction velocity of u_τ=0.28 m s⁻¹ for the new vertical dimension, the pressure gradient is adjusted to ${\partial }_x{p}_{\mathrm{\infty }}/\mathit{\rho }=\mathrm{1.0453}\times {\mathrm{10}}^{-\mathrm{4}}$ m s⁻². For the time step, we take Δt=4 s, which is the coarsest time step that still respects a CFL number of 0.8. The simulation setup is summarized in Table 6.

Table 6 Case description for the case from Sect. 5.

5.2 Results and discussion

The results for cases (a), (b) and (c) are summarized in Table 7. By tailoring the control model and through a sensible choice for the receding-horizon parameters as described above, for case (c) we achieve a real-time factor of 1.33 with only a minor decrease in power gain compared to cases (a) and (b) (1.31 versus 1.39 and 1.38 respectively). Apparently, just five L-BFGS-B iterations already suffice to extract an adequate control action.

Table 7 Gain factor versus real-time factor for the optimal control cases. Error bars are also shown in the brackets.

It is important to note that the setup used here was tailored to the fully aligned TCRWP, and therefore the results may not directly translate to other farm configurations. Nevertheless, the observations suggest that using coarse-grid LES for real-time wind farm control is a viable approach, if another order of magnitude (factor of 10) in computational or algorithmic speedup can be found (since, for practical wind farm control, the optimization process must be at least twice as fast as real time to perform state estimation within the computational window). However, as CPU computing power continues to advance, and since this first investigation already achieves near-real-time speed, a factor of 10 still remains within reach. Further speedups may also be attained through GPU-accelerated computing. Furthermore, it is worth mentioning that the SP-Wind code can also be enhanced. For instance, a new version of the code is currently being developed that employs 3D domain decomposition for the spatial parallelization (as opposed to the 2D domain decomposition used in this work). This upgrade will be necessary to handle even bigger optimization cases in real time. Very recently; Janssens and Meyers (2023) proposed a multiple shooting algorithm for large-scale optimal control cases, such as the ones considered here. The additional speedup due to the temporal parallelization in that case may potentially narrow the gap towards achieving actual, practical wind farm control in real time.

6 Conclusions

In the current paper, we investigated the influence of the grid resolution of the LES-based control model and receding-horizon parameters on the performance of the controller, both in terms of power gain and computational cost. To that end, we defined a set of optimal control cases with varying optimization horizons and control update times, as well as a fine-grid LES emulator model, applied to the TotalControl Reference Wind Power Plant. For each case, we defined three grid resolutions for the LES-based control model and performed a complete optimal control simulation on each of the grids.

Regarding the receding-horizon parameters, on the one hand, the results indicate that the optimization horizon should be long enough to take into account turbine–wake interactions over subsequent turbine rows. In that case, upstream turbines are misaligned to steer away the wakes from downstream turbines in a quasi-static way, resulting in significant power gains compared to the uncontrolled reference case and steady yaw control case. Downstream turbines are yawed as well, but to a lesser extent and in a more complex pattern, mostly in response to the local unsteadiness of the flow. On the other hand, the control update times should be short enough to mitigate end-of-time effects, i.e. to discard controls near the end of the optimization window that are affected by the finiteness of the optimization window. Taking this into account, optimal control case 4 (with T=350 s and T=50 s, i.e. the longest horizon and shortest update time) consistently produces the highest power gains – up to η_P=1.51 on grid level 1. Furthermore, it must be noted that all optimal control cases improve on the uncontrolled reference case in terms of power gain, except case 10 (T=50 s, T_A=50 s) due to the short optimization horizon. Moreover, cases 1–7 (where T≥250 s) also outperform a steady yaw control case that was optimized based on the Bastankhah wake model combined with a recursive wake merging method (see Sood and Meyers, 2022). Obviously, increasing the optimization horizon and decreasing the control update time both increase the real-time factor.

Somewhat surprisingly, the grid resolution has no significant impact on the performance of the controllers, at least not in the range of resolutions considered here. Sometimes refining the grid resolution results in higher gains, whereas sometimes the gains decrease. This may be attributed to the existence of many local optima in the wind farm optimization problem, as well as the model uncertainty, since the finest grid in this work was still more than 3 times coarser than the reference simulation. It is expected that finer grids (with resolutions close to that of the reference) would eventually improve on the coarse models investigated here, as soon as the small scales are represented sufficiently accurately to take them into account in the control action. However, optimizations on these kinds of resolution would be prohibitive due to computational cost and outweigh the potential gains originating from the decreased model mismatch. It must also be noted that we did not investigate the influence of the convergence criterion; therefore, it is possible that somewhat better results can be obtained on each of the grid levels by tweaking the stopping criterion.

In terms of power gains, the coarse-grid control models investigated in the current paper perform surprisingly well, with gains up to and above 40 % if the time horizon is sufficiently long and the update time sufficiently short. This is a powerful result, given the model uncertainty due to the additional prediction step in the time-delayed MPC loop. Regarding the complex physics underlying wind farm control, the results suggest that the large-scale spatial and temporal structures in the wind farm boundary layer (i.e. the ones that are accurately represented by the coarse LES models) suffice to extract an efficient yaw controller. From the viewpoint of real-time LES-based optimal control, this means that high-performance controllers can be obtained at only a fraction of the computational cost by coarsening the grid resolution, potentially bridging the gap between theoretical studies based on LES and practical, real-time wind farm control. Via a proper choice of the receding-horizon parameters and optimal control domain, we achieved a gain of 31 % with a real-time factor of 1.33, i.e. only 1.33 times slower than real time. Even better control signals (in terms of the resulting power gains) may be obtained through an enhanced wake mixing by better exploiting dynamical induction control (Munters and Meyers, 2018a) or via the helical wake structures originating from individual pitch control (Frederik et al., 2020). However, this would require a much finer numerical grid, which currently does not allow for a real-time implementation in LES.

It must be noted that, upon the restriction from the actual wind farm (the fine-grid LES wind farm model) to the coarse resolution of the prediction and control model by the restriction operator, we still assume the entire flow field from the reference simulation is available. In practice, this is not the case; the flow field is only available in the form of discrete measurements, for example, lidar measurements or SCADA data. Consequently, for a practical controller, the flow field must first be reconstructed from these measurements in a state estimator, for example, using Kalman filtering or 4D variational data assimilation (see, e.g., Bauweraerts and Meyers, 2021). Future work should therefore focus on the design of an efficient state estimator tailored to wind farm control, which can then be incorporated into the LES-based wind farm control methodology proposed here. As shown in previous work from Bauweraerts and Meyers (2019), LES-based state estimation may potentially benefit from grid coarsening as well. However, in that case, the question still remains of whether combined LES-based state estimation and control can be fast enough for real-time applications, as well as how the additional reconstruction errors originating from the state estimation will affect the performance of the controllers. This would be the next step towards actual, practical wind farm control.

We also admit that the considered test case (TCRWP) is quite susceptible to high power gains due to the fully aligned turbine configuration. Applying the proposed wind farm controller to other wind farm layouts could be a topic for future research, although it is expected that the controller will still be able to produce competitive gains. To increase the credibility of the proposed control strategy in general, another interesting direction would be to use more accurate turbine models in the wind farm emulator model (e.g., actuator-sector models (ASMs) or actuator-line models (ALMs)), instead of the very basic non-rotating, actuator-disc model used here. In the first place, one could investigate the effect of the additional model mismatch when using an ADM control model in the coarse LES in combination with an ASM or ALM reference model.

Appendix A: Lookup tables for power correction

This appendix lists the correction factors M from Eq. (10) for the disc-averaged velocity for different values of the thrust coefficient ${\widehat{C}}_{\mathrm{T}}^{\prime }$ and yaw angle θ for the three grid resolutions considered in the present paper. The corrections are chosen in such a way that the corrected disc-averaged velocity on the coarse grid equals the disc-averaged velocity on the reference grid for a given thrust coefficient and yaw angle for uniform inflow. The resulting lookup values are summarized in Tabes A1–A3 for grid level 0, grid level 1 and grid level 2 respectively. These tables are constructed based on uniform inflow simulations with U_∞=8 m s⁻¹, prescribed using a fringe region spanning the final 20 % of the simulation domain. All simulations are conducted for a single DTU 10 MW reference turbine using the same domain size $L_x\times L_y\times L_z=\mathrm{26.92}$ D × 13.46 D × 8.41 D, where D=178.3 m is the rotor diameter. For the resolution of the reference, we use $\mathrm{\Delta }x\times \mathrm{\Delta }y\times \mathrm{\Delta }z=\mathrm{13.33}\times \mathrm{13.33}\times \mathrm{6.67}$ m³ as prescribed by Table 1. For a given ${\widehat{C}}_{\mathrm{T}}^{\prime }$ and θ, the corresponding correction factor M is computed by (linearly) interpolating the lookup tables.

Finally, we also briefly examine the performance of the lookup table approach for turbulent inflow. To that end, as an example, we take the setup of case 1 from Tables 2 and 3 on grid level 0. For ${\widehat{C}}_{\mathrm{T}}^{\prime }=\mathrm{2.0}$ and θ=0^∘, Fig. A1 shows the corrected (i.e. using the lookup table approach) and uncorrected power predictions on the coarse grids compared to the fine-grid reference for a horizon of 300 s. With the lookup table approach, the error on the average power prediction is decreased by 7 %, i.e. from an overestimation of 21 % for the uncorrected prediction to an underestimation of 14 % for the corrected one.

Table A1 Lookup table for grid level 0 ($\mathrm{\Delta }x\times \mathrm{\Delta }y\times \mathrm{\Delta }z=\mathrm{80}\times \mathrm{80}\times \mathrm{50}$ m³): correction factor M from Eq. (10) for different thrust coefficient ${\widehat{C}}_{\mathrm{T}}^{\prime }$ and yaw angle θ.

Table A2 Lookup table for grid level 1 ($\mathrm{\Delta }x\times \mathrm{\Delta }y\times \mathrm{\Delta }z=\mathrm{60}\times \mathrm{60}\times \mathrm{37.5}$ m³): correction factor M from Eq. (10) for different thrust coefficient ${\widehat{C}}_{\mathrm{T}}^{\prime }$ and yaw angle θ.

Table A3 Lookup table for grid level 2 ($\mathrm{\Delta }x\times \mathrm{\Delta }y\times \mathrm{\Delta }z=\mathrm{40}\times \mathrm{40}\times \mathrm{25}$ m³): correction factor M from Eq. (10) for different thrust coefficient ${\widehat{C}}_{\mathrm{T}}^{\prime }$ and yaw angle θ.

Figure A1 Power predictions for turbine column C1 without correction (blue) and using the lookup table approach (red) compared to the fine-grid reference for ${\widehat{C}}_{\mathrm{T}}^{\prime }=\mathrm{2.0}$ and θ=0^∘.

Appendix B: Derivation and verification of the temporally discrete adjoint method and gradient

In this appendix, we formulate the temporally discrete adjoint method and derive Eq. (18) for the cost functional gradients. These gradients are also validated via finite differences. The derivation and notation are similar to the one from Yilmaz and Meyers (2019); differences are explicitly formulated throughout the text.

B1 Discretization of the optimal control problem

The discrete adjoint method first discretizes and then linearizes the state equations; next it formulates the discrete adjoint of the linear system (Giles and Pierce, 2000). Using an explicit fourth-order Runge–Kutta discretization with N time steps, the discretization of the optimization problem (Eqs. 1–7) can symbolically be written as

$$\begin{array}{}\text{(B1)}& {\displaystyle }{\min }_{\mathit{\phi },\mathit{q}}{\stackrel{\mathrm{̃}}{\mathcal{J}}}^N=\sum _{n=\mathrm{1}}^N{I}^n\mathrm{s}.\mathrm{t}.{\mathit{q}}_i^n={\mathit{q}}^n+\mathrm{\Delta }t\sum _{j=\mathrm{1}}^{\mathrm{4}}{a}_{ij}\mathit{Y}\left({\mathit{q}}_j^n,{\mathit{\phi }}^n\right)i=\mathrm{1}\mathrm{\dots }\mathrm{4},\text{(B2)}& {\displaystyle }{\mathit{q}}^{n+\mathrm{1}}={\mathit{q}}^n+\mathrm{\Delta }t\sum _{i=\mathrm{1}}^{\mathrm{4}}{b}_i\mathit{Y}\left({\mathit{q}}_i^n,{\mathit{\phi }}^n\right),\end{array}$$

where ${\mathit{q}}^n=[{\mathit{u}}^n,p^n,{\widehat{C}}_{\mathrm{T}}^{{}^{\prime }n},{\mathit{\theta }}^n]$ and ${\mathit{\phi }}^n=[C_{\mathrm{T}}^{{}^{\prime }n},{\mathit{\omega }}^n]$ are respectively the state variables and controls at time instant tⁿ=nΔt. In the Runge–Kutta stage equations, Y denotes the right side of the governing equations consisting of the Navier–Stokes equations, thrust coefficient filter equation and yaw rate equation. I_n denotes the discretized objective functional. Below, we continue the derivation for a general explicit (j<i) fourth-order Runge–Kutta scheme; the simulations in SP-Wind are carried out using classical Runge–Kutta 4, for which the nonzero Butcher tableau coefficients are $a_{\mathrm{21}}=a_{\mathrm{32}}=\mathrm{1}/\mathrm{2}$, a₄₃=1, $b_{\mathrm{1}}=b_{\mathrm{4}}=\mathrm{1}/\mathrm{6}$ and $b_{\mathrm{2}}=b_{\mathrm{3}}=\mathrm{1}/\mathrm{3}$.

As opposed to the basic first-order discretization from Yilmaz and Meyers (2019), in this work the intermediate Runge–Kutta stages are also used in the discretization of the cost functional. In order to achieve this, the continuous cost functional (Eq. 1), here symbolically written as

$$\stackrel{\mathrm{̃}}{\mathcal{J}}(\mathit{\phi },\mathit{q}(\mathit{\phi }\left)\right)=\underset{\mathrm{0}}{\overset{T}{\int }}J(\mathit{\phi },\mathit{q}(\mathit{\phi }\left)\right)\mathrm{d}t,$$

is first rewritten in the form of an ordinary differential equation (ODE),

$$\begin{array}{}\text{(B3)}& {\displaystyle \frac{\mathrm{d}{\mathcal{J}}_t\left(t\right)}{\mathrm{d}t}}=J\left(\mathit{\phi }\right(t),\mathit{q}(t\left)\right){\mathcal{J}}_t\left(\mathrm{0}\right)=\mathrm{0},\end{array}$$

where $\stackrel{\mathrm{̃}}{\mathcal{J}}(\mathit{\phi },\mathit{q}(\mathit{\phi }\left)\right)={\mathcal{J}}_t\left(T\right)$. The ODE in Eq. (B3) is then discretized using Runge–Kutta 4, yielding a new expression for the discretized cost function that is of fourth-order accuracy:

$$\begin{array}{}\text{(B4)}& {\stackrel{\mathrm{̃}}{\mathcal{J}}}^N=\sum _{n=\mathrm{1}}^N{I}^n=\sum _{n=\mathrm{1}}^N\left(\mathrm{\Delta }t\sum _{i=\mathrm{1}}^{\mathrm{4}}{b}_{i}J\left({\mathit{\phi }}^n,{\mathit{q}}_i^n\right)\right).\end{array}$$

Note that compared to Yilmaz and Meyers (2019), the control variables φⁿ in the Runge–Kutta discretization in Eqs. (B1)–(B2) are kept constant over the Runge–Kutta stages, since (in practice) controls are kept constant over the control time step (which we assume here is equal to the discretization time step Δt).

B2 Linearization of the state equation

Analogously to Yilmaz and Meyers (2019), the discretized state equations in Eqs. (B1)–(B2) can be linearized, resulting in a linear system for every time step n:

$$\begin{array}{}\text{(B5)}& {\displaystyle }{\mathbf{K}}^n{\mathbf{L}}^n={\mathbf{M}}^n,\text{(B6)}& {\displaystyle }\begin{array}{rl}& \mathbf{K}=\\ & \left[\begin{array}{cccccc}-\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{1}& -\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{1}& a_{\mathrm{21}}\mathrm{\Delta }t{\mathbf{Y}}_{\mathit{q}}\left({\mathit{q}}_{\mathrm{1}}^n,{\mathit{\phi }}^n\right)& -\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{1}& a_{\mathrm{31}}\mathrm{\Delta }t{\mathbf{Y}}_{\mathit{q}}\left({\mathit{q}}_{\mathrm{1}}^n,{\mathit{\phi }}^n\right)& a_{\mathrm{32}}\mathrm{\Delta }t{\mathbf{Y}}_{\mathit{q}}\left({\mathit{q}}_{\mathrm{2}}^n,{\mathit{\phi }}^n\right)& -\mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{1}& a_{\mathrm{41}}\mathrm{\Delta }t{\mathbf{Y}}_{\mathit{q}}\left({\mathit{q}}_{\mathrm{1}}^n,{\mathit{\phi }}^n\right)& a_{\mathrm{42}}\mathrm{\Delta }t{\mathbf{Y}}_{\mathit{q}}\left({\mathit{q}}_{\mathrm{2}}^n,{\mathit{\phi }}^n\right)& a_{\mathrm{43}}\mathrm{\Delta }t{\mathbf{Y}}_{\mathit{q}}\left({\mathit{q}}_{\mathrm{3}}^n,{\mathit{\phi }}^n\right)& -\mathrm{1}& \mathrm{0}\\ \mathrm{1}& b_{\mathrm{1}}\mathrm{\Delta }t{\mathbf{Y}}_{\mathit{q}}\left({\mathit{q}}_{\mathrm{1}}^n,{\mathit{\phi }}^n\right)& b_{\mathrm{2}}\mathrm{\Delta }t{\mathbf{Y}}_{\mathit{q}}\left({\mathit{q}}_{\mathrm{2}}^n,{\mathit{\phi }}^n\right)& b_{\mathrm{3}}\mathrm{\Delta }t{\mathbf{Y}}_{\mathit{q}}\left({\mathit{q}}_{\mathrm{3}}^n,{\mathit{\phi }}^n\right)& b_{\mathrm{4}}\mathrm{\Delta }t{\mathbf{Y}}_{\mathit{q}}\left({\mathit{q}}_{\mathrm{4}}^n,{\mathit{\phi }}^n\right)& -\mathrm{1}\end{array}\right],\end{array}\text{(B7)}& {\displaystyle }\begin{array}{rl}& \mathbf{L}=\left(\begin{array}{c}\mathit{\delta }{\mathit{q}}^n\\ \mathit{\delta }{\mathit{q}}_{\mathrm{1}}^n\\ \mathit{\delta }{\mathit{q}}_{\mathrm{2}}^n\\ \mathit{\delta }{\mathit{q}}_{\mathrm{3}}^n\\ \mathit{\delta }{\mathit{q}}_{\mathrm{4}}^n\\ \mathit{\delta }{\mathit{q}}^{n+\mathrm{1}}\end{array}\right)\\ & \mathbf{M}=-\left(\begin{array}{c}\mathrm{0}\\ \mathrm{0}\\ {\sum }_{j=\mathrm{1}}^{\mathrm{1}}{a}_{\mathrm{2}j}\mathrm{\Delta }t{\mathbf{Y}}_{\mathit{\phi }}\left({\mathit{q}}_j^n,{\mathit{\phi }}^n\right)\mathit{\delta }{\mathit{\phi }}^n\\ {\sum }_{j=\mathrm{1}}^{\mathrm{2}}{a}_{\mathrm{3}j}\mathrm{\Delta }t{\mathbf{Y}}_{\mathit{\phi }}\left({\mathit{q}}_j^n,{\mathit{\phi }}^n\right)\mathit{\delta }{\mathit{\phi }}^n\\ {\sum }_{j=\mathrm{1}}^{\mathrm{3}}{a}_{\mathrm{4}j}\mathrm{\Delta }t{\mathbf{Y}}_{\mathit{\phi }}\left({\mathit{q}}_j^n,{\mathit{\phi }}^n\right)\mathit{\delta }{\mathit{\phi }}^n\\ {\sum }_{i=\mathrm{1}}^{\mathrm{4}}{b}_i\mathrm{\Delta }t{\mathbf{Y}}_{\mathit{\phi }}\left({\mathit{q}}_i^n,{\mathit{\phi }}^n\right)\mathit{\delta }{\mathit{\phi }}^n\end{array}\right).\end{array}\end{array}$$

B3 Adjoint equations

Denote the adjoint vector for time step n by ${\mathit{N}}^n=\left(\begin{array}{c}{\mathit{q}}^{{*}^n}{\mathit{q}}_{\mathrm{1}}^{{*}^n}{\mathit{q}}_{\mathrm{2}}^{{*}^n}{\mathit{q}}_{\mathrm{3}}^{{*}^n}{\mathit{q}}_{\mathrm{4}}^{{*}^n}{\mathit{q}}^{{*}^{n+\mathrm{1}}}\end{array}\right)$, where ${\mathit{q}}^{{*}^n}=[{\mathit{\xi }}^n,{\mathit{\pi }}^n,{\mathit{\sigma }}^n,{\mathit{\eta }}^n]$ are the adjoint variables associated with the state vector qⁿ. Again in the notation of Yilmaz and Meyers (2019), we define the adjoint variables as the solution of the adjoint equation:

$$\begin{array}{}\text{(B8)}& {\left({\mathit{N}}^n{\mathbf{K}}^n\right)}^T={\mathit{I}}_{\mathit{q}}^n={\left(\begin{array}{cccccc}{I}_{{\mathit{q}}^n}^n& I_{{\mathit{q}}_{\mathrm{1}}^n}^n& I_{{\mathit{q}}_{\mathrm{2}}^n}^n& I_{{\mathit{q}}_{\mathrm{3}}^n}^n& I_{{\mathit{q}}_{\mathrm{4}}^n}^n& I_{{\mathit{q}}^{n+\mathrm{1}}}^n\end{array}\right)}^T.\end{array}$$

As opposed to Yilmaz and Meyers (2019), the partial derivatives of the cost function (Eq. B4) with respect to the Runge–Kutta stages are now nonzero (due to the more elaborate discretization):

$$\begin{array}{}\text{(B9)}& {\displaystyle }& {\displaystyle }{I}_{{\mathit{q}}_i^n}^n=b_i\mathrm{\Delta }t{J}_{\mathit{q}}\left({\mathit{\phi }}^n,{\mathit{q}}_i^n\right)i=\mathrm{1}\mathrm{\dots }\mathrm{4},\text{(B10)}& {\displaystyle }& {\displaystyle }{I}_{{\mathit{q}}^n}^n=I_{{\mathit{q}}^{n+\mathrm{1}}}^n=\mathbf{0}.& {\displaystyle }\end{array}$$

Based on Eq. (B8), the temporally discrete adjoint equations become

$$\begin{array}{}\text{(B11)}& \begin{array}{rl}{\mathit{q}}_i^{{*}^n}& =b_i\mathrm{\Delta }t\left({\mathbf{Y}}_{\mathit{q}}^T\left({\mathit{q}}_i^n,{\mathit{\phi }}^n\right){\mathit{q}}^{{*}^{n+\mathrm{1}}}-J_{\mathit{q}}\left({\mathit{q}}_i^n,{\mathit{\phi }}^n\right)\right)\\ & +\mathrm{\Delta }t\sum _{j=i+\mathrm{1}}^{\mathrm{4}}{a}_{ji}{\mathbf{Y}}_{\mathit{q}}^T\left({\mathit{q}}_i^n,{\mathit{\phi }}^n\right){\mathit{q}}_j^{{*}^n}i=\mathrm{1}\mathrm{\dots }\mathrm{4},\end{array}\end{array}$$

$$\begin{array}{}\text{(B12)}& {\mathit{q}}^{{*}^n}={\mathit{q}}^{{*}^{n+\mathrm{1}}}+\sum _{i=\mathrm{1}}^{\mathrm{4}}{\mathit{q}}_i^{{*}^n}.\end{array}$$

These equations run backward in time, starting from terminal conditions ${\mathit{q}}^{{*}^N}=I_{{\mathit{q}}^N}^N=\mathbf{0}$. Since each Runge–Kutta stage contributes to the (discretized) cost function, the cost function gradient J_q now appears as a source term in each equation. Note that, structurally, the temporally discrete adjoint equations differ from the forward Runge–Kutta scheme, since the method is not self-adjoint.

Due to the linearity of the adjoint equations, the stage equations from Eq. (B12) can be rewritten as follows:

$$\begin{array}{}\text{(B13)}& \begin{array}{rl}{\mathit{q}}_i^{{*}^n}& =\mathrm{\Delta }t{\mathbf{Y}}_{\mathit{q}}^T\left({\mathit{q}}_i^n,{\mathit{\phi }}^n\right)\left(b_i{\mathit{q}}^{{*}^{n+\mathrm{1}}}+\sum _{j=i+\mathrm{1}}^{\mathrm{4}}{a}_{ji}{\mathit{q}}_j^{{*}^n}\right)\\ & -b_i\mathrm{\Delta }t{J}_{\mathit{q}}\left({\mathit{q}}_i^n,\mathit{\phi }\right)i=\mathrm{1}\mathrm{\dots }\mathrm{4}.\end{array}\end{array}$$

To ease the notation in the remainder of the derivation, we introduce the auxiliary variable ${\widehat{\mathit{q}}}_i^{{*}^n}=[{\widehat{\mathit{\xi }}}_i^n,{\widehat{\mathit{\pi }}}_i^n,{\widehat{\mathit{\sigma }}}_i^n,{\widehat{\mathit{\eta }}}_i^n]$,

$$\begin{array}{}\text{(B14)}& {\widehat{\mathit{q}}}_i^{{*}^n}=b_i{\mathit{q}}^{{*}^{n+\mathrm{1}}}+\sum _{j=i+\mathrm{1}}^{\mathrm{4}}{a}_{ji}{\mathit{q}}_j^{{*}^n}i=\mathrm{1}\mathrm{\dots }\mathrm{4},\end{array}$$

such that we arrive at the following adjoint equations,

$$\begin{array}{}\text{(B15)}& {\mathit{q}}_i^{{*}^n}=\mathrm{\Delta }t{\mathbf{Y}}_{\mathit{q}}^T\left({\mathit{q}}_i^n,{\mathit{\phi }}^n\right){\widehat{\mathit{q}}}_i^{{*}^n}-b_i\mathrm{\Delta }t{J}_{\mathit{q}}\left({\mathit{q}}_i^n,\mathit{\phi }\right)i=\mathrm{1}\mathrm{\dots }\mathrm{4},\end{array}$$

which can be solved backwards, starting from i=4.

For the wind farm power optimization problem at hand, the adjoint Jacobians ${\mathbf{Y}}_{\mathit{q}}^T$ and J_q are exactly equal to the ones used in Munters and Meyers (2018b), where the continuous adjoint approach was used. To apply the temporally discrete adjoint equations to the wind farm control problem (Eqs. 1–7), we can therefore recycle these expressions; however, here we must also include the extra dependency of the lookup correction factor in Eq. (10) on the thrust coefficient and yaw angle. As such, we arrive at the following temporally discrete adjoint Runge–Kutta 4 scheme for the Navier–Stokes equations (where here V and $\stackrel{\mathrm{‾}}{V}$ denote the uncorrected and corrected disc-averaged velocities; see Eq. 10):

$$\begin{array}{}\text{(B16)}& {\displaystyle }\begin{array}{rl}{\mathit{\xi }}_i^n/\mathrm{\Delta }t& =\left(-{\left(\mathrm{\nabla }{\mathit{u}}_i^n\right)}^T+\left({\mathit{u}}_i^n\cdot \mathrm{\nabla }\right)\right){\widehat{\mathit{\xi }}}_i^n-\mathrm{\nabla }{\mathit{\pi }}_i^n/\mathit{\rho }\\ & -\mathrm{\nabla }\cdot {\mathit{\tau }}_{\mathrm{sgs}}^{*}\left({\mathit{u}}_i^n,{\widehat{\mathit{\xi }}}_i^n\right)+\sum _{m=\mathrm{1}}^{N_t}{\mathit{f}}_{m,i}^{{*}^n}i=\mathrm{1}\mathrm{\dots }\mathrm{4},\end{array}\text{(B17)}& {\displaystyle }\begin{array}{rl}{\mathrm{\nabla }}^{\mathrm{2}}{\mathit{\pi }}_i^n/\mathit{\rho }& =\mathrm{\nabla }\cdot \left[\left(-{\left(\mathrm{\nabla }{\mathit{u}}_i^n\right)}^T+\left({\mathit{u}}_i^n\cdot \mathrm{\nabla }\right)\right){\widehat{\mathit{\xi }}}_i^n\right.\\ & \left.-\mathrm{\nabla }\cdot {\mathit{\tau }}_{\mathrm{sgs}}^{*}\left({\mathit{u}}_i^n,{\widehat{\mathit{\xi }}}_i^n\right)+\sum _{m=\mathrm{1}}^{N_t}{\mathit{f}}_{m,i}^{{*}^n}\right]\\ & i=\mathrm{1}\mathrm{\dots }\mathrm{4},\end{array}\end{array}$$

$$\begin{array}{}\text{(B18)}& {\mathit{\xi }}^n={\mathit{\xi }}^{n+\mathrm{1}}+\sum _{i=\mathrm{1}}^{\mathrm{4}}{\mathit{\xi }}_i^n,\end{array}$$

where Eq. (B17) represents the Poisson equation obtained from the time-discrete adjoint momentum equations (Eq. B16) for the Runge–Kutta stages. The Poisson equation is solved at every time step and for every stage using a direct method. The adjoint wind farm force in Runge–Kutta stage i in Eqs. (B16)–(B17) then becomes (with m the turbine number)

$$\begin{array}{}\text{(B19)}& \begin{array}{rl}{\mathit{f}}_{m,i}^{{*}^n}=& {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\widehat{C}}_{{\mathrm{T}}_m}^{\prime n}M\left({\widehat{C}}_{{\mathrm{T}}_m}^{\prime n},{\mathit{\theta }}_m^n\right){\stackrel{\mathrm{‾}}{V}}_{m,i}^n\\ & \left(\mathrm{3}{b}_i{\stackrel{\mathrm{‾}}{V}}_{m,i}^n-\mathrm{2}{\widehat{X}}_{m,i}^n\right){\mathcal{R}}_m{\mathit{e}}_{⟂,m}^{n}i=\mathrm{1}\mathrm{\dots }\mathrm{4},\end{array}\end{array}$$

where ${\widehat{X}}_{m,i}^n=\frac{\mathrm{1}}{A_m}{\int }_{\mathrm{\Omega }}{\mathcal{R}}_m{\widehat{\mathit{\xi }}}_i^n\cdot {\mathit{e}}_{⟂,m}^n\mathrm{d}\mathit{x}$ denotes the disc-averaged adjoint velocity based on ${\widehat{\mathit{\xi }}}_i^n$, and ${\stackrel{\mathrm{‾}}{V}}_{m,i}^n$ denotes the corrected disc-averaged forward velocity based on ${\mathit{u}}_i^n$ (see Eq. 10).

Analogously, for the thrust coefficient filter and yaw rate equations, we get (i and j subscripts denote Runge–Kutta stages; m subscripts denote turbine numbers)

$$\begin{array}{}\text{(B20)}& {\displaystyle }\begin{array}{rl}{\mathit{\sigma }}_{m,i}^n/\mathrm{\Delta }t& ={\displaystyle \frac{\mathrm{1}}{\mathit{\tau }}}\left({\widehat{\mathit{\sigma }}}_{m,i}^n-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\stackrel{\mathrm{‾}}{V}}_{m,i}^{n^{\mathrm{2}}}\left(b_i{\stackrel{\mathrm{‾}}{V}}_{m,i}^n-{\widehat{X}}_{m,i}^n\right)A_m\right.\\ & -{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\widehat{C}}_{{\mathrm{T}}_m}^{\prime n}{\stackrel{\mathrm{‾}}{V}}_{m,i}^n{V}_{m,i}^n{\displaystyle \frac{\partial M\left({\widehat{C}}_{{\mathrm{T}}_m}^{\prime n},{\mathit{\theta }}_m^n\right)}{\partial {\widehat{C}}_{\mathrm{T}}^{\prime }}}\\ & \left.\left(\mathrm{3}{b}_i{\stackrel{\mathrm{‾}}{V}}_{m,i}^n-\mathrm{2}{\widehat{X}}_{m,i}^n\right)A_m\right)i=\mathrm{1}\mathrm{\dots }\mathrm{4},\end{array}\text{(B21)}& {\displaystyle }{\mathit{\sigma }}_m^n={\mathit{\sigma }}_m^{n+\mathrm{1}}+\sum _{i=\mathrm{1}}^{\mathrm{4}}{\mathit{\sigma }}_{m,i}^n,\end{array}$$

and

$$\begin{array}{}\text{(B22)}& {\displaystyle }\begin{array}{rl}{\mathit{\eta }}_{m,i}^n/\mathrm{\Delta }t=& -{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\widehat{C}}_{{\mathrm{T}}_m}^{\prime n}{\stackrel{\mathrm{‾}}{V}}_{m,i}^n\left[\underset{\mathrm{\Omega }}{\int }\left(\left(\mathrm{3}{b}_i{\stackrel{\mathrm{‾}}{V}}_{m,i}^n-\mathrm{2}{\widehat{X}}_{m,i}^n\right)\right.\right.\\ & \left.\left.M{\mathit{u}}_i^n-{\stackrel{\mathrm{‾}}{V}}_{m,i}^n\mathit{\xi }\right)\cdot \left({\mathit{e}}_{\left|\right|,m}{\mathcal{R}}_m+{\mathit{e}}_{⟂,m}{\mathcal{D}}_m\right)\mathrm{d}\mathit{x}\right]\\ & -{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\widehat{C}}_{{\mathrm{T}}_m}^{\prime n}{\stackrel{\mathrm{‾}}{V}}_{m,i}^n{V}_{m,i}^n{\displaystyle \frac{\partial M\left({\widehat{C}}_{{\mathrm{T}}_m}^{\prime n},{\mathit{\theta }}_m^n\right)}{\partial \mathit{\theta }}}\\ & \left(\mathrm{3}{b}_i{\stackrel{\mathrm{‾}}{V}}_{m,i}^n-\mathrm{2}{\widehat{X}}_{m,i}^n\right)A_{m}i=\mathrm{1}\mathrm{\dots }\mathrm{4},\end{array}\text{(B23)}& {\displaystyle }{\mathit{\eta }}_m^n={\mathit{\eta }}_m^{n+\mathrm{1}}+\sum _{i=\mathrm{1}}^{\mathrm{4}}{\mathit{\eta }}_{m,i}^n.\end{array}$$

For a more detailed explanation on all the terms and equations, as well as the derivation of the adjoints, the reader is referred to Goit and Meyers (2015) and Munters and Meyers (2017, 2018b).

B4 Adjoint gradients

The total variation of the cost functional follows from the chain rule

$$\begin{array}{}\text{(B24)}& \mathit{\delta }{\stackrel{\mathrm{̃}}{\mathcal{J}}}^N=\sum _{n=\mathrm{1}}^N\mathit{\delta }{I}^n=\sum _{n=\mathrm{1}}^N{\left({\mathbf{L}}^n\right)}^T{\mathit{I}}_{\mathit{q}}^n+I_{\mathit{\phi }}\mathit{\delta }{\mathit{\phi }}^n.\end{array}$$

We can now plug in the adjoint equation (Eq. B8) and then the linearization (Eq. B5), such that

$$\begin{array}{}\text{(B25)}& \begin{array}{rl}\mathit{\delta }{\stackrel{\mathrm{̃}}{\mathcal{J}}}^N& =\sum _{n=\mathrm{1}}^N\mathit{\delta }{I}^n=\sum _{n=\mathrm{1}}^N{\left({\mathbf{L}}^n\right)}^T{\left({\mathbf{K}}^n\right)}^T{\left({\mathit{N}}^n\right)}^T+I_{\mathit{\phi }}\mathit{\delta }{\mathit{\phi }}^n\\ & =\sum _{n=\mathrm{1}}^N{\left({\mathbf{M}}^n\right)}^T{\left({\mathit{N}}^n\right)}^T+I_{\mathit{\phi }}\mathit{\delta }{\mathit{\phi }}^n.\end{array}\end{array}$$

The derived expression hence does not require the forward sensitivity matrix L, which would otherwise have to be determined for every control perturbation, making the computation very expensive. From Eq. (B25), we can derive the gradient in any control direction. In practice, the L-BFGS-B library needs the gradient for every control variable δφⁿ, which amounts to

$$\begin{array}{}\text{(B26)}& \begin{array}{rl}{\displaystyle \frac{\partial {\stackrel{\mathrm{̃}}{\mathcal{J}}}^N}{\partial {\mathit{\phi }}^n}}& =-\mathrm{\Delta }t\sum _{i=\mathrm{1}}^{\mathrm{4}}\left(\sum _{j=\mathrm{1}}^{i-\mathrm{1}}\left(a_{ij}{\mathbf{Y}}_{\mathit{\phi }}^T\left({\mathit{q}}_j^n,{\mathit{\phi }}^n\right){\mathit{q}}_i^{{*}^n}\right)\right.\\ & \left.+b_i\left({\mathbf{Y}}_{\mathit{\phi }}^T\left({\mathit{q}}_i^n,{\mathit{\phi }}^n\right){\mathit{q}}^{{*}^{n+\mathrm{1}}}-J_{\mathit{\phi }}\left({\mathit{q}}_i^n,{\mathit{\phi }}^n\right)\right)\right).\end{array}\end{array}$$

Note that Eq. (B26) is the exact gradient of the discretized objective function in Eq. (B2), which converges to the gradient of the continuous problem in the limit of Δt→0 (Giles and Pierce, 2000).

Applied to the wind farm control problem (Eqs. 1–7), we arrive at the following simplified expressions for the gradient of the wind farm power objective function with respect to the thrust coefficients and yaw rates (i subscripts denote Runge–Kutta stages; m subscripts denote turbine numbers):

$$\begin{array}{}\text{(B27)}& \begin{array}{rl}{\mathrm{\nabla }}_{{\mathit{\phi }}_m^n}{\stackrel{\mathrm{̃}}{\mathcal{J}}}^N& =\left(\begin{array}{c}\frac{\partial {\stackrel{\mathrm{̃}}{\mathcal{J}}}^N}{{\widehat{C}}_{{\mathrm{T}}_m}^{\prime n}}\\ \frac{\partial {\stackrel{\mathrm{̃}}{\mathcal{J}}}^N}{\partial {\mathit{\omega }}_m^n}\end{array}\right)\\ & =\left(\begin{array}{c}-\mathrm{\Delta }t{\sum }_{i=\mathrm{1}}^{\mathrm{4}}\left(b_i{\mathit{\sigma }}_m^{n+\mathrm{1}}+{\sum }_{j=\mathrm{1}}^{i-\mathrm{1}}{a}_{ij}{\mathit{\sigma }}_{m,i}^n\right)\\ -\mathrm{\Delta }t{\sum }_{i=\mathrm{1}}^{\mathrm{4}}\left(b_i{\mathit{\eta }}_m^{n+\mathrm{1}}+{\sum }_{j=\mathrm{1}}^{i-\mathrm{1}}{a}_{ij}{\mathit{\eta }}_{m,i}^n\right)\end{array}\right).\end{array}\end{array}$$

B5 Gradient verification

In this section, Eq. (B27) for the gradient of the objective function, resulting from the temporally discrete adjoint method, is validated via finite differences. For the gradient verification, we consider the same numerical setup from Sect. 3.1. Via finite differences, the Gateaux derivative of the objective function in the direction δφ is approximated as

$$\begin{array}{}\text{(B28)}& \left(\mathrm{\nabla }{\stackrel{\mathrm{̃}}{\mathcal{J}}}^N,\mathit{\delta }\mathit{\phi }\right)\approx {\displaystyle \frac{{\stackrel{\mathrm{̃}}{\mathcal{J}}}^N(\mathit{\phi }+\mathit{\alpha }\mathit{\delta }\mathit{\phi })-{\stackrel{\mathrm{̃}}{\mathcal{J}}}^N\left(\mathit{\phi }\right)}{\mathit{\alpha }}},\end{array}$$

where we set $\mathit{\alpha }={\mathrm{10}}^{-\mathrm{6}}$. To limit computational costs, we only examine the gradients for turbines R1C1 and R8C4 (respectively front and last row in Fig. 6) for a control horizon T=350 s (the longest horizon considered in this work) and for a limited amount of time instants. The baseline controls, resulting adjoint gradients and finite-difference verification are shown in Fig. B1.

Figure B1 Baseline controls for which the gradient is computed and adjoint gradients versus finite-difference verification.

Appendix C: Yaw and induction characteristics

Figures C1–C4 display the time evolution of the optimized filtered thrust coefficients ${\widehat{C}}_{\mathrm{T}}^{\prime }$ and yaw angles θ for the eight turbines in column C1 for respectively grid level 1 and grid level 2. As in Sect. 4.4 for grid level 0, only the turbines in column C1 are shown for cases 1, 4, 5, 7, 8, 9 and 10. Regarding the influence of the receding-horizon parameters, the conclusions are similar to those for grid level 0. Note that for level 1 and level 2, the dips in ${\widehat{C}}_{\mathrm{T}}^{\prime }$ are more prominent than for grid level 0 in Sect. 4.4. However, in each of the cases, yaw control remains the dominant control mechanism.

Figure C1 Time evolution of filtered thrust coefficients ${\widehat{C}}_{\mathrm{T}}^{\prime }$ for turbine column C1 for different optimal control cases for grid level 1.

Figure C2 Time evolution of yaw angles θ for turbine column C1 for different optimal control cases for grid level 1.

Figure C3 Time evolution of filtered thrust coefficients ${\widehat{C}}_{\mathrm{T}}^{\prime }$ for turbine column C1 for different optimal control cases for grid level 2.

Figure C4 Time evolution of yaw angles θ for turbine column C1 for different optimal control cases for grid level 2.

Appendix D: Steady yaw angles from Sood and Meyers (2022)

The steady yaw angles from Sood and Meyers (2022) for the TotalControl Reference Wind Power Plant are shown in Fig.. D1. In their framework, only front-row turbines are significantly yawed to ±30^∘, whereas all other (downstream) turbines remain mostly unyawed. For a comparison with the yaw angles obtained by the proposed LES-based controller, see Sect. 4.5.

Figure D1 Steady yaw angles [^∘] from Sood and Meyers (2022) for the TotalControl Reference Wind Power Plant. Axes in rotor diameter units, with D=178.3 m. Figure adapted from Andersen et al. (2018).