# Stability analysis of reaction-diffusion PDEs coupled at the boundaries with an ODE

Hugo Lhachemi, Christophe Prieur

► To cite this version:

Hugo Lhachemi, Christophe Prieur. Stability analysis of reaction-diffusion PDEs coupled at the boundaries with an ODE. Automatica, 2022, 144, pp.110465. 10.1016/j.automatica.2022.110465. hal-03726804

HAL Id: hal-03726804

https://hal.science/hal-03726804v1

Submitted on 25 Feb 2025

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.

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

Distributed under a Creative Commons Attribution 4.0 International License
# Stability analysis of reaction-diffusion PDEs coupled at the boundaries with an ODE*

Hugo Lhachemi<sup>a</sup>, Christophe Prieur<sup>b</sup>,

## Abstract

This paper addresses the derivation of generic and tractable sufficient conditions ensuring the stability of a coupled system composed of a reaction-diffusion partial differential equation (PDE), with possibly spatially-varying coefficients, and a finite-dimensional linear time invariant ordinary differential equation (ODE). The coupling of the PDE with the ODE is located at both left and right boundary conditions of the reaction-diffusion equation and takes the form of the input and output of the ODE. We investigate the four possible sets of left/right boundary couplings among Dirichlet and Neumann traces. The adopted approach relies on the spectral reduction of the problem by projecting the trajectory of the PDE into a Hilbert basis composed of the eigenvectors of the underlying Sturm-Liouville operator and yields a set of sufficient stability conditions taking the form of LMIs. We propose numerical examples, consisting of an unstable reaction-diffusion equation and an unstable ODE, such that the application of the derived stability conditions ensure the stability of the resulting coupled PDE-ODE system.

**Key words:** Coupled PDE-ODE, stability, reaction-diffusion equation, modal decomposition, LMI

## 1 Introduction

The stability analysis and control of coupled PDE-ODE systems has emerged relatively recently in the literature (and more generally PDEs with dynamical boundary conditions, see e.g. [22]). Such a trend was driven by a certain number of practical applications involving a finite-dimensional dynamics coupled to a phenomenon described by a PDE. This includes, to cite a few, solid-gas interaction of heat diffusion and chemical reaction [27], flexible cranes [14], flexible aircraft [19], drilling mechanisms [5], and power converters connected to transmission lines [12]. PDE-ODE coupling can also arise due to feedback control. Indeed, the PDE can represent the open-loop plant to be controlled while the ODE part gathers controller and actuator dynamics, see e.g. [16,26]. Conversely, the PDE can represent the dynamics of an actuator (e.g., heat or flux sensors) that is embedded into the closed-loop control of a finite-dimensional plant modeled by an ODE.

One of the very first contributions in this field was reported in [16] dealing with the state-feedback stabilization and the observer design of a diffusion PDE cascaded with an ODE via Dirichlet interconnections (see also [26] for the case of Neumann interconnections). Such a problem can be interpreted as a compensation problem of an infinite-dimensional input dynamics [17] and was solved by employing a backstepping control design procedure. This approach was also reported in the case of string equation in [15,26] and was later on applied to other types of PDEs such as beam [30] and linearized Kirtegew-de Vries [2] equations. This backstepping-based procedure for PDE-ODE cascades was then extended to other boundary stabilization problems such as wave PDEs cascaded with MIMO LTI systems [8], a diffusion PDE coupled with an ODE [27], and a diffusion PDE sandwiched between two ODEs [28]. The robustness of certain of these control strategies for the stabilization of PDE-ODE cascades were studied in [16,26,25], particularly for heat equations w.r.t. the diffusion coefficient and the length of the diffusion domain. Other extensions embracing the augmentation of the backstepping transformation with either adaptive or sliding mode control have been investigated in [20,29]. Recently, a different control design strategy using Sylvester equation was proposed in [21] for ODE-PDE and PDE-ODE cascades.

<footnote><sup>a</sup>Université Paris-Saclay, CNRS, CentraleSupélec, Laboratoire des signaux et systèmes, 91190, Gif-sur-Yvette, France</footnote>

<footnote><sup>b</sup>Université Grenoble Alpes, CNRS, Grenoble-INP, GIPSA-lab, F-38000, Grenoble, France</footnote>

<footnote>* Corresponding author H. Lhachemi. The work of C. Prieur has been partially supported by MIAI@Grenoble Alpes (ANR-19-P3IA-0003)</footnote>

<footnote>Email addresses: hugo.lhachemi@centralesupelec.fr
(Hugo Lhachemi), christophe.prieur@gipsa-lab.fr
(Christophe Prieur).</footnote>
The focus of this paper is put on the derivation of generic and tractable sufficient stability conditions ensuring the stability of coupled systems composed of a general reaction-diffusion PDE and an ODE, with various coupling configurations, rather than on the design of a particular controller for a specific setting. Such generic conditions can be used, e.g., to assess the stability of closed-loop systems with finite-dimensional controllers obtained either on a finite-dimensional truncated model of the PDE or via the approximation of infinite-dimensional output feedback controllers such as the ones described in the previous paragraph. The traditional approach for studying the stability of coupled PDE-ODE systems consists of the adequate selection of a Lyapunov functional. At a very high level, the general trend is to build the Lyapunov functional by considering terms related to 1) the energy of the PDE (measured via a relevant norm); 2) the energy of the ODE; 3) the coupling of the PDE-ODE system. Such Lyapunov functionals can be built manually [10,11,18] but can also be obtained numerically by considering very general Lyapunov functional candidates while resorting to numerical methods, such as a sum of square procedure, to obtain an admissible suitable set of parameters [1,13].

In the abovementioned context, a number of contributions have been reported in the recent years to study the stability of coupled PDE-ODE systems with couplings occurring at the boundaries of the PDE. One of the most fruitful approaches relies on the use of Legendre polynomials as a basis of projection for the PDE trajectories. In essence, this consists of the construction of a classical Lyapunov functional accounting for the PDE and ODE parts considered separately while adding a cross quadratic term mixing the state of the ODE with a finite number of coefficients of projection of the PDE trajectory into the basis of Legendre polynomials. Such an approach was reported in [7] for the study of a coupled system composed of a reaction PDE and an ODE. This method was also reported in [6] in the case of a string equation coupled with an ODE, as well as in [3] for the study of input-output stability. Input-output stability properties for coupled PDE-ODE systems using Legendre polynomials-based projections was further investigated in [4]. Recently, the stability of abstract boundary control systems with dynamic boundary conditions and positive underlying $C_0$-semigroups was studied in [9].

In this paper, we study the stability of a generic 1-D reaction diffusion equation coupled at the boundaries with a finite-dimensional ODE. Compared to [7], which was concerned with an open-loop stable constant coefficient diffusion PDE with left and right Dirichlet couplings, our approach allows the consideration of generic reaction-diffusion PDEs that are possibly open-loop unstable and with the four possible sets of left/right boundary couplings among Dirichlet and Neumann traces. The adopted approach relies on a spectral reduction-based method used to build a suitable Lyapunov functional

candidate. We obtain a set of tractable LMI conditions ensuring the exponential stability of the coupled PDE-ODE system. We also show the exponential decrease to zero of the coupling channels. The relevance of these LMI conditions are assessed based on numerical examples associated with PDEs and ODEs that are all unstable.

The rest of the paper is organized as follows. Section 2 describes the notations and reports a number of basic properties for Sturm-Liouville operators. Then the study is split into two parts. Firstly, the case of a Dirichlet trace used as an input for the ODE is investigated in Section 3. Secondly, the case of a Neumann trace used as an input for the ODE is reported in Section 4. Finally, concluding remarks are formulated in Section 5.

## 2 Notation and properties

Spaces $\mathbb{R}^n$ are endowed with the Euclidean norm denoted by $\| \cdot \|$. The associated induced norms of matrices are also denoted by $\| \cdot \|$. $L^2(0, 1)$ stands for the space of square integrable functions on $(0, 1)$ and is endowed with the inner product $\langle f, g \rangle = \int_0^1 f(x)g(x) dx$ and the norm is denoted by $\| \cdot \|_{L^2}$. For an integer $m \ge 1$, the $m$-order Sobolev space is denoted by $H^m(0, 1)$ and is endowed with its usual norm $\| \cdot \|_{H^m}$. For a symmetric matrix $P \in \mathbb{R}^{n \times n}$, $P \succeq 0$ (resp. $P > 0$) means that $P$ is positive semi-definite (resp. positive definite) while $\lambda_M(P)$ (resp. $\lambda_m(P)$) denotes its maximal (resp. minimal) eigenvalue.

Let $p \in C^1([0, 1])$ and $q \in C^0([0, 1])$ with $p > 0$ and $q \ge 0$. For $n_L, n_R \in \{0, 1\}$, let the Sturm-Liouville operator $A_{n_L, n_R} : D(A_{n_L, n_R}) \subset L^2(0, 1) \to L^2(0, 1)$ be defined by $A_{n_L, n_R}f = -(pf)' + qf$ on the domain $D(A_{n_L, n_R}) = \{f \in H^2(0, 1) : f^{(n_L)}(0) = f^{(n_R)}(1) = 0\}$. The operator $A_{n_L, n_R}$ is self-adjoint and its eigenvalues $\lambda_n, n \ge 1$, are simple, non negative, and form an increasing sequence with $\lambda_n \to +\infty$ as $n \to +\infty$. Moreover, the associated unit eigenvectors $\phi_n \in L^2(0, 1)$ form a Hilbert basis and we also have $D(A_{n_L, n_R}) = \{f \in L^2(0, 1) : \sum_{n \ge 1} |\lambda_n|^2 |\langle f, \phi_n \rangle|^2 < +\infty\}$ with $A_{n_L, n_R}f = \sum_{n \ge 1} \lambda_n \langle f, \phi_n \rangle \phi_n$.

Let $p_*, p^*, q^* \in \mathbb{R}$ be such that $0 < p_* \le p(\xi) \le p^*$ and $0 \le q(\xi) \le q^*$ for all $\xi \in [0, 1]$, then it holds [23]:

$$0 \le \pi^2(n-1)^2 p_* \le \lambda_n \le \pi^2 n^2 p^* + q^* \quad (1)$$

for all $n \ge 1$. Assuming further than $p \in C^2([0, 1])$, we have that $\phi_n(0) = O(1)$ and $\phi'_n(0) = O(\sqrt{\lambda_n})$ as $n \to +\infty$ (see [23]). We also have for all $f \in D(A_{n_L, n_R})$, $\langle A_{n_L, n_R}f, f \rangle = \sum_{n \ge 1} \lambda_n \langle f, \phi_n \rangle^2$ hence

$$\sum_{n \ge 1} \lambda_n \langle f, \phi_n \rangle^2 = \int_0^1 p(f')^2 + qf^2 d\xi. \quad (2)$$
This in particular implies that, for any $f \in D(\mathcal{A}_{n_L, n_R})$ and any $\xi \in [0, 1]$, $f(\xi) = \sum_{n \ge 1} \langle f, \phi_n \rangle \phi_n(\xi)$ and $f'(\xi) = \sum_{n \ge 1} \langle f, \phi_n \rangle \phi'_n(\xi)$. For any $\alpha \in (0, 1)$, we introduce the fractional powers of $\mathcal{A}_{n_L, n_R}$ by defining $D(\mathcal{A}_{n_L, n_R}^\alpha) = \{f \in L^2(0, 1) : \sum_{n \ge 1} |\lambda_n|^{2\alpha} |\langle f, \phi_n \rangle|^2 < +\infty\}$ and $\mathcal{A}_{n_L, n_R}^\alpha f = \sum_{n \ge 1} \lambda_n^\alpha \langle f, \phi_n \rangle \phi_n$.

For any $f \in L^2(0, 1)$ we define $\mathcal{P}_N f = f - \sum_{n=1}^N \langle f, \phi_n \rangle \phi_n = \sum_{n \ge N+1} \langle f, \phi_n \rangle \phi_n$

# 3 Dirichlet trace as an input of the ODE

## 3.1 Coupled PDE-ODE systems

We consider in this section one of the two following PDE-ODE systems:

$$z_t(t, \xi) = (pz_{\xi})_{\xi}(t, \xi) - \tilde{q}(\xi)z(t, \xi) \quad (3a)$$

$$z_{\xi}(t, 0) = 0, \quad z(t, 1) = y(t) = Cx(t) \quad (3b)$$

$$\dot{x}(t) = Ax(t) + Bz(t, 0) \quad (3c)$$

$$z(0, \xi) = z_0(\xi), \quad x(0) = x_0 \quad (3d)$$

or

$$z_t(t, \xi) = (pz_{\xi})_{\xi}(t, \xi) - \tilde{q}(\xi)z(t, \xi) \quad (4a)$$

$$z_{\xi}(t, 0) = 0, \quad z_{\xi}(t, 1) = y(t) = Cx(t) \quad (4b)$$

$$\dot{x}(t) = Ax(t) + Bz(t, 0) \quad (4c)$$

$$z(0, \xi) = z_0(\xi), \quad x(0) = x_0 \quad (4d)$$

for $t > 0$ and $\xi \in (0, 1)$. Here $p \in \mathcal{C}^2([0, 1])$ with $p > 0$, $\tilde{q} \in \mathcal{C}^0([0, 1])$, $A \in \mathbb{R}^{n \times n}$, $B \in \mathbb{R}^n$, and $C \in \mathbb{R}^{n \times 1}$ are matrices, $z_0 \in L^2(0, 1)$ and $x_0 \in \mathbb{R}^n$ are initial conditions, and $z(t, \cdot) \in L^2(0, 1)$ and $x(t) \in \mathbb{R}^n$ are the state of the reaction-diffusion PDE and of the ODE at time $t$, respectively. We introduce without loss of generality $q \in \mathcal{C}^0([0, 1])$ and $q_c \in \mathbb{R}$ all such that

$$\tilde{q} = q - q_c, \quad q \ge 0. \quad (5)$$

In the case (4), without loss of generality, we further constrain $q$ such that $q > 0$.

**Remark 1** Systems (3) and (4), as well as systems (18) and (19) that will be described in the next section, can be used to represent a variety of practical situations. For instance the PDE part can stand for a reaction-diffusion process coupled with a finite-dimensional LTI controller materialized by the ODE. Conversely, the ODE part can merge both finite-dimensional LTI plant along with its associated finite-dimensional LTI controller while the PDE part describes the sensor dynamics. This latter situation is similar to the one described in [16] where a controller was designed for a cascaded PDE-ODE system using a backstepping transformation. One of the main motivations for deriving generic stability conditions for

coupled PDE-ODE systems such as (3) and (4) is ultimately when the to-be-implemented finite-dimensional controller is computed either on a finite-dimensional approximation of the PDE or via the approximation of an infinite-dimensional output feedback controller (obtained, e.g., using backstepping control design procedures).

**Remark 2** The PDE-ODE system (3) was studied in [7] in the case of a stable diffusion PDE (i.e., without reaction term) with a constant diffusion coefficient using the projection of the PDE trajectories into a finite subset of Legendre polynomials. In this paper, we propose a spectral reduction-based framework for handling general reaction-diffusion PDEs, possibly unstable with spatially varying coefficient, and for various combinations of couplings selected among Dirichlet and Neumann traces.

## 3.2 Preliminary spectral reduction

We rewrite (3) and (4) under an equivalent PDE-ODE system with homogeneous boundary conditions. Specifically, introducing for (3) the change of variable

$$w(t, \xi) = z(t, \xi) - \xi^2 y(t) \quad (6)$$

we infer that (3) is equivalent to

$$w_t(t, \xi) = (pz_{\xi})_{\xi}(t, \xi) + (q_c - q(\xi))w(t, \xi) \quad (7a) \\ + a(\xi)y(t) + b(\xi)\dot{y}(t)$$

$$w_{\xi}(t, 0) = w(t, 1) = 0, \quad y(t) = Cx(t) \quad (7b)$$

$$\dot{x}(t) = Ax(t) + Bw(t, 0) \quad (7c)$$

$$w(0, \xi) = w_0(\xi), \quad x(0) = x_0 \quad (7d)$$

with $a(\xi) = 2p(\xi) + 2\xi p'(\xi) + (q_c - q(\xi))\xi^2$, $b(\xi) = -\xi^2$, and $w_0(\xi) = z_0(\xi) - \xi^2 y(0)$. Similarly, we introduce for (4) the change of variable

$$w(t, \xi) = z(t, \xi) - \frac{\xi^2}{2} y(t), \quad (8)$$

with which we obtain that (4) is equivalent to

$$w_t(t, \xi) = (pz_{\xi})_{\xi}(t, \xi) + (q_c - q(\xi))w(t, \xi) \quad (9a) \\ + a(\xi)y(t) + b(\xi)\dot{y}(t)$$

$$w_{\xi}(t, 0) = w_{\xi}(t, 1) = 0, \quad y(t) = Cx(t) \quad (9b)$$

$$\dot{x}(t) = Ax(t) + Bw(t, 0) \quad (9c)$$

$$w(0, \xi) = w_0(\xi), \quad x(0) = x_0 \quad (9d)$$

with $a(\xi) = p(\xi) + \xi p'(\xi) + (q_c - q(\xi))\frac{\xi^2}{2}$, $b(\xi) = -\frac{\xi^2}{2}$, and $w_0(\xi) = z_0(\xi) - \frac{\xi^2}{2}y(0)$. After these change of variable, the well-posedness in terms of classical solutions of the above PDE-ODE systems for initial conditions $w_0 \in D(\mathcal{A}_{1, n_R}^{1/2})$ and $x_0 \in \mathbb{R}^n$ is a consequence of [24, Chap. 6.3]. More precisely, setting
$n_R = 0$ in the case (7) and $n_R = 1$ in the case (9), we have for any $w_0 \in D(A_{1,n_R}^{1/2})$ and any $x_0 \in \mathbb{R}^n$ the existence and uniqueness of a classical solution $(w, x) \in C^0(\mathbb{R}_{\ge 0}; L^2(0, 1) \times \mathbb{R}^n) \cap C^1(\mathbb{R}_{>0}; L^2(0, 1) \times \mathbb{R}^n)$ with $w(t, \cdot) \in D(A_{1,n_R})$ for all $t > 0$. Moreover, from the proof of [24, Thm. 6.3.1], we have $A_{1,n_R} w \in C^0(\mathbb{R}_{>0}; L^2(0, 1))$ and $A_{1,n_R}^{1/2} w \in C^0(\mathbb{R}_{\ge 0}; L^2(0, 1))$.

We now introduce in the case (7) [resp. (9)] the Hilbert basis $\{\phi_i : i \ge 1\}$ of $L^2(0, 1)$ formed by eigenvectors of the Sturm-Liouville operator $A_{1,0}$ [resp. $A_{1,1}$]. From now on, the treatment of both (7) and (9) is identical. We define the coefficients of projection:

$$w_i(t) = \langle w(t, \cdot), \phi_i \rangle, \quad a_i = \langle a, \phi_i \rangle, \quad b_i = \langle b, \phi_i \rangle$$

for $i \ge 1$. Considering classical solutions, we obtain that

$$\dot{w}_i(t) = (-\lambda_i + q_c)w_i(t) + a_i C x(t) + b_i C A x(t) \quad (10a) \\ + b_i C \left\{ Ax(t) + B \sum_{j \ge 1} w_j(t) \phi_j(0) \right\}$$

$$\dot{x}(t) = Ax(t) + B \sum_{j \ge 1} w_j(t) \phi_j(0) \quad (10b)$$

for $i \ge 1$. For an arbitrary integer $N \ge 1$, we define

$$W(t) = [w_1(t) \dots w_N(t)]^\top \in \mathbb{R}^N, \\ A_N = \operatorname{diag}(-\lambda_1 + q_c, \dots, -\lambda_N + q_c) \in \mathbb{R}^{N \times N}, \\ B_{a,N} = [a_1 \dots a_N]^\top \in \mathbb{R}^N, \\ B_{b,N} = [b_1 \dots b_N]^\top \in \mathbb{R}^N, \\ C_N = [\phi_1(0) \dots \phi_N(0)] \in \mathbb{R}^{1 \times N},$$

and $\mathcal{R}(t) = \sum_{i \ge N+1} w_i(t)\phi_i(0)$. We infer from (10a) that

$$\dot{W}(t) = (A_N + B_{b,N}C B C_N)W(t) + (B_{a,N}C + B_{b,N}C A)x(t) + B_{b,N}CB\mathcal{R}(t).$$

Combining this latter identity with (10b) while defining

$$X(t) = \begin{bmatrix} W(t) \\ x(t) \end{bmatrix} \in \mathbb{R}^{N+n},$$

we infer that

$$\dot{X}(t) = F X(t) + G \mathcal{R}(t) \quad (12)$$

where

$$F = \begin{bmatrix} A_N + B_{b,N}C B C_N & B_{a,N}C + B_{b,N}C A \\ B C_N & A \end{bmatrix}$$

and

$$G = \begin{bmatrix} B_{b,N} C B \\ B \end{bmatrix}.$$

The residual dynamics, which corresponds to the modes $i \ge N+1$, is characterized by

$$\dot{w}_i(t) = (-\lambda_i + q_c)w_i(t) + a_i C x(t) + b_i C A x(t) \quad (13) \\ + b_i C B C_N W(t) + b_i C B \mathcal{R}(t).$$

We finally define $M_{1,\phi} = \sum_{i \ge N+1} \frac{\phi_i(0)^2}{\lambda_i} < +\infty$ and

$$H = \begin{bmatrix} H_{1,1} & 0 \\ 0 & H_{2,2} \end{bmatrix}$$

with $H_{1,1} = \|P_N b\|_{L^2}^2 C_N^\top B^\top C^\top C B C_N$ and $H_{2,2} = \|P_N a\|_{L^2}^2 C^\top C + \|P_N b\|_{L^2}^2 A^\top C^\top C A$.

### 3.3 Main result

We can now introduce the main result of this section.

**Theorem 3** Let $p \in C^2([0, 1])$ with $p > 0$, $\tilde{q} \in C^0([0, 1])$, $A \in \mathbb{R}^{n \times n}$, $B \in \mathbb{R}^n$, and $C \in \mathbb{R}^{n \times 1}$ be given. Let $q \in C^0([0, 1])$ and $q_c \in \mathbb{R}$ be such that (5) holds with $q \ge 0$ in the case (3) or $q > 0$ in the case (4). Assume that there exist $N \ge 1$, $P > 0$, $\alpha > 2$, and $\beta > 0$ such that $\Theta_1 < 0$ and $\Theta_2 < 0$ where

$$\Theta_1 = \begin{bmatrix} F^\top P + P F + \alpha H & PG \\ G^\top P & \alpha \|P_N b\|_{L^2}^2 (CB)^2 - \beta \end{bmatrix},$$

$$\Theta_2 = \begin{bmatrix} -\lambda_{N+1} + q_c + \frac{\beta M_{1,\phi}}{2} & \sqrt{2\lambda_{N+1}} \\ \sqrt{2\lambda_{N+1}} & -\alpha \end{bmatrix},$$

then there exist constants $\eta$, $M > 0$ such that, for any initial conditions $z_0 \in H^2(0, 1)$ and $x_0 \in \mathbb{R}^n$ such that $z'_0(0) = 0$ and $z_0(1) = C x_0$ [resp. $z'_0(0) = 0$ and $z'_0(1) = C x_0$], the classical solution of (3) [resp. (4)] satisfies

$$\|z(t, \cdot)\|_{H^1}^2 + \|x(t)\|^2 \le Me^{-2\eta t} (\|z_0\|_{H^1}^2 + \|x_0\|^2) \quad (14)$$

with coupling channels such that

* for (3): $z(t, 0)^2 + z(t, 1)^2 \le Me^{-2\eta t} (\|z_0\|_{H^1}^2 + \|x_0\|^2)$;

* for (4): $z(t, 0)^2 + z_\xi(t, 1)^2 \le Me^{-2\eta t} (\|z_0\|_{H^1}^2 + \|x_0\|^2)$

for all $t > 0$.
**Remark 4** The conclusions of Theorem 3 actually hold for any initial conditions such that $w_0 \in D(\mathcal{A}_{1,n_R}^{1/2})$. In the case (3), it is easy to show based on (2) that $D(\mathcal{A}_{1,0}^{1/2}) = \{f \in H^1(0,1) : f(1) = 0\}$. In this setting the conclusions of Theorem 3 hold for any $z_0 \in H^1(0,1)$ and $x_0 \in \mathbb{R}^n$ such that $z_0(1) = Cx_0$. Similarly, since $D(\mathcal{A}_{1,1}^{1/2}) = H^1(0,1)$, the conclusions of Theorem 3 in the case (4) hold for any $z_0 \in H^1(0,1)$ and $x_0 \in \mathbb{R}^n$.

**Proof.** Let $N \ge 1$, $P > 0$, $\alpha > 2$, and $\beta > 0$ such that $\Theta_1 < 0$ and $\Theta_2 < 0$. Hence, there exist $\eta > 0$ such that $\Theta_{1,\eta} \le 0$ and $\Theta_{2,\eta} \le 0$ where

$$ 
\Theta_{1,\eta} = \begin{bmatrix} F^\top P + PF + 2\eta P + \alpha H & PG \\ G^\top P & \alpha \|P_N b\|_{L^2}^2 (CB)^2 - \beta \end{bmatrix}, 
$$
$$ 
\Theta_{2,\eta} = \begin{bmatrix} -\lambda_{N+1} + q_c + \eta + \frac{\beta M_{1,\phi}}{2} \sqrt{2\lambda_{N+1}} & -\alpha \end{bmatrix}. 
$$

Define the Lyapunov functional candidate

$$ V(X, w) = X^\top P X + \sum_{i \ge N+1} \lambda_i \langle w, \phi_i \rangle^2 \quad (15) $$

with $X \in \mathbb{R}^{N+n}$ and $w \in D(\mathcal{A}_{1,n_R})$ where $n_R = 0$ in the case (7) while $n_R = 1$ in the case (9). With the slight abuse of notation $V(t) = V(X(t), w(t))$, the computation of the time derivative along the system trajectories (12) and (13) yields for $t > 0$

$$ 
\begin{aligned}
\dot{V}(t) &= 2X(t)^\top P \dot{X}(t) + 2 \sum_{i \ge N+1} \lambda_i w_i(t) \dot{w}_i(t) \\
&= X(t)^\top (F^\top P + PF) X(t) + 2X(t)^\top P \mathcal{R}(t) \\
&\quad + 2 \sum_{i \ge N+1} \lambda_i (-\lambda_i + q_c) w_i(t)^2 \\
&\quad + 2 \sum_{i \ge N+1} \lambda_i w_i(t) a_i C x(t) + 2 \sum_{i \ge N+1} \lambda_i w_i(t) b_i C A x(t) \\
&\quad + 2 \sum_{i \ge N+1} \lambda_i w_i(t) b_i C B C_N W(t) \\
&\quad + 2 \sum_{i \ge N+1} \lambda_i w_i(t) b_i C B \mathcal{R}(t).
\end{aligned}
 $$

We estimate the four latter series by using Young's inequality. For instance, the first term is estimated as

$$ 
\begin{aligned}
& 2 \sum_{i \ge N+1} \lambda_i w_i(t) a_i C x(t) \\
& \le \sum_{i \ge N+1} \left\{ \frac{1}{\alpha} \lambda_i^2 w_i(t)^2 + \alpha a_i^2 (Cx(t))^2 \right\} \\
& \le \frac{1}{\alpha} \sum_{i \ge N+1} \lambda_i^2 w_i(t)^2 + \alpha \|P_N a\|_{L^2}^2 x(t)^\top C^\top C x(t).
\end{aligned}
 $$

Similarly, we obtain that

$$ 
\begin{aligned}
& 2 \sum_{i \ge N+1} \lambda_i w_i(t) b_i C A x(t) \\
& \le \frac{1}{\alpha} \sum_{i \ge N+1} \lambda_i^2 w_i(t)^2 + \alpha \|P_N b\|_{L^2}^2 x(t)^\top A^\top C^\top C A x(t),
\end{aligned}
 $$

$$ 
\begin{aligned}
& 2 \sum_{i \ge N+1} \lambda_i w_i(t) b_i C B C_N W(t) \le \frac{1}{\alpha} \sum_{i \ge N+1} \lambda_i^2 w_i(t)^2 \\
& \quad + \alpha \|P_N b\|_{L^2}^2 W(t)^\top C_N^\top B^\top C^\top C B C_N W(t),
\end{aligned}
 $$

and

$$ 
\begin{aligned}
& 2 \sum_{i \ge N+1} \lambda_i w_i(t) b_i C B \mathcal{R}(t) \\
& \le \frac{1}{\alpha} \sum_{i \ge N+1} \lambda_i^2 w_i(t)^2 + \alpha \|P_N b\|_{L^2}^2 (CB)^2 \mathcal{R}(t)^2.
\end{aligned}
 $$

The use of the four latter estimates implies that

$$ 
\begin{aligned}
\dot{V}(t) \le & \begin{bmatrix} X(t) \\ \mathcal{R}(t) \end{bmatrix}^\top \begin{bmatrix} F^\top P + PF + \alpha H & PG \\ G^\top P & \alpha \|P_N b\|_{L^2}^2 (CB)^2 \end{bmatrix} \begin{bmatrix} X(t) \\ \mathcal{R}(t) \end{bmatrix} \\
& + 2 \sum_{i \ge N+1} \lambda_i \left( -\lambda_i + q_c + \frac{2\lambda_i}{\alpha} \right) w_i(t)^2
\end{aligned}
\quad (16) $$

for $t > 0$. Since $\mathcal{R}(t) = \sum_{i \ge N+1} w_i(t) \phi_i(0)$, we infer that $\mathcal{R}(t)^2 \le M_{1,\phi} \sum_{i \ge N+1} \lambda_i w_i(t)^2$. This implies for $t > 0$ that

$$ 
\begin{aligned}
\dot{V}(t) + 2\eta V(t) & \le \begin{bmatrix} X(t) \\ \mathcal{R}(t) \end{bmatrix}^\top \Theta_{1,\eta} \begin{bmatrix} X(t) \\ \mathcal{R}(t) \end{bmatrix} + 2 \sum_{i \ge N+1} \lambda_i \Gamma_i w_i(t)^2
\end{aligned}
\quad (17) $$

where $\Gamma_i = - (1 - \frac{2}{\alpha}) \lambda_i + q_c + \eta + \frac{\beta M_{1,\phi}}{2}$. Now, since $\alpha > 2$, we have $\Gamma_i \le \Gamma_{N+1}$ for all $i \ge N+1$. Moreover, combining $\Theta_{2,\eta} \le 0$ and the Schur complement, we infer that $\Gamma_{N+1} \le 0$. Using also $\Theta_{1,\eta} \le 0$, we obtain that $\dot{V}(t) + 2\eta V(t) \le 0$ for all $t > 0$. Since $\mathcal{A}_{1,n_R}^{1/2} w \in C^0(\mathbb{R}_{\ge 0}; L^2(0,1))$, the mapping $t \mapsto V(t)$ is continuous for $t \ge 0$, implying that $V(t) \le e^{-2\eta t} V(0)$ for all $t \ge 0$. We now note from (15) that $\dot{V}(0) \le \lambda_M(P) \|X(0)\|^2 + \sum_{i \ge N+1} \lambda_i \langle w_0, \phi_i \rangle^2$. Noting on one side that $\|X(0)\|^2 \le \|x_0\|^2 + \|w_0\|_{L^2}^2$, while on the other side, from (2), $\sum_{i \ge N+1} \lambda_i \langle w_0, \phi_i \rangle^2 \le \max(\sup |p|, \sup |q|) \|w_0\|_{H^1}^2$, we infer the existence of a constant $M_1 > 0$ such that $V(0) \le M_1 (\|x_0\|^2 + \|w_0\|_{H^1}^2)$. We now consider separately the cases (3) and (4). In the case
(3), since $w(t, 1) = 0$, Poincaré's inequality along with (2) and (15) shows the existence of a constant $M_2 > 0$ such that $\|w(t, \cdot)\|_{H^1}^2 \le M_2 V(t)$. In the case (4), the extra assumption $q > 0$ directly implies from (2) and (15) the existence of a constant $M'_2 > 0$ such that $\|w(t, \cdot)\|_{H^1}^2 \le M'_2 V(t)$. In both cases, we infer the existence of a constant $M_3 > 0$ such that $\|w(t, \cdot)\|_{H^1}^2 + \|x(t)\|^2 \le M_3 (\|x_0\|^2 + \|w_0\|_{H^1}^2)$. The claimed conclusion follows from the change of variable given by either (6) or (8) and the continuous embedding $H^1(0, 1) \subset L^\infty([0, 1])$. □

From the above proof we deduce the following corollary.

**Corollary 5** In the context of Theorem 3, the decay rate $\eta > 0$ of the stability estimate (14) is guaranteed provided the LMI conditions $\Theta_{1,\eta} \preceq 0$ and $\Theta_{2,\eta} \preceq 0$ are feasible.

## 3.4 Numerical illustration

We illustrate the results of Theorem 3 and Corollary 5 for the coupled PDE-ODE system described by (3) with $p = 1, q = 0, q_c = 3$,

$$A = \begin{bmatrix} 0 & -1/4 & -1/5 & 1/5 & 1/6 \\ 1/2 & 1 & -4 & 9/2 & 7/2 \\ -9/4 & -1/2 & -14 & 23 & 16 \\ -1/5 & -1/2 & -11/4 & 1/10 & 5/4 \\ -4/3 & -4/3 & -9 & 9 & 5/2 \end{bmatrix},$$

$$B = [-7/2 \ -3/2 \ -1/10 \ 1/2 \ 1]^T,$$

$$C = [1/10 \ -1/3 \ -4 \ 7/8 \ 7/8].$$

In this case, both PDE and ODE systems are open-loop unstable. Indeed, the dominant eigenvalue of the PDE is located approximately at +0.533 while the matrix A has two unstable eigenvalues located approximately at +1.046 and +0.247. The application of Theorem 3 with $N = 6$ shows the exponential stability of the coupled PDE-ODE system (3). Moreover, the application of Corollary 5 with $N = 30$ shows the exponential stability of the coupled PDE-ODE system with decay rate $\eta = 0.47$. We illustrate this result with a numerical simulation. The numerical scheme consists in the modal approximation of the PDE plant by its 100 dominant modes. The initial condition is set as $w_0(\xi) = -1 + \xi^2$ and $x_0 = [-2 \ 1 \ 2 \ 1 \ 3]^T$. The obtained results are depicted on Fig. 1, confirming the theoretical predictions of Theorem 3 and Corollary 5.

<figure><img src="image_2.png" /><figcaption>(a) State of the reaction-diffusion system $z(t, \xi)$</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>(b) State of the ODE $x(t)$</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>(c) Coupling channels</figcaption></figure>

<caption>Fig. 1. Time evolution of the coupled PDE-ODE system (3)</caption>

# 4 Neumann trace as an input of the ODE

## 4.1 Coupled PDE-ODE systems

We consider in this section one of the two below PDE-ODE systems. The difference comparing to the previously studied coupled systems (3) and (4) is that the PDE enters into the ODE by means of a Neumann trace instead of a Dirichlet trace.

$$z_t(t, \xi) = (pz_{\xi})_{\xi}(t, \xi) - \tilde{q}(\xi)z(t, \xi) \quad (18a)$$

$$z(t, 0) = 0, \quad z(t, 1) = y(t) = Cx(t) \quad (18b)$$

$$\dot{x}(t) = Ax(t) + Bz_{\xi}(t, 0) \quad (18c)$$

$$z(0, \xi) = z_0(\xi), \quad x(0) = x_0 \quad (18d)$$

or

$$z_t(t, \xi) = (pz_{\xi})_{\xi}(t, \xi) - \tilde{q}(\xi)z(t, \xi) \quad (19a)$$

$$z(t, 0) = 0, \quad z_{\xi}(t, 1) = y(t) = Cx(t) \quad (19b)$$

$$\dot{x}(t) = Ax(t) + Bz_{\xi}(t, 0) \quad (19c)$$

$$z(0, \xi) = z_0(\xi), \quad x(0) = x_0 \quad (19d)$$

for $t > 0$ and $\xi \in (0, 1)$. Here $p \in C^2([0, 1])$ with $p > 0$, $\tilde{q} \in C^0([0, 1])$, $A \in \mathbb{R}^{n \times n}$, $B \in \mathbb{R}^n$, and $C \in \mathbb{R}^{n \times 1}$ are
matrices, $z_0 \in L^2(0,1)$ and $x_0 \in \mathbb{R}^n$ are initial conditions, and $z(t, \cdot) \in L^2(0,1)$ and $x(t) \in \mathbb{R}^n$ are the state of the reaction-diffusion PDE and of the ODE at time $t$, respectively. As in the previous section, we introduce without loss of generality a function $q \in C^0([0, 1])$ and a constant $q_c \in \mathbb{R}$ all such that (5) holds.

## 4.2 Preliminary spectral reduction

We transform (18) and (19) into an equivalent PDE-ODE system with homogeneous boundary conditions. Considering the change of variable

$$w(t, \xi) = z(t, \xi) - \xi y(t) \quad (20)$$

we infer that (18) is equivalent to

$$w_t(t, \xi) = (pz_{\xi})_{\xi}(t, \xi) + (q_c - q(\xi))w(t, \xi) \quad (21a) \\ + a(\xi)y(t) + b(\xi)\dot{y}(t)$$

$$w(t, 0) = w(t, 1) = 0, \quad y(t) = Cx(t) \quad (21b)$$

$$\dot{x}(t) = Ax(t) + B(w_{\xi}(t, 0) + y(t)) \quad (21c)$$

$$w(0, \xi) = w_0(\xi), \quad x(0) = x_0 \quad (21d)$$

while, for the same change of variable (20), (19) is equivalent to

$$w_t(t, \xi) = (pz_{\xi})_{\xi}(t, \xi) + (q_c - q(\xi))w(t, \xi) \quad (22a) \\ + a(\xi)y(t) + b(\xi)\dot{y}(t)$$

$$w(t, 0) = w_{\xi}(t, 1) = 0, \quad y(t) = Cx(t) \quad (22b)$$

$$\dot{x}(t) = Ax(t) + B(w_{\xi}(t, 0) + y(t)) \quad (22c)$$

$$w(0, \xi) = w_0(\xi), \quad x(0) = x_0 \quad (22d)$$

both with $a(\xi) = p'(\xi) + (q_c - q(\xi))\xi$, $b(\xi) = -\xi$, and $w_0(\xi) = z_0(\xi) - \xi y(0)$. Note that, after these change of variables, the well-posedness in terms of classical solutions of the above PDE-ODE systems for initial conditions $w_0 \in \cup_{\alpha_0 \in (3/4, 1)} D(A_{0,n_R}^{\alpha_0})$ and $x_0 \in \mathbb{R}^n$ is a consequence of [24, Chap. 6.3]. More precisely, for a given $\alpha_0 \in (3/4, 1)$ and setting $n_R = 0$ in the case (21) and $n_R = 1$ in the case (22), we have for any $w_0 \in D(A_{0,n_R}^{\alpha_0})$ and any $x_0 \in \mathbb{R}^n$ the existence and uniqueness of a classical solution $(w, x) \in C^0(\mathbb{R}_{\ge 0}; L^2(0, 1) \times \mathbb{R}^n) \cap C^1(\mathbb{R}_{>0}; L^2(0, 1) \times \mathbb{R}^n)$ with $w(t, \cdot) \in D(A_{0,n_R})$ for all $t > 0$. Moreover, from the proof of [24, Thm. 6.3.1], we have $A_{0,n_R} w \in C^0(\mathbb{R}_{>0}; L^2(0, 1))$ and $A_{0,n_R}^{\alpha_0} w \in C^0(\mathbb{R}_{\ge 0}; L^2(0, 1))$ hence $A_{0,n_R}^{1/2} w \in C^0(\mathbb{R}_{\ge 0}; L^2(0, 1))$.

We introduce $\{\phi_i : i \ge 1\}$ a Hilbert basis of $L^2(0, 1)$ composed of eigenvectors of either $A_{0,0}$ in the case (21) or $A_{0,1}$ in the case (22). We now treat the two cases simultaneously. To do so, we define the coefficients of projection:

$$w_i(t) = \langle w(t, \cdot), \phi_i \rangle, \quad a_i = \langle a, \phi_i \rangle, \quad b_i = \langle b, \phi_i \rangle$$

for $i \ge 1$. Considering classical solutions, we obtain that

$$\dot{w}_i(t) = (-\lambda_i + q_c)w_i(t) + a_i C x(t) \quad (23a) \\ + b_i C \left\{ (A + BC)x(t) + B \sum_{j \ge 1} w_j(t) \phi_j'(0) \right\}$$

$$\dot{x}(t) = (A + BC)x(t) + B \sum_{j \ge 1} w_j(t) \phi_j'(0) \quad (23b)$$

for $i \ge 1$. For an arbitrary integer $N \ge 1$, we define

$$W(t) = [w_1(t) \dots w_N(t)]^\top \in \mathbb{R}^N, \\ A_N = \operatorname{diag}(-\lambda_1 + q_c, \dots, -\lambda_N + q_c) \in \mathbb{R}^{N \times N},$$

$$B_{a,N} = [a_1 \dots a_N]^\top \in \mathbb{R}^N,$$

$$B_{b,N} = [b_1 \dots b_N]^\top \in \mathbb{R}^N,$$

$$C_N = [\phi_1'(0) \dots \phi_N'(0)] \in \mathbb{R}^{1 \times N},$$

and $\mathcal{R}(t) = \sum_{i \ge N+1} w_i(t) \phi_i'(0)$. We infer from (23a) that

$$\dot{W}(t) = (A_N + B_{b,N}C B C_N)W(t) + (B_{a,N}C + B_{b,N}C(A + BC))x(t) + B_{b,N}CB\mathcal{R}(t).$$

Combining this latter identity with (23b) while defining

$$X(t) = \begin{bmatrix} W(t) \\ x(t) \end{bmatrix} \in \mathbb{R}^{N+n},$$

we infer that

$$\dot{X}(t) = FX(t) + GR(t) \quad (25)$$

where

$$F = \begin{bmatrix} A_N + B_{b,N}C B C_N & B_{a,N}C + B_{b,N}C(A + BC) \\ BC_N & A + BC \end{bmatrix}$$

and

$$G = \begin{bmatrix} B_{b,N}CB \\ B \end{bmatrix}.$$

The residual dynamics, which corresponds to the modes $i \ge N + 1$, is characterized by

$$\dot{w}_i(t) = (-\lambda_i + q_c)w_i(t) + a_i C x(t) \quad (26) \\ + b_i C(A + BC)x(t) + b_i C B C_N W(t) + b_i C B \mathcal{R}(t).$$

We finally define for any $\epsilon \in (0, 1/2]$ the constant
$$M_{2,\phi}(\epsilon) = \sum_{i \ge N+1} \frac{\phi_i'(0)^2}{\lambda_i^{3/2+\epsilon}} < +\infty \text{ and the matrix}$$

$$H = \begin{bmatrix} H_{1,1} & 0 \\ 0 & H_{2,2} \end{bmatrix}$$

with $H_{1,1} = \|\mathcal{P}_N b\|_{L^2}^2 C_N^\top B^\top C^\top C B C_N$ and $H_{2,2} = \|\mathcal{P}_N a\|_{L^2}^2 C^\top C + \|\mathcal{P}_N b\|_{L^2}^2 (A + BC)^\top C^\top C(A + BC)$.

### 4.3 Main result

We can now introduce the main result of this section.

**Theorem 6** Let $p \in C^2([0, 1])$ with $p > 0$, $\tilde{q} \in C^0([0, 1])$, $A \in \mathbb{R}^{n \times n}$, $B \in \mathbb{R}^n$, and $C \in \mathbb{R}^{n \times 1}$ be given. Let $q \in C^0([0, 1])$ and $q_c \in \mathbb{R}$ be such that (5) holds. Assume that there exist $N \ge 1$, $\epsilon \in (0, 1/2]$, $P > 0$, $\alpha > 2$, and $\beta > 0$ such that $\Theta_1 < 0$, $\Theta_2 < 0$, and $\Theta_3 > 0$ where

$$\Theta_1 = \begin{bmatrix} F^\top P + PF + \alpha H & PG \\ G^\top P & \alpha \|\mathcal{P}_N b\|_{L^2}^2 (CB)^2 - \beta \end{bmatrix},$$

$$\Theta_2 = \begin{bmatrix} -\lambda_{N+1} + q_c + \frac{\beta M_{2,\phi}(\epsilon)}{2} \lambda_{N+1}^{1/2+\epsilon} \sqrt{2\lambda_{N+1}} & -\alpha \\ \sqrt{2\lambda_{N+1}} & \alpha \end{bmatrix},$$

$$\Theta_3 = \begin{bmatrix} 1 - \frac{\beta M_{2,\phi}(\epsilon)}{2\lambda_{N+1}^{1/2-\epsilon}} \sqrt{2} & \alpha \end{bmatrix},$$

then there exist constants $\eta$, $M > 0$ such that, for any initial conditions $z_0 \in H^2(0, 1)$ and $x_0 \in \mathbb{R}^n$ such that $z_0(0) = 0$ and $z_0(1) = C x_0$ [resp. $z_0(0) = 0$ and $z_0'(1) = C x_0$], the classical solution of (18) [resp. (19)] satisfies

$$\|z(t, \cdot)\|_{H^1}^2 + \|x(t)\|^2 \le Me^{-2\eta t} (\|z_0\|_{H^1}^2 + \|x_0\|^2) \quad (27)$$

with coupling channels such that

* for (18): $z_\xi(t, 0)^2 + z(t, 1)^2 \le Me^{-2\eta t} (\|A_{0,0} w_0\|_{L^2}^2 + \|z_0\|_{H^1}^2 + \|x_0\|^2)$;

* for (19): $z_\xi(t, 0)^2 + z_\xi(t, 1)^2 \le Me^{-2\eta t} (\|A_{0,1} w_0\|_{L^2}^2 + \|z_0\|_{H^1}^2 + \|x_0\|^2)$

for all $t > 0$.

**Remark 7** For any fixed $\alpha_0 \in (3/4, 1)$, the estimate (27) of Theorem 6 actually hold for any initial conditions such that $w_0 \in D(\mathcal{A}_{0,n_R}^{\alpha_0})$. In this case, the estimate regarding the coupling channels also hold when replacing $\|A_{0,n_R} w_0\|_{L^2}^2$ by $\|A_{0,n_R}^{\alpha_0} w_0\|_{L^2}^2$.

**Proof.** Let $N \ge 1$, $\epsilon \in (0, 1/2]$, $P > 0$, $\alpha > 2$, and $\beta > 0$ such that $\Theta_1 < 0$, $\Theta_2 < 0$, and $\Theta_3 > 0$. Hence, there

exist $\eta > 0$ such that $\Theta_{1,\eta} \le 0$ and $\Theta_{2,\eta} \le 0$ where

$$\Theta_{1,\eta} = \begin{bmatrix} F^\top P + PF + 2\eta P + \alpha H & PG \\ G^\top P & \alpha \|\mathcal{P}_N b\|_{L^2}^2 (CB)^2 - \beta \end{bmatrix},$$

$$\Theta_{2,\eta} = \begin{bmatrix} -\lambda_{N+1} + q_c + \eta + \frac{\beta M_{2,\phi}(\epsilon)}{2} \lambda_{N+1}^{1/2+\epsilon} \sqrt{2\lambda_{N+1}} & -\alpha \\ \sqrt{2\lambda_{N+1}} & \alpha \end{bmatrix}.$$

Considering the Lyapunov functional candidate (15) with $X \in \mathbb{R}^{N+n}$ and $w \in D(\mathcal{A}_{0,n_R})$ where $n_R = 0$ in the case (21) while $n_R = 1$ in the case (22) and adopting the same approach as the one reported in the previous section, the computation of the time derivative along the system trajectories (25) and (26) yields (16) for $t > 0$. Since $\mathcal{R}(t) = \sum_{i \ge N+1} w_i(t) \phi_i'(0)$, we infer that $\mathcal{R}(t)^2 \le M_{2,\phi}(\epsilon) \sum_{i \ge N+1} \lambda_i^{3/2+\epsilon} w_i(t)^2$. This implies that (17) holds for $t > 0$ with $\Gamma_i = -(1 - \frac{2}{\alpha}) \lambda_i + q_c + \eta + \frac{\beta M_{2,\phi}(\epsilon)}{2} \lambda_i^{1/2+\epsilon}$. Since $\epsilon \in (0, 1/2]$, we observe for $i \ge N+1$ that $\lambda_i^{1/2+\epsilon} = \lambda_i/\lambda_i^{1/2-\epsilon} \le \lambda_i/\lambda_{N+1}^{1/2-\epsilon}$ hence $\Gamma_i \le -\left(1 - \frac{2}{\alpha} - \frac{\beta M_{2,\phi}(\epsilon)}{2\lambda_{N+1}^{1/2-\epsilon}}\right) \lambda_i + q_c + \eta$. Using $\Theta_3 > 0$ and Schur's complement, we infer that $\Gamma_i \le -\lambda_{N+1} + q_c + \eta + \frac{\beta M_{2,\phi}(\epsilon)}{2} \lambda_{N+1}^{1/2+\epsilon} + \frac{2\lambda_{N+1}}{\alpha}$ for all $i \ge N+1$. Using now $\Theta_{2,\eta} \le 0$ and Schur's complement, we obtain that $\Gamma_i \le 0$ for all $i \ge N+1$. Combining this result with $\Theta_{1,\eta} \le 0$, we deduce from (17) that $\dot{V}(t) + 2\eta V(t) \le 0$ for all $t > 0$. From now on, the proof of (27) follows from the same arguments than the ones reported in the previous section. Noting that $z(t, 1) = C x(t)$ in the case (18) while $z_\xi(t, 1) = C x(t)$ in the case (19), we only need to establish the exponential decrease of the term $z_\xi(t, 0)$ to conclude the proof. This is done in Appendix by invoking a small gain argument. □

**Corollary 8** In the context of Theorem 6, the decay rate $\eta > 0$ of the stability estimate (27) is guaranteed provided the LMI conditions $\Theta_{1,\eta} \le 0$, $\Theta_{2,\eta} \le 0$, and $\Theta_3 > 0$ are feasible.

### 4.4 Numerical illustration

We illustrate the results of Theorem 6 and Corollary 8 for the coupled PDE-ODE system described by (19) with $p = 1$, $q = 0$, $q_c = 3$,

$$A = \begin{bmatrix} -1/4 & -1/6 & 2 & 1 & 1/12 \\ -3/2 & -3/2 & 5 & 5 & 1/6 \\ 3/2 & -4 & -15/2 & -5 & -1/3 \\ -13/2 & 22 & 22 & -14 & -1/2 \\ 1/7 & -1/2 & -1/2 & 1/5 & -5/2 \end{bmatrix},$$

$$B = \begin{bmatrix} -5/4 & 2/3 & 1/6 & -1/6 & 0 \end{bmatrix}^\top,$$
<figure><img src="image_5.png" /><figcaption>Fig. 2. Time evolution of the coupled PDE-ODE system (19)</figcaption></figure>

$$C = \begin{bmatrix} -2/5 & -5/4 & 3/2 & 1/3 & 1/40 \end{bmatrix}.$$

In this case, both PDE and ODE systems are open-loop unstable. Indeed, the dominant eigenvalue of the PDE is located approximately at +0.533 while the matrix A has one unstable eigenvalue located approximately at +0.393. The application of Theorem 3 with N = 2 shows the exponential stability of the coupled PDE-ODE system (3). Moreover, the application of Corollary 5 with N = 30 shows the exponential stability of the coupled PDE-ODE system with decay rate η = 0.37. We illustrate this result with a numerical simulation. The numerical scheme consists in the modal approximation of the PDE plant by its 100 dominant modes. The initial condition is set as $w_0(\xi) = 5\xi(1 - \xi)^2 \cos(3\pi\xi)$ and $x_0 = \begin{bmatrix} -1 & 1 & -2 & 2 & -1 \end{bmatrix}^T$. The obtained results are depicted on Fig. 2, confirming the theoretical predictions of Theorem 6 and Corollary 8.

## 5 Conclusion

This paper has addressed the topic of assessing the stability of coupled systems composed of a reaction-

diffusion equation and a finite-dimensional linear time-invariant ordinary differential equation (ODE). The considered coupling channels are located at the boundaries of the PDE and consist of the input and output signals of the ODE. The reported sufficient stability conditions take the form of tractable LMIs and have been derived by adopting a spectral reduction-based method. A distinguished feature is that this approach allows to capture all the possible combinations of either Dirichlet or Neumann boundary conditions in terms of coupling channels. Moreover, we have also assessed the exponential decrease to zero of the aforementioned coupling channels, particularly in the case of Neumann boundary couplings. The drawback of the present Lyapunov function based approach is that the derived stability condition are only sufficient, hence may be not satisfied by some stable reaction-diffusion systems. Nevertheless, as illustrated via the reported numerical examples, this method can be successfully applied to assess the exponential stability of coupled PDE-ODE systems for which both the open-loop PDE and ODE plants are exponentially unstable.

## References

[1] Mohamadreza Ahmadi, Giorgio Valmorbida, and Antonis Papachristodoulou. Dissipation inequalities for the analysis of a class of PDEs. Automatica, 66:163–171, 2016.

[2] Habib Ayadi. Exponential stabilization of cascade ODE-linearized KdV system by boundary dirichlet actuation. European Journal of Control, 43:33–38, 2018.

[3] Matthieu Barreau, Frédéric Gouaisbaut, Alexandre Seuret, and Rifat Sipahi. Input/output stability of a damped string equation coupled with ordinary differential system. International Journal of Robust and Nonlinear Control, 28(18):6053–6069, 2018.

[4] Matthieu Barreau, Carsten Scherer, Frédéric Gouaisbaut, and Alexandre Seuret. Integral quadratic constraints on linear infinite-dimensional systems for robust stability analysis. In IFAC World Congress, 2020.

[5] Matthieu Barreau, Alexandre Seuret, and Frédéric Gouaisbaut. Exponential Lyapunov stability analysis of a drilling mechanism. In 2018 IEEE Conference on Decision and Control (CDC), pages 6579–6584. IEEE, 2018.

[6] Matthieu Barreau, Alexandre Seuret, Frédéric Gouaisbaut, and Lucie Baudouin. Lyapunov stability analysis of a string equation coupled with an ordinary differential system. IEEE Transactions on Automatic Control, 63(11):3850–3857, 2018.

[7] Lucie Baudouin, Alexandre Seuret, and Frédéric Gouaisbaut. Stability analysis of a system coupled to a heat equation. Automatica, 99:195–202, 2019.

[8] Nikolaos Bekiaris-Liberis and Miroslav Krstic. Compensating the distributed effect of a wave pde in the actuation or sensing path of MIMO LTI systems. Systems & Control Letters, 59(11):713–719, 2010.

[9] Abed Boulouz, Hamid Bounit, and Said Hadd. Well-posedness and exponential stability of boundary control systems with dynamic boundary conditions. Systems & Control Letters, 147:104825, 2021.
[10] Jean-Michel Coron and Emmanuel Trélat. Global steady-state controllability of one-dimensional semilinear heat equations. *SIAM Journal on Control and Optimization*, 43(2):549–569, 2004.

[11] Jean-Michel Coron and Emmanuel Trélat. Global steady-state stabilization and controllability of 1D semilinear wave equations. *Commun. Contemp. Math.*, 8(04):535–567, 2006.

[12] Jamal Daafouz, Marius Tucnak, and Julie Valein. Nonlinear control of a coupled PDE/ODE system modeling a switched power converter with a transmission line. *Systems & Control Letters*, 70:92–99, 2014.

[13] Aditya Gahlawat and Matthew M Peet. A convex sum-of-squares approach to analysis, state feedback and output feedback control of parabolic PDEs. *IEEE Transactions on Automatic Control*, 62(4):1636–1651, 2016.

[14] Wei He and Shuzhi Sam Ge. Cooperative control of a nonuniform gantry crane with constrained tension. *Automatica*, 66:146–154, 2016.

[15] Miroslav Krstic. Compensating a string PDE in the actuation or sensing path of an unstable ODE. *IEEE Transactions on Automatic Control*, 54(6):1362–1368, 2009.

[16] Miroslav Krstic. Compensating actuator and sensor dynamics governed by diffusion PDEs. *Systems & Control Letters*, 58(5):372–377, 2009.

[17] Miroslav Krstic and Nikolaos Bekiaris-Liberis. Compensation of infinite-dimensional input dynamics. *Annual Reviews in Control*, 34(2):233–244, 2010.

[18] Hugo Lhachemi and Christophe Prieur. Feedback stabilization of a class of diagonal infinite-dimensional systems with delay boundary control. *IEEE Transactions on Automatic Control*, 66(1):105–120, 2020.

[19] Hugo Lhachemi, David Saussié, and Guchuan Zhu. Boundary feedback stabilization of a flexible wing model under unsteady aerodynamic loads. *Automatica*, 97:73–81, 2018.

[20] Jian Li and Yungang Liu. Adaptive stabilization of coupled PDE-ODE systems with multiple uncertainties. *ESAIM: Control, Optimisation and Calculus of Variations*, 20(2):488–516, 2014.

[21] Vivek Natarajan. Compensating PDE actuator and sensor dynamics using Sylvester equation. *Automatica*, 123:109362, 2021.

[22] Serge Nicaise, Julie Valein, and Emilia Fridman. Stability of the heat and of the wave equations with boundary time-varying delays. *Discrete and Continuous Dynamical Systems* Series S, 2(3):559, 2009.

[23] Yury Orlov. On general properties of eigenvalues and eigenfunctions of a Sturm-Liouville operator: comments on "ISS with respect to boundary disturbances for 1-D parabolic PDEs". *IEEE Transactions on Automatic Control*, 62(11):5970–5973, 2017.

[24] Amnon Pazy. *Semigroups of linear operators and applications to partial differential equations*, volume 44. Springer Science & Business Media, 2012.

[25] Ricardo Sanz, Pedro García, and Miroslav Krstic. Robust compensation of delay and diffusive actuator dynamics without distributed feedback. *IEEE Transactions on Automatic Control*, 64(9):3663–3675, 2018.

[26] Gian Antonio Susto and Miroslav Krstic. Control of PDE-ODE cascades with neumann interconnections. *Journal of the Franklin Institute*, 347(1):284–314, 2010.

[27] Shuxia Tang and Chengkang Xie. State and output feedback boundary control for a coupled PDE-ODE system. *Systems & Control Letters*, 60(8):540–545, 2011.

[28] Ji Wang and Miroslav Krstic. Output feedback boundary control of a heat PDE sandwiched between two ODEs. *IEEE Transactions on Automatic Control*, 64(11):4653–4660, 2019.

[29] Jun-Min Wang, Jun-Jun Liu, Beibei Ren, and Jinhao Chen. Sliding mode control to stabilization of cascaded heat PDE-ODE systems subject to boundary control matched disturbance. *Automatica*, 52:23–34, 2015.

[30] Huai-Ning Wu and Jun-Wei Wang. Static output feedback control via PDE boundary and ODE measurements in linear cascaded ODE-beam systems. *Automatica*, 50(11):2787–2798, 2014.

## A End of the proof of Theorem 6

We investigate the exponential decrease of $z_{\xi}(t, 0)$ to zero. Using the change of variable (20) and the identity $y(t) = Cx(t)$, we have $z_{\xi}(t, 0) = w_{\xi}(t, 0) + Cx(t)$ for $t > 0$. Hence, based on (27), we only need to study the term $w_{\xi}(t, 0) = \sum_{j \ge 1} \phi'_{j}(0)w_{j}(t)$. Let $\alpha_{0} \in (3/4, 1)$ be such that $w_{0} \in D(A^{\alpha_{0}}_{0,n_{R}})$. Let $N_{0} \ge 1$ and $\kappa > 0$ be such that $-\lambda_{n} + q_{c} \le -\eta - \kappa$ for all $n \ge N_{0} + 1$. Consider an arbitrary integer $M \ge N_{0}$. Then we have

$$ \begin{aligned} |w_{\xi}(t, 0)| &\le \sum_{j \ge 1} |\phi'_{j}(0)||w_{j}(t)| \\ &\le C_{M} \sqrt{\sum_{j=1}^{M} w_{j}(t)^2 + C_{\phi} \sqrt{\sum_{j \ge M+1} \lambda_{j}^{2\alpha_{0}} w_{j}(t)^2}} \\ &\le C_{M} \|w(t, \cdot)\|_{L^2} + C_{\phi} \sqrt{\sum_{j \ge M+1} \lambda_{j}^{2\alpha_{0}} w_{j}(t)^2} \end{aligned} $$

where $C_{M} = \sqrt{\sum_{j=1}^{M} \phi'_{j}(0)^2}$ and $C_{\phi} = \sqrt{\sum_{j \ge 1} \frac{\phi'_{j}(0)^2}{\lambda_{j}^{2\alpha_{0}}}} < \infty$. Based again on (20) and (27) we only need to study the term $S_{M}(t) = \sum_{j \ge M+1} \lambda_{j}^{2\alpha_{0}} w_{j}(t)^2$. To do so, we integrate (23a) for $i \ge N_{0} + 1$ and direct estimations give

$$ \begin{aligned} |\lambda_i^{\alpha_0} w_i(t)| &\le e^{(-\lambda_i + q_c)t} |\lambda_i^{\alpha_0} w_i(0)| \\ &+ (|a_i||C|| + |b_i||CA_e||) \int_0^t \lambda_i^{\alpha_0} e^{(-\lambda_i + q_c)(t-s)} \|x(s)\| \mathrm{d}s \\ &+ |b_i||CB|| \int_0^t \lambda_i^{\alpha_0} e^{(-\lambda_i + q_c)(t-s)} |w_{\xi}(s, 0)| \mathrm{d}s \end{aligned} $$

with $A_e = A + BC$. For any function $f \in L_{\text{loc}}^{\infty}(\mathbb{R})$, we note that

$$ \begin{aligned} & \int_0^t \lambda_i^{\alpha_0} e^{(-\lambda_i + q_c)(t-s)} |f(s)| \mathrm{d}s \\ &= e^{-\eta t} \int_0^t \lambda_i^{\alpha_0} e^{(-\lambda_i + q_c + \eta)(t-s)} \times e^{\eta s} |f(s)| \mathrm{d}s \\ &\le \frac{\lambda_i^{\alpha_0}}{\lambda_i - q_c - \eta} e^{-\eta t} \sup_{s \in [0,t]} e^{\eta s} |f(s)|. \end{aligned} $$
Introducing $\mu = \mu(\alpha_0) = (q_c + \eta + \kappa)^{\alpha_0}/\kappa$, we have
$\lambda_i^{\alpha_0}/(\lambda_i - q_c - \eta) \le \mu$ for all $i \ge N_0 + 1$. Combining the
latter estimates, we infer that

$$
\begin{aligned}
& |\lambda_i^{\alpha_0} w_i(t)| \le e^{-\eta t} |\lambda_i^{\alpha_0} w_i(0)| \\
& + \mu (|a_i| \|C\| + |b_i| \|CA_e\|) e^{-\eta t} \sup_{s \in [0,t]} e^{\eta s} \|x(s)\| \\
& + \mu C_M |b_i| \|CB\| e^{-\eta t} \sup_{s \in [0,t]} e^{\eta s} \|w(s, \cdot)\|_{L^2}^2 \\
& + \mu C_\phi |b_i| \|CB\| e^{-\eta t} \sup_{s \in [0,t]} e^{\eta s} \sqrt{S_M(s)}
\end{aligned}
$$

for all $i \ge N_0 + 1$ and all $t \ge 0$. The use of Young's
inequality and summing for $i \ge M + 1 \ge N_0 + 1$ we
obtain that

$$
\begin{aligned}
S_M(t) &\le 5e^{-2\eta t} S_M(0) \\
&\quad + 5\mu^2 (\|\mathcal{P}_M a\|_{L^2}^2 \|C\|^2 + \|\mathcal{P}_M b\|_{L^2}^2 \|CA_e\|^2) \\
&\quad \times e^{-2\eta t} \sup_{s \in [0,t]} e^{2\eta s} \|x(s)\|^2 \\
&\quad + 5\mu^2 C_M^2 \|\mathcal{P}_M b\|_{L^2}^2 \|CB\|^2 e^{-2\eta t} \sup_{s \in [0,t]} e^{2\eta s} \|w(s, \cdot)\|_{L^2}^2 \\
&\quad + 5\mu^2 C_\phi^2 \|\mathcal{P}_M b\|_{L^2}^2 \|CB\|^2 e^{-2\eta t} \sup_{s \in [0,t]} e^{2\eta s} S_M(s)
\end{aligned}
$$

for all $t \ge 0$, hence

$$
\begin{aligned}
\sup_{s \in [0,t]} e^{2\eta s} S_M(s) &\le 5S_M(0) \\
&\quad + 5\mu^2 (\|\mathcal{P}_M a\|_{L^2}^2 \|C\|^2 + \|\mathcal{P}_M b\|_{L^2}^2 \|CA_e\|^2) \\
&\quad \times \sup_{s \in [0,t]} e^{2\eta s} \|x(s)\|^2 \\
&\quad + 5\mu^2 C_M^2 \|\mathcal{P}_M b\|_{L^2}^2 \|CB\|^2 \sup_{s \in [0,t]} e^{2\eta s} \|w(s, \cdot)\|_{L^2}^2 \\
&\quad + 5\mu^2 C_\phi^2 \|\mathcal{P}_M b\|_{L^2}^2 \|CB\|^2 \sup_{s \in [0,t]} e^{2\eta s} S_M(s).
\end{aligned}
$$

Since $\|\mathcal{P}_M b\|_{L^2} \to 0$ when $M \to +\infty$, we infer the existence of a large enough integer $M \ge N_0 + 1$, independent of the initial condition $w_0 \in D(\mathcal{A}_{0,n_R}^{\alpha_0})$, such that $5\mu^2 C_\phi^2 \|\mathcal{P}_M b\|_{L^2}^2 \|CB\|^2 < 1$. Since all the supremums appearing in the latter estimate are finite (recall that $\mathcal{A}_{0,n_R}^{\alpha_0} w \in C^0(\mathbb{R}_{\ge 0}; L^2(0,1)))$, we obtain the existence of a constant $M_4 > 0$ such that

$$
\sup_{s \in [0,t]} e^{2\eta s} S_M(s) \le M_4 S_M(0) + M_4 \sup_{s \in [0,t]} e^{2\eta s} \|x(s)\|^2 \\
+ M_4 \sup_{s \in [0,t]} e^{2\eta s} \|w(s, \cdot)\|_{L^2}^2
$$

for all $t \ge 0$. Noting that $S_M(0) \le \|\mathcal{A}_{0,n_R}^{\alpha_0} w_0\|_{L^2}^2$, the
claimed conclusion follows from (20) and (27). In the
case, $w_0 \in D(\mathcal{A}_{0,n_R}) \subset D(\mathcal{A}_{0,n_R}^{\alpha_0})$ it can be easily seen
that $\|\mathcal{A}_{0,n_R}^{\alpha_0} w_0\|_{L^2}^2 \le \|w_0\|_{L^2}^2 + \|\mathcal{A}_{0,n_R} w_0\|_{L^2}^2$, which con-
cludes the proof.