# Semidiscrete approximation of the penalty approach to the stabilization of the Boussinesq system by localized feedback control

Mejdi Azaiez, Kévin Le Balc'h

► To cite this version:

Mejdi Azaiez, Kévin Le Balc'h. Semidiscrete approximation of the penalty approach to the stabilization of the Boussinesq system by localized feedback control. 2021. hal-03442394

HAL Id: hal-03442394

https://hal.science/hal-03442394

Preprint submitted on 23 Nov 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# Semidiscrete approximation of the penalty approach to the stabilization of the Boussinesq system by localized feedback control

Mejdi Azaiez* Kévin Le Balc'h†

November 23, 2021

## Abstract

We study the numerical approximation of the stabilization of the semidiscrete linearized Boussinesq system around an unstable stationary state. The stabilization is achieved by internal feedback controls applied on the velocity and the temperature equations, localized in an arbitrary open subset. This article follows the framework of [11], considering the continuous linearized Boussinesq system. The goal is to study the approximation by penalization of the free divergence condition in the semidiscrete case. More precisely, considering infinite time horizon LQR optimal control problem, we establish convergence results for the optimal controls, optimal solutions and Riccati operators when the penalization parameter goes to zero. We then propose a numerical validation of these results in a two-dimensional setting.

**Keywords:** Boussinesq system, penalty method, stabilization, feedback control

**MSC classification:**

## Contents

<table><tr><td>1</td><td>Introduction</td><td>2</td></tr><tr><td>2</td><td>Semidiscrete approximations</td><td>4</td></tr><tr><td>2.1</td><td>Finite dimensional approximations of the systems</td><td>4</td></tr><tr><td>2.2</td><td>Finite-dimensional linear controlled systems</td><td>5</td></tr><tr><td>2.3</td><td>Some background on LQR optimal control problems</td><td>7</td></tr><tr><td>3</td><td>Main result</td><td>7</td></tr><tr><td>4</td><td>Proof of the main result</td><td>8</td></tr><tr><td>4.1</td><td>Uniform bounds, convergence results for the free solution and the input-to-state maps</td><td>9</td></tr><tr><td>4.2</td><td>Uniform null-controllability and stabilization results</td><td>11</td></tr><tr><td>5</td><td>Numerical experiments and assessments</td><td>12</td></tr><tr><td>5.1</td><td>First example</td><td>12</td></tr><tr><td>5.2</td><td>Second example: Flow past a square section cylinder</td><td>24</td></tr></table>

<footnote>*Université de Bordeaux, CNRS, Bordeaux-INP, I2M UMR5295, F-33400, Talence, France. E-mail: azaiez@bordeaux-inp.fr</footnote>

<footnote>†Inria, Sorbonne Université, Université de Paris, CNRS, Laboratoire Jacques-Louis Lions, Paris, France. E-mail: kevin.le-balc-h@inria.fr</footnote>
# 1 Introduction

The optimal control of the Boussinesq system and of its linearization around a stationary state are matters of great interest in various applications fields, such as designing and exploiting energy efficient buildings, (see, for instance, [3], [17] and [5], [13]).

In the recent paper [11], the authors study the approximation by penalization of the free divergence condition for infinite time horizon Linear Quadratic Regulator (LQR) optimal control problem associated to the (continuous) linearized Boussinesq system around a stationary state. From a theoretical point of view, the first interest of considering such an approximation is to write the system as a well-posed control system in the sense of Salomon-Weiss. From a numerical point of view, the second interest is to avoid projection methods to treat the free divergence condition. Indeed, the matrix of the Leray projector is of big size and is not sparse. Note however that projection based efficient numerical methods have been recently developed in the literature, see for instance [2, 14, 9].

In this article, we study the numerical stabilization of the linearized Boussinesq system by localized feedback controls. In particular, we investigate the penalty approach for the semidiscrete linearized Boussinesq system.

Let us first introduce the model that we consider. Let $\Omega$ be a smooth domain contained in $\mathbb{R}^d$ with $d = 2, 3$, let $\mathcal{O}$ be a nonempty open subset contained in $\Omega$, let $\Gamma := \partial\Omega$ and let $n$ be the outer unit normal vector. The incompressible Boussinesq system with Neumann boundary conditions writes as follows

$$ egin{cases} 
\partial_t v - \operatorname{div} \sigma(v, p) + (v \cdot \nabla)v = ye_d + f & \text{in } (0, \infty) \times \Omega, \\
\partial_t y - \alpha \Delta y + v \cdot \nabla y = g & \text{in } (0, \infty) \times \Omega, \\
\operatorname{div} v = 0 & \text{in } (0, \infty) \times \Omega, \\
\sigma(v, p)n = h, \quad \partial_n y = k & \text{on } (0, \infty) \times \Gamma, \\
v(0, \cdot) = v_0, \quad y(0, \cdot) = y_0 & \text{in } \Omega.
\end{cases} 
\quad (1.1) $$

In (1.1), $v$ denotes the fluid velocity, $p$ is the fluid pressure, $y$ is the temperature of the fluid, $\sigma(v, p) = \nu((\nabla v) + (\nabla v)^{\mathrm{tr}}) - pI$ is the Cauchy stress tensor, $\nu > 0$ is the kinematic viscosity of the fluid, $\alpha > 0$ is the heat conductivity of the fluid, $e_d$ is the last vector of the canonical basis of $\mathbb{R}^d$. The terms $f: \Omega \to \mathbb{R}^d$, $g: \Omega \to \mathbb{R}$ describe respectively the influence of internal field forces and heat sources. The boundary term $h: \Gamma \to \mathbb{R}^d$ is an open boundary condition which is an example of non-reflecting outlet boundary, it is known to be efficient for low Reynolds number, see [12]. The Neumann boundary condition $k: \Gamma \to \mathbb{R}$ prescribes the heat flux. We assume that $(v_s, p_s, y_s) \in W^{1,\infty}(\Omega; \mathbb{R})^{d+2}$ is a real-valued solution to the stationary Boussinesq system

$$ egin{cases} 
-\operatorname{div} \sigma(v_s, p_s) + (v_s \cdot \nabla)v_s = y_s e_d + f & \text{in } \Omega, \\
-\alpha \Delta y + v_s \cdot \nabla y_s = g & \text{in } \Omega, \\
\operatorname{div} v_s = 0 & \text{in } \Omega, \\
\sigma(v, p)n = h, \quad \partial_n y = k & \text{on} \Gamma, 
\end{cases} 
\quad (1.2) $$

We consider the linearized incompressible Boussinesq system around $(v_s, p_s, y_s)$

$$ egin{cases} 
\partial_t v - \operatorname{div} \sigma(v, p) + (v \cdot \nabla)v_s + (v_s \cdot \nabla)v = ye_d + \tilde{u}1_{\mathcal{O}} & \text{in } (0, \infty) \times \Omega, \\
\partial_t y - \alpha \Delta y + v_s \cdot \nabla y + v \cdot \nabla y_s = u_{d+1}1_{\mathcal{O}} & \text{in } (0, \infty) \times \Omega, \\
\operatorname{div} v = 0 & \text{in } (0, \infty) \times \Omega, \\
\sigma(v, p)n = 0, \quad \partial_n y = 0 & \text{on } (0, \infty) \times \Gamma, \\
v(0, \cdot) = v_0, \quad y(0, \cdot) = y_0 & \text{in } \Omega.
\end{cases} 
\quad (1.3) $$

In the controlled system (1.3), $(v, p, y): (0, \infty) \times \Omega \to \mathbb{R}^{d+2}$ is the state, while $u := (\tilde{u}, u_{d+1})$: $(0, \infty) \times \mathcal{O} \to \mathbb{R}^{d+1}$ is the control, i.e. a localized field force $\tilde{u}$ and a heat source $u_{d+1}$ that one
can choose in order to modify the dynamics of the velocity $y$, the pressure $p$ and the temperature $y$.

It is worth mentioning that without any extra assumption on $(v_s, y_s)$, (1.3) with control $u = 0$ can be unstable. One way to stabilize (1.3) by feedback controls consists in looking for a control $u$ as a solution to a quadratic minimization problem, called LQR optimal control problem, see Section 2.3 for a very brief introduction to this notion. Such a control will be denoted by $u_{\text{opt}}$ and the associated trajectory, i.e. the solution to (1.3) with $u = u_{\text{opt}}$ by $(v_{\text{opt}}, p_{\text{opt}}, y_{\text{opt}})$.

Now we introduce a perturbed form of (1.3) which may be easier to solve in practice. More precisely, we penalize the free divergence condition. This method has been introduced for the first time in [15] for the Navier-Stokes system. Taking $\varepsilon > 0$ a small parameter, we consider the system

$$ egin{cases} 
\partial_t v - \operatorname{div} \sigma(v, p) + (v \cdot \nabla) v_s + (v_s \cdot \nabla) v = y e_d + \tilde{u} 1_O & \text{in } (0, \infty) \times \Omega, \\
\partial_t y - \alpha \Delta y + v_s \cdot \nabla y + v \cdot \nabla y_s = u_{d+1} 1_O & \text{in } (0, \infty) \times \Omega, \\
\operatorname{div} v + \varepsilon p = 0 & \text{in } (0, \infty) \times \Omega, \\
\sigma(v, p) n = 0, \partial_n y = 0 & \text{on } (0, \infty) \times \Gamma, \\
v(0, \cdot) = v_0, y(0, \cdot) = y_0 & \text{in } \Omega.
\end{cases} 
\quad (1.4) $$

By eliminating the pressure in (1.4), we can reformulate the above equations as the following parabolic system

$$ egin{cases} 
\partial_t v - \operatorname{div} \sigma(v, p) + (v \cdot \nabla) v_s + (v_s \cdot \nabla) v - \frac{1}{\varepsilon} \nabla \operatorname{div} v = y e_d + \tilde{u} 1_O & \text{in } (0, \infty) \times \Omega, \\
\partial_t y - \alpha \Delta y + v_s \cdot \nabla y + v \cdot \nabla y_s = u_{d+1} 1_O & \text{in } (0, \infty) \times \Omega, \\
\sigma(v, p) n = 0, \partial_n y = 0 & \text{on } (0, \infty) \times \Gamma, \\
v(0, \cdot) = v_0, y(0, \cdot) = y_0 & \text{in } \Omega.
\end{cases} 
\quad (1.5) $$

As before, let us denote by $u_{\text{opt},\varepsilon}$ the optimal control for the associated LQR problem to (1.5) and $(v_{\text{opt},\varepsilon}, y_{\text{opt},\varepsilon})$ the corresponding trajectory, i.e. the solution $(v, y)$ of (1.5) with $u = u_{\text{opt},\varepsilon}$.

In [11], for homogeneous Dirichlet boundary conditions instead of homogeneous Neumann boundary conditions that are considered here, the authors roughly prove that

$$ \lim_{\varepsilon \to 0} u_{\text{opt},\varepsilon} = u_{\text{opt}} \quad \text{and} \quad \lim_{\varepsilon \to 0} \begin{bmatrix} v_{\text{opt},\varepsilon} \\ y_{\text{opt},\varepsilon} \end{bmatrix} = \begin{bmatrix} v_{\text{opt}} \\ y_{\text{opt}} \end{bmatrix}, \quad (1.6) $$

see [11, Theorem 3.2] for a precise version of this statement.

The goal of this article consists in obtaining the convergence results (1.6) for semidiscrete approximations of the continuous systems (1.3) and (1.5), then proposing a numerical validation of these results.

The organization of the paper is as follows. In Section 2, we introduce the semidiscrete approximations of the linearized Boussinesq system and recall standard results related to LQR problems. In Section 3, we present the main results of the article, that is to say, convergence results for the optimal control and the optimal trajectories. In Section 4, we give the proofs of the main results by insisting only on the main differences between the proofs in the continuous case and the semidiscrete case. Finally, in Section 5, we proceed to the numerical validation of our results for several examples in a two-dimensional setting.

**Acknowledgements.** The authors have benefited from financial support of the French State, managed by the French National Research Agency (ANR) in the frame of the “Investments for the future” Programme IdEx Bordeaux - SysNum (ANR-10-IDEX-03-02). The last author has been supported by the grant “Numerical simulation and optimal control in view of temperature regulation in smart buildings” of the Nouvelle Aquitaine Region. The last author is partially supported by the Project TRECOS ANR-20-CE40-0009 funded by the ANR
(2021-2024). The authors would like to thank Marius Tucsnak for his advises and discussions contributing to the improvement of the manuscript.

# 2 Semidiscrete approximations

## 2.1 Finite dimensional approximations of the systems

The goal of this part is to introduce the discrete spatial approximations, based on a finite element method, associated to the controlled systems (1.3) and (1.5).

To approximate the systems (1.3) and (1.5) by a finite element method, we introduce finite dimensional subspaces $V_h \subset H^1(\Omega; \mathbb{R}^d)$ for the velocity, $M_h \subset L^2(\Omega; \mathbb{R})$ for the pressure, $Y_h \in H^1(\Omega; \mathbb{R})$ for the temperature and $U_h \subset L^2(\mathcal{O}; \mathbb{R}^{d+1})$ for the control. We denote by $(\phi_i)_{1\le i\le N_v}$ a basis of $V_h$, by $(\psi_i)_{1\le i\le N_p}$ a basis of $M_h$, by $(\theta_i)_{1\le i\le N_y}$ a basis of $Y_h$ and by $(\zeta_i)_{1\le i\le N_u}$ a basis of $U_h$. We set

$$v = \sum_{i=1}^{N_v} V_i \phi_i, \quad p = \sum_{i=1}^{N_p} p_i \psi_i, \quad y = \sum_{i=1}^{N_y} Y_i \theta_i, \quad u = \sum_{i=1}^{N_u} U_i \zeta_i, \qquad (2.1)$$

$$v_0 = \sum_{i=1}^{N_v} V_{0,i} \phi_i, \quad y_0 = \sum_{i=1}^{N_y} Y_{0,i} \theta_i. \qquad (2.2)$$

If we denote by boldface letters the coordinate vectors, we have

$$\mathbf{v} = (V_1, \dots, V_{N_v})^T, \quad \mathbf{p} = (p_1, \dots, p_{N_p})^T, \quad \mathbf{y} = (Y_1, \dots, Y_{N_y})^T, \quad \mathbf{u} = (U_1, \dots, U_{N_u})^T, \qquad (2.3)$$

$$\mathbf{v}_0 = (w_{0,1}, \dots, V_{0,N_v})^T, \quad \mathbf{y}_0 = (Y_{0,1}, \dots, Y_{0,N_y})^T. \qquad (2.4)$$

The finite dimensional approximation of (1.3) consists in finding $v \in H_{\text{loc}}^1([0, \infty); V_h)$, $p \in L_{\text{loc}}^2([0, \infty); M_h)$, $y \in H_{\text{loc}}^1([0, \infty); Y_h)$, such that

$$\frac{d}{dt} \int_{\Omega} (v(t) \cdot \phi) \, dx = a_0^{vv}(v(t), \phi) + b(\phi, p(t)) + a_0^{vy}(y(t), \phi) + \langle \tilde{u}(t), \phi \rangle_{L^2(\mathcal{O})^d}, \quad \forall \phi \in V_h, \qquad (2.5)$$

$$\frac{d}{dt} \int_{\Omega} (y(t)\theta) \, dx = a_0^{yy}(v(t), \theta) + a_0^{yy}(y(t), \theta) + \langle u_{d+1}(t), \theta \rangle_{L^2(\mathcal{O})}, \quad \forall \theta \in Y_h, \qquad (2.6)$$

$$b(v(t), \psi) = 0, \quad \forall \psi \in M_h, \qquad (2.7)$$

where

$$a_0^{vv}(v, \phi) = -\nu \int_{\Omega} \nabla v \cdot \nabla \phi \, dx - \int_{\Omega} [(v_s \cdot \nabla)v + (v \cdot \nabla)v_s] \cdot \phi \, dx, \qquad (2.8)$$

$$b(\phi, p) = \int_{\Omega} (\operatorname{div} \phi)p \, dx, \qquad (2.9)$$

$$a_0^{vy}(y, \phi) = \int_{\Omega} y \phi_d \, dx, \qquad (2.10)$$

$$a_0^{yy}(v, \theta) = - \int_{\Omega} [v \cdot \nabla y_s] \cdot \theta \, dx, \qquad (2.11)$$

$$a_0^{yy}(y, \theta) = -\alpha \int_{\Omega} \nabla y \cdot \nabla \theta \, dx - \int_{\Omega} [v_s \cdot \nabla y] \cdot \theta \, dx. \qquad (2.12)$$
We introduce the stiffness matrices $A_{vv}$, $A_{vp}$, $A_{vy}$, $A_{yv}$, $A_{yy}$; the mass matrices $M_{vv}$, $M_{pp}$, $M_{yy}$; and the matrices of control $B_{uv}$, $B_{uy}$ by

$$
\begin{align*}
A_{vv}^{ij} &= a_0^{vv}(\phi_j, \phi_i), \quad A_{vp}^{ij} = b(\phi_i, \psi_j), \quad A_{vy}^{ij} = a_0^{vy}(\theta_j, \phi_i), \quad A_{yv}^{ij} = a_0^{vy}(\theta_j, \theta_i), \\
M_{zz}^{ij} &= (\phi_i, \phi_j)_{L^2(\Omega)^d}, \quad M_{yy}^{ij} = (\theta_i, \theta_j)_{L^2(\Omega)}, \\
B_{uv}^{ij} &= (\tilde{\zeta}_i, \phi_j)_{L^2(\mathcal{O})^d}, \quad B_{uy}^{ij} = ((\zeta_i)_{d+1}, \theta_j)_{L^2(\mathcal{O})}.
\end{align*}
$$

System (2.5)-(2.7) may be written in the following form

$$
M_{vv} \mathbf{v}'(t) = A_{vv} \mathbf{v}(t) + A_{vp} \mathbf{p}(t) + A_{vy} \mathbf{y}(t) + B_{uv} \tilde{\mathbf{u}}(t), \quad (2.13)
$$

$$
M_{yy} \mathbf{y}'(t) = A_{yy} \mathbf{y}(t) + A_{yy} \mathbf{y}(t) + B_{uy} \mathbf{u}_{d+1}(t), \quad (2.14)
$$

$$
A_{vp}^{T} \mathbf{v}(t) = 0. \tag{2.15}
$$

Similarly for the perturbed form (1.5) the finite dimensional approximation consists in finding
$v \in H_{\text{loc}}^{1}([0, \infty); V_h)$, $y \in H_{\text{loc}}^{1}([0, \infty); Y_h)$, such that

$$
\frac{d}{dt} \int_{\Omega} v(t) \cdot \phi dx = a_{\varepsilon}^{vv}(v(t), \phi) + a_{0}^{vy}(y(t), \phi) + \langle \tilde{u}(t), \phi \rangle_{L^2(\mathcal{O})^d}, \quad \forall \phi \in V_h, \quad (2.16)
$$

$$
\frac{d}{dt} \int_{\Omega} y(t) \theta dx = a_0^{yv}(v(t), \theta) + a_0^{yy}(y(t), \theta) + \langle u_{d+1}(t), \theta \rangle_{L^2(\mathcal{O})}, \quad \forall \theta \in Y_h, \quad (2.17)
$$

where

$$
\begin{align}
a_{\epsilon}^{vv}(v, \phi) &= -\nu \int_{\Omega} \nabla v : \nabla \phi \,dx - \int_{\Omega} [(v_s \cdot \nabla)v + (v \cdot \nabla)v_s] \cdot \phi \,dx - \frac{1}{\epsilon} \int_{\Omega} (\operatorname{div} v)(\operatorname{div} \phi) \,dx. \tag{2.18} \\
&\quad - \nu \int_{\Omega} \nabla v : \nabla \phi \,dx - \int_{\Omega} [(v_s \cdot \nabla)v + (v \cdot \nabla)v_s] \cdot \phi \,dx - \frac{1}{\epsilon} \int_{\Omega} (\operatorname{div} v)(\operatorname{div} \phi) \,dx. \tag{2.18}
\end{align}
$$

and the other bilinear forms have already been defined in (2.10), (2.11) and (2.12).

System (2.16)-(2.17) may be written in the following form

$$
M_{vv} \mathbf{v}'(t) = A_{vv,\epsilon} \mathbf{v}(t) + A_{vy} \mathbf{y}(t) + B_{uv} \tilde{\mathbf{u}}(t), \quad (2.19)
$$

$$
M_{yy} \mathbf{y}'(t) = A_{yy} \mathbf{y}(t) + A_{yy} \mathbf{y}(t) + B_{uy} \mathbf{u}_{d+1}(t). \quad (2.20)
$$

2.2 Finite-dimensional linear controlled systems

The goal of this part is to first reformulate (2.13)-(2.15) into a more suitable differential system,
without the constraint formula (2.15). This will be done thanks to the projection of the dynamics
by the discrete Leray projector. Then, we state a well-posedness result for both finite-dimensional
linear controlled systems.

In the following result whose proof is an easy adaptation of [1, Proposition 3.1], we define
the discrete Leray projector $P_h$ and state some important properties for $P_h$.

**Proposition 2.1.** The projector $P_h$ in $\mathbb{R}^{N_v}$ onto $\ker(A_{vp}^T M_{vv}^{-1})$ parallel to $\operatorname{Im}(A_{vp})$ and the projector $P_h^T$ in $\mathbb{R}^{N_v}$ onto $\ker(A_{vp}^T)$ parallel to $\operatorname{Im}(M_{vv}^{-1} A_{vp})$ are defined by

$$
P_h = I - A_{vp}(A_{vp}^T M_{vv}^{-1} A_{vp})^{-1} A_{vp}^T M_{vv}^{-1}, \quad P_h^T = I - M_{vv}^{-1} A_{vp}(A_{vp}^T M_{vv}^{-1} A_{vp})^{-1} A_{vp}^T. \quad (2.21)
$$

Moreover, we have $P_h A_{vp} = 0$, $P_h M_{vv} = M_{vv} P_h^T$, $M_{vv}^{-1} P_h = P_h^T M_{vv}^{-1}$.
The above result enables us to project the dynamics of (2.13)-(2.15), we then obtain

$$
\begin{align*}
M_{vv} P_h^T \mathbf{v}'(t) &= A_{vv} P_h^T \mathbf{v}(t) + A_{vy} \mathbf{y}(t) + P_h B_{uv} \tilde{\mathbf{u}}(t), \\
M_{yy} \mathbf{y}'(t) &= A_{yy} \mathbf{v}(t) + A_{yy} \mathbf{y}(t) + B_{uy} \mathbf{u}_{d+1}(t),
\end{align*}
$$

and then using that $A_{vp}^T \mathbf{v}(t) = 0$ so $P_h^T \mathbf{v}(t) = \mathbf{v}(t)$, we thus obtain

$$
M_{vv} P_h^T \mathbf{v}'(t) = A_{vv} P_h^T \mathbf{v}(t) + A_{vy} \mathbf{y}(t) + P_h B_{uv} \tilde{\mathbf{u}}(t), \quad (2.22)
$$

$$
M_{yy} \mathbf{y}'(t) = A_{yy} P_h^T \mathbf{v}(t) + A_{yy} \mathbf{y}(t) + B_{uy} \mathbf{u}_{d+1}(t), \quad (2.23)
$$

$$
P_h^T \mathbf{v}(t) = \mathbf{v}(t). \tag{2.24}
$$

Before going further, we define the following matrices

$$
\begin{align*}
P &= \begin{bmatrix} P_h & 0 \\ 0 & I_{N_y} \end{bmatrix}, \quad
M = \begin{bmatrix} M_{vv} & 0 \\ 0 & M_{yy} \end{bmatrix}, \quad
A_0 = \begin{bmatrix} A_{vv} & A_{vy} \\ A_{yv} & A_{yy} \end{bmatrix}, \quad
A_\epsilon = \begin{bmatrix} A_{vv,\epsilon} & A_{vy} \\ A_{yv} & A_{yy} \end{bmatrix}, \\
B &= \begin{bmatrix} B_{uv} \\ B_{uy} \end{bmatrix}, \quad
B_0 = PB,
\end{align*}
$$

and finite-dimensional spaces

$$
H_h = \mathrm{Ker}(A_{vp}^T) \times Y_h, \quad X_h = V_h \times Y_h.
$$

We have the following standard well-posedness results for (2.22)-(2.24) and (2.19)-(2.20).

**Proposition 2.2.** For every $[\begin{matrix} v_0 \\ y_0 \end{matrix}] \in H_h$, $u \in L^2([0, \infty); U_h)$, the Cauchy problem

$$
M \begin{bmatrix} \dot{v} \\ \dot{y} \end{bmatrix} = A_0 \begin{bmatrix} v \\ y \end{bmatrix} + B_0 u, \quad \begin{bmatrix} v(0) \\ y(0) \end{bmatrix} = \begin{bmatrix} v_0 \\ y_0 \end{bmatrix}, \qquad (2.25)
$$

admits an unique solution $[\begin{matrix} v \\ y \end{matrix}] \in C^0([0, +\infty); H_h)$ defined by

$$
\begin{bmatrix} v(t) \\ y(t) \end{bmatrix} = e^{tM^{-1}A_0} \begin{bmatrix} v_0 \\ y_0 \end{bmatrix} + \int_0^t e^{(t-\sigma)M^{-1}A_0} M^{-1} B_0 u(\sigma) \, d\sigma \quad (t \ge 0). \quad (2.26)
$$

Similarly, for every $[\begin{matrix} v_0 \\ y_0 \end{matrix}] \in X_h$, $u \in L^2([0, \infty); U_h)$, the Cauchy problem

$$
M \begin{bmatrix} \dot{v} \\ \dot{y} \end{bmatrix} = A_{\epsilon} \begin{bmatrix} v \\ y \end{bmatrix} + Bu, \quad \begin{bmatrix} v(0) \\ y(0) \end{bmatrix} = \begin{bmatrix} v_0 \\ y_0 \end{bmatrix}, \qquad (2.27)
$$

admits an unique solution $[\begin{matrix} v \\ y \end{matrix}] \in C^0([0, +\infty); X_h)$, defined by

$$
\begin{bmatrix} v(t) \\ y(t) \end{bmatrix} = e^{tM^{-1}A_{\epsilon}} \begin{bmatrix} v_0 \\ y_0 \end{bmatrix} + \int_0^t e^{(t-\sigma)M^{-1}A_{\epsilon}} M^{-1} Bu(\sigma) \, d\sigma \quad (t \ge 0). \quad (2.28)
$$

The proof is standard so we omit it.
## 2.3 Some background on LQR optimal control problems

We recall below some basic facts on LQR problems in a finite dimensional context. Roughly speaking, the main result of the theory stipulates that the optimal control of the cost quadratic function is a feedback controller that can be explicitly determined by solving a Riccati equation.

Let $\mathcal{A} \in \mathbb{R}^{n \times n}$ and $\mathcal{B} \in \mathbb{R}^{n \times m}$. The LQR theory in infinite horizon focus on trajectories of the linear control system $x' = \mathcal{A}x + \mathcal{B}u$, $x(0) = x_0$, which minimize the quadratic cost functional

$$ \mathcal{J}(u; x_0) = \int_0^{+\infty} \left( \|u(t)\|_{\mathbb{R}^m}^2 + \|x(t)\|_{\mathbb{R}^n}^2 \right) dt \quad (u \in L^2(0, +\infty; \mathbb{R}^m)). \quad (2.29) $$

We recall below the well-known notion of stabilizability.

**Definition 2.3.** The pair $(\mathcal{A}, \mathcal{B})$ is stabilizable if there exist $\mathcal{K} \in \mathbb{R}^{m \times n}$, $M > 0$ and $\omega > 0$ such that

$$ \|e^{t(\mathcal{A}-\mathcal{B}\mathcal{K})}\| \leq Me^{-\omega t} \quad (t \geq 0). $$

One of the main results of the LQR theory states that the optimal control associated to (2.29) is given by a feedback law, see for instance [16, Théorème 4.4.5].

**Theorem 2.4.** If the pair $(\mathcal{A}, \mathcal{B})$ is stabilizable then for every $x_0 \in \mathbb{R}^n$, $J(\cdot; x_0)$ admits a unique minimum $u_{\text{opt}}$ given by

$$ u_{\text{opt}}(t) = -\mathcal{B}^*\mathcal{P}e^{t(\mathcal{A}-\mathcal{B}\mathcal{B}^*\mathcal{P})}x_0 \quad (t \in [0, \infty)). \quad (2.30) $$

where $P \in \mathbb{R}^{n \times n}$, is the unique nonnegative self-adjoint solution of the Riccati equation

$$ \mathcal{A}^*\mathcal{P} + \mathcal{P}\mathcal{A} - \mathcal{P}\mathcal{B}\mathcal{B}^*\mathcal{P} + I_n = 0. \quad (2.31) $$

Moreover, for every $x_0 \in \mathbb{R}^n$,

$$ \min_{u \in L^2([0, +\infty); \mathbb{R}^m)} \mathcal{J}(u; x_0) = \langle \mathcal{P}x_0, x_0 \rangle. \quad (2.32) $$

# 3 Main result

In this section, we use the notations introduced in Section 2 and we state the main result of the paper.

Assume that $(A_0, B_0)$ is stabilizable. By Theorem 2.4, let us denote by $\Pi_0 \in \mathcal{L}(H_h)$ the unique nonnegative self-adjoint solution of

$$ A_0^*\Pi_0 + \Pi_0 A_0 - \Pi_0 B_0 B_0^*\Pi_0 + I_{H_h} = 0. \quad (3.1) $$

In (3.1), $I_{H_h}$ denotes the identity operator on $H_h$.

In the same way, assume that $(A_\varepsilon, B)$ is stabilizable. Then let us denote by $\Pi_\varepsilon \in \mathcal{L}(X_h)$ the unique nonnegative self-adjoint solution of

$$ A_{\varepsilon}^*\Pi_{\varepsilon} + \Pi_{\varepsilon}A_{\varepsilon} - \Pi_{\varepsilon}BB^*\Pi_{\varepsilon} + I_{X_h} = 0. \quad (3.2) $$

In (3.2), $I_{X_h}$ denotes the identity operator on $X_h$.

The main result of the article focuses on asymptotic properties of $\Pi_\varepsilon$ in the limit $\varepsilon \to 0$.
**Theorem 3.1.** We assume that

$$ \mathrm{rank}(B) = N_v + N_y. \qquad (3.3) $$

Then, $(A_0, B_0)$ is exponentially stabilizable and $(A_\epsilon, B)$ is uniformly (with respect to $\epsilon$) exponentially stabilizable. Moreover, there exist $C > 0$ and $\omega > 0$, independent of $\epsilon$ such that the solution

$$ \|\Pi_{\epsilon}\|_{\mathcal{L}(X_h; X_h)} \le C \quad (\epsilon \ge 0), \qquad (3.4) $$

$$ \left\| e^{t(A_{\epsilon} - BB^{*}\Pi_{\epsilon})} \right\|_{\mathcal{L}(X_h; X_h)} \le C e^{-\omega t} \quad (\epsilon, t \ge 0), \qquad (3.5) $$

In addition,

$$ \lim_{\epsilon \to 0+} \| \Pi_{\epsilon} P - \Pi_{0} \|_{\mathcal{L}(H_h; X_h)} = 0 \quad (t \ge 0), \qquad (3.6) $$

$$ \lim_{\epsilon \to 0+} \left\| e^{t(A_{\epsilon} - BB^{*}\Pi_{\epsilon})} P - e^{t(A_0 - B_0 B_0^{*}\Pi_0)} \right\|_{\mathcal{L}(H_h; X_h)} = 0 \quad (t \ge 0), \qquad (3.7) $$

$$ \lim_{\epsilon \to 0+} \left\| -B^{*}\Pi_{\epsilon} e^{t(A_{\epsilon} - BB^{*}\Pi_{\epsilon})} P - (-B_{0}^{*}\Pi_{0} e^{t(A_{0} - B_{0}B_{0}^{*}\Pi_{0})}) \right\|_{\mathcal{L}(H_h; U_h)} = 0 \quad (t \ge 0), \qquad (3.8) $$

and the last two convergences are uniform with respect to $t$ on compact intervals.

Before going further, let us make some comments on Theorem 3.1.

* Theorem 3.1 is the analogue discrete version of [11, Theorem 3.2]. It is worth mentioning that the norm-convergence results of the Riccati operators (3.6), the closed-loop semigroups (3.7) and the optimal controls (3.8) are actually equivalent to ponctual convergence results because the functional spaces $H_h$, $X_h$ and $U_h$ are finite dimensional.

* The full-rank condition on $B$, i.e. (3.3) is used for obtaining null-controllability results for the pairs $(A_0, B_0)$ and $(A_\epsilon, B)$, see the proof of Proposition 4.5 below. This condition is checked numerically. It is worth mentioning that [11, Section 7] asserts that in the continuous case, the control matrix $B$ has not to satisfy a full-rank condition. Obtaining the same results of Theorem 3.1 without assuming (3.3) is an interesting open problem.

* Note that Theorem 3.1 proves convergence results for LQR problems associated to the pairs $(A_0, B_0)$ and $(A_\epsilon, B)$ but the proof that we will sketch can be easily adapted to the pairs $(M^{-1}A_0, M^{-1}B_0)$ and $(M^{-1}A_\epsilon, M^{-1}B)$. Therefore, one can recover the expected convergence results for LQR problems associated to the true semidiscrete systems (2.25) and (2.27), that we consider in the numerical validation, see Section 5.

* The convergence results of Theorem 3.1 are not uniform with respect to the parameter $h$ because the constants $C$ that appear in (3.4) and in (3.5) depend a priori on $h$. Obtaining such a uniform result in the spirit of [10] is an interesting open question.

# 4 Proof of the main result

In this section, we give the main points for obtaining the proof of Theorem 3.1. We crucially follow the methodology developed in [11] that we briefly recall here.

One of the key point consists in obtaining convergence results for the free solution, i.e. $e^{tA_\epsilon}P \to e^{tA_0}P$ and for the input-to-state maps i.e. $\int_0^t e^{(t-s)A_\epsilon}Bu(s)ds \to \int_0^t e^{(t-s)A_0}B_0u(s)ds$ as $\epsilon \to 0$. This type of results will be collected in Section 4.1. Then, to prove Theorem 3.1, we first establish uniform bounds (3.4) and (3.5). For this, we prove uniform null-controllability results for the pairs $(A_0, B_0)$ and the pairs $(A_\epsilon, B)$. While this type of results are obtained in the continuous case through observability estimates for the adjoint systems, here we simply
derive these results using the full-rank condition on the matrix B, i.e. (3.3) that implies that the so-called Kalman condition is fulfilled then null-controllability holds true. These results will be proved in Section 4.2. The last point consists in proving the convergence results (3.6), (3.7), (3.8) by following carefully [11, Section 6]. Some details of the proofs will be omitted, we will only insist on the main differences with [11].

## 4.1 Uniform bounds, convergence results for the free solution and the input-to-state maps

The goal of this part is to obtain convergence results for the free solution and the controlled solution, thanks to convergence results for the resolvents and the use of Trotter-Kato theorem. We have the following result.

**Proposition 4.1.** There exist $C \ge 1$ and $\omega \in \mathbb{R}$ such that for every $\varepsilon > 0$

$$ \|e^{tA_{\varepsilon}}\|_{\mathcal{L}(X_h; X_h)} \le C e^{\omega t} \quad (\varepsilon > 0, t \ge 0). \quad (4.1) $$

*Proof.* By integration by parts and Young's inequalities, it is not difficult to see that there exists a positive constant $c > 0$ depending on $\|(v_s, y_s)\|_{W^{1,\infty}(\Omega)^{d+1}}$ (independent of $\varepsilon > 0$), such that

$$ \begin{aligned} & - \langle A_{\varepsilon} \begin{pmatrix} \varphi \\ \zeta \end{pmatrix}, \begin{pmatrix} \varphi \\ \zeta \end{pmatrix} \rangle_{L^2(\Omega)^{d+1}} \\ & \ge \nu \|\nabla \varphi\|_{L^2(\Omega)^d}^2 + \frac{1}{\varepsilon} \|\operatorname{div} \varphi\|_{L^2(\Omega)^d}^2 + \alpha \|\nabla \zeta\|_{L^2(\Omega)^d}^2 - c \|(\varphi, \zeta)\|_{L^2(\Omega)^{d+1}}^2 \quad \left( \begin{bmatrix} \varphi \\ \zeta \end{bmatrix} \in X_h \right). \end{aligned} \quad (4.2) $$

We then easily deduce that

$$ \|e^{t(A_{\varepsilon}-cI)}\|_{\mathcal{L}(X_h; X_h)} \le 1 \quad (\varepsilon > 0, t \ge 0), $$

which gives the expected result. $\square$

We next prove the convergence of the resolvents of $A_{\varepsilon}$ (respectively $A_{\varepsilon}^*$) towards the resolvents of $A_0$ (respectively $A_0^*$).

**Proposition 4.2.** There exists $\lambda_0 > 0$ such that for every $\lambda \in \mathbb{C}$ with $\operatorname{Re}\lambda \ge \lambda_0$, we have

$$ \lim_{\varepsilon \to 0} \|[\lambda I - A_{\varepsilon}]^{-1} - [\lambda I - A_0]^{-1} P\|_{\mathcal{L}(X_h)} = 0, \quad (4.3) $$

$$ \lim_{\varepsilon \to 0} \|[\lambda I - A_{\varepsilon}^*]^{-1} - [\lambda I - A_0^*]^{-1} P\|_{\mathcal{L}(X_h)} = 0. \quad (4.4) $$

*Proof.* We only prove (4.3), since the proof of (4.4) is fully similar.

$$ \text{Let } \begin{bmatrix} f \\ g \end{bmatrix} \in X_h \text{ and} \qquad \lambda_0 > c > 0, \quad (4.5) $$

where the constant $c$ is the one appearing in (4.2). Then, for $\lambda \in \mathbb{C}$ with $\operatorname{Re}\lambda \ge \lambda_0$, by setting

$$ \begin{bmatrix} \varphi_{\varepsilon} \\ \zeta_{\varepsilon} \end{bmatrix} = (\lambda I - A_{\varepsilon})^{-1} \begin{bmatrix} f \\ g \end{bmatrix}, $$
we have

$$
\begin{align}
& \lambda \int_{\Omega} \varphi_{\varepsilon} \cdot \psi \, dx + \nu \int_{\Omega} \nabla \varphi_{\varepsilon} : \nabla \psi \, dx + \int_{\Omega} [(v_s \cdot \nabla)\varphi_{\varepsilon} + (\varphi_{\varepsilon} \cdot \nabla)v_s] \cdot \psi \, dx \nonumber \\
& + \frac{1}{\varepsilon} \int_{\Omega} (\operatorname{div} \varphi_{\varepsilon})(\operatorname{div} \psi) \, dx - \int_{\Omega} \zeta_{\varepsilon}\psi_d \, dx + \lambda \int_{\Omega} \zeta_{\varepsilon}\eta \, dx + \alpha \int_{\Omega} \nabla\zeta_{\varepsilon} \cdot \nabla\eta \, dx \nonumber \\
& + \int_{\Omega} [v_s \cdot \nabla\zeta_{\varepsilon} + \varphi_{\varepsilon} \cdot \nabla y_s] \cdot \eta \, dx = \int_{\Omega} f \cdot \psi \, dx + \int_{\Omega} g\eta \, dx \tag{4.6} \\
& \qquad \qquad \qquad \qquad \qquad \qquad \left( \begin{bmatrix} \psi \\ \eta \end{bmatrix} \in X_h \right). \nonumber
\end{align}
$$

Taking $\begin{bmatrix} \psi \\ \eta \end{bmatrix} = \begin{bmatrix} \varphi_{\varepsilon} \\ \zeta_{\varepsilon} \end{bmatrix}$ in (4.6) and using (4.2) we obtain

$$
\|\varphi_{\varepsilon}\|_{H^1(\Omega)^d}^2 + \frac{1}{\varepsilon} \|\operatorname{div} \varphi_{\varepsilon}\|_{L^2(\Omega)}^2 + \|\zeta_{\varepsilon}\|_{H^1(\Omega)}^2 \le c \left\| \begin{bmatrix} f \\ g \end{bmatrix} \right\|_{X_h}^2 \quad (\varepsilon > 0),
$$

where $c > 0$ is some another constant. The above estimate implies that there exists $\begin{bmatrix} \varphi \\ \zeta \end{bmatrix} \in H_h$
such that

$$
\begin{bmatrix} \varphi_{\epsilon} \\ \zeta_{\epsilon} \end{bmatrix} \to \begin{bmatrix} \varphi \\ \zeta \end{bmatrix} \text{ as } \epsilon \to 0 \text{ in } X_h,
$$

For $\begin{bmatrix} \psi \\ \eta \end{bmatrix} \in H_h$, note that $\psi \in \operatorname{Ker}(A_{vp}^T)$, we can thus pass to the limit in (4.6) to obtain that

$$
\begin{align}
& \lambda \int_{\Omega} \varphi \cdot \psi \, dx + \nu \int_{\Omega} \nabla \varphi : \nabla \psi \, dx + \int_{\Omega} [(v_s \cdot \nabla)\varphi + (\varphi \cdot \nabla)v_s] \cdot \psi \, dx \nonumber \\
& - \int_{\Omega} \zeta \psi_d \, dx + \lambda \int_{\Omega} \zeta \eta \, dx + \alpha \int_{\Omega} \nabla \zeta \cdot \nabla \eta \, dx \nonumber \\
& + \int_{\Omega} [v_s \cdot \nabla \zeta + \varphi \cdot \nabla y_s] \cdot \eta \, dx = \int_{\Omega} f \cdot \psi \, dx + \int_{\Omega} g\eta \, dx, \tag{4.7} \\
& \qquad \qquad \qquad \qquad \qquad \qquad \left( \begin{bmatrix} \psi \\ \eta \end{bmatrix} \in H_h \right), \nonumber
\end{align}
$$

then using

$$
\int_{\Omega} f \cdot \psi \, dx + \int_{\Omega} g\eta \, dx = \int_{\Omega} P \begin{bmatrix} f \\ g \end{bmatrix} \cdot \begin{bmatrix} \psi \\ \eta \end{bmatrix} \, dx \quad \left( \begin{bmatrix} \psi \\ \eta \end{bmatrix} \in H_h \right), \quad (4.8)
$$

we deduce (4.3) from (4.7).

<figure><img src="image_1.png" /></figure>

An important consequence of Proposition 4.2 is the following result.

**Proposition 4.3.** We have

$$
\lim_{\epsilon \to 0} \| e^{tA_\epsilon} P - e^{tA_0} \|_{\mathcal{L}(H_h; X_h)} = 0 \quad (t \ge 0), \quad (4.9)
$$

$$
\lim_{\epsilon \to 0} \|(e^{tA_{\epsilon}})^* P - (e^{tA_0})^*\|_{\mathcal{L}(H_h; X_h)} = 0 \quad (t \ge 0), \quad (4.10)
$$

uniformly with respect to $t$ on compact intervals.

The proof comes from an adaptation of Trotter-Kato's theorem, see [11, Proposition 4.2] for details. Note again that the convergence are norm-convergence operators because the functional spaces are finite-dimensional.

For $u \in L^2([0, +\infty; U_h)$, we define the input-to-state maps

$$
\Xi_t^\epsilon u = \int_0^t e^{(t-s)A_\epsilon} Bu(s)ds \quad \text{and} \quad \Xi_t^0 u = \int_0^t e^{(t-s)A_0} B_0 u(s)ds. \quad (4.11)
$$

We have the following result.
**Proposition 4.4.** With the above notation, we have

$$ \lim_{\epsilon \to 0+} \| \Xi_t^\epsilon u - \Xi_t^0 u \|_{X_h} = 0 \quad (t \ge 0, u \in L^2([0, +\infty; U_h)), \quad (4.12) $$

$$ \lim_{\epsilon \to 0+} \|(\Xi_t^\epsilon)^* u - (\Xi_t^0)^* u\|_{X_h} = 0 \quad (t \ge 0, u \in L^2([0, +\infty; U_h)), \quad (4.13) $$

Moreover, the convergence above is uniform for $t$ in compact intervals.

The proof uses Lax-Phillips semigroups theory and Trotter-Kato's theorem, see [11, Proposition 4.5] for details.

## 4.2 Uniform null-controllability and stabilization results

The goal of this part is to state a uniform null-controllability result for the pair $(A_\epsilon, B)$. Here, the proof is different from the one in the continuous case that has been obtained through Carleman estimates. By using the full-rank condition on $B$, explicit controls are designed.

The main result of this part is the following.

**Proposition 4.5.** For every $\epsilon > 0$, the pair $(A_\epsilon, B)$ is exactly controllable. Moreover, for every $T > 0$, there exists $C > 0$ such that for every $\epsilon > 0$, $\begin{bmatrix} v_0 \\ y_0 \end{bmatrix} \in X_h$, there exists $u_\epsilon \in L^2([0, T]; U_h)$ satisfying

$$ \|u_\epsilon\|_{L^2([0,T];U_h)} \le C \left\| \begin{bmatrix} v_0 \\ y_0 \end{bmatrix} \right\|_{X_h}, \quad (4.14) $$

such that the solution to

$$ \begin{bmatrix} v'_\epsilon(t) \\ y'_\epsilon(t) \end{bmatrix} = A_\epsilon \begin{bmatrix} v_\epsilon(t) \\ y_\epsilon(t) \end{bmatrix} + B u_\epsilon(t), \quad \begin{bmatrix} v_\epsilon(0) \\ y_\epsilon(0) \end{bmatrix} = \begin{bmatrix} v_0 \\ y_0 \end{bmatrix}, \quad (4.15) $$

satisfies

$$ \begin{bmatrix} v_{\epsilon}(T) \\ y_{\epsilon}(T) \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \end{bmatrix}. \quad (4.16) $$

It is a well-known fact that in finite dimension, the exact-controllability (which is equivalent to the null-controllability) is equivalent to the Kalman rank condition

$$ \mathrm{rank}(B|A_\epsilon B| \dots | A_\epsilon^{N_v+N_y} B) = N_v + N_y, \quad (4.17) $$

then does not depend on the time of control $T > 0$. Here, we recall that we have assumed (3.3) so (4.17) trivially holds. So the main difficulty of Proposition 4.5 is to obtain the uniform bound on the control cost, i.e (4.14).

*Proof.* We recall that $B$ is full-rank by (3.3), so there exists $\tilde{B}$ such that $B\tilde{B} = I_{N_v+N_y}$. Let us fix $T > 0$ and $\begin{bmatrix} v_0 \\ y_0 \end{bmatrix} \in X_h$, and let us define the following controlled trajectory for $t \in [0, T]$,

$$ \begin{bmatrix} v_{\epsilon}(t) \\ y_{\epsilon}(t) \end{bmatrix} = \frac{T-t}{T} e^{tA_{\epsilon}} \begin{bmatrix} v_0 \\ y_0 \end{bmatrix}, \quad u_{\epsilon}(t) = -\frac{1}{T} \tilde{B} e^{tA_{\epsilon}} \begin{bmatrix} v_0 \\ y_0 \end{bmatrix}. \quad (4.18) $$

We check easily that

$$ \begin{bmatrix} v'_{\epsilon}(t) \\ y'_{\epsilon}(t) \end{bmatrix} = A_{\epsilon} \begin{bmatrix} v_{\epsilon}(t) \\ y_{\epsilon}(t) \end{bmatrix} + B u_{\epsilon}(t), \quad \begin{bmatrix} v(0) \\ y(0) \end{bmatrix} = \begin{bmatrix} v_0 \\ y_0 \end{bmatrix}, \quad \begin{bmatrix} v_{\epsilon}(T) \\ y_{\epsilon}(T) \end{bmatrix} = 0. \quad (4.19) $$

Finally, (4.14) is an easy consequence of the definition of $u_\epsilon$ and (4.1). □
A similar proof enables us to obtain the null-controllability of the pair $(A_0, B_0)$ in $H_h$.

By following carefully [11, Section 6], we can use Proposition 4.5 to prove uniform stabilization results, stated in (3.4) and (3.5) then convergence results (3.6), (3.7) and (3.8) which leads to the conclusion of the proof of Theorem 3.1

# 5 Numerical experiments and assessments

Prior to any numerical discussion, we briefly describe the finite element discretization we use in our computations for achieving the numerical validation of the method. To alleviate the presentation, we consider that $\Omega$ is a two dimensional domain. Assume that it is polygonal and is the union of a finite number of triangles. Then, we introduce a regular family $(\mathcal{T}_h)_h$ of triangulations of $\Omega$.

* The domain $\Omega$ is the union of all elements of $\mathcal{T}_h$;

* The intersection of two different elements of $\mathcal{T}_h$, if not empty, is a vertex or a whole edge of both of them;

* The ratio of the diameter $h_K$ of any element $K$ of $\mathcal{T}_h$ to the diameter of its inscribed circle is smaller than a constant $\sigma$ independent of $h$.

The mesh-size $h$ is the maximum of the diameters $h_K$. We refer to [4, 6, 7] for the basics of the finite element method.

Furthermore, it is well known that the velocity and pressure approximation spaces used in the spatial discretization are required to satisfy a discrete LBB condition, see [7] for more details concerning this point and more details about the finite element discretization of Navier-Stokes problem. In this work we choose to use the standard $P_2 - P_1$ and $P_2$ to respectively approximate the velocity, pressure and temperature fields. More precisely, we consider the discrete spaces defined by

$$Y_h = \{ \chi_h \in H^1(\Omega); \forall K \in \mathcal{T}_h, (\chi_h)|_K \in \mathcal{P}_2(K) \}, \quad V_h = (Y_h \cap H_0^1(\Omega))^2,$$

and

$$M_h = \{ \chi_h \in L^2(\Omega); \forall K \in \mathcal{T}_h, (\chi_h)|_K \in \mathcal{P}_1(K) \},$$

where $\mathcal{P}_i(K)$ stands for the space of restrictions to $K$ of polynomial functions of degree lower or equal to $i$ on $\mathbb{R}^2$.

## 5.1 First example

We first describe an example to assess and illustrate the main results of the paper, i.e. the convergence results (3.6), (3.7) and (3.8).

Let us introduce first the data used in the numerical experiments. The whole domain (see Fig. 1) is $\Omega = [0, 1]^2$, and the subset of control is

$$\mathcal{O} = \{(x_1, x_2) \in \Omega; (x_1 - 0.5)^2 + (x_2 - 0.5)^2 < 0.01\}. \quad (5.1)$$

We choose as stationary state the velocity and temperature given by

$$\begin{bmatrix} v_s^1(x_1, x_2) \\ v_s^2(x_1, x_2) \\ y_s(x_1, x_2) \end{bmatrix} = \begin{bmatrix} 100x_2(1 - x_2) \\ 0 \\ -100(1 - x_2) \end{bmatrix}, \quad (5.2)$$
<figure><img src="image_2.png" /><figcaption>Figure 1: Domain Ω and subset of control O</figcaption></figure>

and initial condition

$$ \begin{bmatrix} v_0^1(x_1, x_2) \\ v_0^2(x_1, x_2) \\ y_0(x_1, x_2) \end{bmatrix} = \begin{bmatrix} 2 \sin(2\pi x_1) \cos(2\pi x_2) \\ 2 \sin(2\pi x_1) \cos(2\pi x_2) \\ 1 - x_2 \end{bmatrix}. \qquad (5.3) $$

All our computations are realized by means of the free finite element software FreeFem++ (see [8]). The time and space approximation step are fixed to be $\Delta t = 0.01$ and $h = 0.05$

The resolution strategy consists of first proving the numerical evidence that the matrix $A_0$ is unstable. This can be confirmed by studying numerically the presence of eigenvalues with positive real part for the resulting discrete generalized eigenvalue problem

$$ \begin{pmatrix} A_{vv} & A_{vp} & A_{vy} \\ A_{zp}^T & 0 & 0 \\ A_{yv} & 0 & A_{yy} \end{pmatrix} \begin{pmatrix} v \\ p \\ y \end{pmatrix} = \lambda \begin{pmatrix} M_{vv} & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & M_{yy} \end{pmatrix} \begin{pmatrix} v \\ p \\ y \end{pmatrix}. \qquad (5.4) $$

Figure 2 illustrates the computed generalized eigenvalue problem (5.4) solved by using the ARPACK eigenvalue package. The spectrum confirms the unstable nature of the matrix $A_0$.

Our method is mainly based on the penalty method and it is therefore essential to be aware of the quality of the convergence of our results with respect to the parameter of penalization $\epsilon$.

In Figure 3 we represent the errors in a logarithm scale between the Riccati matrices $\Pi_{\epsilon}$ and $\Pi_0$ as a function of $\epsilon$. The numerical behavior is in total agreement with the result given by (3.6). The computation of the Riccati matrices are done thanks to the package “control” in Octave. We estimate a convergence rate of about one, which is in agreement with the usual results of penalty methods.

Figure 4 shows the error in a logarithm scale between the controlled solutions $(v, y)_{\epsilon}$ and $(v, y)$ while the same measurement is given for controls $u_{\epsilon}$ and $u$ in Figure 5, as a function of $\epsilon$. Once again the numerical behaviors are in total agreement with the theory. To solve the evolution problems, the time derivative is approximated thanks to a standard Euler explicit scheme. We estimate a convergence rate of about one.
<figure><img src="image_3.png" /><figcaption>Figure 2: Eigenvalues associated to (5.4)</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3: Illustration of the convergence result (3.6)</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4: Illustration of the convergence result (3.7)</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 5: Illustration of the convergence result (3.8)</figcaption></figure>

We present in the following array the exact computing values of the errors.

<table><thead><tr><th>ε</th><th>||Π<sub>ε</sub>P - Π<sub>0</sub>|| / ||Π<sub>0</sub>||</th><th>||(v, y)<sub>ε</sub> - (v, y)||</th><th>||u<sub>ε</sub> - u||</th></tr></thead><tbody><tr><td>1 e-1</td><td>9.10 e-3</td><td>1.65 e-1</td><td>5.90 e-3</td></tr><tr><td>1 e-2</td><td>9.73 e-4</td><td>1.41 e-2</td><td>7.76 e-4</td></tr><tr><td>1 e-3</td><td>9.80 e-5</td><td>1.46 e-3</td><td>8.03 e-5</td></tr><tr><td>1 e-4</td><td>9.82 e-6</td><td>1.47 e-4</td><td>8.06 e-6</td></tr></tbody></table>

**Remark 5.1.** We observe for the three previous cases a slope that is equal to one. This illustrates numerically the expected speed of convergence with respect to ε. Hence, we naturally conjecture that we actually have

$$ \|\Pi_{\varepsilon}P - \Pi_{0}\|_{\mathcal{L}(H_h; X_h)} = \mathcal{O}(\varepsilon) \quad (t \ge 0), \qquad (5.5) $$

$$ \left\| e^{t(A_{\varepsilon} - BB^{*}\Pi_{\varepsilon})} P - e^{t(A_0 - B_0 B_0^{*} \Pi_0)} \right\|_{\mathcal{L}(H_h; X_h)} = \mathcal{O}(\varepsilon) \quad (t \ge 0), \qquad (5.6) $$

$$ \left\| -B^{*}\Pi_{\varepsilon}e^{t(A_{\varepsilon}-BB^{*}\Pi_{\varepsilon})}P - (-B_{0}^{*}\Pi_{0}e^{t(A_{0}-B_{0}B_{0}^{*}\Pi_{0})}) \right\|_{\mathcal{L}(H_h; U_h)} = \mathcal{O}(\varepsilon) \quad (t \ge 0). \qquad (5.7) $$

We now give some snapshots of the evolution of the horizontal velocity and the temperature during the control process in Figures 6 and 7. We observe the stabilization of the solution. We also observe that it is in the control zone where the values of the solution are the smallest.

To better appreciate the efficiency of the method we continue our study by giving some snapshots of the evolution of the horizontal velocity and the temperature without control in Figures 8 and 9. We remark that the both velocity and temperature do not converge to 0, this is due to the instability of the spectrum of the linearized system.

We then plot the evolution of the controlled and uncontrolled solution.
<figure><img src="image_7.png" /><figcaption>Figure 6: Snapshots of the controlled horizontal velocity</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 7: Snapshots of the controlled temperature</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 8: Snapshots of the uncontrolled horizontal velocity</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 9: Snapshots of the uncontrolled temperature</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 10: Evolution of the controlled and uncontrolled horizontal velocity</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 11: Evolution of the controlled and uncontrolled temperature</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure 12: Domain Ω including another subset of control O</figcaption></figure>

We end this example by considering the effect of the control location. This consists in considering our configuration by just replacing

$$ \mathcal{O} = \{(x_1, x_2) \in \Omega ; (x_1 - 0.25)^2 + (x_2 - 0.75)^2 < 0.01\}. \quad (5.8) $$

To illustrate the effect of the control we will give hereafter We now some snapshots of the evolution of the horizontal velocity and the temperature during the control process in Figures 13 and 14. We observe the stabilization of the solution as well as that it is in the control zone where the values of the solution are the smallest.
<figure><img src="image_14.png" /><figcaption>Figure 13: Snapshots of the controlled horizontal velocity</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Figure 14: Snapshots of the controlled temperature</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Figure 15: Domain $\Omega$ with obstacle and the control zone $\mathcal{O}$</figcaption></figure>

## 5.2 Second example: Flow past a square section cylinder

We now consider a more complex configuration described by a rectangle in which we fix an obstacle. More precisely, we consider the following configuration

$$ \Omega := (0, 1) \times (0, 4) \setminus [23/16, 25/16] \times [7/16, 9/16], \qquad (5.9) $$

and

$$ \mathcal{O} := \{(x_1, x_2) \in \Omega ; (x_1 - 3/8)^2 + (x_2 - 0.5)^2 < 0.01\}. \qquad (5.10) $$

The stationary state $(v_s, y_s)$ and the initial condition $(v_0, y_0)$ are chosen as in (5.2) and (5.3), respectively.

The purpose here is to illustrate our numerical results in a quantitative way. This consists in displaying some snapshots of the evolution of the horizontal velocity and the temperature during the control process in Figures 16 and 17.
<figure><img src="image_17.png" /></figure>

$t = 0$

<figure><img src="image_18.png" /></figure>

$t = 0.5$

<figure><img src="image_19.png" /></figure>

$t = 1$

<figure><img src="image_20.png" /></figure>

$t = 1.5$

<figure><img src="image_21.png" /></figure>

$t = 2$

<figure><img src="image_22.png" /></figure>

$t = 2.5$

<caption>Figure 16: Snapshots of the controlled horizontal velocity</caption>
<figure><img src="image_23.png" /><figcaption>Figure 17: Snapshots of the controlled temperature</figcaption></figure>
References

[1] C. AIRIAU, J.-M. BUCHOT, R. K. DUBEY, M. FOURNIÉ, J.-P. RAYMOND, AND J. WELLER-CALVO, *Stabilization and best actuator location for the Navier-Stokes equations*, SIAM J. Sci. Comput., 39 (2017), pp. B993-B1020.

[2] P. BENNER, M. HEINKENSCHLOSS, J. SAAK, AND H. K. WEICHELT, *Efficient solution of large-scale algebraic riccati equations associated with index-2 daes via the inexact low-rank newton-adi method*, Applied Numerical Mathematics, 152 (2020), pp. 338-354.

[3] J. BORGGAARD, J. A. BURNS, A. SURANA, AND L. ZIETSMAN, *Control, estimation and optimization of energy efficient buildings*, in 2009 American Control Conference, IEEE, 2009, pp. 837-841.

[4] S. C. BRENNER AND L. R. SCOTT, *Mathematical Theory of Finite Element Methods*, Springer, 2008, 12008.

[5] J. A. BURNS, X. HE, AND W. HU, *Feedback stabilization of a thermal fluid system with mixed boundary control*, Comput. Math. Appl., 71 (2016), pp. 2170-2191.

[6] P. CIARLET, *The Finite Element Method for Elliptic Problems*, North Holland, 1978.

[7] V. GIRAULT AND P.-A. RAVIART, *Finite Element Methods for Navier-Stokes Equations*, Springer Series in Computational Mathematics, Springer-Verlag, 1986. Theory & applications.

[8] F. HECHT, *New development in freefem++*, J. Numer. Math., 20 (2012), pp. 251-265.

[9] M. HEINKENSCHLOSS, D. C. SORENSEN, AND K. SUN, *Balanced truncation model reduction for a class of descriptor systems with application to the Oseen equations*, SIAM J. Sci. Comput., 30 (2008), pp. 1038-1063.

[10] S. LABBÉ AND E. TRÉLAT, *Uniform controllability of semidiscrete approximations of parabolic control systems*, Systems Control Lett., 55 (2006), pp. 597-609.

[11] K. LE BALC'H AND M. TUCSNAK, *A penalty approach to the infinite horizon LQR optimal control problem for the linearized Boussinesq system*, ESAIM Control Optim. Calc. Var., 27 (2021), pp. Paper No. 17, 30.

[12] J. LIU, *Open and traction boundary conditions for the incompressible Navier-Stokes equations*, J. Comput. Phys., 228 (2009), pp. 7250-7267.

[13] M. RAMASWAMY, J.-P. RAYMOND, AND A. ROY, *Boundary feedback stabilization of the Boussinesq system with mixed boundary conditions*, Journal of Differential Equations, 266 (2019), pp. 4268-4304.

[14] J.-P. RAYMOND, *Stokes and Navier-Stokes equations with nonhomogeneous boundary conditions*, Ann. Inst. H. Poincaré Anal. Non Linéaire, 24 (2007), pp. 921-951.

[15] R. TEMAM, *Une méthode d'approximation de la solution des équations de Navier-Stokes*, Bull. Soc. Math. France, 96 (1968), pp. 115-152.

[16] E. TRÉLAT, *Contrôle optimal*, Mathématiques Concrètes. [Concrete Mathematics], Vuibert, Paris, 2005. Théorie & applications. [Theory and applications].
[17] U. VAIDYA, R. RAJARAM, AND S. DASGUPTA, *Actuator and sensor placement in linear advection PDE with building system application*, Journal of Mathematical Analysis and Applications, 394 (2012), pp. 213–224.