David Bou

Federico Bribiesca-Argomedo

Jean Auriol

Michael Di Loreto

Florent Di Meglio

Stability Analysis for a Class of Linear 2 × 2 Hyperbolic PDEs Using a Backstepping Transform

In this paper, we develop a sufficient stability condition for a class of coupled first-order linear hyperbolic PDEs with constant coefficients that appear when considering target systems for backstepping boundary control. Using a backstepping transform, the problem is reformulated as a stability problem for a difference equation with distributed delay. Finding the explicit solution to the backstepping kernels, we derive an explicit sufficient condition depending on the plant coefficients. This stability condition is compared to an existing stability result based on a Lyapunov analysis. Both the proposed and existing sufficient conditions are then contrasted in some examples to a (computationally expensive) numerical approximation of a necessary and sufficient condition for exponential stability to illustrate their conservatism.

Index Terms-Hyperbolic Partial Differential Equations, Difference equations, Stability analysis, Backstepping design.

I. INTRODUCTION

Systems of first-order hyperbolic PDEs, often representing systems of conservation and balance laws, have been broadly studied both from the control and analysis perspectives. They appear in many practically relevant engineering models, including those of electric transmission lines [START_REF] Magnusson | Transmission lines and wave propagation[END_REF], hydraulic channels [START_REF] De Halleux | Boundary feedback control in networks of open channels[END_REF], pipelines [START_REF] Rager | Simplified fluid transmission line model for pneumatic control applications[END_REF]- [START_REF] Gugat | Time-delayed boundary feedback stabilization of the isothermal Euler equations with friction[END_REF], heat exchangers [START_REF] Xu | Exponential stability and transfer functions of processes governed by symmetric hyperbolic systems[END_REF] and oil well drilling [START_REF] Aamo | Leak detection, size estimation and localization in pipe flows[END_REF]. Different control approaches have been employed for the stabilization or observation of such systems (see for instance [START_REF] Aamo | Leak detection, size estimation and localization in pipe flows[END_REF], [START_REF] Coron | A strict Lyapunov function for boundary control of hyperbolic systems of conservation laws[END_REF], [START_REF] Santos | Boundary control with integral action for hyperbolic systems of conservation laws: Stability and experiments[END_REF], [START_REF] Vazquez | Backstepping boundary stabilization and state estimation of a 2×2 linear hyperbolic system[END_REF], [START_REF] Coron | Local exponential H 2 stabilization of a 2×2 quasilinear hyperbolic system using backstepping[END_REF] and [START_REF] Lamare | Control of 2×2 linear hyperbolic systems: Backstepping-based trajectory generation and PI-based tracking[END_REF]). These control designs are accompanied by a stability analysis to ensure that the resulting closed-loop system (or, equivalently, the target system in the case of backstepping control) is stable. When in-domain and boundary couplings are both present (both of which can be sources of instability), the stability analysis of these systems can be challenging. Therefore, practical ways to guarantee their stability are required. Some existing stability conditions for systems of conservation and balance laws (in different norms) can be found for instance in [START_REF] Bastin | Stability and boundary stabilization of 1-D hyperbolic systems[END_REF]. Other available results, based on (weighted L 2 ) Lyapunov functions, are [START_REF] Coron | Dissipative boundary conditions for one-dimensional nonlinear hyperbolic systems[END_REF] and [START_REF] Prieur | Lyapunov functions for switched linear hyperbolic systems[END_REF] for the nonlinear and switched hyperbolic cases, respectively. The authors of [START_REF] Tchousso | Exponential stability of distributed parameter systems governed by symmetric hyperbolic partial differential equations using Lyapunov's second method[END_REF] treat the case of symmetric hyperbolic PDEs, and prove the stability of these systems under the hypothesis of dissipation both inside the domain and through the boundary conditions. In [START_REF] Safi | Tractable sufficient stability conditions for a system coupling linear transport and differential equations[END_REF], a sufficient stability condition was developped for Corresponding author: M. Bou Saba (email: david.bou-saba@insa-lyon.fr).

a coupled hyperbolic PDE-ODE system using a Lyapunov functional. Singular perturbation analysis was used in [START_REF] Tang | Tikhonov theorem for linear hyperbolic systems[END_REF] for systems of conservation laws and in [START_REF]Singular perturbation approximation of linear hyperbolic systems of balance laws[END_REF] for hyperbolic systems of balance laws. This method was explored in [START_REF] Tang | Stability analysis of coupled linear ODEhyperbolic PDE systems with two time scales[END_REF] and [START_REF] Tang | Stability analysis of a singularly perturbed coupled ODE-PDE system[END_REF] for the development of sufficient stability conditions for systems of hyperbolic PDEs coupled to an ODE, and in [START_REF] Cerpa | Effect of time scales on stability of coupled systems involving the wave equation[END_REF] for the case of a wave equation coupled to an ODE. Also, in [START_REF] Karafyllis | Input-to-State Stability for PDEs[END_REF], Input-to-State Stability (ISS) bounds for hyperbolic PDEs are used to give sufficient conditions for stability of systems of two coupled hyperbolic PDEs in the sup norm.

Interestingly, systems of first-order hyperbolic PDEs can be transformed to neutral functional equations as in [START_REF] Russell | Neutral FDE canonical representations of hyperbolic systems[END_REF] or to a difference-delay system as in [START_REF]Canonical forms and spectral determination for a class of hyperbolic distributed parameter control systems[END_REF]. In [START_REF] Karafyllis | On the relation of delay equations to first-order hyperbolic partial differential equations[END_REF], the authors established the equivalence between a single first-order hyperbolic PDE and systems described by an integral delay equation by proving that the mapping between the solutions is unique in appropriate spaces. Furthermore, it is easier to analyze the stability when converting the first-order hyperbolic PDEs to a delay form. In [START_REF] Auriol | Delay-robust stabilization of a hyperbolic PDE-ODE system[END_REF], a difference equation with distributed delay was derived from 2 coupled hyperbolic PDEs to achieve a robust stabilization with respect to input delays.

Given this relation between hyperbolic PDEs and delay systems it makes sense to also consider existing stability results for the latter. In the case of difference equations with pointwise and distributed delays, stability conditions can be found in [START_REF] Damak | Stability of linear continuoustime difference equations with distributed delay: Constructive exponential estimates[END_REF] and [START_REF] Melchor-Aguilar | Exponential stability of linear continuous time difference systems with multiple delays[END_REF], and, in the case of delay-differential equations in [START_REF] Hale | Strong stabilization of neutral functional differential equations[END_REF], [START_REF] Glizer | A novel approach to exact slowfast decomposition of linear singularly perturbed systems with small delays[END_REF], [START_REF] Bellman | Differential-difference equations[END_REF], [START_REF] Diekmann | Delay equations: functional-, complex-, and nonlinear analysis[END_REF] and [START_REF] Niculescu | On delay-dependent stability under model transformations of some neutral linear systems[END_REF].

The motivation for this work is to enlarge the class of target systems that can be considered when designing a backstepping controller. In the first available backstepping designs for these systems, boundary reflections in the target system were regarded as undesirable and were eliminated, together with some well-chosen in-domain coupling terms, allowing for a relatively simple stability analysis [START_REF] Vazquez | Backstepping boundary stabilization and state estimation of a 2×2 linear hyperbolic system[END_REF], [START_REF] Di Meglio | Stabilization of a linear hyperbolic system with one boundary controlled transport PDE coupled with n counterconvecting PDEs[END_REF], [START_REF] Hu | Control of homodirectional and general heterodirectional linear coupled hyperbolic PDEs[END_REF], [START_REF] Di Meglio | Stabilization of coupled linear heterodirectional hyperbolic PDE-ODE systems[END_REF]. Given that all of these systems were finite-time stable (and many were minimal-time convergent [START_REF] Auriol | Minimum time control of heterodirectional linear coupled hyperbolic PDEs[END_REF], [START_REF]Two-sided boundary stabilization of two linear hyperbolic PDEs in minimum time[END_REF]) the question of adding degrees of freedom to the controller was not extremely urgent to address. Nevertheless, recent results have made clear that removing reflections often results in closed-loop systems that have zero delay margins [START_REF] Auriol | Delay-robust control design for heterodirectional linear coupled hyperbolic PDEs[END_REF], [START_REF] Auriol | Delay-robust stabilization of a hyperbolic PDE-ODE system[END_REF]. Since it is now known that, in general, some reflections must be preserved, we end up with few target system options for which the stability analysis is simple. At the same time, any new target systems will naturally be exponentially stable instead of finite-time stable, which means that adding degrees of freedom that allow us to alter the rate of convergence and transient responses is more important than ever. It is within this context that we develop a new sufficient condition for the stability analysis of such target systems that is easy to verify and depends explicitly on the system coefficients.

It should be noted that this is not the first time alternate target systems have been studied. As an example, some recent works on boundary control of port-Hamiltonian systems, have also considered alternative internal coupling structures [START_REF] Ramirez | On backstepping boundary control for a class of linear port-hamiltonian systems[END_REF], [START_REF] Macchelli | Boundary control of distributed port-hamiltonian systems via generalised canonical transformations[END_REF], yet these designs rely on canceling boundary reflections and the use of specific strict Lyapunov functions for stability analysis, which impose a particular structure for the retained closed-loop system (it requires some degree of dissipation in the domain).

In this work, we analyze the stability of a class of 2 × 2 linear hyperbolic coupled PDEs. Using a backstepping transformation, the original system is mapped to a target system that can be written as a difference equation with distributed delay. This fact highlights the usefulness of the backstepping technique not only for control design, but for analysis as well. Furthermore, to the best of the authors' knowledge, few easily verifiable conditions exist for the class of systems obtained after this backstepping transform, which contain a distributed delay, and even then they would not be expressed in terms of the original system parameters, adding a supplementary layer of complexity to the analysis.

The contribution of this paper is the introduction of an easy to apply, explicit in terms of the plant parameters, sufficient condition for stability for a class of hyperbolic PDEs that can be naturally obtained through the backstepping procedure when applied to sytems of 2 × 2, counterconvecting, linear hyperbolic PDEs with constant coefficients.

In order to construct this result, we need to derive the explicit expressions of the backstepping kernels, using a similar procedure to that used in [START_REF] Vazquez | Marcum Q-functions and explicit kernels for stabilization of 2×2 linear hyperbolic systems with constant coefficients[END_REF]. Incidentally, combining the results of [START_REF] Vazquez | Marcum Q-functions and explicit kernels for stabilization of 2×2 linear hyperbolic systems with constant coefficients[END_REF] and the present paper yields analytic expressions of the complete set of backstepping kernels (both direct and inverse) which can be useful for analyzing transient behavior of the closed-loop system. These analytic expressions allow us to express the target system as a difference equation with distributed delay for which a necessary and sufficient stability condition can be found in the Laplace domain [START_REF] Vidyasagar | Input-output stability of a broad class of linear timeinvariant multivariable systems[END_REF]. Upperbounding some of the terms in this necessary and sufficient condition we are able to derive a relatively simple sufficient stability condition, which is roughly equivalent (in terms of complexity) to that in [START_REF] Bastin | Stability and boundary stabilization of 1-D hyperbolic systems[END_REF]Theorem 5.4] for the same class of PDEs. These two explicit sufficient conditions obtained by completely different approaches are then compared and shown to not be equivalent, as neither one implies the other.

In order to illustrate the conservatism of the proposed condition, it is compared to a numerical approximation of the necessary and sufficient condition. This approximation is obtained by leveraging the explicit expressions of the backstepping kernels in order to both restrict the region of the complex plane where unstable poles may exist and then to determine a stopping criterion that, for stable systems and with a perfect computer, should result in a finite execution time. In practice, some tolerance is added in order to account for numerical errors and to avoid infinite execution times for unstable systems.

The paper is organized as follows. In Section II, we present the class of systems under consideration. In Section III, we reformulate the original system as a difference equation with distributed delay by means of a backstepping transform and we give the associated necessary and sufficient stability condition in the Laplace domain. In Section IV, an explicit sufficient stability condition in terms of the original system coefficients is proposed. This condition is compared to an existing (and similarly complex) condition obtained by Lyapunov analysis directly and through some examples. A numerical approximation of the necessary and sufficient condition is then constructed. In Section V, this numerical approximation is applied, together with both explicit conditions, to some illustrative examples to assess the conservatism of the proposed condition. Concluding remarks are given in Section VI.

Throughout this paper, we use the following definitions. We define the functional space χ . = L 2 ([0, 1]; R) × L 2 ([0, 1]; R), the space where the PDE states will be defined, with the associated norm

(1) (ϕ, ψ) χ . = 1 0 ϕ 2 (θ)dθ + 1 0 ψ 2 (θ)dθ 1/2
, for any (ϕ, ψ) ∈ χ . We define as well the partial trajectory of a function ϕ as

ϕ [t] (θ) . = ϕ(t + θ), θ ∈ [-τ, 0), with norm in L 2 ([-τ, 0); R), τ > 0 (2) ϕ [t] L 2 ([-τ,0);R) . = 0 -τ ϕ 2 (t + θ)dθ 1/2
for t ≥ τ , and its norm in C 1 ([-τ, 0]; R) will be defined as

(3) ϕ [t] C 1 ([-τ,0];R) . = max θ∈[-τ,0] |ϕ(t + θ)|+ max θ∈[-τ,0] |ϕ t (t + θ)|
for t ≥ τ . ∂ t ϕ will designate the partial derivative of ϕ with respect to the independent variable t.

We define the Banach algebra A equipped with a standard convolution product in the sense of [START_REF] Vidyasagar | Input-output stability of a broad class of linear timeinvariant multivariable systems[END_REF]. A consists of BIBOstable generalized functions of the form (4)

g(t) = g 0 (t) + ∞ i=1 g i δ(t -t i ), with g 0 ∈ L 1 (R + ; R) and ∞ i=1 |g i |< ∞. δ denotes the unit Dirac distribution. The algebra A is equipped with the norm (5) g A . = g 0 L 1 + ∞ i=1 |g i |.
We denote  the set of Laplace transforms of elements in A.  is a Banach algebra with unit. The norm on  satisfies ĝ Â= g A , where ĝ stands for the Laplace transform of g ∈ A.

J n and I n , for n ∈ Z, stand for the Bessel functions and the modified Bessel functions of the first kind, respectively.

Given a set Ω ⊆ R, its characteristic function 1 Ω : R → R is defined as

1 Ω (θ) . = 1 if θ ∈ Ω 0 otherwise. A function f defined on ω ⊆ C is Lipschitz-continuous with Lipschitz constant m on ω, if there exists a constant m ≥ 0 such that |f (z 1 ) -f (z 2 )|≤ m|z 1 -z 2 |, ∀z 1 , z 2 ∈ ω,

II. PROBLEM STATEMENT

We are interested in the stability analysis of a class of 2 × 2 linear first-order hyperbolic coupled PDEs. Systems under consideration are of the form

(6) ∂ t u(t, x) + 1 r 1 ∂ x u(t, x) = σ 1 v(t, x) (7) ∂ t v(t, x) - 1 r 2 ∂ x v(t, x) = σ 2 u(t, x)
for all t ≥ 0, x ∈ [0, 1], where 1 r1 > 0 and 1 r2 > 0 are the transport velocities associated to the two transport equations ( 6) and ( 7) respectively. The coefficients σ 1 and σ 2 represent constant in-domain couplings between the PDEs. The associated initial conditions are [START_REF] Santos | Boundary control with integral action for hyperbolic systems of conservation laws: Stability and experiments[END_REF]

u 0 (•) . = u(0, •) v 0 (•) . = v(0, •), defined in L 2 ([0, 1], R).
In addition, u(t, x) and v(t, x) satisfy the following boundary conditions ( 9)

u(t, 0) = qv(t, 0) v(t, 1) = ρu(t, 1)
for all t > 0, where ρ and q are boundary coupling coefficients such that [START_REF] Coron | Local exponential H 2 stabilization of a 2×2 quasilinear hyperbolic system using backstepping[END_REF] |ρq|< 1.

Note that if |ρq|≥ 1, the system cannot be stable [START_REF] Bastin | Stability and boundary stabilization of 1-D hyperbolic systems[END_REF]Theorem 2.6]. We assume that ρ and q are not both zero. Without loss of generality, we consider that q = 0. The reader will notice first that no diagonal coupling terms are considered in ( 6)- [START_REF] Coron | A strict Lyapunov function for boundary control of hyperbolic systems of conservation laws[END_REF]. This is partly because these terms can be removed using an exponential transform together with the backstepping control design presented in [START_REF] Vazquez | Backstepping boundary stabilization and state estimation of a 2×2 linear hyperbolic system[END_REF], [START_REF] Coron | Local exponential H 2 stabilization of a 2×2 quasilinear hyperbolic system using backstepping[END_REF], and therefore such a target system is always achievable (when considering constant-coefficient equations), and also because, if kept, they prevent the development of the explicit expressions for the backstepping kernels. One could also attempt to use only the exponential change of variables such as that in [START_REF] Lamare | Control of 2×2 linear hyperbolic systems: Backstepping-based trajectory generation and PI-based tracking[END_REF], [START_REF] Bastin | On boundary feedback stabilization of non-uniform linear 2×2 hyperbolic systems over a bounded interval[END_REF] in order to absorb the diagonal coupling terms into the off-diagonal ones at the expense of non-constant coefficients, which pose the same problem. Therefore we only consider the case where diagonal coefficients either do not appear in the original system or, more likely, are compensated by a backstepping-based control.

Our goal is therefore to construct sufficient stability conditions that guarantee the exponential stability of system ( 6)-( 9) in an L 2 sense. That is, for any (u 0 , v 0 ) ∈ χ, there exist µ > 0 and C ≥ 0, such that the solution (u, v) of system ( 6)-( 9) satisfies [START_REF] Lamare | Control of 2×2 linear hyperbolic systems: Backstepping-based trajectory generation and PI-based tracking[END_REF] (u(t,

•), v(t, •)) χ ≤ C (u 0 , v 0 ) χ e -µt , t ≥ 0.
As a first step, we will employ a backstepping transform that leads to a target system where the in-domain couplings are shifted to the boundary. The structure of the target system will then allow us to reformulate the system as a difference equation with distributed delay.

III. DIFFERENCE EQUATION WITH DISTRIBUTED DELAY

In this section, we consider the class of PDEs described by ( 6)-( 9) and establish an equivalent realization as a difference equation. The objective of this section is to give a necessary and sufficient stability condition of the system. To do so, we use a backstepping transform, then the Laplace transform and therefore analyze the BIBO stability [START_REF] Vidyasagar | Input-output stability of a broad class of linear timeinvariant multivariable systems[END_REF] in Laplace domain.

A. Backstepping transform

In order to rewrite the original system in a difference equation form, we employ the following Volterra transformation used in [START_REF] Coron | Local exponential H 2 stabilization of a 2×2 quasilinear hyperbolic system using backstepping[END_REF] (12)

u(t, x) = α(t, x) - x 0 L 11 (x, y)α(t, y)dy - x 0 L 12 (x, y)β(t, y)dy v(t, x) = β(t, x) - x 0 L 21 (x, y)α(t, y)dy - x 0 L 22 (x, y)β(t, y)dy,
where the kernels L ij , i, j ∈ {1, 2} verify the following PDEs in the triangular domain T

. = {(x, y) ∈ [0, 1] × [0, 1], y ≤ x} [10] (13) 1 r 1 ∂ x L 11 (x, y) + 1 r 1 ∂ y L 11 (x, y) = σ 1 L 21 (x, y) 1 r 1 ∂ x L 12 (x, y) - 1 r 2 ∂ y L 12 (x, y) = σ 1 L 22 (x, y) 1 r 2 ∂ x L 21 (x, y) - 1 r 1 ∂ y L 21 (x, y) = -σ 2 L 11 (x, y) 1 r 2 ∂ x L 22 (x, y) + 1 r 2 ∂ y L 22 (x, y) = -σ 2 L 12 (x, y),
with boundary conditions ( 14)

L 12 (x, x) = -σ 1 r 1 r 2 r L 11 (x, 0) = L 12 (x, 0) r 1 qr 2 L 21 (x, x) = σ 2 r 1 r 2 r L 22 (x, 0) = L 21 (x, 0) r 2 r 1 q,
where we denote r . = r 1 + r 2 . Given q = 0, equations ( 13)-( 14) are of the class considered in [10, Theorem A.2]. Accordingly, these equations have a unique solution

L ij , i, j ∈ {1, 2} in C ∞ (T ; R).
To the best of the authors' knowledge, an analytic solution for these equations is not given in the literature, including [START_REF] Coron | Local exponential H 2 stabilization of a 2×2 quasilinear hyperbolic system using backstepping[END_REF]. Since an analytic solution is needed in our analysis, we develop in Appendix A the explicit expression of the backstepping kernels in terms of Bessel functions.

Lemma 1: PDEs (13) with boundary conditions [START_REF] Prieur | Lyapunov functions for switched linear hyperbolic systems[END_REF], have a unique smooth solution, given by ( 63)-(66). The proof is given in Appendix A.

We consider the following target system (15)

∂ t α(t, x) + 1 r 1 ∂ x α(t, x) = 0 ∂ t β(t, x) - 1 r 2 ∂ x β(t, x) = 0,
with initial conditions ( 16)

α 0 (•) . = α(0, •) β 0 (•) . = β(0, •), defined in L 2 ([0, 1], R), and boundary conditions (17) α(t, 0) = qβ(t, 0) (18)
β(t, 1) = ρα(t, 1) + 1 0 N α (y)α(t, y)dy + 1 0 N β (y)β(t, y)dy, where (19)
N α (y) . = L 21 (1, y) -ρL 11 (1, y) N β (y) . = L 22 (1, y) -ρL 12 (1, y).
The backstepping transform ( 12) is a bounded linear operator from χ to χ [START_REF] Coron | Local exponential H 2 stabilization of a 2×2 quasilinear hyperbolic system using backstepping[END_REF], and maps the target system ( 15)-( 18) to the original system (6)- [START_REF] Vazquez | Backstepping boundary stabilization and state estimation of a 2×2 linear hyperbolic system[END_REF]. Notice that the backstepping transform ( 12) is an inverse transform, i.e., it explicitly maps the target system to the original system. It is used instead of the direct transform since it allows us to write the boundary condition [START_REF]Singular perturbation approximation of linear hyperbolic systems of balance laws[END_REF] as function of the states α and β.

Proposition 2: For any initial condition (α 0 , β 0 ) ∈ χ, the target system ( 15)-( 18) has a unique weak solution

(α(t, •), β(t, •)) ∈ χ, for t ≥ 0 and β [t] (•, 1) ∈ L 2 ([-r, 0); R), ∀t ≥ r.
Proof: An initial condition (α 0 , β 0 ) ∈ χ is mapped, by the backstepping transform [START_REF] Bastin | Stability and boundary stabilization of 1-D hyperbolic systems[END_REF], to an initial condition (u 0 , v 0 ) ∈ χ . According to [START_REF] Bastin | Stability and boundary stabilization of 1-D hyperbolic systems[END_REF]Theorem A.4], the Cauchy problem ( 6)-( 9) with this initial condition has a unique solution (u(t, •), v(t, •)) ∈ χ . Since the transform ( 12) is a boundedly invertible operator from χ to χ, we conclude the proof.

Remark 1: Assume that u 0 and v 0 ∈ C 1 ([0, 1]; R), and the initial conditions of the original system verify the following compatibility conditions [START_REF] Tang | Stability analysis of a singularly perturbed coupled ODE-PDE system[END_REF]

u 0 (0) = qv 0 (0) v 0 (1) = ρu 0 (1). u 0 (0) = -q r 1 r 2 v 0 (0) -qr 1 σ 2 u 0 (0) + r 1 σ 1 v 0 (0) v 0 (1) = -ρ r 2 r 1 u 0 (1) -r 2 σ 2 u 0 (1) + ρr 2 σ 1 v 0 (1).
Adapting the proof of [START_REF] Coron | Local exponential H 2 stabilization of a 2×2 quasilinear hyperbolic system using backstepping[END_REF]Theorem A.1] to C 1 ([0, 1]; R), we have that the target system (15) has a unique solution α(t, x) and β(t, x) verifying α(t, •) and

β(t, •) ∈ C 1 ([0, 1]; R), ∀t ≥ 0 and β [t] (•, 1) ∈ C 1 ([-r, 0]; R), ∀t ≥ r.

B. Distributed delay form

Following the method of characteristics, the transport equations [START_REF] Tchousso | Exponential stability of distributed parameter systems governed by symmetric hyperbolic partial differential equations using Lyapunov's second method[END_REF] can be written as the following difference equations [START_REF] Cerpa | Effect of time scales on stability of coupled systems involving the wave equation[END_REF] 21)-( 22) in the terms of the boundary condition [START_REF]Singular perturbation approximation of linear hyperbolic systems of balance laws[END_REF] verified by β(t, 1) and using ( 17) yield ( 23)

α(t, x) = α(t -r 1 x, 0), ∀t ≥ r 1 x (22) β(t, x) = β(t -r 2 (1 -x), 1), ∀t ≥ r 2 (1 -x) ∀x ∈ [0, 1]. Plugging equations (
α(t, 1) = α(t -r 1 , 0) = qβ(t -r, 1), 1 0 N α (y)α(t, y)dy = q 1 0 N α (y)β(t -r 2 -r 1 y, 1)dy = q r 1 r r2 N α θ -r 2 r 1 β(t -θ, 1)dθ, 1 0 N β (y)β(t, y)dy = 1 0 N β (y)β(t -r 2 (1 -y), 1)dy = 1 r 2 r2 0 N β 1 - θ r 2 β(t -θ, 1)dθ.
Finally, the difference equation with distributed delay satisfied by β(t, 1) is ( 24) 25)). N (θ) can be expressed in terms of the Bessel functions (see Appendix A, equation ( 69)). The continuity of N at r 2 is ensured due to the same expression satisfied by N (θ) for θ ∈ [0, r 2 ] (equation (67)) and θ ∈ (r 2 , r] (equation (68)).

β(t, 1) = ρqβ(t -r, 1) + r 0 N (θ)β(t -θ, 1)dθ for t ≥ r, with (25)
N (θ) = 1 r 2 N β 1 - θ r 2 1 [0,r2] (θ) + q r 1 N α θ -r 2 r 1 1 (r2,r] (θ), θ ∈ [0, r]. Notice that N α and N β ∈ C ∞ ([0, 1]; R) (as defined in (19)) and N ∈ C ∞ ([0, r]; R) ⊂ L 1 ([0, r]; R) (as defined in (
Equations ( 21)- [START_REF] Karafyllis | Input-to-State Stability for PDEs[END_REF] represent the solution to the PDE states inside the domain as a function of the solution at the boundaries.

One can check that the stability of β(•, 1) implies in particular the stability of the state β satisfying [START_REF] Karafyllis | Input-to-State Stability for PDEs[END_REF]. In turn, the stability of the state α is ensured at x = 0 by [START_REF] Tang | Tikhonov theorem for linear hyperbolic systems[END_REF] and inside the domain by [START_REF] Cerpa | Effect of time scales on stability of coupled systems involving the wave equation[END_REF]. Furthermore, due to the boundedness and the invertibility of the backstepping transfrom, there is an equivalence between the stability of the original and the target systems.

We also note that for t ∈ [0, r], the regularity of the solution of the target system (15)-( 18) depends on the regularity of the initial conditions α 0 and β 0 . Thus, since the transform ( 12) is boundedly invertible from χ to χ [START_REF] Coron | Local exponential H 2 stabilization of a 2×2 quasilinear hyperbolic system using backstepping[END_REF], any initial condition (u 0 , v 0 ) ∈ χ is mapped to an initial condition (α 0 , β 0 ) ∈ χ. Therefore, (α(t, •),

β(t, •)) ∈ χ for t ∈ [0, r].
Consequently, we will focus our analysis on the boundary stability criterion and give stability conditions for β(•, 1) satisfying the difference equation with distributed delay [START_REF]Canonical forms and spectral determination for a class of hyperbolic distributed parameter control systems[END_REF].

C. Necessary and sufficient stability condition

We now present a necessary and sufficient stability condition for our system rewritten as a difference equation with distributed delay. Using the Laplace transform, the characteristic equation associated to ( 24) is

(26) φ(s) . = 1 -ρqe -sr - r 0 N (θ)e -sθ dθ = 0.
For the stability analysis, we restrict our study to the halfplane Re(s) ≥ 0. Consider the transfer function

(27) Ĝ(s) . = - r 0 N (θ)e -sθ dθ 1 -ρqe -sr •
With |ρq|< 1, one can verify that the conditions required by [43, Theorem 1] are fulfilled by Ĝ allowing us to apply the stability result therein. In addition, we have that

(28) F (s) . = Ĝ(s) 1 + Ĝ(s) = - r 0 N (θ)e -sθ dθ φ(s)
and ( 29)

1 + Ĝ(s) = φ(s) |1 -ρqe -sr | • Since N (•) ∈ C ∞ ([0, r]; R), r 0 N (θ)e -sθ
dθ is a Finite Impulse Response (FIR) filter, the stability of F (s) is equivalent to the absence of roots of φ(s) in the complex right half-plane. According to [START_REF] Vidyasagar | Input-output stability of a broad class of linear timeinvariant multivariable systems[END_REF]Theorem 1], φ has an inverse in  if and only if [START_REF] Glizer | A novel approach to exact slowfast decomposition of linear singularly perturbed systems with small delays[END_REF] inf

Re(s) ≥0 |1 + Ĝ(s)|> 0,
which, by [START_REF] Hale | Strong stabilization of neutral functional differential equations[END_REF], is equivalent to

(31) inf Re(s) ≥0 | φ(s)|> 0.
Condition ( 31) is necessary and sufficient for the exponential stability of β(•, 1) satisfying (24) in L 2 ((-r, 0]; R) according to [START_REF] Baker | Input-output stability of linear time-invariant systems[END_REF]Remark 1], provided that (α 0 , β 0 ) ∈ χ. The stability of system ( 6)-( 9) in the χ-norm is consequently equivalent to [START_REF] Bellman | Differential-difference equations[END_REF].

This necessary and sufficient condition cannot be easily verified since it requires the evaluation of φ at every point s in the complex right half-plane (RHP). This evaluation is complicated by the distributed delay term (integral term in [START_REF] Auriol | Delay-robust stabilization of a hyperbolic PDE-ODE system[END_REF]) and, therefore, it is not a practical tool to check the stability of the system under consideration.

However, one may proceed by upper-bounding the modulus of the terms in φ. This yields a simpler (only sufficient) condition, that is much easier to apply and depends only on the system coefficients. This is the purpose of Section IV-A.

Another approach to tackle this problem would be to numerically approximate this condition. This is done in Section IV-C for validation purposes. Despite being more computationally expensive, the numerical approximation is less conservative than the explicit criterion. Crucially for the implementation of the algorithm, we only need to verify condition [START_REF] Bellman | Differential-difference equations[END_REF] in a compact domain in the RHP which can be analytically bounded in terms of the system coefficients, as will appear in Section IV-C.

Both alternatives require some sort of solution of N (θ). It is possible to get the backstepping gains without an explicit expression, by using numerical successive approximations. Nevertheless, applying this method is more complicated than using their analytic expressions and could yield more conservative criteria since approximation errors have to be taken into account. Instead, we explicitly solve the hyperbolic PDEs ( 13)-( 14) and obtain the resulting expression of N (θ). This is essential in order to state an explicit stability criterion and to construct the compact domain where the numerical algorithm makes the tests. In addition, it is used to derive the expressions of many parameters required by the algorithm (as will be detailed in section IV-C). Appendix A contains the derivations of the explicit solution to the hyperbolic PDEs satisfied by the backstepping kernel gains, leading to the expressions of N (θ), N (θ) and the upper-bounds of their absolute values.

IV. STABILITY ANALYSIS OF A DIFFERENCE EQUATION WITH DISTRIBUTED DELAY: SUFFICIENT CONDITIONS

In this section, we present the explicit sufficient stability condition for the system under consideration, obtained by relaxing the necessary and sufficient condition previosly presented. This explicit criterion, in terms of the system coefficients is compared to a similarly complex criterion available in the literature and, by numerically approximating the necessary and sufficient condition, an algorithm is constructed. Its task is to verify the validity of (31) inside a compact domain in the RHP.

A. Explicit sufficient condition

We state now a sufficient stability condition of the system under consideration, depending on the sign of σ 1 σ 2 and ρq. In order to simplify the notation, we define the parameters [START_REF] Diekmann | Delay equations: functional-, complex-, and nonlinear analysis[END_REF] a .

= qr 2 σ 2 + ρr 1 σ 1 (33) R . = r 1 r 2 σ 1 σ 2 .
Proposition 3: If the constant parameters of system ( 6)-( 9) satisfy either of the following set of inequalities:

(i) σ 1 σ 2 ≥ 0, ρq ≥ 0 and

|a|+|R| 1 1 + |ρq| - 1 -|ρq| 2 < 1 -|ρq|, (ii) σ 1 σ 2 ≥ 0, ρq < 0 and |a|+|R| 1 + |ρq| 2 < 1 -|ρq|, (iii)
σ 1 σ 2 < 0, ρq ≥ 0 and

|a|I 0 |R| + |R| 1 1 + |ρq| - 1 -|ρq| 2 × I 0 |R| -I 2 |R| < 1 -|ρq|, (iv)
σ 1 σ 2 < 0, ρq < 0 and

|a|I 0 |R| + |R| 1 + |ρq| 2 × I 0 |R| -I 2 |R| < 1 -|ρq|,
then, system ( 6)-( 9) is exponentially stable in the χ-norm in the sense of [START_REF] Lamare | Control of 2×2 linear hyperbolic systems: Backstepping-based trajectory generation and PI-based tracking[END_REF].

Proof: In this proof, we rely on the necessary and sufficient condition [START_REF] Bellman | Differential-difference equations[END_REF], which is fulfilled if [START_REF] Di Meglio | Stabilization of a linear hyperbolic system with one boundary controlled transport PDE coupled with n counterconvecting PDEs[END_REF] sup

Re(s) ≥0 |1 -φ(s)|< 1.
|1 -φ(s)| can be bounded as

(35) |1-φ(s)|= ρqe -sr + r 0 N (θ)e -sθ dθ ≤ |ρq|+ r 0 |N (θ)|dθ.
Using [START_REF] Hu | Control of homodirectional and general heterodirectional linear coupled hyperbolic PDEs[END_REF], a sufficient condition for ( 34) is

(36) r 0 |N (θ)|dθ < 1 -|ρq|.
The crucial point is then the upper-bound of r 0 |N (θ)|dθ given in (84) in Appendix A, depending on the sign of σ 1 σ 2 and ρq. Substituting these expressions into (36), we conclude the proof.

Remark that inequality (36) can be tested by numerically approximating the integral term and using the expression of N (θ) in (69). In the case of spatially varying coefficients σ 1 (x) and σ 2 (x), [START_REF] Di Meglio | Stabilization of coupled linear heterodirectional hyperbolic PDE-ODE systems[END_REF] can also be tested, with N (θ) generated numerically using the successive approximation method.

Following the assumptions of Remark 1, and according to [45, Remark 1], the sufficient condition of Proposition 3 implies the point-wise exponential stability of system [START_REF]Canonical forms and spectral determination for a class of hyperbolic distributed parameter control systems[END_REF]. In other words, there exist µ > 0 and C ≥ 0 such that [START_REF] Auriol | Minimum time control of heterodirectional linear coupled hyperbolic PDEs[END_REF]

|β(t, 1)|≤ C sup -r≤θ≤0 |β(θ, 1)|e -µ(t-r) , t ≥ r.
Let us now compare the presented result to some existing results in the literature, and apply it to some illustrative examples.

B. Comparison of the explicit sufficient condition with other criteria and numerical validation

In order to evaluate the conservatism of the sufficient condition in Proposition 3, we compare it with other sufficient conditions from the literature for the stability of systems of linear first-order hyperbolic PDEs. In particular, we compare it with a Lyapunov-based exponential stability condition for linear systems of balance laws presented in [START_REF] Bastin | Stability and boundary stabilization of 1-D hyperbolic systems[END_REF] and an ISS and small-gain based result in [START_REF] Karafyllis | Input-to-State Stability for PDEs[END_REF].

ρq ≥ 0 |a|+|R| 1 1+|ρq| - 1-|ρq| 2 < 1 -|ρq| ρq < 0 |a|+|R| 1+|ρq| 2 < 1 -|ρq| (Not satisfied) σ 1 = 0, σ 2 = 0 |ρq|< 1 |ρq|< 1 σ 1 σ 2 < 0 ρq ≥ 0 |a|I 0 |R| + |R| 1 1+|ρq| - 1-|ρq| 2 × I 0 |R| -I 2 |R| < 1 -|ρq| ρq < 0 |a|I 0 |R| + |R| 1+|ρq| 2 × I 0 |R| -I 2 |R| < 1 -|ρq| ρ 2 < -σ 2 σ 1 r 2 r 1 < 1 q 2
According to [12, Theorem 5.4], system (6)-( 9) (with constant coefficients σ 1 and σ 2 ) is exponentially stable in the L 2 -norm if there exists a 2 × 2 positive diagonal real

matrix Π . = π 1 0 0 π 2 such that (38)
Σ Π + ΠΣ 0 P 2 < 1 where Σ = 0 -σ 1 -σ 2 0 and P =   0 q π1 π2 r2 r1 ρ π2 π1 r1 r2 0   .
The norm considered in [START_REF]Two-sided boundary stabilization of two linear hyperbolic PDEs in minimum time[END_REF] is the spectral-norm, defined as the square root of the largest eigenvalue of P P . This condition was originally derived for systems with diagonal coupling terms, possibly spatially varying. In case of ( 6)-( 9), condition [START_REF]Two-sided boundary stabilization of two linear hyperbolic PDEs in minimum time[END_REF] reduces to: either

σ 1 = σ 2 = 0 , or (39)
σ 1 σ 2 < 0, and

ρ 2 < - σ 2 σ 1 r 2 r 1 < 1 q 2 •
Otherwise, if σ 1 σ 2 > 0 or only one of them is equal to zero, no conclusion can be made regarding the stability. We present in Table I the sufficient conditions required by Proposition 3 and [12, Theorem 5.4] for the different signs of σ 1 σ 2 and ρq.

One may check that in the case of σ 1 = σ 2 = 0 , system (6)-( 9) is viewed as a system of two linear conservation laws with two reflection terms at the boundaries. In the case of σ 1 σ 2 < 0 , the condition of Proposition 3 and ( 39) are both sufficient and can be applied simultaneously. For the sake of clarity, we give some numerical examples in Table II.

Another criterion for the stability of symmetric hyperbolic PDEs is found in [START_REF] Tchousso | Exponential stability of distributed parameter systems governed by symmetric hyperbolic partial differential equations using Lyapunov's second method[END_REF]. When applying [15, Theorem 2.1] for the class of 2 × 2 constant-coefficients, the required condition is [START_REF] Ramirez | On backstepping boundary control for a class of linear port-hamiltonian systems[END_REF] σ 1 + σ 2 = 0, and

c 1 r 2 + ρ 2 ≤ r 2 r 1 ≤ 1 q 2 -c 2 r 2 q 2 ,
with c 1 ≥ 0 and c 2 ≥ 0 and c 1 and c 2 are not both zero. This criterion, developed under the hypothesis of a dissipation on both equation and the boundary is more restrictive than [12, Theorem 5.4] regarding the stability of the system with respect to the values of σ 1 and σ 2 .

Using recent ISS results and a small-gain property, another sufficient stability condition for a related class of systems of hyperbolic PDEs is given in Chapter 9 of [START_REF] Karafyllis | Input-to-State Stability for PDEs[END_REF]. When applied to system (6)-( 9), this condition requires the existence of a constant K > 0, with (|ρ|+|q|)e -K < 1, such that (41)

e 2K -e K K r 2 |σ 2 | + |ρ| × e 2K -e K K r 1 |σ 1 | + |q| < 1.

This condition is not satisfied for the examples considered in

Table II so an extra column was not added. It will however be compared to the numerical algorithm developed in the next section.

As we can see from these results, the stability of some systems can be tested using the condition required by Proposition 3. Nevertheless, many other cases are still not covered, i.e. we cannot conclude regarding the stability of many plants by applying both sufficient criteria. The conservatism in our approach is a result of bounding the modulus of the exponential terms in φ in [START_REF] Hu | Control of homodirectional and general heterodirectional linear coupled hyperbolic PDEs[END_REF]. An alternative is to consider the necessary and sufficient condition and numerically approximate it. That will be the goal of the next section, where an algorithm is constructed and reveals to be less conservative than both approaches previously compared at a much higher computational cost.

C. Numerical approximation of the Necessary and Sufficient condition

We develop a numerical test of the condition [START_REF] Bellman | Differential-difference equations[END_REF] for the stability of system [START_REF]Canonical forms and spectral determination for a class of hyperbolic distributed parameter control systems[END_REF], based on the explicit expressions of N (θ), N (θ) and the upper-bounds of their absolute values, given in Appendix A. Our objective is to test the absence of zeros (up to a certain precision) in the right half plane of the characteristic equation [START_REF] Auriol | Delay-robust stabilization of a hyperbolic PDE-ODE system[END_REF].

Notice that numerical algorithms for spectrum analysis for time-delay systems as DELAY_PSA developed in [START_REF] Gumussoy | A predictor-corrector type algorithm for the pseudospectral abscissa computation of time-delay systems[END_REF], YALTA in [START_REF] Avanessoff | Yalta: a Matlab toolbox for the H∞-stability analysis of classical and fractional systems with commensurate delays[END_REF] or DDE-BIFTOOL in [START_REF] Engelborghs | Numerical bifurcation analysis of delay differential equations using dde-biftool[END_REF] require explicit characteristic quasipolynomials. However, such analytic representation can not be handled in our problem, and we are restricted to use only a numerical realization for the distributed delay. For this, a numerical algorithm for (pseudo)spectrum of the time-delay plant ( 24) is developed in this section.

This study requires the consideration of the whole RHP. However, having the explicit expression of N (θ) and therefore the behavior of φ, we can restrict our analysis to a compact subset of the RHP, which is the only subset of the RHP where equation ( 26) can have roots as detailed by Lemma 4. This will allow us to construct a finite initial grid for the evaluation of φ which we will then refine as needed. Using the Lipschitz constant of φ the algorithm can conclude on the absence of roots in a region around each evaluation point. Whenever the evaluation of φ is insufficient to guarantee the absence of roots between two adjacent evaluation points, the evaluation grid is refined around those points. This operation is repeated until the algorithm can either guarantee the absence of roots in the initial compact set (in which case we conclude that the system is stable) or a certain tolerance level in the grid size or value of | φ| is reached (in which case we cannot conclude on the stability of the system).

Let the compact subset of C:

C . = {s = z + iy ∈ C s.t. 0 ≤ z ≤ z l and |s|≤ M (z)}
, where we define [START_REF] Vazquez | Marcum Q-functions and explicit kernels for stabilization of 2×2 linear hyperbolic systems with constant coefficients[END_REF]

z l . = |ρq| r + ν M (z) . = |N (r)|e -zr + |N (0)|+ ν z (1 -e -zr ) 1 -|ρq|e -zr , with (43)
N (0) = 1 r (a + R) N (r) = 1 r (a -ρqR).
The values of a and R are given in ( 32) and ( 33

| φ(s)|= 1 -ρqe -sr - r 0 N (θ)e -sθ dθ ≥ 1 -|ρqe -sr |- r 0 N (θ)e -sθ dθ ≥ 1 -|ρq|e -zr -ν r 0 e -zθ dθ ≥ 1 -|ρq|e -zr - ν z + ν z e -zr ≥ 1 -|ρq|e -zr - ν z ≥ 1 - |ρq| zr - ν z = 1 - z l z > 0.
On the other hand, if |s|> M (z), and by integrating by parts the term

| φ(s)|≥ 1 -|ρq|e -zr - |N (r)|e -zr |s| - |N (0)| |s| - ν z|s| (1 -e -zr ) ≥ 1 - M (z) |s| > 0.
Note that M (z) is strictly decreasing with the increase of z and

lim z→0 M (z) = |N (r)|+|N (0)|+νr 1 -|ρq| and lim z→∞ M (z) = |N (0)| .
Let us first construct a grid defined as a set G 0 of nodes in a compact C containing C as depicted in Figure 1. The nodes in the grid are equally distant along the z and y axis with steps denoted by δz and δy respectively. In order to cover all the compact subset C by means of a finite number of nodes, we consider B(s, ) consisting in a closed ball centered at each s ∈ G 0 and having a radius [START_REF] Avanessoff | Yalta: a Matlab toolbox for the H∞-stability analysis of classical and fractional systems with commensurate delays[END_REF] . Proof:

having N ∈ C ∞ ([0, r]; R), φ is Lipschitz- continuous. Furthermore (50) φ(s 1 ) -φ(s 2 ) ≤ φ(s 1 ) -φ(s 2 ) ≤ max s∈ C ∂ ∂z φ(s) , ∂ ∂y φ(s) |s 1 -s 2 |, with (51)
∂ ∂y φ(s) = i ∂ ∂z φ(s) = irρqe -sr + i r 0 θN (θ)e -sθ dθ,
which makes φ Lipschitz-continuous with Lipschitz constant upper-bounded by κ.

In the light of this fact, the algorithm will start by calculating φ(s) ∀s ∈ G 0 . If ∀s ∈ G 0 , (49) is satisfied, then φ(s) > 0 ∀s ∈ C. Consequently, φ(s) > 0 ∀ Re(s) ≥ 0 and the system is stable according to the equivalent condition [START_REF] Bellman | Differential-difference equations[END_REF]. Otherwise, a more detailed grid will be constructed by considering further nodes in the vicinity of every point s ∈ G 0 where ( 49) is not verified (starting by the point where φ has the lowest value) and by dividing the radius of the closed ball in (47) by 2. We will denote a

G k (s) . = s = z + iy ∈ C s.t. z ∈ z - δz 2 k , z, z + δz 2 k and y ∈ ȳ - δy 2 k , ȳ, ȳ + δy 2 k . ( 52
)
In a refined grid G k , φ is calculated at all the newly considered nodes and a test similar to (49) will be executed. Before visiting another node in the original grid G 0 , the algorithm will track all the nodes in a refined grid using a similar technique of a depth-first search. The algorithm will backtrack a node by moving to a higherlevel grid G k-1 if all the nodes of G k verify (53) φ(s) > κ 2 k • The iteration progress is shown in Figure 3. First, we consider a point s in G 0 (blue '×'), such that φ(s) does not satisfy [START_REF] Michiels | Evaluating and approximating FIR filters : An approach based on functions of matrices[END_REF]. Subsequently, G 1 (s) is constructed (red circles). If the lowest value of φ in G 1 (s) does not satisfy (53) (for k = 1), then G 2 ( s) is constructed by adding more nodes in the vicinity of s .

= argmin s∈G 1 (s) φ(s) as shown in pink '+' signs. The algorithm will iterate in the same way to get the nodes pointed in black asterisks and cyan squares in subsequent steps.

The construction of a refined grid G k will always be priority to the node s = argmin s∈G k-1 φ(s) before backtracking the other nodes, in order to reduce the total time needed for a complete test in the case of unstable system.

The same procedure is repeated until one of the situations below occurs:

• The algorithm could visit all the points considered in G 0 (and all the refined grids) and either the inequality [START_REF] Michiels | Evaluating and approximating FIR filters : An approach based on functions of matrices[END_REF] or ( 53) is satisfied for all the created nodes. The system is therefore stable.

• The process is still iterating at a certain point, by exploring more detailed grids until φ reaches a tolerance level ε (a predefined value), then the algorithm will return that there is no guarantee for the stability of the system. A description of the algorithm is detailed in Algorithm 1. Remark 2: Let us notice that the stability analysis of ( 24) can be carried out using the Nyquist stability criterion. Indeed, the associated characteristic equation ( 26) is equivalent to 1 + Ĝ(s) = 0, where Ĝ(s) is defined in [START_REF] Damak | Stability of linear continuoustime difference equations with distributed delay: Constructive exponential estimates[END_REF]. This last equality corresponds to the characteristic equation of a feedback system whose open-loop transfer function is Ĝ(s), which is stable and strictly proper. In order to get stability for this feedback system, it is required that the Nyquist graph of Ĝ(iω) does not encircle -1. Since Ĝ(s) is strictly proper, the graph of Ĝ(iω) has to be generated only over the finite frequency interval ω ∈ [0, ω max ] where ω max is such that | Ĝ(iω)|< 1 for all ω ≥ ω max .

The construction of this frequency interval can be obtained as follows. From (45), we have, for all ω > 0,

| N (iω)|≤ |N (0)|+|N (r)|+νr ω ,
where N (0) and N (r) are defined in (43) while ν is defined in (82). Introducing [START_REF] Damak | Stability of linear continuoustime difference equations with distributed delay: Constructive exponential estimates[END_REF]. The Nyquist graph of Ĝ(iω) requires a numerical approximation of the distributed delay, that can be performed using standard techniques [START_REF] Michiels | Evaluating and approximating FIR filters : An approach based on functions of matrices[END_REF]. This algorithm leads to a simpler numerical test of stability than the one developed in Section IV-C, that can be used to test stability or instability with reduced complexity and a graphical criterion. The algorithm developed in Section IV-C uses instead analytical properties of the characteristic equation for testing stability or instability, which also points to the numerical location of (possible) poles of the plant in the righthalf complex plane. We would like to emphasize the importance of the derivation of the kernel expressions in the Appendix A as they provide the possibility of having the explicit criterion and constructing the algorithm.

ω max = |N (0)|+|N (r)|+νr 1 -|ρq| , we have that, for all ω > ω max , | N (iω)|< 1 -|ρq|, which in turn implies | Ĝ(iω)|< 1 from

V. SIMULATION RESULTS

In order to illustrate the validity of the algorithm, let us take the example of system ( 6)-( 9) with the following parameters: σ 1 = 2.3, σ 2 = -3.5, r 1 = 0.8, r 2 = 1.1, ρ = 0.5 and q = -0.7. One can check that sufficient conditions suggested in Proposition 3, [START_REF] Bastin | Stability and boundary stabilization of 1-D hyperbolic systems[END_REF]Theorem 5.4] and [START_REF] Karafyllis | Input-to-State Stability for PDEs[END_REF] fail to prove the stability of this system. Meanwhile, the algorithm described above is able to visit all the nodes and successfully verify the criterion (53) as depicted in Figure 4.

Furthermore, let us compare the stability domains obtained by the explicit criterion of Proposition 3, condition [START_REF] Auriol | Delay-robust control design for heterodirectional linear coupled hyperbolic PDEs[END_REF] required by [12, Theorem 5.4], condition (41) obtained in [START_REF] Karafyllis | Input-to-State Stability for PDEs[END_REF] and the numerical algorithm. We take the following Algorithm 1: Numerical Approximation Input: system parameters: σ 1 , σ 2 , r 1 , r 2 , ρ and q a tolerance level: ε steps: δz and δy example: r 1 = 0.5, r 2 = 1.1, σ 1 = 1 and q = 1.2. ρ is varied in order to have -1 < ρq < 1 with a step of 0.02 and σ 2 is taken between -1.6 and 1.6, with a step of 0.04. As depicted in Figure 5, if the pair (ρ , σ 2 ) is inside the dashed red curve then the system is stable according to Proposition 3. If it is inside the the solid blue curve, then it is guaranteed to be stable by [START_REF] Bastin | Stability and boundary stabilization of 1-D hyperbolic systems[END_REF]Theorem 5.4]. And if it is inside the dotted magenta curve, it is guaranteed to be stable by [START_REF] Karafyllis | Input-to-State Stability for PDEs[END_REF]. For these values, the stability domain obtained by applying the condition (41) of [START_REF] Karafyllis | Input-to-State Stability for PDEs[END_REF] is included in the domain obtained by Proposition 3. There is no inclusion between Proposition 3 and [START_REF] Bastin | Stability and boundary stabilization of 1-D hyperbolic systems[END_REF]Theorem 5.4]. As expected, the stability domain obtained by the numerical approximation algorithm (shaded area inside the black dash-dotted curve) covers all the other stability domains.

VI. CONCLUSION

In this paper, a stability analysis of a class of 2 × 2 linear first-order hyperbolic coupled PDEs is presented. A backstepping transformation yields a delay form of the original system. Then, the explicit expressions of the kernel gains of the backstepping transform are derived using the Bessel functions of the first kind. These expressions are used to obtain an explicit stability criterion in terms of the constant system coefficients, depending on their sign. The criterion is validated through academic examples. A numerical algorithm that guarantees the stability of the system in a less conservative way is developed. It works by testing if the necessary and sufficient condition is satisfied in a compact space in the right half plane. Interestingly, if we could deal with perfect computers, the algorithm would guarantee the stability with a finite number of steps. Future works will be focused on extending these results to general cases of linear first-order hyperbolic coupled PDEs, for instance, by taking into consideration diagonal in-domain couplings.

VII. ACKNOWLEDGMENTS

The authors would like to thank the reviewers for their particularly insightful and thorough comments that had an important impact in improving the paper.

APPENDIX A PROOF OF LEMMA 1

We derive the exact analytic solution for the kernel gains L ij , i, j ∈ {1, 2} satisfying the hyperbolic PDEs [START_REF] Coron | Dissipative boundary conditions for one-dimensional nonlinear hyperbolic systems[END_REF] with boundary conditions [START_REF] Prieur | Lyapunov functions for switched linear hyperbolic systems[END_REF], as functions of the constant parameters (σ 1 , σ 2 , r 1 , r 2 , ρ and q) of the original system ( 6)- [START_REF] Vazquez | Backstepping boundary stabilization and state estimation of a 2×2 linear hyperbolic system[END_REF]. Notice that the solution given in [START_REF] Vazquez | Marcum Q-functions and explicit kernels for stabilization of 2×2 linear hyperbolic systems with constant coefficients[END_REF] using the Bessel and Marcum Q-functions cannot be applied since the class of systems considered is a 2 × 2 coupled hyperbolic PDEs, in contrast to the 4 × 4 systems that appear in ( 13)- [START_REF] Prieur | Lyapunov functions for switched linear hyperbolic systems[END_REF].

We first find the integral equations satisfied by L ij before applying the successive approximation method that leads to power series expressions of these gains. Bessel functions and modified Bessel functions of the first kind are then used to give their exact analytic solution. Next, we find the analytic expressions of N (θ) and N (θ) and give upper-bounds of their absolute values. These solutions, as well as the upper-bounds, will be used to express the sufficient condition for the stability of the considered system, given in Proposition 3 in Section IV-A, and in the numerical approximation of the necessary and sufficient condition [START_REF] Bellman | Differential-difference equations[END_REF] developed in section IV-C.

A. Integral Equations

Following the method of characteristics, and integrating through the characteristic lines emerging from the boundaries [START_REF] Prieur | Lyapunov functions for switched linear hyperbolic systems[END_REF] where the kernels are defined, we find that L ij , i, j ∈ {1, 2} verify the following integral equations (54)

L 11 (x, y) = - σ 1 q r 2 1 r + r 1 σ 1 y 0 L 21 (ξ + x -y, ξ)dξ + σ 1 q r 1 r 2 (x-y) r 1 r 2 r 0 L 22 ξ r 1 + (x -y) r 1 r , - ξ r 2 + (x -y) r 1 r dξ, L 12 (x, y) = -σ 1 r 1 r 2 r + σ 1 (x-y) r 1 r 2 r 0 L 22 ξ r 1 + x r 1 r + y r 2 r , - ξ r 2 + x r 1 r + y r 2 r dξ, L 21 (x, y) = σ 2 r 1 r 2 r -σ 2 (x-y) r 1 r 2 r 0 L 11 ξ r 2 + x r 2 r + y r 1 r , - ξ r 1 + x r 2 r + y r 1 r dξ, L 22 (x, y) = σ 2 q r 2 2 r -r 2 σ 2 y 0 L 12 (ξ + x -y, ξ)dξ -σ 2 q r 2 r 1 (x-y) r 1 r 2 r 0 L 11 ξ r 2 + (x -y) r 2 r , - ξ r 1 + (x -y) r 2 r dξ.
A general form of these integral equations can be found in [START_REF] Coron | Local exponential H 2 stabilization of a 2×2 quasilinear hyperbolic system using backstepping[END_REF]. Let us write equations (54) in a compact matrix form as

(55) L(x, y) = L 0 + A[L](x, y), with L(x, y) . =     L 11 (x, y) L 12 (x, y) L 21 (x, y) L 22 (x, y)     , L 0 . =      -σ1 q r 2 1 r -σ 1 r1r2 r σ 2 r1r2 r σ 2 q r 2 2 r      .
A is a linear integral operator defined as

(56) A[L] . =     0 0 r 1 σ 1 A 1 σ1 q r1 r2 A 2 0 0 0 σ 1 A 3 -σ 2 A 4 0 0 0 -σ 2 q r2 r1 A 5 -r 2 σ 2 A 1 0 0         L 11 L 12 L 21 L 22     , with (57)
A

1 [L ij (•, •)](x, y) . = y 0 L ij (ξ + x -y, ξ)dξ, A 2 [L 22 (•, •)](x, y) . = (x-y) r 1 r 2 r 0 L 22 ξ r 1 + (x -y) r 1 r , - ξ r 2 + (x -y) r 1 r dξ, A 3 [L 22 (•, •)](x, y) . = (x-y) r 1 r 2 r 0 L 22 ξ r 1 + x r 1 r + y r 2 r , - ξ r 2 + x r 1 r + y r 2 r dξ, A 4 [L 11 (•, •)](x, y) . = (x-y) r 1 r 2 r 0 L 11 ξ r 2 + x r 2 r + y r 1 r , - ξ r 1 + x r 2 r + y r 1 r dξ, A 5 [L 11 (•, •)](x, y) . = (x-y) r 1 r 2 r 0 L 11 ξ r 2 + (x -y) r 2 r , - ξ r 1 + (x -y) r 2 r dξ.

B. Successive approximation and power series expressions

The solution of the integral equation (55) will be generated using the method of successive approximations (see, for instance, [START_REF] Carothers | An overview of the modified Picard method[END_REF] and [START_REF] Krstic | Backstepping boundary control for firstorder hyperbolic PDEs and application to systems with actuator and sensor delays[END_REF]) as

(58) L(x, y) = ∞ n=0 A n [L 0 ](x, y), with A 0 [L 0 ] = L 0 . We denote L ij k (x, y) = A k [L 0 ](x, y), k ∈ N.
We denote as well γ(x, y) . = σ 1 σ 2 (x -y)(r 2 x + r 1 y) and γ(x, y) . = σ 1 σ 2 (x -y)(r 1 x + r 2 y). In order to get the general term of the power series satisfied by each L ij (x, y) as in (58), the first few terms (L ij k (x, y), k = 1 • • • 5) have been expressed. The patterns were then obtained by induction. Adding the consecutive terms of each L ij , i, j ∈ {1, 2}, we get the following power series expressions:

(59)

L 11 (x, y) = - r 2 1 r σ 1 q ∞ m=0 (-1) m m! 2 r 2 1 r 2 r 2 γ(x, y) m + r 2 1 r 2 r 2 σ 1 σ 2 (r 2 x + r 1 y) × ∞ m=0 (-1) m m! (m + 1)! r 2 1 r 2 r 2 γ(x, y) m , L 12 (x, y) = - r 1 r 2 r σ 1 ∞ m=0 (-1) m m! 2 r 1 r 2 2 r 2 γ(x, y) m + q r 1 r 3 2 r 2 σ 1 σ 2 (x -y) × ∞ m=0 (-1) m m! (m + 1)! r 1 r 2 2 r 2 γ(x, y) m , L 21 (x, y) = r 1 r 2 r σ 2 ∞ m=0 (-1) m m! 2 r 2 1 r 2 r 2 γ(x, y) m + r 3 1 r 2 r 2 σ 1 σ 2 q (x -y) × ∞ m=0 (-1) m m! (m + 1)! r 2 1 r 2 r 2 γ(x, y) m , L 22 (x, y) = q r 2 2 r σ 2 ∞ m=0 (-1) m m! 2 r 1 r 2 2 r 2 γ(x, y) m + r 1 r 2 2 r 2 σ 1 σ 2 (r 1 x + r 2 y) × ∞ m=0 (-1) m m! (m + 1)! r 1 r 2 2 r 2 γ(x, y) m .

C. Explicit expressions

Using the Bessel functions of the first kind J 0 and J 1 [START_REF] Kreh | Bessel functions[END_REF], we have that (60)

L 11 (x, y) = - r 2 1 r σ 1 q J 0 2r 1 √ r 2 r γ(x, y) + r 1 √ r 2 r σ 1 σ 2 (r 2 x + r 1 y) x -y × J 1 2r 1 √ r 2 r γ(x, y) , L 12 (x, y) = - r 1 r 2 r σ 1 J 0 2 √ r 1 r 2 r γ(x, y) + q √ r 1 r 2 2 r σ 1 σ 2 (x -y) r 1 x + r 2 y × J 1 2 √ r 1 r 2 r γ(x, y) , L 21 (x, y) = r 1 r 2 r σ 2 J 0 2r 1 √ r 2 r γ(x, y) + r 2 1 √ r 2 qr σ 1 σ 2 (x -y) r 2 x + r 1 y × J 1 2r 1 √ r 2 r γ(x, y) , L 22 (x, y) = q r 2 2 r σ 2 J 0 2 √ r 1 r 2 r γ(x, y) + √ r 1 r 2 r σ 1 σ 2 (r 1 x + r 2 y) x -y × J 1 2 √ r 1 r 2 r γ(x, y) .
Notice that for x = y, the above expressions are well defined and simplify to (61)

L 11 (x, x) = - r 2 1 r σ 1 q + r 2 1 r 2 r σ 1 σ 2 x L 12 (x, x) = - r 1 r 2 r σ 1 L 21 (x, x) = r 1 r 2 r σ 2 L 22 (x, x) = q r 2 2 r σ 2 + r 1 r 2 2 r σ 1 σ 2 x.
Since it can be problematic to numerically manipulate (60) around x = y, we can use some properties of the Bessel functions as [START_REF] Kreh | Bessel functions[END_REF] (62) Remark 3: In the case σ 1 σ 2 < 0, it may be easier to use the modified Bessel functions of the first kind I n , having the relation J n (iξ) = i n I n (ξ), ∀n ∈ Z, ∀ξ ∈ R [START_REF] Kreh | Bessel functions[END_REF].

J n (ξ) = ξ 2n (J n-1 (ξ) + J n+1 (ξ)) ∀n ∈ Z and ∀ξ ∈ R, to get (63) L 11 (x, y) = - r 2 1 r σ 1 q + r 2 1 r 2 r 2 σ 1 σ 2 (r 2 x + r 1 y) × J 0 2r 1 √ r 2 r γ(x, y) + r 2 1 r 2 r 2 σ 1 σ 2 (r 2 x + r 1 y) × J 2 2r 1 √ r 2 r γ(x, y) , (64)
L 12 (x, y) = - r 1 r 2 r σ 1 + q r 1 r 3 2 r 2 σ 1 σ 2 (x -y) × J 0 2 √ r 1 r 2 r γ(x, y) + q r 1 r 3 2 r 2 σ 1 σ 2 (x -y) × J 2 2 √ r 1 r 2 r γ(x, y) , (65)
L 21 (x, y) = r 1 r 2 r σ 2 + r 3 1 r 2 r 2 σ 1 σ 2 q (x -y) × J 0 2r 1 √ r 2 r γ(x, y) + r 3 1 r 2 r 2 σ 1 σ 2 q (x -y) × J 2 2r 1 √ r 2 r γ(x, y) , (66)

APPENDIX B EXPRESSION AND BOUNDS OF N (θ)

A. Expression of N (θ)

The explicit criterion given in Proposition 3 is based on the expressions of N (θ) and the upper-bound of |N (θ)|, θ ∈ [0, r], as function of the original system parameters (σ 1 , σ 2 , r 1 , r 2 , ρ and q). In the sequel, we find these expressions before differentiating N (θ) and upper-bound its absolute value.

For θ ∈ [0, r 2 ], evaluating (64) and (66) at x = 1 and y = 1 -θ r 2 yields (67)

1 r 2 L 22 1, 1 - θ r 2 = q r 2 r σ 2 + R r 2 (r -θ) J 0 2 Rθ(r -θ) r + R r 2 (r -θ)J 2 2 Rθ(r -θ) r , - ρ r 2 L 12 1, 1 - θ r 2 = ρ r 1 r σ 1 -ρq R r 2 θ J 0 2 Rθ(r -θ) r -ρq R r 2 θJ 2 2 Rθ(r -θ) r ,
where R is defined in [START_REF] Niculescu | On delay-dependent stability under model transformations of some neutral linear systems[END_REF]. For θ ∈ (r 2 , r], evaluating (63) and (65) at x = 1 and y = θ -r 2 r 1 yields (68) Plugging these expressions in [START_REF] Tang | Stability analysis of coupled linear ODEhyperbolic PDE systems with two time scales[END_REF] and [START_REF] Karafyllis | On the relation of delay equations to first-order hyperbolic partial differential equations[END_REF] we obtain and a is defined in [START_REF] Diekmann | Delay equations: functional-, complex-, and nonlinear analysis[END_REF].

q r 1 L 21 1, θ -r 2 r 1 = q r 2 r σ 2 + R r 2 (r -θ) J 0 2 Rθ(r -θ) r + R r 2 (r -θ)J 2 2

B. Bounds of |N (θ)| and |N (θ)|

The following calculations will be valid for σ 1 σ 2 ≥ 0. The case where σ 1 σ 2 < 0 will be treated next. We find an upper bound of |N (θ)| using some properties of the Bessel functions of the first kind [START_REF] Kreh | Bessel functions[END_REF] such as (71) J 0 2 h(θ) ≤ 1, ∀θ ∈ [0, r] J 0 2 h(θ) + J 2 2 h(θ) ≤ 1, ∀θ ∈ (0, r).

d (θ) = -R(1 + ρq), h (θ) = R r 1 -2 θ r , J 0 2 h(θ) = 1 2 J -1 2 h(θ) -J 1 2 h(θ) = -J 1 2 h(θ) , J 2 2 h(θ) = 1 2 J 1 2 h(θ) -J 3 2 h(θ) .
Plugging (74) in (73) and using (62) yields

ν =                    1 r (|a|+R) if σ 1 σ 2 ≥ 0, |a| r I 0 |R| + |R| r I 0 |R| -I 2 |R| if σ 1 σ 2 < 0, ν =                            |R| r 2 1 + ρq + |a|+ |R| 2 if σ 1 σ 2 ≥ 0, |R| r 2 1 + ρq + |a|+ |R| 2 I 0 |R| -I 2 |R| + R 2 6r 2 I 2 |R| -I 4 |R| if σ 1 σ 2 < 0.
The values of ν and ν in (82) will be used in the construction of the algorithm in section IV-C. However, in order to be more precise, and give the integral term in the sufficient condi-

  where |•| denotes the modulus of an element in C.

  ) respectively. ν and ν, given in (82), are the upper-bounds of |N (θ)| and |N (θ)| respectively.An illustration of such a compact is shown in Figure1. The parameters chosen for this example are σ 1 = 2, σ 2 = 1.5, r 1 = 0.85, r 2 = 1.1, ρ = 0.7 and q = -0.4.Lemma 4: By definition of C, we have that ∃ > 0, such that ∀s ∈ RHP \ C, φ(s) ≥ .

Figure 1 .

 1 Figure 1. Grid G 0 construction inside C ⊃ C

r 0 N 0 N

 00 (θ)e -sθ dθ to get[START_REF] Baker | Input-output stability of linear time-invariant systems[END_REF] r (θ)e -sθ dθ = -N (r) s e -sr + N

Figure 2 .

 2 Figure 2. Closed ball B(s, )

= 1 2 5 :

 25 δz 2 + δy 2 as shown in Figure 2. Notice that C ⊂ C ⊂ ∪ s∈G 0 B(s, ) . The algorithm is based on the following Lemma. Lemma Let If at a certain point s ∈ G 0 , (49) φ(s) > κ , then φ(s) > 0 ∀s ∈ B(s, ).

Figure 3 .

 3 Figure 3. Example grid after the fourth iteration of the algorithm.

1 13 k ← k + 1 14 15 forall s ∈ G k do 16 Calculate φ(s) 17 end 18 if min s∈G k φ(s) > κ 2 k then 19 if 20 k ← k - 1 21Figure 4 .Figure 5 .

 131151617181920145 Figure 4. Main and refined grids created by the algorithm after a complete for the system with parameters: σ 1 = 2.3, σ 2 = -3.5, r 1 = 0.8, r 2 = 1.1, ρ = 0.5 and q = -0.7.

L 22 (x, y) = q r 2 2 r σ 2 + r 1 r 2 2 r 2 σ 1 σ 2 (r 1 x + r 2 y) × J 0 2 √ r 1 r 2 r γ(x, y) + r 1 r 2 2 r 2 σJ 2 2 √ r 1 r 2 r

 222222222222 1 σ 2 (r 1 x + r 2 y) × γ(x, y) .

2 2

 2 Rθ(r -θ) r .

2 J

 2 2 2 h(θ) , ∀θ ∈ [0, r],

r 2 J

 2 By remarking that |d(θ)|≤ |R|r , we have 0 2 h(θ) + J 2 2 h(θ)

2 .

 2 Using the inequalities verified by the Bessel functions of first kind:J 0 2 h(θ) + J 2 2 h(θ) ≤ 1 and J 0 2 h(θ) + 4 3 J 2 2 h(θ) + 1 3 J 4 2 h(θ) ≤ 1 and by remarking that 1 -2 θ r ≤ 1, we can bound |N (θ)| as (76) |N (θ)|≤ |R| r 2 1 + ρq + |a|+ |R|We should note that when σ 1 σ 2 < 0, |N (θ)| and |N (θ)| are bounded differently. Recalling (69) and according to Remark 3, N (θ) can be written in this case as (77)N (θ) = a r + d(θ) r 2 I 0 2 h(θ) -d(θ) r 2 I 2 2 h(θ) ,where we have defined(78) h(θ) . = |h(θ)|= |R| r 2 θ(r -θ).According to the property of the modified Bessel functions[START_REF] Kreh | Bessel functions[END_REF]: ∀n ∈ Z, I n (θ) -I n+2 (θ) ≥ 0, ∀θ ≥ 0 andI n (•) -I n+2 (•) is a strictly increasing function, |N (θ)| can be h(θ) .|N (θ)| can be bounded using the modified Bessel functions as in (79)We define ν and ν as upper-bounds of |N (θ)| and |N (θ)| respectively, as function of the system parameters, deduced from (72), (76), (79) and (81):(82)

Table I COMPARISON

 I BETWEEN THE CONDITION OF PROPOSITION 3 AND CONDITION[START_REF] Auriol | Delay-robust control design for heterodirectional linear coupled hyperbolic PDEs[END_REF] REQUIRED BY[START_REF] Bastin | Stability and boundary stabilization of 1-D hyperbolic systems[END_REF] THEOREM 5.4] FOR THE STABILITY OF SYSTEM (6)-[START_REF] Vazquez | Backstepping boundary stabilization and state estimation of a 2×2 linear hyperbolic system[END_REF]

			Condition (39)
	case	Condition required by Proposition 3 satisfied if	required by [12, Theorem 5.4]
			satisfied if
	σ 1 σ 2 > 0 or
	only one of
	them is equal
	to zero

Table II VALIDATING

 II EXAMPLES. ALL THESE SYSTEMS ARE STABLEExample of system (6)-(9)(σ 1 , σ 2 , r 1 , r 2 , ρ, q)

		Condition of Proposition 3	Condition (39) required by [12, Theorem 5.4]
	(1.1, 0.4, 1, 1.2, 0.4, -0.5)	Satisfied	Not satisfied (σ 1 σ 2 > 0)
	(-0.8, 0.7, 1, 1.2, 0.4, 0.25)	Satisfied	Satisfied
	(1.3, -0.95, 1.8, 0.44, 0.45, 0.25)	Satisfied	Not satisfied
	(1.3, -1.2, 1.8, 1.5, 0.45, 0.25)	Not satisfied	Satisfied
	(2.3, -3.5, 0.8, 1.1, 0.5, -0.7)	Not satisfied	Not satisfied

Michael Di Loreto obtained his Ph.D degree in control theory from Université de Nantes (France) in 2006. Since 2007, he is an associate professor at Laboratoire Ampère, INSA de Lyon, France. His main research interests are in time-delay systems, distributed parameter systems and linear control theory. Florent Di Meglio is tenured professor at the Centre Automatique et Systèmes of MINES ParisTech, part of PSL Research University. He received his Ph.D. from the same university in Mathematics and Control in 2011, and was a Posdoctoral Researcher at UC San Diego from 2011 to 2012. His current research interests include control and estimation design for hyperbolic PDEs, with application to process control, most notably multiphase flow control and oil drilling.

tion [START_REF] Di Meglio | Stabilization of coupled linear heterodirectional hyperbolic PDE-ODE systems[END_REF] a less conservative approximation, we express |d(θ)| depending on the sign of ρq as


Gopal Raghav

Swaminathan Gopalswamy

Karthikeyan Radhakrishnan

Jérôme Hugues

Julien Delange

Model Based Code Generation for Distributed Embedded Systems

Keywords: AADL, Architecture Driven, Distributed Embedded Software

Embedded systems are becoming increasingly complex and more distributed. Cost and quality requirements necessitate reuse of the functional software components for multiple deployment architectures. An important step is the allocation of software components to hardware. During this process the differences between the hardware and application software architectures must be reconciled. In this paper we discuss an architecture driven approach involving model-based techniques to resolve these differences and integrate hardware and software components. The system architecture serves as the underpinning based on which distributed real-time components can be generated. Generation of various embedded system architectures using the same functional architecture is discussed. The approach leverages the following technologies -IME (Integrated Modeling Environment), the SAE AADL (Architecture Analysis and Design Language), and Ocarina. The approach is illustrated using the electronic throttle control system as a case study.

INTRODUCTION

Embedded systems in ground vehicles are becoming increasingly complex in the functionality they support. Safety and security are very critical. Innovative approaches are needed to develop such systems efficiently without compromising on quality. A growing trend in development of complex embedded systems is the use of model-based development (MBD) techniques. Essentially MBD involves modeling the behavior of the embedded systems to enable simulation of the embedded system performance for various stimuli under various operating conditions. MBD supported by CAE tools facilitates the design of advanced control functionality by enabling early V&V before the mechanical and electronic hardware become available. The current state of MBD technologies is evolved enough to allow embedded software to be automatically generated from the functional models. Such tools and processes facilitate code generation for a single ECU. In practice however, as the number of processors and complexity of algorithms keep growing, there are two critical needs that emerge:

(i) The development framework needs to support modular development of embedded software promoting re-usability. Further, we need to support multiple variants in the implementation of re-usable components. This leads to the idea of an "architecture" becoming the underpinning description of a system, with variant management built around this architecture -Architecture Driven Development (ADD). These issues have already been addressed in previous work, e.g. [START_REF] Gopalswamy | Practical Considerations for the Implementation of Model Based Control System Development Processes[END_REF].

(ii) The second need is related to the fact that the functional model of an application often has a very different architecture from the architecture of the application embedded software. The functional architecture of the system corresponds to the optimum architecture required for control system development, and is concerned with the functional performance of the physical system being controlled. On the other hand, the embedded systems architecture required for any application is concerned with the number of processors, the different tasks and threads and their scheduling, etc. We need an approach that explicitly attempts to resolve such differences.

The focus of this paper is to discuss the second need above. In particular, this paper discusses an approach to extending the current technologies to allow generation of distributed embedded software from functional models, seamlessly reconciling the differences between functional architectures and embedded systems architecture.

In addition to the above consideration, a single functional architecture could often support multiple embedded systems architectures. For example, there could be technological advancements in the hardware used in the system, requiring new hardware architectures even though the functional architecture does not change.

It is also common to reuse software components between different vehicles where the functional architecture might change but the hardware architecture remains the same.

The approach discussed in this paper also addresses such additional practical issues that are of relevance during the development of complex and large embedded systems.

Section 2 discusses the architecture driven development approach for generation of distributed software. Two major activities are described in this section -(a) development of functional models consistent with a system architecture and (b) generation of distributed embedded software. Workflows to perform these activities are discussed. Section 3 discusses enabling tools and technologies that are leveraged in this process. Three major technologies are leveraged in our work -the SAE AADL (Architecture Analysis and Design Language), IME (Integrated Modeling Environment) and Ocarina (AADL toolsuite with code generation facilities). Section 4 describes a case study using the electronic throttle control application. Finally in section 5 we summarize our findings.

ARCHITECTURE DRIVEN APPROACH

The system architecture can be used as the underpinning using which functional models as well as embedded system models and software can be generated. The approach primarily consists of two major activities. First the system architecture needs to be defined and functional models (Executable Specifications) be developed that are consistent with the system architecture. Second the functional models should be integrated with the hardware architecture and embedded software generated from the complete system model.

Generation of Functional Models

The Generation of Functional Models can be captured through four key steps in our proposed approach:

(i) Functional Architecture Definition: One of the first steps is to develop the system architecture. The architecture is the topology of the system and describes the structural hierarchy of the subsystems and their interfaces and connections. Several stakeholders are involved in this step -control engineers, software engineers and managers. Usually this step is performed by the OEM based on product goals and requirements. This architecture can then be used to communicate the requirements to the suppliers of the individual components in the system, who could be a division of the OEM or an outside supplier.

(ii) Organize and Mine Component Functional Models Repository: A step that happens in tandem with the architecture is the development of component models. To enable maximum reuse, modular component models are developed over time by the organization and collected in a repository that is accessible to all authorized developers. The larger the repository of such models the quicker the development of functional applications. However, it is important to have the repository well organized and searchable. In particular, being able to search the repository based on architectural metadata of the component models (such as their interfaces, hierarchy, etc) will dramatically improve the efficiency of application development.

(iii) Associate Component Models to Architecture: Since we want the application models to be consistent with the system architecture, we need to identify component models that fit into the corporate hierarchy. As part of such an association, the ability to search the model repository for architectural metadata should be leveraged to verify architectural consistency.

Both structural and interface consistency must be performed. Additionally constraints specified in the architecture can be combined with metadata inserted into the component models to enable guided searches of the model repository.

(iv) Generate Simulatable Application Models: The final task is the generation of simulatable application models that can be used for confirming the functionality of the control system. As an example, from the system architecture we can compose models in a simulation domain such as Simulink [START_REF]Simulink is a registered trademark of The Mathworks[END_REF]. Some of the steps in the generation of functional models are shown in Figure 1.

Generation of Distributed Embedded Software

The generation of Distributed Embedded Software can be captured in these next steps:

(i) Embedded Systems Architecture Definition: The embedded system architecture is concerned with the processors in the system, the information communicated between them, the processes and threads within each of the processor, and the scheduling of those processes and threads. AADL is a powerful mechanism for describing and communicating such an architecture [START_REF][END_REF]. Thus the embedded systems architecture is developed by the stakeholders consisting of the program managers, and the embedded software and hardware engineering team.

(ii) Reconciliation of Functional Architecture with Embedded Systems Architecture: It is important to recognize that the architecture used to generate functional models is often quite different from the embedded system (hardware) architecture. What is common between the functional architecture and the embedded systems architecture is the set of component functional models associated with either architecture. Thus, in order to reconcile the functional architecture with the embedded systems architecture, a key step is to define the bindings between the component functional models in the functional architecture and the nodes of the embedded systems architecture. When such a binding is done, it is important to ensure that the functional connectivity between the software components as specified in the functional architecture is not broken, while simultaneously the communication requirements of the embedded systems architecture is satisfied.

(iii) Generation of Distributed Embedded Software: The final step in this process is to generate the distributed embedded software based on the embedded systems architecture and the associated component models.

Figure 2: Process Steps in generating (Distributed) Embedded Software for Deployment

This workflow is shown in Figure 2 above. If an automated process exists to generate the distributed embedded software from the architecture, then this provides a powerful methodology to generate different deployment variants of the distributed embedded system, all starting from the same functional architecture. The different deployments could pose different constraints, and these constraints could be accommodated easily by the appropriate binding between the embedded architecture and functional architecture. (see Figure 3 An example of multiple deployments for a given functional architecture could be as below: (i) A lowcost solution dictates that the entire application is deployed within a single processor, accepting the associated schedulability issues, and consequent performance pull-back. (ii) A medium-cost solution that allows for smart actuators and/or smart sensors to run some of the functionality in local processors, enabling faster internal feedback, and consequent improved system performance; the bulk of the functionality still runs in a central processor. (iii) A high-cost system that calls for separate processors to provide monitoring, safety and redundancy functionalities. By using the advocated approach, each of the deployment could be initiated conveniently from a common functional architecture.

below)

ENABLING TECHNOLOGIES

In the previous section we discussed how the system architecture drives the various activities involved in the embedded software generation process. Several enabling technologies are needed to efficiently realize the processes: Architecture Description Methodology: In order to perform reconciliation between the functional and hardware architectures a first requirement is that we need to be able to have standardized descriptions of both the architectures. Such an Architecture Description Language (ADL) needs to allow specification of several different components (functional, software, electronics, sensors and actuators) as well as the communication between them. Along with the architectural view that lists the topology of the system, non-functional properties of relevance (priority, time slot for CAN buses, memory capacity, etc.) need to be specified at the different nodes of the architecture. There is also a need for a well defined mechanism to bind the functional components to the embedded system components while still retaining both their properties. We chose the SAE AADL as the backbone modeling notation because it has sufficient richness in definition to be a good language for architecture driven development of embedded systems. Further the AADL is a tool neutral language and so facilitates exchange of architecture descriptions between different tools.

Architecture Driven Modeling: We need an environment where (i) Architectures are comprehended and (ii) Architecture transformations (from functional architecture to embedded system architecture) can be achieved. Such an environment should allow for linking and management of the relationships between architectures and functional models that are developed in state of the art modeling tools such as Simulink. Such an environment should also provide support for migration of the functional models to embedded system models. In this paper we use the commercial tool IME (Integrated Modeling Environment) to evaluate and demonstrate the approach.

Architecture Driven Embedded Software Creation: We need technologies that will synthesize the embedded software corresponding to the functional models, and the distributed real-time software components. For generation of the embedded software corresponding to the functional models, many of the native modeling tools themselves provide this capability. In this paper we use Simulink as the modeling tool. For generating the real time executive for distributed processor system, in this paper, we use Ocarina to evaluate and demonstrate the approach. Ocarina toolsuite provides support for code generation from AADL models for the real time executables, with appropriate integration of the software autogenerated from functional models. These technologies -AADL, IME and Ocarina facilitate the designer to go from high-level modeling down to code seamlessly, and are described in detail in the following sections.

AADL

The Architecture Analysis and Design Language (AADL) was adopted as an SAE (Society of Automotive Engineers) standard in the year 2004. The AADL is a tool neutral language that can be used to describe the run-time architecture of the embedded system. Model-based analysis for schedulability, safety, security, etc. can be performed using these descriptions [START_REF] Lewis | Multi-Dimensional Model Based Engineering for Performance Critical Computer Systems using the AADL[END_REF]. It naturally provides the ability to describe the architecture of both hardware and software components and data along with their variations with regards to implementations. AADL components are described as component types and component implementations. Component type defines the interface of the component. Component implementation inherits the properties of the component type and describes the sub-components and connections. There can be different component implementations for a single component type definition leading to variant implementations for the same component. In our work the following AADL components are used.

• System -These components can be used to describe the system architecture. They can be used to represent any software or hardware component or a combination of both. • Process -Process components are an abstraction of software responsible for scheduling and for executing threads. They execute in their own memory in a processor. • Thread -Thread components are an abstraction of software responsible for scheduling and executing sub-programs. Thread execution periods and execution times can be set as properties. • Sub-programs -Subprograms represent elementary pieces of code that process inputs to produce outputs. In the AADL only the interfaces are described. The implementations must be provided by the host language. In our work the functional models developed in Simulink provide the implementation.

• Processors -These components are abstractions of hardware and software that schedule and execute processes. Each processor will have its own clock which is the base time for all the components running on the processor. • Buses -Bus components are used to exchange data between hardware components.

[5] provides a complete description of the AADL. The functional architecture can be initially described using the AADL System components. Ports and connections define the component interfaces and communication. Once functional models are associated with the System components the architecture can be migrated to the behavioral modeling domain. The System components can later be translated to embedded components consisting of Processor, Process, Thread and Subprogram components.

Integrated Modeling Environment (IME)

IME is a model management and architecture creation and analysis tool [START_REF]IME[END_REF]. As described in Section 2 the system architecture drives the creation of functional models and generating distributed software. A visualization environment must be provided for the system architecture. Such an environment should be able to exchange architecture descriptions with other tools. Component behavior models developed over time are archived in a model repository for reuse. In order to find the consistent models their architectural information must be extracted and stored. Intelligent queries to search for certified models coupled with consistency checks facilitate selection of models. Once the selections are complete the architecture needs to be migrated to the behavioral modeling domain. As discussed in section 2.2 the functional models must be integrated seamlessly with the hardware architecture descriptions. During this step the environment must facilitate engineers to define the bindings between functional and embedded components. The bindings define the allocation of functional components to the AADL Processes which must be bound to AADL Processors. The AADL Processor components represent the different ECUs. Once the bindings are defined the embedded system architecture must be generated. From this architecture AADL descriptions need to be exported. IME is a tool that can support all the activities described above. The resulting AADL models can then be used in code generation.

Ocarina

Ocarina [START_REF] Lasnier | OCARINA: An Environment for AADL Models Analysis and Automatic Code Generation for High Integrity Applications[END_REF] is a toolsuite developed by the AADL group at Telecom ParisTech. It aims at providing AADL model manipulation, syntactic/semantic analysis, model analysis capabilities (using external tools like Cheddar [START_REF] Cheddar | [END_REF]) or embedded (e.g. generation of Petri Net models, analysis of models using the REAL constraint language).

Besides, Ocarina proposes code generation from AADL models to either C or Ada using the PolyORB-HI family of AADL runtimes. Code for the runtime, and the code generated follow carefully restrictions for the High-Integrity domains as mandated by the space, avionics or automotive domains. Targeted RTOS range from bare boards Ada runtime, real-time executive like RTEMS or RT-Linux, domain-specific OS like POK for avionics systems [START_REF] Delange | Code Generation Strategies for Partitioned Systems[END_REF] or native platforms for rapid prototyping. Ocarina supports both AADLv1.0 and AADLv2. Contrary to many MDE tools, Ocarina relies on its own internal meta-model engine, closer to a compiler AST. This allows for a wide range of internal optimizations, allowing to process large models quickly and efficiently. Ocarina is available under the GPL license and runs on most operating systems (Windows, Mac OS X, Linux). Releases, examples of application-level models, case studies as well as documentation are available on our AADL portal [10]. Ocarina has been successfully tested in the European project IST-ASSERT [11], led by the European Space Agency, and the French R&D project Flex-eWare [12], led by Thales.

Their partners, as part of advanced technology transfer projects, are currently testing Ocarina further, on industry-size case studies. One notable feature of Ocarina is the ability to include seamlessly any functional notations as implementation of blocks such as subprograms or threads. Instead of writing C code, the designer may insert a SCADE or Simulink functional blocks. Then, the Ocarina code generator will generate all the required glue code to integrate C code generated from Simulink Real-Time workshop or SCADE's kcg C code into AADL runtime code. This suppresses a tedious and error prone integration work, and allows the user to focus only on the behavior of its system instead of bothering with low-level implementation details. This feature allows for a natural bridge between a high level architecture based tool like IME and simulation on real hardware, in a distributed setup, using Ocarina generated code. We illustrate this in the next section.

CASE STUDY

In this section we describe the application of the architecture driven approach to the Electronic Throttle Control (ETC) system. ETC replaces a mechanical system consisting of a linkage between the accelerator pedal to the throttle plate. In the mechanical system the vehicle operator directly regulates the engine airflow by adjusting the position of the throttle plate via the accelerator pedal. At idle speed conditions, the airflow bypasses the throttle plate and is regulated with an Idle Air Control (IAC) valve. In the ETC system the throttle plate is actuated electronically. The desired throttle plate position (setpoint) is determined based on the pedal position as well as other inputs and operating conditions. A primary benefit of ETC is it enables system designers to incorporate throttle control into other vehicle functions, such as cruise control and vehicle stability control. ETC is considered a safety critical system. As a result a considerable portion of ETC functionality is in place for redundancy and safety monitoring. Figure 4 shows the functional architecture of one ETC system as visualized within IME. The main components in the architecture are the core controller, actuator, sensors and the plant. The core controller consists of three important functionssafety monitor, manager and servo control. In this case study we will consider three deployment scenarios for the embedded software.

Scenario1 -Single ECU: We consider the case in which a supplier is assigned the responsibility of developing the actuator system. The supplier develops a first version of the actuator based on the specifications from the OEM. The product development cycle requires virtual integration of the entire system much in advance of the availability of actual physical hardware. Therefore the supplier develops the actuator driver models and the plant models and delivers the closed loop actuator system model to the OEM. The OEM integrates the actuator models into the bigger ETC system models. The next step for the OEM is to deploy all the controller components on a single target ECU in order to eventually generate the control software.

Scenario2 -Single ECU + Smart Actuator: As an alternative to the scenario 1 above, the supplier offers a technologically advanced actuator in which the driver software is tightly integrated with the physical actuator. In turn, the OEM can offer superior functionality on some of its product lines, without disturbing the basic hardware architecture. However, in this case, the actuator driver software executes on a dedicated ECU. When the supplier delivers the actuator driver and plant models to the OEM for integration into the ETC system model, the OEM should ensure that code generation does not include the drivers with the main ECU. Scenario3 -Two ECUs + Smart Actuator: The OEM wants to reuse the same physical system and controls architecture for an advanced defense application, where safety criticality and redundancy are highly prioritized. In order to accommodate this modification the OEM wants to deploy the safety monitor component on a separate target ECU; the core functionality runs on the main ECU; the drivers run on the processor with the actuator. In all the above scenarios the functional architecture remains the same but the deployment architecture is different. We now discuss the workflow and tool chain to support the deployment in the above scenarios. Figure 5 depicts the steps involved in the generation of functional models from the system architecture. The system architecture can be developed using AADL authoring tools such as OSATE (Open Source AADL Tool Environment) [5] and imported into IME. The model repository can be mined and consistent functional models can be selected and associated with the system architecture. The functional Simulink models of the system can then be generated for further analysis. The next step is to deploy the functional models on the target. Since the functional architecture is different from the target architecture and also the target architecture is different in the three scenarios considerable re-architecting of the functional architecture must be done. The workflow for allocating the functional software components to the target is shown in figure 6. The functional system architecture can be annotated with the binding information. The binding information is different in the three deployment scenarios. The binding information must include the processor, process and thread in which each component must execute. Then the embedded architectures are generated from the functional architecture. Figures 7, 8 and 9 illustrate the embedded architectures for scenarios 1, 2 and 3 respectively. Figure 9 shows the detailed hierarchy of the embedded system. The embedded system consists of two ECUs -ECU1 and ECU2. A single process executes in ECU1 which manages a single thread. This thread executes the monitor component. The connections to the plant are not shown in the figure. In all the three target architectures the communication between the software components is maintained as specified in the functional architecture.

CONCLUSIONS

In this work we discussed the approaches and specific technologies that enable the generation of distributed embedded software from functional models. An architecture driven approach facilitates the reconciliation of the functional and embedded architectures and their integration. The specific tools and technologies also support modular development and reuse of software components. These improve the efficiency of the engineering activities. Further studies needs to be performed with regards to deploying the embedded software on the target hardware and real-time testing of the system.

Figure 1 :

 1 Figure 1: Process Steps to Generating (executable) Application Models from Functional Architecture

Figure 3 :

 3 Figure 3: A single Functional Architecture can support Multiple Deployment Architectures

Figure 4 :

 4 Figure 4: ETC Closed Loop System

Figure 5 :

 5 Figure 5: Generation of functional models consistent with system architecture.

Figure 6 :

 6 Figure 6: Generation of embedded system architectures from functional architectures

Figure 7 :Figure 8 :Figure 9 :

 789 Figure 7: Illustrates scenario 1 -Controller and actuator drivers are executed in a single ECU

Figure 10 :

 10 Figure 10: Communication between the monitor and manager components executing in separate ECUs.

Figure 11 :

 11 Figure 11: Structure of a subprogram