Mirko Fiacchini
email: mirko.fiacchini@gipsa-lab.fr

Irinel-Constantin Morȃrescu
email: constantin.morarescu@univ-lorraine.fr

Stability analysis for systems with asynchronous sensors and actuators

The paper provides computation-oriented necessary and sufficient conditions for the global exponential stability of linear systems with asynchronous sensors and actuators. Precisely, we focus on continuous-time linear systems whose state undergoes finite jumps referred to as impulsions. The impulsions are of two types: those related to input updates and those related to measurement updates. We assume that impulsions of each type occur periodically but the periods may be different and the clocks at the sensor and at the actuator are not synchronized. We first show that the analysis can be reduced to a finite time domain. Based on that, we provide necessary and sufficient conditions for the global exponential stability of systems belonging to the class under study. An example illustrates numerically the proposed results.

I. INTRODUCTION

The analysis of impulsive systems is motivated by sampled data 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], [START_REF] Hespanha | Lyapunov conditions for input-to-state stability of impulsive systems[END_REF], electronic circuits with diodes [START_REF] Vasca | A new perspective for modeling power electronics converters: complementarity framework[END_REF], [START_REF] Schumacher | Complementarity systems in optimisation[END_REF] or mechanical systems with impacts and dry friction [START_REF] Facchinei | Finite-dimensional Variational Inequalities and Complementarity ProblemsResearch[END_REF], [START_REF] Filippov | Differential equations with discontinuous righthand sides[END_REF], [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF]. Stability and performances of impulsive systems received a lot of attention as illustrated by the following recent works [START_REF] Beker | Fundamental properties of reset control systems[END_REF], [START_REF] Nesić | Stability properties of reset systems[END_REF], [START_REF] Banos | Reset times-dependent stability of reset control systems[END_REF] and the references therein. Stability and stabilization of linear impulsive systems with quasi-periodic impulsions have been considered in [START_REF] Hetel | Stabilization of linear impulsive systems through a nearly-periodic reset[END_REF]. Furthermore, we proposed in [START_REF] Fiacchini | Constructive necessary and sufficient condition for the stability of quasi-periodic linear impulsive systems[END_REF] a constructive methodology to get necessary and sufficient conditions for the stability of this class of systems.

As pointed out in [START_REF] Wakaiki | Stabilization of systems with asynchronous sensors and controllers[END_REF], in real applications we can encounter situations where the local clocks at the sensor and at the actuator are not synchronized. As mentioned there, clocks synchronization over networks has fundamental limitations [START_REF] Freris | Fundamental limits on synchronizing clocks over network[END_REF] and GPS synchronization may be vulnerable against malicious attacks [START_REF] Jiang | Spoofing gps receiver clock offset of phasor measurement units[END_REF]. In this context [START_REF] Wakaiki | Stabilization of systems with asynchronous sensors and controllers[END_REF] addresses the problem of how large the clock offset between actuator and sensor can be, without compromising the existence of a stabilizing linear time-invariant controller.

In this work we consider a dynamic linear controller that regulates a linear plant. The sensor measurements and the control inputs can be updated periodically but with different sampling period. In this framework, we give computation oriented conditions which are necessary and sufficient for stability of the closed loop. The results provide also a preliminary analysis to the more realistic case in which the sampling periods are not only different, but uncertain within bounded known intervals. The conditions presented in this paper are then necessary for robust stability and the analysis method is extendable, in our opinion, to the uncertain case.

The problem is formulated in terms of global stability analysis of a subclass of hybrid systems studied in [START_REF] Goebel | Hybrid Dynamical Systems[END_REF]. The main contribution of this work is that for this specific subclass we are able to highlight necessary and sufficient conditions characterizing the global uniform exponential stability. Furthermore, our main results provide a numerical methodology to check these necessary and sufficient conditions.

The paper is organized as follows. Section II presents the problem statement and the framework assumptions. Preliminary results allowing to reduce stability analysis to a finite time-domain are provided in Section III. The main results on the constructive necessary and sufficient conditions for the global exponential stability of the system at hand, are developed in Section IV. We illustrate the methodology with one numerical example in Section V before giving some concluding remarks.

II. PROBLEM FORMULATION

In this paper, we consider a continuous-time linear plant

ẋp (t) = A p x p (t) + B p u(t), y(t) = C p x p (t), (1)
where x p (t) ∈ R n p is the state, u(t) ∈ R m the control input and y(t) ∈ R n y the plant output. The plant is in closed loop with a dynamic controller of the form

ẋk (t) = A k x k (t) + B k y(t), u(t) = C k x k (t), (2)
where x k (t) ∈ R n k is the controller's state. We investigate the scenario where the actuator and the sensor are implemented on a digital platform with zero-order-hold devices. Moreover the updates of the control signal and the measured output are done periodically but independently, i.e. we don't assume clocks synchronization and same period for the input and output of the plant. In order to formalize that, we consider T 1 , T 2 ∈ R + the sampling periods of u and y, respectively and we define the following sequences

T 1 = {t i } i∈N : t i+1 = t i + T 1 , t 0 ∈ [0, T 1 ], ∀i ∈ N , T 2 = {t j } j∈N : t j+1 = t j + T 2 , t 0 ∈ [0, T 2 ], ∀ j ∈ N .
Therefore, closed-loop dynamics (1)-( 2) becomes

       ẋp (t) = A p x p (t) + B p v(t), y(t) = C p x p (t), ẋk (t) = A k x k (t) + B k w(t), u(t) = C k x k (t), (3)
where v and w are defined as

v(t) = u(t j ), ∀t ∈ [t j ,t j+1 ), ∀t j ∈ T 1 , w(t) = y(t h ), ∀t ∈ [t h ,t h+1 ), ∀t h ∈ T 2 .
In other words, considering the state vector x(t) = (x p (t), w(t), x k (t), v(t)), the closed-loop dynamics can be formulated as a continuous-time linear system affected by impulses at time instants t j ∈ T 1 ∪ T 2 . The continuous-time dynamics is defined by:

       ẋp (t) = A p x p (t) + B p v(t), ẇ(t) = 0, ẋk (t) = A k x k (t) + B k w(t), v(t) = 0, ∀t ∈ R + \ (T 1 ∪ T 2 ) (4)
while the impulses are given by one of the following equations:

       x p (t j ) = x p (t - j ), w(t j ) = w(t - j ), x k (t j ) = x k (t - j ), v(t j ) = C k x k (t - j ), ∀ j ∈ T 1 (5)
and

       x p (t j ) = x p (t - j ), w(t j ) = C p x p (t - j ), x k (t j ) = x k (t - j ), v(t j ) = v(t - j ), ∀ j ∈ T 2 (6)
where, for any function f (t) we introduce the notation f (t -) = lim τ→0,τ<0 f (t + τ). Summarizing, the overall asynchronously sampled system can be modelled as the following linear impulsive system:

   ẋ(t) = A c x(t), ∀t ∈ R + \ (T 1 ∪ T 2 ), x(t) = A 1 x(t -), ∀t ∈ T 1 , x(t) = A 2 x(t -), ∀t ∈ T 2 , (7)
where x ∈ R n is the system state, i.e. x = (x p , w, x k , v) with n = n p + m + n k + p, and A c , A 1 and A 2 are appropriate matrices defined from ( 4)- [START_REF] Fiacchini | Constructive necessary and sufficient condition for the stability of quasi-periodic linear impulsive systems[END_REF]. Dynamics [START_REF] Filippov | Differential equations with discontinuous righthand sides[END_REF] can be also obtained in the framework of multi-agent systems that update asynchronously their inputs. We will refer to jump (or impulse) of type 1 (resp. 2) if the discrete-time dynamics is defined by the matrix A 1 (resp. A 2 ). Remark 1:

• Throughout this paper we suppose that T 1 ≤ T 2 . This would not entail any loss of generality, since a simple redefinition of the system structure could be done if necessary.

• We don't exclude the case where t ∈ T 1 ∩ T 2 . In this case the state-jump is defined by both A 1 and A 2 . Furthermore, the jump at t ∈ T 1 ∩ T 2 is defined either as x(t) = A 1 A 2 x(t -), or x(t) = A 2 A 1 x(t -) in order to ensure the outer semi-continuity of the jump map (see [START_REF] Goebel | Hybrid Dynamical Systems[END_REF] for details).

A hybrid model of the systems (7) can be directly given by introducing two additional variables τ 1 and τ 2 taking into account the time passed from the last jump of type 1 or 2, respectively, flow and jump maps and sets. That is, denoting

z = (x, τ 1 , τ 2 ) ∈ R n+2 define    ẋ(t) = A c x(t), τ1 (t) = 1, if z ∈ C, τ2 (t) = 1, z + ∈ G(z), if z ∈ D (8)
where

C = R n × [0, T 1 ] × [0, T 2 ], D = R n × {T 1 } × {T 2 } and G(z) = (A 1 x, 0, τ 2 ), if τ 1 = T 1 , (A 2 x, τ 1 , 0), if τ 2 = T 2 . (9)
Notice that, while the continuous-time dynamics is given by a function, G(z) is a set-valued map, as τ 1 = T 1 and τ 2 = T 2 implies that both jumps occurs concomitantly yielding

z + ∈ {(A 1 A 2 x, 0, 0), (A 2 A 1 x, 0, 0)} if τ 1 = T 1 , τ 2 = T 2 .
In this work we deal with the case of commensurate sampling periods. Therefore, throughout the paper, we suppose that the following assumption holds true. Assumption 1:

There exist n 1 , n 2 ∈ N such that T 1 n 1 = T 2 n 2 and n 1 /n 2 is irreducible.
This first assumption imposes the commensurability of T 1 and T 2 and introduces the minimal pair n 1 , n 2 ∈ N defining the ratio T 2 /T 1 . Under this assumption we investigate the global exponential stability of system [START_REF] Filippov | Differential equations with discontinuous righthand sides[END_REF] regardless the initialization of (τ

1 (0), τ 2 (0)) ∈ [0, T 1 ] × [0, T 2 ].
In other words we want to solve the following equivalent problem.

Problem 1: Find necessary and sufficient conditions under which the compact set

A = {0} × [0, T 1 ] × [0, T 2 ], (10)
is globally exponentially stable with respect to the hybrid dynamics [START_REF] Freris | Fundamental limits on synchronizing clocks over network[END_REF].

The main difference with respect to [START_REF] Goebel | Hybrid Dynamical Systems[END_REF] is that we give constructive necessary and sufficient condition that can be numerically checked.

Remark 2: Considering the case of commensurate sampling periods with different clock initializations

(τ 1 (0), τ 2 (0)) ∈ [0, T 1 ] × [0, T 2 ]
is a preliminary analysis for the more general problem in which the sampling periods are not constant but they belong to given bounded intervals. Notice that the general case encloses the one of fixed but not commensurate sampling periods. Moreover, the problem studied here leads to conditions that result to be necessary for the non-constant periods case.

III. ANALYSIS OF THE TIME DOMAIN

In this section we show that the analysis can be reduced to some finite time domains included in [0, T 1 ) × [0, T 2 ). First, notice that, since we give global asymptotic stability results, the analysis can always be reduced to the case in which the first jump does not occur at t = 0. This is equivalent to consider that the following assumption holds.

Assumption 2: Suppose that

τ 1 (0) ∈ [0, T 1 ) and τ 2 (0) ∈ [0, T 2 ).
Throughout the rest of the paper we denote τ = (τ 1 (0), τ 2 (0)).

A. Geometrical illustration of the jumps sequence

Definition 1: Consider the system ( 8)-( 9) with τ

1 (0) = τ1 ∈ [0, T 1 ), τ 2 (0) = τ2 ∈ [0, T 2 ). Let T 1 , T 2 satisfy Assump- tion 1 and D = T 1 n 1 = T 2 n 2 , n D = n 1 + n 2 , T ∈ [0, D]. We define T 1,T ( τ) = {t ∈ R + : (t -τ1 ) mod T 1 = 0, ∧ (t -τ1 ) ≤ T } ∈ [0, T ]}, T 2,T ( τ) = {t ∈ R + : (t -τ2 ) mod T 2 = 0, ∧ (t -τ2 ) ≤ T } ∈ [0, T ]}, n 1,T ( τ) = |T 1,T ( τ)|, n 2,T ( τ) = |T 2,T ( τ)|, n T ( τ) = n 1,T ( τ) + n 2,T ( τ), T T ( τ) = {{t k } k∈N n T : t k ≤ t k+1 , ∀k ∈ N n T -1 ∧ t k ∈ T 1,T ∪ T 2,T , ∀k ∈ N n T }, T ( τ) = lim T →+∞ T T ( τ). ( 11
)
Roughly speaking: T i,T ( τ) is the sequence of instants of jumps of type i in the interval [0, T ] and n i,T ( τ) its length, with i ∈ N 2 ; T T ( τ) the overall sequence of jump instants. Notice that two jump instants can coincide, i.e.t k could be equal with t k+1 .

Let us illustrate the set T D ( τ) (i.e. T T ( τ) in ( 11) when T = D) with the following example.

Example 1: Consider the system ( 8) with ( 9) and focus on the dynamics of the variables τ 1 and τ 2 , that are not affected by the value of x. Let us suppose for this illustration that τ 1 (0) = τ 2 (0) = 0. We consider that Assumption 1 holds with T 1 /T 2 = n 2 /n 1 = 7/3, say T 1 = 7a and T 2 = 3a with a > 0. Thus, as can be noticed by graphical inspection in Figure 1, left, the sequence of the first 8 jumps whose time is denoted as t i with i ∈ N 8 , is {(t 1 , 2), (t 2 , 2), (t 3 , 1), (t 4 , 2), (t 5 , 2), (t 6 , 1), (t 7 , 2), (t 8 , 2)} [START_REF] Jiang | Spoofing gps receiver clock offset of phasor measurement units[END_REF] while at time t 9 both the first and the second jump occur, that means that the subsequence (12) might be followed either by {(t 9 , 1), (t 10 , 2)} or {(t 9 , 2), (t 10 , 1)} with t 9 = t 10 . The values of t i are {t i } i∈N 10 = {3a, 6a, 7a, 9a, 12a, 14a, 15a, 18a, 21a, 21a} [START_REF] Khatib | Verification and synthesis of timing contracts for embedded controllers[END_REF] and

τ 1 (21a) = τ 2 (21a) = 0 = τ 1 (0) = τ 2 (0)
. This yields that, at the time instant which is the minimal common multiplier of n 1 and n 2 , the partial state (τ 1 , τ 2 ) recovers its initial value. Notice that t h 1 +1 = t h 1 + 7a and t h 2 +1 = t h 2 + 3a as required, for all h 1 ∈ T 1,D and h 2 ∈ T 2,D . Thus, the jump sequence is given by the concatenation of the subsequence [START_REF] Jiang | Spoofing gps receiver clock offset of phasor measurement units[END_REF] and the jump instants are

T = {t i } i∈N = {3a + 21ak, 6a + 21ak, 7a + 21ak, 9a + 21ak, 12a + 21ak, 14a + 21ak, 15a + 21ak, 18a + 21ak, 21a + 21ak, 21a + 21ak} k∈N .
An analogous graphical illustration of the sequence of jumps is illustrated in Figure 1, right. Since Assumption 1 is satisfied, the set in the space of (τ 1 , τ 2 ) ∈ R 2 can be constructed by replicating the basic set [0, 7a] × [0, 3a] to obtain a square in R 2 , in particular [0, 21a] × [0, 21a]. Notice that the jump instants are given by the points in (τ 1 , τ 2 ) ∈ R 2 such that τ 1 = 7ak and τ 2 = 3ak for all k ∈ N. Moreover, analogous considerations hold for the case of (τ 1 (0), τ 2 (0)) = (0, 0).

B. Finite time-domain reduction

The sets and sequences in Definition 1 permit to determine a cyclic behavior of the hybrid system (8) with initial time variable τ.

Proposition 1: Given the system ( 8) and ( 9) satisfying Assumption 1, the sequence of jumps is Corollary 1: Given the system (8) and ( 9) satisfying Assumption 1, for every τ ∈ [0, T 1 ) × [0, T 2 ), the set T D ( τ) contains at most one element t k such that t k = t k+1 .

T ( τ) = {T D ( τ), T D ( τ), . . .}. ( 14
Proof: Suppose by contradiction that, there are more than one element in T D ( τ) such that t k = t k+1 and consider two of them. Denote them t j and t i with j < i, and n1 and n2 the jumps of type T 1 and T 2 respectively, that occur in (t j ,t i ]. By hypothesis, n1 ≥ 1 and n2 ≥ 1. Since t jt i must be smaller than D, then n1 < n 1 and n2 < n 2 . Thus it follows that t it j = T 1 n1 = T 2 n2 with n1 < n 1 and n2 < n 2 , that contradicts Assumption 1.

Consider first the case in which the simultaneuos jump occurs at one jump instant, denote it t j with j ∈ N n D . The jump map takes the set value {(A 1 x, 0, τ 2 )∪(A 2 x, τ 1 , 0)} and the following jump instant is t j+1 = t j . Thus the discrete-time dynamics of x is given by the following difference inclusion

x(k + 1) ∈ A D ( τ)x(k), (15)
if t(0) = τ, where we denoted x(kD) as x(k), with slight abuse of notation, and

A D ( τ) = A (12) D ( τ) ∪ A (21) D ( τ) with A (12) D ( τ) e A c δ n D +1 n D ∏ k∈ j+2 A q k e A c δ k A 2 A 1 e A c δ j j-1 ∏ k∈1 A q k e A c δ k A (21) D ( τ) e A c δ n D +1 n D ∏ k∈1 A q k e A c δ k A 1 A 2 e A c δ j j-1 ∏ k∈1 A q k e A c δ k (16)
where δ k = t kt k-1 , with τ 0 ( τ) = 0 and δ n D +1 ( τ) = min{ τ1 , τ2 } and q k ∈ {1, 2}. In the case when simultaneous jumps do not occur before D the jump map [START_REF] Goebel | Hybrid Dynamical Systems[END_REF] takes values that are singletons. Then the dynamics of the state x at the time instants kD with k ∈ N is [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF] with

A D ( τ) = A (0) D ( τ) e c δ n D +1 n D ∏ k∈1 A q k e A c δ k . ( 17
)
Indeed, the dynamics reduces in this case to a difference equation rather than an inclusion. Then, summarizing, the discrete-time dynamics are given by [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF] with

A D = A (12) D ∪ A (21) D if ∃ j ∈ N n D -1 : t j = t j+1 , A (0) D , otherwise, (18)
where the dependence on τ has been avoided for notational simplicty.

Of course, A must be Schur for all A ∈ A D ( τ) for the discrete-time system [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF] to stay bounded and converge to the origin. Notice that such condition is only necessary, the necessary and sufficient one will be discussed in the sequel. The fact that A D ( τ) contains Schur or not Schur matrices might depend on the initial value τ, as illustrated in the example below.

Example 2: Consider the system (7) with matrices

A c = 0 -2π 2π 0 , A 1 = a -ε 0 0 a -1 , A 2 = a -1 0 0 a -ε with a > 1, ε > 0 such that a -ε > 1 and T 1 = T 2 = 0.5.
Recall that e A c t = cos(2πt)sin(2πt) sin(2πt) cos(2πt) .

For initial condition τ = (0.25, 0) one has

A D ( τ) = A 2 e A c 0.25 A 1 e A c 0.25 = -a -2 0 0 -(a -ε) 2 ,
which is not Schur from a-ε > 1. On the other hand, setting τ = (0, 0) and noticing that A 1 A 2 = A 2 A 1 , one gets

A D ( τ) = A 2 A 1 e A c 0.5 = -a -1 (a -ε) 0 0 -a -1 (a -ε) ,
whose eigenvalues are a -1 (aε) < a -1 a = 1 and then the matrix is Schur. An interesting issue is whether all A ∈ A D ( τ) are Schur for every possible initial condition τ ∈ [0, T 1 ) × [0, T 2 ). The property below is used in the following.

Proposition 2: Given the system ( 8) and ( 9) satisfying Assumptions 1 and 2, the spectrum of the matrices A (0)

D ( τ + d1), A (12)
D ( τ + d1) and A

D ( τ) = A (0) D-d ( τ + d1)A d ( τ). More- over A (0) D ( τ + d1) = A (0) d ( τ + D1)A (0) D-d ( τ + d1) and, since A (0) d ( τ + D1) = A (0) d ( τ), then σ A (0) D ( τ) = σ A (0) D-d ( τ + d1) A (0) d ( τ) = σ A (0) d ( τ) A (0) D-d ( τ + d1) = σ A
D ( τ) = σ A (0) D ( τ + kD1) = σ A (0) D ( τ + kD1 + d) = σ A (0) D ( τ + d1) .
Geometrically, it means that analyzing the stability property of A D ( τ) implies analyzing the same property for every (τ 1 , τ 2 ) along the line starting from τ.

From Proposition 2, the analysis can be limited to the points, for instance, such that either τ 1 = 0 and τ 2 ≥ 0 or τ 2 = 0 and τ 1 ≥ 0.

Proposition 3: Given the system (8) and ( 9) satisfying Assumption 1, the matrices A with A ∈ A D ( τ) are Schur for every τ ∈ [0, T 1 ) × [0, T 2 ), if and only if they are Schur for every τ ∈ [0, T 1 ) × 0 or for every τ ∈ 0 × [0, T 2 ).

Proof: Necessity is clear. For sufficiency, suppose that A is Schur for every A ∈ A D ( τ) for all τ ∈ [0, T 1 ) × 0, the case τ ∈ 0×[0, T 2 ) is analogous. For every τ ∈ [0, T 1 )×[0, T 2 ) there exists d ∈ (0, T 2 ] such the jump of type 2 occurs at time

d. Then τ 2 + d = kT 2 and τ + 1d ∈ [0, T 1 ) × 0 with k ∈ N.
Thus, from Proposition 2, since the matrices of A D ( τ) and A D ( τ + 1d) have the same spectrum and the latter are Schur by assumption then every A ∈ A D ( τ) is Schur.

C. Further reduction of the time-domain

Hereafter it will be proved that the stability analysis can be further reduced to some subset of [0, T 1 ) × 0.

Lemma 1: Given T 1 and T 2 such that Assumption 1 holds, for every initial condition τ(0) = τ of (8) with τ ∈ [0, T 1 )×0 and every a ∈ N it follows that

   τ(aT 2 ) = τ if a = mn 2 , |τ(aT 2 ) -τ| ≥ T 1 n 2 otherwise. ( 19
)
with m ∈ N.

Proof: First notice that τ 2 (aT 2 ) = τ 2 (0) = τ2 for all a ∈ N then it is sufficent to consider |τ 1 (aT 2 ) -τ1 |, in particular the condition on a such that

|τ 1 (aT 2 ) -τ1 | < T 1 n 2 , or equivalently τ1 - T 1 n 2 < τ 1 (aT 2 ) < τ1 + T 1 n 2 . ( 20
)
Since τ 1 (aT 2 ) = ( τ1 + aT 2 )mod (T 1 ) and

T 2 = T 1 n 1 /n 2 , from Assumption 1, then (20) holds only if there exists b ∈ N such that τ1 - T 1 n 2 + bT 1 < τ1 + a n 1 n 2 T 1 < τ1 + T 1 n 2 + bT 1 ⇔ -1 < an 1 -n 2 b < +1. (21)
Being

an 1 -n 2 b an integer, then condition (21) is satisfied if and only if an 1 -n 2 b = 0, that is if a = n 2 (b/n 2 )
. Then (20) cannot hold if there is not m = (b/n 2 ) ∈ N such that a = mn 2 , which proves the second part in [START_REF] Wakaiki | Stabilization of systems with asynchronous sensors and controllers[END_REF]. Moreover, if a = mn 2 with m ∈ N then τ 1 (aT 2 ) = τ 1 (mn 2 T 2 ) = τ 1 (mn 1 T 1 ) and then τ 1 (aT 2 ) = τ1 , which terminates the proof. Geometrically, the Lemma 1 means that, if the trajectory of τ is considered in the extended space, as in Figure 2, right, then at the time of the first n 2 -1 jump of type 2, τ 1 is not in the open neighborhood of witdh T 1 /n 2 around τ1 . Moreover, at the n 2 -th jump of type 2, one has t(n 2 T 2 ) = τ.

Remark 3: From Lemma 1 it follows that the only initial condition τ ∈ [0, T 1 /n 2 ) × 0 for which simultaneous jumps occur is τ = (0, 0). Proposition 4: Given T 1 and T 2 such that Assumption 1 holds, for every initial condition τ(0) = τ of [START_REF] Freris | Fundamental limits on synchronizing clocks over network[END_REF]

with τ ∈ T 1,2 (k) = [(k -1)T 1 /n 2 , kT 1 /n 2 )×0, for k ∈ N n 2 , one has that τ(aT 2 ) /
∈ T 1,2 for all a ∈ N n 2 -1 and τ(n 2 T 2 ) = τ. Moreover, the jump sequence associated with any τ(0) ∈ T 1,2 (k) is the same.

The matrix A D ( τ), as in Definition 1, can be written in the following form

A D ( τ) =    A (12) D ( τ) ∪ A (21) D ( τ) if τ = (0, 0), ∏ k∈N m R k e A k δ , otherwise, (22)
for all τ = (δ , 0) ∈ [0, T 1 /n 2 ) × 0, with m ≤ n D , for appropriate matrices R k , A k ∈ R n×n for k ∈ N m . The following example shows how to build the matrix representation of A D ( τ) as in ( 22).

Example 3: Consider the system ( 8) and ( 9) satisfying Assumption 1, with T 1 = 3 and T 2 = 7 and τ = (δ , 0) ∈ [0, 1) × 0. Then, n 1 = 7, n 2 = 3 and D = 21. The graphical representation of the trajectory of the states τ 1 and τ 2 is depicted in Figure 2, right. From geometrical inspection we have that, for δ = 0, one has

A (12) D ( τ) = A 2 A 1 e A c 3 A 1 e A c 3 A 1 e A c A 2 e A c 2 •A 1 e A c 3 A 1 e A c 2 A 2 e A c A 1 e A c 3 A 1 e A c 3 , A (21) D ( τ) = A 1 A 2 e A c 3 A 1 e A c 3 A 1 e A c A 2 e A c 2 •A 1 e A c 3 A 1 e A c 2 A 2 e A c A 1 e A c 3 A 1 e A c 3 , (23)
while

A (0) D ( τ) = A 2 e A c δ A 1 e A c 3 A 1 e A c 3 A 1 e A c e -A c δ A 2 e A c 2 e A c δ •A 1 e A c 3 A 1 e A c 2 e -A c δ A 2 e A c e A c δ A 1 e A c 3 A 1 e A c 3 e -A c δ ,
for δ ∈ (0, 1), that can be rewritten as in ( 22) with m = 6 and

A 1 = -A c , A 2 = A c , A 3 = -A c , A 4 = A c , A 5 = -A c , A 6 = A c , R 1 = A 1 e A c 3 A 1 e A c 3 , R 2 = A 2 e A c , R 3 = A 1 e A c 3 A 1 e A c 2 , R 4 = A 2 e A c 2 , R 5 = A 1 e A c 3 A 1 e A c 3 A 1 e A c , R 6 = A 2 .

IV. STABILITY ANALYSIS: NECESSARY AND SUFFICIENT CONDITIONS

Before providing the main result of the paper, concerning the constructive necessary and sufficient conditions for global exponential stability, a property is given that is functional for the objective.

Definition 2: The convex compact set Ω ⊆ R n with 0 ∈ int Ω is a λ -contractive set, with λ ∈ [0, 1), for the difference inclusion x + ∈ Āx with Ā ⊆ R n×n if AΩ ⊆ λ Ω for all A ∈ Ā.

The following result comes from the literature, [START_REF] Molchanov | Criteria of asymptotic stability of differential and difference inclusions encounterd in control theory[END_REF], [START_REF] Blanchini | Nonquadratic Lyapunov functions for robust control[END_REF]. Lemma 2: The difference inclusion x + ∈ Āx is globally exponentially stable if and only if there exists a polytopic λ -contractive set.

The set of symmetric positive definite matrices in R n×n is denoted S n .

Theorem 1: The compact set A defined in [START_REF] Hespanha | Lyapunov conditions for input-to-state stability of impulsive systems[END_REF] is globally exponentially stable for the system (8) with (9) satisfying Assumption 1 with decay rate λ ∈ [0, 1), if and only if there exits P

i : (0, T 1 /n 2 ) → S n for i ∈ N m+1 such that e A T i δ R T i P i+1 (δ )R i e A i δ ≤ P i (δ ), ∀i ∈ N m , (24)
and P 1 (δ ) ≤ λ 2 P m+1 (δ ), with τ = (δ , 0), for all δ ∈ (0, T 1 /n 2 ), with R k , A k , for k ∈ N m as in ( 22) and a λcontractive polytope for Ā = A

D ((0, 0)) ∪ A

D ((0, 0)). Proof: Due to space limitations we do not provide the proof here.

The condition given in Theorem 1 is necessary and sufficent for stability of system (8) and ( 9), for every initialization of (τ 1 , τ 2 ) ∈ [0, T 1 ) × [0, T 2 ), but it concerns an infinite number of LMI conditions. Since our aim is to provide a constructive method, a computation-oriented equivalent result is given hereafter.

Proposition 5: The compact set A defined in ( 10) is globally exponentially stable for the system (8) with [START_REF] Goebel | Hybrid Dynamical Systems[END_REF] satisfying Assumption 1 with decay rate λ ∈ [0, 1), if and only if for all ν ∈ (λ , 1) there exist: µ ∈ R; J ∈ N; P i, j ∈ S n and α i, j ∈ R for i ∈ N m+1 and j ∈ N J such that

e A T i δ j R T i P i+1, j R i e A i δ j ≤ µ 2 P i, j , ∀i ∈ N m , j ∈ N J , A T i P i, j + P i, j A i ≤ α i, j P i, j , ∀i ∈ N m , j ∈ N J , max{1, e √ α i, j ∆ } ≤ µ -1 , ∀i ∈ N m , j ∈ N J , δ j+1 = δ j + ∆, ∀ j ∈ N J , (25)
with ∆ = T 1 /(n 2 J); and P 1, j ≤ ν 2 P m+1, j , with R k , A k , for k ∈ N m as in ( 22) and a λ -contractive polytope for Ā = A

D ((0, 0)) ∪ A

D ((0, 0)). Proof: The proof consists in demonstrating that (24) and (25) are equivalent. An analogous result has been proved in [START_REF] Fiacchini | Constructive necessary and sufficient condition for the stability of quasi-periodic linear impulsive systems[END_REF].

V. NUMERICAL EXAMPLE

Although a direct comparison is not possible since the problems are not the same, we inspire our example from the literature, in particular 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], [START_REF] Khatib | Verification and synthesis of timing contracts for embedded controllers[END_REF]. The plant dynamics is given by

A p = 0.1 1 0 -0.1 , B p = 0 0.1 , C p = 1 0 ,
which is analogous to an example treated 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], [START_REF] Khatib | Verification and synthesis of timing contracts for embedded controllers[END_REF] but we replaced the pole in 0 with one in 0.1 and the state is supposed to be not accessible.

A k = A p + B p K + LC p , B k = -L, C k = K,
where L is the observer gain designed such that A p + LC p is Hurwitz. The overall system is exponentially stable in closed loop. The matrices A c , A 1 and A 2 in (7) are obtained by using the numerical values given above in equations ( 4)- [START_REF] Fiacchini | Constructive necessary and sufficient condition for the stability of quasi-periodic linear impulsive systems[END_REF]. We consider the asynchronous sampling times T 1 = 3a and T 2 = 7a and try to maximize a such that the resulting sampled-data system is exponentially stable for every possible initialization of the clocks. We found that one among the matrices A D ( τ) has a spectral radius of 0.9989 for a = 0.254 and of 1.0076 for a = 0.255. Applying our method we found that condition (25) with a = 0.254 is satisfied with J = 4000 and the maximal contraction obtained is 0.9999. Moreover, since A 1 and A 2 commute, than A and then testing that one of them is Schur, which is the case, is sufficient. Finally we run a simulation with a = 0.254 to test the convergence of the system, see Figure 3.

VI. CONCLUSIONS

In this paper we presented a constructive method to check whether a linear closed-loop system with asynchronous input and output sampling is globally exponentially stable. The approach is based on a computation-oriented necessary and sufficient condition leading to a finite set of convex constraints. The result presented is also a preliminary step for the analysis of systems with uncertain sampling periods.

Fig. 1 .

 1 Fig. 1. Graphical construction of the jump sequence.

  ) Proof: The result follows directly from the construction of the structures defined in[START_REF] Hetel | Stabilization of linear impulsive systems through a nearly-periodic reset[END_REF]. First notice that T 1,D ( τ) and T 2,D ( τ) are the time instants in the interval (0, D] at which the first and the second type impulses occur, respectively. Thus, T D ( τ) is the sequence of jump time instants in [0, D]. Moreover, from Assumption 1, the value of τ 1 , and τ 2 , after a time interval of length D are τ1 and τ2 , respectively, and n 2 jumps of type 1 and n 1 jumps of type 2 have occurred, see Figure2, left. Therefore, the subsystem whose state is (τ 1 , τ 2 ) returns to the initial value after n 2 jumps 1 and n 1 jumps 2 after D, and then the same happens after 2D, and kD with k ∈ N in general.

Fig. 2 .

 2 Fig. 2. Graphical construction of the jump sequence

D

  ( τ + d1) does not depend on d ∈ R.Proof: We consider in the proof the matrix A (0) D ( τ), analogous reasonings hold for the matrices A

D

  ( τ) as in[START_REF] Nesić | Stability properties of reset systems[END_REF]. Suppose first that d ∈ [0, D]. By construction we have that A (0)

D

  ( τ + d1) , where the second equality follows from the fact that σ (AB) = σ (BA) for every A, B ∈ R n×n . Finally, if d / ∈ [0, D], there exists k ∈ Z such that kD ≤ d ≤ (k + 1)D which means that d = (d -kD) ∈ [0, D] and then, applying iteratively the result above, one has σ A (0)

Fig. 3 .

 3 Fig. 3. Time evolution of the system ouput with T 1 = 3 • 0.254 = 0.7620 and T 2 = 7 • 0.254 = 1.7780

  The dynamical output-

	8
			2.6
			2.4
	6
			2.2
			2
	4
			1.8
			1.6
	2		1.4
			97	98	99	100	101	102	103	104	105	106
	0
	-2
	-4
	0	100	200		300				400	500

This work was partially funded by the ANR project COMPACS -"Computation Aware Control Systems", ANR-13-BS03-004.


Flora Lucassou

Florian Koch

What uses of piezometric network's data in hard-rock aquifers ? Example of Brittany's network 10 years after its implementation.

What uses of piezometric network's data in hard-rock aquifers ? Example of Brittany's network 10 years after its implementation.

Mougin B.

(1)* , Lucassou F. (1) , Koch F. (1) (1) b.mougin@brgm.fr, BRGM Direction régionale Bretagne, Rennes, France For a long time ago considered as a region devoid of big aquifers, the Brittany region (in France) has a piezometric network since 10 years contrary to other regions equipped since the 1970s. Among the 52 monitoring stations, 48 are settled in hard-rock aquifers (4 others follow alluviums or tertiary aquifers). 10 years after its implementation, piezometric network's basic data of Brittany are useful for various valorisations and plenty of different uses.

References (in French)

HAPPY 10th BIRTHDAY

Brittany's regional piezometric network « 0400000020 -RRESOUPBRE » in French National Data Base on Groundwater ADES

  Agence de l'eau Loire-Bretagne (décembre 2013) -État des lieux du bassin Loire-Bretagne établi en application de la Directive cadre européenne sur l'eau Lucassou F., Mougin B. avec la collaboration de Stollsteiner P. (2015a) -Essai d'élaboration d'indicateurs piézométriques pour la gestion quantitative AEP dans le département des Côtes d'Armor. Rapport final BRGM/RP-64123-FR Lucassou F. (2015b) -Avis hydrogéologique suite aux inondations par remontée de nappe survenues le 12 février 2014 sur la commune de Tréogat (Finistère), dans le cadre d'une demande de reconnaissance de l'état de catastrophe naturelle. Rapport d'expertise CATNAT. Rapport BRGM/RP-65017-FR Mougin B., collaboration Jégou J-P. et Marmu B. (2012) -SILURES Suivi -Gestion du réseau de surveillance piézométrique régional de Bretagne -Valorisation des données des 52 piézomètres et annuaire 2011 -Rapport final -BRGM/RP-60673-FR Observatoire de l'Eau en Bretagne (Janvier 2012) -Tableau de bord interactif Eau souterraine. Variation du niveau des nappes Vernoux J.F. (2015) -Réseau de référence piézométrique pour le suivi de l'impact du changement climatique sur les eaux souterraines. Rapport final BRGM/RP-64858-FR Abstract n°1565 (Session 8.03 Hard-Rock aquifers)

  Knowledge of hard-rock aquifers : dynamics and slowness, cycles of data variation, memory effect, overflowing, time of water table's drop, sea-level, ground signals, and climate change observation(Vernoux et al., 2015)

	Brittany's piezometric network is constituted by 52 boreholes
	(average depth 35 m) equipped with automatic acquisition of
	groundwater level (1 measure/hour) and with data's remote transmission. Their location has been carefully cogitated : it allows to cover main geological and hydrogeological formations existing in Brittany and takes into account weather	Evaluation of groundwater bodies quantitative status (Agence de l'eau Loire-Bretagne, 2013)
	conditions which are very contrasted from west to east and
	between coasts and the inland. After check and validation,
	data are stored (1 measure/day) and put at the public's disposal twice a month on the French National groundwater Data Base called ADES (Mougin et al., 2012).	Bulletins of groundwater's level situations (regional, basin, national)

Elaboration of piezometric indicators dedicated to the departmental management of drinking water (Lucassou et al., 2015a) Contribution of groundwater to rivers flows Diagnosis of floods by water table rise (Lucassou et al., 2015b) Interactive dashboard (Observatoire de l'Eau en Bretagne, 2012)


Masoud Najafi
email: masoud.najafi@inria.fr

Damiano Dipenta
email: damiano.di-penta@renault.com

Karim Bencherif
email: karim.ben-cherif@renault.com

Michel Sorine
email: michel.sorine@inria.fr

D Di Penta

Modeling and simulation of a reformate supplied PEM fuel cell stack, application to fault detection

I. INTRODUCTION

Proton exchange membrane (PEM) fuel cells are the main type of fuel cell developed for ground vehicle applications.

Such applications require a fuel cell control, fault detection, and diagnosis. The three major control subsystem loops in a fuel cell system are applied to the air/fuel supply regulation, the water management, and the heat management [START_REF] Pukrushpan | Modeling and control for PEM fuel cell stack system[END_REF].

In this paper we concentrate on the modeling and simulation of the thermal effects and the temperature distribution in a fuel cell stack.

Existing stack thermal models as in [START_REF] Dipenta | A Reduced Fuel Cell Stack Model for Control and Fault Diagnosis[END_REF] consider the stack as a lumped mass and do not describe thermal distributions in the stack. Typically, these models use a single average cell temperature or stack coolant outlet temperature, as an indicator of the stack thermal condition. However, in a stack, each cell have its own electrical and thermal characteristics different from its neighbors. Furthermore, existing models do not explicitly account for the effect of the end-cells of the stack. The thermal mass of the end-cells draws heat from the cells at both ends of the stack and affects the stack temperature distribution. And above all, the model obtained for the stack is often a large stiff differential equation that is difficult or is not appropriate for real-time applications such as fault detection and process control.

One of the main difficulties in the fuel cell control system is the need for an on-board control and fault detection device to supervise the operation of the system through reliable strategies. These strategies are based on real-time analysis of on-board sensor signals by checking their coherence with mathematical models. The model of each cell includes electrochemical, electrical, and thermal phenomena. Consequently, the large range of transients involved makes the model stiff. The stack model is, then, a large stiff DAE that is too slow for real-time applications. In this paper, we develop a model for fuel cell stack, then we apply the orthogonal collocation method to reduce the model and simulate it with only a few cells. Using the reduced model of stack will be well suited for controllers and diagnosis algorithms.

The one-dimensional, fuel cell stack thermal model developed in this paper considers the features which form an analysis tool for thermal characteristic and fault detection of a PEM fuel cell stack. The fuel cell stack model has been developed with the Modelica language [START_REF] Fritzson | Principles of Object-Oriented Modeling and Simulation with Modelica 2[END_REF]. The Modelica models are then simulated in the Scilab/Scicos (www.scilab.org) environment. The developed models of the fuel cell stack consists of N (50, 100) cell, electrically in series, fed with hydrogen-rich reformate and compressed air.

II. FUEL CELL MODEL

In this section, the model of a PEM fuel cell stack is described. Fig. anode of the PEMFC stack whereas the cathode is fed by compressed air. The reformer breaks down hydrocarbons or hydrogen rich fuels into hydrogen, carbon dioxide, carbon monoxide and water. One of the main problems of reforming hydrocarbons is the production of carbon monoxide which is a poison for the stack but this effect can be mitigated by introducing an air bleed into the anode (see Fig. 1). The model which includes electro-chemical, electrical and thermal phenomena is used to build a model for the fuel cell stack. The internal structure of the fuel cell will not be explained here, interested users are referred to [START_REF] Yi | An Along-the-Channel Model for Proton Exchange Membrane Fuel Cells[END_REF][START_REF] Mench | An Introduction to Fuel Cells and Related Transport Phenomena[END_REF][START_REF] Dipenta | A Reduced Fuel Cell Stack Model for Control and Fault Diagnosis[END_REF].

Electrodes of each cell in the stack are supplied by a mass fluid flow q in k at the temperature T in k , the pressure P in k and the mass composition Y in i,a , where i ∈ {h, co, n, v, air}, and Y in i,c , where i ∈ {o, n, v}. We suppose that there is no liquid water entering the cell electrodes. The electrochemical reaction produces a current noted I and voltage V cell .

The solid temperature is noted T s and lastly the cell is cooled by a coolant flow q in cool at the temperature T in cool . We define The input variable vector U in k and output variables vector U out k of the fuel cell channel k ∈ {a, c, cool} are illustrated in Fig. 2, with

U in a = T in a q in a P in a λ in a Y in h Y in co Y in air Y in va U out a = [T out a P out a Y out h Y out wa ] U in c = T in c q in c P in c Y in o Y in vc λ in c U out c = [T out c P out c Y out o Y out wc ] U in cool = q in cool T in cool U out cool = T out
← U in a ← U in cool ← U in c V cell

Ts

Electrical axis

I θ co

M embrane

Anode Cathode E l e c t r i c a l a x i s

U out c ← U out cool ← U out a ←
FIG. 2: Fuel cell stack and input and output variables in the fuel cell model

Mass balance equations

Mass balance equations are obtained from mass conservation laws. We apply the mass conservation law on hydrogen, oxygen and water in electrodes. Since the procedure to get these equations are identical, we perform it on hydrogen species, the others are found in a similar way. Let m h be the hydrogen mass in the anode. Thus, we have

dm h dt = q in h -q out h -M h I 2F = q in a Y in h -q out a Y out h -M h I 2F
where M h I 2F indicates the conversion of H 2 to H + ions and is obtained from the Faraday law. We suppose that the hydrogen mass in anode m h is equal to Y out h ρ a V a . Assuming a dilute concentration, i.e., q out a ≈ q in a and ρ a slightly constant, we get the differential equation ( 1) for the conservation of hydrogen mass in the anode. In a similar way, the oxygen mass balance equation, i.e., (2), the water mass balance equations in electrodes, i.e., (3) are obtained.

dY out h dt = q in a ρ a V a (Y in h -Y out h ) - M h I 2ρ a F V a ( 1
)
dY out o dt = q in c ρ c V c (Y in o -Y out o ) - M o I 4ρ c F V c ( 2
)
dY out w,k dt = q in k ρ k V k (Y in v,k -Y out w,k ) + α k M w I 2ρ k F V k (3)
The value of the input mass flows q in k k ∈ {a, c}, are determined as a function of the current I, of the inlet mass fraction of hydrogen and oxygen, noted Y in h and Y in o , and of the stoichiometries in the electrodes, noted λ a and λ c , as given in ( 4) and [START_REF] Mench | An Introduction to Fuel Cells and Related Transport Phenomena[END_REF].

q in a = λ a M h I 2F Y in h ( 4
)
q in c = λ c M o I 4F Y in o ( 5
)
where M h and M o are the molar mass of the H 2 and O 2 respectively.

As we will see in the next section, the cell voltage is a function of the amount of reactants inside the electrodes.

Thus we define the mean partial pressures of H 2 at the anode and of O 2 at the cathode by

P h = M a 2M h (Y in h P in a + Y out h P out a ) (6)
P o = M c 2M o (Y in o P in c + Y out o P out c ) (7)
where the average molar mass in the cathode and the anode, neglecting the mass fraction of CO and air bleed, will be:

M a = Y in h M h + Y in v,c M w + Y in n,a M n,a -1 M c = Y in o M o + Y in v,c M w + Y in n,c M n,c -1
where the inlet mass fractions of neutral gas are given by:

Y in n,a = 1 -Y in h -Y in w,a (8)
Y in n,c = 1 -Y in o -Y in w,c (9)
where the inlet pressures are obtained through a simple pressure drop laws:

P in k = P out k + K dp,k q in k P out k
Finally, we suppose that all the CO and the air bleed is consumed in the the anode:

P co = M a 2M co Y in co P in a ( 10
)
P air = M a 2M air Y in air P in a ( 11
)

Fuel cell voltage

The electrical voltage of a cell can by modeled by the difference between the two electrode potentials and overpotentials, i.e.,

V cell = E r -η a -η c -η ohm (12)
We assume that water is produced in liquid phase, so E r , the reversible potential, corresponding to the liquid water formation is:

E r = 1 2F ∆H 0 l -∆S 0 l T s = E 0 - ∆S 0 l 2F T s + RT s 2F ln P o P h
where E 0 , ∆H 0 l , and ∆S 0 l are the thermo-neutral potential, the variations in standard enthalpy, and standard entropy of liquid water formation, respectively.

The ohmic loss is the voltage loss due to the internal resistance R ohm of the stack cells. Supposing that the membranes are optimally hydrated, we can consider that these resistances are constant, thus the ohmic losses can be computed as given:

η ohm = R ohm J ( 13
)
where J is the current density defined by I = JA act where A act is the active area of a cell.

The anode overpotential η a is mostly linked to CO poisoning is expression was obtain by Di Penta et al. [START_REF] Dipenta | Fuel Cell Carbon Monoxide Poisoning Estimation and Control with Air Bleed Injection[END_REF] :

η a = 2RT s F sinh -1 J 2(k eh θ h + k ec θ co + k res ) ( 14
)
where T s , R, and F are the cell solid temperature, the gas constant, and the Faraday constant. σ is the density of sites on the catalyst layer. The coverage fractions of carbon monoxide and hydrogen carbon monoxide on the anode catalyst sites noted θ co and θ h are given by

σ dθ co dt = k f c P co - k f o k oc θ co P air k oc θ co + k oh θ h + b f o k f o (1 -θ h -θ co ) -b f c k f c θ co - k ec θ co J k eh θ h + k ec θ co + k res (15) θ h = k f h P h k f h P h + b f h k f h (1 -θ co ) - J k f h P h + b f h k f h (16) with b f c = b 0 f c exp -E bf c RT s
where P co , P air , P h , are the mean partial pressures of CO, air bleed, and H 2 at the anode which were defined ( 6), ( 10) and [START_REF] Pereira | Modeling of Dynamic Systems[END_REF].

Contrary to the anode, the cathode is not a place of competing reactions. Nevertheless, the cathode is at the origin of most of the activation and diffusion losses. The cathode overpotential is given by:

η c = 2RT s F sinh -1 k dif f P o J 2k act (k dif f P o -J) (17)
where

k act = k 0 act exp -E kact RT s , k dif f = k 0 dif f exp -E kdif f RT s
where P o is defined in [START_REF] Dipenta | Fuel Cell Carbon Monoxide Poisoning Estimation and Control with Air Bleed Injection[END_REF].

Fuel cell thermal model

We suppose that solid part of the fuel cell is isothermal and the temperature of fluids at the outlet of the cell is the same as the temperature of the fluids inside the cell. We have also supposed that the specific heat capacity of elements are constant. We have considered four thermal equations in each cell, i.e., an equation for the solid mass equation, two equations for anode and cathode gases, and an equation for the coolant. The fuel cell solid temperature T s is obtained with the first law of thermodynamics, i.e.,

m s C p,s dT s dt = W trans + W prod ( 18
)
where m s and C p,s are the mass and the specific heat capacity of the solid, respectively. W trans which is the heat transfer to the solid is obtained as the total thermal energy transfered from existing species in the electrodes to the solid. Using the Newton law of cooling, we get

W trans = h a S a (T out a -T s ) + h c S c (T out c -T s ) + h cool S cool (T out cool -T s ) (19)
where h a S a , h c S c , and h cool S cool are the heat transfer coefficients between the solid and the gas in the anode, the gas in the cathode, and the coolant, respectively. As shown in Fig. 3, it is assumed that the chemical power of hydrogen is shared between electrical power and the heat generation from reactions. Supposing that reactions take place on

Heat power

Electrical power

E 0 Er V cell , Cell Voltage (V) J, Current density (mA/cm 2 )
FIG. 3: Thermal balance in the fuel cell the catalyst (solid), we obtain the heat generation as a function of the current density J and the cell voltage V cell , as given in 20.

W prod = (E 0 -V cell )I (20)
Differential equations describing thermal behavior of fluids in electrodes are

V k ρ k C p,k dT out k dt = W k,evap/cond + W k,trans , k ∈ {a, c} (21)
where W k,evap/cond is the heat transfered from evaporation/condensation which is obtained from

W k,evap/cond = q k L v (Y in wk -Y out vk ), k ∈ {a, c} (22)
where L v is the latent heat of water. W k,evap/cond is the heat transfered from the solid to gas is obtained from

W k,trans = q in k C p,k (T in k -T out k ) + h k S k T s -T out k , k ∈ {a, c} (23)
Finally, the thermal behavior of the coolant the is given by:

m cool C p,cool dT out cool dt = q in cool C p,cool T in cool -T out cool + h cool S cool (T s -T out cool ) (24)

Fuel cell complete model

The complete model of the fuel cell is an ODE with eight state variables. Defining state and input vectors as

X = [Y out h , Y out o , Y out w,a , Y out w,c , θco, T out a , T out c , T out cool ] T U = [I, Y in h , Y in co , Y in air , Y in o , Y in w,a , Y in w,c
, q in a , q in c , q in cool , P we can summarize the model as follows

m s C p,s dTs dt = h a S a T out a -T s + h c S c T out c -T s + h cool S cool T out cool -T s + (E 0 -V cell )I dX dt = A(X, Ts, U ) V cell = G(X, Ts, U )
This model has been implemented and the experimental data are used to identify and validate the model parameters, see [START_REF] Dipenta | A Reduced Fuel Cell Stack Model for Control and Fault Diagnosis[END_REF]. We have used the Modelica language [START_REF] Fritzson | Principles of Object-Oriented Modeling and Simulation with Modelica 2[END_REF][START_REF] Najafi | New features and new challenges in modeling and simulation in Scicos[END_REF] to model the fuel cell. The Modelica model is then simulated in Scicos [START_REF] Campbell | Modeling and simulation Scilab/Scicos[END_REF]. In this simulation, we have used the input variables which are given in table I. In Fig. 4 the polarization curve of the fuel cell which we obtained is plotted. Although, this assumption is reasonable for many cases, there are situations where this assumption does not hold; e.g., a stack with some faulty cells, or a fuel cell stack whose cells have different thermal or chemical characteristics.

Based on an evaluation of existing stack thermal models in the literature, we build a thermal model for the stack.

We assume a fuel cell stack composed of cells electrically in series with thermal and electrical contacts. The stack model is a 1-D model in which fuel cells are independent and each cell can only influence the others via thermal conduction. Although a cell is composed of several layers, we assume that the cells are lumped plates of width δ and mass m s . Each cell is in thermal contact with two adjacent cells and the input of all cells are identical except the coolant electrodes. In fact, here we assume that the cooling circuit does not pass through all cells. Generally, five or even more cells are cooled down by a single cooling plate. In Fig. 5 the cooling channel layout of our fuel cell stack is given. Thus, Q in cool , the flow rate of the coolant in the stack channels is defined as

q in cool (i) =     
q in cool,0 if i = 5, 10, 15, ...

else

where q in cool,0 is the coolant flow in a 'cooled' cell, see Fig. 5.

00 00 00 00 00 00 00 00 00 00 In the stack, each cell affects the two adjacent cells through thermal conduction. Thus, the model of i th cell in the stack is

               ms,iCs,i dTs,i dt = haSa T out a,i -Ts,i + hcSc T out c,i -Ts,i + h cool S cool (T out cool,i -Ts,i) + (E 0 -G(Xi, Ts,i, Ui))I + ψi dXi dt = A(Xi, Ts,i, Ui) V cell,i = G(X i , T s,i , U i ) (25)
where ψ i is defined as

ψ i =              D δ (Ts,2 -Ts,1) + h inf S inf (T inf -Ts,1) i = 1 D δ (Ts,N-1 -Ts,N ) + h inf S inf (T inf -Ts,N ) i = N D δ (T s,i-1 + T s,i+1 -2T s,i ) 1 < i < N (26)
The terminal voltage of the stack is

V stack = N i=1 V cell,i = N i=1 G(X i , T s,i , U i )
The simulation result of a stack composed of 100 cells is given in Fig. 6. For this simulation we have used the parameters given in the table I. In Fig. 6, the temperatures of the first twenty cells of the fuel cell stack (i.e., T s,i , i = 1, • • • , 20) are given for the time interval, time=[0, 70] seconds. Fig. 7 gives the value of T s,i , i = 1, • • • , 20 at time=70 sec. In this simulation, the initial cell temperatures are 300K. A typical fuel cell stack is composed of at least 50 fuel cells. Modeling such a stack usually ends up with a large stiff differential equation. Consequently, the simulation of the model will be slow. In some applications such as fault detection where the output of a physical system is compared with the output of a model, we need a real-time simulation. Thus, a compromise between low model complexity and high solution accuracy should be found. Therefore, a functional approximation method is necessary to obtain a lower model size as well as an accurate enough modeling and simulation.

IV. MODELING OF THE FUEL CELL STACK AS A 1D PDE

The dynamics of a single fuel cell is described by an ODE and the fuel cell stack model (25,26) developed in section III is composed of N cells in series with thermal and electrical contacts. This model is too large and too slow to be used in control applications. In order to come up with a solution for this problem and obtain a reduced model, the stack model is transformed into a partial differential equation (PDE). For that, we assume that the stack is a uniform mass whose thermal behavior is given by the following general PDE and two Robin's type boundary conditions. For the sake of simplicity, in the rest of the paper we assume that cooling channels pass through all cells. Thus, the model of stack (25,26) can be considered as the PDE (27) discretized at N distinct points.

                         ∂X ∂t = A(X , Ts, U ) c0 ∂T s ∂t = c1 ∂ 2 T s ∂z 2 + c2(T out a -Ts) + c3(T out c -Ts) + c4(T out cool -Ts) + (E 0 -G(X , Ts, U ))I 0 = c5 ∂T s ∂z + c2(T out a -Ts) + c3(T out c -Ts) + c4(T out cool -Ts) + c7(T inf -Ts) + (E 0 -G(X , Ts, U ))I| z = 0 0 = c6 ∂Ts ∂z + c2(T out a -Ts) + c3(T out c -Ts) + c4(T out cool -Ts) + c7(T inf -Ts) + (E 0 -G(X , Ts, U ))I| z = L (27)
where T s , and X are the temperature and state variables of the uniform mass as a function of z and relates to T s,i and X i via

     T s,i ≈ T s (z)| z=(i-1/2)δ , i = 1, • • • , N X i ≈ X (z)| z=(i-1/2)δ , i = 1, • • • , N
The coefficients c j , where i ∈ 0 • • • 7, are unknowns that can be obtained if we discretize the PDE via the finite difference method and compare the resulting ODE with original stack model (25,26). The coefficients can then be computed as given in the following table.

c0 = msCs c1 = Dδ c2 = haSa c3 = hcSc c4 = h cool0 S cool c5 = D c6 = -D c7 = h inf S inf
The obtained PDE can be solved instead of simulating the stack. But the simulation of a PDE is not a trivial task and needs a discretization of the space at several points to have a reasonable accuracy. There are, however, some spectral methods to discretize a PDE using only a small number of discretization points which will be explained in the next section.

V. THE COLLOCATION METHOD

The collocation method, which is a useful method for solving PDEs, was developed originally as a stable, predictable, and simple way to implement pseudo-spectral technique. Because of its reliability, it has become a standard method for solving boundary-value problems by polynomial trial function expansions. This method allows discretizing a PDE with a selection of only 3 to 7 points in the region, comparing to 10 to 20 points in finite difference method [START_REF] Villadsen | Solution of Differential Equation Models by Polynominal Approximation[END_REF][START_REF] Rice | Applied Mathematics and Modeling for Chemical Engineers[END_REF][START_REF] Pereira | Modeling of Dynamic Systems[END_REF].

The formulation of this method is based on choosing a set of trial functions from an orthogonal polynomial sequence, with the discretization points computed as the roots of the polynomial next in the sequence. Similar to the finite difference method, in collocation methods N c points are selected which partition the region into N c + 1 segments.

Then, Φ(z, t) is approximated by a Lagrange interpolation polynomial Φ(z, t) of order N c + 1 using z i as interpolation points. Then, the approximated partial derivatives at z i (corresponding derivatives of Φ(z, t)) are inserted in the PDE to satisfy the PDE only at collocation points z i . To implement this method in a 1-D region with N c collocation points, Φ and its partial derivatives are approximated as linear combinations of basis functions L i , i.e.,

Φ(z, t) ≈ Φ(z, t) = Nc i=1 Φi(t)Li(z) (28a) ⇒ ∂Φ(z, t) ∂z ≈ ∂ Φ(z, t) ∂z = Nc i=1 Φi(t) dLi(z) dz (28b) ⇒ ∂ 2 Φ(z, t) ∂z 2 ≈ ∂ 2 Φ(z, t) ∂z 2 = Nc i=1 Φi(t) d 2 L i (z) dz 2 (28c) ⇒ ∂Φ(z, t) ∂t ≈ ∂ Φ(z, t) ∂t = N c i=1 d Φi(t) dt Li(z) (28d)
L i (z) are the N c + 1 linearly independent basis functions of the N th c order Lagrange polynomial. They are determined by the N c + 1 collocation points z i as follows:

L i (z) = Nc j=1, i =j z -z j z i -z j = Nc k=1 a k z k (29)
The coefficients a k can be computed from the known points z i . Hence, the derivatives of the Lagrange-Polynomial (29) can be obtained easily as given in (30).

dL i (z) dz = Nc k=1 ka k z k-1 (30a)
d 2 L i (z) dz 2 = Nc k=2 k(k -1)a k z k-2 (30b)
Since the PDE is to be satisfied only at the collocation points and the Lagrange polynomials (29) have the property (31). Thus the derivatives of the Lagrange polynomials (31) at the collocation points are obtained to be used to compute (28).

     L i (z j ) = 1 for i = j L i (z j ) = 0 for i = j (31)
The interpolation with the Lagrange polynomials tends to provide oscillating curves yielding a very bad approximation, if the collocation points (z i ) are not well chosen. In order to avoid this problem, the basis functions of the Lagrange polynomial should be a set of orthogonal functions in the interval 0 ≤ ξ ≤ 1 with respect to some weighting function w(ξ) assuming a normalized ξ-domain with L = 1. When the interpolation points ξ i are chosen as roots of Jacobi polynomials, the Lagrange polynomials are orthogonal in the sense of

1 0 w(ξ)L i (ξ)L j (ξ)dξ for i = j, w(ξ) = (1 -ξ) τ1 ξ τ2 (32)
The Jacobi polynomial of order N is defined as:

P (τ1,τ2) N = N i=0 (-1) N -i ζ(N, i)ξ i
The function ζ(N, i) can be computed recursively as follows:

ζ(N, i) = N -i+1 i N +i+τ1+τ2 i+τ2 ζ(N, i -1) ζ(N, 0) = 1
where the parameters τ 1 and τ 2 can be used to influence the position of collocation points (ξ i ). This results from the weighting function w(ξ) = (1ξ) τ1 ξ τ2 . If τ 1 is small, the concentration of collocation points will be toward ξ = 1 and when τ 2 is small, the concentration of collocation points will be around ξ = 0. The concentration will be uniform in the region, if τ 1 = τ 2 .

VI. APPLICATION

As a simple application, the orthogonal collocation method is used to discretize the PDE (27) using six collocation points. The simulation result is given in Fig. 8. In this figure, the temperature of the stack is given as a function of time and z. In order to compare the simulation result, the original fuel cell, composed of 50 cells, was simulated with the same parameters and input vector (U ). The result of two simulations is given in Fig. 9. In this figure, the temperature of stack via two method is given. The initial temperature is 293K and the steady state temperature is about 420K. The simulation shows that two results are nearly identical. In order to demonstrate the efficiency of the collocation method in simulation of cell cell stacks, we simulate a fault in a fuel cell. Consider a fuel cell stack composed of 50 cells, and suppose that due to a fault, cooling circuit in cell numbered 24, 25, 26, and 27 is blocked and these cells are not cooled down properly. We have modeled such a fault by replacing coolant flow rate (q in cool (i)) in the faulty cells in the stack in (24,25,26) by equation (33).

q in cool (i) =      0.5q in cool,0
for i = 24, 25, 26, 27 q in cool,0 other i (33) Similar to the previous example, two methods are used; A model composed of 50 cells with a fault, and another model with orthogonal collocation method. In the collocation method, we assume that the stack is composed of three regions and we discretize each region with four points (N c = 4). So, the total number of collocation points is twelve.

The τ 1 and τ 2 of each domain are chosen so that the concentration of the points around the faulty region become higher. In Fig. 10 and 11, the temperature of the cell is given as a function of z and time. The simulation shows a very good accordance between the simulation of the stack of 50 cells and the simulation of the stack simulated with only 12 cells.

VII. CONCLUSION

The reaction-diffusion model of fuel cell stack which may be composed of at least 50 cells is modeled and simulated with a small number of cells. We have used the orthogonal collocation method to sample the stack at only a few points, which is quite lower than the number of cells in the stack. The main principle of the orthogonal collocation method is that the material and energy balances remains unchanged and they are still satisfied exactly at collocation points. Therefore the model of stack via the collocation method becomes smaller and easier to simulate. There are several advantages when using the reduced model: the model of the fuel cell stack is no longer discrete, which allows a uniform modeling. Then, this method reduces the size of the process model in terms of number of variables and balances leading to significant computational savings.

FIG. 4 :

 4 FIG. 4: Polarization curve of the fuel cell, with and without air bleed

FIG. 5 :

 5 FIG. 5: Q in cool in the fuel cell stack

FIG. 6 :FIG. 7 :

 67 FIG. 6: Temperature distribution in the stack of 100 fuel cells

FIG. 8 :

 8 FIG.8: Simulation of a fuel cell stack with the orthogonal collocation method using four discretization points

  FIG. 9: Temperature distribution (in time and z) in a fuel cell stack using two simulation methods

FIG. 10 :

 10 FIG. 10: Simulation result of a faulty fuel cell stack with multi-domain collocation method and standard method of the stack of 50 cells

NomenclatureAin=

  act = the active membrane area bi= kinetics constants i ∈ {f h, f c, f o} Ce= specific heat capacity of elements e ∈ { s, a, c, cool } D= heat transfer coefficients between two adjacent cells E i = activation energies i ∈ {bf c, kact, kdif f } E 0 = thermo-neutral potential (1.48 volts) F = Faraday constant (96450 coulombs/mol) h k S k = heat transfer coefficients between electrodes and gas k ∈ {a, c} h cool S cool = heat transfer coefficient between solid and coolant h inf S inf = heat transfer coefficient between end-cells and the environment I, J= current of the cell and the current density (I = JA act)k i = kinetics constants i ∈ {f h, eh, f c, ec, f o, oh, oc, res} k 0 i = kinetics constants i ∈ {act, dif f } L= the fuel cell stack length M f =molar mass of fluids f ∈ { h, o, co, w, air } M n,k = average molar mass of neutral gases in electrodes k ∈ {a, c} P out k = output pressure in electrodes k ∈ {a, c} q in cool = inlet coolant mass flow R= gas constant (8.314 J/(mol K)) R ohm = ohmic resistance of a cell T in k = inlet fluid temperature at the electrodes k ∈ {a, c} T in cool = inlet coolant temperature Ts= solid temperature T out k = outlet fluid temperature at the electrodes k ∈ {a, c} T out cool = outlet coolant temperature T inf = environment temperature V cell = cell voltage V k = volume of electrode channels k ∈ {a, c} Y in h = anode inlet mass fraction of H2 Y in co = anode inlet mass fraction of CO Y in air = anode inlet mass fraction of air bleed Y in o = cathode inlet molar fraction of O2 Y in v,k = inlet mass fraction of vapor water at electrodes k ∈ {a, c} Y out h = anode outlet mass fraction of H2 Y out o = anode outlet mass fraction of O 2 in cathode Y out w,k = outlet mass fraction of water at the electrodes k ∈ {a, c} Greek ∆S 0 l = change in entropy of liquid water formation under the standard condition ∆H 0 l = change in enthalpy of liquid water formation under the standard condition δ= The thickness of a single fuel cell ρ k = average mass density of gas in electrodes k ∈ {a, c} λ k = inlet stoichiometric ratio at the electrodes k ∈ {a, c} α k = water transfer coefficient through membrane in electrodes k ∈ {a, c} and αa = 1 -αc θ co = anode catalyst CO poisoning rate β k = pressure drop parameter in electrodes k ∈ {a, c} Input to the cell out= Output from the cell


Huy-Dinh Vu

Jacques Renault
email: jacques.renault@univ-rennes1.fr

Thierry Roisnel

Nicolas Gouault

Philippe Uriac

Synthesis of 3-substituted indolizidines from amino-ynones derivatives

Keywords: amino-ynone, indolizidine, methanesulfonic acid, zinc chloride

Two types of indolizidine derivatives can be synthesized from amino-ynones. In the presence of methanesulfonic acid, Meyer-Schuster rearrangement furnished vinylogues of indolizidin-3-ones whereas treatment with zinc chloride gave stable 3-alkynyl-1,5,6,7,8,8a-hexahydroindolizine-zinc complexes that were reduced to 3-alkynylindolizidines.

In previous publications, 1,2 the reactivity of ynones bearing N-Boc protected primary amines I in the presence of electrophilic reagents was reported. Treatment with methanesulfonic acid (CH 3 SO 3 H), gave pyrrolidine exocyclic vinylogous amides via Meyer-Schuster rearrangement, 1 whereas zinc chloride (ZnCl 2 ) furnished acetylenic cyclic imines. In both cases a cyclic carbinolamine intermediate was formed, which underwent either a Meyer-Schuster rearrangement with protic acid 1 or dehydration with Lewis acid. 2 In this paper, we present an extension of this study to amino-ynones having a secondary Boc-protected amine. Thus, when applied to the piperidine derivatives II, reaction with CH 3 SO 3 H or ZnCl 2 resulted in formation of vinylogues of indolizidin-3-ones III or very stable cyclic enamine-zinc complexes.

Reduction of these complexes led to the indolizidines IV, which are frequently found as key structural components of numerous bioactive alkaloids (Fig. 1). 3 The starting amino-ynones were easily prepared from racemic pipecolic acid. The method included transformation of the carboxylic acid to an aldehyde, Wittig chain elongation, alkene hydrogenation 4 and conversion to the racemic Weinreb amide 1. Reaction of 1 with acetylides 5 allowed isolation of compounds 2 to 5 in good yields (> 80% from Weinreb amide 1, Scheme 1). Compounds 2-5 were then submitted to reaction with either methanesulfonic acid or zinc chloride.

Five equivalents of methanesulfonic acid were used to perform the N-Boc deprotection. The mixture was then made alkaline with sodium or potassium carbonate, leading to formation of the indolizidin-3one vinylogues 6-8 through a Meyer-Schuster rearrangement 1 (Scheme 2). The very poor yield observed for compound 6 (R = H) either under our conditions or those reported by Georg (HCl in dioxane then Na 2 CO 3 ) 5b is presumably due to its instability; numerous degradation products were observed by TLC analysis. Contrary to a previous report, 5b that described formation of seven-membered enaminones under similar conditions, we obtained indolizidine vinylogous amides 6-8, whose derivatives are seldom reported. 6 Structural assignment was made on the basis of several key analytical features: the 1 H NMR spectrum of 6 showed a resonance signal at 9.38 ppm, corresponding to an aldehydic proton that

N O R Boc N N O Boc CH 3 OCH 3 a or b 1 2 (R = H, 82%) 3 (R = C 3 H 7 , 85 %) 4 (R = C 6 H 5 , 89 %) 5 (R = , 87%) OCH3 N O R Boc 2-4 N 6 (R = H, 15%) 7 (R = C 3 H 7 , 90 %) 8 (R = C 6 H 5 , 86 %) O H R a N O R Seven membered enaminone
would not be observed in a seven member enaminone, which correlates closely with structurally related compounds found in the literature 7 (Fig. 2).

Figure 2. 1 H NMR data of compounds related to 6-8.

The 1 H-coupled 13 C NMR spectrum of 6 (R = H) showed a doublet (J = 158.9 Hz) at 187.0 ppm, in accordance with an aldehydic carbonyl. NOESY correlation between the H5 and the ethylenic H

confirmed the E configuration of the C=C double bond. Finally, the X-ray crystal data of compound 8 (Fig. 3) were in accordance with these findings. 8 (Comins et al. 7a ) deshielded (δC3 163 ppm). [START_REF] Cook | Enamines: Synthesis, Structure, and Reactions[END_REF] Elemental analysis performed with 11 was in accordance with a M 3 L 2 zinc complex structure (see SI). As enamine reduction needs a protic medium in order to form an iminium salt, [START_REF] Cook | Enamines: Synthesis, Structure, and Reactions[END_REF] we attempted direct reduction of these complexes using NaBH 4 10 in ethanol. The 3alkynylindolizidines 12-14 were obtained as a 60/40 mixture of diastereomers, with the major isomers (12a-14a) having H3 and H8a configured cis. This configuration was assigned on the basis of NOESY correlations between H3, H8 and the axial H5 in compound 14a (See SI). Also of note was shielding of the equatorial H5 by the acetylenic substituent in the major isomers (Scheme 3). The yields of the isolated compounds were strongly affected by alkyne substituent R and isolation was complicated by difficult chromatographic separation. The best result (63% for 14) was obtained with a bulky 5methoxynaphtyl group that facilitated chromatographic purification.

Figure 1 .

 1 Figure 1. Overall aim of this work.

Scheme 1 .

 1 Scheme 1. Synthesis of starting materials 2 to 5. Reagents and conditions: (a) H-C≡CMgBr, THF,

Scheme 2 .

 2 Scheme 2. Synthesis of amide vinylogues 6-8. Reagents and conditions: (a) CH 3 SO 3 H (5 equiv), CH 2 Cl 2 /CH 3 OH, RT, 4 h then Na 2 CO 3 , 12 h.

Figure 3 .

 3 Figure 3.ORTEP drawing of the X-ray crystal structure of 8.

Scheme 3 .Scheme 4 .

 34 Scheme 3. Synthesis of zinc complexes 9-11 and their reduction to 12-14. Reagents and conditions:

Acknowledgments

The authors are deeply grateful to Fabian Lambert (CRMPO) for X-ray analyses, to the Ministry of Education and Training of Vietnam for financial support and to Daniel Best for his contribution to this work.

Supplementary data

Supplementary data associated with this article can be found, in the online version, at http….