Lyapunov stability analysis of a string equation coupled with an ordinary differential system

# Lyapunov stability analysis of a string equation coupled with an ordinary differential system

Matthieu Barreau, Alexandre Seuret, Frédéric Gouaisbaut and Lucie Baudouin M. Barreau, A. Seuret, F. Gouaisbaut and L. Baudouin are with LAAS - CNRS, Université de Toulouse, CNRS, UPS, France. This work is supported by the ANR project SCIDiS contract number 15-CE23-0014.
###### Abstract

This paper considers the stability problem of a linear time invariant system in feedback with a string equation. A new Lyapunov functional is proposed using augmented states which enriches and encompasses the classical functionals of the literature. It results in tractable stability conditions expressed in terms of linear matrix inequalities. This methodology follows from the application of the Bessel inequality to the projections over the Legendre polynomials. Numerical examples illustrate the potential of our approach through three scenari: a stable ODE perturbed by the PDE, an unstable open-loop ODE and an unstable closed-loop ODE stabilized by the PDE.

String equation, Ordinary differential equation, Lyapunov functionals, LMI.

## I Introduction

This paper presents a novel approach to assess stability of a heterogeneous system composed of the interconnection of a partial differential equation (PDE), more precisely a damped string equation, with a linear ordinary differential equation (ODE). While the topic of stability and control of PDE systems has a rich literature between applied mathematics [9, 18] and automatic control [20]; the stability analysis (and the control) of such a coupled system belongs to a recent research area. To cite a few related results, one can refer to [5, 6, 25] where an ODE is interconnected with a transport equation, to [27] for a heat equation, [15, 7] for the wave equation and [29] for the beam equation.

Generally, the PDE is viewed as a perturbation to be compensated for instance using a backstepping method as proposed in [16], where infinite dimensional controllers are provided to cope with the undesirable effect of the PDE. Another interesting point of view relies on the converse approach: the ODE system can be seen as a finite dimensional boundary controller for the PDE (see [12, 22, 23]). A last strategy describes a robust control approach, aiming at characterizing the robustness of the PDE-ODE interconnection [13].

In the present paper, we consider a damped string equation, i.e. a stable one-dimensional wave equation which is connected at its boundary to a stable or unstable ODE. The proposed method to assess stability is inspired by the recent developments on the stability analysis of time-delay systems based on Bessel inequality and Legendre polynomials [26]. Since time-delay systems represent a particular class of systems coupling a transport PDE with a classical ODE system (see for instance [3]), the main motivation of this work is to show how this methodology can be adapted to a larger class of PDE/ODE systems as demonstrated with the heat equation in [2].

Compared to the literature on coupled PDE/ODE systems, the proposed methodology aims at designing a new Lyapunov functional, integrating some cross-terms merging the ODE’s and the PDE’s usual terms. This new class of Lyapunov functional encompasses the classical notion of energy usually proposed in the literature by offering more flexibility. Hence, it allows to guarantee stability for a larger set of systems, for instance, unstable open-loop ODE and, for the first time to the best of our knowledge, even an unstable closed-loop ODE; meaning that the PDE helps for the stabilization.

The paper is organized as follows. The next section formulates the problem and provides some general results on the existence of solutions and equilibrium. In Section 3, after a modeling phase inspired by the Riemann coordinates, a generic form of Lyapunov functional is introduced, and its associate analysis leads to a first stability theorem. Then, in Section 4, an extension using Bessel inequality is provided. Finally, Section 5 discusses the results on three examples. The last section draws some conclusion and perspectives.

Notations: is the closed set and . is a multi-variable function from to . The notation stands for . We also use the notations and for the Sobolev spaces: and particularly . The norm in is . For any square matrices , the operations ‘He’ and ‘diag’ are defined as follows: and . A symmetric positive definite matrix of belongs to the set or we write more simply .

## Ii Problem Statement

We consider the coupled system described by

 ˙X(t) =AX(t)+Bu(1,t), t⩾0, (1a) utt(x,t) =c2uxx(x,t), x∈Ω,t⩾0, (1b) u(0,t) =KX(t), t⩾0, (1c) ux(1,t) =−c0ut(1,t), t⩾0, (1d) u(x,0) =u0(x), ut(x,0)=v0(x), x∈Ω, (1e) X(0) =X0, (1f)

with the initial conditions and such that equations (1c) and (1d) are respected. They are then called “compatible” with the boundary conditions. and are time-invariant matrices of appropriate size.

###### Remark 1

When no confusion is possible, parameter may be omitted and so do the domains of definition.

This system can be viewed as an interconnection in feedback between a linear time invariant system (1a) and an infinite dimensional system modeled by a string equation (1b). The latter is a one dimension hyperbolic PDE, representing the evolution of a wave of speed and amplitude . To keep the content clear, is assumed to be a scalar but the calculations are done as if it was a vector of any dimension. The measurement we have access to is the state at which is the right extremity of the string and the control is a Dirichlet actuation (equation (1c)) because it affects directly the state and not its derivative. Another boundary condition must be added. It is defined at by . This is a well-known damping condition when (see for example [17]). As presented in [4], we find this kind of systems for instance when modeling a drilling mechanism. The control is then given at one end and the measurement is done at the other end.

More generally, this system can be seen either as the control of the PDE by a finite dimensional dynamic control law generated by an ODE [8] or on the contrary the robustness of a linear closed loop system with a control signal conveyed by a damped string equation. On the first scenario, both the ODE and the PDE are stable and the stability of the coupled system is studied. The second case corresponds to an unstable but stabilizable ODE connected to a stable PDE. To sum up, this paper focuses on the stability analysis of closed-loop coupled system (1) with a potentially unstable closed-loop ODE but a stable PDE. This differs significantly from the backstepping methodology of [15] which aims at designing an infinite dimensional control law ensuring the stability of a cascaded ODE-PDE system with a closed-loop stable ODE.

### Ii-a Existence and regularity of solutions

This subsection is dedicated to the existence and regularity of solutions to system (1). We first introduce for . We consider the classical norm on the Hilbert space :

 ∥(X,u,v)∥2H=|X|2n+∥u∥2+c2∥ux∥2+∥v∥2.

This norm can be seen as the sum of the energy of the ODE system and the one of the PDE.

###### Remark 2

A more natural norm for space would be which is equivalent to . The norm used here makes the calculations easier in the sequel.

Once the space is defined, we model system (1) using the following linear unbounded operator :

 T(Xuv)=(AX+Bu(1)vc2uxx) and
 D(T)={(X,u,v)∈H2,u(0)=KX,ux(1)=−c0v(1)}.

This operator is said to be dissipative with respect to a norm if its time-derivative along the trajectories generated by is strictly negative. The goal of this paper is then to find an equivalent norm to which allows us to refine the dissipativity analysis of . This equivalent norm is derived from a general formulation of a Lyapunov functional, whose parameters are chosen using a semi-definite programming optimization process.
Beforehand, from the semi-group theory, we propose the following result on the existence of solutions for (1).

###### Proposition 1

If there exists a norm on for which the linear operator is dissipative with non singular, then there exists a unique solution of system (1) with initial conditions . Moreover, the solution has the following regularity property: .

Proof : This proof follows the same lines than in [21]. Applying Lumer-Phillips theorem (p103 from [28]), as the norm is dissipative, it is enough to show that for all with , where is the range operator. Let , we want to show that for this system, there exists for which the following set of equation is verified:

 λX−AX−Bu(1) =r, (2a) λu(x)−v(x) =g(x), (2b) λv(x)−c2uxx(x) =h(x), (2c)

for all and a given . Equations (2b), (2c) give:

 ∀x∈(0,1),  u(x)=k1exp(λc−1x)+k2exp(−λc−1x)+G(x)

where . are constants to be determined. Using the boundary condition , we get:

 ∀x∈(0,1),  u(x)=2k1sinh(λc−1x)+KXe−λc−1x+G(x),

Taking its derivative at the boundary we get:

 ux(1)=2λc−1k1cosh(λc−1)−λc−1KXe−λc−1+G1,

with known. We also have , leading to with and . Then using (2a), we get:

 (λIn−(A+BKf(λc−1)))X=r+BG2.

Since when and is non singular, there exists such that is non singular and

 ∀λ∈(0,λmax),det(λIn−(A+BKf(λc−1)))≠0.

Then, there is a unique for a given . We immediately get that is in . Then for , . The regularity property falls from Lumer-Phillips theorem.

### Ii-B Equilibrium point

An equilibrium of system (1) is such that , i.e. it verifies the following linear equations:

 0 =AXe+Bue(1), (3a) 0 =c2∂xxue(x),x∈(0,1), (3b) ve(x) =0,x∈(0,1), (3c) ue(0) =KXe, (3d) ∂xue(1) =0. (3e)

Using equation (3b), we get as a first order polynomial in but in accordance to equation (3e), is a constant function. Then, using equation (3d), we get . That leads to: . We obtain the following proposition:

###### Proposition 2

An equilibrium of system (1) verifies , , . Moreover, if is not singular, system (1) admits a unique equilibrium .

## Iii A First Stability Analysis Based on Modified Riemann Coordinates

This part is dedicated to the construction of a Lyapunov functional. We introduce therefore a new structure based on variables directly related to the states of system (1).

### Iii-a Modified Riemann coordinates

The PDE considered in system (1) is of second order in time. As we want to use some tools already designed for first order systems, we propose to define some new states using modified Riemann coordinates, which satisfy a set of coupled first order PDEs and diagonalize the operator. Let us introduce these coordinates, defined as follows:

 χ(x)=[ut(x)+cux(x)ut(1−x)−cux(1−x)]=[χ+(x)χ−(1−x)].

The introduction of such variables is not new and the reader can refer to articles [24, 4] or [10] and references therein about Riemann invariants. and are eigenfunctions of equation (1b) associated respectively to the eigenvalues and . Therefore, using , the previous equation leads to a transport PDE:

 ∀t≥0,∀x∈Ω,χt(x,t)=cχx(x,t). (4)
###### Remark 3

The norm of the modified state can be directly related to the norm of the functions and . Indeed simple calculations and a change of variables give:

 ∥χ∥2=2(∥ut∥2+c2∥ux∥2). (5)
###### Remark 4

This manipulation does not aim at providing an equivalent formulation for system (1) but at identifying a manner to build a Lyapunov functional for system (1).

The second step is to understand how the extra-variable interacts with ODE (1a). Hence using (1c), we notice:

 ˙X=AX+B(u(1)−u(0)+KX),=(A+BK)X+B∫10ux(x)dx.

To express the last integral term using , we note that:

 2c∫10ux(x)dx=∫10χ+(x)dx−∫10χ−(x)dx.

This expression allows us to rewrite the ODE system as where and . The extra-state follows the dynamics:

 ˙X0=c∫10χx(x)dx=c[χ(1)−χ(0)].

The ODE dynamic can then be enriched by considering an extended system where is viewed as a new dynamical state:

 ˙X0=[A+BK~B02,n02]X0+[0n,2cI2](χ(1)−χ(0)), (6)

with . Hence, associated to the original system (1), we propose a set of equation (4)-(6). They are linked to system (1) but enriched by extra dynamics aiming at representing the interconnection between the extended finite dimensional system and the two transport equations. Nevertheless, these two systems are not equivalent. The transport equation gives trajectories of and but can be defined within a constant. The second set of equations just induces a formulation for a Lyapunov functional candidate which is developed in the subsection below.

### Iii-B Lyapunov functional and stability analysis

The main idea is to rely on the auxiliary variables satisfying equations (4) and (6) to define a Lyapunov functional for the original system (1). The associated Lyapunov function of ODE (6) is a simple quadratic term on the state , with . It introduces automatically a cross-term between the ODE and the original PDE through . Hence, the auxiliary equations of the previous paragraph shows a coupling between a finite dimensional LTI system and a transport PDE. For the latter, inspired from the literature on time-delay systems [3, 10], we provide a Lyapunov functional:

 V(u)=∫10χ⊤(x)(S+xR)χ(x)dx,

with . The use of the modified Riemann coordinates enables us to consider full matrices and . As the transport described by the variable is going backward, is multiplied by . Thereby, we propose a Lyapunov functional for system (1) expressed with the extended state variable :

 V0(X0,u)=X⊤0P0X0+V(u). (7)

This Lyapunov functional is actually made up of three terms:

• A quadratic term in introduced by the ODE;

• A functional for the stability of the string equation;

• A cross-term between and described by the extended state .

The idea is that this last contribution is interesting since we may consider the stability of system (1) with an unstable ODE, stabilized thanks to the string equation. At this stage, a stability theorem can be derived using the Lyapunov functional .

###### Theorem 1

Consider the system defined in (1) with a given speed , a viscous damping with initial conditions . Assume there exist and such that the following LMI holds:

 Ψ0=He(Z⊤0P0F0)−c~R0+c(H⊤0(S+R)H0−G⊤0SG0)≺0 (8)

where

 (9)

Then, there exists a unique solution to system (1) and it is exponentially stable in the sense of i.e. there exist such that the following estimate holds for :

 ∥(X(t),u(t),ut(t))∥2H⩽γe−δt∥(X0,u0,v0)∥2H. (10)
###### Remark 5

The LMI includes a necessary condition given by with , which is . This inequality is guaranteed if and only if the matrix has its eigenvalues inside the unit cycle of the complex plan, i.e. , which is consistent with the result on exponential stability of [14].

### Iii-C Proof of Theorem 1

The proof of stability is presented below.

#### Iii-C1 Preliminaries

As a first step of this proof, an inequality on is presented below.

###### Lemma 1

For , the following inequality holds:

 ∥u∥2⩽2∥ux∥2+2|u(0)|2.

Proof : Since , Young and Jensen inequalities imply that for all :

 u(x)2=(∫x0us(s)ds+u(0))2⩽2∫x0u2s(s)ds+2|u(0)|2.

The proof of Theorem 1 consists in explaining how the LMI condition presented in the statement implies that there exist a functional and three positive scalars and such that the following inequalities hold:

 ε1∥(X,u,ut)∥2H⩽V(X,u)⩽ε2∥(X,u,ut)∥2H,˙V(X,u)⩽−ε3∥(X,u,ut)∥2H. (11)

The next steps aim at proving (11) in order to obtain the convergence of the state to the equilibrium.

#### Iii-C2 Well-posedness

If the conditions of Theorem 1 are satisfied, then the inequality holds where . After some simplifications, we get , for some matrix depending on , and . This strict inequality requires that is non singular and, in light of Propositions 1 and 2, the problem is indeed well-posed and is the unique equilibrium point. Furthermore, note that since is not necessarily symmetric, then matrix does not have to be Hurwitz.

#### Iii-C3 Existence of ε1

Conditions and mean that there exists , such that for all :

 P0⪰ε1diag(In+2K⊤K,02),S+xR ⪰ S⪰ε12+c22c2I2.

 V0(X0,u)⩾ε1(|X|2n+|KX|2+2+c22c2∥χ∥2)+∫10χ⊤(x)(S+xR−ε12+c22c2I2)χ(x)dx⩾ε1(|X|2n+|KX|2+2+c22c2∥χ∥2).

Using boundary condition (1c) and equality (5), it becomes

 V0(X0,u)⩾ε1(|X|2n+∥u∥2+∥ut∥2+c2∥ux∥2)+2ε1c2∥ut∥2+ε1(2∥ux∥2+2|u(0)|2−∥u∥2).

Then, we apply Lemma 1 to ensure that the last term is positive. It follows that , which ends the proof of existence of .

#### Iii-C4 Existence of ε2

Since and , there exists such that for :

 P0⪯diag(ε2In,ε24I2),S+xR ⪯ S+R⪯ε24I2.

From equation (7), we get:

 V0(X0,u)⩽ε2(|X|2n+14X⊤0X0+14∥χ∥2)+∫10χ⊤(x)(S+xR−ε24I2)χ(x)dx⩽ε2(|X|2n+12∥χ∥2) (12)

where we have used , which is a result of Jensen’s inequality. The proof of the existence of ends by using (5) so that we get:

 V0(X0,u)⩽ε2(|X|2n+∥ut∥2+c2∥ux∥2)⩽ε2∥(X,u,ut)∥2H.

#### Iii-C5 Existence of ε3

Differentiating in (7) along the trajectories of system (1) leads to

Our goal is to express an upper bound of thanks to the extended vector defined as follows:

 ξ0=[X⊤X⊤0ut(1)cux(0)]⊤. (13)

Let us first concentrate on . Equation (4) yields:

 ˙V(u)=2c∫10χ⊤x(x,t)(S+xR)χ(x,t)dx. (14)

Integrating by parts the last expression leads to:

 ˙V(u)=c(χ⊤(1)(S+R)χ(1)−χ⊤(0)Sχ(0)∫10χ⊤(x)Rχ(x)dx−∫10χ⊤(x)Rχ(x)dx). (15)

Then we note that , , , with defined in (13) and the matrices above in (9). We get and which results in the following expression for :

 ˙V0(X0,u)=ξ⊤0(He(Z⊤0P0F0)+cH⊤0(S+R)H0−cG⊤0SG0)ξ0−c∫10χ⊤(x)Rχ(x)dx. (16)

Then, using the definition of given in (8), the previous expression can be rewritten as follows:

 ˙V0(X0,u)=ξ⊤0Ψ0ξ0+cX⊤0RX0−c∫10χ⊤(x)Rχ(x)dx. (17)

Since and , there exists such that:

 R⪰ ε32c2+c2c2I2, (18a) Ψ0⪯ −ε3diag(In+2K⊤K,2+c22c2I2,02). (18b)

Using (18b) and the boundary condition , equation (17) becomes:

 ˙V0(X0,u)⩽−ε3(|X|2n+2|u(0)|2+2+c22c2∥χ∥2)+cX⊤0(R−ε32c2+c2c2I2)X0−c∫10χ⊤(x)(R−ε32c2+c2c2I2)χ(x)dx,

so that we get by application of Jensen’s inequality:

 ˙V0(X0,u)⩽−ε3(|X|2n+2|u(0)|2+2+c22c2∥χ∥2). (19)

The most important part of the proof lies in the following trick. Since (5) holds, we get:

 ˙V0(X0,u)⩽−ε3∥(X,u,ut)∥2H−ε32c2∥ut∥2−ε3(2|u(0)|2+2∥ux∥2−∥u∥2).

Moreover, Lemma 1 ensures that the last term of the previous expression is negative so that we have , which concludes this proof of existence.

#### Iii-C6 Conclusion

Finally, there exist such that (11) holds for a functional . Hence defines an equivalent norm to and is dissipative. It means, according to Propositions 1 and 2, that there exists a unique solution to system (1) in . Equation (11) also brings: and

 ∀t>0,∥(X(t),u(t),ut(t))∥2H⩽ε2ε1e−ε3ε2t∥(X0,u0,v0)∥2H,

which shows the exponential convergence of all the trajectories of system (1) to the unique equilibrium . In other words, the solution to system (1) is exponentially stable.

###### Remark 6

It is also worth noting that LMI (8) can be transformed to extend this theorem to uncertain ODE systems subject to polytopic-type uncertainties for instance.

## Iv Extended Stability Analysis

In the previous analysis, we have proposed an auxiliary system presented in (4)-(6) helping us to define a new Lyapunov functional for system (1). The notable aspect is that the term appears naturally in the dynamics of system (1). In light of the previous work on integral inequalities in [26], this term can also be interpreted as the projection of the modified state over the set of constant functions in the sense of the canonical inner product in . One may therefore enrich (6) by additional projections of over the higher order Legendre polynomials, as one can read in [26, 3] in the context of time-delay systems. The family of shifted Legendre polynomials, denoted and defined over by with , form an orthogonal family with respect to the inner product (see [11] for more details).

### Iv-a Preliminaries

The previous discussion leads to the definition of the projection of any function in on the family :

 ∀k∈N,Xk:=∫10χ(x)Lk(x)dx,

An augmented vector is naturally derived for any :

 XN=[X⊤X⊤0⋯X⊤N]⊤. (20)

Following the same methodology as in Theorem 1, this specific structure suggests to introduce a new Lyapunov functional, inspired from (7), with :

 VN(XN,u)=X⊤NPNXN+V(u). (21)

In order to follow the same procedure, several technical extensions are required. Indeed, the stability conditions issued from the functional are proved using Jensen’s inequality and an explicit expression of the time derivative of . Therefore, it is necessary to provide an extended version of Jensen’s inequality and of this differentiation rule. These technicals steps are summarized in the two following lemmas.

###### Lemma 2

For any function and symmetric positive matrix , the following Bessel-like integral inequality holds for all :

 ∫10χ⊤(x)Rχ(x)dx⩾N∑k=0(2k+1)X⊤kRXk. (22)

This inequality includes Jensen’s inequality as the particular case , suggesting that this lemma is an appropriate extension and should help to address the stability analysis using the new Lyapunov functional (21) with the augmented state defined in (20).

The proof of Lemma 2 is based on the expansion of the positive scalar where can be interpreted as the approximation error between and its orthogonal projection over the family .

The next lemma is concerned by the differentiation of .

###### Lemma 3

For any function , the following expression holds for any in :

where

 (23)

with if and otherwise.

Proof : The proof of this lemma is presented in appendix because of its technical nature.

### Iv-B Main result

Taking advantage of the previous lemmas, the following extension to Theorem 1 is stated:

###### Theorem 2

Consider system (1) with a given speed , a viscous damping and initial conditions . Assume that, for a given integer , there exist and such that inequality

 ΨN=He(Z⊤NPNFN)−c~RN+c(H⊤N(S+R)HN−G⊤NSGN)≺0 (24)

holds, where

 (25)

and where matrices , and are given in (23).

Then, the coupled infinite dimensional system (1) is exponentially stable in the sense of norm and there exist and such that energy estimate (10) holds.

###### Remark 7

Remark 5 also applies for this theorem and it means that must be strictly positive. In other words, these theorems cannot ensure the stability of the interconnection if the PDE is undamped.
Also note that Theorem 2 with leads exactly to the same conditions as presented in Theorem 1.

###### Remark 8

This methodology introduces a hierarchy in the stability conditions inspired from what one can read in [26] in the case of time-delay systems. More precisely, the sets

 CN={c>0 s.t. ∃PN∈Sn+2(N+1)+,S,R∈S2+,ΨN≺0}

representing the parameters for which the LMI of Theorem 2 is feasible for a given system (1) and for a given , satisfy the following inclusion: . In other words, if there exists a solution to Theorem 2 at an order , then there also exists a solution at any order . The proof is very similar to the one given in [26]. We can proceed by induction with and a sufficiently small . Then, . The calculations are tedious and technical and we do not intend to give them in this article.

### Iv-C Proof of Theorem 2

The proof of dissipativity follows the same line as in Theorem 1 and consists in proving the existence of positive scalars and such that the functional verifies the inequalities given in (11).

#### Iv-C1 Well-posedness

Using a similar reasoning to Theorem 1, a necessary condition for LMI (24) to be verified is that is non singular. Then, according to Propositions 1 and 2, the problem is well-posed and is the unique equilibrium.

#### Iv-C2 Existence of ε1

It strictly follows the proof in Theorem 1 and is therefore omitted.

#### Iv-C3 Existence of ε2

Since and are definite positive matrices, there exists such that:

 PN⪯diag(ε2In,ε24diag{(2k+1)In}k∈(0,N)),(S+xR)⪯S+R ⪯ ε24I2,∀x∈(0,1).

Then, from equation (21), we get:

 VN(XN,u)⩽ε2|X|2nN∑k=0+ε24(N∑k=0(2k+1)X⊤kXk+∥χ∥2)⩽ε2(|X|2n+12∥χ∥2).

While the first inequality is guaranteed by the constraint , for all the second estimate results from the application of Bessel inequality (22). Therefore, following the same procedure as in the proof of Theorem 1 after equation (12), there indeed exists such that .

#### Iv-C4 Existence of ε3

Differentiating in time defined in (21) along the trajectories of system (1) leads to:

The goal here is to find an upper bound of using the following extended state: . Using equation (15) and Lemma 3, we note that where matrices are given in (