Leonardo F Toso
email: leonardo-felipe.toso@student-cs.fr

Giorgio Valmorbida
email: giorgio.valmorbida@centralesupelec.fr

Lyapunov Function computation for Periodic Linear Hybrid Systems via Handelman, Polya and SoS approaches: A comparative study

Keywords: Periodic linear hybrid systems, Lyapunov stability analysis, polynomial methods, convex optimization, semidefinite programming

We propose a method for the stability analysis of linear hybrid systems with periodic jumps. The method relies on the solution to polynomial inequalities based on the Handelman decomposition. Compared to existing approaches, such as sum-of-squares (SoS) and Polya's theorem, the proposed method reduces the computation time to obtain stability certificates.

INTRODUCTION

Hybrid Dynamical Systems allow to model systems composed of continuous-time processes interconnected to discrete-time digital devices. The additional complexity introduced by combining continuous and discrete-time dynamics requires the development of mathematical tools for the stability analysis to go beyond the analysis of purely continuous-or discrete-time systems [START_REF] Goebel | Hybrid dynamical systems[END_REF][START_REF] Goebel | Hybrid dynamical systems: Modeling, stability, and robustness[END_REF]; [START_REF] Cox | Hybrid internal models for spline[END_REF]. The particular case studied in this paper is the class of linear hybrid systems with periodic jumps.

A Lyapunov-based approach for stability of linear hybrid systems with periodic jumps is to impose the decrease of the Lyapunov function after the continuous-time dynamics (flows) followed by the the discrete-time evolution (jumps). This approach allows us to formulate stability analysis and the design of stabilizing periodic control laws based on the solution to differential matrix inequalities, as illustrated in Galeani et al. (2015). When using a polynomial parametrization of the functions in the differential matrix inequalities these inequalities become polynomial matrix inequalities which are required to be non-negative in a compact interval.

In a different context, the stability of nonlinear polynomial systems can be carried out by searching polynomial Lyapunov functions. The existence of these polynomial functions is a necessary and sufficient condition for local stability (over bounded sets) of exponentially stable systems, [START_REF] Peet | Exponentially stable nonlinear systems have polynomial lyapunov functions on bounded regions[END_REF]. Handelman's theorem introduced by Handelman et al. (1988) has been used as an alternative to assess the non-negativity of polynomial expressions. Handelman decomposition exploits the fact that positive polynomials can be recast in linear combinations of elements from its corresponding Handelman basis. Some works exploiting this decomposition in [START_REF] Kamyar | Constructing Piecewise-Polynomial Lyapunov Functions for Local Stability of Nonlinear Systems Using Handelman's Theorem[END_REF]; [START_REF] Briat | Robust stability and stabilization of uncertain linear positive systems via integral linear constraints: L1-gain and L2-gain characterization[END_REF]; [START_REF] Sassi | Linear relaxations of polynomial positivity for polynomial Lyapunov function synthesis[END_REF] illustrate the possibility of using Handelman theorem for the stability analysis of linear and nonlinear systems as well as to compute inner estimates to regions of attraction. Moreover, Polya's theorem can be exploited to analyze the non-negativity of polynomials Kamyar et al. (2013); Kamyar and Peet (2013). However, as [START_REF] Kamyar | Constructing Piecewise-Polynomial Lyapunov Functions for Local Stability of Nonlinear Systems Using Handelman's Theorem[END_REF] points out, the growth of the associated optimization problems on the parameters of the polynomial can prevent the application of these methods. Indeed, to construct a polynomial Lyapunov function given by a polynomial of degree six for a system composed of ten variables, the Polya's theorem approach requires to set up a semidefinite programming with ≈ 10 8 variables and ≈ 10 5 constrains.

In this paper, we search for polynomial parametrization on Handelman's bases, Polya's theorem and Sum-of-Squares decompositions, to conclude about the non-negativity of a polynomial parametrization of differential matrix inequalities.

The paper is structured as follows. After some preliminaries on definitions and notations related to the Handelman decomposition that are established in Section 2, the stability analysis of linear hybrid systems with periodic jumps is addressed in Section 3 through Lyapunov-based stability conditions recast in terms of differential matrix inequalities. Numerical examples are provided in Section 4 to illustrate the potential of the proposed approach. Finally, some conclusions and potential future work are drawn in Section 5.

Notation:

The identity matrix of dimension n is I n and the the matrix of zeros of dimension n × m is 0 n×m . The set symmetric and positive definite matrices is denoted by S n ≻0 . Periodic linear hybrid systems solutions are piecewise absolutely continuous functions X (., .) that are almost everywhere differentiable and satisfy a differential equation (flow dynamics) Ẋ = f (X ) when (t, k) ∈ T is such that t ∈ (kτ, (k + 1)τ ), and satisfy a difference equation (jump dynamics) X + = g(X ) when (t, k) ∈ T is such that t = (k + 1)τ . For brevity, the notation x [k] , x <k> may be adopted in place of x(kτ, k) and x((k + 1)τ, k) to refer to the value at the beginning and at the end, respectively, of the k th flow time interval.

PRELIMINARIES

In this section, we give conditions to check the nonnegativity of a polynomial. These conditions will be instrumental in the stability conditions of next section. Definition 1. (Scalar Polynomial): A scalar and monovariable polynomial p : R → R of degree n is defined by p(θ) = p n θ n + p n-1 θ n-1 + . . . + p 1 θ + p 0 , where p i for all i = 1, . . . , n correspond to the polynomial coefficients. A polynomial p(θ) is said to be strictly positive within an interval [START_REF] Grünbaum | Convex Polytopes[END_REF]) An n-dimensional convex polytope P ∈ R n composed of a set of K vertices V := {v j ∈ R, j = 1, . . . , K} is defined by

I = {θ ∈ R : θ ∈ [0, 1]}, if p(θ) > 0 for all θ ∈ I. Definition 2. (Convex Polytope,
P := {θ ∈ R : θ = K j=1 β j v j , K j=1 β j = 1, ∀j = 1, . . . , K}.
The interval I is an convex polytope and it can be expressed as

P := {θ ∈ R : e ⊤ i θ + f i ≥ 0, ∀ i = 1, 2}, (1)
with e 1 = 1, f 1 = 0, e 2 = -1 and f 2 = 1 being affine coefficients of the convex polytope P defined within I. Definition 3. (Handelman basis) The Handelman basis of degree d m associated to the convex polytope in (1) is defined by

B(P, dm) :={ρ(θ) : ρ(θ) = 2 j=1 (e ⊤ j θ + f j ) α j , α = (α 1 , α 2 ) ∈ N 2 , |α| 1 ≤ dm}, ={1, θ, (1 -θ), θ 2 , θ(1 -θ), (1 -θ) 2 , . . . , (1 -θ) dm , θ dm }.
From the above definitions, we formulate conditions for the non-negativity in a polytope P of a polynomial expressed in terms of the Handelman basis associated to P.

Handelman Decomposition

The following theorem exploits the Handelman basis to provide a condition for the non-negativity of a polynomial in a set P. Theorem 1. Assume that a polynomial p(θ) : R → R is strictly positive on the convex set (1), that is, p(θ) > 0, ∀θ ∈ P. There exists a degree d m , and a set of coefficients c l ∈ R ≥0 , such that p is expressed as

p(θ) := M l=1 c l K j=1 (e ⊤ j θ + f j ) αj ,
where M corresponds to the number of elements in B(P, d m ).

From the above theorem, we can conclude that if a given polynomial p(θ) is strictly positive on the set P, it can be decomposed as a sum of non-negative polynomials on P multiplied by non-negative coefficients. Given a polynomial expressed, or decomposed, in the Handelman basis of a polytope P with non-negative coefficients, it is clear that the polynomial is non-negative in P. See [START_REF] Handelman | Representing polynomials by positive linear functions on compact convex polyhedra[END_REF] for more details on the proofs of Theorem 1.

Remark 1. The existence of Handelman decomposition is only guaranteed if the polynomial is strictly positive on P. ⋆

Example 1 (Scalar Polynomial) Consider the following strictly positive polynomial in a single varaible p 1 (θ) : R → R,

p 1 (θ) = 2θ 2 - 1 4 θ + 1,
where θ ∈ P = [0, 1]. The polynomial p 1 (θ) can be decomposed into the following polynomial expression in terms of the Handelman basis

p 1 (θ) = c 1 + c 2 θ + c 3 (1 -θ) + . . . + c M (1 -θ) dm .
We then look for conditions on the coefficients c l for the non-negativity of the polynomial following the conditions of Theorem 1, obtained from the identity

c 1 + c 2 θ + c 3 (1 -θ) + . . . + c M (1 -θ) dm = 2θ 2 - 1 4 θ + 1,
The coefficients of the decomposition can then be computed by solving the corresponding system of linear equations with non-negativity constraints on the coefficients c l .

For the case d m = 3, we obtain

c 2 -c 3 + 2c 5 + c 8 -3c 9 + c 10 = - 1 4 , c 4 + c 5 + c 7 -2c 8 + 3c 9 -c 10 = 2, c 1 + c 3 + c 6 + c 9 = 1, c 6 -c 7 -c 8 -c 9 = 0.
The above set of equations and the non-negativity of c l , ∀l = 1, . . . , 10 are taken as constraints of a linear programming (LP). The LP is then solved by standard LP solvers (e.g. Gurobi Optimization, LLC (2021), MOSEK ApS (2019), Cplex (2009)). We illustrate solutions of the corresponding LPs of the above example in Table 1 for different values of the degree d m obtained with MOSEK ApS (2019).

Example 2 (Polynomial Matrix) Consider a polynomial symmetric matrix To test whether F 1 (θ) > 0, ∀θ ∈ I by using the Theorem 1, we express the positive definiteness of the matrix F 1 (θ) by the equivalent condition f p (θ,

F 1 : R → S 2×2 F 1 (θ) = 3 + 3θ θ -1 θ -1 2 .
x 1 , x 2 ) > 0 ∀x ∈ R 2 \ {0}, ∀θ ∈ I with f p (θ, x 1 , x 2 ) = x 1 x 2 ⊤ 3 + 3θ θ -1 θ -1 2 x 1 x 2 = 3x 2 1 (1 + θ) + 2x 2 x 1 (θ -1) + 2x 2 2 .
(3)

Since the above expression is homogeneous on x, it suffices to verify f p (θ,

x 1 , x 2 ) > 0 ∀x ∈ {x ∈ R 2 : |x j | ≤ 1 ∀i = 1, 2}
, which is a set containing the unit ball. Therefore, by using Theorem 1, the positive-definiteness of the F 1 (θ) on the polytope (1), can be concluded through the verification of the positivity of f p (θ, x 1 , x 2 ) on the following convex polytope

P m = {(θ, x 1 , x 2 ) ∈ R 3 : θ ≥ 0, θ -1 0, -1 < x 1 < 1, -1 < x 2 < 1}.
By constructing a Handelman basis for the polytope P m , we can proceed as in the previous example to obtain a Handelman decomposition to f p (θ, x 1 , x 2 ) via the solution of an LP. However, we can notice that to recast a polynomial matrix of dimension m in terms of a scalar expression, is necessary to add m auxiliary variables, increasing the complexity and solving time of the corresponding LP and making the approach impracticable for a high-dimensional polynomial matrices.

To avoid the increase of the number of variables describing the polynomial, we can formulate the Handelman decomposition by solving a semidefinite programming (SDP), as detailed in the following theorem. Theorem 2. Assume that a polynomial symmetric matrix F (θ) ∈ S n×n is positive definite on the convex set (1), that is F (θ) ≻ 0, ∀θ ∈ P. There exists a degree d m , and a set of positive definite matrices C l ∈ S n ≥0 , such that F (θ) can be expressed as

F (θ) = M l=1 C l K j=1 (e ⊤ j θ + f j ) αj .
where M corresponds to the number of elements in B(P, d m ).

See [START_REF] Handelman | Representing polynomials by positive linear functions on compact convex polyhedra[END_REF] for more details on Handelman decomposition.

Example 2 cont. (Polynomial Matrix) By exploiting the Theorem 2 it is possible to find positive definite matrices C l ∈ S 2 ≻0 for l = 1, . . . , M , such that F 1 (θ) is decomposed into the following expression,

F 1 (θ) = 1.
C 6 θ(1 -θ),
for d m = 2. Where the matrices C l were obtained by solving the associated SDP with MOSEK ApS (2019).

Sum-of-Squares (SoS) Decomposition

The sum-of-squares approach replaces conditions for a polynomial to be non-negative by the sufficient condition that the polynomial is decomposed into a sum of squared polynomials. Theorem 3. (Sum-of-Squares decomposition) A polynomial p(θ) : R → R is a sum-of-squares of polynomials if and only if there exist a positive semidefinite matrix G, called Gram matrix, such that,

p(θ) = z(θ) ⊤ Gz(θ),
where z(θ) = [1, θ, . . . , θ dm ] is the vector of monomials of degree up to d m .

Therefore, if we can find a vector of monomials z(θ) and a positive semidefinite matrix G, the non-negativity of p(θ) is guaranteed.

For the examples in the previous section, the scalar polynomial p 1 (θ) can be decomposed as a sum-of-squares as follows

p 1 (θ) = 1 θ ⊤ 1 -0.125 -0.125 2 1 θ
= p 2 11 (θ) + p 2 12 (θ), with p 11 (θ) = 0.9987 -0.05183θ and p 12 (θ) = 1.4133θ -0.05183. Also, the expression (3) related to the polynomial matrix F 1 (θ) can be lower bounded as f p (θ, x 1 , x 2 ) ≥ x ⊤ F lb (θ)x. The eigenvalues λ 1 and λ 2 of F lb (θ) are illustrated in Figure 1.

It thus appears evident that the polynomial p 1 (θ) is not negative and the polynomial matrix F 1 (θ) is positive semidefinite within the interval θ ∈ [0, 1]. See (Parrilo, 2000, Section 4.2) for more details on sum-of-squares programming.

In the rest of the paper, we look for a local certificate of the non-negativity of a polynomial p(θ) by searching for an SoS decomposition in the interval I by verifying whether p(θ) -(1 -θ)θg(θ), (4) is an SoS polynomial with g(θ) an SoS polynomial. If there exists an SoS decomposition of p(θ) -(1 -θ)θg(θ), it is immediate that p(θ) is not negative for all θ ∈ [0, 1], namely, the set where (1 -θ)θ ≥ 0.

Polya's Theorem

Polya's theorem, allows to check whether a polynomial p(θ) is non-negative as follows.

θ ∈ {θ ∈ R N | θ i ≥ 0, N i=1 θ i = 1}
, then there exists an integer d m , such that all coefficients of the polynomial

(θ 1 + θ 2 + . . . + θ N ) dm p(θ)
are strictly positive.

We refer the reader to [START_REF] Pólya | Inequalities[END_REF] for the proofs of Theorem 4.

Since, in this paper, we study scalar polynomials in the interval θ ∈ [0, 1] on variable θ, we shall introduce two variables θ 1 and θ 2 , thus N = 2 to describe the set [0, 1] as above. Namely we use

{(θ 1 , θ 2 ) ∈ R 2 | θ 1 ≥ 0, θ 2 ≥ 0, θ 1 = 1 -θ 2 }
and set θ = θ 1 to limit its value to the interval [0, 1]. This way, the conditions for the positivity of the scalar example given by p 1 (θ) above are that all the coefficients of the polynomial

p1 (θ) = (θ 1 + θ 2 ) dm (2θ 2 1 - 1 4 θ 1 + 1), (5)
are positive for some integer d m . For the matrix example above we generalize Polya's theorem to the verification of the positive semi-definiteness of all matrices multiplying the monomials in the polynomial matrix expression

F1 (θ) = (θ 1 + θ 2 ) dm F 1 (θ 1 ) = (θ 1 + θ 2 ) dm 3 + 3θ 1 θ 1 -1 θ 1 -1 2 for some degree d m .

PROBLEM STATEMENT

In this section, a class of linear hybrid systems with periodic time domains is introduced, and its stability analysis is formulated through differential matrix inequalities based upon Lyapunov-based conditions detailed in [START_REF] Goebel | Hybrid dynamical systems: Modeling, stability, and robustness[END_REF]; [START_REF] Cox | Hybrid internal models for spline[END_REF]Galeani et al. (2015).

Linear Hybrid Systems with periodic jumps

Consider a particular class of linear hybrid systems with all solutions defined on the same periodic hybrid time domain given by T := {(t, k) : t ∈ [kτ, (k + 1)τ ], k ∈ N}, with τ > 0 given and t denoting the current value of continuous-time and k the number of jumps already produced.

Consider the linear hybrid system, defined by

ẋ = Ax(t, k) + Bu F , (6a) x + = Ex <k> + F u J , (6b)
where x ∈ R n corresponds to the state of the linear hybrid system (6). Moreover, A ∈ R n×n , E ∈ R n×n , with u F ∈ R mc denoting the continuous-time control input and u J ∈ R m d representing the an impulsive control input.

From the theory of linear hybrid systems formulated in [START_REF] Goebel | Hybrid dynamical systems: Modeling, stability, and robustness[END_REF] (see Chapter 3), the stability conditions for the origin of the system can be established by a Lyapunov function V (x) : R n → R ≥0 satisfying positivity constraints and the decrease along the trajectories of the system. In addition, in the hybrid context imposing decrease of V (x) with respect to the flow dynamics (6a) as well as the jump dynamics (6b) with the same time-invariant quadratic function V (x) = x ⊤ P x, leads in general to restrictive constraints, as detailed in (Zaccarian et al., 2011, Section 4). In the following, the stability properties of system ( 6) are studied as timeinvariant quadratic Lyapunov functions defined for the equivalent discrete-time system.

Stability Analysis

As detailed in Galeani et al. ( 2015) the stability analysis of ( 6) can be characterised by the variation of V (x) = x ⊤ P x over one period, that is

∆ τ (V (x)) = (k+1)τ kτ V (x)dξ ∆c(V (x)) + x ⊤ [k+1] -x ⊤ ⟨k⟩ P x ⟨k⟩ ∆ d (V (x)) , (7)
where ∆ c (V (x)) and ∆ d (V (x)) correspond to the variation of V (x) over the continuous-time dynamics and after a jump, respectively. Moreover, setting θ = ξτ , the expression ( 7) can be re-written as follows

∆ τ (V (x)) = 1 0 2x ⊤ τ P ẋ + x ⊤ [k+1] P x [k+1] -x ⊤ ⟨k⟩ P x ⟨k⟩ dθ.
thus transforming the integral into the interval [0, 1] and assuming u F = u J = 0, we have the following theorem. Theorem 5. (Stability Analysis, Galeani et al. ( 2015)) Consider the linear hybrid system (6) its origin is asymptotically stable if there exists a quadratic Lyapunov function V (x) = x ⊤ P x with P ∈ S n×n >0 , matrix-valued functions R(θ), H(θ) : [0, 1] → S n×n and a matrix W ∈ R n×n , such that

d dθ R(θ) + R ⊤ (θ)A + A ⊤ R(θ) ≤ -ϵI n×n , (8a)
R(0) R(0)E E ⊤ R(0) -R(1) ≤ 0 n×n , (8b)
hold for all θ ∈ [0, 1], with R(θ) = P + H(θ), W = P and H(0) = 0.

See Galeani et al. (2015)([Propositions 2.1 and 3.1] and the references therein) for more details on the proof of Theorem 5. With the above theorem, the conditions for stability of linear hybrid systems with periodic jumps are expressed as the set of inequalities.

A solution to (8), namely matrices P and H can be obtained by imposing these matrices to be polynomials and by further imposing the inequalities to be strict. The resulting polynomial inequality can then be solved by computing the coefficients of a Handelman decomposition, following Theorem 2. It is important to notice that to exploit Theorem 2 we impose strict inequality to (8).

Remark 2. Note that the computation of V (x) = x ⊤ P x detailed in this section can also be addressed by solving the equation E ⊤ e A ⊤ τ P e Aτ E -P + Q = 0, where Q is any positive definite matrix. However, the formulation using inequalities makes it possible to cope with uncertainties, as for instance, modeled by polytopes.

⋆

NUMERICAL EXAMPLES

Numerical examples illustrate the proposed strategy to solve the matrix inequalities (8). The computation time to obtain the Handelman decomposition that gives the quadratic Lyapunov function is compared against the sumof-squares programming and Polya's theorem results. For this numerical example, the toolbox CVX, [START_REF] Grant | CVX: MATLAB Software for Disciplined Convex Programming[END_REF][START_REF] Grant | CVX: MATLAB Software for Disciplined Convex Programming[END_REF]SOSTOOLS, Papachristodoulou et al. (2021), along with the solver SeDuMi [START_REF] Sturm | Using SeDuMi 1.02, a MATLAB toolbox for optimization over symmetric cones[END_REF] were used to generate the results presented below with tolerances ϵ = 10 -6 .

For this purpose, consider a linear hybrid system with periodic jumps given by

ẋ = -1 0 0 1 x(t, k) + 1 0 u F , (9a)
x

+ = 2 0 1 1 2 x ⟨k⟩ , (9b)
u F = [-44.57 -30.05] x(t, k), (9c)
where x ∈ R 2 and a stabilizing control u F is given by the above expression.

It appears evident that, in open-loop, x 1 is continuoustime exponentially stable and discrete-time unstable, whereas x 2 is flow unstable and jump exponentially stable. Using the stabilizing continuous-time control law u F , we obtain the converging time histories of the system states, as illustrated in Figure 2.

Figure 3 shows the comparison between the time required to solve polynomial relaxations of the differential matrix inequalities in (8) for each approach: the green line concerns the results obtained by using the Handelman decomposition, whereas the blue and red line correspond to the solving time needed if we exploit the sum-of-squares and Polya's approaches, respectively. This figure illustrates the impact of increasing the maximum polynomial degree d m on the solving time of the corresponding semidefinite programming for each method. We can note that by using the sum-of-squares and Polya's approaches, the computation time required to solve the inequalities ( 8) is considerably larger than the time required by the Handelman decomposition for the same maximum polynomial degree. The improvement in solving time is even more relevant when a higher maximum polynomial degree d m is considered. Moreover, Figure 4 compares the Lyapunov function V (x) = x ⊤ P x as a function of time obtained for each one of those approaches pointing out that the solutions do not coincide, and that the monotonic decay of such solutions happens after a pair of a flow and a jump. All computations were performed on a DELL L8TQN4S laptop with Intel 2.7 GHz i7-7500U CPU and 16 GB RAM. Overall, the average improvement in solving time using Handelman decomposition rather than sum-of-squares and Polya's theorem was ≈ 79.79% and ≈ 52.58%, respectively.

CONCLUSIONS AND FUTURE WORK

The stability analysis of linear hybrid systems with periodic jumps was considered. By exploiting the fact that non-negativity analysis of polynomials can be recast in terms of Handelman decomposition, a method was proposed to solve the differential matrix inequalities related to the stability analysis of such linear hybrid systems. Numerical examples illustrate the potential of the proposed method, showing a considerable reduction in computation time compared to other approaches allowing to verify polynomial non-negativity. The potential of the proposed method should be fully exploited for systems of larger dimension for which existing methods, such as Polya's theorem in Kamyar et al. (2013) and sum-of-squares programming in [START_REF] Powers | Positive polynomials and sums of squares: Theory and practice[END_REF], are still of limited application.

Fig. 1 .

 1 Fig. 1. Eigenvalues of the lower bound matrix F lb (θ) within the interval θ ∈ [0, 1]. Theorem 4. If a polynomial p : R N → R is positive for all θ ∈ {θ ∈ R N | θ i ≥ 0,

Fig. 2 .

 2 Fig. 2. Time histories of the state of system (9) in closedloop with the proposed stabilizing controller u F .

Fig. 3 .

 3 Fig. 3. Comparison between the solving time needed to solve (8) by using Handelman decomposition, sumof-squares and Polya's theorem considering different values of d m .

Fig. 4 .

 4 Fig. 4. Comparison between the Lyapunov functions obtained to solve (8) by using Handelman decomposition, sum-of-squares and Polya's theorem considering d m = 2.

Table 1 .

 1 Handelman coefficients c l for different values of maximum polynomial degree d m .

	dm	c1	c2	c3	c5	c6	c7	c10	c14	c19
	2	0	0.75	1.00	2.00	0	0	0	0	0
	3	0.88	0	0.13	0	0.13	1.88	0	0	0
	4	0.67	1.08	0	0	0.33	0	1.00	0.33	0
	5	0.67	1.08	0	0	0	0.33	0	1.00	0.33


B Khaddoumi

H Rix

M Fereniec

NEW MATHEMATICAL METHODS IN ANALYSIS OF HIGH RESOLUTION ELECTROCARDIOGRAHY MAPPING

Recent mathematical methods, dealing with shape analysis and shape clustering are presented. The advantage of taking signal shape into account is shown in an application to High Resolution ECG Mapping: the spatial dispersion of the ECG waves recorded with a 64 lead device, and mainly the T-wave, is well correlated with the presence of Myocardial Infarcts.

Introduction

High resolution electrocardiography (HR ECG) and body surface potential mapping (BSPM) are non-invasive methods used for evaluation of electric activity of the heart. In this paper we demonstrate the usefulness of a 64-lead high-resolution ECG measurement, combined with shape analysis (SA), for providing pertinent information from the ECG. The aim of the paper is to show, that the spatial variability of ECG waves can be significantly estimated and may assess departures from the healthy state. In particular, in the case of MI patients threatened by ventricular arrhythmia compared to healthy subjects, the presence of a necrosis generates local alterations of the myocardium, leading to inhomogeneous repolarization. So we could expect to observe an increased spatial variability of the shapes, on T-waves, at the body surface.

The examination was carried out in the electrically shielded room using the high-resolution ECG measurement system. The system consists of 64 low noise amplifiers with 16-bit A/D converters (BIOSEMI, the Netherlands). Digital signals were transformed to the serial optical format and then were transferred to the computer via an optical fiber. The data acquisition was controlled by the LabView measurement software. To improve the signal-to-noise ratio the cross-correlation averaging and filtering methods were applied to 64 signals obtained from the lead position on the torso according to the University of Amsterdam lead system 2. Theory Before using the notion of signal shape in signal processing, we have first to define shape equality, then to provide a measure of shape difference and finally to show how to make statistics on shapes.

Shape equality

A shape is, for us, a class of equivalence of signals. The equivalence is associated to the affine group acting both on the left and on the right of the set of signals we are interested in. Thus, signals s(t) and v(t) are the same shape if and only if they are linked by:

) 1 ( 0 , 0 ) ( ) ( 0 > > + + = k C t t ks t v α α
In the application to ECG, the baseline is subtracted. Therefore we will assume C = 0 and equation (1) will be replaced by:

) 2 ( 0 , 0 ) ( ) ( 0 > > + = k t t
∫ ∫ ∫ = = d c t c b a t a d v d v t V d s d s t S τ τ τ τ τ τ τ τ 4 ( ) ( ) ( ) ( ; ) ( ) ( ) ( ∫
s(t) and v(t) being positive, S(t) and V(t) are increasing functions (on the supports of s(t) and v(t) respectively) which can be linked by an increasing function φ(t), modeling the change in shape between the two signals :

S V or V S e i t V S t t S t V o o 1 1 1 1 . . ) 5 ( )) ( ( ) ( )) ( ( ) ( - - - - = = = ⇔ = ϕ ϕ ϕ ϕ Practically, in the numerical form, the interval [0, 1] is sampled by M equidistant values i y : M to i y i 1 1 0 = < <
. Solving (using e.g. linear interpolation) the equation:

) 6 ( ) ' ( ) ( i i i t V t S y = = gives a set of couples ) ' , ( i i t t linked by ) ( ' i i t t ϕ = . In case of equal shape signals, i.e. if equation (2) is true, 0 ) ( t t t + = α ϕ
. The departure of function φ from linearity measures the shape difference. For example, fitting the least mean square line to the data ) ' , ( i i t t , the root mean square of the residue may be such a measure.

Averaging signals

If we want to make shape clustering, i.e. to cluster the signals according to their shapes, we need an averaging technique in addition to the possibility of measuring shape differences. The first constraint imposed to such a technique is to preserve shape when averaging equal shape signals. In the following, we show that the classical signal averaging is not a good candidate and then recall Integral Shape Averaging (ISA) and Corrected Integral Shape Averaging (CISA).

Classical Signal Averaging * Ideal case

This technique, also known as ensemble averaging is well fitted no noise reduction when the following hypotheses hold: equal shape, equal width and perfectly aligned signals. In addition the zero mean added noise and the signals are assumed to be independent. So, the signals to average can be modeled as:

1 ( ) ( ) ( ) (7 ) 1 , ( ) 1 (
) ( ) ( ) ( ) ( ) (8 )
= + = = = = + ≈
∑ So, for N sufficiently large and assuming the hypotheses, signal average gives a good estimate of s(t) within an amplitude factor. In the following this factor will be taken equal to 1.

* With jitter

When the alignment process is not able to completely suppress all the jittering effect, the model becomes:

( ) ( ) ( ) 1 ( 9
) i i i i
x t k s t d n t i to N

= -+ =

When N goes to infinity, the arithmetic average does not estimate the signal s(t) but its convolution with the probability density function (pdf) of the random delay:

( ) * ( ) (10)
( ) ( ) ( ) 1 ( 11
) i i i i i x t k s t d n t i to N α = - + =
It has been shown [START_REF] Rix | Averaging Signal with Random Time Shift and Time Scale Fluctuations[END_REF], assuming α i is a realisation of the random variable A, independent of the delay D, that (10) must be relaced by:

0 1 ( ) ( )( )( ) (12) A D t x t f s f d t τ τ τ ∞ = * ∫ Equation (12)
, where f A is the pdf of the random variable A, means that signal averaging can be viewed as a linear time-shift invariant filter followed by a linear time-scale invariant filter [START_REF] Rix | Averaging Signal with Random Time Shift and Time Scale Fluctuations[END_REF]. Each filter modifies the shape of the signal.

Integral Shape Averaging (ISA)

Let s i (t), i = to N, be N positive signals and S i (t) their normalized integrals. The ISA signal is defined [START_REF] Rix | Averaging Signal with Random Time Shift and Time Scale Fluctuations[END_REF][START_REF] Boudaoud | Integral Shape Averaging and Structural Average Estimation: A Comparative Study[END_REF] by the inverse of its normalized integral:

The ISA signal is: When all the signals s i (t) are the same shape, i.e. they are all shifted and scaled versions of the same signal, the interesting property of ISA is to give a mean signal wiyh the same shape, shifted and scaled by the average delay and the average scale respectively. In this case, the shape of the ISA signal is invariant through affine transforms of the signals to average. An illustration of ISA is given in Fig. 1. In this example the positive and the negative parts of the signals are processed separately.

( ) ( ) (14)
In the case of signals s i (t) with different shapes, ISA still provides a mean signal, but the invariance property for the affine group is not true. The algorithm of Corrected Integral Shape Averaging (CISA), proposed in [START_REF] Boudaoud | Corrected Integral Shape Averaging applied to Obstructive Sleep Apnoea detection from the Electrocardiogram[END_REF], makes it possible to extract an average shape from any set of signals.

Corrected Integral Shape Averaging

Let us call Γ CISA the normalized integral of the CISA signal γ CISA . Each S i can be linked to Γ CISA by an inceasing function:

To separate affine transformations from shape distortion, the function φ i is modeled as:

1 1 1 1 ( ) ( ) 0 1 ( 13
) N i i S y S y y N - - = = ≤ ≤ ∑ % ( ) ( ( )) (15) i CISA i

S t t

Γ ϕ =

The CISA signal is obtained by a minimization algorithm using the distance defined below.

1/2 1 1 2 1 2 1 2 0 ( , ) (( ) ( )) (17)
CISA CISA d s s n n y dy

-   = - Γ     ∫ o
The distance between s 1 and s 2 is null if the non linear distortions n 1 and n 2 are the same.

1 2 1 1 1 1 arg min ( , ) (18)
N N CISA i i i CISA i i d d s A S dt N µ µ γ - - = =     = =           ∑ ∑ o
For more details on the algorithm, see [START_REF] Boudaoud | Corrected Integral Shape Averaging applied to Obstructive Sleep Apnoea detection from the Electrocardiogram[END_REF].

Application to High Resolution ECG Maping

Shape analysis been applied to the ECG waves recorded from a BSPM device of 64 leads.

The aim was to link the spatial variability of some ECG wave, mainly the T-wave, to the presence of a myocardial infarct (MI) in patients with or without implanted cardioverter defibrillator (ICD). The application of the tools described in section 2 has been published with a first set of healthy and MI patients without ICD in [START_REF] Khaddoumi | Body Surface ECG Signal Shape Dispersion[END_REF]. In this paper we recall the principle of the procedure and the results of the application to ICD patients.

Procedure

In each column of leads, each lead is associated to the absolute value of an average wave (in order to work on positive signals). The ISA signal is computed in each column and the shape dispersion is measured by the mean shape difference of each lead in the column, using DFM, to this ISA signal. This reference is probably suboptimal, comparing to CISA, but the difference is slight and the computing time is significantly shorter with ISA.

Results

In Fig. 2, the shape dispersion averaged on each set of subjects, is shown column by column, for each ECG wave. The difference between the upper curve (non healthy) and the lower (healthy) is obvious for the T wave and specially in the precordial region. In order to show the clustering power of shape dispersion, we applied the Karhunen-Loève Transform to the two populations. The data were vectors X p (associated to each person p) of 61 components (number of active electrodes) representing the shape difference between the ECG wave recorded at the corresponding lead and the ISA Signal of its column. Let us call V 1 (respectively V 2 ) the main eigen vector of the healthy group (respectively the patient group), a 1 (p) (respectively a 2 (p)) the scalar product of X p with V 1 (respectively V 2 ). The discriminant parameter for the clustering was:

i i i i i i i i i i i i t v A t v t t m t t v t t n t A t t ϕ ω α (16) ( ) ( ( )) ( ) ( ) ( ) ( ) ( ) ( )
d(p) = a 1 (p) -a 2 (p)
The results, using T-wave shape dispersion from our first data [START_REF] Khaddoumi | Spatial shape variability of ECG waves[END_REF], have been confirmed by new data containing patients with ICD. In Fig. 3, as expected, the index d(p) is positive (except for one) for the healthy group and negative (except for three) for the patient group.

Conclusion

The presented results show that High Resolution ECG Mapping, associated to signal shape annalysis using tools fitted to statistics on shapes, makes it possible to extract valuable information from spatial dispersion of the ECG waveforms. Further works, including more data, will attempt to cluster MI patients with ICD (i.e. having an history of tacchycardia or ventricular fibrillation) from MI patients witout ICD.

=*

  With jitter and scale fluctuations Another drawback of the classical signal averaging technique, which is never mentioned in any alignment process, is the possibility, even if shape can be assumed constant, of fluctuations of the duration of the signals, modeled by time-scale fluctuations. Let us model these fluctuations, according to:

Fig. 1 .

 1 Fig.1. Illustration of Integral Shape Averaging.The signals on the left are shifted and scaled versions of the same signal; on the right this signal is completely recovered by ISA;the classical averaged signal is given for comparison.

1 0

 1

Fig. 2 .

 2 Fig. 2. Shape dispersion using ISA as the reference in each column (a) T-Waves, (b) QRS Complexes, (c) P-Waves from [Khaddoumi et al. 2006]

Fig. 3 .

 3 Fig. 3. Clustering with T-wave shape. The index d(p) for the patient group (on the left) and the healthy group (on the right)


Aurora Rossi

Pierluigi Crescenzi

Emanuele Samuel Deslauriers-Gauthier

Emanuele Natale

Natale Hyperbolic

Hyperbolic Model Captures Temporal Small Worldness of Brain Dynamics

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Hyperbolic Model Captures Temporal Small Worldness of Brain Dynamics

Aurora Rossi, Pierluigi Crescenzi, Samuel Deslauriers-Gauthier, Emanuele Natale

Brain activity can be represented as a complex network

Regular Complex Random

Linear Regression and Bandpass Filtering to remove patient, respiratory and cardiac movements

Clean BOLD signal of a voxel

Results

The Hyperbolic model not only is close to real data but it has the same behaviour

• 1 3 c = 1 ••

 11 clustering coefficient: average over all the nodes of the local clustering coefficient, measures the tendency of the neighbors of a node to form a clique c = 0 c = temporal path length: average of the fastest paths between all pairs of nodes Small Worldness: combination of high clustering coefficient (local property) and short temporal path length (global reachability property) • Scale Freeness: small number of nodes have a degree that exceeds the average one BOLD signal of a voxel Pipeline What is the best model that simulates the functional connectivity of the brain and can be used as a null model? minutes minutes t = 1 fastest path length from the green to the orange nodes = 3