Neural networks in feedback for flow analysis, sensor placement and
  control
This work presents a novel methodology for analysis and control of nonlinear
fluid systems using neural networks. The approach is demonstrated on four
different study cases being the Lorenz system, a modified version of the
Kuramoto-Sivashinsky equation, a streamwise-periodic 2D channel flow, and a
confined cylinder flow. Neural networks are trained as models to capture the
complex system dynamics and estimate equilibrium points through a Newton
method, enabled by backpropagation. These neural network surrogate models
(NNSMs) are leveraged to train a second neural network, which is designed to
act as a stabilizing closed-loop controller. The training process employs a
recurrent approach, whereby the NNSM and the neural network controller (NNC)
are chained in closed loop along a finite time horizon. By cycling through
phases of combined random open-loop actuation and closed-loop control, an
iterative training process is introduced to overcome the lack of data near
equilibrium points. This approach improves the accuracy of the models in the
most critical region for achieving stabilization. Through the use of L1
regularization within loss functions, the NNSMs can also guide optimal sensor
placement, reducing the number of sensors from an initial candidate set. The
datasets produced during the iterative training process are also leveraged for
conducting a linear stability analysis through a modified dynamic mode
decomposition approach. The results demonstrate the effectiveness of
computationally inexpensive neural networks in modeling, controlling, and
enabling stability analysis of nonlinear systems, providing insights into the
system behaviour and offering potential for stabilization of complex fluid
systems.
Tarcísio Déda, William Wolf, Scott Dawson
# Neural networks in feedback for flow analysis, sensor placement and control

###### Abstract

This work presents a novel methodology for analysis and control of nonlinear fluid systems using neural networks. The approach is demonstrated on four different study cases being the Lorenz system, a modified version of the Kuramoto-Sivashinsky equation, a streamwise-periodic 2D channel flow, and a confined cylinder flow. Neural networks are trained as models to capture the complex system dynamics and estimate equilibrium points through a Newton method, enabled by backpropagation. These neural network surrogate models (NNSMs) are leveraged to train a second neural network, which is designed to act as a stabilizing closed-loop controller. The training process employs a recurrent approach, whereby the NNSM and the neural network controller (NNC) are chained in closed loop along a finite time horizon. By cycling through phases of combined random open-loop actuation and closed-loop control, an iterative training process is introduced to overcome the lack of data near equilibrium points. This approach improves the accuracy of the models in the most critical region for achieving stabilization. Through the use of L1 regularization within loss functions, the NNSMs can also guide optimal sensor placement, reducing the number of sensors from an initial candidate set. The datasets produced during the iterative training process are also leveraged for conducting a linear stability analysis through a modified dynamic mode decomposition approach. The results demonstrate the effectiveness of computationally inexpensive neural networks in modeling, controlling, and enabling stability analysis of nonlinear systems, providing insights into the system behaviour and offering potential for stabilization of complex fluid systems.

## 1 Introduction

Neural networks (NNs) are powerful tools for data-driven mathematical modelling, finding a range of applications in the field of fluid mechanics such as flow analysis, optimization and control (Gad-el Hak, 1996; Lee _et al._, 1997; Kutz, 2017; Brunton _et al._, 2020). Recent advances in experimental data sampling, the availability of data from increasingly complex simulations, and the widespread availability of libraries that facilitate the design and training of NNs has encouraged their rapid development and application for the purposes of data compression, model order reduction, flow analysis, turbulence modelling, and feature extraction, all leveraging the power of NNs as general nonlinear function approximators. In the field of turbulence modeling, deep neural networks have been used to improve the prediction capabilities of RANS closure
models (Ling _et al._, 2016; Maulik _et al._, 2019). In the area of flow estimation and reconstruction, Milano & Koumoutsakos (2002) presented one of the first uses of neural networks in fluid mechanics for reconstructing near-wall turbulent flow fields. More recently, for example, Fukami _et al._ (2019) employed convolutional neural networks (CNNs) to accurately reconstruct turbulent flows from coarse flowfield images. In the context of reduced order modeling, Lui & Wolf (2019) combined deep neural networks and proper orthogonal decomposition (POD) to build accurate models of unsteady flows involving transient features, rewriting the Navier-Stokes equations as a set of nonlinear ordinary differential equations. More recently, Miotto & Wolf (2023) employed vision transformers and CNNs to extract quantities of interest such as aerodynamic coefficients and skin friction from flowfield visualizations. These authors demonstrated that the neural network models based on input images were effective in interpolating and extrapolating between flow regimes, offering a promising alternative for sensors in experimental campaigns and for building surrogate models of complex unsteady flows.

In the field of flow control, different techniques based on training NNs are found in the literature to either build a model for enabling control design (Morton _et al._, 2018; Bieker _et al._, 2020) or to directly infer a suitable control law for a given task. Training of NN-based models for control design was conducted by Bieker _et al._ (2020), who applied limited sensor data supported by delay coordinates, and Morton _et al._ (2018), in an approach that modelled the dynamics of a cylinder flow with models trained to learn state mappings to span a Koopman invariant subspace. Regarding model-free approaches, reinforcement learning (RL) can be highlighted as a promising tool for automatic training of control strategies with applications to drag reduction, for example. In their work, Rabault _et al._ (2019), Ren _et al._ (2021) and Castellanos Garcia de Blas _et al._ (2022) developed RL strategies to train controllers able to reduce drag in a confined cylinder flow through a set of measurements from velocity probes. Their sensor/actuation setups are employed to test the techniques proposed in the current work. Similarly, Li & Zhang (2022) considered RL strategies to suppress vortex shedding utilising stability analyses and unstable equilibrium computation of the underlying system to obtain improved results. Approaches based on RL have also been successfully applied to complex turbulent flows in both experimental (Fan _et al._, 2020) and computational (Guastoni _et al._, 2023) setups, being able, in some cases, to stabilise the flows (Sonoda _et al._, 2023). Comprehensive reviews of RL for flow control can be found in Vignon _et al._ (2023) and Viquerat _et al._ (2022).

Models based on sparse identification of nonlinear dynamics (SINDy) (Brunton _et al._, 2016) have been combined with order reduction techniques, such as POD, as an alternative to black box NN models, being able to find simpler and more intuitive models based on a library of interpretable functions. The approach can be applied for building nonlinear models with control inputs, which can be leveraged for closed-loop control applications (Kaiser _et al._, 2018). On the topic of flow optimization through open-loop control, genetic algorithms can be a candidate tool for parameter search. In such case, studies have been performed to find the best active control setup to improve a fitness function that aims to reduce drag in a bluff body (Zigunov _et al._, 2022), to maximise the thrust vectoring angle in a supersonic jet (Zigunov _et al._, 2023), and to find the best open-loop rotation values in order to reduce drag and symmetrise the wake of the flow past a cluster of cylinders (Raibaudo _et al._, 2020).

Traditional feedback control design approaches often rely on well established control theory, whose literature provides techniques that can ensure optimality, robustness, adaptability or the achievement of specific performance characteristics such as settling (convergence) time. Such linear control methods have been applied for the control of a variety of flow problems (Sipp & Schmid, 2016), such as cavities (Rowley & Williams, 2006; Barbagallo _et al._, 2009), bluff body wakes (Illingworth, 2016; Flinois & Morgans, 2016), and unsteady aerodynamic systems (Brunton _et al._, 2014; Herrmann _et al._, 2022; Sedky _et al._, 2023). However, the broad application of such approaches in fluid mechanics is limited by the fact that they typically rely on linear
approximations of the dynamics near desired equilibrium points. Furthermore, building data-driven linear models of flows is also hindered by the fact that signals sampled from fluid flows can be governed by nonlinear phenomena such as chaos, limit cycles, multiple equilibria, switching, and hysteresis. In the case of unstable flows (more specifically globally unstable flows), datasets containing measurements at (approximately) linear regions near equilibrium can be totally absent. Alternatively, nonlinear control techniques, which include sliding-mode control (Shi _et al._, 2019; Baek _et al._, 2019), feedback linearization (Alyoussef & Kaya, 2019), backstepping control (Zulu & John, 2016), extremum seeking control (Deda & Wolf, 2022) and switched control (Egidio _et al._, 2022), can be applied to the task of controlling nonlinear systems, although stabilization relative to a natural equilibrium point may still require a suitable model with some level of accuracy near equilibrium. Differentiable nonlinear models built from data can be leveraged for model predictive control (MPC), which employs real-time optimization to a cost function within a finite horizon. Although MPC still presents similar limitations near the linear regions for some flows -such as difficulties in stabilizing the system- it has been successfully applied to different flow control problems (Bieker _et al._, 2020; Arbabi _et al._, 2018; Deda _et al._, 2023; Morton _et al._, 2018). For example, near stabilization of the flow past a cluster of cylinders is achieved by Maceda _et al._ (2021).

The determination of optimal or advantageous locations to place sensors to enable full state reconstruction has been explored through a variety of methods, including via gappy proper orthogonal decomposition (Willcox, 2006), a pivoted QR decomposition of an identified set of basis functions (Manohar _et al._, 2018), data-driven dynamical models with Kalman filters (Sashittal & Bodony, 2021; Graff _et al._, 2023), and shallow decoder neural networks (Williams _et al._, 2022). As well as state reconstruction, the problem of sensor placement arises in flow control problems, where the placement of both sensors and actuators are important, and can be optimised and analysed using linear control theory (Chen & Rowley, 2011; Jin _et al._, 2022; Freire _et al._, 2020). The sensor placement methodology employed in the present work differs from these previous methods, and utilises a modified neural network loss function that incorporates the L1 norm of node (and thus physical location) weights. The use of L1-norm-based loss functions has been applied to find sparse solutions in a range of flow analysis problems, such as to find a reduced set of DMD modes to represent flow dynamics (Jovanovic _et al._, 2014), and to find spatially-localised resolvent modes (Skene _et al._, 2022; Lopez-Doriga _et al._, 2023).

In the present work, we present a novel approach that enables nonlinear control design by utilizing neural network architectures to model both the system and control law. To achieve stabilization of complex systems exhibiting nonlinear dynamics, a neural network controller (NNC) is trained in closed loop with a neural network surrogate model (NNSM), which is built from data obtained from a real plant (Deda _et al._, 2023). The methodology is proposed as a candidate solution to overcome the aforementioned lack of measurements near the linear region of unstable nonlinear systems. The training framework also allows for optimal sensor placement by automatic selection from a set of candidate probes. Here, through penalization of model inputs, a sparse configuration is aimed that exclude the less relevant sensors. Furthermore, a byproduct consisting of datasets that contain measurements from approximately linear regions of the state space is leveraged for the data-driven estimation of equilibrium points--which are used as control setpoints--and for conducting a data-driven linear stability analysis. This novelty allows for approaching both problems without the need to modify the solver and enables the construction of linearised flow models involving complex geometries, as well as other nonlinear systems. Moreover, it potentially enables equilibrium estimation and stability analysis through data gathered from experimental applications. The capabilities of the proposed methodology are demonstrated through analysis and control of four different plants, namely the Lorenz system, a modified version of the Kuramoto-Sivashinsky equation, a streamwise-periodic 2D channel flow, and a confined cylinder flow.
The remaining of this work is organized as follows: Section 2 describes the tools and techniques introduced and employed along the work, while Section 3 describes in detail the proposed iterative approach for training the NNs. In Section 4, we describe the different systems/plants to which the techniques are applied. Section 5 then presents the results for each choice of plant. Finally, the benefits and limitations of the approaches presented are discussed in Section 6. Finally, the Appendix A describes the hyperparameters for building and training the neural networks for each case studied.

## 2 Methodology

### Neural network surrogate models

In this work, NNSMs are employed to represent the dynamics of complex systems, such as fluid flows. Given a fixed time step \(\Delta t\), we are interested in discrete systems of the form

\[\mathbf{x}_{k+1}=F(\mathbf{x}_{k},\mathbf{\mu}_{k})\, \tag{1}\]

where \(\mathbf{x}\) is the vector containing the system states, \(\mathbf{\mu}\) is the vector of control inputs and the subscripts represent the discrete time steps such that time \(t=k\Delta t\). The states vector \(\mathbf{x}\) consists of a set of sensor measurements from which the dynamics can be inferred using data-driven methods. This set of of measurements consists of a subset of the true plant full set of states.

From a dataset containing measurements of both \(\mathbf{x}\) and \(\mathbf{u}\) from the plant, for multiple time steps, a NNSM is trained to approximate \(F(\mathbf{x},\mathbf{u})\). A schematic of the neural network, with its inputs and outputs, is presented in figure 1. Here, the symbol \(\tilde{F}\) is used to denote the function represented by the NNSM, and thus a good fit would imply \(\tilde{F}\approx F\). The general approach is similar to that described in Deda _et al._ (2023). The network is composed of fully connected hidden layers which employ a rectified linear unit (ReLU) activation function, whereas a linear function is used for the output layer. The data for training and prediction are normalised so the average for each input and output is 0 and the standard deviation is 1. The loss function used to train this NNSM is discussed in section 2.2. We also point out that the present NNs have a simple complexity, being kept at a maximum of two hidden layers for the cases studied in this work (with further details provided in Appendix A).

### Sensor placement

Another idea explored and applied in the current work is related to the problem of sensor placement. There are several reasons to incorporate sensor placement into the proposed methodology. First, using a smaller number of sensors reduces the total number of training parameters in the NNSM, which in turn reduces the amount of training data required. Second, real-world flow control applications typically use limited real-time measurements, and selecting a small number of sensors that are optimal for control has practical advantages. Lastly, the incorporation of sensor placement into the NN training process allows for optimal sensor locations to be identified without relying on any prior knowledge or intuition of the system dynamics.

The approach proposed here consists of searching for a subset of relevant sensors among an initial set of candidates. To do so, we implement a layer of trainable parameters that are used to weight the NNSM inputs. These weights are penalised using an L1 regularization. The main objective is the production of a sparse layer excluding unnecessary measurements. Hence, after training, the weights below a given threshold (in absolute values) are truncated to 0. For the initialization, each weight value is set to 1.

To avoid these weights getting too low by the cost of growing hidden layer weights, an L2 regularization is applied to the hidden layers. This ensures that the input weights can only decrease if they are irrelevant to estimate the output, which is still the full set of states in a future
time instant. Figure 1 illustrates such configuration. The loss function

\[\mathcal{L}=\frac{1}{n}\mathbf{g}\cdot\mathbf{g}+r_{2}\mathbf{w}_{\text{h}}\cdot\mathbf{w}_{ \text{h}}+r_{1}\sum\lvert\mathbf{w}_{\text{s}}\rvert\, \tag{2}\]

is used for evaluation of the network convergence, where \(\mathbf{g}\) is the array of differences between the labelled data and the NNSM output, \(\mathbf{w}_{\text{h}}\) is the stacked weights array for the hidden layers, \(\mathbf{w}_{\text{s}}\) is the array of weights for the sparsity layer, \(n\) is the total number of candidate sensors, and \(r_{1}\) and \(r_{2}\) are the L1 and L2 regularization factors. In practice, \(r_{1}\) can be adjusted such that larger values tend to provide sparser configurations at the cost of possibly reducing the accuracy of the model. The parameter \(r_{2}\) can be increased to avoid the hidden layer weights to grow indefinitely larger to compensate for the sparsity layer weights becoming small.

In practice, the final model can be viewed as a function of a reduced vector \(\mathbf{x}_{\mathbf{r}}\) containing \(n_{r}\leq n\) values, corresponding to the nonzero components of the sparsity weights, \(\mathbf{w}_{\text{s}}\). Although this is not done in the current work, pruning could also be conducted to reduce the complexity of the neural network by excluding unnecessary connections (Zahn et al., 2022). Regardless of pruning, it is now possible to evaluate the NNSM as a function with the following form

\[\mathbf{x}_{k+1}=\hat{F}_{r}\left(\mathbf{x}_{rk},\mathbf{u}_{k}\right)\, \tag{3}\]

so the full state vector is reconstructed from the reduced set of measurements. Similarly, we can define a function that predicts the components of the state corresponding to identified sensor locations,

\[\mathbf{x}_{rk+1}=\tilde{F}_{rr}\left(\mathbf{x}_{rk},\mathbf{u}_{k}\right). \tag{4}\]

In short, \(\hat{F}_{r}\) and \(\hat{F}_{rr}\) are evaluated with the reduced set of states, but the former outputs the full set of states while the latter outputs the reduced set. Both \(\hat{F}_{r}\) and \(\hat{F}_{rr}\) are used for different goals in the present work, as it will be clear in the following sections.

### Estimation of equilibrium state

To tackle equilibrium estimation, a linearization of the NNSM is obtained through backpropagation. This linearization approach was first introduced by Deda et al. (2023) and, here, an improved technique is proposed. In the first step, the reduced states matrix \(\mathbf{A}_{\mathbf{rr}\left(n_{r}\times n_{r}\right)}\) is obtained

Figure 1: Schematic of the NNSM where the red (blue) nodes are related to weights penalised by the L2 (L1) regularization. A ReLU function is employed in the nodes with a white mark, while a linear function is used for the other nodes.
for the time invariant system of the form

\[\mathbf{d_{rk+1}} =\mathbf{A_{rr}}\mathbf{d_{rk}}+\mathbf{B_{r}}\mathbf{u}_{k}\, \tag{5}\] \[\mathbf{d_{rk}} =\mathbf{x_{rk}}-\mathbf{x_{r}}^{*}\, \tag{6}\]

where \(\mathbf{A_{rr}}(n_{r}\times n_{r})\) and \(\mathbf{B_{r}}(n_{r}\times m)\) are constant matrices obtained from the Jacobian of \(\tilde{F}_{rr}\) evaluated at a given equilibrium operating condition \(\mathbf{x_{r}}^{*}\) and \(\mathbf{u}^{*}\) as follows:

\[\nabla\tilde{F}_{rr}=\begin{bmatrix}\mathbf{A_{rr}}&\mathbf{B_{r}}\end{bmatrix}= \begin{bmatrix}\frac{\partial\tilde{F}_{rr}}{\partial\mathbf{x}}&\frac{\partial \tilde{F}_{rr}}{\partial\mathbf{u}}\end{bmatrix}=\begin{bmatrix}\frac{\partial \tilde{F}_{rr}}{\partial\mathbf{x}}&\frac{\partial\tilde{F}_{rr}}{\partial\mathbf{u} }&\frac{\partial\tilde{F}_{rr}}{\partial\mathbf{u}_{1}}&\dots&\frac{\partial \tilde{F}_{1}}{\partial\mathbf{u}_{m}}\\ \vdots&\ddots&\vdots&\vdots&\ddots&\vdots\\ \frac{\partial\tilde{F}_{rr}}{\partial\mathbf{x}_{1}}&\dots&\frac{\partial \tilde{F}_{rr}}{\partial\mathbf{x}_{nr}}&\frac{\partial\tilde{F}_{rr}}{\partial \mathbf{u}_{1}}&\dots&\frac{\partial\tilde{F}_{nr}}{\partial\mathbf{u}_{m}}\end{bmatrix}. \tag{7}\]

Here, \(\tilde{f}_{i}\) is the \(i\)-th element in the output of \(\tilde{F}\).

In this work, \(\mathbf{u}^{*}=0\) is considered for calculation of a natural equilibrium point of the system, i.e., the equilibrium point related to the uncontrolled system. To compute \(\nabla\tilde{F}_{rr}\), an initial guess \(\mathbf{x_{r}}_{0}^{*}\) is used. The first estimate for the matrix \(\mathbf{A_{rr}}\) is obtained by taking the first \(n_{r}\) columns of \(\nabla\tilde{F}_{rr}\), while the remaining columns compose \(\mathbf{B_{r}}\) (see equation 7). The Newton method is employed to search for an equilibrium point for \(\tilde{F}_{rr}\).

Since \(\mathbf{u}^{*}=0\), the system of the form

\[\mathbf{x_{rk+1}}=\mathbf{A_{rr}}\mathbf{x_{rk}}+\mathbf{b_{r}} \tag{8}\]

is considered, where \(\mathbf{b_{r}}\) is an invariant vector. In this context, \(\mathbf{b_{r}}\) is introduced to add a bias that shifts equilibrium from the origin. For a given equilibrium point \(\mathbf{x_{r}}^{*}\),

\[\mathbf{x_{r}}^{*}=\mathbf{A_{rr}}\mathbf{x_{r}}^{*}+\mathbf{b_{r}}\, \tag{9}\]

which gives

\[\mathbf{x_{r}}^{*}=(I-\mathbf{A_{rr}})^{-1}\mathbf{b_{r}}. \tag{10}\]

The first guess for \(\mathbf{b_{r}}_{0}\) is obtained by evaluating

\[\mathbf{b_{r}}_{0}=\tilde{F}_{rr}(\mathbf{x_{r}}^{*},\mathbf{0})-\mathbf{A_{rr}}\mathbf{x_{r}}^{*} _{0}\, \tag{11}\]

which comes from the linear approximation of \(\tilde{F}_{rr}\). With \(\mathbf{b_{r}}\), \(\mathbf{x_{r}}^{*}\) can be updated by evaluating equation 10, which can subsequently be used to update \(\mathbf{b_{r}}\) via equation 11. The iterative process can be repeated to obtain an estimate of the equilibrium point. If desired, the equilibrium point can be mapped to the full set of sensors by evaluating

\[\mathbf{x}^{*}=\tilde{F}_{r}(\mathbf{x_{r}}^{*},\mathbf{u}_{k}). \tag{12}\]

The complete procedure is synthesised in Algorithm 1. Rather than iterating for a fixed number of steps, one could alternatively employ a stopping criteria based on the size of the update, or estimated error. It is important to mention that this approach is better suitable when the spectrum of \(\mathbf{A_{rr}}\) does not contain eigenvalues at 1+0j (integrator poles). In these cases, there is a continuous space of possible equilibrium points, any of which can be found by applying the Newton method. An alternative method for refining the equilibrium estimate, which is based on dynamic mode decomposition with control (DMDc) and which we observe can give improved results, will be described in Section 2.5.

### Neural network controller

Neural network controllers (NNC) can be employed for the task of controlling complex nonlinear systems (Deda _et al._, 2023). A nonlinear control law of the form

\[\mathbf{u}_{k}=K(\mathbf{x_{rk}}) \tag{13}\]
is implemented as a machine learning model, where the reduced set of states \(\mathbf{x}_{rk}\) are real-time measurements used to evaluate the control signal. A recurrent network is built for the training task as shown in figure 2. The model is trained such that a set of initial conditions is brought closer to the equilibrium point \(\mathbf{x}_{0}^{*}\) along iterations within a fixed discrete horizon \(n_{h}\).

The training data consists of a set of measurements \(\mathbf{X}_{0}=[\mathbf{x}_{r0,1}~{}\mathbf{x}_{r0,2}~{}\dots~{}\mathbf{x}_{r0,p}]\) containing \(p\) samples of the \(n_{r}\)-state system at an initial time \(k=0\). Note that the subscript \(0\) refers to the initial condition in the context of the finite-time closed-loop recurrent network, not being related to the simulation initial condition. The following loss function is targeted for minimization:

\[\mathcal{L} =\frac{1}{p\,n_{h}}\left(\frac{\mathbf{w}_{\mathbf{\theta}}}{n}\cdot\sum_ {j=1}^{p}\sum_{i=1}^{n_{h}}\mathbf{\phi}_{i,j}\circ\mathbf{\phi}_{i,j}+\frac{\mathbf{w}_{ \mathbf{u}}}{m}\cdot\sum_{j=1}^{p}\sum_{i=0}^{n_{h}-1}\mathbf{u}_{i,j}\circ\mathbf{u}_{i,j }\right)\,, \tag{14}\] \[\mathbf{\phi}_{i,j} =\mathbf{\tilde{x}}_{i,j}-\mathbf{x}^{*}\, \tag{15}\]

where \(p\) is the size of the training dataset. Here, the weight vector \(\mathbf{w}_{\mathbf{\theta}}\) has all elements equal to \(1/n_{r}\), which averages the result of the sum using an equal weighting of the states. Similarly, \(\mathbf{w}_{\mathbf{u}}\) has all elements equal to \(w_{\mathbf{u}}\), a scalar hyperparameter used to penalise the control inputs.

The hidden NNC layers are implemented with ReLU activation, except for the output layer that has a sigmoid activation, which is normalised to limit the control effort within the desired range. Once the training process is finished, the NNC block is used independently of the training loop to evaluate the control law in real time through equation 13. In the present work, all of the NNCs are implemented with a single hidden layer containing \(8\) nodes, which means they are computationally inexpensive as required for real-time applications. Appendix A further describes the hyperparameters employed.

### Linear stability analysis

Linear stability theory can be employed to elucidate mechanisms underlying the generation and amplification of flow instabilities. Here, we aim to employ the proposed neural network models to conduct linear stability analyses of unsteady flows, thereby illuminating physical mechanisms associated with unstable equilibrium points. In particular, with such an approach, it will be possible to compute the frequencies and growth rates of the most unstable modes and their associated eigenfunctions.

In an unsteady flow, a Reynolds decomposition allows the splitting of the flow states \(\mathbf{x}(\mathbf{x}_{\mathbf{c}},t)\) in a base flow \(\mathbf{x}^{*}(\mathbf{x}_{\mathbf{c}})\), here obtained at the equilibrium point, plus a time-dependent fluctuation component \(\mathbf{x}^{\prime}(\mathbf{x}_{\mathbf{c}},t)\). In the present notation, \(\mathbf{x}_{\mathbf{c}}\) represents the spatial coordinates. If the fluctuations are sufficiently small, the equations can be linearised about the base flow and the Navier-Stokes equations can then be written in the following linear form

\[\frac{\partial\mathbf{x}^{\prime}}{\partial t}=\mathcal{A}\mathbf{x}^{\prime}\, \tag{16}\]
where the matrix \(\mathbf{\mathcal{A}}=\mathbf{\mathcal{A}}(\mathbf{x}^{*})\) is a linear operator. For a system that is homogeneous in the streamwise (\(x_{c}\)) and spanwise (\(z_{c}\)) directions, the evolution of linear disturbances by equation 16 can be investigated by a direct analysis of the operator \(\mathbf{\mathcal{A}}\) with the transformation

\[\mathbf{x}^{\prime}(x_{c},y_{c},z_{c},t)=\int_{-\infty}^{\infty}\int_{-\infty}^{ \infty}\int_{-\infty}^{\infty}\hat{\mathbf{x}}(\alpha,y_{c},\beta,\omega)e^{ \mathrm{i}(\alpha x_{c}+\beta z_{c}-\omega t)}\;d\alpha\;d\beta\;d\omega\;, \tag{17}\]

where both the streamwise and spanwise wavenumbers \(\alpha,\beta\in\mathbb{R}\), and the frequency \(\omega\in\mathbb{C}\).

In discrete form, and under the transformation of equation 17, equation 16 can be written as

\[-\mathrm{i}\omega\hat{\mathbf{v}}=\mathbf{\mathcal{A}}\hat{\mathbf{v}}\;. \tag{18}\]

In this case, the linear operator becomes \(\mathbf{\mathcal{A}}=\mathbf{\mathcal{A}}(\mathbf{x}^{*},\alpha,\beta)\), and it can be analysed for each separate wavenumber pair \((\alpha,\beta)\) as an eigenvalue problem

\[\mathbf{V}\mathbf{\Lambda}=\mathbf{\mathcal{A}}\mathbf{V}\;. \tag{19}\]

In this equation, the columns \(\hat{\mathbf{v}}_{j}\) of \(\mathbf{V}\) are the eigenvectors of \(\mathbf{\mathcal{A}}\), and the eigenvalues \(\lambda_{j}=-\mathrm{i}\omega_{j}\) are the corresponding diagonal entries of \(\mathbf{\Lambda}\). Here the frequency and growth rate are the real and imaginary parts of \(\omega_{j}\), respectively.

There are several ways in which the identified NNSMs and NNCs can be leveraged to conduct stability analyses of the underlying dynamical systems. Most directly, we could utilise the matrix \(\mathbf{\mathcal{A}}_{\mathbf{tr}}\) arising from the linearization of the NNSM obtained through backpropagation (equation 7). Rather than linearising the global nonlinear model, we can alternatively identify a linear model directly from data taken near the equilibrium point, which is the method used in the present work. This approach amounts to performing a linear regression inspired in DMDc (Proctor _et al._, 2016). The idea comes from the fact that the controlled system can be slightly perturbed, producing a rich set of data near the equilibrium point, where the dynamics of the system tends to be approximately linear.

Given a dataset \(\mathbf{X}=[\mathbf{d}_{0}\dots\mathbf{d}_{p-1}]\) containing \(p\) consecutive full state measurements, a dataset \(\mathbf{U}=[\mathbf{u}_{0}\dots\mathbf{u}_{p-1}]\) containing the control inputs history, and a dataset \(\mathbf{X}^{\prime}=[\mathbf{d}_{1}\dots\mathbf{d}_{p}]\) containing

Figure 2: Schematic of the NNC training where the initial dataset is propagated through recurrent evaluations of the NNSM and NNC. Only the NNC weights and biases are updated during training in order to bring the states \(\mathbf{\chi}_{\mathbf{r}1\leqslant i\leqslant n_{h},j}\) closer to \(\mathbf{\chi}_{\mathbf{r}^{*}}\).
full state measurements with a unit shift, matrices \(\mathbf{A}\) and \(\mathbf{B}\) can be inferred as follows:

\[\mathbf{G} =\mathbf{X}^{\prime}\mathbf{\Omega}^{\dagger}\, \tag{20}\] \[\mathbf{G} \approx\begin{bmatrix}\mathbf{A}&\mathbf{B}\end{bmatrix}\,\] (21) \[\mathbf{\Omega} =\begin{bmatrix}\mathbf{X}\\ \mathbf{U}\end{bmatrix}. \tag{22}\]

As seen, matrix \(\mathbf{\Omega}\) is built from \(\mathbf{X}\) and \(\mathbf{U}\). By calculating its Moore-Penrose inverse \(\mathbf{\Omega}^{\dagger}\), it is possible to compute \(\mathbf{G}\), from which approximations for \(\mathbf{A}\) and \(\mathbf{B}\) can be extracted.

Note that \(\mathbf{X}\) is composed of the deviation of the states from their estimated equilibrium values, as shown in equation 6. Admitting that the estimation of this equilibrium \(\mathbf{X}^{*}\) is not exact, the proposed linear regression can be contaminated. To account for such imperfect estimation, a modification to this method is proposed, where a matrix \(\mathbf{O}=[1\ldots 1]_{(1\times p)}\) containing ones is employed. This allows for the modified procedure:

\[\mathbf{H} =\mathbf{X}^{\prime}\mathbf{\Psi}^{\dagger}\, \tag{23}\] \[\mathbf{H} \approx\begin{bmatrix}\mathbf{A}&\mathbf{B}&\mathbf{c}\end{bmatrix}\,\] (24) \[\mathbf{\Psi} =\begin{bmatrix}\mathbf{X}\\ \mathbf{U}\\ \mathbf{O}\end{bmatrix}. \tag{25}\]

For systems with no integrator poles, this modification allows for the estimation of a system of the form

\[\mathbf{d}_{k+1}=\mathbf{A}\mathbf{d}_{k}+\mathbf{B}\mathbf{u}_{k}+\mathbf{c}\, \tag{26}\]

where \(\mathbf{c}\) is a constant vector that displaces the equilibrium point of the new states \(\mathbf{d}\) from \(0\). Beyond providing a better estimation of \(\mathbf{A}\) to conduct stability analysis, the addition of the constant \(\mathbf{c}\) allows for the correction of the previous estimate of \(\mathbf{X}^{*}\). Although this is not done in the current work, it could also be used, for example, to compensate control loop measurements to improve NNC results. Since the NNC is trained to control the NNSM and, furthermore, \(\mathbf{x}^{*}\) is an equilibrium point for the NNSM, \(\mathbf{c}\) could be used to adapt the NNC loop for the real plant.

To reduce computational costs, two different sets of matrices are used based on the sparse configuration obtained as described in section 2.2:

* The first is obtained by using the reduced states dataset \(\mathbf{X}_{\mathbf{r}}=[\mathbf{d}_{\mathbf{r}0}\ldots\mathbf{d}_{\mathbf{r}p-1}]\). Replacing \(\mathbf{X}\) by \(\mathbf{X}_{\mathbf{r}}\), the algorithm produces a matrix \(\mathbf{A}_{\mathbf{r}(n\times n_{r})}\) which maps the reduced set of sensors to the full set at the next timestep, as well as \(\mathbf{B}\) and \(\mathbf{c}\).
* The second set of matrices is obtained by also using \(\mathbf{X}_{\mathbf{r}}\) with the further substitution of \(\mathbf{X}^{\prime}\) by a reduced matrix \(\mathbf{X}_{\mathbf{r}^{\prime}}\). This produces the same matrices \(\mathbf{A}_{\mathbf{r}(n_{r}\times n_{r})}\) and \(\mathbf{B}_{\mathbf{r}(n_{r}\times m)}\) shown in section 2.3--with potentially better accuracy--as well as the reduced bias vector \(\mathbf{c}_{\mathbf{r}}\).

A stability analysis can be conducted by finding the eigenvectors \(\mathbf{V}_{\mathbf{r}i}\) and eigenvalues \(\lambda_{i}\) from the square matrix \(\mathbf{A}_{\mathbf{r}\mathbf{r}}\), where \(i=1,2,\ldots,n_{r}\). Since each vector \(\mathbf{V}_{\mathbf{r}i}\) is an eigenvector of \(\mathbf{A}_{\mathbf{r}\mathbf{r}}\), each \(\mathbf{A}_{\mathbf{r}\mathbf{r}}\mathbf{V}_{\mathbf{r}i}\) is also an eigenvector of \(\mathbf{A}_{\mathbf{r}\mathbf{r}}\), which allows the mapping

\[\mathbf{V}_{\mathbf{r}(n\times 1)}=\mathbf{A}_{\mathbf{r}(n\times n_{r})}\,\mathbf{V}_{\mathbf{r}i(n_{ r}\times 1)} \tag{27}\]

to estimate the eigenvectors across all sensor states. This mapping is employed in this work for visualization of modes projected over the full set of candidate sensors.

## 3 Iterative training

A key aspect of the proposed methodology is the development of a training methodology that ensures that both the NNSM and NNC are trained with sufficient quantities of data, particularly
in the regions of state space where model accuracy is important for control effectiveness. To train and improve the neural network models employed in this work, an iterative training algorithm is proposed. Four different modes of data sampling are considered for training the NNSM and the NNC. Each mode produces data that are compiled for training both networks. They are:

* _Uncontrolled sweep_ mode: the system is perturbed by an open-loop control signal with no closed-loop control;
* _Controlled sweep_ mode: the system is controlled in closed loop at the same time that it is perturbed by an open-loop signal;
* _Release_ mode: the plant runs with no control and no open-loop perturbation, allowing for better visualization of the next mode and to test whether the system can be successfully controlled from an initially uncontrolled state;
* _Control_ mode: the system is controlled in closed loop with no open-loop perturbation.

All open-loop perturbations applied in this work are random staircase signals with constant step length.

The _uncontrolled sweep_ is the first mode employed so as to produce data with control input history. A number of samples is gathered from the perturbed plant. These first samples are used to compute the means and standard deviations used for data normalization of the neural networks inputs and outputs. This data is also used for training the NNSM, estimating its equilibrium point \(\boldsymbol{x_{r}}^{*}\) and training the NNC, all for the first time. The NNSM sparsity layer has its negligible weights truncated according to a threshold (further described in Appendix A), which accounts for identifying and updating the relevant sensor locations.

The second mode to be run is the _release_ mode, which allows the system to achieve its natural (uncontrolled) behaviour, i.e., a limit cycle or a chaotic attractor, for example. The _control_ mode is applied after the system completely or partially returns to its natural operation. When the NNC is turned on, the dynamics of the controlled system drives the plant states to work under a new behaviour, which can tend to a new limit cycle, a new chaotic attractor or even to an equilibrium point - which can either be or not be desirably close to the estimated equilibrium, \(\boldsymbol{\chi}^{*}\).

Data is gathered in each of these three steps, building a database that grows along the iterative process. After they are run, the process is repeated for a number of iterations, but instead of doing an _uncontrolled sweep_, the _controlled sweep_ mode is applied in every iteration beyond the first one. The idea is to leverage the fact that the controlled system should have a tendency to get closer to the equilibrium point, and by perturbing the controlled system appropriately, a larger set of data is obtained around this equilibrium. When a new iteration begins, the amplitudes of the open-loop perturbations are reduced by a constant factor \(0<\alpha\leq 1\), so that data is gathered closer to the actual plant equilibrium. At each iteration, the NNSM and the NNC are retrained with new data collected from the full database and the equilibrium point is updated using the retrained NNSM. Algorithm 2 synthesises the steps described, and a schematic showing the phase portrait for the Lorenz system is also presented in figure 3. Each iteration in this figure depicts 1000 measurements sampled during the _sweep_ modes (_uncontrolled_ for iteration 1, _controlled_ for subsequent iterations). As the iterations advance, the closed-loop control tends to bring states closer to the origin, which is an equilibrium point for the Lorenz system. This enhances the estimate of the equilibrium point, and it can be seen in the smaller limits of the axes, as shown for the higher iterations. Also, beginning in iteration 7, since the perturbations are reduced between iterations, they become unable to shoot the controlled states out of the main stable orbit around the now-stabilised equilibrium.
```
1:for\(i=1\) to \(n\)do
2:if i==1 then
3: Run uncontrolled sweep mode to gather data
4: Compute parameters for normalization
5:else
6: Reduce the amplitude of open-loop perturbations by a factor \(\alpha\)
7: Run controlled sweep mode to gather data around the estimated equilibrium
8:endif
9: Update the database with the new data
10: Train NNSM with the updated database
11: Truncate negligible weights in the sparsity layer
12: Estimate equilibrium from NNSM
13: Train NNC with the updated database
14: Run release mode to allow the system to reach its natural behaviour
15: Gather data in release mode
16: Run control mode to drive the system to a new behaviour and test control effectiveness
17: Gather data in control mode
18:endfor
```

**Algorithm 2** Data Sampling Algorithm

Figure 3: Phase portraits of the Lorenz system for data sampled during the _sweep_ modes at different iterations. Note that for iterations 7 to 9 the axes limits are considerably smaller. The red dot indicates the origin, which is an equilibrium point for the Lorenz system.
For the stability analysis, only the data produced by the last _sweep_ mode is used since it may contain more accurate data near the equilibrium point, hopefully ruled by approximately linear dynamics. One of the advantages of using such type of data is enabling DMD-based techniques to identify structures corresponding to eigenvectors of the equillibirum-linearised system, even for systems exhibiting strong nonlinearities.

## 4 Study cases

This section outlines the four systems that will be used to test and validate the proposed methodology.

### Lorenz equations

The first study case for this work is the Lorenz system, whose dynamics are described by the following set of nonlinear ordinary differential equations

\[\dot{x} = \sigma(y-x)\, \tag{1}\] \[\dot{y} = \rho x-y-xz+u\,\] (2) \[\dot{z} = -\beta z+xy\, \tag{3}\]

where \(\sigma=10\), \(\beta=8/3\) and \(\rho=28\) are chosen for the present analysis. With this set of parameters, the uncontrolled system (\(u(t)\)=0) behaves as a chaotic attractor. The goal is to apply the techniques proposed in this work to find an equilibrium point; to stabilise the system, bringing the states towards this point; and to compute a linear model to obtain the system poles. In this problem we do not conduct sensor placement.

### Modified Kuramoto-Sivashinsky equation

The Kuramoto-Sivashinsky (KS) equation is given by the following nonlinear partial differential equation

\[\frac{\partial v}{\partial t}+v\frac{\partial v}{\partial x_{c}}=-\frac{1}{R} \left(P\frac{\partial^{2}v}{\partial x_{c}^{2}}+\frac{\partial^{4}v}{ \partial x_{c}^{4}}\right)\, \tag{4}\]

where \(R\) is equivalent to the Reynolds number in a fluid flow, and \(P\) represents a balance between energy production and dissipation. This system has been used to emulate, for example, combustion instabilities in flame fronts (Sivashinsky 1977; Kuramoto 1978) and hydrodynamic instabilities in shear flows (Fabbiane _et al._ 2014). In this context, \(x_{c}\) is the spatial coordinate and \(v(x_{c})\) is the variable of interest.

It is clear that, for any constant \(c\), \(v(x_{c})=c\) corresponds to an equilibrium point, i.e., \(\partial v/\partial x_{c}=0\). Since we are interested--among other goals--in computing an equilibrium point for this system, we propose a modified version of the KS equation given by

\[\frac{\partial v}{\partial t}+v\frac{\partial v}{\partial x_{c}}=-\frac{1}{R} \left(P\frac{\partial^{2}v}{\partial x_{c}^{2}}+\frac{\partial^{4}v}{ \partial x_{c}^{4}}\right)-\frac{Q}{L}\int_{0}^{L}\left(v(x_{c})-V\right)dx_{c }+\sum_{i=0}^{m}B_{i}(x_{c})u_{i}. \tag{5}\]

The new term \(-Q/L\int_{0}^{L}(v(x_{c})-V)\ dx_{c}\) corresponds to a spatial average of the difference between \(v(x_{c})\) and a chosen equilibrium \(V\). We make use of this modification to ensure that the equilibrium point \(v(x_{c})=V\) is unique for \(u_{i}(t)=0\). A globally unstable configuration is set by choosing \(R=0.25\), \(P=0.05\), \(Q=0.0005\), \(V=0.2\) and \(L=60\) with periodic boundary conditions. The control term \(\sum_{i=0}^{m}B_{i}(x_{c})u_{i}\) comes similarly to the approach proposed by Fabbiane _et al._ (2014), where \(m\) is the number of control inputs, \(B_{i}(x_{c})\) is a window function for each actuator and \(u_{i}\) represents the control signals.

The equation is discretised in the present simulation with \(\Delta x_{c}=1\) and \(\Delta t=0.025\). Sampling
for control is made each 400 timesteps, thus giving a control timestep \(\Delta t_{c}=10\). Three actuators are employed by using Gaussian window functions centered at \(x_{c}=10\), \(x_{c}=30\) and \(x_{c}=50\). Sampling of \(v\) is conducted at each of the 60 discrete grid points. Figure 4 shows a sensor/actuation schematic. In this test case, we conduct equilibrium estimation, sensor placement, control and stability analysis.

### Periodic 2D channel flow

Another problem we consider is streamwise-periodic channel flow. We employ numerical simulations using an open source Nek5000 CFD code previously applied for RL control (Li & Zhang, 2022) to solve the incompressible Navier-Stokes equations for a Reynolds number \(\mathrm{Re}=8000\), based on the half channel height \(h\) and the velocity at the channel center line \(U\). The channel length is set to \(L=36h\) so that two unstable modes appear, at streamwise wavelengths L/5 and L/6. A total of \(m=8\) control inputs \(u_{1},\ldots,u_{8}\) modulate the velocities of 8 pairs of actuators. Each pair works in opposition to ensure a zero-net mass flux, being weighted by a parabolic window function as depicted in figure 5.

The distribution of candidate sensors is presented in figure 6. The chosen locations are concentrated near the wall to capture velocity fluctuations from hydrodynamic instabilities. A vertical line of 50 sensors is also employed for comparison with stability analysis computed from the solution of the Orr-Sommerfeld equation. In total, 332 candidate sensor positions are chosen, where horizontal and vertical velocities \(u(x_{c},y_{c})\) and \(v(x_{c},y_{c})\) are sampled, totalling 664 signals. For this channel problem, we are interested in equilibrium estimation, sensor placement, control, and stability analysis.

### Confined cylinder flow

Finally, the study of a 2D confined flow past a cylinder is presented. The initial configuration of sensors and actuators is the same as that proposed by Rabault _et al._ (2019). For this task, a Nek5000 setup (Li & Zhang, 2022) is employed considering a Reynolds number \(\mathrm{Re}=150\) based on the cylinder diameter \(D\) and maximum velocity \(U\) of the inflow profile. The initial configuration of sensors is equivalent to the one presented by Rabault _et al._ (2019), depicted in

Figure 4: Actuation and sensor setup for controlling the modified KS equation. The coloured dots represent the control input profiles \(B_{i}\) at each of the 60 positions measured initially (candidate sensors).

Figure 5: Actuation setup for controlling the 2D channel flow. Each pair of actuators works in opposition, ensuring zero-net mass flux.
figure 7, where 153 locations are chosen, totalling 306 velocity measurements for \(u\) and \(v\). No-slip wall boundary conditions are applied at the top and bottom limits of the domain.

The actuation scheme is presented in figure 8. A single control input modulates minijets in opposition, thus providing zero-net mass flux. For this problem, we are interested in finding an equilibrium point (only to estimate a control setpoint), sensor placement and flow control.

## 5 Results

### Lorenz equations

To assess the performance of the proposed techniques, results for the low order Lorenz system are presented. Starting with the problem of estimating an equilibrium point for this system, figure 9 presents the values of each state at each iteration. Two plots compare this computation using both the backpropagation-based linearization of the neural network (as described in section 2.3) and the alternative version employing the modified DMDc method proposed in section 2.5. In both cases, a convergence towards equilibrium at \(x=y=z=0\) is seen. It is clear that, for the initial iterations, better approximations are obtained using backpropagation (figure 9(a)). The fact that DMD based techniques rely on data sampled from approximately linear systems compromise their performance when the plant is strongly nonlinear. As discussed by Deda _et al._ (2023),

Figure 8: Actuation scheme applied to the cylinder flow. Blowing/suction devices in opposition are modulated by a single control input.

Figure 6: Initial candidate sensor positions for the channel flow. The dimensions are distorted for better visualization.

Figure 7: Initial candidate sensor positions for the confined cylinder flow. The box represents the spatial domain used in the simulation.
the linearization of a nonlinear model at the point of interest can provide better results when nonlinearities cannot be neglected. On the other hand, for the later iterations, the DMD based approach is able to provide more accurate results, since the closed-loop control is able to maintain states near equilibrium, i.e., operating almost as a linear system (figure 9(b)).

The stabilization of the controlled Lorenz system is depicted in figure 10, where signals are presented against time for iterations 1 and 4. For the former case, depicted in figure 10(a), it is still not possible to obtain a stabilizing controller. This might be due to poor estimation of the dynamics near the equilibrium point obtained for such iteration. On the other hand, at the iteration 4 shown in figure 10(b), the trained NNC is able to fully stabilise the plant, bringing the states close to \(x=y=z=0\). As seen in figure 9(a), the equilibrium point for iteration 4 is estimated at \(x\approx 0.14\), \(y\approx-0.01\) and \(z\approx 0.53\), which is close enough to allow for the NNC to stabilise the plant (though subsequent iterations have much more accurate estimates).

Finally, a modal stability analysis is conducted, where poles are computed for the controlled Lorenz system. The modified DMDc approach is employed to data obtained from a single iteration of the _sweep_ mode. Similarly to what is observed in figure 9(b), later iterations tend to provide better approximations of equilibria. Figure 11 presents results for iterations 4 and 9. For comparison, the Lorenz system equations are linearised analytically, and the obtained poles \(s_{i}\) are discretised so that \(z_{i}=e^{s_{i}\Delta t}\). These poles are presented as "ground truth" in figure 11. As also shown in figure 3, at iteration 4, the perturbed controlled system is driven to nonlinear orbits. At iteration 9, sampled states are concentrated much closer to equilibrium, thus providing better approximations through linear regression.

### Modified Kuramoto-Sivashinsky equation

For the modified KS equation, two different training approaches are presented. They only differ in terms of application of the sparsity layer. One of the cases does not include the sparsity layer on the neural network, (i.e., the full set of sensors is employed), whereas the other case includes a sparsity layer to reduce the number of sensors required. In this second case, the total number of

Figure 9: Equilibrium point estimation along iterations. The estimation through backpropagation is shown in (a), where more accurate estimates are made at first iterations. The fixed estimation through the DMDc variant is shown in (b), where better approximations are found in the last iterations. In this case, the red dots are hidden behind the blue ones.
sensors used to infer a model for the system dynamics, as well as to conduct closed-loop control, is reduced from 60 to 9, as represented in figure 12. In both cases, the trained process is conducted along 10 iterations of _sweep_, _release_ and _control_. We observe in figure 12 that the sparse sensors are selected to be approximately evenly distributed throughout the domain, with three sensors located across each of the three Gaussian functions defining the actuators, and with one sensor at or slightly downstream of each of the actuator peaks.

Figure 13 shows the calculation of the equilibrium points along \(x_{c}\) found in each case using Newton's method and the modified DMDc. It is possible to observe in figure 13 (a) that the modified DMDc is able to enhance the accuracy of the equilibrium point when the full set of sensors is employed. However, as a trade-off for reducing the number of sensors, figure 13 (b) shows that the estimate of the equilibrium is worsened with sparse sensors. Despite this, the

Figure 11: Poles of the discrete-time Lorenz system obtained by linearization of the controlled system compared to ground truth positions. (a) At iteration 4, the approximation is still not satisfactory, which is explained by the nonlinear dynamics present in the data employed. (b) At iteration 9, the estimated poles are found considerably closer to the ground truth solution.

Figure 10: _Sweep_, _release_ and _control_ modes (_s.m._, _r.m._ and _c.m._, respectively) for iterations (a) 1 and (b) 4. The trained NNC is unable to stabilise the system in iteration 1, in contrast to the trained control system in iteration 4, where stabilization is achieved. This can be seen by comparing both control modes (_c.m._).
trained NNC is still able to stabilise the system, as exemplified in figure 14, where the waves are controlled until convergence.

For the modal stability analysis of the modified KS equation, the computed poles are shown in figure 15. The ground truth poles are obtained by linearizing equation 4.5 with a discretization performed using fourth-order centered finite difference schemes for all spatial derivatives. The poles are found by computing the eigenvalues \(s_{i}\) of the matrix built from the approximation (with periodic boundary conditions) and transformed to their discrete version \(z_{i}=e^{s_{i}\Delta t}\). The map obtained with the proposed technique is fairly accurate, with good preservation of frequency and growth rates of the eigenvalues closest to the unit circle. Since the sparse case comprises only 9 sensors, only 9 eigenvalues are found.

### Periodic 2D channel flow

The 2D channel flow case is studied with dynamic sensor placement. The final distribution of probes after 19 iterations of training (with _sweep_, _release_, and _control_ modes) is presented in figure 16. It is possible to notice the preference for horizontal velocity probes (\(u\)) rather than vertical (\(v\)) ones. The sensor placement is able to achieve a reduction from 664 to 181 signals, 148 and 33 correspond to \(u\) and \(v\) velocity probes, respectively. The spatial distribution of these chosen sensors shows that almost all of the near-wall \(u\) sensors are included, possibly highlighting the

Figure 12: Sparse configuration obtained in training for the case with sensor placement, where 9 sensors are kept for the KS equation.

Figure 13: Estimated equilibrium points of the KS equation found with (a) full and (b) sparse set of sensors. The true equilibrium is \(v(x_{c})=0.2\).
importance of near-wall information for estimating and controlling (with wall-based actuation) the system dynamics.

With feedback of the signals probed at the locations obtained through training, the NNC is able to completely stabilise the flow. A comparison of the uncontrolled and controlled flows is presented in figure 17, where \(v\)-velocity contours are shown. As expected, the actuation and the flow field do not converge to the exact equilibrium, since its estimate is not perfect, especially considering that, for control, we employ the approximation obtained through Newton's method (see section 2.5). The controlled flow is presented in two ways: 1) using the same contour range of the uncontrolled flow, and 2) with a more saturated contour range that allows one to notice some

Figure 14: Space-time map showing the growth (decay) of disturbances in the KS equation before (after) turning the controller for the (a) full and (b) sparse set of sensors. The black vertical line highlights the instant when the control is turned on.

Figure 15: Poles of the discrete-time modified KS system obtained by linearization of the controlled system compared to the ground truth solution. The (a) full and (b) sparse sensor approaches are compared. Most of the ground truth poles are concentrated close to the origin, making them harder to visualise.
small residual velocity fluctuations. When employing the same scale as that of the uncontrolled case, the bias becomes nearly imperceptible. For the uncontrolled case, the snapshot represents a frame captured from an unsteady flow, where the flow structures are transported from left to right, while the controlled case is simply a steady state.

The present flow includes two unstable modes. These are computed through the eigendecomposition of the associated Orr-Sommerfeld operator, which is discretised using a pseudo-spectral method with Chebyshev polynomials utilizing the toolbox by Weideman & Reddy (2000). The wavelengths of the unstable modes obtained by the Orr-Sommerfeld equation correspond to 1/6 (mode 1) and 1/5 (mode 2) of the channel length. This is in agreement with the wavelengths of the dominant structures observed in simulations of the uncontrolled nonlinear system, as shown in figure 17(a).

By applying the modified DMDc approach to the data obtained near equilibrium, two unstable pairs of modes are found, corresponding to the two unstable modes of the true linearised system. For this case, the discrete poles \(z_{i}\) found with the modified DMDc are converted to a continuous version \(s_{i}=\ln{(z_{i})}/\Delta t\), so that one can directly infer the corresponding growth rates and frequencies. A comparison between the true and estimated unstable eigenvalues is shown in Table 1, where close agreement is observed. In these results, the imaginary part corresponds to the oscillation frequency while the real component provides the growth rate. The higher relative error found for the growth rate when compared to the frequency can be explained by the transient time scales being too slow compared to the discretization time \(\Delta t\). Since \(\Delta t\) needs to be small enough to capture the signals without significant aliasing, and the timescale of the oscillations are considerably faster than the growth rate, the discrete time horizon to capture the long-term growth becomes too large. The present approach is also able to provide the associated eigenvectors. In figure 18, the absolute values of the \(u\) and \(v\)-velocity eigenvectors related to

Figure 16: Probe locations employed with sparse sensing, where the number of probes is reduced from a total of 664 to 181, from which 148 are used for the horizontal velocity component, while 33 are used for the vertical one. The grey dots are deactivated by the sparsity layer.

Figure 17: Comparison of uncontrolled and controlled flows through visualization of \(v\)-velocity contours (see Movie 1). For the controlled case, a more saturated contour range is also shown in the bottom plot, and one can see that the actuation does not converge to zero due to an imperfect estimation of the equilibrium.
each pair of eigenvalues are plotted. Here, these eigenvectors are mapped to the vertical line of sensors for comparing with the ground truth values. Despite displaying a slight asymmetry for the proposed method, which could likely be improved by adding more sensor points, the comparisons are in good agreement.

### Confined cylinder flow

As was the case in section 5.3, this system is modelled and controlled while also incorporating a sparsity layer in the NNSM to select sensors. In doing so, the NNSM training process for the confined cylinder flow selects a total of 47 sensors from the candidate set of 306 sensors: 25 and 18 points for \(u\) and \(v\) velocity components, respectively. While the number of sensors is significantly reduced, as presented in figure 19, the NNC approach with iterative training is able to successfully stabilise the this flow. To provide an overview of the training process, the lift coefficient over the cylinder is presented in figure 20 for the first 3 iterations. In the first iteration, the NNSM and the NNC are trained using only open-loop data. For this iteration, the lift fluctuations are only slightly reduced when flow control is turned on. For the following iterations, the models are trained using both open and closed-loop data, which improves the results, as more data is available near the equilibrium point, where the system is approximately linear.

While the stabilization is not completely achieved by the end of the second iteration, a reduction in the main oscillations can be observed. By the third iteration, complete stabilization is achieved through control.

Important features observed in the controlled flow are shown in figure 21, which shows the convergence history of the drag and lift coefficients, as well as of the control input, during the control stage of the third iteration. The present flow control approach is able to suppress both the drag and lift coefficient fluctuations. At the same time, the mean drag is also reduced. The complete

\begin{table}
\begin{tabular}{c c c} \hline \hline  & Mode 1 & Mode 2 \\ Ground truth & 1.8601e-3 \(\pm\)/2.6403e-1 & 6.8185e-4 \(\pm\)/2.0206e-1 \\ Modified DMDc & 2.0545e-3 \(\pm\)/2.6396e-1 & 7.4956e-4 \(\pm\)/2.0214e-1 \\ \hline \hline \end{tabular}
\end{table}
Table 1: Eigenvalues for unstable modes 1 and 2 for 2D channel flow. Colours correspond to modes plotted in figure 18.

Figure 18: Unstable modes for the 2D channel flow. Ground truth and estimated results are compared in terms of absolute values for the eigenvectors mapped to the vertical line of sensors.
stabilization also ensures that only a minimum effort is required to keep the system operating with small oscillations and drag losses. Residual efforts are due to imperfect estimation of equilibrium, and may also be required for compensating eventual perturbations. Figure 22 shows \(u-\)velocity contours for the uncontrolled and controlled cases. In the former, vortex shedding develops along the wake, inside the plane channel. On the other hand, a steady solution is obtained for the controlled case.

To contextualize these observations, we now discuss prior control results for this confined cylinder configuration. The most direct comparison can be made with Li and Zhang (2022), who also attempted to stabilize the flow. They proposed an approach for the stabilization problem using reinforcement learning to train neural networks capable of minimizing the vortex shedding energy, which is translated to a reward function that takes the velocity fluctuations at selected probe locations. Such calculations could be done by either subtracting the mean flow or the equilibrium flow found through selective frequency damping (SFD) from the measured velocity probes. By using the former approach, convergence was not observed, whereas using the equilibrium states from SFD enabled stabilization. The results, however, required manual positioning of probes through a heuristic approach informed by considering the wavemaker region obtained from considering the region of overlap between the leading direct and adjoint eigenmodes. Furthermore, the prior computation of the equilibrium through SFD is a dedicated step that requires additional

Figure 19: Demonstration of the application of the sparsity layer for the confined cylinder flow. The total number of measurements is reduced from 306 to 47 (25 and 18 probes for the horizontal and vertical velocity components, respectively). The grey dots are deactivated by the present L1 regularization.

Figure 20: Evolution of lift coefficient for 3 iterations of the _sweep_ (blue), _release_ (yellow), and _control_ (red) stages for the confined cylinder flow.
simulation-intrusive computations. By contrast, our approach enables automatic equilibrium (and eigenmode) computations, which do not require the use of modified numerical solvers. Li & Zhang (2022) found that successful control of the confined cylinder flow resulted in a small, nonzero mean control signal (i.e., control jet flowrate), which we also observe here. In our case, this is on the order of \(\tilde{u}\) = -6e-4, which is too small to observe on the right subplot of figure 21.

It is more difficult to make a direct comparison with other prior control studies of this system, primarily due to the different control objectives considered. The current work aims for stabilization, whereas other prior investigations considering this confined cylinder flow (Rabault _et al._, 2019; Varela _et al._, 2022) focused on drag minimization. As we holistically approach the control problem proposing methods that also provide estimates of the equilibrium, the stabilization task is directly enabled through NNC training. For this confined cylinder flow, a drag reduction is seen through stabilization, the lift and drag oscillations are suppressed, and the control input required to keep the controlled system at the unstable equilibrium approaches zero. While we find a drag reduction of 2.6% for the fully-stabilized flow, larger drag reduction--5.7% at \(\mathrm{Re}=100\) and 21.6% at \(\mathrm{Re}=200\)--are reported in Varela _et al._ (2022). This suggests that the optimal solution for drag reduction is not the flow at its equilibrium. This can be verified, for example, in the strong overshoot seen in \(C_{D}\) (Fig. 21), where the maximum drag reduction is reached before the flow is stabilized. In fact, in iteration 1, the control system produces an average drag reduction of 3.4%, lower than the stabilized flow at iteration 3 of 2.6%.

Figure 21: History of drag and lift coefficients, and control effort for the third iteration in the control mode. Gray lines show uncontrolled solution.

Figure 22: Comparison of uncontrolled and controlled \(u\)-velocity fields (see Movie 2). The former is represented by a snapshot taken at the flow natural limit cycle. For the controlled case, the snapshot is captured after stabilization is achieved.
## 6 Conclusions

In the present work, neural networks are applied to perform as surrogate models of nonlinear dynamical systems, and also as controllers for stabilizing such systems. The proposed approach is first tested with the Lorenz set of ordinary differential equations, as well as with a modified version of the Kuramoto-Sivashinsky partial differential equation. Then, flow control is demonstrated for a periodic 2D channel flow, for which a data-driven linear stability analysis is also conducted and, finally, for a more complex case of a confined cylinder flow. The trained neural networks provide adequate surrogate models that achieve all the present results with a single hidden layer, except for the confined cylinder case, which employs two layers. In the same fashion, the neural networks obtained to control the investigated systems are also single-layer models with only a few neurons. Despite the simple individual neural network architectures, through a recurrent training strategy, the present framework is able to obtain effective controllers. Results also support the efficacy of neural networks for learning complex system dynamics, as well as to obtain important flow features, such as equilibrium points, and leading eigenmodes of the linearised system about these fixed points. Although these data-driven models are trained as black boxes, relevant information can be extracted by using backpropagation.

In this work, several different nonlinear unstable systems are investigated, and their dynamics do not naturally converge to equilibrium points. This behaviour hinders data sampling that represents the plant dynamics close to such points, where linear approximations should be able to inform the main eigenvalues and linear modes. As a solution for this problem, we propose the application of NNC with an iterative training of the models. As well as achieving stablization, this new method proves to be an efficient approach for equilibrium computation, sensor placement, and modal stability analysis. As demonstrated by results, the iterative approach is key for bringing the system closer to a progressively better estimation of the equilibrium point. In all the cases presented, the first iteration is unable to provide accurate results, specially for equilibrium computation and modal analysis. Even for the control task, stabilization is not achieved, for example, in the cylinder flow case. The improvement achieved with the iterative process puts the presented NN approach as an efficient and accurate candidate for flow stabilization and data-driven flow analysis, with the possibility of enabling limited sensor allocation.

The first plant studied with the proposed methodology is the Lorenz system. With the chosen numerical parameters, its dynamics results in a chaotic attractor featuring 3 different equilibrium points. While not explicitly shown, it was found that the initial guess for the equilibrium point has an important role in the iterative process, which means in general that if initial guesses are too far from an equilibrium of interest, the algorithm may converge to another point, which might not be desirable. In the case presented in the results section, the initial choice at \(x=y=z=0.1\) converges near the origin after a few iterations. Results from this first test indicate that the iterative training of the NNSM and NNC successfully overcomes the lack of data sampled near the equilibrium. As the control strategy improves and the perturbations are reduced in amplitude, datasets containing more measurements near equilibrium are built, improving the quality of linear approximations and the control design aiming stabilization. Linear regression through DMD-based algorithms is also enabled for a considerably broader range of plants, taking advantage of data produced by controlled nonlinear systems that are brought closer to linear operation through feedback control. The modification to the DMDc method used in this work also helps extend the analyses to systems whose equilibrium points are unknown or poorly estimated, allowing for a model fitting that also enhances estimates of equilibrium bias.

The present methodology is next applied to control a more complex problem, a modified version of the Kuramoto-Sivashinsky equation. This modified equation is proposed in order to force a single equilibrium point instead of a continuous space of possibilities. This test case allows the study of the effectiveness of sensor placement through L1 regularization of the inputs. Again,
equilibrium estimation and feedback control obtain satisfactory results, even with a reduction in the number of sensors from 60 to 9. Although this reduction slightly worsens the estimation of equilibrium and, therefore, adds further errors to control convergence, tuning the L1 regularization through the choice of a single hyperparameter leads to a configuration with an intermediate number of sensors at reduced accuracy cost. Even so, the application of Newton's method with linearization of the NNSM through backpropagation is able to provide an equilibrium point estimation that allows for stabilization and, therefore, to obtain the approximately linear dataset subsequently used to identify the dominant eigenmodes of the linearized system.

Application of the NNC is then applied to the Navier-Stokes equations. As a first case, stabilization of a 2D streamwise-periodic channel flow is sought, where the capability of conducting a linear stability analysis is also explored for a flow configuration with two unstable modes. The L1 regularization is applied in this case and, from a total of 664 measured signals, training is able to reduce sensing to 181 probes. Stabilization is achieved with 8 pairs of actuators in opposition through NNC, which enables sampling of data close to equilibrium, i.e., of signals described by approximately linear dynamics. The data is used for performing a linear regression through a modified DMDc. The unstable eigenvalues and eigenmodes computed by this technique show good agreement with those obtained directly from the solution of the linear Orr-Sommerfeld operator. The different magnitudes of the growth rate and frequency of the unstable modes result in different levels of accuracy for their estimation, with the growth rates being close to 0. In this sense, the frequencies present lower relative errors than the growth rates due to the large discrete time intervals associated with the latter for the chosen time step. The natural growth rate corresponds to a very small component of the measured signal compared to the effect of the uncontrolled frequency.

Finally, control of a confined cylinder flow is also presented. The setup chosen with actuators in opposition is the same proposed in different studies found in literature (Rabault _et al._, 2019; Li & Zhang, 2022). Here, we assess the ability of the neural network model to stabilise a complex system through the iterative training framework. Even for this more complex case, the NNC is built from a single layer containing only 8 nodes. Results show that flow control suppresses the lift and drag oscillations from vortex shedding and also reduces the mean drag by bringing the flow to its equilibrium point. Not only is the system successfully controlled, but additionally the low complexity of the NNC suggests that there is considerable room for applying the technique to considerably more complex fluid systems without prohibitive costs. For this confined cylinder configuration, prior work (Li & Zhang, 2022) found that RL control performance could be improved through equilibrium computation and linear stability analysis in the vicinity of this equilibrium. Here, we achieve similarly successful control results without needing to perform these auxiliary computations explicitly, but where the results of such computations are available as a byproduct of our modeling and control methodology.

In each of the examples featuring partial differential equations across a spatial domain, we have utilized a sparsity layer in the NNSM to reduce the number of sensor measurements (and thus model inputs). While the resulting sensor locations are chosen to optimise a given cost function, they are not necessarily unique solutions, in the sense that alternative sensor locations may be obtained for different realizations (with different random input signals and randomly initialized NN weights) of the training procedures. This is clear from observing that the identified sensors do not follow the same symmetry properties as the geometries of the problem. That being said, the chosen sensor locations do reveal some aspects of the underlying physics. For the modified KS equation, the sensors are distributed throughout the domain, in accordance to the translation invariance of the problem (aside from the actuator locations). For the 2D channel flow, the sensor locations indicate the importance of near-wall streamwise velocity measurements, with a sparser set of measurements away from the wall. The asymmetry in the chosen locations of the wall-normal velocity probe locations is perhaps due to the fact that the instabilities in this system
consist of mode shapes where the vertical velocity component extends across the whole domain, so taking measurements at both the upper and lower walls may be unnecessary. In the confined cylinder flow, we find that sensors both close to the cylinder, and downstream in the wake are chosen, indicating the importance of measurements in both regions for suppressing the natural vortex shedding behavior. While beyond the scope of the present work, it could be interesting to compare the identified sensor locations with those obtained from alternative methods (Manohar et al., 2018; Sashittal and Bodony, 2021; Williams et al., 2022; Graff et al., 2023).

The modeling and control framework utilised here is similar in principle to that used in classical linear control theory, with separate models for the plant and controller connected in feedback. While nonlinear control is substantially more complex, here we show that relatively simple neural networks can be used to develop effective nonlinear models for both the plant and controller, when trained with a method that promotes the generation of large quantities of near-equilibrium data. Future work could also utilise linearization of the plant and/or controller models for linear control design, with better authority over the plant behavior through well studied control design techniques that can ensure robustness, optimality or specific pole placement.

While it was demonstrated that the proposed methodology could be used to identify modal linear amplification mechanisms, future work could extend these methods to study nonmodal amplification, such as transient growth and resolvent analysis. Data-driven implementations of such analyses have been implemented previously (Herrmann et al., 2021), though for non-normal systems it may be more difficult to obtain accurate results without adjoint operators/data, which may require further modifications to the methodology presented here.

We lastly emphasise that the methods and results are obtained from non-intrusive methods in the sense that only the nonlinear system (with appropriate control inputs) has been used. This means that the methodology could in principle be applied (for both control and flow physics analysis) in an experimental setting. This, as well as the application to systems presenting higher degree of complexity, remains the subject of future work.

## Acknowledgments

The authors acknowledge the financial support received from Fundacao de Amparo a Pesquisa do Estado de Sao Paulo, FAPESP, under Grants No. 2013/08293-7 and 2021/06448-0. The first author is supported by FAPESP PhD scholarships No. 2019/19179-7 and 2022/00469-8, which are also acknowledged. STMD acknowledges support from NSF Award 2238770. The computational resources used in this work were provided by CENAPAD-SP (Project 551), and by LNCC via the SDumont cluster (Project SimTurb).

## Declaration of Interests

The authors report no conflict of interest.

## Appendix A Hyperparameters

In this appendix, the chosen values for the hyperparameters in each case are presented. They are organised in table 2. Regarding the iterative training, the "number of iterations" correspond to how many times the proposed training steps are conducted. It is expected that at each iteration the perturbed plant controlled by NNC tend to produce more data near equilibrium. The "number of steps in _sweep_", "number of steps in _release_" and "number of steps in _control_" correspond to the number of measurements sampled in each mode at each iteration.

The open loop control signal is parameterised by a "control input saturation value", which defines the possible uniform interval of random values it can assume. The staircase signal has
a step width given by the "number of time steps for each stair step (open loop signal)". Also, the "maximum amplitude decay at each iteration" parameter makes sure the open-loop signal is reduced in amplitude at each training iteration.

From each mode (_sweep_, _release_ and _control_), data is sampled for training, but the dataset size is limited to a "maximum number of points sampled for training NNSM". This ensures that the training cost does not increase beyond a certain point. A "NNSM learning rate" is defined for the ADAM optimization algorithm, and the model is updated through a "number of epochs at each training iteration for the NNSM". The "number of neurons at NNSM hidden layers" hyperparameter defines the number of layers, which is given by the number of elements in the array of numbers, and the number of neurons in each layer, which is the value at each position. Note that, except for the cylinder case in table 2, all NNSM models are trained with a single layer. For assessment of possible overfitting, an "NNSM training set ratio" as a fraction of sampled data is also chosen, which allows for the comparison of losses between a training set and a testing set of samples. The "hidden layers L2 regularization" parameter is the \(r_{2}\) variable defined in this work. "Use sparsity layer" is a Boolean that determines if sparse sensor placement is used. If so, we can set the "sparsity layer L1 regularization" parameter (\(r_{1}\)) and the "sparsity layer truncation tolerance" that determines the absolute value below which the input layer weights are truncated to zero. The input layer is organised according to the "number of control inputs" and the the "number of states" measured.

Similarly a "maximum number of points sampled for training NNC" is also defined to avoid excessive complexity, as well as the "NNC learning rate" for ADAM and the "number of epochs at each training step for NNC". The "training finite horizon length" in number of steps (\(n_{h}\)) is also presented. For the NNC, an "NNC training set ratio" as a fraction of sampled data is also chosen. The "number of neurons at NNC hidden layers" is always the same in all cases, which consists of a single layer containing 8 neurons. "NNC control saturation" stands for the limits of control effort that NNC can provide. It is implemented as an output layer consisting of weighted sigmoid functions. The "control input weight in loss function" corresponds to \(\mathbf{w_{u}}\).

Finally, regarding equilibrium computation, the "initial guess for equilibrium point" is presented for each case, corresponding to the first guess used in the first training iteration. At each of these iterations, the estimation of equilibrium is updated along a "Number of steps for the Newton method".
