# Flow control by a hybrid use of machine learning and control theory

Takeru Ishize
Department of Mechanical Engineering, Keio University
<EMAIL_ADDRESS>
&Hiroshi Omichi
Department of Mechanical Engineering, Keio University;
Department of Mechanical and Aerospace Engineering, UCLA
<EMAIL_ADDRESS>
&Koji Fukagata
Department of Mechanical Engineering, Keio University
<EMAIL_ADDRESS>

###### Abstract

Flow control has a great potential to contribute to the sustainable society
through mitigation of environmental burden. However, high dimensional and
nonlinear nature of fluid flows poses challenges in designing efficient
control laws. This paper aims to propose a hybrid method (i.e., machine
learning and control theory) for feedback control of fluid flows. We propose a
partially nonlinear linear-system extraction autoencoder (pn-LEAE), which
consists of convolutional neural networks-based autoencoder (CNN-AE) and a
custom layer to extract a low-dimensional latent dynamics. This pn-LEAE
basically extracts a linear dynamical system so that the modern control theory
can easily be applied, but at the same time, it is designed to capture a
nonlinear development of the latent dynamics. We demonstrate the effectiveness
of the linear system extracted by the pn-LEAE, as well as the designed control
law’s effectiveness for a flow around a circular cylinder at the Reynolds
number of ${\rm Re}_{D}=100$. This is the first attempt utilizing CNN-AE for
linearization of fluid flows involving transient development to design a
feedback control law.
Keywords Machine learning, Autoencoder, Linear system extraction, Reduced
order modeling, Flow control

##  1 Introduction

Flow control has always been a significant focus in fluid mechanics, aiming at
improving flow characteristics such as reducing drag, enhancing mixing, and
minimizing noise, and the ability to manipulate fluid flow in a versatile
manner is not only a matter of scientific interest but also holds substantial
engineering importance. Flow control strategies can be broadly classified into
passive flow control (PFC) and active flow control (AFC). PFC typically
involves geometric modifications within the initial setup that does not
require external energy input in its operation. In contrast, AFC employs
actuators and requires external energy input in its operation. AFC is gaining
increased attention due to its potential for delivering more pronounced
effects over a broader operational range than PFC. AFC often relies on control
laws rooted in heuristic concepts such as the opposition control (Choi et
al.,, 1994) and superposition of traveling waves (Karniadakis and Choi,, 2003;
Ricco et al.,, 2021; Fukagata et al.,, 2024). At the same time, classical and
modern control theories have also been integrated into fluid mechanics,
providing a robust foundation for AFC controller design (Abergel and Temam,,
1990; Roussopoulos,, 1993; Kim and Bewley,, 2007).

There are primarily two categories of control theories in AFC: model-based
methods and model-free methods. Model-free methods derive control laws from
data in an unsupervised manner, enabling their application without the need
for modeling phenomena (Brunton and Noack,, 2015). These methods find
extensive use across numerical and experimental systems. In recent times,
machine learning (ML) techniques, including reinforcement learning and genetic
programming, have been used to discover control laws, demonstrating
significant potential for AFC (Rabault et al.,, 2019; Ren et al.,, 2019).
Nevertheless, ML methods are often treated as black boxes, making it
challenging to elucidate the underlying control logic. In contrast, model-
based methods derive control laws from system dynamics using mathematical
theories, encompassing linear and stochastic models, among others. These
methods have been applied to numerous flow control problems (Kim and Bewley,,
2007; Kasagi et al.,, 2009). However, modeling flow phenomena remains
difficult due to their high dimensionality and strong nonlinearity, and their
complexity has prevented the establishment of universal control methods.
Therefore, there is a strong need to develop a novel framework to address
these intricate issues.

To address these challenges, the utilization of ML techniques is under
extensive investigation. ML methods are data-driven and are renowned for their
adaptability across diverse cases, without necessitating mathematical models
of the target systems. Machine learning techniques have found broad
applications in the field of fluid mechanics (Brunton et al., 2020b, ; Fukami
et al.,, 2020; Brunton,, 2021), including state estimation (Milano and
Koumoutsakos,, 2002; Tong et al.,, 2022; Farazmand and Saibaba,, 2023),
reduced-order modeling (Taira et al.,, 2020), enhancement in numerical
simulations and experiments (Kochkov et al.,, 2021; Vinuesa and Brunton,,
2022; Vinuesa et al.,, 2023), turbulence modeling (Duraisamy et al.,, 2019),
and super-resolution reconstruction (Fukami et al.,, 2023) to name a few.
Among these, reduced-order modeling has gained significant attention because
of its capability to efficiently handle large datasets. In particular, the
autoencoder (AE) (Hinton and Salakhutdinov,, 2006), which maps high
dimensional fluid data into a low-dimensional feature space, is widely
utilized. AE is an unsupervised model that can be trained using the same data
before and after the bottleneck structure to extract low-dimensional features.
Hasegawa et al., (2020) applied an AE-based ROM to flows around cylinders at
various Reynolds numbers, demonstrating its applicability. San and Maulik,
(2018) introduced an extreme learning-based ROM for quasi-stationary
geophysical turbulent flows and demonstrated its superiority in terms of
stability when compared to the proper orthogonal decomposition (POD) (Lumley,,
1967)-based ROM. Various other studies have reported the use of AE-based ROMs
(Brenner et al.,, 2019; Brunton et al., 2020a, ; Fukagata,, 2023).
Consequently, the utilization of the AE-based approach for high-dimensional
fluid data has demonstrated the wide applicability of latent variables within
the extracted low-dimensional space.

However, the challenge of strong nonlinearity remains unresolved, and the low-
dimensional latent vectors obtained also exhibit nonlinearity. For instance,
Fukami et al., (2021) low-dimensionalized a flow around a circular cylinder at
the Reynolds number at ${\rm Re}_{D}=100$ using a Convolutional Neural Network
(LeCun et al.,, 1998) (CNN)-based AE (CNN-AE) and extracted the ordinary
differential equation (ODE) governing the low-dimensional latent dynamics
using SINDy (Brunton et al.,, 2016). While the dimension of the latent space
was reduced to two, the obtained ODE eligible of capturing the transient flow
was a set of fifth-order nonlinear equations, which is difficult to apply the
conventional control theory. Therefore, the present study aims to develop a
machine learning-based method, referred to as a Linear System Extracting
Autoencoder (LEAE), to extract a linear dynamical system in the latent space.
Although the structure of the LEAE is mostly similar to the conventional CNN-
AE, the LEAE is designed to output the flow field at the next time step
($t+\Delta t$) from that at the current time step ($t$). For that purpose, the
LEAE incorporates a custom layer which performs a time integration of the
latent dynamics based on the Crank-Nicholson scheme.

While the LEAE is expected to simplify the problem by mapping the high-
dimensional flow field governed by a nonlinear dynamics to a low-dimensional
latent variables governed by a linear dynamics, always assuming linear systems
poses difficulties when the nonlinearity plays an essential role, such as in
the cases of transient flows. To handle this issue, we also propose an
enhanced version of LEAE, called a Partially Nonlinear Linear System
Extraction Autoencoder (pn-LEAE), which can account for a nonlinear variation
of the orbital radius of the latent variables. This pn-LEAE, focusing on
continuously changing flow fields, extracts a system capable of representing
the time evolution of latent variables. In the present study, as a proof of
concept, we consider two types of flows: 1) transient flow and 2) steady
periodic flow around a circular cylinder at $Re_{D}=100$. This kind of low-
Reynolds number flow around a circular cylinder has long been considered in
the studies of feedback control aiming at suppression of vortex shedding (Park
et al.,, 1994); however, it has also been revealed that the linear feedback
control does not work at ${\rm Re}_{D}>80$ due to disappearance of the gain
window (Illingworth et al.,, 2014).

The paper is organized as follows. In Section 2, the method for preparing the
dataset, the details of the proposed LEAE, and the method of flow control
using the LEAE are introduced. In Section 3, we assess the linear system
extracted through the proposed procedure and demonstrate its effectiveness.
Then, we design a control law based on the extracted linear system and
investigate its performance. Finally, conclusions are drawn in Section 4.

##  2 Methods

###  2.1 Dataset

We consider the two-dimensional (2D) flow around a circular cylinder as a
basic example. Training data for machine learning (ML) model is prepared by
direct numerical simulation (DNS) (Kor et al.,, 2017), whose governing
equations are the incompressible continuity and Navier-Stokes equations, i.e.,

$\displaystyle{{\mbox{\boldmath$\nabla$}}}\cdot{\mbox{\boldmath$u$}}$
$\displaystyle=0,$ (1)
$\displaystyle\frac{\partial{\mbox{\boldmath$u$}}}{\partial
t}+{{\mbox{\boldmath$\nabla$}}}\cdot({\mbox{\boldmath$uu$}})$
$\displaystyle=-{{\mbox{\boldmath$\nabla$}}}p+\frac{1}{{\rm{Re}}_{D}}\nabla^{2}{\mbox{\boldmath$u$}}.$
(2)

Here, ${{\mbox{\boldmath$u$}}}=\\{u,v\\}$ and $p$ are the velocity vector and
pressure, respectively, nondimensionalized by the fluid density $\rho^{*}$,
the diameter of cylinder $D^{*}$, and the uniform velocity $U_{\infty}^{*}$.
The asterisk $(\cdot)^{*}$ denotes dimensional quantities, and the Reynolds
number is defined as ${\rm Re}_{D}=U^{*}D^{*}/\nu^{*}$, where $\nu^{*}$ is the
kinematic viscosity.

The left side of Figure 1 shows the computational domain, which is defined as
a rectangle with the streamwise $(x)$ and transverse $(y)$ dimensions of
$25.6D$ and $20.0D$, respectively. The center of the cylinder is located at a
distance of $9D$ from the inflow boundary. The present DNS code adopts the
second-order central finite difference method on the Cartesian coordinates,
and the no-slip condition on the cylinder surface is imposed using the
immersed boundary method (Kor et al.,, 2017). The number of computational
cells in the $x$ and $y$ directions is
$\left(N_{x},N_{y}\right)=\left(1024,800\right)$. The time step is set to
$\Delta t_{\rm DNS}=2.5\times 10^{-3}$. Free-slip conditions are applied to
the upper and lower boundaries, a uniform velocity $U_{\infty}=1$ is imposed
at the inflow boundary, and a convective boundary condition is utilized at the
outflow boundary. For more details of the present DNS and its verification and
validation, readers are referred to Kor et al., (2017).

As depicted on the right side of Figure 1, the locations for blowing and
suction are determined following Park et al., (1994), specifically at
$\theta=$$$ from the stagnation point. Time-varying blowing and suction is
applied in the radial direction with $v_{s}(t)=m(t)U_{\infty}$, where $m(t)$
represents the magnitude of the blowing and suction. Blowing is applied to the
upper slot and suction to the lower slot when $m(t)>0$, while the opposite
holds when $m(t)<0$. The tangential distribution of the velocity within the
blowing/suction slot is assumed Gaussian, such that the amplitude decreases to
$1\%$ of its peak at an angle of $\pm$$$ from the center. To construct ML
models, we utilize the velocity components $\left(u,v\right)$ within the
region delineated by the red line in Figure 1, which allows us to concentrate
on the flow around the cylinder. The dimensions of the instantaneous data are
$\left(384\times 192\times 2\right)$, i.e., $384\times 192$ computational
cells in $x$, $y$ directions and two velocity components ($u$, $v$) on each
computational point. As the dataset for the ML model, we generate two
different time series datasets of flow fields. The first dataset captures the
periodic wake behind the cylinder with steady blowing and suction. The
magnitude of blowing and suction is expressed as $m(t)=0,\pm 0.1,\pm 0.3,\pm
0.5,\pm 0.7,\pm 1.0$. Note that $m(t)=0$ is the base flow without blowing and
suction. We utilize data within the time range of $t=[0,100]$, with a time
interval of $\Delta t_{\rm data}=2.5\times 10^{-1}$. The second dataset
focuses on the transient evolution of the cylinder wake. To emphasize the
wake’s development, we employ data within the time span of $t=[35,350]$, with
a time interval of $\Delta t_{\rm data}=2.5\times 10^{-2}$.

Figure 1: Computational domain of DNS (black), the domain used for machine
learning (red), and a schematic of blowing and suction (right).

###  2.2 Partially linear, linear system extraction autoencoder (pn-LEAE)

####  2.2.1 Autoencoder

Figure 2: Schematic of the autoencoder.

Before delving into the methods we propose, it is important to provide an
explanation on the standard autoencoder. Typically, an autoencoder is used to
extract essential features from extensive datasets. The diagram of a typical
autoencoder is depicted in Figure 2. The autoencoder consists of two
components: an encoder $(\cal{F}_{\rm enc})$ and a decoder $(\cal{F}_{\rm
dec})$. The encoder’s function is to transform the input data $q$ into a
different representation called latent variables or latent vectors, denoted as
$r$. Meanwhile, the decoder’s role is to expand this representation to
reconstruct the original data, $\hat{{\mbox{\boldmath$q$}}}$, from the latent
vector $r$. These relationships can be expressed as

${\mbox{\boldmath$r$}}=\cal{F}_{\rm
enc}({\mbox{\boldmath$q$}})\in\mathbb{R}^{\mathit{n}_{\mathit{r}}},\\\ $ (3)
$\hat{{\mbox{\boldmath$q$}}}=\cal{F}_{\rm
dec}({\mbox{\boldmath$r$}})\in{\mathbb{R}}^{\mathit{n}_{\mathit{q}}},$ (4)

where $n_{r}$ represents the dimension of the latent vector, and $n_{q}$
represents the dimension of the original data. The encoder and decoder are
designed using neural network architectures such as Multi-Layer Perceptron
(Rumelhart et al.,, 1986) (MLP) and Convolutional Neural Network (LeCun et
al.,, 1998) (CNN). Their internal weights $w$ are fine-tuned during training
to ensure that the output data $\hat{{\mbox{\boldmath$q$}}}$ match the input
data $q$, i.e.,

$\displaystyle{{\mbox{\boldmath$w$}}}={\rm
argmin}_{{{\mbox{\boldmath$w$}}}}{\parallel{\mbox{\boldmath$q$}}-\hat{{\mbox{\boldmath$q$}}}\parallel}_{2}.$
(5)

If the original data is successfully reconstructed, it implies that the latent
vector $r$ can be seen as an alternative representation of the original data.
When the dimensions of the latent vectors are smaller than those of the
original data, the autoencoder can be employed for developing a reduced-order
model (ROM). A typical autoencoder only reduces the size of the data and does
not consider the temporal relationship of time-series data unless
intentionally designed so (Fukami et al.,, 2019).

####  2.2.2 Linear ordinary differential equation (LODE) layer

Let us explain the details of the linear ordinary differential equation layer
(LODE layer) we introduce in the present study. This customized layer contains
trainable parameters to numerically integrate a linear ordinary differential
equation (ODE). We are dealing with a system governed by a linear ODE
represented as

$\displaystyle\dot{{\mbox{\boldmath$r$}}}={{\mbox{\boldmath$A$}}}{{\mbox{\boldmath$r$}}}+{{\mbox{\boldmath$B$}}}{{\mbox{\boldmath$m$}}},$
(6)

where ${{\mbox{\boldmath$r$}}}\in\mathbb{R}^{\mathit{n}_{\mathit{r}}}$ is the
state of the system and
${{\mbox{\boldmath$m$}}}\in\mathbb{R}^{\mathit{n}_{\mathit{m}}}$ is the
control input. For instance, time integration of the left-hand-side of
equation (6) with the Euler explicit scheme yields

$\displaystyle{{\mbox{\boldmath$r$}}}(t+\Delta t)=\Delta
t({{\mbox{\boldmath$A$}}}{{\mbox{\boldmath$r$}}}(t)+{{\mbox{\boldmath$B$}}}m(t))+{{\mbox{\boldmath$r$}}}(t).$
(7)

The LODE layer calculates $r(t+\Delta t)$ from $r(t)$ and $\Delta t$ by the
operations like equation (7), although in practice we will use the Crank-
Nicolson scheme instead of the Euler explicit scheme, as will be explained
later, In the LODE layer, $A$ and $B$ are set as trainable parameters and are
optimized by training. Assuming that the dynamics of the latent vector is
described by a single oscillator, the matrix $A$ can be expressed as

${{\mbox{\boldmath$A$}}}=\begin{bmatrix}0&-\omega\\\ \omega&0\\\
\end{bmatrix},$ (8)

where the angular frequency $\omega$ is a trainable parameter. The training
process is divided into two steps. First, only the matrix
${{\mbox{\boldmath$A$}}}\ (\omega)$ is trained without control inputs. Then,
only the matrix $B$ is trained using the flow field with control inputs, while
$A$ is fixed. Eventually, the linear ODE represented by $A$ and $B$ serves as
the governing equations for the temporal evolution of $r$.

####  2.2.3 LODE for varying orbital radius

A previous study (Murata et al.,, 2020) demonstrated the 2D flow around a
circular cylinder at ${\rm Re}_{D}=100$ can be mapped to a binary latent
vector. Moreover, the acquired latent vectors were found to exhibit periodic
oscillations. Building on this background, we hypothesize that the governing
equation for the low-dimensional features of the present flow around a
circular cylinder can be simplified to a single oscillator, i.e., equation
(8). However, linear systems face challenges when applied to nonlinear flow
development. To address this, we introduce an enhanced Linear ODE (LODE) layer
that can handle varying orbital radius.

The idea for this enhanced LODE is borrowed from the error in the growth rate
of the radius when the non-conservative time integration schemes are used. For
instance, let us consider the oscillator at the unit angular frequency and the
state $r$ at time $t$ given as

${{\mbox{\boldmath$A$}}}=\begin{bmatrix}0&-1\\\ 1&0\\\
\end{bmatrix},\;\;{{\mbox{\boldmath$r$}}}(0)=\begin{bmatrix}\cos{\theta}\\\
\sin{\theta}\\\ \end{bmatrix}.$ (9)

The state $r$ and the orbital radius $R$ at $t+\Delta t$ computed using the
Crank-Nicholson method are

$\displaystyle{{\mbox{\boldmath$r$}}}(t+\Delta t)$
$\displaystyle=\frac{1}{4+\Delta t^{2}}\begin{bmatrix}(4-\Delta
t^{2})\cos{\theta}-4\Delta t\sin{\theta}\\\ 4\Delta t\cos{\theta}+(4-\Delta
t^{2})\sin{\theta}\\\ \end{bmatrix},$ (10) $\displaystyle R(t+\Delta t)$
$\displaystyle=\parallel{{\mbox{\boldmath$r$}}}(t+\Delta t)\parallel_{2}=1.$
(11)

In this case, the orbital radius $R$ is unchanged, consistent with its
mathematical property. On the other hand, when the Euler explicit method is
used, we can easily compute that the orbital radius changes to $R(t+\Delta
t)=\sqrt{1+\Delta t^{2}}$. It can also easily be shown that the desired rate
of change in the orbital radius, i.e., $\gamma\simeq(R(t+\Delta
t)-R(t))/\Delta t$, can be obtained by time integration using the Euler
explicit method by replacing the time step $\Delta t$ by a new parameter
$\beta$ defined as

$\beta=-\frac{2\gamma}{R(t)^{2}-\gamma^{2}}.$ (12)

Finally, by adding this term representing the growth rate to the Crank-
Nicholson method, we obtain a time integration scheme that can freely
accommodate the evolution of the state, i.e.,

$\displaystyle{{\mbox{\boldmath$r$}}}(t+\Delta t)$
$\displaystyle={{\mbox{\boldmath$r$}}}(t)+\Delta
t{{\mbox{\boldmath$A$}}}\left(\frac{{{\mbox{\boldmath$r$}}}(t+\Delta
t)+{{\mbox{\boldmath$r$}}}(t)}{2}\right)-\beta{{\mbox{\boldmath$A$}}}\left(\frac{{{\mbox{\boldmath$r$}}}(t+\Delta
t)-{{\mbox{\boldmath$r$}}}(t)}{2}\right),$ (13)

which can be explicitly expressed as

$\displaystyle{{\mbox{\boldmath$r$}}}(t+\Delta t)$
$\displaystyle=(2{{\mbox{\boldmath$I$}}}-(\Delta
t+\beta){{\mbox{\boldmath$A$}}})^{-1}(2{{\mbox{\boldmath$I$}}}+(\Delta
t-\beta){{\mbox{\boldmath$A$}}}){{\mbox{\boldmath$r$}}}(t).$ (14)

In the training process, $\gamma$ is set as a trainable parameter during flow
development. According to equation (12), $\beta$ is non-zero and computed by
equation (12), while $\beta=0$ in steady periodic state. More concretely, the
time integration performed by the LODE layer is summarized as

$\begin{cases}{{\mbox{\boldmath$r$}}}(t+\Delta
t)&=(2{{\mbox{\boldmath$I$}}}-(\Delta
t+\beta){{\mbox{\boldmath$A$}}})^{-1}(2{{\mbox{\boldmath$I$}}}+(\Delta
t-\beta){{\mbox{\boldmath$A$}}}){{\mbox{\boldmath$r$}}}(t)\hskip
28.45274pt(t<t_{b}),\\\ {{\mbox{\boldmath$r$}}}(t+\Delta
t)&=(2{{\mbox{\boldmath$I$}}}-\Delta
t{{\mbox{\boldmath$A$}}})^{-1}(2{{\mbox{\boldmath$I$}}}+\Delta
t{{\mbox{\boldmath$A$}}}){{\mbox{\boldmath$r$}}}(t)\hskip 81.09052pt(t_{b}\leq
t),\end{cases}$ (15)

where $t_{b}$ is the prescribed boundary between the transient and the steady
state.

####  2.2.4 Overview of partially nonlinear LEAE

Figure 3: Schematic of the partially nonlinear linear system extraction
autoencoder (pn-LEAE).

In this section, we overview the partially nonlinear LEAE (pn-LEAE), which is
built by combining the methods introduced in Sections 2.2.1, 2.2.2 and 2.2.3.
The schematic of the pn-LEAE is shown in Figure 3. The model is simultaneously
trained for reconstruction and time evolution, and the weights in the two
decoders are shared. The operations with time evolution by the LODE layer are
expressed as

$\displaystyle\begin{cases}{{\mbox{\boldmath$r$}}}(t)&=\mathcal{F}_{\rm
enc}({{\mbox{\boldmath$q$}}}(t)),\\\ {{\mbox{\boldmath$r$}}}(t+\Delta
t)&=\mathcal{F}_{\rm ODE}({{\mbox{\boldmath$q$}}}(t),\Delta t,m(t)),\\\
\hat{{\mbox{\boldmath$q$}}}(t+\Delta t)&=\mathcal{F}_{\rm
dec}({{\mbox{\boldmath$r$}}}(t+\Delta t)),\\\
\hat{{\mbox{\boldmath$q$}}}(t)&=\mathcal{F}_{\rm
dec}({{\mbox{\boldmath$r$}}}(t)),\end{cases}$ (16)

where $\Delta t$ is the time interval, $m(t)$ is the control input (a scalar
in this case). The encoder and decoder are trained concurrently to compress
low-dimensional features into values that align with a linear system, while
also improving reconstruction accuracy.

Furthermore, based on the assumption that the governing equation within the
latent space adheres to a single oscillation equation in the case of steady
periodic flow, it is expected that performance can be improved by imposing
constraints on its trajectory radius. Therefore, we add a regularization term
to the loss function,

$L_{\mathrm{d}is}=\sum(\max{|{{\mbox{\boldmath$r$}}}(t)|}-\min{|{{\mbox{\boldmath$r$}}}(t)|}),$
(17)

where $L_{\rm dis}$ is distance from origin. This constraint is applied to the
steady periodic flow and is performed on each of the mini-batches separated
for the training process.

The network structure of the present pn-LEAE is basically similar to the CNN-
AE used by Murata et al., (2020), with the differences that the aforementioned
LODE layer is inserted in the bottleneck part and that the latent vector is
decoded by two decoders for two consecutive time steps as illustrated in
Figure 3. For further details on the encoder and decoder parts, readers are
referred to Murata et al., (2020).

###  2.3 Control design for linear system

We employ the Linear Quadratic Regulator (LQR) (Anderson and Moore,, 2014) to
formulate a state feedback control law for the extracted linear systems. LQR
is a methodology used to devise an optimal control strategy for systems
governed by linear ODEs while minimizing costs. Now, let us elucidate the
process of designing a state feedback control law for the system outlined in
equation (6).

The quadratic cost function $J$ for the latent vector $r$ and the control
input $m$ is defined as

$J=\int_{0}^{\infty}\left[{{\mbox{\boldmath$r$}}}^{T}(t){{\mbox{\boldmath$Q$}}}{{\mbox{\boldmath$r$}}}(t)+{{\mbox{\boldmath$m$}}}^{T}(t){{\mbox{\boldmath$H$}}}{{\mbox{\boldmath$m$}}}(t)\right]dt,$
(18)

where $Q$ and $H$ is positive definite matrices with size $(n_{r}\times
n_{r})$ and $(n_{m}\times n_{m})$, respectively. The control input $m$ to
minimize the $J$ is derived as

${{\mbox{\boldmath$m$}}}(t)=-{{\mbox{\boldmath$H$}}}^{-1}{{\mbox{\boldmath$BPr$}}}(t),$
(19)

where $P$ is the target matrix obtained by solving the Riccati equation,

${{\mbox{\boldmath$A$}}}^{T}{{\mbox{\boldmath$P$}}}+{\mbox{\boldmath$PA$}}+{{\mbox{\boldmath$Q$}}}-{\mbox{\boldmath$PBH$}}^{-1}{{\mbox{\boldmath$B$}}}^{T}{{\mbox{\boldmath$P$}}}=0.$
(20)

Note that, in the present study, both the system matrix $A$ and the control
matrix $B$ are identified by machine learning using the LEAE.

##  3 Results and discussion

###  3.1 Linear system extraction

Figure 4: Dynamics of the latent vector: $(a)$ time trace of each latent
variable, drawn with the boundary between transient and steady periodic
states, $t_{b}$; $(b)$ trajectory of the latent vector. Reference represents
the dynamics of latent vector simply extracted by the CNN-AE, and pn-LEAE
represents those computed by time integration of the system identified by the
pn-LEAE. Figure 5: Comparison between the reference DNS and the flow fields
reproduced by the pn-LEAE: $(a)$ $t=40$ (transient process); $(b)$ $t=100$
(steady periodic state). Figure 6: Dynamics of the latent vector in the
presence of steady blowing and suction: $(a)$ time trace of each latent
variable, $m(t)=0.3$; $(b)$ trajectory of the latent vector, $m(t)=0.3$; $(c)$
time trace of each latent variable, $m(t)=1.0$; $(d)$ trajectory of the latent
vector, $m(t)=1.0$. Figure 7: Comparison between the reference DNS and the
flow fields reproduced by the pn-LEAE in the presence of steady blowing and
suction: $(a)$ $m(t)=0.3$; $(b)$ $m(t)=1.0$.

In this section, we focus on the linear system extraction by the proposed pn-
LEAE and discuss its performance. First, we demonstrate its applicability to a
flow field involving a transition process. The linear ODE governing the latent
dynamics obtained by the present pn-LEAE is

$\displaystyle\dot{{\mbox{\boldmath$r$}}}(t)=\begin{bmatrix}0&-1.066\\\
1.066&0\\\ \end{bmatrix}{{\mbox{\boldmath$r$}}}(t),$ (21)

with the growth rate in the transient period ($t<t_{b}$) of
$\gamma=1.332\times 10^{-3}$.

Figure 4 shows the predicted results of the temporal evolution of the latent
vectors using the extracted ODEs, along with the boundary distinguishing
between the time intervals where the growth rate is considered or not. The
prediction results align reasonably well with the reference, and the model can
accommodate both the evolution of the latent vectors and the steady periodic
oscillation. In particular, focusing on the period from the transient state to
the steady periodic state, we observe that the evolution of the latent space
comes to a halt, and a gradual transition takes place towards the steady
periodic behavior. These results also suggest that the behavior of the latent
vectors can be described by a simple linear ODE as expressed in equation (21)
with the help of the additional trainable parameter $\gamma$ introduced.

The flow field reconstructed by the decoder from the predicted latent vectors
is depicted in Figure 5. Here, we compare two states: $(a)$ the transient
state and $(b)$ the steady state. Both flow fields closely align with the DNS
data, highlighting the suitability of the temporal evolution of the latent
vector as described by the extracted linear system. This finding is a crucial
indicator of the model’s effectiveness.

Then, we investigate the applicability of the method to flow fields with
prescribed constant control inputs, i.e., $m(t)={\rm const.}$. Here, we
extract a linear system using equation (7) while fixing the weights of matrix
$A$. In other words, the training optimizes only matrix $B$ with respect to
the control inputs. The obtained linear system with the constant control
inputs is

$\displaystyle\dot{{\mbox{\boldmath$r$}}}(t)=\begin{bmatrix}0&-1.066\\\
1.066&0\\\ \end{bmatrix}{{\mbox{\boldmath$r$}}}(t)+\begin{bmatrix}-0.04039\\\
0.01019\\\ \end{bmatrix}{m}(t).$ (22)

Figure 6 shows the prediction results of the latent vectors based on the flow
field under two different control input magnitudes: $m(t)=0.3$ and $1.0$.
While a slight phase shift is observed for both control inputs, the temporal
change in the orbital radius is well captured, showcasing the effectiveness of
the present linear system. The phase-shift problem is an inherent challenge
arising from the accumulation of small errors within the linear system and
remains an important subject for future research.

In Figure 7, we observe the reconstructed fields alongside the reference data
obtained from DNS. Both predictions demonstrate the response of to the steady
blowing and suction, indicating that the model effectively capture the flow
development, including control effects, through a straightforward addition
within the latent space.

###  3.2 Design of control laws using the derived linear system

Figure 8: Control effect when the LQR regulator is applied in the latent
space: $(a)$ $H=0.1$; $(b)$ $H=0.3$ Figure 9: Flow fields reconstructed from
the latent dynamics controlled by the LQR regulator, i.e., equation (22).

In this section, we elucidate the method used to derive control laws from the
linear system obtained in Section 3.1. We also employ optimal control,
following the process described in Section 2.3. Specifically, we design a
state feedback control law for the linear system described by equation (22)
using LQR. It is essential to note that the rank of the controllability
matrix, denoted as $U_{b}$, is equal to $2$, which corresponds to $n_{r}$ in
this system, indicating its controllability. The matrices for LQR are set to
${{\mbox{\boldmath$Q$}}}={{\mbox{\boldmath$I$}}}$ and
${{\mbox{\boldmath$H$}}}=0.1{{\mbox{\boldmath$I$}}}$. Solving equation (20) by
using the Schur decomposition method, we can obtain $P$ as

$\displaystyle{{\mbox{\boldmath$P$}}}=\begin{bmatrix}10.31&0.1857\\\
0.1857&9.450\\\ \end{bmatrix}.$ (23)

Hence, the optimal control input to minimize cost function ${J}$ is determined
by equation (19) as,

$\displaystyle
m(t)=\begin{bmatrix}0.3737&-0.2457\end{bmatrix}{{\mbox{\boldmath$r$}}}(t).$
(24)

Before implementing this control law in DNS, we first apply this feedback
control law to the linear system to confirm its performance within the latent
space. The control results for $r$ are displayed in Figure 8($a$). The state
converges to the origin as a result of the control intervention. Although not
shown here, the control input $m(t)$ is found to fall within the range of
$[-1.0,1.0]$, which is within the maximum and minimum values of the control
input seen in the training data. For comparison, the result obtained with a
different weight for the control input, i.e.,
${{\mbox{\boldmath$H$}}}=0.3{{\mbox{\boldmath$I$}}}$ is also presented in
Figure 8($b$). The result is basically similar to the case of
${{\mbox{\boldmath$H$}}}=0.1{{\mbox{\boldmath$I$}}}$, but with a faster decay.

Figure 9 shows the flow fields reconstructed from the time histories of $r$.
It is evident that this behavior in the latent space corresponds to the
suppression of vortex shedding in the physical space. As a result, we can also
expect that this control law will successfully suppress vortex shedding when
applied to DNS.

###  3.3 Application of designed feedback control law to DNS

Figure 10: Schematic diagram of feedback control using the extracted control
law. Figure 11: Flow fields computed by DNS with the present feedback control
law. Figure 12: Time traces computed by DNS: $(a)$ lift coefficient $C_{L}$
without control (Reference) and with designed feedback control (LQRFeedback);
$(b)$ control input $m(t)$.

Using the derived control law (24), the control input in the DNS is expressed
as

$\displaystyle m(t)=\begin{bmatrix}0.3737&-0.2457\end{bmatrix}{\cal F}_{\rm
enc}({{\mbox{\boldmath$u$}}}(t)),$ (25)

where $m(t)$ represents the intensity of blowing and suction in the physical
space, and ${{\mbox{\boldmath$u$}}}(t)$ is the velocity field computed by DNS.
We apply this feedback control law to the DNS, keeping the DNS configuration
consistent with Section 2.1.

The procedure of the feedback control applied to the DNS is presented in
Figure 10. DNS is performed in the physical space with the control input
updated at prescribed time intervals ($\lambda$ time steps). In the latent
space, the flow field is mapped to latent vectors using the trained encoder,
and the control input is determined by applying the control law (25).
Repeating these procedures is expected to suppress vortex shedding. A crucial
aspect of this method is the connection between the latent space and the real
space through a control law represented by binary values.

Figure 11 shows the DNS results with the control input computed by the
procedure described above. The time interval for applying feedback control was
set to $\lambda=100$ time steps, which corresponds to approximately 1/25 of
the vortex shedding period. Over time, there is a noticeable reduction in
vortex shedding, indicating the success of the feedback control with the
designed control law. The control effect is aldo evident from the lift
coefficient depicted in Figure 12 $(a)$. Nevertheless, the issue of imperfect
control effectiveness persists. We are confronted with the challenge that the
control effect plateaus and does not further improve. The ultimate objective
of this study is to completely eliminate vortex shedding. as observed in the
latent space, also in the physical space. Analyzing the physical quantities
and control inputs presented in Figure 12, it becomes apparent that after
$t=80$, both remain within a specific range. To address this challenge, it
would be necessary to enhance the diversity of training data when training
autoencoders. However, this approach is not feasible due to computational
resource limitations. In the future, research efforts should focus on devising
strategies to improve control effectiveness while operating under various
constraints.

##  4 Conclusions

In this study, we introduced a Linear System Extraction Autoencoder (LEAE),
which is a linear system extraction model based on an autoencoder, designed to
address fluid flow phenomena characterized by high dimensionality and
nonlinearity. The proposed LEAE has shown that a low-Reynolds number flow
around a circular cylinder at ${\rm Re}_{D}=100$ can be effectively described
by a simple equation representing single oscillation. Furthermore, we
presented a partially nonlinear LEAE (pn-LEAE), an enhanced version of LEAE
capable of handling transient and feedback-controlled flows. A pn-LEAE has
demonstrated its ability to generalize and adapt to flow evolution by
incorporating a growth rate of the orbital radius. Finally, we designed a
control law based on the linear systems extracted by LEAE and provided
evidence of its effectiveness. This outcome underscores the capacity of the
designed control law to bridge the gap between latent space and real space.

## Acknowledgments

This work was supported by JSPS KAKENHI Grant Number 21H05007. The authors
acknowledge Mr. Shoei Kanehira (Keio University) for fruitful discussion and
comments.

## References

  * Abergel and Temam, (1990) Abergel, F. and Temam, R. (1990), “On some control problems in fluid mechanics,” Thoer. Comput. Fluid Dyn., Vol. 1, pp. 303–325.
  * Anderson and Moore, (2014) Anderson, B.D.O. and Moore, J.B. (2014), Optimal Control: Linear Quadratic Methods, Dover Publications, New York, pp. 35–37.
  * Brenner et al., (2019) Brenner, M.P., Eldredge, J.D. and Freund, J.B. (2019), “Perspective on machine learning for advancing fluid mechanics”, Physical Review Fluids, Vol. 4, Article No. 100501.
  * Brunton, (2021) Brunton, S. L. (2021), “Applying machine learning to study fluid mechanics”, Acta Mechanica Sinica, Vol. 37, pp. 1718–1726.
  * Brunton and Noack, (2015) Brunton, S.L. and Noack, B.R. (2015), “Closed-loop turbulence control: Progress and challenges”, Applied Mechanics Reviews, Vol. 67, Article No. 050801.
  * (6) Brunton, S.L., Hemati, M.S. and Taira, K. (2020a), “Special issue on machine learning and data-driven methods in fluid dynamics”, Theoretical and Computational Fluid Dynamics, Vol. 34, pp. 333–337.
  * (7) Brunton, S.L., Noack, B.R. and Koumoutsakos, P. (2020b) “Machine learning for fluid mechanics”, Annual Review of Fluid Mechanics, Vol. 52, pp. 477–508.
  * Brunton et al., (2016) Brunton, S.L., Proctor, J.L. and Kutz, J.N. (2016) “Discovering governing equations from data by sparse identification of nonlinear dynamical systems”, Proc. Natl Acad. Sci., Vol. 113, pp. 3932-–3937.
  * Choi et al., (1994) Choi, H., Moin, P. and Kim, J. (1994), “Active turbulence control for drag reduction in wall-bounded flows”, Journal of Fluid Mechanics, Vol. 262, pp. 75–110.
  * Duraisamy et al., (2019) Duraisamy, K., Iaccarino, G. and Xiao, H. (2019) “Turbulence modeling in the age of data”, Annual Review of Fluid Mechanics, Vol. 51, pp. 357–377.
  * Farazmand and Saibaba, (2023) Farazmand M. and Saibaba, A.K. (2023), “Tensor-based flow reconstruction from optimally located sensor measurements”, Journal of Fluid Mechanics, Vol. 962, Article No. A27.
  * Fukagata, (2023) Fukagata, K. (2023), “Reduced order modeling of fluid flows using convolutional neural networks,” J. Fluid Sci. Technol., Vol. 18, Article No. JFST0002.
  * Fukagata et al., (2024) Fukagata, K., Iwamoto, K. and Hasegawa, Y. (2024), “Turbulent drag reduction by streamwise traveling waves of wall-normal forcing,” Annu. Rev. Fluid Mech., Vol. 56, pp. 45-66.
  * Fukami et al., (2020) Fukami, K., Fukagata, K. and Taira, K. (2020), “Assessment of supervised machine learning for fluid flows”, Theoretical and Computational Fluid Dynamics, Vol. 34, pp. 497–519.
  * Fukami et al., (2023) Fukami, K., Fukagata, K. and Taira, K. (2023), “Super-resolution analysis via machine learning: A survey for fluid flows”, Theoretical and Computational Fluid Dynamics, Vol. 37, pp. 421–444.
  * Fukami et al., (2021) Fukami, K., Murata, T., Zhang, K. and Fukagata, K. (2021), “Sparse identification of nonlinear dynamics with low-dimensionalized flow representations”, J. Fluid Mech., Vol. 926, Article No. A10.
  * Fukami et al., (2019) Fukami, K., Nabae, Y., Kawai, K. and Fukagata, K. (2019), “Synthetic turbulent inflow generator using machine learning”, Physical Review Fluids, Vol. 4, Article No. 064603.
  * Hasegawa et al., (2020) Hasegawa, K., Fukami, K., Murata, T. and Fukagata, K. (2020), “CNN-LSTM based reduced order modeling of two-dimensional unsteady flows around a circular cylinder at different reynolds numbers”, Fluid Dynamics Research, Vol. 52, Article No. 065501.
  * Hinton and Salakhutdinov, (2006) Hinton, G.E. and Salakhutdinov, R.R. (2006), “Reducing the dimensionality of data with neural networks”, Science, Vol. 313, pp. 504–507.
  * Illingworth et al., (2014) Illingworth, S.J., Naito, H. and Fukagata, K. (2014), “Active control of vortex shedding: An explanation of the gain window”, Physical Review E, Vol. 90, Article No. 043014.
  * Karniadakis and Choi, (2003) Karniadakis, G.E. and Choi, K.-S. (2003). “Mechanisms on transverse motions in turbulent wall flows,” Annu. Rev. Fluid Mech., Vol. 35, pp. 45–62.
  * Kasagi et al., (2009) Kasagi, N., Suzuki, Y. and Fukagata K. (2009), “Microelectromechanical systems–based feedback control of turbulence for skin friction reduction”, Annual Review of Fluid Mechanics, Vol. 41, pp. 231–251.
  * Kim and Bewley, (2007) Kim, J. and Bewley, T.R. (2007), “A linear systems approach to flow control”, Annual Review of Fluid Mechanics, Vol. 39, pp. 383–417.
  * Kochkov et al., (2021) Kochkov, D., Smith, J.A., Alieva, A., Wang, Q., Brenner, M.P. and Hoyer, S. (2021), “Machine learning-accelerated computational fluid dynamics”, Proceedings of the National Academy of Sciences, Vol. 118, Article No. e2101784118.
  * Kor et al., (2017) Kor, H., Badri Ghomizad, M. and Fukagata, K. (2017), “A unified interpolation stencil for ghost-cell immersed boundary method for flow around complex geometries”, Journal of Fluid Science and Technology, Vol. 12, Article No. JFST0011–JFST0011.
  * LeCun et al., (1998) LeCun, Y., Bottou, L., Bengio, Y. and Haffner, P. (1998), “Gradient-based learning applied to document recognition”, Proceedings of the IEEE, Vol. 86, pp. 2278–2324.
  * Lumley, (1967) Lumley, J.L. (1967), “The structure of inhomogeneous turbulent flows Atmospheric Turbulence and Radio Wave Propagation, Nauka, Moscow, pp. 166–178.
  * Milano and Koumoutsakos, (2002) Milano, M. and Koumoutsakos, P. (2002), “Neural network modeling for near wall turbulent flow”, Journal of Computational Physics, Vol. 182, pp. 1–26.
  * Murata et al., (2020) Murata, T., Fukami, K. and Fukagata, K. (2020), “Nonlinear mode decomposition with convolutional neural networks for fluid dynamics”, Journal of Fluid Mechanics, Vol. 882, Article No. A13.
  * Park et al., (1994) Park, D.S., Ladd, D.M. and Hendricks, E.W. (1994), “Feedback control of von Kármán vortex shedding behind a circular cylinder at low reynolds numbers”, Physics of Fluids, Vol. 6, pp. 2390–2405.
  * Rabault et al., (2019) Rabault, J., Kuchta, M., Jensen, A., Réglade, U. and Cerardi, N. (2019), “Artificial neural networks trained through deep reinforcement learning discover control strategies for active flow control”, Journal of Fluid Mechanics, Vol. 865, pp. 281–302.
  * Ren et al., (2019) Ren, F., Wang, C. and Tang, H. (2019), “Active control of vortex-induced vibration of a circular cylinder using machine learning”, Physics of Fluids, Vol. 31, Article No. 093601.
  * Ricco et al., (2021) Ricco, P., Skote, M. and Leschziner M.A. (2021), “A review of turbulent skin-friction drag reduction by near-wall transverse forcing”, Prog. Aerosp. Sci., Vol. 123, Article No. 100713.
  * Roussopoulos, (1993) Roussopoulos, K. (1993), “Feedback control of vortex shedding at low Reynolds numbers”, Journal of Fluid Mechanics, Vol. 248, pp. 267–296.
  * Rumelhart et al., (1986) Rumelhart, D.E., Hinton, G.E. and Williams, R.J. (1986), “Learning representations by back-propagating errors”, Nature, Vol. 323, pp. 533–536.
  * San and Maulik, (2018) San O. and Maulik, R. (2018), “Extreme learning machine for reduced order modeling of turbulent geophysical flows”, Physical Review E, Vol. 97, Article No. 042322.
  * Taira et al., (2020) Taira, K., Hemati, M.S., Brunton, S.L., Sun, Y., Duraisamy, K., Bagheri, S., Dawson, S. and Yeh, C.-A. (2020), “Modal analysis of fluid flows: Applications and outlook”, AIAA Journal, Vol. 58, pp. 998–1022.
  * Tong et al., (2022) Tong, W., Wang, S. and Yang, Y. (2022), “Estimating forces from cross-sectional data in the wake of flows past a plate using theoretical and data-driven models”, Physics of Fluids, Vol. 34, Article No. 111905.
  * Vinuesa and Brunton, (2022) Vinuesa, R. and Brunton, S.L. (2022), “Enhancing computational fluid dynamics with machine learning”, Nature Computational Science, Vol. 2, pp. 358–366.
  * Vinuesa et al., (2023) Vinuesa, R., Brunton, S.L. and McKeon, B.J. (2023), “The transformative potential of machine learning for experiments in fluid mechanics”, Nature Reviews Physics, Vol. 5, pp. 536–545.