Mohammad Al Khatib
email: mohammad.alkhatib@l2s.centralesupelec.fr

Antoine Girard
email: antoine.girard@l2s.centralesupelec.fr

Thao Dang
email: thao.dang@imag.fr

Stability Verification and Timing Contract Synthesis for Linear Impulsive Systems using Reachability Analysis $

Keywords: stability analysis, parameter synthesis, linear impulsive systems, difference inclusions, reachability analysis

This paper deals with stability analysis for a class of linear impulsive systems subject to a timing contract specifying bounds on the time between two consecutive impulses. We consider the problem of stability verification, which consists in proving stability for a particular timing contract, and the problem of timing contract synthesis, which consists in synthesizing a set of timing contracts that guarantee the stability of the linear impulsive system. Our approach is based on a reformulation using parameterized difference inclusions. We derive theoretical necessary and sufficient conditions for stability based on the propagation of a set by the system dynamics. For linear impulsive systems, this allow us to design a stability verification algorithm using reachability analysis. We then propose an approach to timing contract synthesis, which exploits the monotonicity of stability with respect to timing contract parameters to design an algorithm based on adaptive sampling of the parameter space. Several examples are provided, which allow us to compare our algorithm with several existing techniques, and show the effectiveness of our approach.

Introduction

Impulsive dynamical systems form a class of hybrid systems, which model processes that evolve continuously and undergo instantaneous changes at discrete time instants. Applications of impulsive dynamical systems include sampled-data control systems [START_REF] Briat | Convex conditions for robust stability analysis and stabilization of linear aperiodic impulsive and sampled-data systems under dwell-time constraints[END_REF], networked control systems [START_REF] Donkers | Stability analysis of networked control systems using a switched linear systems approach[END_REF] or multiagent systems [START_REF] Bragagnolo | Reset strategy for consensus in networks of clusters[END_REF]. In this paper, we deal with stability analysis of linear impulsive systems subject to a timing contract specifying bounds on the time between two consecutive impulses. We consider two problems of interest. The first problem under study is stability verification, which consists in proving stability of a linear impulsive system for a particular timing contract. We then consider the timing contract synthesis problem, which consists in synthesizing a set of timing contracts that guarantee stability of the system.

We use a reformulation of the linear impulsive systems in the general framework of difference inclusions. Then, for a fairly large class of difference inclusions, we establish necessary and sufficient conditions for stability. These conditions are based on the successive images of a set under the dynamics of the difference inclusion, and generalize some previous conditions on the stability of discrete-time switched systems [START_REF] Lin | Stability and stabilizability of switched linear systems: a survey of recent results[END_REF][START_REF] Athanasopoulos | Alternative stability conditions for switched discrete time linear systems[END_REF]. For linear impulsive systems, these conditions allow us to design a stability verification algorithm using reachability analysis. Then, we take advantage of previous work [START_REF] Guernic | Reachability analysis of linear systems using support functions[END_REF], which provides an efficient and accurate algorithmic scheme to compute the reachable sets for linear impulsive systems. We then propose an approach to timing contract synthesis, which exploits the monotonicity of stability with respect to timing contract parameters to design an algorithm based on adaptive sampling of the parameter space, borrowing ideas from techniques for approximating the Pareto front of a monotone multi-criteria optimization problem [START_REF] Legriel | Approximating the pareto front of multi-criteria optimization problems[END_REF][START_REF] Tendulkar | Mapping and Scheduling on Multi-core Processors using SMT Solvers[END_REF].

Some results presented in this paper appeared in preliminary form in [START_REF] Khatib | Stability verification of nearly periodic impulsive linear systems using reachability analysis[END_REF] for stability verification and in [START_REF] Khatib | Verification and synthesis of timing contracts for embedded controllers[END_REF] for parameter synthesis. This paper widen their applicability by considering a more general formulation using difference inclusions, and provides technical details that were omitted in the previous works due to space limitation. Also, a new characterization of stability based on convexification of the difference inclusion is established (Theorem 7), which provides new insights regarding the conservatism of the proposed approach.

The paper is organized as follows. In Section 2, the problems under study are formulated and the relation between linear impulsive systems and difference inclusions is formally established. Section 3 addresses the stability verification problem. We provide new theoretical necessary and sufficient conditions for stability of a class of difference inclusions, which allow us to perform stability verification for linear impulsive systems using reachability analysis. Section 4 addresses the timing contract synthesis problem by combining the stability verification algorithm and adaptive sampling of the parameter space. In Section 5, examples are used to compare our technique with existing ones and to show the effectiveness of our approach.

Related work. Several approaches have been developed in the literature for stability analysis of linear impulsive systems. A non-exhaustive list is given in Table 1. From the modeling perspective, the problem can be tackled using difference inclusions, time-delay systems or hybrid systems. On the computational side, most of the approaches are based on semi-definite programming using either Linear Matrix Inequalities (LMI) or Sum Of Squares (SOS) formulations. This makes a clear distinction with our approach which relies on reachability analysis. [START_REF] Fiacchini | Constructive necessary and sufficient condition for the stability of quasi-periodic linear impulsive systems[END_REF] is more closely related to our approach since it relies on the computation, using backward reachability analysis, of sets which are contracting between two successive impulses. In comparison, our approach uses forward reachabiltiy analysis and we compute sets that may need several impulses before contracting. In Section 5, we will provide comparisons on numerical examples between our approach and those that are listed in Table 1.

Table 1: Some of the existing approaches for stability analysis of linear impulsive systems with description of the modeling and computational approaches.

Modeling

Computation References Difference inclusions LMI [START_REF] Bauer | Networked control systems toolbox: robust stability analysis made easy[END_REF][START_REF] Hetel | Stabilization of linear impulsive systems through a nearly-periodic reset[END_REF][START_REF] Seuret | A novel stability analysis of linear systems under asynchronous samplings[END_REF] Notations. Let R, R + 0 , R + , N, N + denote the sets of reals, nonnegative reals, positive reals, nonnegative integers and positive integers, respectively. For

I ⊆ R + 0 , let N I = N ∩ I.
Let . be a norm on R n , and let B denote the associated unit ball. Given a real matrix A ∈ R n×n , A is the norm of A induced by the norm . . Given S ⊆ R n and a real matrix A ∈ R n×n , the set AS = {x ∈ R n : (∃y ∈ S : x = Ay)}; for a ∈ R, aS = (aI n )S where I n is the n × n identity matrix. The convex hull of S is denoted by conv(S). We denote the set of all subsets of R n by 2 R n . We denote by B 0 (R n ) the set of bounded subsets of R n containing 0 in their interior. For any

S ∈ B 0 (R n ), there exist c, c ∈ R + such that cB ⊆ S ⊆ cB. For p, p ∈ R d , p ≤ p if and only if p i ≤ p i , i = 1, . . . , d.

Problem formulation

This paper mainly deals with stability analysis and timing contract synthesis for linear impulsive systems. We use a general formulation based on difference inclusions and later show how linear impulsive systems can be embedded in this framework.

Difference inclusions

We consider discrete-time dynamical systems modeled by the following difference inclusion:

z k+1 ∈ Φ({z k }), k ∈ N (1)
where z k ∈ R n is the state of the system, and Φ : 2 R n → 2 R n is a set-valued map. Stability for systems of the form (1) is considered in the following sense:

Definition 1. System (1) is globally exponentially stable (GES) if there exists (C, ε) ∈ R + × (0, 1) such that for all trajectories (z k ) k∈N of (1), we have

z k ≤ Cε k z 0 , ∀k ∈ N. (2)
The first problem considered in the paper is that of verifying the stability of a difference inclusion of the form (1). We make the following assumptions on the map Φ. Assumption 1. For all S ⊆ R n , λ ∈ R + 0 , the following assertions hold:

(i) Φ(S) = z∈S Φ({z}); (ii) Φ(λS) ⊆ λΦ(S); (iii) if S is bounded, then Φ(S) is bounded.
Under item (i) of Assumption 1, for all S, S ⊆ R n , it follows that Φ(S ∪ S ) = Φ(S) ∪ Φ(S ). Also, if S ⊆ S , then Φ(S) ⊆ Φ(S ). We define the iterates of Φ as Φ 0 (S) = S for all S ⊆ R n , and Φ k+1 = Φ • Φ k for all k ∈ N. Let (z k ) k∈N be a trajectory of (1) such that z 0 ∈ S, then under item (i) of Assumption 1, for all k ∈ N, Φ k (S) is the set of all possible values of z k .

For some results of the paper, the following additional assumption related to the convexity of the map Φ is needed:

Assumption 2. For all S ⊆ R n , Φ(conv(S)) ⊆ conv(Φ(S)).
Then, the stability verification problem can be formulated as follows:

Problem 1 (Stability verification). Under Assumptions 1 and 2, verify that system (1) is GES.

The second problem considered in the paper deals with parameter synthesis for parameterized difference inclusions:

z k+1 ∈ Φ p ({z k }), k ∈ N (3)
where z k ∈ R n is the state of the system, p ∈ P ⊆ R d is a parameter, and for all p ∈ P, Φ p : 2 R n → 2 R n is a set-valued map. For a given value of the parameter p ∈ P, (3) is a difference inclusion of the form [START_REF] Khatib | Stability verification of nearly periodic impulsive linear systems using reachability analysis[END_REF]. The parameter synthesis problem consists in synthesizing a subset of parameters, which guarantee the stability of difference inclusion [START_REF] Athanasopoulos | Alternative stability conditions for switched discrete time linear systems[END_REF]. For that purpose, we shall make the following assumption on the monotonicity of Φ p with respect to parameter p: Assumption 3. For all p, p ∈ P with p ≤ p and S ⊆ R n , Φ p (S) ⊆ Φ p (S).

The previous assumption yields the following result:

Lemma 1. Let p, p ∈ P with p ≤ p ; under Assumption 3, if (3) is GES for parameter p, then (3) is GES for parameter p .

Proof. Let (z k ) k∈N be a trajectory of (3) for parameter p , then for all k ∈ N, z k+1 ∈ Φ p ({z k }) ⊆ Φ p ({z k }). Thus, (z k ) k∈N is also a trajectory of (3) for parameter p. Thus, if (2) holds for all trajectories of (3) for parameter p, it also holds for all trajectories of (3) for parameter p .

We will also make assumptions related to the boundedness and monotonicity of the set of parameters P:

Assumption 4. P = D ∩ C where D = [p m 1 , p M 1 ] × • • • × [p m d , p M d ]
and C ⊆ R d satisfies the following property:

∀p ∈ C, ∀p ∈ R d , p ≤ p =⇒ p ∈ C.
Then, the parameter synthesis problem can be formulated as follows:

Problem 2 (Parameter synthesis). Under Assumptions 3 and 4, synthesize a set P * ⊆ P such that for all p ∈ P * , (3) is GES.

Linear impulsive systems

In this section, we introduce a class of linear impulsive systems and show how they can be embedded in the formalism of the difference inclusion (3). We consider systems of the form:

ẋ(t) = A c x(t), ∀t ∈ (t k , t k+1 ], k ∈ N, (4) x
(t + k ) = A d x(t k ), k ∈ N, (5)
where (t k ) k∈N are the impulse instants, x(t) ∈ R n is the state of the system, A c and A d are n × n matrices and x(t + ) = lim τ →0,τ >0 x(t + τ ). Moreover, the sequence of impulse instants (t k ) k∈N satisfies a timing contract given by

t 0 = 0, t k+1 -t k ∈ [T , T ], k ∈ N (6)
where 0 < T ≤ T represent lower and upper bounds on the period between two consecutive impulses. Stability for linear impulsive systems of the form (4-6) is considered in the following sense: Definition 2. System (4-6) is globally uniformly exponentially stable (GUES) if there exists (C, λ) ∈ R + × R + such that, for all sequences (t k ) k∈N verifying (6) the solutions of (4-5) verify

x(t) ≤ Ce -λt x(0) , ∀t ∈ R + 0 . (7)
The stability verification problem is formulated as follows:

Problem 3 (Stability verification). Given A c , A d ∈ R n×n , 0 < T ≤ T , verify that system (4-6) is GUES.

The second problem we consider is that of timing contract synthesis. For that purpose we define a set of admissible timing contract parameters:

T = [T m , T M ] × [T m , T M ] ∩ {(T , T ) ∈ R 2 | T ≤ T }
where 0 < T m < T M and 0 < T m < T M . The timing contract synthesis problem can then be formulated as follows Problem 4 (Timing contract synthesis). Given A c , A d ∈ R n×n , and T , synthesize a set T * ⊆ T such that for all (T , T ) ∈ T * , system (4-6) is GUES.

In the remaining of the section, we show that Problems 3 and 4 can be reduced to Problems 1 and 2 respectively. To rewrite the linear impulsive system (4-6) as a parameterized difference inclusion of the form (3), we first define the notion of reachable set: Definition 3. Given a continuous-time dynamical system

ẋ(t) = Ax(t), t ∈ R + 0 , x(t) ∈ R n the reachable set on [t, t ] ⊆ R + 0 from the set S ⊆ R n is R A [t,t ] (S) = τ ∈[t,t ] e τ A S.
Let us consider a trajectory x of the linear impulsive system (4-6), and let us define the sequence z k = x(t k ) for all k ∈ N. Then, the dynamics of the sequence (z k ) k∈N is characterized by a difference inclusion of the form (3) with Φ p given by

∀S ⊆ R n , Φ p (S) = R Ac [T ,T ] (A d S), (8)
where the parameter p is defined as p = (T , -T ). The set of parameters P is then derived from T as

P = [T m , T M ] × [-T M , -T m ] ∩ {(p 1 , p 2 ) ∈ R 2 | p 1 ≤ -p 2 }. (9)
The following proposition establishes the equivalence between stability of systems (4-6) and (3).

Proposition 2. System (4-6) is GUES if and only if system (3) is GES with p = (T , -T ) and Φ p given by [START_REF] Donkers | Stability analysis of networked control systems using a switched linear systems approach[END_REF].

Proof. Let us assume that (4-6) is GUES, and let us consider a trajectory (z k ) k∈N of (3) with p = (T , -T ) and Φ p given by [START_REF] Donkers | Stability analysis of networked control systems using a switched linear systems approach[END_REF]. Then, by definition of Φ p there exists a sequence (τ k ) k∈N with τ k ∈ [T , T ] such that z k+1 = e Acτ k A d z k , for all k ∈ N. Let x be the trajectory of (4-6) corresponding to the initial condition x(0) = z 0 and the sequence of impulse instants t 0 = 0, t k+1 = t k + τ k . Then, it follows that x(t k ) = z k for all k ∈ N. By remarking that t k ≥ kT , for all k ∈ N, the stability of (4-6) gives

z k = x(t k ) ≤ Ce -λt k x(0) ≤ C(e -λT ) k z 0
which proves that system (3) is GES. Now let us assume that system (3) is GES with p = (T , -T ) and Φ p given by ( 8), and let us consider a trajectory x of (4-6). By construction, the sequence defined by z k = x(t k ), for all k ∈ N is a trajectory of system (3) with p = (T , -T ). Then, the stability of ( 1) gives for all k ∈ N,

x(t k ) = z k ≤ Cε k z 0 = Cε k x(0) . Now, let t ∈ R + , let k ∈ N be such that t ∈ (t k , t k+1 ], then t -t k ≤ T and k ≥ t/T -1. Since ε ∈ (0, 1) it follows x(t) ≤ e Ac (t-t k ) A d x(t k ) ≤ e Ac (t-t k ) A d Cε k x(0) ≤ e Ac T A d ε (t/T )-1 x(0) = e Ac T A d ε e ln(ε) T t x(0) .
Since ε ∈ (0, 1), system (4-6) is GUES.

We now show that the map Φ p and the parameter set P satisfy the assumptions formulated in Section 2.1.

Proposition 3. Let Φ p and P be given by ( 8) and ( 9). Then, for all p ∈ P, Φ p satisfies Assumptions 1 and 2. Moreover, Assumptions 3 and 4 hold.

Proof. Let us prove that the different assumptions hold. Assumption 1(i) : From the definition of reachable set, we have

Φ p (S) = τ ∈[T ,T ] e τ A A d S = τ ∈[T ,T ] z∈S e τ A A d {z} = z∈S τ ∈[T ,T ] e τ A A d {z} = z∈S Φ p ({z}).
Assumption 1(ii) : From the definition of reachable set, we have

Φ p (λS) = τ ∈[T ,T ] e τ A A d λS = λ τ ∈[T ,T ] e τ A A d S = λΦ p (S). Assumption 1(iii) : Let S ⊆ R n , then let z ∈ Φ p (S), there exists z ∈ S and τ ∈ [T , T ] such that z = e τ Ac A d z. Then, z ≤ e T Ac A d z .
Hence, if S is bounded, so is Φ p (S).

Assumption 2 : Let z ∈ Φ p (conv(S)), then there exist z ∈ conv(S) and τ ∈ [T , T ] such that z = e τ Ac A d z. Since z ∈ conv(S), there exist x, y ∈ S and λ ∈ [0, 1] such that z = λx + (1 -λ)y. Then, by linearity

z = λe τ Ac A d x + (1 -λ)e τ Ac A d y.
By remarking that e τ Ac A d x ∈ Φ p (S) and e τ Ac A d y ∈ Φ p (S), it follows that z ∈ conv(Φ p (S)). Thus, Φ p (conv(S)) ⊆ conv(Φ p (S)).

Assumption 3 : Let p, p ∈ P with p ≤ p . Then, p = (T , -T ) and

p = (T , -T ) gives [T , T ] ⊆ [T , T ]. It follows that Φ p (S) = τ ∈[T ,T ] e τ A A d S ⊆ τ ∈[T ,T ] e τ A A d S = Φ p (S). Assumption 4 : P = D ∩ C where D = [T m , T M ] × [-T M , -T m ] and C = {(p 1 , p 2 ) ∈ R 2 | p 1 ≤ -p 2 }. Moreover, let p ∈ C and p ∈ R 2 with p ≤ p, then p 1 ≤ p 1 ≤ -p 2 ≤ -p 2 which implies that p ∈ C.
It follows from Propositions 2 and 3 that Problems 3 and 4 can be reduced to Problems 1 and 2. Therefore, in the following sections, we develop algorithms to solve Problems 1 and 2.

Stability verification

This section proposes a framework to solve the stability verification problem. More precisely, it presents theoretical necessary and sufficient conditions for stability of system (1). Then, an algorithm is proposed to solve Problem 1.

Necessary and sufficient conditions for stability

The following result characterizes the stability of system (1) in terms of the map Φ.

Theorem 4. Let S ∈ B 0 (R n ), under Assumption 1, the following statements are equivalent:

(a) System (1) is GES; (b) There exists (k, j, ρ) ∈ N + × N [0,k-1] × (0, 1) such that Φ k (S) ⊆ ρΦ j (S); (c) There exists (k, ρ) ∈ N + × (0, 1) such that Φ k (S) ⊆ ρ k-1 j=0 Φ j (S). Proof. It is obvious that (b) =⇒ (c). Hence, it is sufficient to prove that (a) =⇒ (b) and (c) =⇒ (a). (a) =⇒ (b): We prove that there exists (k, ρ) ∈ N + × [0, 1) such that Φ k (S) ⊆ ρS. This is a special case of (b) when j = 0. Since S ∈ B 0 (R n ), then there exist c, c ∈ R + such that cB ⊆ S ⊆ cB. Let (z k ) k∈N be a trajectory of (1) with z 0 ∈ S, then z 0 ≤ c. Under item (i) of Assumption 1, for all k ∈ N, Φ k (S) represents all the possible values of z k . Since (1) is GES, then there exist C ∈ R + and ε ∈ (0, 1) such that z k ≤ Cε k z 0 ≤ Cε k c. This gives us for all k ∈ N, Φ k (S) ⊆ Ccε k B ⊆ C c c ε k S. For k sufficiently large, C c c ε k < 1 and therefore (b) holds. (c) =⇒ (a): Let ε = ρ 1 k ; since ρ ∈ (0, 1) then for all j ∈ N [0,k-1] , ρ ≤ ε k-j and Φ k (S) ⊆ ρ k-1 j=0 Φ j (S) ⊆ k-1 j=0 ε k-j Φ j (S). ( 10
)
Let S = k-1 j=0 ε -j Φ j (S), then using items (i) and (ii) of Assumption 1:

Φ(S ) = Φ k-1 j=0 ε -j Φ j (S) = k-1 j=0 Φ(ε -j Φ j (S)) ⊆ k-1 j=0 ε -j Φ j+1 (S) = k-2 j=0 ε -j Φ j+1 (S) ∪ ε -k+1 Φ k (S).
Making a change of index in the union and using [START_REF] Fridman | A refined input delay approach to sampled-data control[END_REF] yield

Φ(S ) ⊆ k-1 j=1 ε -j+1 Φ j (S) ∪ ε -k+1 k-1 j=0 ε k-j Φ j (S) ⊆ ε k-1 j=0 ε -j Φ j (S) = εS . (11)
Let us remark that S ⊆ S , then cB ⊆ S . In addition, since S is bounded, from item (iii) of Assumption 1, S is bounded and there exists c ∈ R + such that S ⊆ c B. Now consider a trajectory (z k ) k∈N of (1), then z 0 ∈ z 0 B ⊆ z 0 c S . Items (i) and (ii) of Assumption 1 and ( 11) give for all k ∈ N

z k ∈ Φ k z(0) c S ⊆ z(0) c ε k S ⊆ z(0) c ε k c B.
In other words, it holds for all k ∈ N,

z k ≤ c c ε k z 0 .
Since ε ∈ (0, 1) , system (1) is GES.

Theorem 4 shows the existence of a generally non-convex contracting set S , with respect to the system (1) whenever the latter is GES. In addition, when Assumption 2 holds, it is possible to show the existence of a convex contracting set as well.

Corollary 5. Let S ∈ B 0 (R n ), under Assumptions 1 and 2, system (1) is GES if and only if there exists (k, ε)

∈ N + × (0, 1) such that Φ( Ŝ) ⊆ ε Ŝ, where Ŝ = conv k-1 j=0 ε -j Φ j (S) .
Proof. For sufficiency, we assume that there exists (k, ε) ∈ N + × (0, 1) such that Φ( Ŝ) ⊆ ε Ŝ. Following the same steps after [START_REF] Fridman | Robust sampled-data stabilization of linear systems: an input delay approach[END_REF] in the proof of Theorem 4, we conclude that (1) is GES. For necessity, we assume that (1) is GES. Then, from the proof of Theorem 4, there exists (k, ε) ∈ N + × (0, 1) such that S = k-1 j=0 ε -j Φ j (S) satisfies Φ(S ) ⊆ εS . Let Ŝ = conv(S ), then using Assumption 2 we have:

Φ( Ŝ) = Φ (conv (S )) ⊆ conv (Φ (S )) ⊆ conv (εS ) = εconv (S ) = ε Ŝ.
Thus, with the addition of Assumption 2, the stable system (1) admits a convex contracting set Ŝ. In that case, we can further prove that a characterization of the stability of (1) can be given in terms of a convexified version of the set valued-map Φ. Let us consider the set-valued map Φ : 2

R n → 2 R n given by ∀S ⊆ R n , Φ(S) = conv(Φ(S)).
The images of Φ are convex sets and for all S ⊆ R n ; Φ(S) ⊆ Φ(S). The iterates of Φ are defined similarly to those of Φ. Let us also define the dynamical system associated to the set-valued map Φ:

z k+1 ∈ Φ({z k }), k ∈ N. ( 12
)
Let us state some properties of the map Φ: Lemma 6. Let Assumptions 1 and 2 hold. For all S, S ⊆ R n , λ ∈ R + 0 , the following assertions hold:

(i) if S ⊆ S , then Φ(S) ⊆ Φ(S ); (ii) Φ(λS) ⊆ λ Φ(S); (iii) if S is bounded, then Φ(S) is a bounded; (iv) Φ(conv(S)) = Φ(S).
Proof. Let us prove the different assertions.

(i) : From item (i) of Assumption 1, we have that S ⊆ S implies Φ(S) ⊆ Φ(S ). Therefore, conv(Φ(S)) ⊆ conv(Φ(S )).

(ii) : From item (ii) of Assumption 1, conv(Φ(λS)) ⊆ conv(λΦ(S)) = λconv(Φ(S)).

(iii) : From item (iii) of Assumption 1, if S is bounded then Φ(S) and thus conv(Φ(S)) are bounded.

(iv) : From the first item of the Lemma, S ⊆ conv(S) gives Φ(S) ⊆ Φ(conv(S)). Then, from Assumption 2,

Φ(conv(S)) = conv(Φ(conv(S))) ⊆ conv(conv(Φ(S))) = conv(Φ(S)) = Φ(S).
The previous result shows that items (ii) and (iii) of Assumption 1 are transferred from Φ to Φ. This is not the case of item (i) of Assumption 1, where only a weaker property can be stated for Φ (item (i) in the lemma). In particular, for k ∈ N + , Φk (S) generally contains values that are not reachable by any trajectory (z k ) k∈N of ( 12) with z 0 ∈ S. On the other hand, Assumption 2 gives a stronger property for Φ than for the original map Φ (item (iv) in the lemma).

We can now prove the following result which shows equivalence between stability of systems ( 1) and ( 12) and gives a characterization in terms of the set-valued map Φ.

Theorem 7. Let S ∈ B 0 (R n ), under Assumptions 1 and 2, the following statements are equivalent: 12) is GES.

(a) System (1) is GES; (b) There exists (k, j, ρ) ∈ N + × N [0,k-1] × (0, 1) such that Φk (S) ⊆ ρ Φj (S); (c) There exists (k, ρ) ∈ N + ×(0, 1) such that Φk (S) ⊆ ρconv k-1 j=0 Φj (S) ; (d) System (
Proof. Obviously (b) =⇒ (c). Moreover (d) =⇒ (a), since all trajectories of (1) are also trajectories of [START_REF] Herceg | Multi-Parametric Toolbox 3.0[END_REF]. Hence, it is sufficient to prove that (a) =⇒ (b) and (c) =⇒ (d).

(a) =⇒ (b): We prove that there exists (k, ρ) ∈ N + × (0, 1) such that (b) is true for j = 0. We have from Corollary 5 that there exist (i, ε) ∈ N + ×(0, 1) such that Φ( Ŝ) ⊆ ε Ŝ, where Ŝ = conv i-1 j=0 ε -j Φ j (S) . Then,

Φ( Ŝ) = conv(Φ( Ŝ)) ⊆ conv(ε Ŝ) = ε Ŝ. ( 13
)
Also S ∈ B 0 (R n ) implies, from item (iii) of Assumption 1 that Ŝ is bounded.

Then, there exists c ∈ R + such that Ŝ ⊆ cS. Let us remark that S ⊆ Ŝ, then, from [START_REF] Hetel | Stabilization of linear impulsive systems through a nearly-periodic reset[END_REF] and items (i) and (ii) of Lemma 6, for all k ∈ N,

Φk (S) ⊆ Φk ( Ŝ) ⊆ ε k Ŝ ⊆ cε k S.
Since ε ∈ (0, 1), then for k sufficiently large it becomes true that cε k < 1 which allows us to conclude.

(c) =⇒ (d): Let ε = ρ 1 k ; since ρ ∈ (0, 1) then for all j ∈ N [0,k-1] , ρ ≤ ε k-j and Φk (S) ⊆ ρconv k-1 j=0 Φj (S) ⊆ conv k-1 j=0 ε k-j Φj (S) . ( 14
)
Let Ŝ = conv( k-1 j=0 ε -j Φj (S)), then by item (iv) of Lemma 6, items (i) and (ii) of Assumption 1, we have

Φ( Ŝ ) = Φ k-1 j=0 ε -j Φj (S) = conv Φ k-1 j=0 ε -j Φj (S) ⊆ conv k-1 j=0 ε -j Φ( Φj (S)) ⊆ conv k-1 j=0 ε -j Φj+1 (S) ⊆ conv k-2 j=0 ε -j Φj+1 (S) ∪ ε -k+1 Φk (S) .
Making a change of index in the union and using [START_REF] Guernic | Reachability analysis of hybrid systems with linear continuous dynamics[END_REF] yield

Φ( Ŝ ) ⊆ εconv k-1 j=0 ε -j Φj (S) = ε Ŝ .
Let us remark that S ⊆ Ŝ , moreover, since S is bounded then from item (iii) of Lemma 6, Ŝ is bounded. It follows that there exist c ∈ R + , c ∈ R + such that c B ⊆ Ŝ ⊆ c B. Following the same steps after [START_REF] Fridman | Robust sampled-data stabilization of linear systems: an input delay approach[END_REF] of the proof of Theorem 4, one concludes that (12) is GES.

Remark 1. The results in this section can be applied to stability analysis of discrete-time switched linear systems of the form

x k+1 = A i k x k where i k ∈ N [1,N ]
, by defining the associated set-valued map Φ(S) = N i=1 A i S. In particular, by Theorem 7, we can recover the result in [17, Proposition 1] stating the equivalence between stability of the switched system and of the difference inclusion x k+1 ∈ conv({A 1 , . . . , A n })x k . Also, the stability characterizations established in [3, Theorem 1 and Corollary 2] for discrete-time switched linear systems can be obtained directly from Theorems 4 and 7, respectively.

An algorithm for stability verification

In this section, we present an algorithm for verifying the stability of system (1).

A sufficient condition for stability

The maps Φ and Φ involved in Theorems 4 and 7 can be impractical to compute exactly. This is the case for instance with linear impulsive systems, which involve the computation of the reachable set of a linear system on a time interval. In that case, we may use an over-approximation Φ : 2 R n → 2 R n , which is easier to compute and satisfies the following assumption: Assumption 5. For all S ⊆ R n , the following assertions hold:

(i) Φ(S) ⊆ Φ(S); (ii) if S is bounded then Φ(S) is bounded.
The iterates of Φ are defined similarly to those of Φ. We now derive sufficient conditions for stability of system (1) based on Φ.

Corollary 8. Under Assumptions 1 and 5, if there exist

S ∈ B 0 (R n ) and (k, i, ρ) ∈ N + × N [0,k-1] × (0, 1) such that Φ k (S) ⊆ ρΦ i (S), then system (1)
is GES. Then, by items (i) and (ii) of Assumption 1, and item (i) of Assumption 5, we have

Proof. Φ k (S) ⊆ ρΦ i (S) ⊆ ρ k-1 j=0 Φ j (S) ⊆ k-1 j=0 ε k-j Φ j (
Φ(S ) = Φ k-1 j=0 ε -j Φ j (S) ⊆ k-1 j=0 ε -j Φ(Φ j (S)) ⊆ k-1 j=0 ε -j Φ(Φ j (S)) = k-1 j=0 ε -j Φ j+1 (S).
Then, following the same steps as in [START_REF] Fridman | Robust sampled-data stabilization of linear systems: an input delay approach[END_REF], we can show that Φ(S ) ⊆ εS . Following the same lines as in the proof of Theorem 4 after (11) and using item (ii) of Assumption 5, one concludes that (1) is GES.

Let us remark that if the images of Φ are convex sets, then Φ(S) ⊆ Φ(S) ⊆ Φ(S). In such a case, in regards of Theorem 7, the only conservatism introduced by Corollary 8 is due to the over-approximation of Φ(S).

Algorithm

We propose a method to solve Problem 1 based on the sufficient condition given in Corollary 8. The stability verification algorithm consists of an initialization step and a main loop. In the initialization step, we compute an initial set S 0 ∈ B 0 (R n ), which is then propagated in the main loop using the map Φ to check the stability condition given by Corollary 8.

The choice of the initial set is important in order to try to minimize the value of the integer k such that the stability condition given by Corollary 8 holds. One approach to choose this set for the particular case of linear impulsive systems is given in Section 3.3.1. The function computing S 0 is denoted by init(Φ).

In the main loop, the initial set is propagated using the map Φ. The stability condition given by Corollary 8 is checked after each iteration. If the condition is verified then system (1) is proved GES and the algorithm returns true. We impose a maximum number of iterations k max . If that number of iterations is reached then the algorithm fails to prove stability and returns unknown. The overall method is then summarized by the Algorithm 1: Algorithm 1. Stability verification function is GES(Φ) input: Φ output: true if system (1) is proved GES, unknown otherwise parameter: k max ∈ N + 1: S 0 :=init(Φ); compute initial set 2: for k = 1 to k max do 3: The proposed approach above induces conservativeness due to the overapproximation of the map Φ and to the limited number of iterations. Consequently, it is possible that some stable systems (1) cannot be verified by the algorithm. On the other hand, if the maps Φ or Φ can be effectively computed then these can replace Φ in Algorithm 1, and for any initial set S 0 there exists a value for k max ∈ N such that the algorithm returns true if and only if system (1) is GES.

S k :=Φ(S k-1 ); set propagation 4: if ∃(i, ρ) ∈ N [0,k-1] × (0,

Case of linear impulsive systems

In this section, we give the practical details regarding the implementation of Algorithm 1 for the linear impulsive system (4-6).

We use sets given by polytopes of R n , which can be defined as the intersection of a finite number of closed half-spaces, that is S = {x ∈ R n : Hx ≤ b} where H ∈ R m×n , b ∈ R m and the vector of inequalities is interpreted componentwise.

Initial set computation

In order to try to minimize the integer k such that the stability condition given by Corollary 8 holds, we define the initial set as a symmetric polytope S 0 = {x ∈ R n : H 0 x ≤ b 0 }, which is a common contracting set for L ∈ N + linear discrete-time-invariant systems sampling the dynamics of the linear impulsive system (4-6). More precisely, S 0 satisfies for some ρ ∈ (0, 1):

∀j ∈ N [1,L] , e T j Ac A d S 0 ⊆ ρS 0 , where T j = T + (j -1)(T -T ) L . (15)
S 0 can be computed using a backward iterative method in an analogous way as done in [START_REF] Blanchini | Ultimate boundedness control for uncertain discretetime systems via set-induced Lyapunov functions[END_REF] and [START_REF] Fiacchini | Constructive necessary and sufficient condition for the stability of quasi-periodic linear impulsive systems[END_REF].

Reachability analysis

We now need to define the set valued map Φ, satisfying Assumption 5 with Φ given for all S ⊆ R n by

Φ(S) = R Ac [T ,T ] (A d S) = R Ac [0,T -T ] (e T Ac A d S).
For the computation Φ, we use efficient and accurate algorithms presented in [START_REF] Guernic | Reachability analysis of hybrid systems with linear continuous dynamics[END_REF], for over-approximating the reachable set of a linear system. For that purpose, let us introduce some notations. Given a real matrix A ∈ R n×n , |A| is the matrix whose elements are the absolute values of the elements of A.

For a set S ⊆ R n , the interval hull of S is the smallest n-dimensional interval containing the set S and is denoted by (S). The symmetric interval hull of S is the smallest symmetric (with respect to 0) n-dimensional interval containing S and is denoted by (S). Given S, S ⊆ R n , the Minkowski sum of S and S is S ⊕S = {x+x : x ∈ S, x ∈ S }. The Hausdorff distance between S and S is d(S, S ) = max(sup

x∈S inf x ∈S x -x , sup x ∈S inf x∈S x -x ). Theorem 9. [14] For δ ∈ R + , A ∈ R n×n and a set S ⊆ R n , let R A [0,δ] (S) = N i=1 R A [(i-1)h,ih] (S) (16)
where N ∈ N + , h = δ/N is the time step, and R A [(i-1)h,ih] (S) is defined by the recurrence equation:

R A [0,h] (S) = conv(S, e hA S) ⊕ 1/4 h (S), R A [ih,(i+1)h] (S) = e hA R A [(i-1)h,ih] (S), i ∈ N [1,N -1] , with h (S) = (|A| -1 (e h|A| -I) (A(I -e hA )S)) ⊕ (|A| -2 (e h|A| -I -h|A|) (A 2 e hA S)). Then, R A [0,δ] (S) ⊆ R A [0,δ] (S). If S is bounded, d(R A [0,δ] (S), R A [0,δ] (S)) = O(h).
Let us remark that h (S) in the previous theorem can be defined even if |A| is not invertible (see [14, page 56]). From the previous result, it appears that an over-approximation of Φ or Φ can be given by

Φ(S) ⊆ R Ac [0,T -T ] (e T Ac A d S), Φ(S) ⊆ conv R Ac [0,T -T ] (e T Ac A d S) . ( 17
)
These over-approximation is given by the union or convex hull of N sets which may be quite impractical for subsequent manipulations. For that reason, they will be over-approximated by a single polytope. Given a matrix H ∈ R m×n , let H i , i ∈ N [1,m] denote the row vectors of H. For a set S ⊆ R n , let us define the polytope Γ H (S) = {x ∈ R n : Hx ≤ b} where b i = sup x∈S H i x, i ∈ N [1,m] . In other words, Γ H (S) is the smallest polytope whose facets directions are given by H and containing S. Let us remark that if S is bounded and if 0 is in the interior of conv({H 1 , . . . , H m }), then Γ H (S) is bounded. In addition, if S is convex, then it can be approximated arbitrarily close by Γ H (S) by taking a sufficient number of facets directions H 1 , . . . , H m . The over-approximation of Φ and Φ is then given as follows:

Corollary 10. Let the matrix H ∈ R m×n , such that 0 is in the interior of conv({H 1 , . . . , H m }). Let Φ be given by

Φ(S) = Γ H conv R Ac [0,T -T ] (e T Ac A d S) , ( 18
)
where R Ac [0,T -T ] (e T Ac A d S) is computed as in Theorem 9. Then, Φ satisfies Assumption 5.

Proof. By ( 17) and ( 18), we have that for all S ⊆ R n , Φ(S) ⊆ Φ(S). If S is bounded, then conv R Ac [0,T -T ] (e T Ac A d S) is bounded. Furthermore, since 0 is in the interior of conv({H 1 , . . . , H m }), Φ(S) is bounded. Φ(S) defined by ( 18) is a convex set. Then, it follows from Φ(S) ⊆ Φ(S) that Φ(S) ⊆ Φ(S). Let us remark that the computation of Φ(S) is fairly simple in practice using an implementation based on support functions [START_REF] Guernic | Reachability analysis of linear systems using support functions[END_REF]. Indeed, if S is a polytope, then using the properties of support functions, the computation of Φ(S) reduces to solving a set of linear programs.

Remark 2. A normal issue arising after proposing an over-approximation Φ is to know how far the sufficient condition proposed by Corollary 8 is from being sufficient. In regards of Theorem 7, this is related to the distance between Φ(S) and Φ(S). First of all, from Theorem 9, it appears that, by choosing the time step h small enough, Φ(S) and Φ(S) can be approximated arbitrarily close by R

Parameter synthesis

In this section, we propose an approach for solving the parameter synthesis problem (Problem 2). We use the monotonicity properties stated in Assumptions 3 and 4 to design an algorithm which synthesizes a set of parameters ensuring stability of the parameterized difference inclusion (3) by adaptively sampling the parameter space.

An approximation result based on monotonicity

For p ∈ P, let us define the following property:

Stab(p) ≡ (3) is GES with parameter p.
Solving Problem 2 is equivalent to compute (a subset of) the set P o defined by P o = {p ∈ P : Stab(p)}.

We further define the following set

D o = {p ∈ D : (p / ∈ C) ∨ ((p ∈ C) ∧ Stab(p))} .
One can easily check that the following relation holds:

P o = C ∩ D o . (19)
Hence, from the previous equality, we can solve Problem 2 by computing (a subset of) the set D o . Moreover, D o satisfies the following monotonicity property:

Proposition 11. Under Assumptions 3 and 4, for all p, p ∈ D, the following implications hold:

((p ≤ p ) ∧ (p ∈ D o )) =⇒ p ∈ D o . ((p ≤ p ) ∧ (p / ∈ D o )) =⇒ p / ∈ D o .
Proof. Let us assume p ≤ p and p ∈ D o . There are two cases: This proves the first implication. For the second implication, it is sufficient to check that

((p ≤ p ) ∧ (p ∈ D o )) =⇒ p ∈ D o ≡ ¬(p ≤ p ) ∨ (p / ∈ D o ) ∨ (p ∈ D o ) ≡ ((p ≤ p ) ∧ (p / ∈ D o )) =⇒ p / ∈ D o .
The previous property is instrumental for computing a subset of D o since it allows us to state the following theorem:

, it follows that C * is a solution to Problem 2 and C o ⊆ C ∩ D.

An algorithm for parameter synthesis

Theorem 12 shows that it is possible to compute an under and overapproximation of the set D o by sampling the parameter space D. In this section, we use this property to design a parameter synthesis algorithm. Similar algorithms have been used in [START_REF] Legriel | Approximating the pareto front of multi-criteria optimization problems[END_REF][START_REF] Tendulkar | Mapping and Scheduling on Multi-core Processors using SMT Solvers[END_REF] for computing an approximation of the Pareto front of a monotone multi-criteria optimization problem. Indeed, this latter problem can be tackled by computing an under and over-approximation of a set satisfying a monotonicity property similar to that of Proposition 11.

Algorithm 2 computes an under-approximation D and an over-approximation D of the set D o by sampling iteratively the parameter space D. Lines 9 to 14 describe the main loop of the parameter synthesis algorithm. At any time of the execution, D ⊆ D o ⊆ D holds. We pick a sample p ∈ D\D which is the unexplored parameter region lying in the over-approximation of D o but not in its under-approximation. If p ∈ D o (or if p / ∈ D o ), then we update the under-approximation D (or the over-approximation D) according to Theorem 12. The algorithm stops when the Hausdorff distance between the D and D becomes smaller than ε. Of course, the choice of the sample p ∈ D \ D, at line 10, is crucial for the efficiency of the algorithm. In our implementation of the algorithm, we use the selection criteria proposed in [START_REF] Legriel | Approximating the pareto front of multi-criteria optimization problems[END_REF] which consists in choosing the sample that will produce the fastest decrease of the Hausdorff distance d(D, D). In [START_REF] Tendulkar | Mapping and Scheduling on Multi-core Processors using SMT Solvers[END_REF] an alternative selection criteria based on multiscale grid exploration was proposed. It is important to note that Algorithm 2 needs testing if the samples p ∈ D o , which requires checking the condition Stab(p). In our implementation, this is done using Algorithm 1 since we assume furthermore that Φ p satisfies Assumptions 1 and 2. If it returns true, then we can consider that Stab(p) holds. If it returns unknown, we treat the sample as if Stab(p) is false. As a consequence, in practice it may be the case that D is not an over-approximation of D o . However, it always holds that D ⊆ D o and therefore the set returned by Algorithm 2 is always a valid solution to Problem 2. Note that the property Stab(p) need not be checked using Algorithm 1 but one can use any of existing stability verification algorithm.

Applications and Numerical Results

We implemented the algorithms presented in this paper in Matlab using the Multi-Parametric Toolbox [START_REF] Herceg | Multi-Parametric Toolbox 3.0[END_REF]. All reported experiments were realized on a desktop with i7 4790 processor of frequency 3.6 GHz and a 8 GB RAM.

An academic example

Example 1. The following example is taken from [START_REF] Hetel | Stabilization of linear impulsive systems through a nearly-periodic reset[END_REF] proposing an LMI-based approach to verify stability of a linear impulsive system. Consider system (4-6) with

A c =   0 -3 1 1.4 -2.6 0.6 8.4 -18.6 4.6   , A d =   1 0 0 0 1 0 0 0 0   (20)
As noted in [START_REF] Hetel | Stabilization of linear impulsive systems through a nearly-periodic reset[END_REF], the matrix i∈N [START_REF] Khatib | Stability verification of nearly periodic impulsive linear systems using reachability analysis[END_REF][START_REF] Blanchini | Ultimate boundedness control for uncertain discretetime systems via set-induced Lyapunov functions[END_REF] (e T i Ac A d ) has eigenvalues outside the unit circle for T 1 = 0.515 and T i = 0.1, for i ∈ N [START_REF] Khatib | Verification and synthesis of timing contracts for embedded controllers[END_REF][START_REF] Blanchini | Ultimate boundedness control for uncertain discretetime systems via set-induced Lyapunov functions[END_REF] . As a result, we can consider that if T = 0.1, the value 0.515 is an upper bound for admissible values of T . For T = 0.1, stability could be proven up to T = 0.3 following the LMI approach in [START_REF] Hetel | Stabilization of linear impulsive systems through a nearly-periodic reset[END_REF], and up to T = 0.375 following the set based approach in [START_REF] Fiacchini | Constructive necessary and sufficient condition for the stability of quasi-periodic linear impulsive systems[END_REF]. Results obtained using Algorithm 1 with several parameter setups are reported in Table 2. In this example, parameter setups B and D lead to less conservative results than the mentioned approaches since stability is verified at least up to T = 0.5. Moreover, with parameter setup D, the verified value T = 0.514 is tight, since it is very close to the known upper-bound 0.515. [START_REF] Seuret | A novel stability analysis of linear systems under asynchronous samplings[END_REF] for several values of parameters L (number of subsystems chosen to find the initial set S 0 ) and k max (maximum number of iterations of Algorithm 1) with N = 100 (number steps used in reachability analysis): for T = 0.1, maximum value of T for which stability could be proved; T CP U is the computation time in seconds; i, k are the index values for which the stability condition [START_REF] Guernic | Reachability analysis of linear systems using support functions[END_REF] and [START_REF] Liu | Stability of linear systems with general sawtooth delay[END_REF]. Figure 1 shows the polytopes S 0 and S 1 computed by Algorithm 1 using parameter setup C for T = 0.1 and T = 0.5. The inclusion of S 1 in S 0 proves the stability of the linear impulsive system. We now consider the timing contract synthesis problem for system [START_REF] Seuret | A novel stability analysis of linear systems under asynchronous samplings[END_REF]. We used Algorithm 2 with parameters ε = 0.01, T m = T m = 10 -4 and T M = T M = 0.59. Stability for sampled values of the timing contract parameter space was verified using Algorithm 1 were L = 2, k max = 30 and the number of steps used in reachability analysis was N = 100. Figure 2 shows the region T * of timing contract parameters, for which GUES is guaranteed. The total number of parameter samples was 179 with a computation time of about 67 seconds.

S k ⊆ int(S i ) is verified; m is such that H = H 0 ∈ R m×3 in computing

Sampled-data systems

We consider the problem of verifying stability of aperiodic sampled-data control systems. These systems are given under the form:

ż(t) = Az(t) + Bu(t), ∀t ∈ R + u(t) = Kz(t k ), ∀t ∈ (t k , t k+1 ], k ∈ N (21)
where t k+1 -t k is a variable sampling interval bounded in [T , T ], z(t) ∈ R n is the state of the system, u(t) ∈ R m is the control input computed aperiodically at instants t k , and K ∈ R m×n is the feedback gain. This problem can be rewritten in the form (4-6), with:

A c = A B 0 0 , A d = I n 0 K 0 , x(t) = z(t) u(t k ) (22)
with I n as the n × n identity matrix, and with the same timing contract parameters T and T .

Example 2. This sampled data system is taken from [START_REF] Briat | Convex conditions for robust stability analysis and stabilization of linear aperiodic impulsive and sampled-data systems under dwell-time constraints[END_REF], which compares results of LMI or SOS based approaches for stability analysis of linear impulsive systems. Consider the state space plant model given by [START_REF] Seuret | Stability analysis of sampled-data systems using Sum of Squares[END_REF] with

A = 0 1 0 -0.1 , B = 0 0.1 , K = -3.75 -11.5 . (23)
After rewriting the problem in the form of (4-6) with matrices defined as in [START_REF] Tendulkar | Mapping and Scheduling on Multi-core Processors using SMT Solvers[END_REF], we set T = 10 -5 . For this system, we can check numerically that the matrix e τ Ac A d is Schur for τ ∈ [0, 1.7294] and has eigenvalues outside the unit circle for larger values of τ . Thus, we know that 1.7294 is an upper bound for the maximal value of T guaranteeing stability. Table 3 reports the maximum value of T , for which stability could be verified by our approach and other existing methods, as reported in [START_REF] Briat | Convex conditions for robust stability analysis and stabilization of linear aperiodic impulsive and sampled-data systems under dwell-time constraints[END_REF]. The results obtained by our approach are similar to the least conservative result reported in [START_REF] Briat | Convex conditions for robust stability analysis and stabilization of linear aperiodic impulsive and sampled-data systems under dwell-time constraints[END_REF], which was obtained by the method presented in [START_REF] Seuret | Stability analysis of sampled-data systems using Sum of Squares[END_REF]. More precisely, stability could be proven up to T = 1.7294 using Algorithm 1 with parameters L = 2, k max = 1 and N = 1011. This shows the tightness of our approach since we know that the system becomes unstable for T > 1.7294 Note that a matrix H = H 0 ∈ R 10×3 is used in [START_REF] Liu | Stability of linear systems with general sawtooth delay[END_REF] and the computation time was 0.1207 seconds.

Example 3. The second sampled-data control system is also taken from [START_REF] Briat | Convex conditions for robust stability analysis and stabilization of linear aperiodic impulsive and sampled-data systems under dwell-time constraints[END_REF], with:

A = 0 1 -2 0.1 , B = 0 1 K = 1 0 . ( 24
)
We set T = 0.4. Note that the system becomes unstable for T = 1.889 since the matrix i∈N [START_REF] Khatib | Stability verification of nearly periodic impulsive linear systems using reachability analysis[END_REF][START_REF] Khatib | Verification and synthesis of timing contracts for embedded controllers[END_REF] (e T i Ac A d ) has eigenvalues outside the unit circle for T 1 = 0.4 and T 2 = 1.889. Results obtained by our approach and by several others are also reported in Table 3. Our approach has better results than the existing ones since it was able to verify stability for the system up to T = 1.888, instead of T = 1.828 for the method presented in [START_REF] Seuret | Stability analysis of sampled-data systems using Sum of Squares[END_REF]. Again, our approach appears to be quite tight since the maximal value of T for which stability was verified is very close to the known upper bound 1.889. Algorithm 1 was used with parameters L = 2, k max = 30 and the number of time steps used for the over-approximation of the reachable set is N = 100. The stability condition S k ⊆ int(S i ) was verified for k = 14 and i = 12. Also, a matrix H = H 0 ∈ R 18×3 is used in [START_REF] Liu | Stability of linear systems with general sawtooth delay[END_REF] and the computation time was 0.824 seconds. We now consider the timing contract synthesis problem for the sampleddata system given by matrices (24). We used Algorithm 2 with parameters ε = 0.01, T m = T m = 0.2101 and T M = T M = 2.0201. Parameters of Algorithm 1 for stability verification were L = 2, k max = 30 and N = 100. Figure 3 shows the region T * of timing contract parameters (T , T ), for which GUES is guaranteed. The total number of parameter samples was 83 with a computation time of about 77 seconds.

Example 4. We consider the state space plant model of a batch reactor [START_REF] Donkers | Stability analysis of networked control systems using a switched linear systems approach[END_REF], with a static feedback, given by [START_REF] Seuret | Stability analysis of sampled-data systems using Sum of Squares[END_REF]

and one of the two feedback gains K 1 = 0.13 0.02 0.07 -0.18 1.21 0.28 0.48 -0.06 , K 2 = 0.41 -0.45 0.38 -0.59 1.65 -0.2 0.91 -0.59 . We rewrite the problem in the form of a 6-dimensional impulsive system (4-6). Then, we apply Algorithm 1 to check stability of the impulsive system. We compare our results to those obtained using the NCS toolbox [START_REF] Bauer | Networked control systems toolbox: robust stability analysis made easy[END_REF] in Table 4. After setting T = 0.01, we report the maximal value of T for which stability has been verified. Note that we conducted an extra experiment labelled "(exp1)" to compare the results in terms of CPU time after fixing the same values of T in Algorithm 1 and the NCS toolbox. The NCS toolbox uses three different approximation methods to embed the timing uncertainty (Jordan Normal Form (JNF), Cayley Hamilton, and Gridding and Norm Bounding (GNB)), so the experiments are conducted as follows: we search for the maximum value of T that guarantees stability by running experiments using the three approximation methods. Then we report the computation time for the experiment in which we obtained this maximum value. In case the maximum bound could be obtained by more than one experiment, we report the CPU time corresponding to the fastest in terms of computation. Stability for system (25) with the feedback gains K 1 and K 2 is guaranteed using the GNB approximation with 65 gridpoints. Parameter setups used by Algorithm 1 are summarized by Table 5. It is clear, for the systems at hand, that our method is competitive with the NCS toolbox in terms of CPU time and tightness, since Algorithm 1 yields better results for system (25) with the feedback gain K 1 and has quite similar results for the same system with controller K 2 . Notice that in this example the dimension of the problem increased where matrices H = H 0 ∈ R 88×6 and H = H 0 ∈ R 80×6 were used in computing [START_REF] Liu | Stability of linear systems with general sawtooth delay[END_REF] for the former and latter results respectively. For gain K 2 , we could not obtain a better value for T since the Matlab implementation of Algorithm 1 ran into numerical problems when increasing the parameters L or k max or N . L N k max K 1 Algorithm 1(exp1) 3 200 100 Algorithm 1(exp2) 5 100 100 K 2 Algorithm 1 4 100 100

Conclusion

In this work, we proposed a new approach to stability analysis of linear impulsive systems using reachability analysis. We have also shown that the stability verification algorithm can be combined with adaptive sampling of the parameter space for synthesizing timing contracts. The effectiveness of our approach has been shown on several examples, where it has been compared with previously existing techniques. Applications to aperiodic sampled data control systems have been shown. For this class of systems, the timing contract synthesis is of particular interest since it can be used by control and software engineers to derive requirements that must be met by the real-time implementation of a control law. As future work, it would be interesting to handle the problem of controller synthesis given a timing contract, and to co-synthesize the controller and the timing contract parameters.

S) where ε = ρ 1 k

 1 . Similar to the second part of the proof of Theorem 4, let S = k-1 j=0 ε -j Φ j (S).

Ac [ 0

 0 ,T -T ] (e T Ac A d S) and conv R Ac [0,T -T ] (e T Ac A d S) , respectively. Then, conv R Ac [0,T -T ] (e T Ac A d S) can be approximated arbitrarily close by Γ H conv(R Ac [0,T -T ] (e T Ac A d S)) by considering a sufficient number of approximation directions H i . Thus, it follows that, by choosing appropriately the time step h and the matrix H, Φ(S) can be approximated arbitrarily close by Φ(S).

1 .

 1 Following Assumption 4, if p / ∈ C, then p / C and therefore p ∈ D o . 2. If p ∈ C and Stab(p), then either p / ∈ C which implies p ∈ D o , or p ∈ C. In the latter case Stab(p ) follows from Assumption 3, then p ∈ D o .

Theorem 12 .= M 1 j=1 2 j=1

 1212 Let p 1 , . . . , p M 1 ∈ D o , and p 1 , . . . , p M 2 ∈ D \ D o and let D {p ∈ D : p j ≤ p}, D = D \ M {p ∈ D : p ≤ p j }. Then, D ⊆ D o ⊆ D. Moreover, C * = C ∩ D is a solution to Problem 2 and C o ⊆ C ∩ D. Proof. D ⊆ D o ⊆ D is a direct consequence of Proposition 11. Then, from

Lines 1 to 8

 8 correspond to the initialization of these approximations by testing the lower bound p m = (p m 1 , . . . , p m d ) and the upper bound p M = (p M 1 , . . . , p M d ) of the set D. If p m ∈ D o , then by Theorem 12, C ∩ D is a solution to Problem 2. Note that in that case, all timing-contract parameters in C ∩ D guarantee the stability of (1). If p m / ∈ D o , then D \ {p m } is an overapproximation of D o . Similarly, if p M / ∈ D o , then by Theorem 12, D o = ∅. Note that in that case, no timing-contract parameters in C ∩ D can guarantee the stability of (1). If p M ∈ D o , then {p M } is an under-approximation of D o .

Algorithm 2 . 10 :

 210 Timing contract synthesis function TC Synth(Φ p , D, C) input: Φ p , D, C output: C * ⊆ C ∩ D such that for all p ∈ C * , (3) is GES. parameter: ε ∈ R + 1: if p m ∈ D o then 2: return C ∩ D; 3: else D := D \ {p m }; 4: end if 5: if p M / ∈ D o then 6:return ∅; 7: else D := {p M }; 8: end if 9: while d(D, D) > ε do main loop Pick p ∈ D \ D; select next sample 11: if p ∈ D o then D := D ∪ {p ∈ D : p ≤ p }; 12: else D := D \ {p ∈ D : p ≤ p}; 13: end if 14: end while 15: C * := C ∩ D; 16: return C * ;

Figure 1 :

 1 Figure 1: Polytopes S 0 and S 1 computed by Algorithm 1 using parameter setup C for system (20) with T = 0.1 and T = 0.5; S 1 is strictly included in S 0 .

Figure 2 :

 2 Figure 2: Timing contract synthesis for system (20): region T * of timing contract parameters, for which stability is guaranteed.

Figure 3 :

 3 Figure 3: Timing contract synthesis for System (24) in the (T , T ) domain.

Table 2 :

 2 Results of Algorithm 1 on system

Table 3 :

 3 Maximum value of T for which stability of systems (23) and (24) could be proved by our approach and several existing methods, as reported in[START_REF] Briat | Convex conditions for robust stability analysis and stabilization of linear aperiodic impulsive and sampled-data systems under dwell-time constraints[END_REF].

		System (23) System (24)
		T	T	T	T
	[7]	10 -5 1.7279 0.4	1.827
	[11]	10 -5	0.869	-	-
	[19]	10 -5	1.113	-	-
	[10]	10 -5	1.695	-	-
	[18]	10 -5	1.695	-	-
	[20]	10 -5	1.723 0.4	1.251
	[21]	10 -5 1.7294 0.4	1.828
	Algorithm 1 10 -5 1.72941 0.4	1.888

Table 4 :

 4 Results of Algorithm 1 for system (25) with feedback gains K 1 and K 2 .

	T	T	T CP U (s)

Table 5 :

 5 Parameter setup for Algorithm 1 for system (25) with feedback gains K 1 and K 2 .

$ This work was supported by the Agence Nationale de la Recherche (COMPACS project ANR-13-BS03-0004).