Multifidelity multiobjective optimization for wake-steering strategies

Wake steering is an emerging wind power plant control strategy where upstream turbines are intentionally yawed out of perpendicular alignment with the incoming wind, thereby “steering” wakes away from downstream turbines. However, trade-offs between the gains in power production and fatigue loads induced by this control strategy are the subject of continuing investigation. In this study, we present a multifidelity multiobjective optimization approach for exploring the Pareto front of trade-offs between power and loading during wake steering. A large eddy simulation is used as the high-fidelity model, where an actuator line representation is used to model wind turbine blades and a rainflow-counting algorithm is used to compute damage equivalent loads. A coarser simulation with a simpler loads model is employed as a supplementary low-fidelity model. Multifidelity Bayesian optimization is performed to iteratively learn both a surrogate of the low-fidelity model and an additive discrepancy function, which maps the low-fidelity model to the high-fidelity model. Each optimization uses the expected hypervolume improvement acquisition function, weighted by the total cost of a proposed model evaluation in the multifidelity case. The multifidelity approach is able to capture the logit function shape of the Pareto frontier at a computational cost only 30 % that of the single-fidelity approach. Additionally, we provide physical insights into the vortical structures in the wake that contribute to the Pareto front shape.

1 Introduction

As wind energy systems have matured, plant-level control has emerged as a new paradigm, where groups of turbines are controlled in coordination to maximize collective power production. This is in contrast to more traditional control strategies, where individual turbines are controlled to maximize their own power production. A potentially promising form of such plant-level control is “wake steering”, where upstream wind turbine yaw positions are intentionally misaligned from the incoming wind, “steering” the wake away from downstream turbines. A counter-rotating pair of vortices is generated by the lateral thrust of the wind turbine rotor, which is determined by the yaw offset direction (Fleming et al., 2018; Martínez-Tossas and Branlard, 2020). This allows the performance of wind power plants to be improved by diverting wakes away from downstream turbines

It is speculated that wake steering may produce more power while inducing less total fatigue on all turbines when compared to the baseline strategy of aligning each turbine with the incoming wind (Howland et al., 2019; Hulsman et al., 2020). However, very few studies have quantified the trade-offs between power and damage. Hulsman et al. (2020) used an actuator line model to train polynomial chaos surrogates for optimization of a weighted sum of power and damage equivalent loads. Yin et al. (2020) present a multiobjective genetic algorithm for the maximization of power and minimization of total thrust. Damiani et al. (2018) performed a detailed analysis of a single wind turbine, finding that negative yaw offsets tended to increase fatigue loading more than positive yaw offsets, but cautioned that this result could not be generalized given the dependence on incident conditions. Other studies have provided additional evidence supporting the notion that positive yaw misalignment results in less edgewise loading of downstream turbine blades than the corresponding negative yaw misalignment strategy (Zalkind and Pao, 2016; López et al., 2020). Wang et al. (2020) demonstrated the potential of individual pitch control to alleviate loads induced from intentionally offsetting the turbine yaw. Van Dijk et al. (2017) used the FLORIS and CCBlade tools to examine trade-offs between power produced and the edgewise and flapwise fatigue loading induced through wake steering. Lin and Porté-Agel (2020) utilize a large eddy simulation (LES) framework to construct the Pareto set between power and flapwise bending moment loading through a comprehensive parameter sweep. While these studies all provide insights into the trade-offs between power and loading, there is still a need for an efficient optimization algorithm to quantify these trade-offs using computationally intensive simulations.

Despite its promise, plant-level control via wake steering involves complex physics and is challenging to model. Engineering wake models have dubious accuracy when predicting fatigue loading, which higher-fidelity models predict more accurately (Rinker et al., 2021). In this study, we propose a multifidelity multiobjective optimization framework to address this challenge and explore trade-offs between power and loading in wake-steering strategies. In practice, power and loading will likely be optimized in real time using a singular weighted objective. The relative weights may be decided upon by exploring trade-offs between power and loading using multiobjective optimization to estimate the Pareto frontier. When searching for the Pareto set, an efficient algorithm must balance exploration and exploitation. Several models have been developed that may be used to study the effects of control strategies with various levels of mathematical detail and real-world accuracy (i.e., fidelity) (Annoni et al., 2018; Martínez-Tossas et al., 2019; Hulsman et al., 2020).

Multifidelity optimization exploits the correlation between low- and high-fidelity models to reduce the overall computational cost of optimization. For instance, Andersson and Imsland (2020) present a real-time modifier adaptation approach for wake-steering design, where a Gaussian process (GP) is used to iteratively learn the difference between observed operational data and the predictions of an engineering wake model. Ariyarit and Kanazaki (2017) present a two-objective bifidelity approach that iteratively builds a GP discrepancy function. Huang et al. (2006) and Rajnarayan et al. (2008) employ an augmented expected improvement formulation, including three factors to account for the correlation between the low- and high-fidelity models, the observed error, and the cost ratio between the low- and high-fidelity models. It is not always clear when a proposed low-fidelity model is appropriate for use in multifidelity optimization, though it is common to assess candidate low-fidelity models by measuring their correlation with the high-fidelity model (Giselle Fernández-Godino et al., 2019).

The novelty of the present study is the application of this multifidelity technique to wind energy systems, resulting in new insights into wake-steering flow physics. The present approach uses the low-fidelity model to first explore the full parameter space and then iteratively builds the low- and high-fidelity model surrogates to gain the most improvement in the Pareto front per model evaluation costs. While this framework is similar to that presented by Ariyarit and Kanazaki (2017) and Andersson and Imsland (2020), the exact framework outlined here is new, and this is the first demonstration of any such approach in the context of wind energy systems.

2 Optimization framework

A Bayesian framework for multiobjective multifidelity optimization is presented. Throughout this section, we assume that minimization of functions is the objective of the optimization procedure (as opposed to maximization).

This study employs GP models to approximate power and loading dynamics. A GP is a collection of random variables, any finite number of which have a joint Gaussian distribution (Rasmussen and Williams, 2006). The simulated power and loads, f_i, are approximated using individual Gaussian process surrogate models, g, which are defined as

$$\begin{array}{}\text{(1)}& g_i\left(\mathit{\gamma }\right)\sim \mathcal{G}\mathcal{P}\left[{\mathit{\mu }}_i\left(\mathit{\gamma }\right),k_i(\mathit{\gamma },{\mathit{\gamma }}^{\prime })\right],\end{array}$$

where γ is a vector of proposed yaw angles (with dimension equal to the number of turbines), γ^′ is an arbitrary vector of yaw angles, μ_i(γ) is the GP mean function, $k_i(\mathit{\gamma },{\mathit{\gamma }}^{\prime })$ is the GP kernel covariance function, and the index i refers to power (i=1) or loading (i=2) objectives.

We perform Bayesian inference on functions by conditioning the GP on a set of observed input–output pairs ${\mathcal{D}}_i=\left\{\mathbf{\Gamma },{\mathbf{Y}}_i\right\}$, where $\mathbf{\Gamma }=\left[{\mathit{\gamma }}^{\left(\mathrm{1}\right)},{\mathit{\gamma }}^{\left(\mathrm{2}\right)},\mathrm{\dots }\right]$ is a matrix of simulated yaw offsets, and ${\mathbf{Y}}_i=\left[f_i\left({\mathit{\gamma }}^{\left(\mathrm{1}\right)}\right),f_i\left({\mathit{\gamma }}^{\left(\mathrm{2}\right)}\right),\mathrm{\dots }\right]$ is a vector of observations of simulated power or loading outputs. We use the scikit-learn Gaussian process implementation (Pedregosa et al., 2011), which is a well-validated open-source project. The power and loading outputs are normalized to have zero mean and unity variance during a preprocessing step. After conditioning on 𝒟_i, we obtain Gaussian distributions at test locations γ^* with the posterior estimates of the mean,

$$\begin{array}{}\text{(2)}& {\mathit{\mu }}_i\left({\mathit{\gamma }}^{*}|{\mathcal{D}}_i\right)={\mathit{k}}_{i,*}^T{\mathit{K}}_i^{-\mathrm{1}}{\mathbf{Y}}_i,\end{array}$$

and variance,

$$\begin{array}{}\text{(3)}& {\mathit{\sigma }}_i^{\mathrm{2}}\left({\mathit{\gamma }}^{*}|{\mathcal{D}}_i\right)=k_i\left({\mathit{\gamma }}^{*},{\mathit{\gamma }}^{*}\right)-{\mathit{k}}_{i,*}^T{\mathit{K}}_i^{-\mathrm{1}}{\mathit{k}}_{i,*},\end{array}$$

where ${\mathit{k}}_{i,*}=k_i\left(\mathbf{\Gamma },{\mathit{\gamma }}^{*}\right)$ is a vector and ${\mathit{K}}_i=k_i\left(\mathbf{\Gamma },\mathbf{\Gamma }\right)$ is a matrix.

The kernel covariance function encodes prior knowledge about structural properties of the underlying signal, such as smoothness, periodicity, and stationarity. In this study, we employ an anisotropic radial basis function kernel (Rasmussen and Williams, 2006) for k_i, given by

$$\begin{array}{}\text{(4)}& k_i(\mathit{\gamma },{\mathit{\gamma }}^{\prime })=\exp \left[{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left(\sum _{j=\mathrm{1}}^{\dim \left(\mathit{\gamma }\right)}{\displaystyle \frac{({\mathit{\gamma }}_j-{\mathit{\gamma }}_j^{\prime }{)}^{\mathrm{2}}}{l_{ij}^{\mathrm{2}}}}\right)}^{\mathrm{1}/\mathrm{2}}\right],\end{array}$$

where the correlation scale, l_ij, is estimated by maximizing the log-marginal likelihood function (Pedregosa et al., 2011).

2.1 Single-fidelity approach

In Bayesian optimization, an acquisition function is defined to maximize a metric representing both exploration and exploitation (Shahriari et al., 2015). A popular acquisition function for single objective optimization is the expected improvement (EI), which is the expected value of an improvement function (Zhan and Xing, 2020) with respect to the predicted uncertainty of a GP. This acquisition function is employed by the efficient global optimization (EGO) algorithm (Jones et al., 1998). The improvement function, I, quantifies the improvement in the objective function for a new evaluation, as compared to the best sampled objective, and is zero if the new objective does not outperform all of the previously sampled points. This results in

$$\begin{array}{}\text{(5)}& I\left(\mathit{\gamma }\right)=max\left[f^{*}-f\left(\mathit{\gamma }\right),\mathrm{0}\right],\end{array}$$

where f^* is the minimum sampled value and f(γ) is the function value, which is generally unknown and must be predicted by a GP.

There is a range of potential outcomes from sampling a new point, and the GP framework conveniently estimates this uncertainty. These uncertainties are used to compute the expected value of the improvement function. It is important that the improvement function contains the maximum function; otherwise, there would be no exploration of regions of larger uncertainty. Other acquisition functions available include the knowledge gradient (Ghoreishi and Allaire, 2018), expected quantile improvement (He et al., 2017; Picheny et al., 2013), improved expected improvement (Qin et al., 2017), entropy search (Hennig and Schuler, 2012), and minimization of the predictor (Andersson and Imsland, 2020).

The expected improvement may be extended to a multiobjective context. This is done by introducing a hypervolume function, H, which measures the volume of a given Pareto front, A, using a reference point, r. The expected hypervolume improvement (EHVI), introduced by Emmerich et al. (2006), is the multiobjective counterpart to the expected improvement acquisition function used in the EGO algorithm and is given by

$$\begin{array}{}\text{(6)}& \text{EHVI}\left(\mathit{\gamma },\mathit{g}\right)={\mathbb{E}}_{\mathit{g}\left(\mathit{\gamma }\right)}\left\{\text{HVI}\left[\mathit{f}\left(\mathit{\gamma }\right)\right]\right\}.\end{array}$$

Here, 𝔼_g(γ)[⋅] represents the expectation with respect to a normal distribution, g(γ), and is expressed as

$$\begin{array}{}\text{(7)}& \begin{array}{rl}& {\mathbb{E}}_{\mathit{g}\left(\mathit{\gamma }\right)}\left\{\text{HVI}\left[\mathit{f}\left(\mathit{\gamma }\right)\right]\right\}=\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\text{HVI}\left(\right[P,L\left]\right)\\ & {\displaystyle \frac{\mathrm{1}}{{\mathit{\sigma }}_{\mathrm{1}}\left(\mathit{\gamma }\right)\sqrt{\mathrm{2}\mathit{\pi }}}}{e}^{-\frac{\mathrm{1}}{\mathrm{2}}{\left[\frac{P-{\mathit{\mu }}_{\mathrm{1}}\left(\mathit{\gamma }\right)}{{\mathit{\sigma }}_{\mathrm{1}}\left(\mathit{\gamma }\right)}\right]}^{\mathrm{2}}}{\displaystyle \frac{\mathrm{1}}{{\mathit{\sigma }}_{\mathrm{2}}\left(\mathit{\gamma }\right)\sqrt{\mathrm{2}\mathit{\pi }}}}{e}^{-\frac{\mathrm{1}}{\mathrm{2}}{\left[\frac{L-{\mathit{\mu }}_{\mathrm{2}}\left(\mathit{\gamma }\right)}{{\mathit{\sigma }}_{\mathrm{2}}\left(\mathit{\gamma }\right)}\right]}^{\mathrm{2}}}\mathrm{d}P\mathrm{d}L,\end{array}\end{array}$$

where μ₁ and σ₁ are the mean and standard deviation, respectively, of the powers modeled by g₁, and μ₂ and σ₂ are the mean and standard deviation, respectively, of loads modeled by g₂. The hypervolume improvement indicator (HVI) function is the multiobjective counterpart to the improvement indicator function. It is given as

$$\begin{array}{}\text{(8)}& \text{HVI}\left(\right[P,L\left]\right)=H\left(A\cup \left\{[P,L]\right\}\right)-H\left(A\right),\end{array}$$

where A is an estimated Pareto frontier; $A\cup \mathit{\left\{}\right[P,L\left]\mathit{\right\}}$ is a new Pareto set that potentially includes [P,L]; and H is a hypervolume function, which measures the volume spanned by the Pareto set of objective functions relative to a reference point, r, which must not be dominated by $A\cup \left\{[P,L]\right\}$. The Pareto frontier is defined as the set of all function values that are not strictly dominated by other function values. The formal definition is

$$\begin{array}{}\text{(9)}& \begin{array}{rl}A=\mathit{\{}{\mathit{y}}^{\prime }\in & \left\{\mathit{y}\in {\mathbb{R}}^{\dim \left(\mathit{f}\right)}:\mathit{y}=\mathit{f}\left(\mathit{\gamma }\right),\mathit{\gamma }\in {\mathrm{\Omega }}_{\mathit{\gamma }}\right\}\\ & :\mathit{\{}\mathit{y}{}^{\prime \prime }\prec {\mathit{y}}^{\prime },\mathit{y}{}^{\prime \prime }\ne {\mathit{y}}^{\prime }\mathit{\}}=\mathrm{\varnothing }\mathit{\}},\end{array}\end{array}$$

where Ω_γ is the set of allowable yaw offsets, ≺ denotes Pareto dominance (that is, if $\mathit{y}{}^{\prime \prime }\prec {\mathit{y}}^{\prime }$, then $y{{}^{\prime \prime }}_i\le y_i^{\prime }$ for all values of i and $y{{}^{\prime \prime }}_i<y_i^{\prime }$ for at least one value of i; Voorneveld, 2003), and ∅ is the empty set.

The hypervolume, H, measures the extent of the Pareto set as the volume of the Pareto-dominated space bounded by a reference point, r, namely,

$$\begin{array}{}\text{(10)}& H\left(A\right)=\text{Vol}\left(\left\{\mathit{y}\in {\mathbb{R}}^{\dim \left(\mathit{f}\right)}\left|{\mathit{y}}^{\prime }\in A\prec \mathit{y}\right.\text{and}\mathit{y}\prec \mathit{r}\right\}\right).\end{array}$$

In practice, the Pareto set is computed by filtering a set of discrete inputs so that only non-dominated points remain,

$$\begin{array}{}\text{(11)}& A=\mathit{\{}{\mathit{y}}^{\prime }\in \mathbf{Y}:\left\{\mathit{y}{}^{\prime \prime }\prec {\mathit{y}}^{\prime },\mathit{y}{}^{\prime \prime }\ne {\mathit{y}}^{\prime }\right\}=\mathrm{\varnothing }\mathit{\}},\end{array}$$

where Y is a matrix of observed function values. This filters out observed samples that are Pareto dominated by other observed samples.

Although these ideas may also be extended to more objectives, assuming two objectives simplifies the problem, in this case, the observed Pareto set is defined as $A\approx ({\mathit{y}}^{\mathrm{1}},{\mathit{y}}^{\mathrm{2}},\mathrm{\dots },{\mathit{y}}^n$) such that $y_{\mathrm{1}}^{\mathrm{1}}<y_{\mathrm{1}}^{\mathrm{2}}<\mathrm{\dots }<y_{\mathrm{1}}^n$. The hypervolume is estimated using rectangular quadrature as

$$\begin{array}{}\text{(12)}& H\left(A\right)\approx \sum _{i=\mathrm{1}}^{n-\mathrm{1}}(y_{\mathrm{1}}^{i+\mathrm{1}}-y_{\mathrm{1}}^i)(r_{\mathrm{2}}-y_{\mathrm{2}}^i)+(r_{\mathrm{1}}-y_{\mathrm{1}}^n)(r_{\mathrm{2}}-y_{\mathrm{2}}^n),\end{array}$$

where n is the number of points in the given Pareto set and r₁ and r₂ are the components of the reference point. In this study, the EHVI is approximated through Monte Carlo simulation by

$$\begin{array}{}\text{(13)}& \begin{array}{rl}& \text{EHVI}(\mathit{\gamma },\mathit{g})\\ & \approx {\displaystyle \frac{\mathrm{1}}{N_{\mathrm{s}}}}\sum _{k=\mathrm{1}}^{N_{\mathrm{s}}}\left[H\left(A\cup {\left\{\mathcal{N}\left[\mathit{\mu }\right(\mathit{\gamma }),\mathit{\sigma }(\mathit{\gamma }\left)\right]\right\}}^{\left(k\right)}\right)-H\left(A\right)\right],\end{array}\end{array}$$

where N_s is the number of Monte Carlo samples and {𝒩[μ(γ),σ(γ)]}^(k) is draw k from the GP model of power and loading, g.

Once the EHVI is estimated, it must be maximized. This is not necessarily trivial, as the EHVI computation is complicated and difficult to vectorize, and the cost of the optimization grows exponentially with the number of design variables. The EHVI optimum may be determined using a grid search, random sampling, direct optimization, or surrogate-based optimization. While a grid search is the most comprehensive option, the latter approaches are more computationally efficient for high-dimensional design inputs.

2.2 Multifidelity approach

The multifidelity approach introduces computationally cheaper but lower fidelity representations of the high-fidelity model, which allow for greater control between exploration and exploitation in the Bayesian optimization. Samples of the low-fidelity model are adaptively refined throughout the optimization as a cheap means for exploration of the high-fidelity function space. Throughout this section, we assume a known hierarchy of model fidelities, $(f_k^{\mathrm{1}},f_k^{\mathrm{2}},\mathrm{\dots },f_k^N)$, where $f_k^{\mathrm{1}}$ is the lowest-fidelity model of power/loading, N is the number of different fidelity models, and $f_k^N$ is the highest-fidelity model of power/loading. When a model is evaluated at a point, γ, we assume that all lower-fidelity models will also be evaluated at this point.

The lowest-fidelity model, $f_k^{\mathrm{1}}$, is approximated using a GP, $g_k^{\mathrm{1}}$, resulting in the following output distribution:

$$\begin{array}{}\text{(14)}& g_k^{\mathrm{1}}\left(\mathit{\gamma }\right)\sim \mathcal{N}\left[{\mathit{\mu }}_k^{\mathrm{1}}\left(\mathit{\gamma }\right),{\mathit{\sigma }}_k^{\mathrm{1}}\left(\mathit{\gamma }\right)\right],\end{array}$$

where ${\mathit{\mu }}_k^{\mathrm{1}}$ and ${\mathit{\sigma }}_k^{\mathrm{1}}$ are the mean and standard deviations, respectively, associated with the lowest-fidelity power/loading model. Higher-fidelity models, $f_k^i\left(\mathit{\gamma }\right)$, are approximated using additive discrepancy functions that map the next lowest fidelity function, $f_k^{i-\mathrm{1}}\left(\mathit{\gamma }\right)$, to $f_k^i\left(\mathit{\gamma }\right)$:

$$\begin{array}{}\text{(15)}& f_k^i\left(\mathit{\gamma }\right)\approx f_k^{i-\mathrm{1}}\left(\mathit{\gamma }\right)+{\mathit{\delta }}_k^i\left(\mathit{\gamma }\right)\forall i>\mathrm{1},\end{array}$$

where ${\mathit{\delta }}_k^i\left(\mathit{\gamma }\right)$ is the discrepancy function associated with objective k and fidelity i,

$$\begin{array}{}\text{(16)}& {\mathit{\delta }}_k^i\left(\mathit{\gamma }\right)\sim \mathcal{N}\left[{\mathit{\mu }}_k^i\left(\mathit{\gamma }\right),{\mathit{\sigma }}_k^i\left(\mathit{\gamma }\right)\right]\forall i>\mathrm{1},\end{array}$$

where ${\mathit{\mu }}_k^i$ and ${\mathit{\sigma }}_k^i$ are the mean and standard deviations, respectively, modeled by the discrepancy function GP associated with fidelity i.

New GPs are defined to extend the EHVI to a multifidelity context. No matter which fidelity is to be sampled next, the ultimate goal is to minimize the highest fidelity functions, so each GP is constructed to predict the high-fidelity outputs. However, GPs associated with lower-fidelity models should not take into account uncertainties associated with higher-fidelity models, as these uncertainties will not be collapsed if the lower-fidelity model is sampled. Sampling the highest-fidelity model must take all sources of surrogate uncertainty into account, as a high-fidelity model evaluation will be associated with sampling all lower-fidelity models at the same point. So, new GP models, hⁱ, are constructed to predict the high-fidelity output while encoding different uncertainty information. The GPs associated with each fidelity are defined as

$$\begin{array}{}\text{(17)}& h_k^i\left(\mathit{\gamma }\right)\sim \mathcal{N}\mathit{\left\{}\sum _{j=\mathrm{1}}^N{\mathit{\mu }}_k^j\right(\mathit{\gamma }),\sqrt{\sum _{j=\mathrm{1}}^i{\left[{\mathit{\sigma }}_k^j\left(\mathit{\gamma }\right)\right]}^{\mathrm{2}}}\mathit{\}}.\end{array}$$

Putting the above formulations together, it is natural to define a multifidelity multiobjective acquisition function as the ratio of EHVI per evaluation cost:

$$\begin{array}{}\text{(18)}& J(\mathit{\gamma },i)=-{\displaystyle \frac{\text{EHVI}\left(\mathit{\gamma },{\mathit{h}}^i\right)}{\sum _{j=\mathrm{1}}^i{C}_j}},\end{array}$$

where J is the optimization acquisition function, which is a function of a set of proposed yaw offsets and model fidelity, to be minimized with respect to yaw offsets, γ, and model fidelity, i, in each optimization iteration, and C_i is the computational cost associated with model i.

In this study, we examine the bifidelity case (N=2), where

$$\begin{array}{}\text{(19)}& h_k^{\mathrm{1}}\left(\mathit{\gamma }\right)\sim \mathcal{N}\left[{\mathit{\mu }}_k^{\mathrm{1}}\left(\mathit{\gamma }\right)+{\mathit{\mu }}_k^{\mathrm{2}}\left(\mathit{\gamma }\right),{\mathit{\sigma }}_k^{\mathrm{1}}\left(\mathit{\gamma }\right)\right]\end{array}$$

and

$$\begin{array}{}\text{(20)}& h_k^{\mathrm{2}}\left(\mathit{\gamma }\right)\sim \mathcal{N}\mathit{\left\{}{\mathit{\mu }}_k^{\mathrm{1}}\right(\mathit{\gamma })+{\mathit{\mu }}_k^{\mathrm{2}}(\mathit{\gamma }),\sqrt{{\left[{\mathit{\sigma }}_k^{\mathrm{1}}\left(\mathit{\gamma }\right)\right]}^{\mathrm{2}}+{\left[{\mathit{\sigma }}_k^{\mathrm{2}}\left(\mathit{\gamma }\right)\right]}^{\mathrm{2}}}\mathit{\}}.\end{array}$$

The bifidelity workflow is visualized in Fig. 1. The EHVIs associated with the low-fidelity (h¹) and high-fidelity (h²) GP models are maximized. Then, the EHVI per cost is maximized with respect to fidelity in the comparative step, which corresponds to minimizing J(γ,i). If the EHVI per cost associated with evaluating the low-fidelity model is greater than the EHVI per cost associated with evaluating the high-fidelity model, the low-fidelity model is evaluated. Otherwise, the high-fidelity model is evaluated. In the figure, h¹ is represented by GP^LF, h² is represented by GP^HF, C₁ is represented by C_LF, and C₂ is represented by C_HF.

Figure 1 Workflow visualization for the bifidelity optimization case. h¹ is represented by GP^LF, h² is represented by GP^HF, C₁ is represented by C_LF, and C₂ is represented by C_HF.

3 Numerical approach

This section outlines the numerical approaches used in this study. Section 3.1 presents the flow modeling framework employed. Section 3.2 introduces the specific power and loading objective functions used in this study. Section 3.3 presents the approach used to maximize the multiobjective acquisition function. Section 3.3.1 outlines the sampling approach, and Sect. 3.3.2 presents a correlation analysis used to determine the low-fidelity loading proxy.

3.1 Flow modeling

We use the WindSE framework (National Renewable Energy Laboratory, 2021) to model flow within the wind power plant. We investigate a two-turbine case, with a single wind direction and speed, where the turbines are spaced seven rotor diameters apart and the wind direction is such that the front turbine directly wakes the back turbine. The large turbine spacing was chosen to ensure that solutions associated with optimal power production were inside of the boundaries of allowable yaw offsets. When turbines are spaced tightly, it is common for the optimal power to be associated with the largest allowable yaw offset of the front turbine, which is a less challenging optimization case. The inflow boundary is modeled using a logarithmic profile with a hub height wind speed of 7.5 m s⁻¹. The top, side, and outflow boundaries are specified as no-stress boundaries, and the ground is specified as a no-slip boundary. We consider turbine representations of the IEA 3.4 MW reference turbine (Bortolotti et al., 2019), with hub heights of 120 m and rotor diameters of 130 m. The turbine blades are represented as actuator lines with 15 force nodes. This analysis does not consider the turbine nacelle or tower, and the turbine blades are modeled as being rigid. There is no turbine controller, and the rotor speed and blade pitch angles are modeled as constants. The domain is represented with a $\mathrm{2260}\times \mathrm{2000}\times \mathrm{520}$ m³ mesh, corresponding to a 301 s flow-through time. The mesh is refined near the center of the domain and where the turbines are located. A target Courant–Friedrichs–Lewy condition of 0.98 is specified. All simulations were initiated with the same atmospheric conditions.

The simulations solve the filtered conservation of mass and Navier–Stokes equations given by

$$\begin{array}{}\text{(21)}& {\displaystyle }\mathrm{\nabla }\cdot \left(\mathit{\rho }\stackrel{\mathrm{̃}}{\mathit{u}}\right)=\mathrm{0},\text{(22)}& {\displaystyle \frac{D\stackrel{\mathrm{̃}}{\mathit{u}}}{Dt}}=-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}\mathrm{\nabla }\stackrel{\mathrm{̃}}{p}+\left({\displaystyle \frac{\mathit{\mu }}{\mathit{\rho }}}+{\mathit{\nu }}_{\mathrm{t}}\right){\mathrm{\nabla }}^{\mathrm{2}}\stackrel{\mathrm{̃}}{\mathit{u}}+\mathcal{F},\end{array}$$

where $D/Dt=\partial /\partial t+\stackrel{\mathrm{̃}}{\mathit{u}}\cdot \mathrm{\nabla }$ is the material derivative, $\stackrel{\mathrm{̃}}{\mathit{u}}$ is the velocity, t is time, x is the spatial location, ∇ is the spatial gradient, ρ is the density, $\stackrel{\mathrm{̃}}{p}$ is the pressure, μ is the dynamic viscosity, ν_t is the turbulent viscosity, and ℱ is the turbine forcing. The density is specified as ρ=1 kg m⁻³, and the dynamic viscosity is specified as $\mathit{\mu }=\mathrm{1.8}\times {\mathrm{10}}^{-\mathrm{5}}$ kg m⁻¹ s⁻¹. The turbulent viscosity, ν_t, is modeled using the Smagorinsky–Lilly LES model as

$$\begin{array}{}\text{(23)}& {\mathit{\nu }}_{\mathrm{t}}=C_{\mathrm{s}}^{\mathrm{2}}{\mathrm{\Delta }}^{\mathrm{2}}\left|\mathbf{S}\right|,\end{array}$$

where C_s=0.17, Δ is the grid cell size, and S is the strain rate tensor given by

$$\begin{array}{}\text{(24)}& \mathbf{S}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[\mathrm{\nabla }\stackrel{\mathrm{̃}}{\mathit{u}}+{\left(\mathrm{\nabla }\stackrel{\mathrm{̃}}{\mathit{u}}\right)}^T\right],\end{array}$$

and $\left|\mathbf{S}\right|={\left(\mathrm{2}\mathbf{S}:\mathbf{S}\right)}^{\mathrm{1}/\mathrm{2}}$.

The turbine forcing is computed as

$$\begin{array}{}\text{(25)}& \begin{array}{rl}\mathcal{F}\left(\mathit{x}\right)=& \sum _{k=\mathrm{1}}^{\dim \left(\mathit{\gamma }\right)}\sum _{b=\mathrm{1}}^{\mathrm{3}}\sum _{j=\mathrm{1}}^{N_{\mathrm{nodes}}}{\mathit{f}}_{kbj}\left(\mathit{x}\right){\displaystyle \frac{\mathrm{1}}{{\mathit{\pi }}^{\mathrm{3}/\mathrm{2}}{\mathit{ϵ}}^{\mathrm{3}}}}\\ & \exp \left[-{\displaystyle \frac{\left|{\mathit{d}}_{kbj}\right(\mathit{x}){|}^{\mathrm{2}}}{{\mathit{ϵ}}^{\mathrm{2}}}}\right],\end{array}\end{array}$$

where N_nodes is the number of actuator nodes per blade; ϵ is the characteristic width of the actuator forces (specified as 2 m in this study); d_kbj(x) is the distance from node j associated with blade b and turbine k; and f_kbj(x) is the actuator force on node j associated with blade b and turbine k, which is computed based on the wind speed, angle of attack, and the reference airfoil lift and drag coefficients, as well as a tip loss correction, as described by Allen et al. (2022).

Low- and high-fidelity models were developed for this study using the WindSE framework. A Cartesian discretization of the computation domain is specified, where the grid is refined twice in the wake region as well as near the turbine rotors. Each high-fidelity simulation is run to 1200 s using Taylor–Hood elements (Ern and Guermond, 2004). The power and loading results only use information from the final 600 s of simulation time. These time parameters were justified by comparing power and loading computed over time intervals of 600–900 and 900–1200 s, resulting in relative differences of only 2.6 % for power and 4.2 % for loading when $\mathit{\gamma }=(\mathrm{15}{}^{\circ },\mathrm{0}{}^{\circ })$. These time parameters were chosen to obtain a reasonably efficient optimization problem for the present demonstration, accepting the possibility that the optimization results could be slightly different if different time parameters were chosen. For example, using a time interval from 1200 to 1800 s, instead of the present 600–1200 s, resulted in changes to the hypervolumes of the Pareto fronts discovered by the single-fidelity and multifidelity optimizations of less than 0.5 %. The low-fidelity model was selected using the same grid as the high-fidelity model but runs to 400 simulation seconds and uses piecewise linear elements. The power is averaged after 300 s and the loading is estimated using the front and back turbine moments past 300 s. This cut-in time and the total low-fidelity model simulation time were selected to avoid initial transient effects while minimizing the computational cost of the low-fidelity simulation. Using eight processors, the high-fidelity model was measured to take 5.4 h to run, and the low-fidelity model was measured to take 0.24 h to run. This corresponds to a cost ratio of approximately 0.05. The flow fields produced by the low- and high-fidelity models are compared in Fig. 2.

Figure 2 Time-averaged velocity magnitude fields at the turbine hub height associated with the low- and high-fidelity models. Panel (a) shows the flow field associated with the low-fidelity model, and panel (b) shows the flow field associated with the high-fidelity model. In both cases, the front turbine is offset by 30^∘, and the back turbine has no yaw offset. Brighter colors correspond to faster velocity magnitudes. Turbine positions are shown with dark lines.

3.2 Objective functions

The objective of the optimization is to minimize negative power, −P, and loading, L, with respect to each turbine yaw offset such that the yaw offset angles, γ, are bounded between −30 and 30^∘.

Power and loading are quantified using the actuator line model results, discarding an initial transient period. Power is computed as the average total power after the initial transient period. While there are several ways to quantify loading, this study provides a demonstration of the optimization framework by summarizing the time history of the flapwise bending moment of one blade in the front and back turbines after the same initial transient period. Using the high-fidelity model, loading is computed as the sum of damage equivalent loads (DELs) (International Electro-technical Commission, 2015) associated with the front and back turbine flapwise bending moments. The power and flapwise bending moment are computed from the actuator force as

$$\begin{array}{}\text{(26)}& {\mathcal{P}}^k\left(t\right)=\mathit{\omega }\sum _{b=\mathrm{1}}^{\mathrm{3}}\sum _{j=\mathrm{1}}^{N_{\mathrm{nodes}}}{r}_j\left({\mathit{f}}_{kbj}\cdot {\widehat{\mathit{n}}}_{\mathit{\theta },k}\right)\end{array}$$

and

$$\begin{array}{}\text{(27)}& M^k\left(t\right)=\sum _{j=\mathrm{1}}^{N_{\mathrm{nodes}}}{r}_j\left({\mathit{f}}_{k\mathrm{1}j}\cdot {\widehat{\mathit{n}}}_{n,k}\right),\end{array}$$

where 𝒫^k(t) is the power associated with turbine k, M^k(t) is the flapwise bending moment associated with turbine k, ω is the angular speed of the rotor (which is a constant 11.6 rotations per minute in this study), r_j is the radial location associated with node j, ${\widehat{\mathit{n}}}_{n,k}$ is the unit vector orientated outward from the rotor plane associated with turbine k, and ${\widehat{\mathit{n}}}_{\mathit{\theta },k}$ is the unit vector oriented in the direction of rotation of turbine k.

The average power production of each turbine is computed as

$$\begin{array}{}\text{(28)}& {\mathcal{P}}^{k,\mathrm{avg}}={\displaystyle \frac{\mathrm{1}}{t_{\mathrm{f}}-t_{\mathrm{0}}}}\underset{t_{\mathrm{0}}}{\overset{t_{\mathrm{f}}}{\int }}{\mathcal{P}}^k\left(t\right)\mathrm{d}t,\end{array}$$

where 𝒫^k,avg is the average power associated with turbine k, t₀ is the initial time considered, and t_f is the final time of the data set. The total power, measured in megawatts, is computed as the sum of the powers produced by each wind turbine,

$$\begin{array}{}\text{(29)}& P=\sum _{k=\mathrm{1}}^{\dim \left(\mathit{\gamma }\right)}{\mathcal{P}}^{k,\mathrm{avg}}.\end{array}$$

Each DEL is computed using the rainflow-counting algorithm as

$$\begin{array}{}\text{(30)}& \text{DEL}\left(M^k\right)={\left(\sum _{i=\mathrm{1}}^{\mathrm{100}}{R}_i^m{\displaystyle \frac{c_i}{\mathrm{\Delta }t}}\right)}^{\mathrm{1}/m},\end{array}$$

where i loops through each cycle found using the rainflow-counting algorithm, R_i is a load range, c_i is the number of cycles associated with the ith range bin of the moment load spectrum, Δt is the time elapsed in seconds, and m is the Wöhler exponent. In this study, m is set as 10, and R_i and c_i are computed using the fatpack Python package (Frøseth and Capponi, 2021), utilizing 100 loading bins. The loading objective is computed as the sum of the flapwise bending moment DELs associated with the front and back turbines as

$$\begin{array}{}\text{(31)}& {\widehat{L}}^{\mathrm{HF}}=\sum _{k=\mathrm{1}}^{\dim \left(\mathit{\gamma }\right)}\text{DEL}\left(M^k\right).\end{array}$$

The loading objective, L, is normalized to be negative and on a similar scale to power, as

$$\begin{array}{}\text{(32)}& L=\widehat{L}/{\mathrm{10}}^{\mathrm{7}}-\mathrm{10},\end{array}$$

where $\widehat{L}$ is the load prior to normalization, which is computed in newton meters. This normalization was chosen to ensure that power and loading are of similar scale and that both are negative. Because the EHVI is an area produced by the two objectives, we do not expect different values in this scaling function to affect the results of maximizing the acquisition function, provided that all sampled objective values are always less than the associated reference value.

3.3 Optimization implementation

Here we use a simple optimization approach for simplicity of demonstration. In each iteration, the l correlation scale parameter in Eq. (4) is selected based on the maximum likelihood function (Pedregosa et al., 2011), with a lower bound of 5^∘ and an upper bound of 30^∘. The reference point, r, in Eq. (12) is specified as $(\mathrm{0}{}^{\circ },\mathrm{0}{}^{\circ })$. In our formulation, we minimize J(γ,l) using a grid-based search for the maximum value. The grid is evenly spaced with 31 inputs per yaw offset dimension. Individual grids are considered for all values of l. The EHVI is computed using Monte Carlo sampling with 1200 samples taken from the GP. The EHVI may also be computed through numerical quadrature (Emmerich et al., 2011; Hupkens et al., 2015). When dealing with inputs of larger dimension, the EHVI may be maximized using an optimization algorithm, such as a genetic algorithm or the EGO approach. After the optimization, the Pareto set was refined using B-spline interpolation in SciPy (Virtanen et al., 2020).

3.3.1 Initial sampling

Initial sampling points are selected using a heuristic approach, where an assumed kernel is used to progressively minimize the standard deviation of the predictor. The simplest approach to initializing the optimization procedure is to randomly sample the low-fidelity and discrepancy functions. Random initial sampling may affect the optimization results, so a deterministic and symmetric sampling strategy is used as a test case. An isotropic kernel is used with a correlation scale of 10^∘. The GP model is initialized with the point $(\mathrm{0}{}^{\circ },\mathrm{0}{}^{\circ })$. A 100×100 grid of inputs ranging between −30 and 30^∘ degrees is used to find the next point that minimizes uncertainty in the predicted variance. This process was repeated iteratively to generate 100 points. These points were used to naively estimate the optimal power and loading and Pareto hypervolume as a reference. The first 5 points were used as the initial high-fidelity samples, and the first 20 points were used as the initial low-fidelity samples. This heuristic sampling approach is also used to generate 100 samples for use in a correlation analysis and as a naive, baseline approach to searching for the Pareto set.

3.3.2 Low-fidelity loading model

We considered several different low-fidelity model forms for loading and selected the one with the highest correlation to the high-fidelity model, as that is known to result in the best multifidelity performance. We used 100 samples obtained using the heuristic sampling method described in Sect. 3.3.1 to test the correlation. In practice, this correlation test would not be part of the optimization procedure, and reasonably accurate models would be identified based on past experience and/or expert opinion. The correlation analysis revealed a correlation of 0.976 between the low- and high-fidelity power predictions. Using the DEL as the load proxy in the low-fidelity model yielded a low correlation between the low- and high-fidelity models. We explored other potential loading proxies – applying the proxy to both the front and back turbine moment histories and then summing the results – and the results of each correlation analysis are presented in Table 1. We used the proxy associated with the highest-measured correlation, namely

$$\begin{array}{}\text{(33)}& \begin{array}{rl}{\widehat{L}}^{\mathrm{LF}}=& {\mathit{\mu }}_t\left[M^{\mathrm{front}}\left(t\right)\right]+\mathrm{5}{\mathit{\sigma }}_t\left[M^{\mathrm{front}}\left(t\right)\right]\\ & +{\mathit{\mu }}_t\left[M^{\mathrm{back}}\left(t\right)\right]+\mathrm{5}{\mathit{\sigma }}_t\left[M^{\mathrm{back}}\left(t\right)\right],\end{array}\end{array}$$

where μ_t and σ_t are the mean and standard deviation, respectively, with respect to time. The DEL is purposely replaced with lower-order moment functions to avoid the influence of the spurious oscillations caused by the low-fidelity loading model. The DEL is essentially a high-order moment, which is especially susceptible to these oscillations. The lower-order moments in Eq. (33) were less susceptible to the spurious oscillations, which is why larger correlations are observed.

Table 1 Correlations observed between high-fidelity DEL and different loading proxies of the low-fidelity model using 100 heuristic samples.

4 Results and discussion

4.1 Pareto set computation

The convergence of the single-fidelity and multifidelity optimization approaches is compared in Fig. 3. The dashed lines show the hypervolume, best-sampled load, and best-sampled power found from 100 sampled points using the heuristic sampling approach, which the single-fidelity and multifidelity approaches both outperformed. The EHVI associated with the multifidelity approach was generally lower than the EHVI associated with the single-fidelity approach. The multifidelity approach took less than a third as much total time to estimate the optimal power and loading compared to the single-fidelity approach. The optimal power is achieved with 26^∘ offset in the front turbine and 2^∘ offset in the back turbine. The optimal loading is achieved with 22^∘ offset in the front turbine and −30^∘ offset in the back turbine. We also performed several shorter optimizations as part of the development process using random initial samples, confirming that the multifidelity approach consistently determined the correct hypervolume in fewer iterations than the single-fidelity approach.

Figure 3 Convergence history of the single-fidelity and multifidelity approaches. The left plots show the EHVI, hypervolume, and best-observed power and loading. The right plots show the yaw configurations associated with the best-observed power and loading.

Figure 4 compares the Pareto sets found using 50 equivalent high-fidelity evaluations using the single-fidelity and multifidelity approaches. The multifidelity approach captures five Pareto points, and the single-fidelity approach captures four Pareto points. The estimated Pareto sets are very similar, although the single-fidelity approach captures more of the Pareto set close to the optimal power and the multifidelity approach captures more of the Pareto set close to the optimal loading. The results of the single-fidelity and multifidelity approaches are combined to show a single Pareto set, which has a shape similar to a logit function.

Figure 4 Sampled inputs and outputs associated with power greater than 3 MW and loads less than 3 MN m. Points associated with the single-fidelity approach are shown using blue square markers, and points associated with the multifidelity approach are shown using green circular markers. A Pareto set constructed from the single-fidelity and multifidelity results is highlighted with hollow circles, where darker (magenta) circles correspond to Pareto points with lower loads and lower powers.

Although the primary goal of the present study is to develop and demonstrate a multifidelity multiobjective optimization framework, once several of the points in the Pareto set are identified the set can be further refined using a grid search. As a demonstration of this additional step, the Pareto set resulting from the combination of the single-fidelity and multifidelity approaches was interpolated to create refinement samples using B-spline interpolation (Virtanen et al., 2020) with 10 interpolation points. To pick up more of the Pareto set, this interpolation was offset in the γ₁ direction by −2^∘, −1^∘, 1^∘, and 2^∘ when creating the input refinement set.

The Pareto set resulting from these additional refinement samples is visualized in Fig. 5. The resulting Pareto set has three more points than the Pareto set found combining the single-fidelity and multifidelity approaches. The optimization algorithm did not originally fill in these points because they reside in a relatively flat portion of the Pareto set (i.e., $\mathrm{d}P/\mathrm{d}L$ is small), where adding points would not be expected to increase the Pareto set hypervolume based on the rectangular quadrature employed in this study. Adding these points increased the hypervolume of the discovered Pareto set by a very small amount, on the order of 0.002 %. Even when using such a refined set of inputs, there are several points where the Pareto set jumps from one yaw position to another; attention would be needed for an operator to control the turbines to operate along the Pareto set.

Figure 5 Sampled inputs and outputs associated with power greater than 3 MW and loads less than 3 MN m. Points associated with the single-fidelity approach are shown using blue square markers, and points associated with the multifidelity approach are shown using green circular markers. Refinement points are shown as black triangles. A Pareto set constructed from the single-fidelity, multifidelity, and refinement samples is highlighted with hollow circles, where darker (magenta) circles correspond to Pareto points with lower loads and lower powers.

The multifidelity approach was successful in quantifying the trade-offs between loading and power and was shown to be more efficient than its single-fidelity counterpart. From the presented results, we find that loading may be reduced by 4 % while only reducing the optimal power by 0.3 %. The accuracy of the single-fidelity and multifidelity GP models is quantified using a leave-one-out analysis in the Appendix A. Table 2 shows the power and front and back turbine DELs associated with several strategies. Slight adjustments to the back turbine angle result in substantial differences in the back turbine loading. These small changes in yaw position adjust the turbine thrust away from the flow, reducing the total thrust imparted on the back turbine.

Table 2 Observed power and loads for various yaw configurations.

4.2 Flow physics insights

Figure 6 shows the flow fields associated with neutral, −30^∘, and +30^∘ yaw offsets in the front turbine, with the back turbine aligned with the wind direction. When the front turbine is offset, two structures are produced: a pair of counter-rotating vortices and a coherent structure that is drawn from the boundary layer. The direction of vortex rotation is determined by the direction of thrust the turbine imparts on the incoming air. Induced vortices generally rotate in the opposite direction from the blades that generated them, and the location of the vortices is determined by their rotational direction and the direction of blade rotation. The upper vortex associated with the positive yaw offset is lower in elevation than the upper vortex produced by the negative yaw offset. The upper vortex also drifts less in the cross-flow direction when using the positive yaw offset than when using a negative yaw offset. The bottom vortex drifts similarly in both the positive and negative offset cases. All this amounts to a larger and more extreme velocity deficit encroaching on the back turbine when using the negative yaw offset, rather than the positive yaw offset, resulting in more loading and less power.

Figure 6 Flow fields associated with the extreme and neutral offsets in the front turbine, viewed from upstream. The Δx term indicates the distance downstream from the front turbine in terms of rotor diameters. Brighter colors show faster velocities. Streamlines show the direction of the cross-flow and vertical velocity components. In each plot, the vertical and cross-flow location of the back turbine is shown as a white circle. The turbines rotate clockwise when viewed from upstream. A diagram is shown on the left depicting the direction of positive yaw offset when viewing the turbine from above.

Time-averaged flow fields associated with the optimal power, $\mathit{\gamma }=(\mathrm{26}{}^{\circ },\mathrm{2}{}^{\circ })$, and optimal loading, $\mathit{\gamma }=(\mathrm{22}{}^{\circ },-\mathrm{30}{}^{\circ })$, solutions are compared in Fig. 7, which shows vertical slices of the flow field before the flow reaches the back turbine. The lesser front turbine yaw offset angle in the $\mathit{\gamma }=(\mathrm{22}{}^{\circ },-\mathrm{30}{}^{\circ })$ case results in less lateral movement of the wake, and the wake structure has greater overlap with the back turbine than in the $\mathit{\gamma }=(\mathrm{26}{}^{\circ },\mathrm{2}{}^{\circ })$ case. The stronger vortical motion resulting from the $\mathit{\gamma }=(\mathrm{26}{}^{\circ },\mathrm{2}{}^{\circ })$ case results in the boundary layer being convected further inwards. This boundary layer structure also impacts the back turbine less in the $\mathit{\gamma }=(\mathrm{22}{}^{\circ },-\mathrm{30}{}^{\circ })$ case because of the reduced back turbine projected area. As the wake convects past the back turbine, additional vorticity is added to the flow. In the $\mathit{\gamma }=(\mathrm{22}{}^{\circ },-\mathrm{30}{}^{\circ })$ case, the boundary layer structure appears to be pushed back down by the rotation of the bottom vortex.

Figure 7 Time-averaged flow fields associated with the optimal power (left) and loading (right) found by the optimization, viewed from upstream, six rotor diameters away from the front turbine and one rotor diameter upstream of the downstream turbine. Brighter colors indicate faster velocity magnitudes. In each plot, the vertical location and cross-flow location of the back turbine are shown as a white ellipse.

Figure 8 shows time histories of the flapwise bending moment associated with the front and back turbines for various wake-steering strategies. The spikes in the back turbine loading history are caused by the wake impacting the back turbine. When $\mathit{\gamma }=(-\mathrm{30}{}^{\circ },\mathrm{0}{}^{\circ })$, the back turbine shows greater downward spikes in loading than when $\mathit{\gamma }=(\mathrm{30}{}^{\circ },\mathrm{0}{}^{\circ })$, because of the greater velocity deficit discussed above. The $\mathit{\gamma }=(\mathrm{26}{}^{\circ },\mathrm{0}{}^{\circ })$ case yields smaller downward spikes associated with the back turbine loading than when $\mathit{\gamma }=(\mathrm{30}{}^{\circ },\mathrm{0}{}^{\circ })$, because the strength of the counter-rotating vortices is such that the structure convected from the boundary layer does not impact the back turbine as adversely. The $\mathit{\gamma }=(\mathrm{22}{}^{\circ },\mathrm{0}{}^{\circ })$ offset case has larger downward spikes in the back turbine loading than the $\mathit{\gamma }=(\mathrm{26}{}^{\circ },\mathrm{0}{}^{\circ })$ offset case, because the latter steers the wake further away from the back turbine. When the back turbine is offset to −30^∘, the back turbine thrust and associated moments are generally reduced. With this extreme back turbine yaw offset, there is less variation in the back turbine loading when the front turbine is offset by 22^∘ than when it is offset by 26^∘, because the former case results in greater variation of velocity across the back turbine rotor plane. These results are specific to the specified spacing between turbines and atmospheric conditions.

Figure 8 Loading histories associated with different yaw offset strategies (values of γ).

5 Conclusions

This paper has demonstrated a multifidelity multiobjective optimization approach for wake-steering strategies. Actuator line simulations were carried out using the WindSE tool and using a coarser simulation as the low-fidelity model. The high-fidelity loading was characterized as the sum of flapwise bending moment DELs on blades on the front and back turbines. Due to oscillations in the low-fidelity simulations, characterizing the low-fidelity loading with a DEL resulted in a relatively low correlation between the low- and high-fidelity loading predictions, so a different low-fidelity surrogate was developed with a higher correlation.

The multifidelity multiobjective optimization approach was effective in exploring the trade-offs between loading and power when developing a wake-steering design. Convergence was achieved in the multifidelity optimization case after approximately 30 % as many equivalent high-fidelity model evaluations as in the single-fidelity case. Future work should apply this approach and a low-fidelity loading function to more complex wind plant layouts to confirm their effectiveness. Exploring the solutions in the final Pareto sets guided insights into the fundamental flow physics. Given the specified turbine spacing and atmospheric conditions, a positive front turbine yaw offset is more effective at reducing loading and increasing power than a negative yaw offset because the counter-rotating vortices associated with the negative front turbine yaw offset produce a greater velocity deficit in the downstream wake. The boundary layer is convected by the counter-rotating vortices, adversely affecting loading, and this may be avoided using less-extreme front turbine yaw offsets. Slightly modifying the back turbine yaw offset reduced loading by 4 % and only reduced power by 0.3 %. Greater offsets in the back turbine also led to less overall loading, with significantly less power generation.

It is well known that yaw position errors can adversely affect the performance of wake-steering strategies. This is especially true when it comes to turbine loading. A 30^∘ yaw offset is already an aggressive strategy, and unfavorable yaw position errors may result in even more aggressive yaw offsets in practice. Yaw offset errors are generally extreme in lower wind speeds, which is when wake-steering strategies are most efficient at increasing power. Previous work has examined the potential of considering yaw error uncertainties in the wake-steering optimization problem (Quick et al., 2017, 2020). The multifidelity optimization approach presented in this paper could conceivably be extended to optimization under uncertainty, using the final GP models to propagate yaw position uncertainty, and potentially even modifying the EHVI definition to include uncertainty information. Incorporating uncertainty will likely change the shape of the discovered Pareto front.

A drawback of the presented approach is that it requires sequential high-fidelity model evaluations. In practice, it is often feasible to evaluate a high-fidelity model several times in parallel, and the greatest expense is the time needed to run the optimization. This framework may be extended to allow for parallel function evaluations. A simple approach is to use predictions of the GP as stand-ins for future model evaluations, iteratively using these points to construct the next iteration of the GP and the associated EHVI (Ginsbourger et al., 2010). Yang et al. (2019) propose dividing the input space into separate regions for parallelization of EHVI optimization. Another intuitive approach could be to include refinement points during each iteration. Refinement points could be selected using the Pareto set predicted by the GP models or interpolated along the observed Pareto set. Care should also be taken when applying this method to ensure convergence of the Pareto set with respect to convergence of the underlying simulation.

In future work, this framework can be applied to a larger array of turbines using more realistic control strategies with different turbine spacings and atmospheric conditions. While considering more turbines presents additional complications in maximizing the EHVI, we anticipate there will be even greater cost savings from the multifidelity approach as the number of turbines increases. Additionally, the framework can be extended to allow for optimization under uncertainty, as it is not realistic to assume perfect control of wind turbine yaw positions. Finally, the framework can incorporate more lower-fidelity models and be combined with layout optimization to realize the full benefits of multifidelity multiobjective wake-steering optimization.

Appendix A

A leave-one-out analysis was performed to assess the accuracy of the single-fidelity and multifidelity GP models. For each point considered, the model was trained using the remaining points available, excluding the point of interest. Then, the accuracy of the prediction was quantified by comparing it to the observed value. The single-fidelity and multifidelity approaches were analyzed using the data associated with the optimization case study. When assessing the accuracy of the multifidelity model, the low- and high-fidelity samples associated with each point considered in the leave-one-out analysis were removed. Figure A1 shows results of the leave-one-out analysis associated with the single-fidelity GP. Figure A2 shows the results of the leave-one-out analysis associated with the multifidelity GP.

Figure A1 Results of the single-fidelity leave-one-out analysis. Panel (a) shows the leave-one-out prediction errors associated with power and loading, and the points are colored by the sum of both errors. The same points are plotted in panel (b), showing their associated γ₁ and γ₂ values.

Figure A2 Results of the multifidelity leave-one-out analysis. Panel (a) shows the leave-one-out prediction errors associated with power and loading, and the points are colored by the sum of both errors. The same points are plotted in panel (b), showing their associated γ₁ and γ₂ values.

Based on these results, both GPs served as reasonably accurate surrogates. Many of the sampled errors are less than 0.1 MW and 0.1 MN m, particularly in the region of the discovered Pareto set, which correspond to 3 % of the maximum power and 6 % of the minimum loading, respectively. The multifidelity approach yielded higher maximum errors and lower minimum errors than the single-fidelity approach. This analysis focused on the final results of the optimization, and we generally expect the leave-one-out errors to shrink as the optimization progresses.