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

HAL
open science

# Input / Output Stability of a Damped String Equation coupled with Ordinary Differential System

Matthieu Barreau, Frédéric Gouaisbaut, Alexandre Seuret, Rifat Sipahi

► To cite this version:

Matthieu Barreau, Frédéric Gouaisbaut, Alexandre Seuret, Rifat Sipahi. Input / Output Stability of a Damped String Equation coupled with Ordinary Differential System. International Journal of Robust and Nonlinear Control, 2018, 28 (18), pp.6053-6069. 10.1002/rnc.4357 . hal-01690626v2

HAL Id: hal-01690626

https://laas.hal.science/hal-01690626v2

Submitted on 19 Sep 2018

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.

Public Domain
# Input / Output Stability of a Damped String Equation coupled with Ordinary Differential System*

September 19, 2018

Matthieu Barreau¹*, Frédéric Gouaisbaut¹, Alexandre Seuret¹, and Rifat Sipahi².

¹ LAAS - CNRS, Université de Toulouse, CNRS, UPS, France (barreau,fgouaisb,seuret@laas.fr)

² Department of Mechanical and Industrial Engineering, 334 Snell Engineering Center, 360 Huntington Avenue, Northeastern University, Boston, MA 02115 USA (rifat@coe.neu.edu)

## Abstract

The input/output stability of an interconnected system composed of an ordinary differential equation and a damped string equation is studied. Issued from the literature on time-delay systems, an exact stability result is firstly derived using pole locations. Then, based on the Small-Gain theorem and on the Quadratic Separation framework, some robust stability criteria are provided. The latter follows from a projection of the infinite dimensional system states onto an orthogonal basis of Legendre polynomials. Numerical examples comparing these results with the ones in the literature are presented along with demonstrations of the effectiveness of the developed robust stability criteria.

**Keywords:** String equation; Frequency approach; Robust stability; Coupled ODE/PDE; Quadratic Separation.

## 1 Introduction

In a number of complex systems, including cyberphysical systems, one particular interconnection that often arises includes an interaction between a finite and an infinite dimensional system, see for example, a flexible crane system [21], a drilling pipe [7], temperature control [46], vibration in structures [24, 22], MEMS[15] among many others[4]. In the presence of such interactions however, studying stability is a challenge since one cannot separately treat each system, but must consider them both in an ensemble[1, 38]. This paper is aimed at addressing this stability problem, specifically of an interconnection between an ordinary differential equation (ODE) and a partial differential equation (PDE). We consider here only the string equation, which falls into the one-dimension hyperbolic systems. Most of the work on the stability of this class of PDEs has been conducted considering a Lyapunov functional. The main techniques can be classified in two categories. Of the first kind is the backstepping methodology for boundary controlled PDEs. This approach has been utilized by many researchers, see a summary in[29, 28]. The other method relies on semi-group theory and uses an energy argument to conclude on stability. This method has also been widely studied, see, for example[21, 47, 34, 37].

Noticing that the PDEs related to transport or string can be treated as a time-delay system (TDS), another approach considered in the study of stability of such systems relies on their input/output stability[17]. The derivation of a transfer function and the study of the characteristic equation forms the core of this approach toward assessing its input/output stability. On the wave equation, further results can be found in[18, 34, 16, 49] considering an input/output map or also using the port-Hamiltonian approach[48, 9]. On the other hand, in the cited references, the focus is on the PDE system, while the study of a coupled ODE/PDE system requires additional developments.

<footnote>*This work is supported by the ANR project SCIDIS contract number 15-CE23-0014.</footnote>
To address the stability problem of the coupled ODE/PDE system at hand, here we first make use of the comparison between TDS and string equation to obtain a transfer function representing the interconnected dynamics. An exact stability criterion similar to the one developed in [39] is next adopted for this transfer function based on the pole locations of the system characteristic equation on the imaginary axis of the complex plane. Even if this approach provides exact results for stability, robustness and synthesis problems must be separately dealt with. To this end, a complementary approach is to use robust stability analysis tools as for example in [13] which considers an event-trigger control. In the same direction, the Small-Gain theorem [11] and Quadratic Separation [41, 2] provide, as we show, efficient frameworks to reveal the stability conditions for the interconnected ODE/PDE system in a parameter region of interest.

The outline of the paper is as follows. Section 2 states the problem and provides some explanations regarding its physical interpretation. Section 3 deals with the input/output stability of the coupled system whose transfer function is derived and properties are deduced. Based on this, an exact stability criterion is developed in Section 4. Section 5 is dedicated to two tools used in robust stability analysis. More particularly, the Small-Gain theorem provides a very conservative but simple stability test while Quadratic Separation is developed in two stages: a first and then an extended analysis using projections of the infinite dimensional states of the ODE/PDE system. Finally, Section 6 presents examples and comparisons with other techniques.

The main contributions of this paper are as follows. The formulation of the ODE/PDE interconnection in the form of a TDS enables one to obtain some fundamental properties of the ODE/PDE system stability, and furthermore an algorithm for an exact stability criterion can be constructed as a baseline. Moreover, here, we propose a simple robust stability test based on Small-Gain Theorem from which many properties of the coupled system are deduced. Finally, the reformulation of the interconnected system in the Quadratic Separation framework enables a refinement of the previous stability results thanks to the Bessel inequality and Legendre polynomials. This new formulation shows that a new robust stability criterion can be obtained, which needs less computational resources than the one developed in [3].

**Notations:** The following notations are used: $\mathbb{R}^+ = [0, +\infty)$, $\mathbb{C}^+ = \{z \in \mathbb{C} / \Re(z) \ge 0, z \ne 0\}$, $L^2 = L^2([0, 1]; \mathbb{R})$, $H^n = \{z \in L^2; \forall m \le n, \frac{\partial^m z}{\partial x^m} \in L^2\}$ and $i$ is the imaginary number. $L^2$ is equipped with the norm $\|z\|^2 = \int_0^1 |z(x)|^2 dx = \langle z, z \rangle$.

For any square matrices $A$ and $B$, let $\text{diag}(A, B) = [\begin{matrix} A & 0 \\ 0 & B \end{matrix}]$. A symmetric positive definite matrix $P \in \mathbb{S}_+^n \subset \mathbb{R}^{n \times n}$ is written as $P \succ 0$. For a matrix $A \in \mathbb{R}^{m \times p}$, then $A(1:N, :)$ (or, $A(:, 1:N)$) is the sub-matrix of $A$ with the first $N$ lines (or, columns respectively). $A^\perp$ is the nullspace of $A$ and $*$ denotes the transconjugate. $0_{m,p}$ and $1_{m,p}$ are respectively the null and full of ones matrices with $m$ lines and $p$ rows and $0_m = 0_{m,m}$, $1_m = 1_{m,m}$. To keep the presentation concise, the following notations are used for partial derivatives of a function $u$ on an appropriate space: $u_t(x, t) = \frac{\partial u}{\partial t}(x, t)$, $u_{tt}(x, t) = \frac{\partial^2 u}{\partial t^2}(x, t)$, $u_x(x, t) = \frac{\partial u}{\partial x}(x, t)$ and $u_{xx}(x, t) = \frac{\partial^2 u}{\partial x^2}(x, t)$.

# 2 Problem Statement

## 2.1 Problem formulation

The input/output stability of the following interconnected system is studied:

$$ \dot{X}(t) = AX(t) + B(u(1,t) + r(t)), \qquad t \ge 0, \qquad (1a) $$

$$ u_{tt}(x,t) = c^2 u_{xx}(x,t), \qquad x \in [0,1], t \ge 0, \qquad (1b) $$

$$ u(0,t) = KX(t), \qquad t \ge 0, \qquad (1c) $$

$$ u_x(1,t) = -c_0 u_t(1,t), \qquad t \ge 0, \qquad (1d) $$

$$ u(x, 0) = u_0(x), \qquad x \in [0, 1], \qquad (1e) $$

$$ u_t(x, 0) = v_0(x), \qquad x \in [0, 1], \qquad (1f) $$

$$ X(0) = X_0, \qquad (1g) $$

with the initial conditions $X_0 \in \mathbb{R}^n$ and $(u_0, v_0) \in H^2 \times H^1$ such that equations (1c) and (1d) are respected. Here, $(X_0, u_0, v_0)$ are compatibles with the boundary conditions and $r$ is the input function.

System (1) represents the interconnection of a linear time invariant system with a string equation. Here, the state $u$ denotes the amplitude of the wave which belongs to a functional space and consequently
is of infinite dimension. We assume that the state $u(x,t)$ belongs to $\mathbb{R}$ but the analysis can be extended without difficulty to $\mathbb{R}^m$, $m > 1$. The two subsystems are connected through a Dirichlet boundary condition, that is, the output $u(1, \cdot)$ and the input (eq. (1c)) depend on the state $u$ at different positions but not on the derivative of $u$.

For the system to be well-posed, another boundary condition is also needed. In the literature, the choice is to use a boundary damping with equation (1d). Indeed, it has been shown in [30] for example that this condition ensures the stability of the wave equation itself if $c_0 > 0$ and the case $c_0 = 0$ removes the damping.

The wave equation operator associated with the above described boundary conditions is known to be diagonalizable (for more information on semi-group definitions and properties, the reader can refer to [33, 47] and the references therein). Once diagonalized, it can be expressed as the composition of two transport equations, one going forward and another backward. Boundary condition (1d) implies a reflection of the forward wave with a coefficient $\alpha = \frac{1-cc_0}{1+cc_0}$ as explained in [32]. Enforcing $c_0 > 0$ implies $|\alpha| < 1$ and consequently the energy of the wave is decreasing. These details will be the key in developing our results later in the paper.

It is critical to note that the string equation considered here behaves like a communication channel with a string dynamics. Note that most of the systems with a string equation however consider a col-located control [35, 23] which then yields a very different problem in nature compared to the one in the proposed interconnection. The purpose of this paper is to study the input/output stability of this inter-connected system, where the definition of input/output stability for infinite dimensional systems follows from Definition 9.1.1 of [11] as:

**Definition 1** System (1) is said to be **input/output stable** if for all energy bounded input $r$, the energy of the output $Y = KX$ is also bounded.

## 2.2 Existence and Uniqueness of the Solution

Before going further, the existence of a solution to system (1) with $r = 0$ is proven. To do so, a step by step procedure is proposed. To ease the reading, the following sets are defined:

$$ 
\begin{aligned}
\mathbb{I}_n &= [nc^{-1}, (n+1)c^{-1}), \quad n \in \mathbb{N}, \\
\mathbb{H} &= \{(X, u, v) \in \mathbb{R}^n \times \mathbb{H}^2 \times \mathbb{H}^1\}, \\
\mathbb{D} &= \{(X, u, v) \in \mathbb{H} \text{ s.t. } u(0) = KX, u_x(1) = -c_0v(1)\},
\end{aligned} 
\quad (2) 
$$

and using the same technique as in [32, 14] for $(X_0, u_0, v_0) \in \mathbb{D}$, a solution for $t \in \mathbb{I}_0$ is given by:

$$ u(1,t) = u_0(1-ct) + \frac{\alpha}{2}(u_0(1-tc) + u_0(1)) + \frac{1+\alpha}{2} \int_{1-ct}^{1} v_0(s)ds. \quad (3) $$

On $\mathbb{I}_0$, we notice that $u(1, \cdot) \in H^2(\mathbb{I}_0)$ depends only on the initial conditions. Then, $X$ on $I_0$ is the solution of the differential equation (1a), with initial condition $X(0) = X_0$. It is well known that this solution reads:

$$ \forall t \in \mathbb{I}_0, \quad X(t) = e^{At}X_0 + \int_0^t e^{A(t-s)}Bu(1,s)ds. \quad (4) $$

Then, $X \in C^1(\mathbb{I}_0)$ and consequently, $u(0, \cdot)$ has the same regularity property. On the domain $\Gamma_0 = \mathbb{I}_0 \times [0, 1]$, with the boundary conditions defined previously and using the formulas of [32], Equation (1b) has a unique solution $(x, t) \in \Gamma_0 \mapsto u(x, t) \in H^2$ and for $t \in \mathbb{I}_0$, $u_t(\cdot, t) \in H^1$. With all these considerations, it is possible to find $(X(c^{-1}), u(\cdot, c^{-1}), u_t(\cdot, c^{-1})) \in \mathbb{D}$ and then repeat the same procedure for $\mathbb{I}_1$ leading to the existence of a solution. We then get $(X, u, u_t) \in L^2(0, +\infty, \mathbb{H})$. Moreover, $u, u_t, u_x, u_{tt}$ and $u_{xx}$ are $L^2$ on each compact set of $\mathbb{R}^+$. As the system is linear, energy consideration leads to the uniqueness property. If $r$ is not null, then we can calculate the steady-state $(X_f, u_f)$ and then the previous analysis holds for $(X - X_f, u - u_f)$.

**Remark 1** The set of solutions $\mathbb{H}$ is smaller than the one in [3]. The main difference comes from considering the strong solution here while the weak solution will belong to the same space as in [3].
<figure><img src="image_2.png" /><figcaption>Figure 1: Block diagram for the system described by equation (1).</figcaption></figure>

# 3 Input / Output Analysis

One possible way to study the input/output stability of system (1) is to study its characteristic equation, which can be obtained by applying the Laplace transform. This standard approach can be found in [12, 36].

## 3.1 Laplace transform of the wave equation

The notion of Laplace transform extends easily from its traditional definition for finite dimensional systems as noted in [12]. For this, the variables need to be $L^2$ on each compact set of $\mathbb{R}^+$, which is the case in (1) as discussed earlier. Since we are not concerned with the transient response of the dynamics, we assume the initial conditions to be zero, i.e., $u^0 = v^0 = 0$. In Laplace domain, equation (1b) can then be transformed into:

$$ \forall x \in [0, 1], \qquad s^2 U_{xx}(x, s) = c^2 U(x, s), \qquad (5) $$

with $s \in \mathbb{C}$ is the Laplace variable and $U(x, \cdot)$ is the Laplace transform of $u(x, \cdot)$. The solution defined on $x \in [0, 1]$ is given by:

$$ U(x, s) = C_1(s)e^{\frac{s}{c}x} + C_2(s)e^{-\frac{s}{c}x}, \qquad (6) $$

where $C_1$ and $C_2$ are space-independent transfer functions to be determined using the boundary conditions. Once calculated using equations (1c) and (1d), the transfer function for the string equation is obtained as:

$$ W(x, s) = \frac{U(x, s)}{U(0, s)} = \frac{e^{-\frac{s}{c}x} + \alpha e^{\frac{s}{c}(x-2)}}{1 + \alpha e^{-2\frac{s}{c}}}, \qquad (7) $$

with $\alpha = \frac{1-cc_0}{1+cc_0}$ and $x \in [0, 1]$. The transfer function from $U(0, s)$ to $U(1, s)$ then reads:

$$ W(1, s) = W(s) = \frac{U(1, s)}{U(0, s)} = \frac{(1 + \alpha)e^{-s/c}}{1 + \alpha e^{-2s/c}} = \frac{2e^{-s/c}}{1 + cc_0 + (1 - cc_0)e^{-2s/c}}. \qquad (8) $$

Notice that there are infinitely many poles of $W(x, s)$. These poles are independent of $x \in [0, 1]$ and their real part $\frac{c}{2} \log|\alpha|$ is strictly negative if $cc_0 \neq 1$ and $c_0 > 0$.

**Remark 2** When $c_0 = 0$, $W$ has infinitively many poles on the imaginary axis. Then Corollary 9.1.4 from [11] applies and there does not exist any finite-dimensional controller which exponentially stabilizes the wave equation. That is why only the case $c_0 > 0$ is considered.

## 3.2 Transfer function for the coupled system

Considering the block diagram in Figure 1, the transfer function of the finite dimensional system is:

$$ \mathcal{H}(s) = \frac{Y(s)}{U(1, s) + R(s)} = K(sI - A)^{-1}B = \frac{N(s)}{D(s)} = \frac{l_{n-1}s^{n-1} + l_{n-2}s^{n-2} + \cdots + l_0}{m_n s^n + m_{n-1}s^{n-1} + \cdots + m_0}, \qquad (9) $$

where $N$ and $D$ are real polynomials in $s$ with $m_n \neq 0$ and $R$ is the Laplace transform of the reference signal $r$ in (1).
Using transfer function $\mathcal{W}$ developed previously, we then obtain the following closed-loop system transfer function:

$$ \mathcal{F}(s) = \frac{Y(s)}{R(s)} = \frac{\mathcal{H}(s)}{1 - \mathcal{H}(s)\mathcal{W}(s)} = \frac{N(s)(1 + \alpha e^{-2s/c})}{(1 + cc_0)c_{eq}(s, c)}, \quad (10) $$

with

$$ \begin{aligned} c_{eq}(s, c) &= (1 + cc_0 + (1 - cc_0)e^{-2s/c}) D(s) - 2N(s)e^{-s/c} \\ &= m_n (1 + cc_0 + (1 - cc_0)e^{-2s/c}) e^{-2s/c} s^n + \sum_{k=0}^{n-1} ((1 + cc_0 + (1 - cc_0)e^{-2s/c}) m_k - 2l_k) s^k. \end{aligned} \quad (11) $$

The stability of the dynamics is determined by the location of the zeros of (11) on the complex plane as we discuss in detail in the following subsection.

### 3.3 Input / Output stability

A time-delay system is of *neutral type* when its highest order derivative is affected by a delay term $\tau$. Given a characteristic equation, this feature appears as the highest power of the Laplace variable $s$ multiplying the time delay operator $e^{-\tau s}$. Similarly, in state-space representation, neutral property is easily detected by the presence of the derivative of the state term being affected by delays. For a complete treatment of such systems, see details on neutral functional differential equations (NFDE) in [38, 27].

Inspecting the corresponding transfer function $\mathcal{F}$ or the characteristic equation (11) of system (1), it is easy to see that this system exhibits the neutral type property. That is, from an input/output approach, system (1) is classified as a neutral system. Moreover, we can re-write system (1) in the standard NFDE form by collecting the derivative of the states on the left hand side:

$$ \dot{X}(t) + \alpha \dot{X}(t - 2c^{-1}) = AX(t) + \frac{2}{1 + cc_0} BKX(t - c^{-1}) + \alpha AX(t - 2c^{-1}) + Br(t) + \alpha Br(t - c^{-1}), \quad (12) $$

where from [20, 40] the difference operator is defined as

$$ D(X)(t) = X(t) + \alpha X(t - 2c^{-1}), $$

for all $t \in \mathbb{R}^+$.

**Remark 3** One can see that for $cc_0 = 1$, we have $\alpha = 0$ then system (1) is no longer a neutral type system. It is of retarded type [5]. This distinction is critical when studying the stability of TDS. It will also help build various example cases and enable a comparison with the stability tests already available for TDS.

When it comes to studying the stability of NFDEs, contrary to functional differential equation of retarded type, inspecting the pole locations of the corresponding characteristic equation is not sufficient to guarantee input/output stability. This is due to the well-known small-delay phenomenon, also known as small $\tau$-stabilizability. When an NFDE is small $\tau$-stabilizable, then this guarantees that the poles of the NFDE behave in a continuum as the delay value changes from 0 to $0^+$. In other words, the stability properties of the NFDE are preserved when the delay is infinitesimally increased past zero. Whenever this property holds, then studying the pole locations of an NFDE allows one to conclude on the stability of an NFDE (see [5, 19] for more information).

**Definition 2** An NFDE system is said to be **small $\tau$-stabilizable** if the difference operator $D(X)(t)$ is stable.

This definition comes from [20, 40]. According to Theorem 12.5.1 and Corollary 12.5.1 from [20], a necessary and sufficient condition for the difference operator to be stable is that $|\alpha|$ is strictly less than 1.

**Proposition 1** System (1) is neutral and small $\tau$-stabilizable if and only if $c_0 > 0$ and $cc_0 \neq 1$.

Finally, since the system is small $\tau$-stabilizable under a given condition, the input/output stability can be discussed.
**Corollary 1** If $(X_0, u_0, v_0) \in \mathbb{D}$ and all the poles of $\mathcal{F}$ as defined in equation (10) are on the left half of the complex plane, then system (1) is input/output stable. The equilibrium points $(X_e, u_e, v_e) \in \mathbb{H}$ are such that $(A + BK)X_e = 0, u_e = KX_e$ and $v_e = 0$.

Moreover, if the system $\Sigma(A, B, K)$ is stabilizable and detectable such that $A + BK$ is not singular, then system (1) is asymptotically stable.

**Proof:** If system (1) is small $\tau$-stabilizable, then it is input/output stable if all the poles of $\mathcal{F}$ are with a strictly negative real part (Theorem 9.9.1 from [19]). The link between asymptotic and input/output stability comes from the pole/zero cancellation. If the system $\Sigma(A, B, K)$ is stabilizable and detectable, then, the pole/zero simplification in $\mathcal{F}$ are with negative real parts. It means that the non-observable or non-controllable parts of system (1) are stable. Then, all states will be indeed converging. The study of equilibrium points has been established in Proposition 2 of [3] and if $A + BK$ is not singular, state X converges toward zero and the system is asymptotically stable. $\square$

Now that we established some key properties regarding the characteristics of (1), the stability analysis can be pursued. We consider two main approaches, one using pole locations and another with robust stability tools.

# 4 A frequency domain argument on stability

Once the characteristic equation is established, different methodologies can be applied to assert the input/output stability of system (1). Various techniques can be adopted for this purpose, see for example [44], including those where $c_{eq}(\cdot, c)$ in (11) has the coefficients explicitly depending on the delay $c^{-1}$. This problem has been investigated in some studies [6, 26]. In this part, we focus on the exact treatment of this problem coming from the time-delay system literature. Specifically, a pole location argument is taken based on frequency domain tools in order to assess stability. We study the stability problem using the Cluster Treatment of Characteristic Roots (CTCR) methodology originally proposed in [39]. Let $\tau = c^{-1}$ and regrouping the terms by their delay dependence, characteristic equation (11) becomes:

$$c_{eq}(s, \tau) = a_0(s, \tau) + a_1(s, \tau)e^{-\tau s} + a_2(s, \tau)e^{-2\tau s}, \quad (13)$$

with $a_0(s, \tau) = (1 + \frac{c_0}{\tau}) D(s)$, $a_1(s, \tau) = -2N(s)$ and $a_2(s, \tau) = (1 - \frac{c_0}{\tau}) D(s)$. We then claim that the system described by transfer function $\mathcal{F}$ can switch from stable to unstable behavior, or vice-versa, for a given delay only if (11) has a pole on the imaginary axis [39].

CTCR starts by exhaustively detecting all the imaginary axis poles $s = i\omega$, along with their corresponding delay values (Proposition 1 of [39]). Next, CTCR identifies that each pole $s = i\omega$ has a unique crossing direction over the imaginary axis for all the delays creating this crossing (Proposition 2 of [39]). Knowing the number of poles at $\tau = 0$, which is trivial to assess, it is then possible to use the information regarding crossing directions and at which delay values such crossings occur to track the pole locations across the imaginary axis. With this idea, it becomes possible to count the number of poles on the right-half plane for a given delay value $\tau > 0$. Whenever there are no unstable poles for certain delays, we then state that the system at hand is input/output stable for those delays.

CTCR framework has already been demonstrated on neutral systems [40] while respecting the small $\tau$-stabilizability property. In the following, we briefly summarize this framework from the cited studies and also point out the main differences introduced while studying the particular system (1). Following the process described in [39, 40], the Rekasius substitution [42] is defined as:

$$e^{-\tau s} := \frac{1 - Ts}{1 + Ts}, \quad \tau \in \mathbb{R}^{+}, \quad T \in \mathbb{R}, \quad s = i\omega, \quad (14)$$

which is an exact substitution of exponential terms when $s = i\omega$. This substitution is different than Padé approximation since in general $T \neq \tau/2$. Next, substituting (14) into the characteristic equation (11) and expanding by $(1 + Ts)^2$, which does not bring any artificial imaginary poles, we obtain a transformed characteristic equation:

$$\bar{c}_{eq}(s, T) = \left( 1 + \frac{2c_0}{\tau}Ts + T^2s^2 \right) D(s) - N(s)(1 - T^2s^2),$$
which is nothing but a multinomial. Moreover the imaginary poles $s = i\omega$ of the original characteristic equation and this multinomial are identical [39] and hence one can alternatively compute the imaginary poles $s = i\omega$ from this multinomial, which is a much easier task. To this end, first build the coefficients $\{b_k(T, \tau)\}_{k \in (0, n+2)}$, such that:

$$ \bar{c}_{eq}(s, T) = \sum_{k=0}^{n+2} b_k(T, \tau) s^k. \qquad (15) $$

Next notice that the presence of delay $\tau$ as a coefficient in (15) is not a standard form, which may complicate the computation of imaginary poles and the application of Proposition 2 in [39] as coefficients influenced by delays may also rule out certain periodicity properties. This issue is removable in the particular problem at hand. We consider the following manipulation: as $c_0$ is always divided by the delay $\tau$ in $b_k$, if one defines a new positive variable $c_1 = c_0\tau^{-1}$, then $b_k$ depends only on $c_1$ and not on $\tau$ anymore. Working at a given strictly positive $c_1$ removes the dependence of $b_k$ in $\tau$. Since they are now independent, the methodology still applies. This manipulation also shows that $c_1$ is a variable of interest when it comes to studying the parametric stability boundaries of damped waves.

Using $b_k$ and $a_k$ as defined in equations (13) and (15), the CTCR methodology provides the boundaries of the parameter space in which the input-output stability of system (1) holds<fnref n="1" />. This method is summarized as follows for a given $c_1 = c_0\tau^{-1}$:

1. Using $\bar{c}_{eq}$, find the roots $s = i\omega$ corresponding to $T \in \mathbb{R}$, e.g., by using Routh's array. There are only a finite number of such solutions (Proposition 1 of [39]);

2. For each $T \in \mathbb{R}$ obtained previously, we already have $\omega$ where a crossing on the imaginary axis exists. For each $\omega$, there is a root tendency, indicating the unique direction of the crossing independent of the delays creating that crossing (Proposition 2 of [39]).

3. Then, using the inverse transformation of Rekasius transformation in (14), it is possible to find all the delays ($\tau_\ell$) corresponding to each pair of $(T, \omega)$ and falling in an interval from 0 up to a target delay value $\tau_{max}$. Sorting these delays in ascending order, and starting with the number of unstable roots for $\tau = 0$, the number of unstable roots for a delay $\tau < \tau_{max}$ can be accounted by observing the root tendency property of the crossings.

4. The stability areas of system (1) for a given $c_1 = c_0\tau^{-1}$ are when no unstable poles are detected. Since $\tau$ is known at this point, $c_0$ can be recovered: $c_0 = c_1\tau$.

We point out that the coefficients $b_k$ depend linearly on $c_1$ and then the roots of $\bar{c}_{eq}(\cdot, T)$ vary continuously subject to $c_1$. In other words, the delays $\tau$ for which there is a crossing vary continuously relative to $c_1$. The border of each stability area on the map $(c_1, \tau)$ is consequently continuous. This is one of the arguments for considering the mapping $(c_1, \tau)$ instead of $(c_0, c)$.

# 5 Robust Stability Analysis

The previous stability analysis results are exact, however, we still need to develop efficient tools to deal with robustness issues with respect to, for example, an uncertainty on $A$ or $c$. Here, we aim at obtaining tools not only for the specific problem (1), but for a general interconnection dynamic. That is why tools coming from the robust analysis are also considered. In [3], the authors built a Lyapunov functional to ensure the exponential stability and the main drawback was an important number of decision variables making the treatment computationally demanding for large scale systems. This section aims at providing analysis tools with similar results but of lower computational complexity.

Since the studied system is an interconnection between two subsystems, here, the stability of the interconnection is stated under some conditions on each subsystem. Two tools coming from the robust analysis are considered: Small Gain theorem and Quadratic Separation.

<footnote n="1">¹Here, we provide the general framework of CTCR for an exact stability analysis. Handling degenerate cases requires care as was demonstrated in [45]. Such special cases, for a slightly different $c_1$, will however not arise.</footnote>
## 5.1 Small-Gain Theorem

We consider the block diagram of Figure 1 where the wave equation is treated as a disturbance. The transfer functions of the disturbance and the plant are borrowed from Section 3. The following stability criteria is a direct application of the Small-Gain Theorem.

**Theorem 1** Let system (1) with A Hurwitz with its $H_\infty$ less than 1. The system is input-output stable for $(c, c_0)$ if the following condition is satisfied:

$$ \|\mathcal{H}\|_{\infty} = \max_{w \in \mathbb{R}^{+}} |\mathcal{H}(iw)| < cc_0. \qquad (16) $$

**Remark 4** Note that the condition $c_0 > 0$ (and consequently $|\alpha| < 1$) is equivalent to $\|\mathcal{W}\|_{\infty}$ bounded.

**Proof:** Beforehand, the infinity norm of the disturbance is computed. Some calculations lead to $\|\mathcal{W}\|_{\infty} = \frac{2}{(1+cc_0) \min_{\omega \ge 0} |1+\alpha e^{-2i\omega/c}|}$. For $|\alpha| < 1$, the function $\omega \mapsto 1 + \alpha e^{-2i\tau\omega}$ is inside the circle on the complex plane centered at 1 and of radius $|\alpha| < 1$. The minimum of $\|\mathcal{W}\|_{\infty}$ is $\max((cc_0)^{-1}, 1)$. Then, using the small gain theorem for infinite dimensional systems as stated in Theorem 9.1.7 by [11], the interconnected system is asymptotically stable if $\|\mathcal{H}\|_{\infty}\|\mathcal{W}\|_{\infty} < 1$ and each subsystem is stable. The second condition is ensured if $\|\mathcal{H}\|_{\infty} < 1$ and $c_0 > 0$. Considering the case $cc_0 \le 1$ leads to $\|\mathcal{W}\|_{\infty} = (cc_0)^{-1}$, it then follows that robust stability is ensured. If $cc_0 > 1$, then the stability is ensured if $\|\mathcal{H}\|_{\infty} < 1$, which is true by assumption. $\square$

This theorem is quite conservative, mostly because the Small-Gain Theorem provides only a sufficient condition. However the resulting stability test is simple and some properties can be deduced.

**Proposition 2** If $\|\mathcal{H}\|_{\infty} < 1$, then there exists a function $c_0 \mapsto c_{min}(c_0)$ where $c_{min}(c_0) \le \|\mathcal{H}\|_{\infty}c_0^{-1}$. That leads to three properties:

1. Since $\|\mathcal{H}\|_{\infty} < 1$, A + BK is stable,

2. For a given $c_0 > 0$, system (1) is stable for all $c \ge c_{min}(c_0)$,

3. $\lim_{c_0 \to +\infty} c_{min}(c_0) = 0$.

To reduce the conservatism introduced in this subpart, another framework is proposed.

## 5.2 Quadratic Separation - Preliminary result

The Small-Gain Theorem ensures the exponential stability of an interconnected system composed of a disturbance and a plant both stable. To decrease the conservatism and consider a broader class of interconnected systems, the Quadratic Separation (QS) framework can be used. This framework has been originally proposed in [25] and it studies the well-posedness of a closed-loop system made up of an unknown disturbance and a plant.

We describe the methodology of QS to provide a preliminary result on system (1) and then we extend this stability analysis to a more general case. QS states the well-posedness of a generic system described in Figure 2, where $\nabla$ is called the uncertainty matrix and belongs to a set $\mathcal{W}$. The well-posedness is defined as follows:

**Definition 3** The interconnected system described on Figure 2 is well-posed with respect to the norm $\|\cdot\|$ if

$$ \exists \gamma > 0, \forall \nabla \in \mathcal{W}, \forall \bar{\omega}, \bar{z}, \quad \left\| \begin{bmatrix} \omega \\ z \end{bmatrix} \right\| \le \gamma \left\| \begin{bmatrix} \bar{\omega} \\ \bar{z} \end{bmatrix} \right\|, \qquad (17) $$

where $\bar{\omega}$ and $\bar{z}$ are the references.

In our case here, we consider the disturbance to be an operator acting on two signals w and z in Laplace domain. In other words, the following system is obtained:

$$ \left\{ \begin{array}{l} \Omega(s) = \nabla(s)Z(s), \\ \mathcal{E}Z(s) = \mathcal{A}\Omega(s), \end{array} \right. \qquad (18) $$
<figure><img src="image_3.png" /><figcaption>Figure 2: Feedback system for quadratic separation. $\mathcal{A}$ and $\mathcal{E}$ are matrices and $\nabla \in \mathbb{W}$ is the uncertainty. $\bar{\omega}$ and $\bar{z}$ are the references.</figcaption></figure>

with $\mathcal{E}$ full column rank, $Z$ and $\Omega$ are the Laplace transformations of $z$ and $\omega$, respectively and $\mathcal{A}$ is a description of the system, its explicit form is given after. The well-posedness of system (18) is assessed in Theorem 1 and Corollary 2 from [41]. Following this formulation, the next theorem is stated.

**Theorem 2** The system described by Figure 2 and Equation (18) is well-posed if and only if there exists a real matrix of appropriate dimension $\Theta = \Theta^\top$ such that:

$$ \forall s \in \mathbb{C}^+, \quad \begin{bmatrix} I \\ \nabla(s) \end{bmatrix}^* \Theta \begin{bmatrix} I \\ \nabla(s) \end{bmatrix} \preceq 0, \qquad (19) $$

$$ \left[ \mathcal{E} - \mathcal{A} \right]^{\top} \Theta \left[ \mathcal{E} - \mathcal{A} \right]^{\perp} \succeq 0, \qquad (20) $$

where $\nabla^*$ is the trans-conjugate of $\nabla$ and $\mathcal{E}$ is full-rank.

The previous theorem has been adapted to our system considering $\nabla$ is an operator depending only on the Laplace variable $s$ such that $\mathbb{W} = \{\nabla(s)|s \in \mathbb{C}^+\}$. In this case, the well-posedness of system (18) implies its input/output stability. Indeed, the well-posedness of system (18) with $\nabla \in \mathbb{W}$ ensures that there is a unique solution and $\omega$ and $z$ are bounded by the exogenous signals $\bar{\omega}$ and $\bar{z}$. Then, there are no poles with a strictly positive real part and consequently, the system is input/output stable.

Now, we can transform the block diagram in Figure 1 into a suitable form to apply Theorem 2. Consider the following signals:

$$ z(t) = \begin{bmatrix} \dot{X}^\top(t) & KX(t) & KX(t-\tau) & K\dot{X}(t) \end{bmatrix}^\top, \quad \omega(t) = \begin{bmatrix} X^\top(t) & KX(t-\tau) & u(1,t) & K(X(t)-X(t-\tau)) \end{bmatrix}^\top, \qquad (21) $$

For the equivalence between (18) and transfer function (10), $\nabla, \mathcal{A}$ and $\mathcal{E}$ are defined as follows:

$$ \forall s \in \mathbb{C}^+, \nabla(s) = \text{diag} \left( s^{-1} I_n, e^{-\tau s}, \delta(s), \delta_0(s) = \frac{1-e^{-\tau s}}{s} \right), \quad \mathcal{E} = \begin{bmatrix} I_n & 0_{n,1} & 0_{n,1} & 0_{n,1} \\ 0_{1,n} & 1 & 0 & 0 \\ 0_{1,n} & 0 & 1 & 0 \\ -K & 0 & 0 & 1 \\ 0_{1,n} & 1 & -1 & 0 \end{bmatrix}, \quad \mathcal{A} = \begin{bmatrix} A \\ K \\ 0_{1,n} \\ 0_{1,n} \\ 0_{1,n} \end{bmatrix} \qquad (22) $$

Note that $\mathcal{E}$ is full column rank.

**Remark 5** The disturbance $\delta$ is related to the neutral part of system (1). The disturbance component $\delta_0$ is related to Jensen inequality (see for instance [17]), a widely used inequality in the analysis of time-delay systems.
We next need to propose a structure for the real-valued separator $\Theta$ such that inequality (19) always holds for $\nabla \in \mathbb{W}$. For $\Theta = \begin{bmatrix} \Theta_{11} & \Theta_{12} \\ \Theta_{12}^{\top} & \Theta_{22} \end{bmatrix}$, the following structure is proposed:

$$ \Theta_{11} = \operatorname{diag}(0_n, -Q, R(1 - \alpha^2)\gamma^2, -\tau^2 S), \quad \Theta_{12} = \operatorname{diag}(-P, 0, -R\gamma, 0), \quad \Theta_{22} = \operatorname{diag}(0_n, Q, R, S), \qquad (23) $$

with $P \in \mathbb{S}_+^n$ and $Q, R, S \in \mathbb{R}^+$. This selection is not new and is used in the examples presented in [25, 41]. With this specific structure, for all $s \in \mathbb{C}^+$, the following holds:

$$ \begin{bmatrix} I \\ \nabla(s) \end{bmatrix}^* \Theta \begin{bmatrix} I \\ \nabla(s) \end{bmatrix} = \operatorname{diag} \{-2P\Re(e(s^{-1})), Q(|e^{-\tau s}|^2 - 1), \underbrace{R((1-\alpha^2)\gamma^2 - 2\gamma\Re(\delta(s)) + |\delta(s)|^2)}_{\delta_{-1}(s)}, S(|\delta_0(s)|^2 - 1)\} \qquad (24) $$

The third diagonal block can be written differently: $\delta_{-1}(s) = R(|\delta(s) - \gamma|^2 - \alpha^2\gamma^2)$. Equation (22) implies that $\delta$ is inside a circle, its center is $\gamma = \frac{1+\alpha}{1-\alpha^2}$ and its radius is $\frac{(1+\alpha)|\alpha|}{1-\alpha^2} = |\alpha|\gamma$, guaranteeing that $\delta_{-1}(s) \le 0$.

Noticing also that $\forall s \in \mathbb{C}^+, |e^{-\tau s}| \le 1$, $\left|\frac{1-e^{-\tau s}}{s}\right| \le \tau$, we get inequality (19). All these considerations lead to the following stability theorem:

**Theorem 3** If there exist $P \in \mathbb{S}_+^n$ and $Q, R, S \in \mathbb{R}^+$ such that LMI (20) holds for $\Theta$ defined in (23), then system (1) is input-output stable.

## 5.3 Quadratic Separation - Extended stability analysis

### 5.3.1 Motivations and main theorem

The studies of time-delay systems during the last few years focused on reducing the conservatism of stability theorems (see [2] for example). Indeed, it has been showed in [2, 43, 8] that Jensen's inequality leads to conservative results. The idea developed here is to enrich $\Omega$ and $Z$ wisely to improve Theorem 3. Following this idea, we aim at capturing the infinite dimensional behavior of this system, described in equation (7), by adding new signals into the framework.

Indeed, quadratic separation clearly shows that the addition of more information will lead to a smaller kernel of $[\mathcal{E} - \mathcal{A}]$ followed by an improvement of the stability criterion at a price of a more complex $\mathcal{A}$, $\mathcal{E}$ and $\nabla$. Following the methodology described in [43], the new signals are projections of the infinite dimensional state $u$. This state is projected onto the orthogonal basis of shifted Legendre polynomials $\{\mathcal{L}_k\}_{k \in [0,N]}$. Some useful properties of these polynomial are reminded in the sequel. For more information, the reader can refer to [10].

**Theorem 4** For a given $N \in \mathbb{N}$, if there exist $P_N \in \mathbb{S}_+^{n+N}$ and $Q, R, S \in \mathbb{R}^+$ such that the LMI:

$$ \begin{bmatrix} \mathcal{E}_N & -\mathcal{A}_N \end{bmatrix}^{\perp \top} \begin{bmatrix} \Theta_{N,1} & \Theta_{N,2} \\ \Theta_{N,2}^{\top} & \Theta_{N,3} \end{bmatrix} \begin{bmatrix} \mathcal{E}_N & -\mathcal{A}_N \end{bmatrix}^{\perp} \succeq 0, \qquad (25) $$

holds for

$$ \Theta_{N,1} = \operatorname{diag}(0_{n+N}, -Q, R(1 - \alpha^2)\gamma^2, -\tau^2 S), \quad \Theta_{N,2} = \operatorname{diag}(-P_N, 0, -R\gamma, 0_{1,N+1}), \quad \Theta_{N,3} = \operatorname{diag}(0_{n+N}, 0, 0, 0) \qquad (26) $$
<table><caption>Table 1: Number of variables for QS and LS for $n = 4$ and an order $N$.</caption><thead><tr><th></th><th>Lyapunov Functional</th><th>Quadratic Separation</th></tr></thead><tbody><tr><td></td><td><div style="text-align: center;"><div style="border-bottom: 1px solid black; padding-bottom: 2px;">n<sup>2</sup> + n</div><div style="padding-top: 2px;">2</div></div> + 2(N<sup>2</sup> + n) + 5N + Nn + 9</td><td><div style="text-align: center;"><div style="border-bottom: 1px solid black; padding-bottom: 2px;">n<sup>2</sup> + n + N<sup>2</sup> + N</div><div style="padding-top: 2px;">2</div></div> + Nn + 3</td></tr><tr><td>N = 0</td><td>27</td><td>13</td></tr><tr><td>N = 2</td><td>61</td><td>24</td></tr><tr><td>N = 5</td><td>142</td><td>48</td></tr></tbody></table>

$$ \mathcal{E}_N = \begin{bmatrix} I_n & 0_{n,N} & 0_{n,1} & 0_{n,1} & 0_{n,1} \\ 0_{N,n} & I_N & 0_{N,1} & 0_{N,1} & 0_{N,1} \\ 0_{1,n} & 0_{1,N} & 1 & 0 & 0 \\ 0_{1,n} & 0_{1,N} & 0 & 1 & 0 \\ -K & 0_{1,N} & 0 & 0 & 1 \\ 0_{N+1,n} & 0_{N+1,N} & 1_{N+1,1} & -1_{N+1} & 0_{N+1,1} \end{bmatrix}, \quad \mathcal{A}_N = \begin{bmatrix} A & 0_{n,N} & 0_{n,1} & B & 0_{n,N+1} \\ 0_{N,n} & 0_N & 0_{N,1} & 0_{N,1} & \tilde{I}_N(1:N,;) \\ K & 0_{1,N} & 0 & 0 & 0_{1,N+1} \\ 0_{1,n} & 0_{1,N} & 1 & 0 & 0_{1,N+1} \\ 0_{1,n} & 0_{1,N} & 0 & 0 & 0_{1,N+1} \\ 0_{N+1,n} & L_{N(\cdot,1:N)} & 0_{N+1,1} & 0_{N+1,1} & \tilde{I}_N \end{bmatrix}, \qquad (27) $$

with

$$ \begin{align} 1_N &= \begin{bmatrix} (-1)^0 & \dots & (-1)^k & \dots & (-1)^{N-1} \end{bmatrix}^T, & \tilde{I}_N &= \operatorname{diag} \left( \left\{ \frac{1}{\sqrt{2k+1}} \right\}_{k \in [0,N]} \right), \notag \\ L_N &= [\ell_{ij}]_{i,j \in [0,N]}, & \ell_{ik} &= \begin{cases} 0, & \text{if } k \ge i, \\ (2k+1)(1-(-1)^{k+i})c, & \text{otherwise,} \end{cases} \tag{28} \end{align} $$

then system (1) is input/output stable.

**Remark 6** The case $N = 0$ leads to Theorem 3. Theorem 4 also introduces a hierarchy of stability conditions. In other words, if system (1) is proven to be exponentially stable using LMI (25) for a given $N = N_0$, then for all $N \ge N_0$, LMI (25) also assesses the same stability[43].

Compared to classical stability analysis approaches using the Lyapunov Stability (LS) obtained with the Lyapunov functional of [3], this method also uses an LMI solver. It usually results in LMIs with less decision variables than other techniques, which is critical when handling high dimensional systems. Indeed, the difference between LS and Quadratic Separation (QS) is the number of variables. This is due to the two dimensions state extension in LS. The double state extension leads to an increase of the number of variables and then to a slower computation. Table 1 and Figure 3 show that considering $N = 2$ with LS has more decision variables than considering QS at $N = 5$.
<figure><img src="image_4.png" /><figcaption>Figure 3: Number of variables for QS and LS for $n = 4$ and an order $N$.</figcaption></figure>

### 5.3.2 Proof of Theorem 4

The proof of the theorem is divided into 2 parts. First, we will introduce some tools used in the sequel and then the proof itself is derived.

**Legendre Polynomials and frequency-bessel inequality** Before going further, the tools used in the sequel are introduced. For the projections, the scalar product is the canonical inner product of the Hilbert space $L^2([-\tau, 0], \mathbb{C})$ and is denoted $\langle f, g \rangle = \int_{-\tau}^{0} f^*(\theta)g(\theta)d\theta$ for $f, g \in L^2([-\tau, 0], \mathbb{C})$. $\| \cdot \|_2$ is the norm given by the previous inner-product. The shifted Legendre polynomials are defined as follows:

$$ \mathcal{L}_k(u) = (-1)^k \sum_{l=0}^{k} (-1)^l \binom{k}{l} \binom{k+l}{l} \left(\frac{u+\tau}{\tau}\right)^l. \qquad (29) $$

Legendre polynomials have been chosen because of their interesting properties (see [43]). First, because the evaluation of each Legendre polynomial at its boundaries $0$ and $-\tau$ is simple. Secondly, as it is a polynomial basis, the following differentiation rule applies for $N \in \mathbb{N}$, $x \in (-\tau, 0)$:

$$ [\dot{\mathcal{L}}_0(x) \dots \dot{\mathcal{L}}_N(x)]^\top = L_N [\mathcal{L}_0(x) \dots \mathcal{L}_N(x)]^\top, \qquad (30) $$

where $L_N$ is given in (28). Finally, the Legendre polynomials family is orthogonal with respect to the inner product $\langle \cdot, \cdot \rangle$ and it is possible to use Bessel inequality. This inequality is the key part to build a separator and a version adapted to this problem is proposed below:

**Lemma 1** Let $\tau > 0$ and $N \in \mathbb{N}$, then the following inequality holds:

$$ \forall s \in \mathbb{C}^{+}, \quad \delta_{N}^{*}(s)\delta_{N}(s) \le \tau^{2}, \qquad (31) $$

$$ \text{with } \delta_N(s) = \sqrt{\tau} \left[ \left\langle e^{\theta s}, \frac{\mathcal{L}_0(\theta)}{\|\mathcal{L}_0\|_2} \right\rangle \dots \left\langle e^{\theta s}, \frac{\mathcal{L}_N(\theta)}{\|\mathcal{L}_N\|_2} \right\rangle \right]^{\top}. $$

**Proof:** Let $s \in \mathbb{C}^{+}$, Bessel inequality applied to function $\theta \mapsto \exp(\theta s)$ gives: $\sum_{k=0}^{N} \frac{1}{\|\mathcal{L}_k\|_2^2} |\langle e^{\theta s}, \mathcal{L}_k(\theta) \rangle|^2 \le \|e^{\theta s}\|_2^2$.

An identification with $\delta_N$ leads to: $\delta_N^*(s)\delta_N(s) \le \tau |\int_{-\tau}^0 e^{\theta(s+s^*)} d\theta|$. As $s \in \mathbb{C}^{+}$, the right hand side is bounded by $\tau^2$ and that ends the proof. $\square$

**Quadratic Separation formulation of system (1)** In order to ease the comparison with time-delay systems, we introduce a new variable $\theta$ defined as follows: $\theta = -\tau x = -\frac{x}{c}$. If the neutral part is
not taken into consideration, a new infinite dimensional state can be defined: $\tilde{U}(\theta, s) = U(-c\theta, s)\delta^{-1}(s)$. Using equation (7), its expression is then:

$$ \tilde{U}(\theta, s) = \frac{e^{\theta s} + \alpha e^{-(\theta+2\tau)s}}{1 + \alpha} KX(s), \qquad (32) $$

with $t > 0$ and $\theta \in (-\tau, 0)$. There are several starting points to enrich $\Omega$, and the one we choose in this paper is to consider only the projection of the first part of the infinite dimensional state, i.e. $\theta \mapsto e^{s\theta} KX(s)$. We justify our choice in Remark 8.

Let $N \ge 0$ be the number of projections used. The purpose is now to build $\Omega_N$ and $Z_N$ such that:

$$ \begin{aligned} \Omega_N(s) &= \nabla_N(s)Z_N(s), \\ \mathcal{E}_N Z_N(s) &= \mathcal{A}_N \Omega_N(s), \end{aligned} \qquad (33) $$

with $\mathcal{E}_N$ full column-rank where $\Omega_N$ and $Z_N$ are extended versions of $\Omega$ and $Z$. The projections for $k \in [0, N]$ are introduced as follows:

$$ \begin{aligned} \chi_k(s) &= \langle e^{\theta s}, \mathcal{L}_k(\theta) \rangle KX(s), & \chi_N(s) &= \begin{bmatrix} \chi_0(s) & \cdots & \chi_{N-1}(s) \end{bmatrix}^{\top}, \\ \nu_k(s) &= s\sqrt{2k+1}\chi_k(s), & \mathcal{V}_N(s) &= \begin{bmatrix} \nu_0(s) & \cdots & \nu_N(s) \end{bmatrix}^{\top}. \end{aligned} \qquad (34) $$

The new state in Laplace domain is composed of $X(s)$ and the projections of $\theta \mapsto KX(s)e^{\theta s}$ on the orthogonal basis of Legendre polynomials using the projections defined in equation (34). Noting that $\mathcal{L}_k(-\tau) = (-1)^k$ and $\mathcal{L}_k(0) = 1$, the derivation rule in equation (30) and an integration by parts give the following result:

$$ s\chi_k(s) = \left(1 - (-1)^k e^{-\tau s} - \int_{-\tau}^{0} e^{\theta s} \mathcal{L}'_k(\theta) d\theta\right) KX(s) = \left(1 - (-1)^k e^{-\tau s}\right) KX(s) - \sum_{i=0}^{k-1} \ell_{ik}\chi_i(s), \qquad (35) $$

for $k \in [0, N-1]$ and $\ell_{ik}$ as defined in equation (28). Noting that $||\mathcal{L}_k||^2 = \tau(2k+1)^{-1}$, $\mathcal{V}_N$ can be expressed as follows: $\mathcal{V}_N(s) = \delta_N(s) s KX(s)$. The new state vector for the quadratic separation are:

$$ z_N(t) = \begin{bmatrix} \dot{X}^\top(t) & \dot{X}_N^\top(t) & KX(t) & KX(t-\tau) & K\dot{X}(t) \end{bmatrix}^\top, \qquad \omega_N(t) = \begin{bmatrix} X^\top(t) & X_N^\top(t) & KX(t-\tau) & u(1,t) & V_N^\top(t) \end{bmatrix}^\top \qquad (36) $$

With the previous signals, equalities (33) hold for $\mathcal{A}_N$, $\mathcal{E}_N$ defined in the theorem and $\nabla_N(s) = \text{diag}(s^{-1}I_{n+N}, e^{-\tau s}, \delta(s), \delta_N(s))$. Now, we need to find a separator $\Theta_N$. The state extension is based on projections on an orthogonal basis such that the Bessel inequality holds. The result presented in Lemma 1 guarantees that $\Theta_N$ proposed in the theorem is a solution to LMI (19).

**Remark 7** One of the main problem of this method comes from the inclusion of $\delta$ into a disk. This may not be a convenient bound for systems with a high reflexion coefficient $\alpha$. So for $c$ and $c_0$ smaller or $c$ and $c_0$ larger, weaker results are expected.

**Remark 8** ((Another state extension)) We justify the choice of projecting only part of the state $\tilde{U}(\cdot, t)$. First, in equation (32), the state $\tilde{U}$ at a given $\theta$ and $t$ is made up of two contributions. The first one is a result of a wave going forward $\tilde{U}_1(\theta, s) = e^{s\theta} KX(s)$ and another one going backward $\tilde{U}_2(\theta, s) = e^{-s(\theta+2\tau)} KX(s)$. We decided here to project only $\tilde{U}_1$. Another option is to consider the two components $\tilde{U}_1$ and $\tilde{U}_2$ independently and, then, to enrich the state by two projections at each order. This leads to an increased number of variables and, unfortunately, similar performances.

Indeed, to keep the same number of variables, one solution would be to project the whole state $U(x, s)$ on the same basis of Legendre polynomials. Considering $\theta = -\tau x$, we get the following equations for $\theta \in (-\tau, 0)$:

$$ \tilde{U}(\theta, s) = U(\theta, s)\delta^{-1}(s) = \left(e^{s\theta} + \alpha e^{-s(\theta+2\tau)}\right) U(0, s). \qquad (37) $$
Using the dot product defined earlier, we get:

$$ \frac{\chi_k(s)}{U(0, s)} = \left\langle \frac{\tilde{U}(\theta, s)}{U(0, s)}, \mathcal{L}_k(\theta) \right\rangle = \langle e^{s\theta}, \mathcal{L}_k(\theta) \rangle + \alpha e^{-s\tau} \int_{-\tau}^{0} e^{-s(\theta+\tau)} \mathcal{L}_k(\theta) d\theta = \langle e^{s\theta}, \mathcal{L}_k(\theta) \rangle + \alpha e^{-s\tau} \langle e^{s\theta}, \mathcal{L}_k(-\theta - \tau) \rangle. \quad (38) $$

Noting that the shifted Legendre polynomials are symmetric/antisymmetric relative to $-\frac{\tau}{2}$ which means we have $\mathcal{L}_k(-\theta - \tau) = (-1)^k \mathcal{L}_k(\theta)$ and hence we get:

$$ \chi_k(s) = \langle e^{s\theta}, \mathcal{L}_k(\theta) \rangle (1 + (-1)^k \alpha e^{-s\tau}) U(0, s). \quad (39) $$

The problem comes from the Bessel-like inequality in (31) and now reads:

$$ \sum_{k=0}^{N} (2k + 1) \left| \frac{\chi_k(s)}{U(0, s)} \right|^2 \le \tau^2 (1 + |\alpha|)^2. \quad (40) $$

This inequality is not an optimal bound because it will not become an equality for $N \to \infty$, which was the case previously.

# 6 Examples

In this part, we aim at presenting the difference between the exact stability area and the one obtained using robust stability theorems. The same system has also been studied in [3] using a Lyapunov-based stability approach. It also uses a hierarchy of LMI conditions, so the notation “LS, $N = i$” refers to the stability obtained using the other methodology for an order $i$. A comparison of efficiency between all the methods is also presented.

The estimation of the stability area of system (1) is provided on some chosen examples with different behaviors. We consider the interconnection of a stable wave equation ($c_0 > 0$) interconnected with different ODE systems. First, the interconnection with a stable finite dimension LTI system is proposed. Then, an unstable ODE is interconnected. We then aim at proving that there exist unstable systems stabilized thanks to the wave equation. To finish, we analyze a system known to possess stability pockets for $cc_0 = 1$. The LMI solver used in the examples is “sdpt3” with Yalmip[31].

## 6.1 A and $A + BK$ Hurwitz with $\|H\|_{\infty} < 1$

This first example is borrowed from [3] where the $K$ matrix has been slightly modified to get $\|H\|_{\infty} < 1$ such that the Small-Gain Theorem can be applied. System (1) is proposed with the following matrices:

$$ A = \begin{bmatrix} -2 & 1 \\ 0 & -1 \end{bmatrix}, \quad B = \begin{bmatrix} 1 \\ 1 \end{bmatrix}, \quad K = \begin{bmatrix} 0 & -\frac{20}{21} \end{bmatrix}. \quad (41) $$

It is easy to verify that $A$ and $A + BK$ are indeed Hurwitz. Moreover, the infinity norm of the open-loop system is less than 1, then Theorem 1 applies and $\lim_{c_0 \to \infty} c_{min}(c_0) = 0$. Here we are after computing the $c_{min}$, which is the smallest $c$ such that for all $c > c_{min}$, the system is input/output stable. While for the Small-Gain theorem, we have proven this bound, the remaining approaches will reveal the bounds through computations.

The results with quadratic separation (QS) and the method with Lyapunov-based stability (LS) at order $N = 0$ developed in [3] are superposed in the chart of Figure 4. As expected, the small gain theorem provides the worst estimation of the stability area but with stronger properties. The QS and LS approaches provide similar results. The QS does not use an extended state so the results for high $c_0$ are further from the CTCR curve than the one obtained with LS.

Finally, the result obtained with CTCR shows a non-continuous behavior for small $cc_0$ and this observation discourages us to use the $(c_0, c)$ chart to estimate the stability area of system (1). This is not a numerical issue, but rather indicates some singularities regarding system poles crossing the imaginary axis for certain
<figure><img src="image_5.png" /><figcaption>Figure 4: Minimum wave speed $c_{min}$ as a function of $c_0$ for system (1) to be stable. The values for A, B and K are given by equation (41).</figcaption></figure>

values of $c_0$. This can be remedied by viewing the stability on a different domain of parameters. Indeed, an appropriate choice of system coordinates is $(cc_0, c^{-1})$. In order to make a comparison with time-delay systems, the delay $c^{-1}$ needs to be considered as one of the axes. It is natural to add $cc_0$ as the other variable of interest as it pilots the behavior of the system. Then, on the line $cc_0 = 1$, we get a time-delay system (see the state-space representation (12) and Remark 3), and it is possible to compare the results with the literature. From this point, the mapping $(cc_0, c^{-1})$ is used.

## 6.2 A and $A + BK$ are not Hurwitz

The following system is borrowed from [2] where

$$A = \begin{bmatrix} 0 & 1 \\ -2 & 0.1 \end{bmatrix}, \quad B = \begin{bmatrix} 0 \\ 1 \end{bmatrix}, \quad K = \begin{bmatrix} 1 & 0 \end{bmatrix}. \qquad (42)$$

It has been shown in the cited study that for $c_1 = cc_0 = 1$, it is stable for a sufficiently large delay. This is a very different theoretical case compared to the previous one. For $\tau = 0$, uncontrolled system has two unstable poles, which means that a simple output feedback cannot stabilize the system. This is evident from both A and $A + BK$ being unstable. Notice first that this system for $cc_0 = 1$ has been studied in [2], and CTCR recovers the same exact stability area.

QS at $N = 0$ does not provide meaningful results, but, as we can see in Figure 5, higher orders of QS detect a stability area which approaches the exact results by CTCR as the order increases; see the stability regions around $c_1 = 1$. Despite the state augmentation, it is not possible to recover the stability area and computations show that for $cc_0 < 0.5$, QS at any order does not assess stability. This can be a consequence of considering an inappropriate bound on $\delta$.

## 6.3 An example with stability pockets

For this last example, the system is taken from [17]. It is known to possess multiple stable intervals (pockets) along the delay axis for $c_1 = cc_0 = 1$. We investigate whether or not QS and the methodology of [3] can detect these pockets. The system matrices are given by:

$$A = \begin{bmatrix} 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ -11 & 10 & 0 & 0 \\ 5 & -15 & 0 & -0.25 \end{bmatrix}, \quad B = \begin{bmatrix} 0 \\ 0 \\ 1 \\ 0 \end{bmatrix}, \quad K = \begin{bmatrix} 1 \\ 0 \\ 0 \\ 0 \end{bmatrix}^T. \qquad (43)$$
<figure><img src="image_6.png" /><figcaption>Figure 5: Stability areas for system (1) with matrices defined in (42) and $c_1 = cc_0$ and $\tau = c^{-1}$ obtained using CTCR and Theorem 4. The hatched area is the exact unstable area, the remaining area is the exact stable area as detected by CTCR. The colored areas are stated as stable for Theorem 4 at an order $N = 3$. The color scale represents the stability area depending on the order $N$ considered. The areas marked as $S_i$ is the stable area up to an order $i$.</figcaption></figure>

We can compare the efficiency of LS in Figure 6a. The hierarchy property can be seen and the stability pockets are indeed recovered as the order increases. Figure 6b shows the stability result with QS at an order $N$. QS and LS at the same order provide similar results but with lower decision variables for QS, see Table 1 and Figure 3.

Moreover, for small $c_1$, the reflection coefficient is not close to 0 and then, as noted in Remark 7, QS should have a worse estimation. That is why we get poor result on the left side compared to LS. But if $c_1$ is close to 1, results similar to LS are obtained using QS with a subsequently lower number of decision variables. That means at the same number of decision variables, QS can detect more stability pockets around $c_1 = 1$ but LS detects a wider stability area. When $4 \le N \le 7$, the stability area detected by QS is included in the one ensured by LS at the price of more decision variables. Not matter the order, QS never detects the whole stability area.

# 7 Conclusion

In this paper, we studied using an exact method the input/output stability in a parameter region of interest for a system interconnecting an ODE with a string-type PDE equation. We also proposed two robust stability results: a simple one obtained with the Small-Gain theorem and another one using Quadratic Separation. The latter is based on a hierarchy of LMI conditions and is more conservative than some other approaches; but it proposes a subsequently reduced number of decision variables. A perspective would be to enhance the stability results for Quadratic Separation while keeping its low computational burden. Future work will be aimed at determining which approximation introduces the most conservatism in order to reduce the gap between the Lyapunov Stability analysis and the one made with Quadratic Separation.

# References

[1] M Arcak, C Meissen, and A Packard. *Stability of Interconnected Systems*. Springer, 2016.
<figure><img src="image_7.png" /><figcaption>(a) Theorem 2 from [3]: Lyapunov functional</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(b) Theorem 4: Quadratic Separation</figcaption></figure>

Figure 6: Stability areas for system (1) with matrices defined in (43) and $c_1 = cc_0$ and $\tau = c^{-1}$ obtained using CTCR, LS and QS. The hatched area is the exact unstable area, the white area corresponds to the unknown area according to the corresponding theorems but known to be stable with CTCR. The color scale represents the stability area depending on the number of variables considered. To enable comparison, $N$ varies between 0 and 3 for LS and between 0 and 7 for QS. The areas marked as $S_i$ is the stable area up to an order $i$.
[2] Y Ariba, F Gouaisbaut, A Seuret, and D Peaucelle. Stability analysis of time-delay systems via Bessel inequality: A quadratic separation approach. *International Journal of Robust and Nonlinear Control*, 28(5):pp 1507-1527, March 2018.

[3] M Barreau, A Seuret, F Gouaisbaut, and L Baudouin. Lyapunov stability analysis of a string equation coupled with an ordinary differential system. *IEEE Transactions on Automatic Control*, 1, 2018, DOI:10.1109/TAC.2018.2802495.

[4] G. Bastin and J.-M. Coron. *Stability and boundary stabilization of 1-d hyperbolic systems*, volume 88. Springer, 2016.

[5] R E Bellman and K L Cooke. *Differential-difference equations*. Academic Press, 1963.

[6] E Beretta and Y Kuang. Geometric stability switch criteria in delay differential systems with delay dependent parameters. *SIAM Journal on Mathematical Analysis*, 33(5):1144–1165, 2002.

[7] D Bresch-Pietri and M Krstic. Output-feedback adaptive control of a wave PDE with boundary anti-damping. *Automatica*, 50(5):1407–1415, 2014.

[8] C Briat. Convergence and equivalence results for the jensen's inequality - application to time-delay and sampled-data systems. *IEEE Transactions on Automatic Control*, 56(7):1660–1665, 2011.

[9] F Califano, A Macchelli, and C Melchiorri. Stability analysis of repetitive control: The port-hamiltonian approach. (1):1894–1899, 2017.

[10] R Courant and D Hilbert. *Methods of mathematical physics*, volume 1 of *Wiley Classics Library*. John Wiley & Sons, Inc, New York, 1989.

[11] R F Curtain and H J Zwart. *An introduction to infinite-dimensional linear systems theory*, volume 21 of *Texts in Applied Mathematics*. Springer-Verlag, New York, 1995.

[12] G Doetsch. *Introduction to the theory and application of the Laplace transformation*. Springer, 2012.

[13] N Espitia, A Girard, N Marchand, and C Prieur. Event-based control of linear hyperbolic systems of conservation laws. *Automatica*, 70:275–287, 2016.

[14] L C Evans. *Partial Differential Equations*. Graduate studies in mathematics. American Mathematical Society, 2010.

[15] G Flores. Dynamics of a damped wave equation arising from mems. *SIAM Journal on Applied Mathematics*, 74(4):1025–1035, 2014.

[16] P Grabowski and F M Callier. Boundary control systems in factor form: Transfer functions and input-output maps. *Integral Equations and Operator Theory*, 41(1):1–37, Mar 2001.

[17] K Gu, J Chen, and V L Kharitonov. *Stability of time-delay systems*. Springer, 2003.

[18] B Z Guo and C Z Xu. The stabilization of a one-dimensional wave equation by boundary feedback with noncollocated observation. *IEEE Transactions on Automatic Control*, 52(2):371–377, Feb 2007.

[19] J K Hale and S M V Lunel. *Introduction to functional differential equations*. Applied mathematical sciences. Springer-Verlag, New York, Berlin, Heidelberg, 1977.

[20] Hale, J K. *Theory of functional differential equations*. Number vol 3 in Applied Mathematical Sciences Series. Springer Verlag, 1977.

[21] W He and S S Ge. Cooperative control of a nonuniform gantry crane with constrained tension. *Automatica*, 66:146–154, 2016.

[22] W He and S Zhang. Control design for nonlinear flexible wings of a robotic aircraft. *IEEE Transactions on Control Systems Technology*, 25(1):351–357, Jan 2017.

[23] W He, S Zhang, and S S Ge. Adaptive control of a flexible crane system with the boundary output constraint. *IEEE Transactions on Industrial Electronics*, 61(8):4126–4133, Aug 2014.

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

[25] T Iwasaki and S Hara. Well-posedness of feedback systems: insights into exact robustness analysis and approximate computations. *IEEE Transactions on Automatic Control*, 43(5):619–630, May 1998.

[26] C Jin, K Gu, S. I Niculescu, and I Boussaada. Stability analysis of systems with delay-dependent coefficients. *Phd Thesis*, IEEE Access, 2018.
[27] V Kolmanovskii and A Myshkis. *Introduction to the theory and applications of functional differential equations*, volume 463. Springer, 2013.

[28] M Krstic. Dead-time compensation for wave/string pdes. 133:4458-4463, 12 2009.

[29] M Krstic. *Delay compensation for nonlinear, adaptive, and PDE systems*. Springer, 2009.

[30] J Lagnese. Decay of solutions of wave equations in a bounded region with boundary dissipation. *Journal of Differential Equations*, 50(2):163–182, 1983.

[31] J Löfberg. Yalmip : a toolbox for modeling and optimization in matlab. *Computer Aided Control Systems Design*, 2004 IEEE International Symposium, pages 284–289, Sept 2004.

[32] T Louw, S Whitney, A Subramanian, and H Viljoen. Forced wave motion with internal and boundary damping. *Journal of applied physics*, 111:14702–147028, 01 2012.

[33] Z H Luo, B Z Guo, and Ö Morgül. *Stability and stabilization of infinite dimensional systems with applications*. Springer, 2012.

[34] T Meurer and A Kugi. Tracking control design for a wave equation with dynamic boundary conditions modeling a piezoelectric stack actuator. *International Journal of Robust and Nonlinear Control*, 21(5):542–562, 2011.

[35] Ö Morgül. A dynamic control law for the wave equation. *Automatica*, 30(11):1785–1792, 1994.

[36] Ö Morgül. Stabilization and disturbance rejection for the wave equation. *IEEE Transactions on Automatic Control*, 43(1):89–95, 1998.

[37] Ö Morgül. An exponential stability result for the wave equation. *Automatica*, 38(4):731–735, 2002.

[38] S I Niculescu. *Delay effects on stability: a robust control approach*, volume 269. Springer, 2001.

[39] N Olgac and R Sipahi. An exact method for the stability analysis of time-delayed linear time-invariant (LTI) systems. *IEEE Transactions on Automatic Control*, 47(5):793–797, 2002.

[40] N Olgac and R Sipahi. A practical method for analyzing the stability of neutral type LTI-time delayed systems. *Automatica*, 40(5):847–853, 2004.

[41] D Peaucelle, D Arzelier, D Henrion, and F Gouaisbaut. Quadratic separation for feedback connection of an uncertain matrix and an implicit linear transformation. *Automatica*, 43(5):795–804, 2007.

[42] Z. V Rekasius. A stability test for systems with delays. *Proceedings of the joint automatic control conference*, TP9-A, 1980.

[43] A Seuret and F Gouaisbaut. Hierarchy of LMI conditions for the stability analysis of time delay systems. *Systems and Control Letters*, 81:1–7, July 2015.

[44] R Sipahi, S I Niculescu, C T Abdallah, W Michiels, and K Gu. Stability and stabilization of systems with time delay. *IEEE Control Systems*, 31(1):38–65, 2011.

[45] R Sipahi and N Olgac. Degenerate cases in using the direct method. *ASME 2003 International Design Engineering Technical Conferences and Computers and Information in Engineering Conference*, pages 2201–2210, 2003.

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

[47] M Tucnak and G Weiss. *Observation and control for operator semigroups*. Springer, 2009.

[48] A van der Schaft and D Jeltsma. Port-hamiltonian systems theory: An introductory overview. *Foundations and Trends in Systems and Control*, 1(2-3):173–378, 2014.

[49] J Woźniak and B Niesterowicz. Input-output stability analysis of a slowly rotating rotor with friction on one end. *25th Mediterranean Conference on Control and Automation (MED)*, (1):847–851, July 2017.