Effect of wind turbine response time on optimal dynamic induction control of wind farms

In this work, we extend recent research efforts on induction-based optimal control in large-eddy simulations of wind farms in the turbulent atmospheric boundary layer. More precisely, we investigate the effect of wind turbine response time to requested power setpoints on achievable power gains. We do this by including a time-filtering of the thrust coefficient setpoints in the optimal control framework. We consider simulation cases restricted to underinduction compared to the Betz limit, as well as cases that also allow overinduction. Optimization results show that, except for the most restrictive underinductive slow-response case, all cases still yield increases in energy extraction in the order of 10% and more.


Introduction
In large-scale wind farms, turbulent wake interactions between turbines can dramatically reduce power extraction efficiencies in downstream rows. More specifically, in currently operational wind farms, in which turbines individually aim to maximize their power extraction, deficits of up to 40% have been reported. Previous efforts to alleviate wake effects and mitigate downstream power deficits by static downrating of upstream rows have often been unsuccessful, i.e. losses in upstream rows are not fully compensated in downstream rows (see, e.g., Refs [1,2]). On the other hand, recent studies on dynamic induction control of wind turbines, taking into account their unsteady interaction with the atmospheric boundary layer (ABL), have shown gains of 16% and 7% in aggregate power extraction for fully-developed and spatially-developing aligned wind farms respectively. [3,4] The latter studies investigate axial induction-based power-maximizing optimal control of wind-farm boundary layers, considering the individual wind turbines as flow actuators. By dynamically controlling the turbines in a coordinated manner, their energy extraction can be regulated with the aim of optimally influencing the turbulent ABL flow for increased global farm power extraction, e.g. through enhanced wake mixing and kinetic energy transport.
A key feature of the induction factors obtained by dynamic optimal control is the presence of strong temporal variability in response to the unsteady turbulent flow dynamics of the ABL. Although technically feasible, such variable induction factors tend to promote power and thrust variability on short timescales as shown in [3,4], and hence contribute significantly to increased fatigue loading of the wind turbines. In this work, we investigate the effect of wind turbine response time on optimal dynamic induction control of wind farms. We do this by time-filtering the setpoints dictated by the optimal control framework, which allows us to assess whether significant gains in energy extraction can still be obtained for smoothed time variations in wind turbine dynamics. Moreover, in addition to reducing temporal fluctuations in turbine thrust forces, we also study the effect of reducing their amplitudes by limiting the maximal thrust coefficient. This work hence contributes towards actual implementation of cooperative windfarm controllers by quantifying the effects of adding increasingly more practical constraints to the optimal control framework.
The paper is structured as follows: Section 2 elaborates on the applied methodology for optimal control of wind-farm boundary layers. Section 3 subsequently describes the cases considered and their numerical setup. Section 4 discusses optimization results and their dependence on wind turbine response time and maximal thrust coefficient. Section 5 finally summarizes the main contributions of this work and provides suggestions for future research.

Methodology
Simulations are performed using the OSP-Wind optimization platform, which has been developed at KU Leuven over recent years [3,5,6]. We maximize wind farm-aggregated energy extraction by minimizing the cost functional J in the following PDE-constrained optimization problem: (1) The set of constraints for the optimization problem in (1) consists of the filtered Navier-Stokes momentum and continuity equation (i.e. the state equations), a time-filter equation on the control inputs C T,i , and associated bound constraints on the latter.
The optimal control problem is solved in a reduced formulation, i.e. instead of treating the Navier-Stokes state equations as a constraint, they are satisfied explicitly throughout the optimization process by means of large-eddy simulations (LES). In these equations, the system state variablesũ andp denote the filtered velocity and pressure field in the ABL. Furthermore, ∇p ∞ is the driving background pressure gradient, f i are the forces enacted on the flow by the N t wind turbines, and τ sgs is the subgrid-scale stress, for which a Smagorinsky model is employed. The state equations are discretized using a fourth order energy-conservative finite-difference scheme in the vertical direction, combined with a Fourier pseudo-spectral discretization in the horizontal directions (see, e.g. Refs. [6,7] for a detailed description of the solver). The latter choice of discretization entails doubly periodic boundary conditions in the lateral directions, while the top and bottom domain boundaries are treated with a free slip and wall stress boundary condition respectively. Time integration is performed using a fourth order accurate explicit Runge-Kutta scheme. Turbulent inflow conditions are generated by running a separate precursor simulation in which shifted periodic boundary conditions are applied to alleviate spanwise locking of turbulent structures and hence accelerate statistical convergence of time-averaged simulation results [8]. Near the domain outlet, the flow field is nudged towards the generated inflow conditionsũ p through the addition of forcing terms with a smooth masking function λ in a so-called fringe region [9].
The turbine thrust forces are calculated using a classical actuator disk model (ADM) as is the geometric rotor footprint on the LES grid, and e ⊥ denotes the unit vector perpendicular to the rotor plane. The disk-based thrust coefficients C T,i serve as the control variables. The time-filter equation in (1) applies a one-sided exponential time filter with a time constant τ to the control variables, resulting in an ordinary differential equation that is discretized using a backward Euler scheme. By varying τ we can hence control the turbine response time to variations in the requested thrust setpoints. In this way, we can assess how much power can still be gained from optimal control while limiting potentially harmful effects due to increased structural loading variability. Note that, instead of simply time-filtering the optimal controls a posteriori, which would lead to sub-optimal results, the time-filtering is taken into account explicitly during the optimization process.
The optimization problem is solved using the quasi-Newton L-BFGS-B algorithm of Byrd et al. [10], which has been shown to outperform the nonlinear conjugate-gradient algorithms that is traditionally used for PDE-constrained optimization problems [11]. Note that, due to the high control space dimensionality (C T ∈ R m , with m = N t × N ∆t ≈ 10 4 − 10 5 ) computing the cost functional gradient with respect to the control variables using classical finite differences is computationally prohibitive. Instead, the gradient is identified using the continuous adjoint method. In this method, an additional set of partial differential equations is solved, in such a way that the calculation of the cost functional gradient becomes tractable. More specifically, by satisfying the following adjoint system, the gradient of the cost functional with respect to the control parameters ∇ C T J can be evaluated as The adjoint equations are formulated by postulating vanishing derivatives of the Lagrangian of (1) with respect to the state variables. The adjoint variables (ξ, π, σ) are hence the Lagrange multipliers associated with the momentum, continuity, and time-filter equations in (1) respectively. The resulting derivation is straightforward, although lengthy. Therefore, the reader is referred to [12,13] for the derivation of the adjoint formulations of the standard transport terms in the state equations on the one hand, and to [3,4] for the case-specific terms regarding the ADM forces f * i , subgrid-scale stress τ * sgs , and fringe forcing λξ on the other. The continuous adjoint equations are discretized with the same numerical schemes as used for the forward state equations.
The optimal wind-farm control is realized using a receding-horizon model predictive control framework as illustrated in Figure 1: at t = 0 s, wind-farm operation is optimized over a time window T , yielding a set of optimal controls C * T . These controls are subsequently employed in a flow advancement simulation over a time horizon T A . The latter is chosen as a compromise between mitigation of finite-horizon end effects on the one hand (favoring T A T ) and computational cost (favoring T A = T ) on the other. The resulting flow field at t = T A is then Receding horizon model-predictive control strategy for optimal wind-farm control employed in this paper.
used as the initial condition for a new optimization window and the procedure begins anew. The global time horizon T tot over which the optimal control is considered is then constructed by stringing the mentioned windows together. The specific parameters used in this work are described in the following section.

Case description and numerical setup
We consider an aligned wind farm consisting of 12 × 6 wind turbines, spaced apart by 6 rotor diameters in both axial and transversal directions. The simulation domain consists of a rectangular box of 10 × 3.6 × 1 km 3 , in which the final 10% of the domain is reserved for the aforementioned fringe region. The governing equations are solved on a simulation grid of 384×256×144 grid points. Time advancement is performed using a constant timestep ∆t = 0.75 s, for which a Courant-Friedrichs-Lewy number of around 0.4 is attained. The flow is forced through the domain using a driving pressure gradient of ∂ x p ∞ = 2.5 × 10 −4 m s −2 , resulting in an upstream velocity of approximately 8.5 m s −1 at hub height. A typical snapshot of the LES velocity field is shown in Figure 2, illustrating the complex and unsteady nature of the interaction between the turbulent ABL and wind turbine wakes throughout the wind farm. The optimization is performed using the aforementioned L-BFGS-B method, where we retain 5 correction pairs to approximate the inverse Hessian required by the quasi-Newton algorithm. We limit the optimization procedure to 60 iterations, leading to approximately 120 PDE evaluations (LES or adjoint), with a total simulation walltime of about 24 hours per window on 320 Ivy Bridge processor cores. The total amount of PDE evaluations depends on whether the Newton step, initially proposed by the L-BFGS-B algorithm, satisfies the strong Wolfe conditions (see, e.g., [14]). If not, a line search algorithm will require additional PDE evaluations until these conditions are fulfilled. In practice however, we find that a line search is only rarely necessary. The controls are optimized over a prediction horizon of T = 240 s, after which they are applied in a flow advancement simulation of T A = 120 s. In total, we optimize wind-farm operation over 15 time windows, resulting in a total time horizon T tot = 1800 s, or 30 minutes. Simulation and optimization parameters are summarized in Table 1.
We consider a reference case (R) in which turbines greedily maximize individual power by setting their thrust coefficients to the Betz-optimal value C T = 2. Using momentum theory, the relationship C P = 64C T /(C T + 4) 3 can be derived, which has 2 solutions for every desired power coefficient C P < 16/27 (i.e. the Betz limit) : an underinductive (C T < 2), and an overinductive (C T > 2) thrust coefficient. We define a set of optimally controlled cases, differentiated by whether or not overinduction is allowed, i.e. with C max T = 3 and 2 respectively, and by the turbine response time τ of 0, 5, or 30 s. This leads to a total of 6 optimal control cases, as summarized in Table 1. Optimized cases are denoted by C<X>t<Y>, where <X> and <Y> represent C max T and τ respectively, e.g. C3t30 for the case with C max T = 3 and τ = 30 s. Note that, although the optimal control problem in (1) is a general nonlinear program for which a multitude of local minima can exist, the application of global optimization methods is well outside the range of present-day computational resources. Tests with varying starting points however indicate that, even though the algorithm converges to a different set of (locally) optimal controls, the overall cost functional decrease does not differ significantly. Therefore, throughout the current work, we limit ourselves to a single starting point with a steady C T = 2 for all turbines in all optimal control cases.
Simulations are initialized as follows: firstly, a fully-developed high-Reynolds number turbulent boundary layer is generated by simulating a perturbed logarithmic flow over a physical timespan of approximately 36 hours. Subsequently, wind turbines are inserted in the simulation, and the flow field is advanced in time until the effects of any wind-farm start-up transients have subsided in the main domain. This wind-farm boundary layer in statistically stationary state is used as initial condition for all simulation cases, and the corresponding time is set as t = 0 s in the remainder of the manuscript.

Optimization results and discussion
In this section, we discuss the results of the optimal control cases in comparison to the greedy reference case, with special attention towards the effect of wind turbine response time τ and maximal thrust coefficient C max T . Firstly, wind farm power extraction and turbine thrust forces are investigated. Afterwards, a concise discussion of flow features is presented. Figure 3 illustrates energy extraction over the entire operation time of 30 minutes, normalized by the values in the reference case R. Figure 3a shows that power gains range between 2 and

Simulation cases
20%. The overall higher yield of all cases except C2t30, as compared to the 7% obtained in Goit, Munters and Meyers [4] is mainly attributed to a better convergence through the use of a better optimization method with more iterations. As expected, reducing C max T and/or increasing response time τ reduces the potential for power gains through dynamic induction control. Interestingly, the overinductive and slow response case C3t30 shows power gains similar to the underinductive but fast response cases C2t0 and C2t5. This shows that, for the considered wind farm, power can be increased significantly in the range of 10% without introducing either short-timescale variability or large amplitude thrust forces. However, C2t30 shows that this is not possible when trying to avoid both. Figure 3b shows the energy extraction per row, normalized by the first row of the reference case. It can be seen that the largest gains are generally achieved in rows 2 and 3, at the cost of slightly downrating the first row. Contrastingly, C3t30 shows a minor increase in first row power. Aside from the first row of turbines, the last row extracts the most energy in all controlled cases, since its turbines should not take into account downstream neighbors and hence maximize their own power. Figure 4 illustrates the dynamic behavior for all cases over 30 minutes of operation time. Figures 4a and 4c depict a front-row turbine thrust coefficient setpoint C T and its time-filtered counterpart C T for the overinductive and underinductive cases respectively. It can be seen  from these figures that C T varies wildly between its upper and lower bound in response to the turbulent dynamics of the boundary layer, and that time-filtering greatly reduces temporal variability in C T . Figures 4b and 4d show the total power extraction compared to the reference power. We see for all cases that optimal control leads to an increase in temporal variations of extracted power, and that time-filtering significantly reduces potentially harmful variability on shorter timescales. The periodic variations on longer timescales that can be observed for all cases except C2t30 are caused by finite-horizon optimization effects. We note that C2t30 hardly differs from the reference case R, with C T ≈ 2 for virtually all turbines throughout the entire time horizon. This indicates that, for underinductive control of the considered windfarm setup, variations on shorter timescales are essential for significant power increase through optimal control. Further analysis of the optimized controls and the physical mechanisms behind increasing power extraction is beyond the scope of this work and is subject of current research. Figure 5 shows the dynamic behavior of the thrust force f t for a turbine in row 6 of the wind farm, where turbulence levels at hub height approach a fully-developed state. Figure 5a illustrates the time series of f t , normalized by the time-averaged values of the same turbine in the reference case. It can be seen that case C3t0 exhibits high frequency force variations and peak forces regularly exceeding values twice as large as the time average reference f t,R . Cases C3t5 and C3t30 show that increasing response time τ not only filters out high frequency force fluctuations, but also leads to a reduction of the largest force amplitudes for this wind turbine. A further decrease in amplitudes is achieved by reducing C max T to 2. Figure 5b shows the power spectra of the time series illustrated in 5a. Low frequency behavior (≤ 10 −2 Hz) is similar for all cases and tends to a -1 slope, indicating that the slowest force variations are governed by the background atmospheric dynamics. Around 10 −2 Hz, a small -5/3 range is observed. At higher frequencies, a slope between -3 and -4 characterizes the controlled cases, while the static reference case decays much more quickly. In this range, the turbine reponse time to variations in control setpoint directly translates to the frequency content of the thrust force enacted by the turbine. It is important to note that, not only should fatigue be taken into account when considering additional cyclic turbine loading due to dynamic wind farm control, also natural frequencies of relevant deformation modes should be avoided at all cost. For example, the natural frequency associated with the tower fore-aft deformation mode of the NREL 5 MW turbine is approximately 0.3 Hz [15]. It can be seen from Figure 5 that cases C3t0 and C2t0   show significant spectral energy in this range. These results indicate that a detailed analysis of turbine structural loading, through the addition of a multi-physics aero-elastic turbine model in the governing equations is an essential step towards actual rollout of fast-response dynamic wind farm controllers. This is subject of ongoing work.
To conclude the discussion on the results of the optimal control simulations, we provide a qualitative view of the boundary layer flow through the wind farm in Figure 6. Figures 6a and  6b illustrate the time-averaged streamwise velocity through the rotor centerline as a function of streamwise distance, averaged over all 6 columns. It is shown that all cases exhibit higher velocities at the turbine disks starting from row 2, which directly leads to the observed increase in energy extraction discussed above. The difference in disk velocity for row 1 is very small, corresponding with the minor differences in first row power shown previously in Figure 3. In addition, overinductive cases with C max T = 3 generate deeper near-wakes behind the first row and rows 8 to 12, while this is not the case for underinductive cases with C max T = 2, where velocities are higher over the entire streamwise extent of the domain. Overall, the figure shows that, for the overinductive case, trough-to-peak recovery between near and far wake is enhanced significantly, while for underinductive cases wake recovery remains roughly constant.
This observation is also supported by the turbulence kinetic energy (TKE) levels at hub height throughout the wind farm, as illustrated in Figures 6c and 6d. For C max T = 3, higher turbulence energy is observed between every row of turbines, consistent with an increase in turbulent mixing leading to better wake recovery. Interestingly, even the slow response case C3t30 shows a significant increase in TKE compared to the reference case, even though turbine thrust forces and power extraction vary smoothly in time. For C max T = 2, TKE levels are very similar to the static reference case.

Conclusions
In this work we continue the application of optimal dynamic induction control for power maximization in LES of aligned wind farms started in earlier work [3,4]. More specifically, we discuss the sensitivity of achievable gains in energy extraction with respect to wind turbine response times, with the aim of reducing detrimental effects associated with fast-varying optimized thrust coefficients obtained by the latter studies. We consider simulation cases restricted to underinduction, as well as cases that also allow overinduction. In summary, the results obtained in this work show that cooperative wind farm control still yields significant power gains after eliminating potentially structurally harmful characteristics of the dynamic control schemes, i.e. strong temporal variations and large thrust forces. For instance, cases C3t30 and C2t5 show increases in energy extraction compared to greedy control strategies of 8 and 10% respectively, which would still involve very significant economic profits. Spectral analysis of thrust forces on a wind turbine indicate the need for multi-physics optimization to accurately model additional fatigue loading and avoid the resonance excitation of deformation modes. Near-future work will focus on identifying the main mechanisms behind the observed increases in energy extraction, and translating the computationally expensive optimal controller to practical suboptimal control schemes.