# Data-Driven Wind Turbine Wake Modeling via Probabilistic Machine Learning

S. Ashwin Renganathan
Mathematics & Computer Science
<EMAIL_ADDRESS>
&Romit Maulik
Mathematics & Computer Science
<EMAIL_ADDRESS>
&Stefano Letizia
Department of Mechanical Engineering
The University of Texas at Dallas,
Dallas, TX
<EMAIL_ADDRESS>
&Giacomo Valerio Iungo
Department of Mechanical Engineering
The University of Texas at Dallas,
Dallas, TX
<EMAIL_ADDRESS>
Argonne National LaboratoryArgonne National Laboratory

###### Abstract

Wind farm design primarily depends on the variability of the wind turbine wake
flows to the atmospheric wind conditions, and the interaction between wakes.
Physics-based models that capture the wake flow-field with high-fidelity are
computationally very expensive to perform layout optimization of wind farms,
and, thus, data-driven reduced order models can represent an efficient
alternative for simulating wind farms. In this work, we use real-world light
detection and ranging (LiDAR) measurements of wind-turbine wakes to construct
predictive surrogate models using machine learning. Specifically, we first
demonstrate the use of deep autoencoders to find a low-dimensional _latent_
space that gives a computationally tractable approximation of the wake LiDAR
measurements. Then, we learn the mapping between the parameter space and the
(latent space) wake flow-fields using a deep neural network. Additionally, we
also demonstrate the use of a probabilistic machine learning technique,
namely, Gaussian process modeling, to learn the parameter-space-latent-space
mapping in addition to the epistemic and aleatoric uncertainty in the data.
Finally, to cope with training large datasets, we demonstrate the use of
variational Gaussian process models that provide a tractable alternative to
the conventional Gaussian process models for large datasets. Furthermore, we
introduce the use of active learning to adaptively build and improve a
conventional Gaussian process model predictive capability. Overall, we find
that our approach provides accurate approximations of the wind-turbine wake
flow field that can be queried at an orders-of-magnitude cheaper cost than
those generated with high-fidelity physics-based simulations.

##  1 Introduction

Understanding and modeling wind farm flows still represent major challenges
for wind farm designers and operators. Important aspects, such as the
interaction of the atmospheric boundary layer with wind turbines [1, 2] and
prediction of the wake morphology and their superposition [3], are far from
being fully understood, in spite of having been the object of numerous
insightful scientific studies.

The term turbine wake refers to the low momentum and highly turbulent region
located downstream of an operating wind turbine, which is a direct consequence
of the extraction of kinetic energy from the incoming wind field. Significant
power losses [4, 5, 6] and enhanced fatigue loads [7, 8] were documented for
turbines impinged by upstream wakes. The study of turbine wakes through
numerical simulations is encumbered with difficulties due to the high Reynolds
number flow (which entails a great span of length and time scales involved)
[3], the unsteadiness of the inflow conditions [9, 10], and the relevant role
of atmospheric stability [11]. Recently, large-eddy simulations (LES) have
become a well-established tool for the simulation of wind farm flows [12, 13,
14]. However, their computational costs are still prohibitive for large scale
tasks, such as for layout optimization and real-time performance diagnostics.

Reynolds-averaged Navier-Stokes [15, 16] and engineering wake models (e.g.
[17, 18, 19]) provide faster and simpler alternatives to LES, since they do
not explicitly solve the unsteady small-scale turbulent eddies. However, while
the former are oftentimes affected by inaccuracy due to the turbulence closure
[15], the latter require a thorough calibration to achieve a satisfactory
agreement with experimental data [20].

The above-mentioned challenges have spurred the interest of wind energy
scientists in the experimental characterization of the wind farm flow. In
particular, the improvements of remote sensing instruments, such as wind Light
Detection and Ranging (LiDAR), have promoted the proliferation of field
experimental campaign investigating the wakes of utility-scale wind turbines
(e.g., see [9, 21, 20]). These studies have highlighted the great complexity
and sensitivity to environmental conditions on the characteristics and
downstream evolution of wind turbine wakes.

Following up on past work [22], in this work, LiDAR measurements of individual
wakes generated by utility-scale wind turbines under broad ranges of
atmospheric and wind conditions are leveraged to develop data-driven machine
learning models to enable accurate predictions of the wake velocity field for
prescribed wind/atmospheric conditions, while requiring computational costs as
low as for empirical wake engineering models. Specifically, we explore the
application of deep neural networks (DNN) for efficient data reduction via
autoencoders, as well as to build predictive models via multilayer
perceptrons. Furthermore, we also explore the combination of the data
reduction via the DNN and Gaussian process (GP) models to develop predictive
probabilistic models of the wake flow. This way, we show how probabilistic
machine learning can be leveraged to learn predictive wake models from noisy
and incomplete measurement data. To summarize, the contributions of this
article are:

  1. 1.

A novel convolutional autoencoder framework to obtain low-dimensional
embeddings of wind LiDAR measurements for wind-turbine wakes. In this manner,
we develop a compressed—and hence tractable—representation of the high-
dimensional LiDAR data.

  2. 2.

A DNN to learn the mapping between the input parameter space and the latent
space derived from the convolutional autoencoder. Therefore, we learn the
mapping between the input parameters and the spatially distributed wake
measurements, via two levels of DNNs–one for data reduction and one for
prediction. This serves as a cheap-to-evaluate surrogate model to predict wind
turbine wake velocity fields.

  3. 3.

Additionally, similar to item 2, we propose the use of probabilistic machine
learning model via GP regression to learn latent-space to wake flow field
mapping, which can simultaneously learn the noise in the data. To address the
tractability issues associated with GP regression with _big_ data, we also
investigate the use of _variational_ and active learned GP regression.

Overall, this work brings together state-of-the-art data-driven machine
learning and the classic field of wind energy to build cheap and reliable
surrogate models that can be leveraged to perform exploratory studies. fig. 1
provides a high-level summary of the work.

Figure 1: Summary of the overall methodology.

The remainder of this article is organized as follows. In section 2, we
present the details of our experimental setup and data collection via LiDAR.
In sections section 3 and section 4 we present the theoretical details of our
machine learning methods. We present the results and associated discussion in
section 5, followed by concluding remarks and an outlook on future work in
section 6.

##  2 Experimental data and collection methods

The wind LiDAR measurements used in this work were collected in the period
from August 2015 to March 2017 at a wind farm located in North Texas 111The
original LiDAR data used in this work are available upon reasonable request
from the fourth author, who may be contacted at<EMAIL_ADDRESS>The wind farm includes 25 identical wind turbines with a nameplate capacity of
2.3 MW. The rotor diameter is $d=127$ m and the height of the hub is 89 m
above the ground. The local topographic map provided by [23] with a resolution
of 100 m shows that 95% of the terrain within the farm has a slope lower than
$3^{\circ}$, which allows to rule out effects due to the terrain on the flow.

A WindCube 200S scanning pulsed Döppler LiDAR was installed at turbine 11 (see
fig. 2) and probed the wakes stemming from turbines 01 to 06 during the
occurrence of southerly winds. Plan Position Indicator (PPI) scans were
scheduled targeting the wake of the turbine with the best line-of-sight
alignment with the LiDAR, which enables achieving an optimal spatio-temporal
sampling resolution.

Figure 2: Map of the experimental site. On the left, is the topographic map
and farm layout, where the green arcs indicate wind sectors where no relevant
wake interactions are expected. On the right, it is the directional histogram
of wind speed and turbulence intensity based on the meteorological data of the
80-meters tower at the "Met" location. $U_{norm}$ refers to the hub-height
wind speed normalized by the rated wind speed of the turbines (11 m s-1). Grey
sectors are likely affected by wakes.

Furthermore, meteorological and SCADA data were continuously collected for the
whole duration of the campaign in the form of 10-minute mean and standard
deviation of wind speed, wind direction, temperature, atmospheric pressure,
active power, RPM, and blade pitch angle. The list of input parameters is
presented in table 1; for a more detailed description of the experimental site
and LiDAR scanning strategy please refer to [5, 11].

The wind rose of the site (fig. 2) indicates high likelihood of southerly
winds and a negligible directional dependence of the turbulence intensity at
hub-height, which confirms the findings of previous studies [5, 11] on the
prevalent thermally-driven cycle of the atmospheric boundary layer caused by
the diurnal variation of the solar irradiance.

The LiDAR data undergo a quality control which excludes all the points
characterized by a carrier-to-noise ratio lower than $-25$ dB. Subsequently,
the wake data are realigned with the wind direction, which is estimated as the
10-minute moving averaged wake direction. This also allows estimating the
horizontal equivalent velocity as:

$u_{eq}\sim\frac{u_{\text{LOS}}}{\cos(\theta-\theta_{w})\cos\beta},$ (1)

where $u_{\text{LOS}}$ is the line-of-sight velocity measured by the LiDAR,
$\theta_{w}$ is the wind direction, $\theta$ and $\beta$ are the azimuth and
elevation angles of the LiDAR, respectively. The vertical variability of the
incoming wind due to wind shear is corrected by normalizing the equivalent
velocity by the vertical undisturbed velocity profile. After the outlined
post-processing, 6,654 quality-control re-aligned and non-dimensional LiDAR
scans are made available for the following analysis. For more technical
details on this procedure, the interested reader shall refer to [11].

Parameter | Range | Description
---|---|---
SCADA_WS [m/s] | $[2.92,15.22]$ | Mean Hub-height wind speed recorded by the SCADA
MET_WS_80m [m/s] | $[3.9,15.22]$ | Mean Hub-height wind speed recorded by the met-tower
SCADA_TI [-] | $[0.04,0.36]$ | Hub-height wind turbulence intensity recorded by the SCADA
MET_BulkRichardson [-] | $[-0.01,0.01]$ | Bulk Richardson number recorded by met-tower
SCADA_Power [kW] | $[58.8,2423]$ | Mean power capture recorded by the SCADA
SCADA_RPM [RPM] | $[7.07,16.95]$ | Mean rotor rotational speed recorded by the SCADA
SCADA_Pitch [deg] | $[-2,80]$ | Mean blade pitch angle recorded by the SCADA
Table 1: List of input parameters and their ranges.

##  3 Deep learning for parametric flow prediction

In the following section, we introduce our deep neural network architectures
for establishing a viable emulation strategy for data obtained from LiDAR
measurements.

###  3.1 Convolutional autoencoder

Autoencoders are neural networks that learn a new representation of the input
data, usually with lower dimensionality. The initial layers, called the
_encoder_ , map the input $\mathbf{x}\in\mathbb{R}^{n}$ to a new
representation $\tilde{\mathbf{x}}\in\mathbb{R}^{k}$ with $k\ll n$. The
remaining layers, called the _decoder_ , map $\tilde{\mathbf{x}}$ back to
$\mathbb{R}^{n}$ with the goal of reconstructing $\mathbf{x}$. The objective
is to minimize the reconstruction error. Autoencoders are unsupervised; the
data $\mathbf{x}$ is given, but the representation $\tilde{\mathbf{x}}$ must
be learned.

More specifically, we use autoencoders that have convolutional layers. In a
convolutional layer, instead of learning a matrix that connects all $m$
neurons of layer’s input to all $n$ neurons of the layer’s output, we learn a
set of filters that are convolved with regions of the layer’s input. Suppose a
one-dimensional (1-d) convolutional layer has filters of length $L_{f}$, then,
each of the layer’s output neurons corresponding to a specific filter
$\mathbf{f^{i}}$ is connected to a patch of $L_{f}$ of the layer’s input
neurons. In particular, a 1-d convolution of filter $\mathbf{f^{i}}$ and patch
$\mathbf{p}$ is defined as
$\mathbf{f^{i}}\ast\mathbf{p}=\sum_{j}f^{i}_{j}p_{j}$ (where $f^{i}_{j}$
corresponds to the stencil coefficient in the filter for index $j$). In other
words, convolutional neural networks identify stencil values $f^{i}_{j}$ that
obtain coherent translationally invariant features relevant to a particular
function approximator. Then, for a typical 1-d convolutional layer, the
layer’s output neuron $y_{ij}=\varphi(\mathbf{f^{i}}\ast\mathbf{p_{j}}+B_{i})$
where $\varphi$ is an activation function, and $B_{i}$ are the entries of a
bias term. As $\mathbf{j}$ increases, patches are shifted by stride $s$. For
example, a 1-d convolutional layer with a filter $\mathbf{f^{0}}$ of length
$m_{f}=3$ and stride $s=1$ could be defined so that $y_{0j}$ involves the
convolution of $\mathbf{f^{0}}$ and inputs $\mathbf{j}-1,\mathbf{j}$, and
$\mathbf{j}+1$. To calculate the convolution, it is common to add zeros around
the inputs to a layer, which is called _zero padding_. In the decoder, we use
deconvolutional layers to return to the original dimension. These layers
upsample with nearest-neighbor interpolation.

Two-dimensional convolutions are defined similarly, but each filter and each
patch are two-dimensional. A 2-d convolution sums over both dimensions, and
patches are shifted both ways. For a typical 2-d convolutional layer, the
output neuron $y_{hij}=\varphi(\mathbf{f^{h}}\ast\mathbf{p_{ij}}+B_{h})$.
Input data can also have a “channel” dimension, such as Red/Green/Blue values
for images. The convolutional operator sums over channel dimensions, but each
patch contains all of the channels. The filters remain the same size as
patches, so they can have different weights for different channels. It is
common to follow a convolutional layer with a _pooling_ layer, which outputs a
sub-sampled version of the input. In this paper, we specifically use max-
pooling layers. Each output of a max-pooling layer is connected to a patch of
the input, and it returns the maximum value in the patch.

Autoencoders have recently become popular for the nonlinear dimensionality
reduction of datasets extracted from several high dimensional systems. These
have been motivated by the extraction of coherent structures that parameterize
low-dimensional embeddings in manifolds [24, 25, 26], and the utilization of
these embeddings for efficient surrogate models of nonlinear dynamical systems
[27, 28, 29, 30, 31, 32, 33, 34]. In this work, we utilize convolutional
autoencoders to identify low-dimensional representations of experimentally
collected data for building parameter-observation maps where the former are
obtained through meteorological and wind turbine data and the latter are LiDAR
measurements collected in the wake generated by wind turbines.

###  3.2 Multilayered perceptron (MLP)

One technique to obtain a mapping from the meteorological and turbine datasers
and the latent space embeddings of the convolutional autoencoder is through
the use of a multilayered perceptron (MLP) architecture, which is a subclass
of feedforward artificial neural network. A general MLP consists of several
neurons arranged in multiple layers. These layers consist of one input and one
output layer along with several hidden layers. Each layer (with the exception
of an input layer) represents a linear operation followed by a nonlinear
activation that allows for great flexibility in representing complicated
nonlinear mappings. This may be expressed as

$\displaystyle\mathcal{L}^{l}\left(\boldsymbol{t}^{l-1}\right):=\boldsymbol{w}^{l}\boldsymbol{t}^{l-1}+\boldsymbol{b}^{l},$
(2)

where $\boldsymbol{t}^{l-1}$ is the output of the previous layer, and
$\boldsymbol{w}^{l},\boldsymbol{b}^{l}$ are the weights and biases associated
with that layer. The output $\boldsymbol{t}^{l}$, for each layer may then be
transformed by a nonlinear activation, such as rectified linear activation:

$\displaystyle\eta(a)=\text{ReLU}(a)=\max(a,0).$ (3)

For our experiments, the inputs $\boldsymbol{t}^{0}$ are in $\mathbb{R}^{d}$
(i.e., $d$ is the number of inputs) and the outputs $\boldsymbol{t}^{K}$ are
in $\mathbb{R}^{k}$ (i.e., $k$ is the number of outputs). The final map is
given by

$\begin{split}F:\mathbb{R}^{d}\mapsto\mathbb{R}^{k}&,\quad\boldsymbol{t}^{0}\mapsto\boldsymbol{t}^{K}=F(\boldsymbol{t}^{0};(\boldsymbol{w},\boldsymbol{b})),\\\
\text{where}\\\
F(\boldsymbol{t}^{0};\boldsymbol{w},\boldsymbol{b})=&\eta^{K}\left(\mathcal{L}^{K}\circ\eta^{K-1}_{\text{act}}\circ\mathcal{L}^{K-1}\circ\ldots\circ\eta^{1}_{\text{act}}\circ\mathcal{L}^{1}\right)(\boldsymbol{t}^{0})\end{split}$
(4)

is a complete representation of the neural network and where $\mathbf{w}$ and
$\mathbf{b}$ are a collection of all the weights and the biases of the neural
network. These weights and biases, lumped together as
$\phi=\\{\mathbf{w},\mathbf{b}\\}$, are trainable parameters of our map, which
can be optimized by examples obtained from a training set. The supervised
learning framework requires for this set to have examples of inputs in
$\mathbb{R}^{N_{ip}}$ and their corresponding outputs $\mathbb{R}^{N_{op}}$.
This is coupled with a cost function $\mathcal{C}$, which is a measure of the
error of the prediction of the network and the ground truth. Our cost function
is given by

$\displaystyle\mathbb{C}=\frac{1}{|T|}\sum_{(\boldsymbol{\theta},\tilde{\mathbf{x}})\in
T}\left\|\tilde{\mathbf{x}}-F(\boldsymbol{\theta};\phi)\right\|^{2}$ (5)

with $|T|$ indicates the cardinality of the training data set given by

$T=\left\\{\left(\boldsymbol{\theta}_{i},\tilde{\mathbf{x}}_{i}\right):\tilde{\mathbf{x}}_{i}=f\left(\boldsymbol{\theta}_{i}\right)\right\\}.$

and where $f\left(\boldsymbol{\theta}_{i}\right)$ are examples of the true
targets obtained from the compressed training data using the autoencoder
introduced in the previous sections. Gradients of this cost function can then
be used in an optimization framework to obtain the best weights and biases,
given the training data. Finally, the trained MLP may be used for uniformly
approximating any continuous function on compact domains [35, 36], provided
$\eta(x)$ is not polynomial in nature.

##  4 Gaussian process regression

We now review the preliminaries of GP models. Our primary interest in the use
of GP models stems from its promise of offering enhanced data efficiency in
emulation compared to DNNs [37] as well as in sequential decision-making [38].
Although not pursued in this work, GPs offer greater potential in emulating
complex functions when combined with DNNs; e.g., see [39, 40]. GP models
provide a probabilistic approximation to an unknown function
$f(\boldsymbol{\theta})$. Specifically, $f(\boldsymbol{\theta})$ is assumed to
take the form of a GP, where each realization (or sample path) is a function.
This _prior_ assumption on the function can then be combined with the
probability of actual observations conditional on the prior (a.k.a, the
_likelihood_), using Bayes’ rule [41]. In this section, we provide a brief
overview of the theory behind GP models and highlight the difference between
exact and approximate inference, the latter finding applications in the
presence of large datasets.

###  4.1 Exact Gaussian process regression

We begin by placing a GP prior assumption on the unknown function, that is,
$f(\boldsymbol{\theta})\sim\mathcal{GP}(\mu(\boldsymbol{\theta}),k(\boldsymbol{\theta},\boldsymbol{\theta}^{\prime}))$,
where $\mu(\boldsymbol{\theta})$ is a mean function and
$k(\boldsymbol{\theta},\boldsymbol{\theta}^{\prime})$ is a covariance function
(or _kernel_), and $\boldsymbol{\Theta}\in\mathcal{T}$. We assume that we have
noisy observations of $f(\boldsymbol{\theta})$ of the form

$y_{i}=f(\boldsymbol{\theta})+\epsilon_{i},\quad i=1,\ldots,n$

where $\epsilon_{i}$ represents the observation noise. We assume
$\epsilon_{i}$ to take an independent and identically distributed normal
distribution with zero mean and a variance of $\sigma_{\epsilon}^{2}$, that
is, $\epsilon_{i}\sim\mathcal{N}(0,\sigma_{\epsilon}^{2})$. Furthermore, we
assume a Gaussian likelihood that defines the probability of observing the
data given the GP assumption.

Let $\mathbf{y}_{n}=[y_{1},\ldots,y_{n}]^{\top}$ denote the vector of noisy
observations of $f$ at
$\boldsymbol{\Theta}=[\boldsymbol{\theta}_{1},\ldots,\boldsymbol{\theta}_{n}]^{\top}$,
and $\mathcal{D}_{n}=\\{(\boldsymbol{\theta}_{i},y_{i}),~{}i=1,\ldots,n\\}$,
then, applying Bayes’ rule, the posterior distribution [42] is given by

$\begin{split}f(\boldsymbol{\theta})|\mathcal{D}_{n}\sim&\mathcal{GP}(\mu_{n}(\boldsymbol{\theta}),\sigma_{n}^{2}(\boldsymbol{\theta})),\\\
\textnormal{where~{}}\mu_{n}(\boldsymbol{\theta})=&\mathbf{k}^{\top}\mathbf{K}^{-1}(\mathbf{y}_{n}-\mu(\boldsymbol{\Theta}))\\\
\sigma_{n}^{2}(\boldsymbol{\theta})=&k(\boldsymbol{\theta},\boldsymbol{\theta})-\mathbf{k}^{\top}\mathbf{K}^{-1}\mathbf{k}.\end{split}$
(6)

In (6), $\mu_{n}$ and $\sigma_{n}^{2}$ are the mean and variance of the
posterior distribution, $\mathbf{K}\subset\mathcal{K}$ is the $n\times n$
covariance matrix with
$\mathbf{K}_{ij}=k(\boldsymbol{\theta}_{i},\boldsymbol{\theta}_{j}),\forall\boldsymbol{\theta}_{i},\boldsymbol{\theta}_{j}\in\boldsymbol{\Theta}$
and $\mathcal{K}$ being the cone of symmetric positive definite matrices,
$\mu({\boldsymbol{\Theta}})=[\mu(\boldsymbol{\theta}_{1}),\ldots,\mu(\boldsymbol{\theta}_{n})]^{\top}$,
and
$\mathbf{k}=[k(\boldsymbol{\theta},\boldsymbol{\theta}_{1}),\ldots,k(\boldsymbol{\theta},\boldsymbol{\theta}_{n})]^{\top}$.

The emulation properties of the GP are driven by the choice of the mean and
covariance functions. In this work, we standarize the observations
$\mathbf{y}_{n}$ such that they have a mean of zero; that is we set
$\mu(\boldsymbol{\theta})=0$ in the prior assumption. On the other hand, we
use a covariance function from the Matern class [43, 44, 45], given by

$k(\boldsymbol{\theta},\boldsymbol{\theta}^{\prime})=\frac{2^{1-\nu}}{\Gamma(\nu)}\left(\frac{\sqrt{2\nu}\|\boldsymbol{\theta}-\boldsymbol{\theta}^{\prime}\|}{\ell}\right)^{\nu}K_{\nu}\left(\frac{\sqrt{2\nu}\|\boldsymbol{\theta}-\boldsymbol{\theta}^{\prime}\|}{\ell}\right),$
(7)

with positive parameters $\nu$ and $\ell$, where $K_{\nu}$ is a modified
Bessel function, $\Gamma()$ is the Gamma function and $\|\cdot\|$ denotes the
Euclidean distance. The parameter $\nu$ (which we set to $3/2$) controls the
differentiability of the sample paths of the GP and is fixed, whereas $\ell$
is a _lengthscale_ parameter that controls the rate of change of the sample
paths in $\mathcal{T}$ and is a hyperparameter. The GP hyperparameter set
$\Omega=\\{\ell,\sigma_{\epsilon}^{2}\\}$ is estimated via a maximum
likelihood estimation (MLE) procedure frequently followed in fitting GP models
[42, 46].

The fact that the posterior distribution of the function given in (6) is
available in closed-form, the inference of such a model is called _exact_
inference. One of the main bottlenecks of the exact inference is that the
computation of the posterior mean and variance in (6) involves the inversion
of the matrix $\mathbf{K}$, whose computational cost scales as
$\mathcal{O}(n^{3})$, and hence gets expensive as $n$ increases. This
motivates the _approximate_ inference technique for GPs, namely, variational
GP regression, which trades some accuracy for large gains in computational
efficiency in fitting GP models when $n$ is large.

###  4.2 Approximate Gaussian process regression

In addition to the cost of estimating hyperparameters of the exact GP
involving $\mathcal{O}(n^{3})$ floating point operations, the prediction with
the GP costs $\mathcal{O}(n)$ and $\mathcal{O}(n^{2})$ operations,
respectively, for the posterior mean and variance computation. Sparse GP
models [47] circumvent this overhead by identifying a subset $m<n$ of the
training points, resulting in _reduced_ computational costs that scale as
$\mathcal{O}(m^{2}n)$, $\mathcal{O}(m)$, and $\mathcal{O}(m^{2})$ for fitting,
predicting mean, and predicting variance, respectively.

The choice of the subset of $m$ inducing points can be treated as another
hyperparameter and estimated by maximizing the marginal likelihood, just like
in the exact GP model; this results in an extended set of hyperparameters
$\\{\Omega,\bar{\boldsymbol{\Theta}}\\}$, where
$\bar{\boldsymbol{\Theta}}\in\mathbb{R}^{m\times d}$ are the _inducing
points_.

While sparse GPs bring down the cost of inverting the covariance matrix
$\mathbf{K}$ and predictions with the GP, the number of hyperparameters
increases (due to the addition of the inducing points), thereby making
inference computationally more expensive. To overcome this, we use variational
inference [48], which provides a tractable alternative to approximate unknown
probability densities in Bayesian models. Below, we briefly provide an
overview of sparse GP models and variational inference.

We now present a sparse model that is computationally tractable in terms of
inference and prediction with GPs. The sparsity arises because we consider a
sparse dataset $\bar{D}$ of size $m<n$ with inducing inputs
$\bar{\boldsymbol{\Theta}}=\\{\bar{\boldsymbol{\theta}_{i}},~{}i=1,\ldots,m\\}\subset\mathcal{T}$.
These inducing inputs can either be a subset of the training inputs or can be
randomly sampled from $\mathcal{T}$ [47]. Let $u=f(\bar{\boldsymbol{\theta}})$
be the inducing outputs, which are sampled from the same prior on the true
function $f$. In this work, we treat $\bar{\boldsymbol{\Theta}}$ as
hyperparameters and estimate them from maximizing the marginal likelihood.

Let the prior on the inducing outputs be given as

$p(\mathbf{u}|\bar{\boldsymbol{\theta}})\sim\mathcal{N}(0,\mathbf{K}_{mm}),$
(8)

where $\mathbf{K}_{mm}$ is the matrix of covariances between the inducing
inputs. This prior follows from the assumption that the inducing outputs also
behave like the latent variable $f$ which has a Gaussian prior. Therefore,
$\mathbf{u}$ and $\mathbf{f}$ have a joint Gaussian distribution [49] given by

$\begin{split}p(\mathbf{f},\mathbf{u})=&p(\mathbf{f|\mathbf{u}})p(\mathbf{u})\\\
=&\mathcal{N}(\mathbf{K}_{nm}\mathbf{K}^{-1}_{mm}\mathbf{u},\mathbf{K}_{nn}-\mathbf{Q}_{nn})\times\mathcal{N}(\mathbf{0},\mathbf{K}_{mm}),\end{split}$
(9)

where
$\mathbf{Q}_{nn}=\mathbf{K}_{nm}\mathbf{K}^{-1}_{mm}\mathbf{K}_{nm}^{\top}$.
To estimate the extended hyperparameter set
$\\{\Omega,\bar{\boldsymbol{\Theta}}\\}$, we first need to define a marginal
likelihood. In this case, the marginal likelihood is marginalized over
$\mathbf{f}$ and $\mathbf{u}$, that is,

$\begin{split}p(\mathbf{y}|\boldsymbol{\Theta},\Omega,\bar{\boldsymbol{\Theta}})=&\int
p(\mathbf{y}|\mathbf{u})p(\mathbf{u})d\mathbf{u}\\\ \end{split}$ (10)

Note that the density $p(\mathbf{y}|\mathbf{u})$ still involves inverting the
matrix of size $n\times n$ and hence is expensive to compute. Therefore, this
density is approximated using the evidence lower bound (ELBO) [48] as

$p(\mathbf{y}|\mathbf{u})\geq\mathbb{E}_{p(\mathbf{f}|\mathbf{u})}\left[p(\mathbf{y}|\mathbf{f})\right].$
(11)

Substituting (11) in (10), the lower bound on the marginal log likelihood can
be approximated as [49]

$p(\mathbf{y}|\boldsymbol{\Theta},\Omega,\bar{\boldsymbol{\Theta}})\geq\textnormal{log}\mathcal{N}(\mathbf{0},\mathbf{K}_{nm}\mathbf{K}^{-1}_{mm}\mathbf{K}_{nm}^{\top}+\sigma^{2}_{\epsilon}\mathbf{I})-\frac{1}{2}\textnormal{tr}(\mathbf{K}_{nn}-\mathbf{Q}_{nn}).$
(12)

Equation 12 can be maximized with respect to $\bar{\boldsymbol{\Theta}}$ and
$\Omega$ to estimate the hyperparameters, where the bound in (12) costs
$\mathcal{O}(nm^{2})$ for computation.

###  4.3 Active learning for Gaussian process regression

Whereas variational GPs provide an approximation to exact GP regression,
another approach to improving tractability of exact GP models is active data
selection [50, 51]. Specifically, given a training data set $\mathcal{D}_{n}$,
the data set is adaptively augmented as

$\begin{split}\mathcal{D}_{n+i}&:=\mathcal{D}_{n+i-1}\bigcup\\{(\boldsymbol{\theta}_{n+i},y_{n+i})\\},~{}i=1,\ldots,m\end{split}$
(13)

where $m$ is the number of adaptive model building steps. Furthermore, each
adaptive step can select a batch of $q$ training points jointly; when $q=1$ we
call the approach _sequential_ active learning and when $q>1$, _batch-
sequential_ active learning. At the end of the active learning process, the GP
is trained with a total of $n+mq$ training points. The main objective of
active learning is to choose points judiciously such that they are _optimal_
in the sense of improving model fit.

In this work, we choose points that are optimal in reducing the overall
uncertainty about the GP model. That is, we select points such that

$\boldsymbol{\theta}_{n+1}=\textnormal{arg}\max_{\boldsymbol{\theta}^{\prime}\in\mathcal{T}}~{}\int_{\mathcal{T}}-\sigma_{n+1}^{2}(\boldsymbol{\theta}^{\prime})~{}d\boldsymbol{\theta}^{\prime}$

where $\sigma_{n+1}^{2}(\boldsymbol{\theta})$ is the posterior variance of the
GP, having observed the $(n+1)$th point, and we introduce the negative sign to
solve a maximization problem. Essentially, we treat the integrated posterior
variance as a measure of uncertainty over our domain $\mathcal{T}$ and seek to
choose training data that are optimal in minimizing this uncertainty.

A fundamental issue with the above equation is that the posterior variance
$\sigma_{n+1}^{2}(\boldsymbol{\theta})$ is unknown until we actually commit to
$\boldsymbol{\theta}_{n+1}$ and choose a training point there. To circumvent,
we simulate the choice of $\boldsymbol{\theta}_{n+1}$ via the GP trained with
data $\mathcal{D}_{n+1}$, and choose $\boldsymbol{\theta}_{n+1}$ as the point
that reduces the _expected_ uncertainty:

$\begin{split}\boldsymbol{\theta}_{n+1}=&\textnormal{arg}\max_{\boldsymbol{\theta}^{\prime}\in\mathcal{T}}~{}-\int_{\mathcal{Y}}\int_{\mathcal{T}}\sigma_{n}^{2}(\boldsymbol{\theta})|\mathcal{D}_{n}\cup\\{\boldsymbol{\theta}^{\prime},y(\boldsymbol{\theta}^{\prime})\\}~{}d\boldsymbol{\theta}dy\\\
=&\textnormal{arg}\max_{\boldsymbol{\theta}^{\prime}\in\mathcal{T}}~{}-\mathbb{E}_{y\sim
Y_{n}}\left[\int_{\mathcal{T}}\sigma_{n}^{2}(\boldsymbol{\theta})|\mathcal{D}_{n}\cup\\{\boldsymbol{\theta}^{\prime},y(\boldsymbol{\theta}^{\prime})\\}~{}d\boldsymbol{\theta}\right].\end{split}$
(14)

Essentially, (14) seeks to find the point in $\mathcal{T}$ that likely leads
to the least overall uncertainty, if the corresponding training point was
chosen and the model updated.

Similarly, the batch-sequential active design, to select
$\Theta=\\{\boldsymbol{\theta}_{1},\ldots,\boldsymbol{\theta}_{q}\\}$, is
performed by choosing points as

$\boldsymbol{\theta}_{n+1:q}=\textnormal{arg}\max_{\Theta^{\prime}}~{}-\mathbb{E}_{y\sim
Y_{n}}\left[\int_{\mathcal{T}}\sigma_{n}^{2}(\boldsymbol{\theta})|\mathcal{D}_{n}\cup\\{\Theta^{\prime},y(\Theta^{\prime})\\}~{}d\Theta\right].$

The choice of $q>1$ particularly has advantages when the training data are
generated by running expensive computer simulations, which can be evaluated
synchronously in parallel. In the cases such as the present work, where we
seek to actively select training data from an existing set, batch-sequential
selection results in fewer hyperparameter training steps, which in the case of
exact GPs scales as $\mathcal{O}(n^{3})$. However, in terms of the improvement
in model fit, it is not obvious what choice of $q$ is the best. Therefore, we
perform a simple sensitivity study to investigate the effect of the choice of
$q$ on the model fit.

Finally, to assess model fit, we evaluate the GP posterior mean on a hold-out
test set
$\\{\boldsymbol{\theta}_{i},f(\boldsymbol{\theta}_{i})\\},~{}i=1,\ldots,n_{\textnormal{test}}$,
and compute the log root mean squared error (RMSE), defined as

$\log(RMSE)=\log\left\\{\frac{1}{n_{\textnormal{test}}}\sum_{i=1}^{n_{\textnormal{test}}}[\mu(\boldsymbol{\theta}_{i})-f(\boldsymbol{\theta})_{i}]^{2}\right\\}^{1/2}.$
(15)

We note that the log(RMSE) is computed at each step of the active learning
process with the GP model trained with the training data selected up to the
previous step. Furthermore, the log(RMSE) is independently computed for each
of the latent space outputs using the corresponding GP model.

##  5 Results

We now demonstrate our proposed data-driven machine learning approaches toward
predicting the wake velocity field, using data generated from a scanning wind
LiDAR. We begin by first discussing the latent space reconstruction accuracy
of the data, purely from the convolutional autoencoder. The reconstruction
results provide a visual estimate of the trade-off between compression due to
the autoencoder and loss/retention of information. With the latent space
representation available from the autoencoder compression, we then proceed to
learn the mapping between the inputs (operating conditions) and the latent
space, via the machine learning models. Finally, the predicted wake fields are
"decoded" via the decoder, as shown in fig. 1.

###  5.1 Compression accuracy

We first examine the ability of the convolutional autoencoder to effectively
compress the LiDAR observations to a suitable latent space. fig. 3 shows the
ability of the autoencoder to reconstruct observed data, with only four latent
dimensions, through its bottleneck neural architecture. The figures
demonstrate that despite the drastic dimensionality reduction (i.e., 2501 to
4), the reconstruction accuracy has not been compromised significantly.
However, the larger goal here is that we want to be able to generalize to a
similar reconstruction error everywhere in the parameter space. For this, we
introduce the results of the methods to obtain parameter-output maps, which we
present next.

Figure 3: The compressive effectiveness of the convolutional autoencoder on
the LiDAR data set for 3 test examples. The left column shows the true
behavior of the wake for a set of testing data (unseen during model fit) and
the right column shows reconstruction from a 4-dimensional latent space by the
autoencoder.

###  5.2 Parameteric reconstruction

The latent space $\tilde{\mathbf{x}}\in\mathbb{R}^{k}$ provides a very concise
encoding of the high-dimensional velocity field since $k\ll d$. Therefore, we
learn the function
$f:\mathbf{\theta}\in\mathbb{R}^{p}\rightarrow\tilde{\mathbf{x}}\in\mathbb{R}^{k}$.
We learn $f$ via MLP and GP regression (with exact inference). Furthermore, we
also use VGP regression (approximate inference) to improve the tractability of
GP models for large datasets. Finally, we also show the performance of
choosing training data via active learning for the exact GP.

The original LiDAR dataset has a dimensionality $d=2501$ which is reduced to
$k=4$, via the convolutional autoencoder. fig. 4 shows the actual-vs-predicted
plot of each of the four latent space dimensions predicted via all of the
three machine learning models: MLP, GP, and VGP; note that we also show
results of fitting GP with active learning used for training data selection.
Firstly, the plots show that the predictive accuracy for all the machine
learning models are somewhat similar. Secondly, there are outliers in the
dataset—as in, for example, the lower-left and upper-right corners of
subfigure $(a)$, where the prediction accuracy is poor. We attribute this
primarily to the fact that the LiDAR measurements are corrupted by noise,
whose structure is unknown and not captured by the models. Even though the GP
models do resolve the noise in the dataset, our model makes simplifying
assumptions such as independent and identically distributed noise, which might
not necessarily provide a realistic model of the noise (although they are
relatively computationally more tractable). Furthermore, the raw measurements
have missing elements, which are imputed via a local interpolation,
which—although inevitable—is expected to bias the dataset. Given these
characteristics of our real-world dataset, including more sophistication, such
as heteroscedastic noise variance, could potentially overfit the data; we
reserve those approaches for future work.

(a) Dimension 1

(b) Dimension 2

(c) Dimension 3

(d) Dimension 4

Figure 4: Fitting a parameteric map between MET data and the latent space
representation extracted by the convolutional autoencoder on the LiDAR
snapshots. The different colors indicate different latent space parameteric
mapping techniques. Similar results are obtained across different methods.

We show the predicted flow field with our machine learning models, for two
unseen parameters in fig. 5. In addition to the true flow-field obtained from
experimental observations, and the _best possible_ reconstruction—via
_autoencoding_ the true flow-field, the prediction accuracy is more-or-less
uniform across the various models, and visually presents a very close match to
the true flow field. Note that, the machine learning models at best can only
emulate the reconstruction decoded via the autoencoder and hence comparison
against the _Decoded-exact_ plot is most appropriate.

(a) Example 1

(b) Example 2

Figure 5: Two examples of parameteric reconstruction ability where information
from the MET data is used to obtain a latent-space representation for a test
data point. Following this, the decoder of the autoencoder is used to
reconstruct in physical space. We show results for the true test snapshots,
their exact reconstructions using only the autoencoder, and various
parameteric predictions in latent space followed by use of the decoder.

In fig. 6, we also show a _worst-case_ prediction via our machine learning
models. These plots correspond to the outliers identified in fig. 4. Overall,
we noticed that total number of such outliers fall roughly within 10% of the
overall dataset. For the specific wake measurement shown in fig. 6, we notice
that an anomalous speed-up is observed on the side of the wake for negative
values of the radial position. This flow feature can be either due to flow
interaction with side wind turbines or to large coherent flow structures
typically present in the atmospheric boundary layer. Even though this kind of
wake realizations are realistic, their occurrence can be relatively low and,
thus, not captured from the training dataset.

Figure 6: Potential issue with an outlier - showing limited parameteric map
expressiveness across various latent space interpolation methods. Note that
outliers were limited to around 10% of the dataset.

We have used the sequential GP model fitting via AL is another approach to
improve the tractability of (exact) GP regression for large datasets.
Essentially, we choose the most _relevant_ subset of the training data set,
instead of using the entire data set or randomly sampling from it. fig. 7
shows the log RMSE of the GP—based on a held out test data set of 1781
points—with sequential increments in the training data. Recall that the points
are chosen sequentially to reduce the average uncertainty about the GP in the
entire input domain $\mathcal{X}$, and fig. 7 shows the resulting prediction
error (log RMSE) variation with training data set size. We start the GP fit
with a randomly selected 50 points (from the full data set of 5000 points) and
sequentially add 100 more points. We also show the impact of the choice of the
number of points selected ($q$) at each step by varying it between $q=1$ to
$q=8$. We repeat this exercise 20 times for independently chosen random
starting points and plot the mean and $\pm 1$ standard deviation. Finally, we
also show the prediction error due to selecting 250 points randomly in _one
shot_ (i.e., no sequential point selection). The results show that sequential
point selection results in a smaller prediction error, regardless of the
choice of $q$, for the first 3 latent space dimensions. For the fourth latent
space dimension ($\tilde{\mathbf{x}}_{4}$), the $q=1$ still outperforms the
one-shot selection. It is worth noting that the one-shot selection of points
still has 100 points more than what was supplied to the AL GP. The reason for
the AL GP outperforming the one-shot GP is because training points are more
judiciously chosen–in this case, they are chosen specifically to minimize the
uncertainty in the GP about its own prediction. Furthermore, it can be shown
[46] that the average uncertainty in the GP is equivalent to the mean-squared
prediction error (MSPE) of the posterior mean of the GP, and therefore in
effect the AL chooses points to minimize the overall prediction error.

Figure 7: One-shot vs sequential GP fit. Horizontal lines are the RMSE with a
one-shot GP fit with 500 uniformly random points; solid and dashed lines are
the mean and $\pm 1\sigma$, respectively. Sequential fit outperforms one-shot
fit despite fewer training data.

The spatial distribution of points is visualized in fig. 8, which is a
scatterplot matrix of all possible combinations of the inputs listed in table
1. In this figure, the blue symbols indicate the full training data set (5000
points) and the red symbols are the points selected via AL. Notice that the
red points show a better spread and hence coverage of the design space
compared to the full data set. This is further emphasized by the kernel
density plots shown along the diagonal of the scatterplot matrix, where the AL
shows a larger variance which is indicative of the fact that points are more
spread out. Overall, we see that the choosing points via AL another way toward
building a more tractable GP model with large data sets, without compromising
on the predictive accuracy.

Figure 8: Selected points via active learning. Blue represents the full
training data set (5000 points) and red represents the points selected via
active learning (150 points).

fig. 9 shows a probability density plot of the prediction accuracy for all the
machine learning models we employed in this work. Note that the densities have
similar shapes indicating that the overall predictive accuracy for all the
used models is somewhat similar, as mentioned previously. It is also worth
noting that the AL GP still has the lowest accuracy amongst all the three GP
approaches presented; this can be appreciated by considering the $R^{2}$ value
in fig. 4 and/or the sorted error and kernel density plot of the error shown
in fig. 10. Our hypothesis for this behavior is that these plots show just one
realization of the GP (as opposed to fig. 7 where we show the average of 20
independent realizations). In general, the predictive accuracy of the GP is
sensitive to the choice of starting points, but fig. 7 proves that on average,
the AL GP outperforms a GP with one-shot point selection. Despite this fact,
we still show just one realization of the AL GP in fig. 4 and fig. 10 (despite
showing a relatively poor predictive accuracy) because this is more
representative of a realistic scenario.

(a) Neural network

(b) Gaussian process

(c) Variational Gaussian Process

(d) Active Learning

Figure 9: Density plots of data-driven predictions in latent space followed by
reconstruction through the decoder for (a) a fully connected neural network
(b) a Gaussian process (c) a sparse variational Gaussian process and (d) a
reconstruction using a Gaussian process trained by active learning.

(a) Sorted Error

(b) Kernel density estimate of MSE

Figure 10: The sorted root-mean-squared error on the testing data set after
parameteric latent space prediction and reconstruction in physical space
(left). Kernel density estimate of the mean-squared errors for different
parametric maps (right). The RMSEs are competitive across a large portion of
the test data set with certain outliers (approximately 200 data points).

##  6 Conclusions

Accurate predictions of wind-turbine wake flows are crucial to estimate power
efficiency of wind farms, power losses due to wake interactions, and optimal
design of wind farm layout. High-fidelity large eddy simulations, which can
resolve the turbulent flow fields with high accuracy, are prohibitively
expensive from the computational-cost standpoint to tackle the foregoing wind
energy problems. Scanning Doppler wind LiDARs provide detailed measurements of
the wake flow field generated from utility-scale wind turbines; however,
experimental data can be incomplete and noisy.

In this work, we develop data-driven machine learning models for prediction of
wind turbine wakes by leveraging wind LiDAR measurements. Specifically, we
develop several data-driven models from high-dimensional LiDAR data and
analyze their comparative performance in terms of parametric prediction. In
this regard, we first use deep neural networks to achieve a drastically
compressed ($2501\rightarrow 4$) latent-space representation of the high-
dimensional data. Then, we use multilayered perceptrons and Gaussian processes
to learn the input parameter-latent-space map. Results show that the
predictive capability of all the machine learning models are somewhat similar.
However, GP regression with exact inference resulted in the least RMSE (best
prediction). Furthermore, to address the well-known tractability issues with
exact-inference Gaussian processes, we also propose variational sparse
Gaussian processes as well as active learned Gaussian processes. These two
approaches, while resulting in a significant saving in the computation
necessary for model inference, are observed to make only a marginal trade in
the accuracy. Overall, once trained, all of the machine learning model takes
only $O(1)$ second of wall-clock time to evaluate.

This work brings together state-of-the-art machine learning, to develop data-
driven models of real-world data pertinent to wind energy. We show that robust
and accurate models can be developed despite the data being noisy and
incomplete. With the developed predictive models, sensitivity of the wake flow
field to the input parameters can be analyzed in real-time, compared to the
unrealistic cpu-hours necessary for high-fidelity simulations. In the future,
we hope to streamline the whole framework for data-driven wake modeling via
automated neural architecture search and improved uncertainty quantification
(e.g., using heteroscedastic noise models for Gaussian processes).
Furthermore, developing hybrid probabilistic machine learning models such as
deep GP models, which combine multilayered perceptrons and Gaussian processes,
is another focus of future work.

## Acknowledgement

This material is partially based upon work supported by the U.S. Department of
Energy (DOE), Office of Science, Office of Advanced Scientific Computing
Research, under Contract DE-AC02-06CH11357. This research was funded in part
and used resources of the Argonne Leadership Computing Facility, which is a
DOE Office of Science User Facility supported under Contract DE-
AC02-06CH11357. SAR acknowledges the support by Laboratory Directed Research
and Development (LDRD) funding from Argonne National Laboratory, provided by
the Director, Office of Science, of the U.S. Department of Energy under
contract DE-AC02-06CH11357. This research has been partially funded by a grant
from the National Science Foundation CBET Fluid Dynamics, award number
1705837. Pattern Energy Group is acknowledged to provide access to the wind
farm for the LiDAR experiment and wind farm data.

The submitted manuscript has been created by UChicago Argonne, LLC, Operator
of Argonne National Laboratory ("Argonne”). Argonne, a U.S. Department of
Energy Office of Science laboratory, is operated under Contract No. DE-
AC02-06CH11357. The U.S. Government retains for itself, and others acting on
its behalf, a paid-up nonexclusive, irrevocable worldwide license in said
article to reproduce, prepare derivative works, distribute copies to the
public, and perform publicly and display publicly, by or on behalf of the
Government. The Department of Energy will provide public access to these
results of federally sponsored research in accordance with the DOE Public
Access Plan (http://energy.gov/downloads/doe-public-access-plan).

## References

  * [1] Javier Sanz Rodrigo, Roberto Aurelio Chávez Arroyo, Patrick Moriarty, Matthew Churchfield, Branko Kosović, Pierre Elouan Réthoré, Kurt Schaldemose Hansen, Andrea Hahmann, Jeffrey D. Mirocha, and Daran Rife.  Mesoscale to microscale wind farm flow modeling and evaluation.  Wiley Interdisciplinary Reviews: Energy and Environment, 6(2), 2017\.
  * [2] P Veers, K Dykes, E Lantz, S Barth, CL Bottasso, O Carlson, A Clifton, J Green, P Green, H Holttinen, D Laird, V Lehtomäki, JK Lundquist, J Manwell, M Marquis, C Meneveau, P Moriarty, X Munduate, M Muskulus, J Naughton, L Pao, J Paquette, J Peinke, A Robertson, J. S Rodrigo, AM Sempreviva, JC Smith, A Tuohy, and R Wiser.  Grand challenges in the science of wind energy.  Science, 366(6464), 2019.
  * [3] Fernando Porté-agel, Majid Bastankhah, and Sina Shamsoddin.  Wind-Turbine and Wind-Farm Flows : A Review.  Boundary-Layer Meteorology, 2019.
  * [4] RJB Barthelmie, SCP Pryor, STF Frandsen, KSH Hansen, JGS Schepers, K Rados, W Schelz, A Neubert, LE Jensen, and S Neckelmann.  Quantifying the Impact of Wind Turbine Wakes on Power Output at Offshore Wind Farms.  J Atmos Ocean Technol, pages 1302–1317, 2010.
  * [5] S El-Asha, L Zhan, and GV Iungo.  Quantification of power losses due to wind turbine wake interactions through SCADA , meteorological and wind LiDAR data.  Wind Energy, 20(June):1823–1839, 2017.
  * [6] Alessandro Sebastiani, Antonio Segalini, Francesco Castellani, and Giorgio Crasto.  Data analysis and simulation of the Lillgrund wind farm.  Wind Energy, (November):1–15, 2020.
  * [7] Matthew J. Churchfield, Sang Lee, John Michalakes, and Patrick J. Moriarty.  A numerical study of the effects of atmospheric and wake turbulence on wind turbine dynamics.  Journal of Turbulence, 13:1–32, 2012.
  * [8] Davide Conti, Nikolay Dimitrov, and Alfredo Peña.  Aeroelastic load validation in wake conditions using nacelle-mounted lidar measurements.  Wind Energy Science, 5(3):1129–1154, 2020.
  * [9] Giacomo Valerio Iungo and Fernando Porté-Agel.  Volumetric lidar scanning of wind turbine wakes under convective and neutral atmospheric stability regimes.  Journal of Atmospheric and Oceanic Technology, 31(10):2035–2048, 2014.
  * [10] Stefano Letizia, Lu Zhan, and Giacomo Valerio Iungo.  LiSBOA: LiDAR Statistical Barnes Objective Analysis for optimal design of LiDAR scans and retrieval of wind statistics. Part I: Theoretical framework.  Atmospheric Measurement Techniques, 14:2065–2093, 2021.
  * [11] Lu Zhan, Stefano Letizia, and Giacomo Valerio Iungo.  Lidar measurements for an onshore wind farm: Wake variability for different incoming wind speeds and atmospheric stability regimes.  Wind Energy, 23(3):501–527, 2020.
  * [12] D Mehta, A H Van Zuijlen, B Koren, J G Holierhoek, and H Bijl.  Large Eddy Simulation of wind farm aerodynamics: A review.  J. of Wind Engineering and Industrial Aerodynamics, 133:1–17, 2014\.
  * [13] S Breton, J Sumner, J N Sørensen, K S Hansen, S Sarmast, and S Ivanell.  A survey of modelling methods for high-fidelity wind farm simulations using large eddy simulation.  Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 375(20160097), 2017.
  * [14] C. Santoni, E. J. Garcia-Cartagena, U. Ciri, L. Zhan, G. V. Iungo, and S. Leonardi.  One-way mesoscale-microscale coupling for simulating a wind farm in North Texas: Assessment against SCADA and LiDAR data.  Wind Energy, 23(3):691–710, 2020.
  * [15] B. Sanderse, S.P. van der Pijl, and B. Koren.  Review of computational fluid dynamics for wind turbine wake aerodynamics.  Wind Energy, 14:799–819, 2011.
  * [16] Giacomo Valerio Iungo, Vignesh Santhanagopalan, Umberto Ciri, Francesco Viola, Lu Zhan, Mario A Rotea, and Stefano Leonardi.  Parabolic rans solver for low-computational-cost simulations of wind turbine wakes.  Wind Energy, 21(3):184–197, 2018.
  * [17] N O Jensen.  A note on wind generator interaction.  Technical report, Risø, Roskilde, Denmark, 1983.
  * [18] Sten Frandsen, Rebecca Barthelmie, Sara Pryor, Ole Rathmann, Søren Larsen, Jørgen Højstrup, and Morten Thøgersen.  Analytical modelling of wind speed deficit in large offshore wind farms.  Wind Energy, 9(1-2):39–53, 2006.
  * [19] Majid Bastankhah and Fernando Porté-Agel.  A new analytical model for wind-turbine wakes.  Renewable Energy, 70:116–123, 2014.
  * [20] Lu Zhan, Stefano Letizia, and Giacomo Valerio Iungo.  Optimal tuning of engineering wake models through lidar measurements.  Wind Energy Science, 5(4):1601–1622, 2020.
  * [21] Ewan Machefaux, Gunner C Larsen, Tilman Koblitz, Niels Troldborg, Mark C Kelly, Abhijit Chougule, Kurt Schaldemose Hansen, and Javier Sanz Rodrigo.  An experimental and numerical study of the atmospheric stability impact on wind turbine wakes.  Wind Energy, 19(10):1785–1805, 2016.
  * [22] Romit Maulik, Vishwas Rao, S Ashwin Renganathan, Stefano Letizia, and Giacomo Valerio Iungo.  Cluster analysis of wind turbine wakes measured through a scanning doppler wind lidar.  In AIAA Scitech 2021 Forum, page 1181, 2021.
  * [23] U.S.G.S.  U.S. Geological Survey Website, 2017.
  * [24] Kai Fukami, Taichi Nakamura, and Koji Fukagata.  Convolutional neural network based hierarchical autoencoder for nonlinear mode decomposition of fluid field data.  Physics of Fluids, 32(9):095110, 2020.
  * [25] Takaaki Murata, Kai Fukami, and Koji Fukagata.  Nonlinear mode decomposition with convolutional neural networks for fluid dynamics.  Journal of Fluid Mechanics, 882, 2020.
  * [26] Bernhard Vennemann and Thomas Rösgen.  A dynamic masking technique for particle image velocimetry using convolutional autoencoders.  Experiments in Fluids, 61(7):1–11, 2020.
  * [27] Francisco J Gonzalez and Maciej Balajewicz.  Deep convolutional recurrent autoencoders for learning low-dimensional feature dynamics of fluid systems.  arXiv preprint arXiv:1808.01346, 2018.
  * [28] Romit Maulik, Bethany Lusch, and Prasanna Balaprakash.  Reduced-order modeling of advection-dominated systems with recurrent neural networks and convolutional autoencoders.  Physics of Fluids, 33(3):037106, 2021.
  * [29] Youngkyu Kim, Youngsoo Choi, David Widemann, and Tarek Zohdi.  A fast and accurate physics-informed neural network reduced order model with shallow masked autoencoder.  arXiv preprint arXiv:2009.11990, 2020.
  * [30] Pin Wu, Siquan Gong, Kaikai Pan, Feng Qiu, Weibing Feng, and Christopher Pain.  Reduced order model using convolutional auto-encoder with self-attention.  Physics of Fluids, 33(7):077107, 2021.
  * [31] M Cheng, F Fang, CC Pain, and IM Navon.  An advanced hybrid deep adversarial autoencoder for parameterized nonlinear fluid flow modelling.  Computer Methods in Applied Mechanics and Engineering, 372:113375, 2020.
  * [32] S Ashwin Renganathan, Yingjie Liu, and Dimitri N Mavris.  Koopman-based approach to nonintrusive projection-based reduced-order modeling with black-box high-fidelity models.  AIAA Journal, 56(10):4087–4111, 2018.
  * [33] S Ashwin Renganathan.  Koopman-based approach to nonintrusive reduced order modeling: Application to aerodynamic shape optimization and uncertainty propagation.  AIAA Journal, 58(5):2221–2235, 2020.
  * [34] S Ashwin Renganathan, Romit Maulik, and Vishwas Rao.  Machine learning for nonintrusive model order reduction of the parametric inviscid transonic flow past an airfoil.  Physics of Fluids, 32(4):047110, 2020.
  * [35] George Cybenko.  Approximation by superpositions of a sigmoidal function.  Mathematics of control, signals and systems, 2(4):303–314, 1989\.
  * [36] Andrew R Barron.  Universal approximation bounds for superpositions of a sigmoidal function.  IEEE Transactions on Information theory, 39(3):930–945, 1993.
  * [37] S Ashwin Renganathan, Romit Maulik, and Jai Ahuja.  Enhanced data efficiency using deep neural networks and gaussian processes for aerodynamic design optimization.  Aerospace Science and Technology, 111:106522, 2021.
  * [38] S Ashwin Renganathan, Jeffrey Larson, and Stefan M Wild.  Lookahead acquisition functions for finite-horizon time-dependent bayesian optimization and application to quantum optimal control.  arXiv preprint arXiv:2105.09824, 2021.
  * [39] Dushhyanth Rajaram, Tejas G Puranik, Ashwin Renganathan, Woong Je Sung, Olivia J Pinon-Fischer, Dimitri N Mavris, and Arun Ramamurthy.  Deep gaussian process enabled surrogate models for aerodynamic flows.  In AIAA Scitech 2020 Forum, page 1640, 2020.
  * [40] Dushhyanth Rajaram, Tejas G Puranik, S Ashwin Renganathan, WoongJe Sung, Olivia Pinon Fischer, Dimitri N Mavris, and Arun Ramamurthy.  Empirical assessment of deep gaussian process surrogate models for engineering problems.  Journal of Aircraft, 58(1):182–196, 2021.
  * [41] Andrew Gelman, John B Carlin, Hal S Stern, and Donald B Rubin.  Bayesian data analysis.  Chapman and Hall/CRC, 1995.
  * [42] Christopher K Williams and Carl Edward Rasmussen.  Gaussian processes for machine learning, volume 2.  MIT press Cambridge, MA, 2006.
  * [43] Michael L Stein.  Interpolation of spatial data: some theory for kriging.  Springer Science & Business Media, 2012.
  * [44] Noel Cressie and Hsin-Cheng Huang.  Classes of nonseparable, spatio-temporal stationary covariance functions.  Journal of the American Statistical Association, 94(448):1330–1339, 1999.
  * [45] Bertil Matérn.  Spatial variation, volume 36.  Springer Science & Business Media, 2013.
  * [46] Thomas J Santner, Brian J Williams, William I Notz, and Brain J Williams.  The design and analysis of computer experiments, volume 1.  Springer, 2003.
  * [47] Edward Snelson and Zoubin Ghahramani.  Sparse gaussian processes using pseudo-inputs.  Advances in neural information processing systems, 18:1257–1264, 2005.
  * [48] David M Blei, Alp Kucukelbir, and Jon D McAuliffe.  Variational inference: A review for statisticians.  Journal of the American statistical Association, 112(518):859–877, 2017.
  * [49] James Hensman, Alexander Matthews, and Zoubin Ghahramani.  Scalable variational gaussian process classification.  In Artificial Intelligence and Statistics, pages 351–360. PMLR, 2015.
  * [50] David A Cohn, Zoubin Ghahramani, and Michael I Jordan.  Active learning with statistical models.  Technical report, MASSACHUSETTS INST OF TECH CAMBRIDGE ARTIFICIAL INTELLIGENCE LAB, 1995.
  * [51] David A Cohn, Zoubin Ghahramani, and Michael I Jordan.  Active learning with statistical models.  Journal of artificial intelligence research, 4:129–145, 1996.