Efficient computations for linear feedback control problems for target velocity matching of Navier-Stokes flows via POD and LSTM-ROM

1.   Introduction

The problem of fluid flow control has long been the focus of fluid mechanics research and is still an important research topic for many researchers. One of the most important fluid control problems is target velocity matching or tracking issues for fluid flow. The issue of target velocity matching or tracking problems for fluid flows has been a very active area of research e.g. [1,15,17,18,19,21,22,29,31] as well as others. The computation of optimally controlled flows based on solving an optimality system is expensive in terms of CPU time and memory space. This is due to the fact that such an approach involves a coupled system of state and adjoint equations with initial and final conditions. The system has to be solved on the entire space-time domain and cannot be solved by marching in time. Therefore, a feedback control design and/or reduced order modeling (ROM) is required for efficient and real time computing.

Feedback control is a less expensive solution technique that gives qualitatively similar results to optimal control. Feedback control can effect good velocity tracking and at the same time the solution can be obtained step-by-step in time at the cost of a single flow solve. Mathematical theories and approximation techniques for optimal control problems for the Navier-Stokes equations have been developed in various areas of fluid dynamics; see, e.g., [1,17,21,22]. Some numerical methods for solving control problems for unsteady flows have been proposed and tested; see, e.g., [17,18,19,21,22] and references therein. Many of these deal with distributed controls and take different approaches to limiting the size of the control. Feedback controls for the Navier-Stokes equations have also been the subject of previous, mostly computational, studies; see, e.g., [11,16] and references therein. In this article, we consider desired states ${\boldsymbol U}$ in ${\bf L}^2$ which are not necessarily divergence-free, which is a feature not seen in previous work.

Proper orthogonal decomposition (POD) is a common technique for extracting the dominant mode that contributes the most to the energy of the entire system [6]. POD combined with Galerkin projection (GP) have been used for many years to formulate ROMs for dynamic systems [5,8,9,32,24]. In such ROMs, the full-order set of equations is projected onto a reduced space, resulting in a dynamic system (modal coefficients) of much lower order than the full order model (FOM). However, in many situations there is a discrepancy between the governing equation and the observed system. This can be due to an approximation of the underlying phenomenon, incorrect parameterization, or insufficient information about the source terms contained in the system. This is especially evident in fluid flow systems where many complex phenomena and source terms interact in different ways. A common reduced order modeling (ROM) development procedure may be described by the following tasks:

1. Reduced space and basis identification.

2. Nonlinear dynamical system evolution in the reduced space.

3. Reconstruction in full-order space for assessments.

Machine learning (ML) tools have had considerable success in the fluid mechanics community, identifying basic structures and mimicking dynamics [3,7,25,26]. However, modeling with ML, especially deep learning, has faced strict opposition from both academia and industry alike because it can produce non-physical results due to its black box nature and its lack of interpretability and generalization [10,23]. In the course of conducting this study, we also experienced these points seriously. Even if the input and output data used to train the ML algorithm is physically accurate, the quantity interpolated with the ML approach can deviate significantly from a physically accurate solution. However, the fluid mechanics problem solving using deep learning is considered a promising research in the future, and we think it should be continued. A perspective on machine learning for advancing fluid dynamics can be found in a recent review articles [7] and references therein.

The plan of the rest of the paper is as follows. In Section 2, we define and discuss the linear feedback control. We derive a time-space discretized version of the feedback control of Navier-Stokes equations. In Section 3, we explain the POD reduced basis and Galerkin Projection ROM. In Section 4, we introduce a closure model using deep neural networks. Finally, in section 5 some numerical results will be given.

2.   A feedback control problem for Navier-Stokes equations

Let ${\Omega}$ be a bounded open set. We shall use the standard function spaces and their norms; for details see [2,13]. For any nonnegative integer $m$, we define the Sobolev spaces $H^m ({\Omega})$ by

where $D^{\alpha} u$ denotes the weak (or distributional) partial derivative and $\alpha$ is a multi-index, $|\alpha| = \sum_i \alpha_i$. Note that $H^0({\Omega}) = L^2({\Omega})$. We equip $H^m({\Omega})$ with the norm

The usual inner product associated with $H^m({\Omega})$ will be denoted by $(\cdot , \cdot)_m$. Let $H^m_0({\Omega})$ denote the closure of $C^{\infty}_0({\Omega})$ under the norm $|| \cdot ||_m$ and $H^{-m}_0({\Omega})$ be the dual space of $H^m_0({\Omega})$.

For vector-valued functions, we define the Sobolev space ${\bf H}^m({\Omega})$ (in all cases, boldface indicates vector-valued) by

where ${\boldsymbol u} = (u_1,u_2)$ and its associated norm by

We also define particular subspace

We introduce the solenoidal spaces

The spaces ${\bf V}$ and ${\bf W}$ are closure of ${\boldsymbol {{\mathcal V}}}$ in ${\bf H}^1_0({\Omega})$ and ${\bf L}^2 ({\Omega})$, respectively. All subspaces are equipped with the norms inherited from the corresponding underlying spaces.

In order to define a weak form of the Navier-Stokes equations, we introduce the continuous bilinear form

where $D({\boldsymbol u}) = (1/2) (\nabla {\boldsymbol u}+ \nabla {\boldsymbol u}^T)$ and the trilinear form

A weak formulation of the Navier-Stokes equations is given by given ${\boldsymbol g} \in {\bf H}^{-1}({\Omega})$, find $({\boldsymbol u}, p) \in {\bf H}^1_0({\Omega}) \times L^2_0({\Omega})$ satisfying

with initial velocity ${\boldsymbol u}_0$, where $\langle \cdot, \cdot \rangle$ denotes the duality paring between ${\bf H}^1_0({\Omega})$ and ${\bf H}^{-1}_0({\Omega})$.

We consider the target velocity fields ${{\boldsymbol U}} \in {\mathcal U}_{ad}$, where ${\mathcal U}_{ad}$ is defined by

We will consider the target velocity matching problem in which we want to minimize

To do this, we define the following control and cost functional. Let ${\boldsymbol f} \in L^2((0,T) ; {\bf L}^2({\Omega}))$ denote the distributed control. Given $T$, we define the functional

The minimization of the $\int_0^T \int_{{\Omega}} | {\boldsymbol u} - {\boldsymbol U} |^2 d{\Omega} dt$ is the real goal of velocity matching problem, which is to keep the solution ${\boldsymbol u}$ to ${\boldsymbol U}$ over $(0,T) \times {\Omega}$. The $\int_0^T \int_{{\Omega}} \frac{\beta}{2}| {\boldsymbol f} |^2 d{\Omega} dt$ term is introduced in order to bound the control function and to prove the existence of the optimal control (regularization term). The $\int_{{\Omega}} | {\boldsymbol u}(T) - {\boldsymbol U}(T) |^2 d{\Omega}$ term is necessary in order to keep the solution ${\boldsymbol u}$ to ${\boldsymbol U}$ near the time $T$. An optimal control problem is defined by

seek $({\boldsymbol u}, {\boldsymbol f}) \in L^2((0,T);{\bf H}^1_0({\Omega})) \times L^2((0,T);{\bf L}^2({\Omega}))$ such that the cost functional (6) is minimized subject to constraints

with initial velocity ${\boldsymbol u}_0 \in {\bf H}^1_0({\Omega})$.

Using the Lagrange multipliers method, one can obtain the optimality system: seek ${\boldsymbol u} \in L^2((0,T);{\bf H}^1_0({\Omega})), p \in L^2((0,T);{\bf L}^2_0({\Omega})), {\boldsymbol w} \in L^2((0,T);{\bf H}^1_0({\Omega}))$ and $r \in L^2((0,T);{\bf L}^2_0({\Omega}))$ such that

with initial velocity ${\boldsymbol u}_0 \in {\bf H}^1_0({\Omega})$ and the adjoint equations

with final velocity ${\boldsymbol w}(T,{\boldsymbol x}) = \delta({\boldsymbol u}(T)-{\boldsymbol U}(T))$ and the optimality condition

Several treatments of similar optimal control problems can be found in the literature, most notably in [1,17]. The numerical treatment of the velocity tracking problem is also an outstanding problem and other algorithms have been proposed. For example, a quasi-optimal control has been studied in [21,22].

Remark 1. Thus the optimality system is a system of the coupled nonlinear partial differential equations (8)-(10). We already mentioned that solving an optimality system is very expensive in terms of CPU time and memory space. This is due to the fact that such an approach involves a coupled system of state and adjoint equations with initial and final conditions. The system has to be solved on the entire space-time domain and cannot be solved by marching in time. A gradient algorithm is needed and so the final algorithm is complex, involving multiple flow solutions, and the convergence may be slow.

2.1.   Feedback control

From (10), we see that a control ${\boldsymbol f}$ is related linearly with the adjoint variable ${\boldsymbol w}$ which is a solution of the adjoint equations (9) with ${\boldsymbol f} = - (1/\beta) {\boldsymbol w}$. Also, we see that ${\boldsymbol w}$ is linearly dependent only on the source term of the adjoint equations (9) $\alpha ({\boldsymbol u} - {\boldsymbol U})$ since the adjoint equations are linear equations. So, we see that

Usually, $\beta$ is very small real number to get a small value of the first term $\int_0^T \int_{{\Omega}} | {\boldsymbol u} - {\boldsymbol U} |^2 d{\Omega} dt$ in (6).

with very large value of $\gamma$. Remember that $\gamma$ is closely related to $\beta$. The smaller $\beta$, the greater control ${\boldsymbol f}$. Increasing control is not a problem in numerical simulation, but it is expensive in actual control problems. Therefore, $\gamma$ can be set according to the actual problem.

Now, we consider the feedback control problem

Remark 2. In [19], the control is achieved by means of a linear feedback law relating the body force to the velocity field, i.e., ${\boldsymbol f} = {\boldsymbol F} - \gamma ({\boldsymbol u}- {\boldsymbol U})$ where ${\boldsymbol U}$ is in the set of the admissible target velocity ${\boldsymbol U}_{ad}$ if ${\boldsymbol U}$ is a divergence free vector field in the set $\{{\boldsymbol c}: {\boldsymbol c} \in C((0,T); {\bf H}^2(\Omega) \cap {\bf H}_0^1(\Omega)) : \partial_t {\boldsymbol c} \in C((0,T); {\bf H}^1(\Omega))\}$. The corresponding body force generated by ${\boldsymbol U}$ is defined as

In this article, ${\boldsymbol U} \in L((0,T);{\bf L}^2(\Omega))$ so ${\boldsymbol F}$ can not be defined and we do not include ${\boldsymbol F}$ in the control. We will compare our results with those in [19]. Our computational result is almost the same as that of [19] depending on the value of $\gamma$ (see Figure 3).

Then, a linear feedback control problem of Navier-Stokes equations can be formulated by

with initial velocity ${\boldsymbol u}_0 \in {\bf H}^1_0({\Omega})$

2.2.   The Navier-Stokes system and its finite element approximation

A typical finite element approximation of (14) is defined as follows: we first choose conforming finite element subspaces ${\mathbb{V}}_h \subset {\bf H}^1(\Omega)$ and $S_h \subset L^2(\Omega)$ and then define ${\mathbb{V}}_{h,0} = {\mathbb{V}}_h\cap {\bf H}^1_0(\Omega)$, and $S_{h,0} = S_h\cap L^2_0(\Omega)$. One then seeks ${\boldsymbol u}_h(t,\cdot) \in {\mathbb{V}}_{h,0}$ and $p_h \in S_{h,0}$ such that

where ${\boldsymbol u}_{0,h}({\boldsymbol x}) \in {\mathbb{V}}_{h,0}$ is an approximation, e.g., a projection, of ${\boldsymbol u}_0({\boldsymbol x})$. Discretization is completed by choosing a time-marching method such as the backward Euler scheme. See, e.g., [13] for details about finite element discretizations of the Navier-Stokes system.

Let $\sigma_N = \{ t_n \}_{n = 0}^N$ be a partition of $[0,T]$ into equal intervals $\Delta t = T/N$ with $t_0 = 0$ and $t_N = T$. Then a time discretization of (15) is given by

Since the trilinear term is nonlinear, a linearization of the discrete equations (16) must be introduced in order to solve it by a iteration method. We solve the following linearized equations:

where $\tilde{c} \left({\boldsymbol u};{\boldsymbol v},{\boldsymbol w} \right)$ is the modified trilinear form (see [12]) to preserve the antisymmetry of the trilinear form $c \left({\boldsymbol u};{\boldsymbol v},{\boldsymbol w} \right)$ on the finite element spaces

Given the velocity ${\boldsymbol u}_h^{(n-1)}$ and an integer $m$, one can generate the sequence $\{ {\boldsymbol u}_h^{(n)}(k),$ $p^{(n)}_h \}$ for $k = 1,2,\dots$ by solving the linear problem with $\sigma = 0$ for $k \leq m$ and $\sigma = 1$ for $k>m$. The algorithm terminates when the maximum value of

is less or equal to $\tau$ on ${\Omega}$. The global convergence of the simple iteration method $(\gamma = 0)$, and the rapid convergence of Newton's method $(\gamma = 1)$ are combined in this numerical algorithm. The value of $m$ is set to be two or three, depending on the time step $\Delta t$.

3.   Reduced order method using proper orthogonal decomposition

Let us briefly introduce the POD method following the formulation of Sirovich [30], Rathinam [28] and Burkardt et. al. [8,9]. POD provides a method for finding the best approximating subspace to a given set of data. Originally POD was used as data representation technique. Proper orthogonal decomposition (POD), also known as Karhunen–Lóeve decomposition or principal component analysis, provides a technique for analyzing multidimensional data. This method essentially provides an orthonormal basis for representing the given data in a certain least squares optimal sense.

3.1.   POD reduced-order bases

Given a discrete set of snapshot vectors $W = \{{{\boldsymbol w}}_n \}_{n = 1}^N$ belonging to ${\mathbb{R}}^J$, where $N<J$, we form the $J\times N$ snapshot matrix ${\mathbb A}$ whose columns are the snapshot vectors ${{\boldsymbol w}}_n$:

Let

where ${\mathbb U}$ and ${\mathbb V}$ are $J\times J$ and $N\times N$ orthogonal matrices, respectively, and $\Sigma = {\rm diag} (\sigma_1,\ldots,\sigma_{\widetilde N})$ with $\sigma_1\ge \sigma_2 \ge \cdots \ge \sigma_{\widetilde N}$ be the singular value decomposition of ${\mathbb A}$. Here, ${\widetilde N}$ is the rank of ${\mathbb A}$, i.e., the dimension of the snapshot set $W$, which would be less than $N$ whenever the snapshot set is linearly dependent. It is well known [14] that if

then

so that also

so that $\sigma_i^2$, $i = 1,\ldots,\widetilde N$, are the nonzero eigenvalues of ${\mathbb A}^T{\mathbb A}$ (and also of ${\mathbb A}{\mathbb A}^T$) arranged in nondecreasing order. Note that the matrix ${\mathbb C} = {\mathbb A}^T{\mathbb A}$ is simply the correlation matrix for the set of snapshot vectors $W = \{{{\boldsymbol w}}_n \}_{n = 1}^N$, i.e., we have that ${\mathbb C}_{mn} = {{\boldsymbol w}}_m^T{{\boldsymbol w}}_n$.

In the reduced-order modeling context, given a set of snapshots $W = \{{{\boldsymbol w}}_n \}_{n = 1}^N$ belonging to ${\mathbb{R}}^J$, the POD reduced-basis of dimension $K \le N <J$ is the set $\{ {\boldsymbol {\phi}}_k \}_{k = 1}^K$ of vectors also belonging to ${\mathbb{R}}^J$ consisting of the first $K$ left singular vectors of the snapshot matrix ${\mathbb A}$. Thus, one can determine the POD basis by computing the (partial) singular value decomposition of the $J\times N$ matrix ${\mathbb A}$. Alternately, one can compute the (partial) eigensystem $\{\sigma_k^2,{\boldsymbol {\psi}}_i\}_{i = 1}^K$ of the $N \times N$ correlation matrix ${\mathbb C} = {\mathbb A}^T{\mathbb A}$ and then set ${\boldsymbol {\phi}}_k = \frac1{\sigma_k} {\mathbb A}{\boldsymbol {\psi}}_k$, $k = 1,\ldots,K$.

The $K$-dimensional POD basis has the obvious property of orthonormality. It also has several other important properties which we now mention. Let $\{{\boldsymbol {s}}_k\}_{k = 1}^K$ be an arbitrary set of $K$ orthonormal vectors in ${\mathbb{R}}^J$ and let $\Pi{\boldsymbol w}$ denote the projection of a vector ${\boldsymbol w} \in {\mathbb{R}}^J$ onto the subspace spanned by that set. Further, let

i.e., ${\mathcal E}$ is the sum of the squares of the error between each snapshot vector ${{\boldsymbol w}}_n$ and its projection $\Pi {{\boldsymbol w}}_n$ onto the span of $\{{\boldsymbol {s}}_k\}_{k = 1}^K$. Then, it can be shown that

In fact, often the POD basis corresponding to a set of snapshots $W = \{{{\boldsymbol w}}_n \}_{n = 1}^N$ and then its relation to the singular value decomposition of the matrix ${\mathbb A}$ or to the eigenvalue decomposition of ${\mathbb A}^T{\mathbb A}$ are derived properties. We note that ${\mathcal E}({\boldsymbol {\phi}}_1,\ldots,{\boldsymbol {\phi}}_K)$ is referred to as the "POD energy" or "error in the POD basis". Also, it can be shown that

i.e., the error in the POD basis is simply the sum of the squares of the singular values corresponding to the neglected POD modes.

3.2.   Galerkin Projection-ROM (GP-ROM)

We now show how a POD basis is used to define a reduced-order model for the Navier-Stokes system. For the sake of brevity, we only discuss the case for which the snapshot set is viewed as a set of finite element coefficient vectors; the case for which the snapshot set is a set of finite element functions proceeds in an almost identical manner.

Let $\{\boldsymbol {\phi}\}_{k = 1}^K$ be a $K$-dimensional POD basis corresponding to the snapshot set $\{{{\boldsymbol w}}_n\}_{n = 1}^N$ and let

We then determine ${{\boldsymbol u}}^K_h(t,\cdot) \in {\bf U}_K$ from the discrete problem

For a given reduced space ${\bf U}_K$, GP-ROM finds the approximation of the velocity field spanned by the low-dimensional space,

where $\left\{ a_k(t) \right\}_{k = 1}^K$ are the sought time-varying coefficients. The GP-ROM can be obtained by projecting the FOM (Full order finite element model) onto the POD space:

for $\ell = 1,\ldots,K$. Equivalently, we have the system of nonlinear ordinary differential equations that determine the coefficient functions $\{a_k(t)\}_{k = 1}^K$:

or

where the Gram matrix ${\mathbb G}$, stiffness matrix ${\mathbb K}$, convection tensor ${\mathbb N}$, and solution vector ${\boldsymbol {a}}(t)$ are respectively given by

for $k,\ell,m = 1,\ldots,K$. Note that all of these matrices and tensors are full; however, since $K$ will be chosen small, this does not cause any computational inefficiencies.

4.   Closure model using least squares projection and LSTM

It is easy to see that the computational cost of GP-ROM (24) is $O(K^3)$, which limits the number of modes to be used in the ROM. Several efforts have been devoted to introduce stabilization and closure techniques to account for the effects of truncated modes on ROM's dynamics. In this study, we utilize a LSTM closure model to improve the accuracy and efficiency for our computations.

The model reduction error ${\mathcal E}^{\mbox{ ROM}}(t,{\boldsymbol x})$ is equal to

The first term ${\mathcal E}_{{{{\boldsymbol U}}_{K}}^{\perp}}$ is the projection error that is a result of neglecting the components of ${\boldsymbol u}^{\mbox{ FOM}}(t,{\boldsymbol x})$ that lie in the space orthogonal to ${\bf U}_K$. This error is represented as ${\mathcal E}_{{{{\boldsymbol U}}_{K}}^{\perp}}$ in Fig. 1. The second term ${\mathcal E}_{{{{\boldsymbol U}}_{K}}}$ is the reduced order modeling error that results from solving a different dynamical system than the original. This error vector lies in the subspace ${\bf U}_K$ spanned by the reduced order basis $\{{\boldsymbol {\phi}}_k \}_{i = 1}^K$. This error represents ${\mathcal E}_{{\boldsymbol {U}}_K}$ in Fig. 1. For more details, one can refer the articles [3,4].

The GP-ROM equations (24) with linear and nonlinear operators can be written as:

where ${\mathcal L}$, ${\mathcal I}$ and ${\mathcal N}$ are the linear, identity and nonlinear operators in (22), respectively. The above equation (27) can be written more explicitly

which can be numerically solved using time-stepping integration,

where $s$ and $\beta_q$ depend on the numerical scheme. Since this problem is a non-stiff problem, we usually use the fourth-order Adams-Bashforth method, e.g., ode113 in Matlab. From (29), we see that

To reduce the error ${\mathcal E}^{\mbox{ ROM}}$, we use a closure model which can be written as:

where $\alpha_k (t)$ is the Least-Squares (LS) projection modal coefficients which can be obtained by:

where the angle parentheses refer to the Euclidean inner product in ${\mathbb{R}}^K$. The LS projection modal coefficients include the hidden physics and its interaction with the dynamical core of the system. For simplicity, we assume that ${C_k}$ depend only on $a_k$ and $\alpha_k$ for each $k$. The correction function $C_k$ can be defined in various way, especially and naturally,

and

which are used in [25] and [26], respectively. In theory, they are different but we can not find the differences in computational performance. So, we will study more on these problem in the later articles. In this article, we adopt the equation (33) to the correction step in each time step $t_n$. One can use any of the suitable machine learning algorithm to learn this correction term $C_k$ such as ResNET and LSTM network, so on. In this article, we employ LSTM neural network algorithm to learn the mapping from LS projection modal coefficients to the correction term

5.   Computational procedure and numerical results

5.1.   Full order finite element solutions

Let $\sigma_N = \{ t_n \}^N_{n = 0}$ be a partition of $[0,T]$ into equal intervals $\Delta t = T/N$ with $t_0 = 0$ and $t_N = T$. The finite element spaces are chosen to be piecewise quadratic for the velocity and linear for the pressure, i.e., the Taylor-Hood finite element pair, based on a rectangular mesh. We use FreeFEM++ [20].

In this test we are interested in the convergence history for the parameters involved and so a simple stationary target velocity ${\boldsymbol U} = (U,V)$ is chosen, where

For the comparison purpose with the results in [19], we set up the initial velocity as follows:

This initial velocity rotates in the opposite direction to the target velocity ${\boldsymbol U}$. For this calculation, we use $\Delta t = 0.0001$, $h = 1/64$ and $T = 1.0$. The number of unknowns (degree of freedom) is 49923.

The flow evolution is given in Fig. 2. The controlled fluid is on the left and the desired flow is on the right. All the figures are normalized. At the beginning, we have a reduction in the magnitude, then we have a change in shape, and finally a change in magnitude again. This evolution is typical of linear feedback control. As shown in Fig. 2, the change in shape is so quick that it is difficult to see the full evolution without using very small time steps. The error $||{\boldsymbol u}(t) - {\boldsymbol U}(t) ||$ rapidly goes to zero and a good match is reached at $t = 0.5.$ To compare computational results for our feedback law with those in [19], we compute $||{\boldsymbol u}(t)-{\boldsymbol U}(t)||_{{\bf L}^2(\Omega)}$ for $\gamma = 10, 30$. Fig. 3 shows the $L^2$ norm $||{\boldsymbol u}(t) - {\boldsymbol U}(t) ||_{{\bf L}^2(\Omega)}$ between the calculated discretized controlled flow (here denoted by ${\boldsymbol u}$ instead of ${\boldsymbol u}_h$) and the exact target flow ${\boldsymbol U}$. We see that our feedback law is almost same effect to that of [19].

5.2.   POD-ROM

Now, we set up the Reynolds number as $Re = 10$ and the initial velocity for our computational experiments as follows:

We generate the data snapshots for $\gamma = [10, 20, 30, 40]$ with a time-step of 0.001 from time $t = 0$ to $T = 0.25$. The $||{\boldsymbol u}(t)-{\boldsymbol U}(t)||_{{\bf L}^2({\Omega})}$ for each cases of $\gamma$ is shown in Fig. 4.

The POD reduced basis is determined from the corresponding snapshot set as described in Section 3. Note that each basis function satisfies the discretized continuity equation, i.e., it is discretely solenoidal. The eight-dimensional POD basis functions are displayed in Fig. 5 (Upper two rows). We see that the difference of basis between those of $\gamma = 10$ and $\gamma = 40$.

We can compute the energy retained by POD basis functions using a relative information content (RIC) formula as given below:

where $\sigma$'s are singular values and ${\tilde N}$ is the total number of singular values. RIC represents the fraction of information of the the total data that can be recovered using $K$ basis functions.

Fig. 6 displays the convergence of relative information content with respect to the number of POD basis functions used to represent the reduced order system for different control number $\gamma$. For the snapshot set determined as described in Section 2.2 with $\gamma = [10, 20, 30, 40]$, the first 8 singular values of the corresponding snapshot matrix are listed in Table 1.

We retained eight basis functions (i.e., $K = 8$) as they captured more than $99.997\%$ of the energy for all control numbers $\gamma = [10,20,30,40]$ (see Fig. 6 and Table 1 for $\gamma = 40$). It can be seen that the singular values for the feedback control problem decreases rapidly. We consider here, a "small-dimensional" POD basis can capture most of the information contained in the snapshot set.

Now, given the desired state ${{\boldsymbol U}}$, we use $K$-dimensional system of of nonlinear ordinary differential equations (24) to determine reduced-order solutions of the Navier-Stokes system. Approximations of solutions ${\boldsymbol {a}}_k^{(n)}$ of the system of ordinary differential equation (27) are determined using a fourth order Adams-Moulton method. To quantify the performance of different frameworks, we defined the root mean squared error (RMSE) between the FOM velocity field and the velocity field computed with different ROM frameworks. The RMSE is defined as:

and

Fig. 7 shows the RMSE$(t)$ over time for four cases $\gamma = [10,20,30,40]$. The LS Projection error is ${\mathcal E}_{{{{\boldsymbol U}}_{K}}^{\perp}}(t,{\boldsymbol x}) = {\boldsymbol u}^{\mbox{ FOM}}(t,{\boldsymbol x}) - {\boldsymbol u}^{\mbox{ Proj}}(t,{\boldsymbol x})$ and GP-ROM error is ${\mathcal E}^{\mbox{ ROM}}(t,{\boldsymbol x}) = {\boldsymbol u}^{\mbox{ FOM}}(t,{\boldsymbol x}) - {\boldsymbol u}^{\mbox{ ROM}}(t,{\boldsymbol x})$. We calculate the LS projection modal coefficients by orthogonally projecting the state variables onto the reduced spaces. We will use the differences between RSME of LS projection and that of GP-ROM in the LSTM-ROM study, which will be studied in the next subsection.

5.3.   LSTM-ROM

In this subsection, we test GP-LSTM closure model for the out-of-sample condition $\gamma = 50$. For testing our GP-LSTM closure model, we have generated four cases in different values of parameters such that $\gamma = [10,20,30,40]$.

The LS projection modal coefficients include the hidden physics and its interaction with the dynamical core of the system. We can then define the correction term as:

A supervised learning framework is applied to model the correction term $C$ with information obtained from FOM and projection data. We train our LSTM neural network to learn the mapping from GP modal coefficients to the correction term. Since GP modal coefficients is used as input features to the LSTM network, the parameter $\gamma$ governing the system's behavior is taken implicitly into account. We train the LSTM network with two hidden layers and 120 cells.

Fig. 8 shows a sketch of Long Short-Term Memory unit and an architecture of LSTM network training. For more details for LSTM architecture, one can refer to the articles [27] and references therein. Programs obtained from GitHub, namely ETC-ROM Master, Hybrid-Modeling Master, UROM Master and mnni-rom Master, are adopted, modified and used in the learning process.

Once the model is trained, we could correct the GP modal coefficients with LSTM-based correction to approximate true projection modal coefficients.

We give the numerical procedure as follows:

1. Data Generation: For $\gamma = 10, 20, 30, 40$ to use in learning algorithm and $\gamma = 50$ to use a test.

2. Basis Construction: For $\gamma = 10, 20, 30, 40$.

3. The velocity modal coefficients $a_k^{(n)}$ of Galerkin projection reduced order modeling calculation for $\gamma = 10, 20, 30, 40$ using each corresponding basis.

4. LSTM Training using $C_k^{(n+1)}$ and $a_k^{(n+1)}$. A summary of the adopted hyper-parameters is presented in Table 2. The training and validation loss during the learning procedure are in Fig. 9.

5. Prediction for $\gamma = 50$ using the basis which is obtained from the snapshot sets with $\gamma = 40$ (see Fig. 10).

6. Reconstruction and compare with GP-ROM and LS Projection (see Fig. 11).

In step 5, interpolation techniques such as Grassman manifold interpolation [3,4,25], or the discrete empirical interpolation method, (DEIM) are typically applied to postprocess the results. For this work, however, such interpolations were not applied, so as to focus on the effects obtained from the Deep Learning technique. Thus, we use the same basis which is obtained from the snapshot sets with $\gamma = 40$. The POD modal coefficients predicted by these approaches are in Figure 10. As shown in the figure 10, it can be seen that the error of GP-LSTM is much smaller than that of GP-ROM.

Finally, we report the "online" computing time for our numerical experiments. We show the computational time as well as the RMSE of reconstructed fields at the final time $t = 0.25$ for the test cases at $\gamma = 50$ in Table 3. All calculations were performed on an iMac 3.6Ghz 8 Core. We observe that computing time of GP-LSTM is about twice that of GP-ROM(8) and the error of GP-LSTM(8) is smaller (about 1/10) than that of GP-ROM(8). If we use as much data as and use interpolations like in other papers, for example ([3,25,26]), one can reduce the error even more. However, in this papers, LSTM is applied in the simplest way, to see the possibility of applying the deep learning technique to the optimal control problem or UQ problem of fluid flows.

6.   Concluding remarks

In this study, a simple feedback rule is developed in order to reduce the amount of computation. It has been shown that our feedback law works very well. It is believed that mathematical proof will be possible without much difficulty. For real-time computation, a numerical experiment was performed by adopting a GP-ROM. In order to increase the GP-ROM's accuracy, the deep learning method, especially the LSTM method, which has been actively developed recently, was studied and applied. The ROM using deep learning such as LSTM performed in this study is considered to be worth continuing research. However, it is difficult to study systematically because the mathematical theory is not supported. We intend to apply GP-LSTM and GP-ResNet methods to the next studies, the optimal control problems and the uncertainty quantification problems of fluid flows.

Acknowledgments

We thank the reviewer for the comments that helped improve the paper.