# Hyperparameter Search using Genetic Algorithm for Surrogate Modeling of
Geophysical Flows

Suraj Pawar 111Ph.D. Candidate, School of Mechanical & Aerospace Engineering,
Oklahoma State University, Stillwater, OK 74078, USA. School of Mechanical &
Aerospace Engineering, Oklahoma State University, Stillwater, OK 74078, USA.
Omer San 222Associate Professor, School of Mechanical & Aerospace Engineering,
Oklahoma State University, Stillwater, OK 74078, USA. School of Mechanical &
Aerospace Engineering, Oklahoma State University, Stillwater, OK 74078, USA.
Gary G. Yen 333Regents Professor, School of Electrical & Computer Engineering,
Oklahoma State University, Stillwater, OK 74078, USA. School of Electrical &
Computer Engineering, Oklahoma State University, Stillwater, OK 74078, USA.

###### Abstract

The computational models for geophysical flows are computationally very
expensive to employ in multi-query tasks such as data assimilation,
uncertainty quantification, and hence surrogate models sought to alleviate the
computational burden associated with these full order models. Researchers have
started applying machine learning algorithms, particularly neural networks, to
build data-driven surrogate models for geophysical flows. The performance of
the neural network highly relies upon its architecture design and selection of
other hyperparameters. These neural networks are usually manually designed
through trial and error to maximize their performance. This often requires
domain knowledge of the neural network as well as the problems of interest.
This limitation can be addressed by using an evolutionary algorithm to
automatically design architecture and select optimal hyperparameters of the
neural network. In this paper, we apply the genetic algorithm to effectively
design the long short-term memory (LSTM) neural network to build the non-
intrusive reduced order model of the sea-surface temperature field.

##  1 Introduction

There are numerous tasks in engineering design, control, and climate modeling
that require solving a partial differential equation (PDE) based forward model
for many ensembles. For example, in sequential data assimilation, an ensemble
of models is run to approximate the covariance matrix [1]. As the desired
level of accuracy increases, the resolution of the spatiotemporal numerical
discretization also increases. This can cause a significant increase in the
computational resources and can become the bottleneck in the outer loop of the
design and forecast cycle. For example in computational fluid dynamics (CFD),
the three-dimensional simulation of turbulent flows [2] is often used to
obtain an exact solution to the Navier-Stokes equation but is less frequently
used in related tasks such as shape optimization [3]. Therefore, there is a
sustained interest in developing a reduced-order model (ROM) that is
computationally much cheaper than the full order model (FOM) and also provides
the solution with a sufficient level of accuracy [4, 5, 6].

In recent years, there is a growing interest in developing non-intrusive ROMs
for physical systems [7, 8, 9]. The non-intrusive refers to the use of only
data to construct the ROM. The non-intrusive ROMs are particularly attractive
for systems where the perfect model is unknown, but there is an abundance of
data available for that system. This situation is commonplace for geophysical
flows. The geophysical system model can be imperfect due to a variety of
reasons such as approximation of subgrid-scale processes due to insufficient
grid resolution, uncertainty in model parameters or incorrect structure of the
model itself [10]. However, there is an abundance of data from remote sensing,
satellite observations, and experimental measurements for geophysical flows.
As a result, there have been several works to harness this data for the
effective forecast of geophysical processes [11, 12].

In this work, we focus on the use of linear dimensionality reduction technique
followed by the machine learning algorithm to evolve the latent space [13, 14,
15]. Specifically, we use proper orthogonal decomposition (POD) to identify
the latent space of the FOM, and then the long short-term memory (LSTM) neural
network is utilized to model the evolution of latent space. One of the main
advantages of this surrogate modeling technique is that it is purely data-
driven and hence is particularly appealing to geophysical flows. The archival
data collected from remote sensing and in-situ measurements can be exploited
to build the surrogate model and then the surrogate model can be employed for
the forecast task or in the data assimilation cycle. As the new data becomes
available, the surrogate model can be retrained using the transfer learning
approach [16] to improve the forecast performance.

Apart from ROM, the use of the deep neural network has significantly increased
for many scientific applications [17, 18, 19]. One of the main challenges with
neural networks is that their performance is highly dependent on their
architecture [20]. Additionally, the neural network has a large number of
hyperparameters that are problem-dependent. Usually, the neural networks are
manually designed through trial and error and this procedure can be time-
consuming. Moreover, rich domain knowledge about the data and neural networks
is required to achieve a promising performance. There are methods like grid
search [21] or a random search [22] to find a good combination of
hyperparameters. But these methods are not scalable as the dimension of the
search space goes up. Hence, there is a growing interest in automating the
design of deep neural network architecture and hyperparameter search that will
allow the user with no specialized knowledge of the neural network to apply to
their problem of interest [23]. In the present study, we use the genetic
algorithm (GA) to optimize the architecture design and hyperparameters of the
LSTM neural network. Readers are referred to a recent survey by Liu et al.
[24] on most recent evolutionary computation-based neural architecture search
methods in light of the core components, to systematically discuss their
design principles and justifications on the design. The LSTM network is used
within the surrogate model design of the NOAA sear-surface temperature (SST)
data set.

The rest of the paper is structured as follows. Firstly, we describe the non-
intrusive ROM methodology and data pre-processing in Section 2. Then, the
details of the proposed algorithm are provided in Section 3. Section 4
discusses the experimental results and their analysis. Finally, in Section 5,
conclusions and future works are outlined.

##  2 Non-intrusive Reduced Order Modeling

###  2.1 Proper orthogonal decomposition

We use the proper orthogonal decomposition (POD) to extract the dominant modes
from the data. We collect the data snapshots
$\mathbf{u}_{1},\mathbf{u}_{2},\dots,\mathbf{u}_{N}$ $\in$ $\mathbb{R}^{M}$ at
different time instances. Here, $M$ is the spatial degree of freedom which is
equal to the total number of grid points, and $N$ is the number of snapshots.
In POD, we construct a set of orthonormal basis functions that optimally
describes the field variable of the system in $L_{2}$-norm. The snapshot data
matrix $\mathbf{A}$ is formed as

$\mathbf{A}=[\mathbf{\tilde{u}}_{1}|\mathbf{\tilde{u}}_{2}|\dots|\mathbf{\tilde{u}}_{N}]\in\mathbb{R}^{M\times
N},$ (1)

where the mean-subtracted (or anomaly) fields $\mathbf{\tilde{u}}_{i}\in$
$\mathbb{R}^{M}$ are computed as follows

$\mathbf{\tilde{u}}_{i}=\mathbf{u}_{i}-\bar{\mathbf{u}},\quad\bar{\mathbf{u}}=\frac{1}{N}\sum_{i=1}^{N}\mathbf{u}_{i},$
(2)

where $\bar{\mathbf{u}}$ is the time-average of the solution field. Once the
snapshot data matrix is formed, we use singular value decomposition (SVD) to
compute the left and right singular vectors of the matrix $\mathbf{A}$. In
matrix form, the SVD can be written as

$\mathbf{A}=\mathbf{W}\mathbf{\Sigma}\mathbf{V}^{T}=\sum_{k=1}^{N}\sigma_{k}\mathbf{w}_{k}\mathbf{v}_{k}^{T},$
(3)

where $\mathbf{W}\in\mathbb{R}^{M\times N}$,
$\mathbf{\Sigma}\in\mathbb{R}^{N\times N}$, and
$\mathbf{V}\in\mathbb{R}^{N\times N}$. The $\mathbf{W}$ and $\mathbf{V}$
contains the left and right singular vectors which are identical to the
eigenvectors of $\mathbf{A}\mathbf{A}^{T}$ and $\mathbf{A}^{T}\mathbf{A}$,
respectively. Also, the square of singular vales are equal to the eigenvalues,
i.e., $\lambda_{k}=\sigma_{k}^{2}$. The vectors $\mathbf{w_{k}}$ (also the
eigenvectors of $\mathbf{A}\mathbf{A}^{T}$) are the POD basis functions and we
denote them as $\phi_{k}$ in this text. The POD basis functions are
orthonormal (i.e., $\langle\phi_{i},\phi_{j}\rangle=\delta_{ij}$, where
$\delta_{ij}$ is the Kronecker delta) and are computed in an optimal manner in
the $L_{2}$ sense [25, 26]. The state of the dynamical system can be
approximated using these POD basis functions as follows

$\mathbf{u}(\mathbf{x},t)=\bar{\mathbf{u}}+\sum_{k=1}^{R}a_{k}(t)\phi_{k}(\mathbf{x}),$
(4)

where $R$ is the number of retained basis functions such that $R<<N$ and
$a_{k}$ are the time-dependent modal coefficients. The POD basis functions
minimize the mean-square error between the field variable and its truncated
representation. The number of retained modes is usually decided based on their
energy content. Using these retained modes, we can form the POD basis set
$\mathbf{\Phi}=\\{\phi_{k}\\}_{k=1}^{R}$ to build the ROM.

###  2.2 Model-free evolution of the ROM

The model-free prediction of the evolution of ROMs has a promising
application, especially for geophysical dynamical systems. For a variety of
tasks such as data assimilation, uncertainty quantification of geophysical
flows, the ensemble of the forward model needs to be run and this can be
computationally expensive. Instead of evolving the forward model, the
surrogate model based on ROMs can be utilized. Also for many geophysical
dynamical systems, the mechanistic description of the dynamical system is
unavailable or incomplete due to coarse grid resolution. However, there is an
abundance of data obtained from local and satellite observations for many past
decades. In recent years, data-driven approaches exploiting the recurrent
neural network (RNN) are proven successful for the model-free prediction of
spatiotemporal chaotic dynamical systems [27, 28].

The time-dependent modal-coefficients given in Equation 4 are obtained by
projecting the mean subtracted field on the POD basis as follows

$a_{k}(t_{i})=\langle\mathbf{\tilde{u}}_{i};\phi_{k}\rangle,$ (5)

where angle parenthesis denotes the inner product of two functions. Therefore,
if we can model the prediction of the evolution of modal coefficients we can
reconstruct the solution field in the future. One of the popular approaches of
time-series prediction is the RNN and has been recently applied to several
studies on non-intrusive ROMs for physical systems [13, 29, 30, 31]. This is
motivated by the fact that RNN takes into account the past history of the
system for future state prediction and the availability of historical data to
train the neural network. The long short-term memory (LSTM) neural network is
one of the most successful variants of the RNNs that solves the vanishing
gradient problem [32]. In this study, we consider an input sequential data
matrix $\mathcal{X}$ and the output sequential data matrix $\mathcal{Y}$. Each
sample of the input training matrix $\mathcal{X}$, i.e., $\mathcal{X}_{n}$ is
constructed as $\left\\{a_{1}^{(k)},\dots,a_{R}^{(k)};\ \dots\
;a_{1}^{(k-\sigma+1)},\dots,a_{R}^{(k-\sigma+1)}\right\\}$ and the
corresponding output sample in output sequential data matrix $\mathcal{Y}$,
i.e., $\mathcal{Y}_{n}$ is
$\left\\{a_{1}^{(k+1)},\dots,a_{R}^{(k+1)}\right\\}$, where $k$ corresponds to
the time index, and $n$ is for the $n$th sample. The parameter $\sigma$ is
called the lookback time-windows. Since we are using the information of only
$\sigma$ past temporally consecutive states as the input, the LSTM can capture
dependencies up to $\sigma$ previous time steps.

###  2.3 Data preparation

Model-free prediction is particularly important for real-world data which is a
manifestation of many multiscale processes that cannot be approximated
accurately with any of the model-based approaches. The dynamics of the flow in
real-life is governed by the coupling with external systems, and this coupling
might be unknown or cannot be modeled. Therefore, the real-world data is not
well-behaved and non-intrusive methods are well suited for these complex
flows. To this end, we investigate the application of non-intrusive ROMs for
the NOAA Optimum Interpolation Sea Surface Temperature (SST) V2 data
set444https://psl.noaa.gov/ (Accessed 2020). This data set consists of the
weekly average sea surface temperature (SST) data snapshots on a one-degree
resolution grid and is produced by combining satellite and local measurements.
The seasonal fluctuations give rise to a strong periodic structure to the
temperature field in this data set.

Each snapshot of the data set has the dimension of $180\times 360$ along the
latitudes and longitudes with one-degree grid resolutions. The data points
corresponding to land are removed by using the mask operation and the
surrogate model is built for the flattened data corresponding to only the
ocean surface. This data is available from October 22, 1981, to June 30, 2018
(i.e., total 1,914 snapshots).

Once the data prepossessing is done, then the snapshot data matrix
$\mathbf{A}$ is formed using the mean-subtracted anomalies of the temperature
field. The number of modes for the ROM, R, is selected based on a relative
information content (RIC) formula as given below

$\text{RIC}(R)=\left(\frac{\sum_{k=1}^{R}\sigma_{k}^{2}}{\sum_{k=1}^{N}\sigma_{k}^{2}}\right).$
(6)

The RIC represents the fraction of information (variance) of the total data
that can be recovered using $R$ basis functions. Figure 1 shows the RIC
percentage for the NOAA SST data set. We fix the number of retained modes to
be $R=8$ that captures around 92% (i.e., $9.182\times 10^{1}$) of the variance
of the data and these modes are sufficient to capture the seasonal trends in
the SST data set. It can also be seen that the increase in the number of modes
after 8 modes gives a very small increase in the total information variance,
as these modes are mainly responsible for capturing small-scale fluctuations.

Figure 1: Percentage of the square of singular values of the snapshot data
matrix $\mathbf{A}$ (equivalent to eigenvalues of $\mathbf{A}\mathbf{A}^{T}$
or $\mathbf{A}^{T}\mathbf{A}$). The number of retained modes for the ROM is
$R=8$.

##  3 Genetic Algorithm for LSTM Network

The performance of the neural network is highly dependent on the design of the
neural network architecture and the selection of other hyperparameters such as
activation function, optimizer, weight initialization, etc. Similar to other
hyperparameter optimization studies [33, 34], we encode the entire LSTM
architecture using different types of small architectures. Figure 2 shows
different types of small architectures that act as the building blocks of the
main LSTM network. In each of these small architectures, skip-connection is
employed as it allows training a deep neural network without the vanishing
gradient problem [35, 36]. The final LSTM network is constructed by ordering
these building blocks. The number of LSTM cells and the number of building
blocks are the two additional parameters related to the design of the neural
network architectures that need to be optimized using the genetic algorithm.
The other hyperparameters to optimize are the type of the optimizer (for
example, SGD and Adam), the learning rate of the optimizer, type of
initialization for weights and biases such as RandomUniform and Xavier
initialization, and the type of activation function (for example, ReLU and
tanh).

Figure 2: Different types of encoded small architectures that act as the
building block of the entire LSTM architecture.

Algorithm 1 lists the different steps of the proposed framework. We start by
initializing the population in a given size. The LSTM network is randomly
constructed for each individual from the predefined building blocks and other
hyperparameters are assigned based on the assignment. The neural network is
evaluated using the three-fold cross-validation on the training dataset. Once
the fitness (i.e., the validation mean squared error in this study) is
evaluated for each individual in the population, they are sorted based on
their fitness values. From this population, a certain percentage of elite
individuals (i.e., the individuals with better fitness values) are retained
and some individuals are selected using the tournament selection for offspring
generation. The offspring are generated from the selected parents using
crossover and mutation operators. The new population is formed by combining
the retained elite individuals from the old population and the newly generated
offspring. The counter for the number of generations is increased by one and
this procedure is continued for the specified number of generations.

Algorithm 1 The proposed genetic algorithm

0: A set of predefined building blocks, population size, maximum number of
generations, dataset for training, elite percentage to retain.

0: Discovered LSTM architecture along with other hyperparameters

1: $P_{0}\leftarrow$ Initialize the population with the given population size.

2: $k\leftarrow 0$

3: for $k=0$ to $the~{}maximum~{}generation~{}number$ do

4: Evaluate the fitness of each individual within the population $P_{k}$;

5: $P_{k}^{E}\leftarrow$ Retain the elite population with the maximum value of
the fitness;

6: $P_{k}^{T}\leftarrow$ Select the population for offspring generation using
tournament selection;

7: $Q_{k}\leftarrow$ Generate offspring from the selected population using the
proposed crossover and mutation operations

8: $P_{k+1}\leftarrow$ $P_{k}^{E}\cup Q_{k}$

9: $k=k+1$

10: end for

11: return the individual with the best fitness form $P_{k}$

Each individual in the population is defined using 6 integer numbers and 1
real number. These numbers represent the hyperparameter of the LSTM network.
The upper and lower limit of these hyperparameters and the type of
distribution used for generating these numbers are reported in Table 1. The
first parameter is the type of the building block of the network and four
different building blocks as shown in Figure 2 that can be used to design the
LSTM network. The second and third parameters are related to how deep the LSTM
network will be and how wide the neural network within the LSTM cell will be.
The other parameters such as activation function, optimizer, learning rate,
and weight initialization are related to the training of the LSTM network.
During the evaluation, the LSTM network corresponding to each individual is
trained and the validation mean squared error (MSE) is assigned as the fitness
for that individual. The population is sorted based on the value of the
fitness in ascending order (lower MSE means better fitness).

Table 1: Hyperaprameters of the LSTM network that are optimized Hyperparameter |  Lower bound |  Upper bound |  Type |  Available options
---|---|---|---|---
Type of building block |  1 |  4 |  Uniform |  Shown in Figure 2
Number of hidden layers |  2 |  8 |  Uniform |  2,3,4,5,6,7,8
Number of LSTM cell in each block |  40 |  160 |  Normal |  Integer between 40 and 160
Optimizer |  1 |  3 |  Uniform |  Adam, RMSprop, SGD
Activation function |  1 |  3 |  Uniform |  ReLU, tanh, LeakyReLU
Weight initialization |  1 |  3 |  Uniform |  RandomNormal, RandomUniform, GlorotNormal
Learning rate of the optimizer |  $10^{-4}$ |  $10^{-2}$ |  Lognormal |  Real number between $10^{-4}$ to $10^{-2}$

The next step is the selection of parents for the offspring generation. We use
tournament selection to select parents. Algorithm 2 shows the procedure of
tournament selection. In tournament selection, a few individuals are selected
randomly from the population and the best individual is selected for offspring
generation.

Algorithm 2 Parents selection

0: The population $P_{k}$, population size $N$, number of elite individuals
$M$, tournament size $T$

0: Population for offspring generations $P_{k}^{T}$

1: $P_{k}^{T}\leftarrow\emptyset$

2: for $i=0$ to $N-M$ do

3: $\mathbf{T}\leftarrow$ Randomly generate $T$ numbers without repetition;

4: $P_{k}^{T}\leftarrow P_{k}^{T}\cup P_{k}[\text{min}(\mathbf{T})]$;

5: end for

6: return The selected population $P_{k}^{T}$

The details of generating the offspring are provided in Algorithm 3. The
offspring population is generated in two stages. The first stage is the
crossover operation (lines 2 -13), and we utilize the uniform crossover
operator. In the uniform crossover operation, two individuals are selected
from the population. Then for every gene of the individual, a random number is
generated, and if this number is greater than the crossover probability then
the corresponding gene at the position $j$ of two individuals is swapped.
During the mutation, the gene at the position $i$ is changed if the random
number is greater than the mutation probability.

Algorithm 3 Offspring generation

0: The selected population $P_{k}^{T}$ containing individual and their
fitness, the probability of crossover operation $p_{c}$, mutation operation
probability $p_{m}$.

0: The offspring population $Q_{k}$

1: $Q_{k}\leftarrow\emptyset$

2: while $i<|P_{k}^{T}|$ do

3: $p_{1},p_{2}\leftarrow$ Select two consecutive individuals from $P_{k}$
(i.e, $2i,2i+1$)

4: for $j$ = $0$ to $|p_{1}|$ do

5: $r\leftarrow$ Randomly generate a number between $[0,1]$;

6: if $r>p_{c}$ then

7: Swap the genes of $p_{1}$ and $p_{2}$

8: end if

9: $j\leftarrow j+1$

10: end for

11: $Q_{k}\leftarrow Q_{k}\cup p_{1}\cup p_{2}$;

12: $i\leftarrow i+1$

13: end while

14: for each individual $p$ in $Q_{k}$ do

15: $r\leftarrow$ Randomly generate a number between $[0,1]$;

16: if $r<p_{m}$ then

17: $i\leftarrow$ Randomly choose a point in $p$;

18: Change the gene at the point $i$ of $p$;

19: end if

20: end for

21: return Offspring population $Q_{k}$

##  4 Numerical Experiments

In this Section, we provide the result of the GA in optimizing the
architecture design and hyperparameters of the neural network. Then, we show
the performance of the optimized LSTM network in emulating and forecasting the
SST. The NOAA SST data set is available from October 22, 1981, to June 30,
2018, which corresponds to 1,914 snapshots. We utilize 70% of the data
selected randomly from the first 1,500 snapshots as the training data set.
During the hyperparameter optimization of the LSTM network, the individual
within the population is trained for 100 epochs and three-fold cross-
validation is used to avoid overfitting. The lookback time window of the LSTM
network is set at $\sigma=8$, which is found to capture the temporal
correlation in modal coefficients. The fitness of each individual is the mean
of the mean squared error (MSE) for the held-out data set from three folds.
The batch size is fixed at 64 for the training. The population size for this
experiment over the number of generations is set at 20. Figure 3 displays the
evolution trajectory of the fitness of individuals within the population with
generations of the GA. The boxplot is utilized to show the statistics of the
population at each generation. Figure 3 also reports the median and minimum
MSE of the population at each generation. As the evolution progresses, the MSE
on the validation data set decreases. The height of the box indicates the
variance of the fitness of the population at each generation, and we can
observe that the variance is substantially low after the first generation. The
MSE decreases sharply from the first generation to the second generation,
which can be attributed to the random initialization of the population at the
beginning of the GA. For the problem investigated in this study, the five to
ten generations seem enough to find the best set of hyperparameters.

Figure 3: The evolutionary trajectory of the proposed algorithm in discovering
the best architecture of the LSTM on the NOAA SST data set.

The optimal neural network architecture discovered by the genetic algorithm is
shown in Figure 4. Once the best architecture and other hyperparameters are
discovered using the GA, the best LSTM network is trained for 1200 epochs with
a batch size of 64. During the deployment of the trained network, the auto-
regressive method is used to forecast the modal coefficients. The initial
condition for modal coefficients for the first 8 time steps (equal to the
lookback of the LSTM) is provided. This information is used to predict the
modal coefficients at the 9th time step. Then the modal coefficients from the
2nd time step to the 9th time step are used to predict the modal coefficients
at the 10th time step. This procedure is repeated up to the final time step,
i.e., the 1,914th-time step. Since we use the auto-regressive deployment of
the trained LSTM network, only the initial condition corresponding to the LSTM
lookback is required. After a few time steps, only the prediction of the LSTM
is used to forecast the modal coefficients in the future. In Figure 5, the
true modal coefficients and the predicted modal coefficients are shown. We can
see that there is a very good agreement between the true and predicted modal
coefficients, especially for the first few modes which are responsible for
capturing large-scale fluctuations and seasonal patterns.

Figure 4: Optimal neural network architecture discovered using the genetic
algorithm. The selected hyperparameters are RandomNormal initialization, Adam
optimizer with the learning rate = 0.0047.  Figure 5: Time series prediction
of the modal coefficients with the best LSTM network discovered by the
proposed GA algorithm.

Figure 6 depicts the true and reconstructed temperature field at two different
times. The temperature field is reconstructed using Equation 4, where the
average temperature field is computed using the first 1,500 snapshots. The
large patterns in the temperature field are accurately captured with the LSTM
based ROM which demonstrates the potential of data-driven ROMs for geophysical
flows. In Figure 7, the $L_{2}$-norm of the difference between the true and
predicted temperature field is shown. We can observe that the error grows
toward the final time, which could be due to the inaccurate prediction of
modal coefficients and the inability of the POD basis functions to capture the
spatial pattern with high accuracy in the extrapolated time region (i.e.,
beyond 1500). One of the remedies to this behavior with the LSTM error piling
up is to train and deploy the network in a non-auto-regressive manner [37].

Figure 6: Sample averaged temperature field in degrees Celsius at time $T=750$
(left) and $T=1500$ (right). Top row represent the true temperature field and
the bottom row represent the reconstructed temperature field using the ROM.
Figure 7: Temporal variation of the $L_{2}$-norm of the difference between the
true and predicted temperature field.

##  5 Concluding Remarks

We presented a genetic algorithm (GA) for automatic architecture search and
hyperparameter optimization of the long short-term memory (LSTM) neural
network for the task of surrogate modeling of geophysical flows. We have
successfully demonstrated the surrogate model built using the optimized LSTM
is able to forecast the sea-surface temperature (SST) field. To facilitate the
training of deeper neural networks, we adopt an encoding strategy in which the
LSTM network is designed using small building blocks that incorporate skip
connections. The optimized LSTM network can predict the seasonal variation in
the SST field accurately without any instability issue for a sufficiently long
period.

We observe that the discrepancy between the true and predicted modal
coefficients is higher for the forecast period compared to the training zone.
This issue can be addressed either using the non-auto-regressive deployment or
transfer learning approach where the LSTM is retrained as the new data becomes
available. In the present study, we assume that the number of LSTM cells in
each hidden layer is constant. This constraint can be removed and the LSTM
network can be optimized with the variable-length encoding strategy for each
individual in the population. The variable-length encoding strategy will let
us design the neural network architecture that employs a different number of
hidden units in each of the LSTM cells and we plan to investigate this in our
future studies. Another challenge with the neural architecture search is the
computational resources required to train the neural network, and this
challenge can be addressed using asynchronous computation.

## Acknowledgments

This material is based upon work supported by the U.S. Department of Energy,
Office of Science, Office of Advanced Scientific Computing Research under
Award Number DE-SC0019290. O.S. gratefully acknowledges their support.
Disclaimer: This report was prepared as an account of work sponsored by an
agency of the United States Government. Neither the United States Government
nor any agency thereof, nor any of their employees, makes any warranty,
express or implied, or assumes any legal liability or responsibility for the
accuracy, completeness, or usefulness of any information, apparatus, product,
or process disclosed, or represents that its use would not infringe privately
owned rights. Reference herein to any specific commercial product, process, or
service by trade name, trademark, manufacturer, or otherwise does not
necessarily constitute or imply its endorsement, recommendation, or favoring
by the United States Government or any agency thereof. The views and opinions
of authors expressed herein do not necessarily state or reflect those of the
United States Government or any agency thereof.

## References

  * Lewis et al. [2006] Lewis, J. M., Lakshmivarahan, S., and Dhall, S., _Dynamic data assimilation: a least squares approach_ , Vol. 104, Cambridge University Press, Cambridge, 2006.
  * Moin and Mahesh [1998] Moin, P., and Mahesh, K., “Direct numerical simulation: a tool in turbulence research,” _Annual Review of Fluid Mechanics_ , Vol. 30, No. 1, 1998, pp. 539–578.
  * Polat and Tuncer [2013] Polat, O., and Tuncer, I. H., “Aerodynamic shape optimization of wind turbine blades using a parallel genetic algorithm,” _Procedia Engineering_ , Vol. 61, 2013, pp. 28–31.
  * Lucia et al. [2004] Lucia, D. J., Beran, P. S., and Silva, W. A., “Reduced-order modeling: new approaches for computational physics,” _Progress in Aerospace Sciences_ , Vol. 40, No. 1-2, 2004, pp. 51–117.
  * Taira et al. [2017] Taira, K., Brunton, S. L., Dawson, S. T., Rowley, C. W., Colonius, T., McKeon, B. J., Schmidt, O. T., Gordeyev, S., Theofilis, V., and Ukeiley, L. S., “Modal analysis of fluid flows: An overview,” _AIAA Journal_ , Vol. 55, No. 12, 2017, pp. 4013–4041.
  * Benner et al. [2015] Benner, P., Gugercin, S., and Willcox, K., “A survey of projection-based model reduction methods for parametric dynamical systems,” _SIAM Review_ , Vol. 57, No. 4, 2015, pp. 483–531.
  * Xiao et al. [2015] Xiao, D., Fang, F., Pain, C., and Hu, G., “Non-intrusive reduced-order modelling of the Navier–Stokes equations based on RBF interpolation,” _International Journal for Numerical Methods in Fluids_ , Vol. 79, No. 11, 2015, pp. 580–595.
  * Hesthaven and Ubbiali [2018] Hesthaven, J. S., and Ubbiali, S., “Non-intrusive reduced order modeling of nonlinear problems using neural networks,” _Journal of Computational Physics_ , Vol. 363, 2018, pp. 55–78.
  * Pawar et al. [2019] Pawar, S., Rahman, S., Vaddireddy, H., San, O., Rasheed, A., and Vedula, P., “A deep learning enabler for nonintrusive reduced order modeling of fluid flows,” _Physics of Fluids_ , Vol. 31, No. 8, 2019, p. 085101.
  * Schneider et al. [2017] Schneider, T., Lan, S., Stuart, A., and Teixeira, J., “Earth system modeling 2.0: A blueprint for models that learn from observations and targeted high-resolution simulations,” _Geophysical Research Letters_ , Vol. 44, No. 24, 2017, pp. 12–396.
  * O’Gorman and Dwyer [2018] O’Gorman, P. A., and Dwyer, J. G., “Using machine learning to parameterize moist convection: Potential for modeling of climate, climate change, and extreme events,” _Journal of Advances in Modeling Earth Systems_ , Vol. 10, No. 10, 2018, pp. 2548–2563.
  * Reichstein et al. [2019] Reichstein, M., Camps-Valls, G., Stevens, B., Jung, M., Denzler, J., and Carvalhais, N., “Deep learning and process understanding for data-driven Earth system science,” _Nature_ , Vol. 566, No. 7743, 2019, pp. 195–204.
  * Mohan and Gaitonde [2018] Mohan, A. T., and Gaitonde, D. V., “A deep learning based approach to reduced order modeling for turbulent flow control using LSTM neural networks,” _arXiv preprint arXiv:1804.09269_ , 2018.
  * Rahman et al. [2019] Rahman, S. M., Pawar, S., San, O., Rasheed, A., and Iliescu, T., “Nonintrusive reduced order modeling framework for quasigeostrophic turbulence,” _Physical Review E_ , Vol. 100, No. 5, 2019, p. 053306.
  * Pawar et al. [2020] Pawar, S., Ahmed, S. E., San, O., and Rasheed, A., “Data-driven recovery of hidden physics in reduced order modeling of fluid flows,” _Physics of Fluids_ , Vol. 32, No. 3, 2020, p. 036602.
  * Tan et al. [2018] Tan, C., Sun, F., Kong, T., Zhang, W., Yang, C., and Liu, C., “A survey on deep transfer learning,” _International conference on artificial neural networks_ , Springer, 2018, pp. 270–279.
  * Mjolsness and DeCoste [2001] Mjolsness, E., and DeCoste, D., “Machine learning for science: state of the art and future prospects,” _Science_ , Vol. 293, No. 5537, 2001, pp. 2051–2055.
  * Brunton et al. [2020] Brunton, S. L., Noack, B. R., and Koumoutsakos, P., “Machine learning for fluid mechanics,” _Annual Review of Fluid Mechanics_ , Vol. 52, 2020, pp. 477–508.
  * Pawar and San [2021] Pawar, S., and San, O., “Data assimilation empowered neural network parametrizations for subgrid processes in geophysical flows,” _Physical Review Fluids_ , Vol. 6, No. 5, 2021, p. 050501.
  * Krizhevsky et al. [2017] Krizhevsky, A., Sutskever, I., and Hinton, G. E., “Imagenet classification with deep convolutional neural networks,” _Communications of the ACM_ , Vol. 60, No. 6, 2017, pp. 84–90.
  * Bergstra et al. [2011] Bergstra, J., Bardenet, R., Bengio, Y., and Kégl, B., “Algorithms for hyper-parameter optimization,” _25th annual conference on neural information processing systems (NIPS 2011)_ , Vol. 24, Neural Information Processing Systems Foundation, 2011.
  * Bergstra and Bengio [2012] Bergstra, J., and Bengio, Y., “Random search for hyper-parameter optimization.” _Journal of Machine Learning Research_ , Vol. 13, No. 2, 2012\.
  * Elsken et al. [2018] Elsken, T., Metzen, J. H., and Hutter, F., “Neural architecture search: A survey,” _arXiv preprint arXiv:1808.05377_ , 2018.
  * Liu et al. [2021] Liu, Y., Sun, Y., Xue, B., Zhang, M., Yen, G. G., and Tan, K. C., “A survey on evolutionary neural architecture search,” _IEEE Transactions on Neural Networks and Learning Systems_ , 2021.
  * Holmes et al. [2012] Holmes, P., Lumley, J. L., Berkooz, G., and Rowley, C. W., _Turbulence, coherent structures, dynamical systems and symmetry_ , Cambridge University Press, Cambridge, 2012.
  * Rowley and Dawson [2017] Rowley, C. W., and Dawson, S. T., “Model reduction for flow analysis and control,” _Annual Review of Fluid Mechanics_ , Vol. 49, 2017, pp. 387–417.
  * Pathak et al. [2018] Pathak, J., Hunt, B., Girvan, M., Lu, Z., and Ott, E., “Model-free prediction of large spatiotemporally chaotic systems from data: A reservoir computing approach,” _Physical Review Letters_ , Vol. 120, No. 2, 2018, p. 024102.
  * Vlachas et al. [2018] Vlachas, P. R., Byeon, W., Wan, Z. Y., Sapsis, T. P., and Koumoutsakos, P., “Data-driven forecasting of high-dimensional chaotic systems with long short-term memory networks,” _Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences_ , Vol. 474, No. 2213, 2018, p. 20170844.
  * Maulik et al. [2020a] Maulik, R., Mohan, A., Lusch, B., Madireddy, S., Balaprakash, P., and Livescu, D., “Time-series learning of latent-space dynamics for reduced-order model closure,” _Physica D: Nonlinear Phenomena_ , Vol. 405, 2020a, p. 132368.
  * Ahmed et al. [2019] Ahmed, S. E., Rahman, S. M., San, O., Rasheed, A., and Navon, I. M., “Memory embedded non-intrusive reduced order modeling of non-ergodic flows,” _Physics of Fluids_ , Vol. 31, No. 12, 2019, p. 126602.
  * Wan et al. [2018] Wan, Z. Y., Vlachas, P., Koumoutsakos, P., and Sapsis, T., “Data-assisted reduced-order modeling of extreme events in complex dynamical systems,” _PLoS ONE_ , Vol. 13, No. 5, 2018, p. e0197704.
  * Hochreiter and Schmidhuber [1997] Hochreiter, S., and Schmidhuber, J., “Long short-term memory,” _Neural Computation_ , Vol. 9, No. 8, 1997, pp. 1735–1780.
  * Xie and Yuille [2017] Xie, L., and Yuille, A., “Genetic CNN,” _Proceedings of the IEEE international conference on computer vision_ , 2017, pp. 1379–1388.
  * Sun et al. [2020] Sun, Y., Xue, B., Zhang, M., Yen, G. G., and Lv, J., “Automatically Designing CNN Architectures Using the Genetic Algorithm for Image Classification,” _IEEE Transactions on Cybernetics_ , 2020.
  * Drozdzal et al. [2016] Drozdzal, M., Vorontsov, E., Chartrand, G., Kadoury, S., and Pal, C., “The importance of skip connections in biomedical image segmentation,” _Deep Learning and Data Labeling for Medical Applications_ , Springer, 2016, pp. 179–187.
  * He et al. [2016] He, K., Zhang, X., Ren, S., and Sun, J., “Deep residual learning for image recognition,” _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2016, pp. 770–778.
  * Maulik et al. [2020b] Maulik, R., Lusch, B., and Balaprakash, P., “Non-autoregressive time-series methods for stable parametric reduced-order models,” _Physics of Fluids_ , Vol. 32, No. 8, 2020b, p. 087115.