G Zheng

Y Orlov

W Perruquetti

J.-P Richard

Finite time observer design for nonlinear impulsive systems with impact perturbation

This paper investigates the observer design problem of nonlinear impulsive systems with impact perturbation. By using the concept of normal form, this paper proposes a full order finite time observer, which guarantees the finite time convergence independent of the impact perturbation. Reduced order observer is additionally developed for the proposed normal form. An illustrative example is given in order to illustrate the capability of the proposed method.

Introduction

Impulsive systems are widely studied in chemical, engineering, and biological fields [START_REF] Lakshmikantham | Theory of Impulsive Differential Equations[END_REF]]) where instantaneous changes of the state variables are admitted at various time instants. Moreover, from the theoretical point of view, some non smooth dynamical systems ([Brogliato(1999)]) and sampleddata systems ([Sun (1993), [START_REF] Sagfors | H ∞ control of multirate sampled-data systems: A state-space approach[END_REF]]) can be treated from an impulsive systems point of view as well.

Fundamental problems for impulsive systems, such as observability, reachability and controllability, have been widely investigated for different types of impulsive systems [START_REF] Lakshmikantham | Theory of Impulsive Differential Equations[END_REF], Yang(2001), [START_REF] Li | [END_REF]]). Particularly, in the case of linear impulsive systems with constant coefficients, the concepts of observability/reachability and controllability proved to be equivalent to those of linear time-invariant systems [START_REF][END_REF]]). Controllability and observability results were developed in [Guan (2002)] for linear impulsive systems where the controlled actuation were available for the continuous-time dynamics at impact time instants only and the impulsive effects were limited to scalings of the state. For the same type of linear impulsive system, [START_REF] Xie | Controllability and observability of a class of linear impulsive systems[END_REF]] used a geometric framework to generalize the results of [Guan (2002)]. For nonlinear systems, finite time observers have been exhaustively discussed in the literature, by using different methods, such as sliding mode technique [START_REF] Perruquetti | A note on sliding observer and controller for generalized canonical forms[END_REF]]), homogeneity ([Perruquetti (2008)]) delay measurement ([Sauvage (2007)]), output injection ( [START_REF] Engel | A continuous-time observer which converges in finite time[END_REF]]) and algebraic methods ( [START_REF] Fliess | Reconstructeurs d'état[END_REF], Barbot (2007)]). How-ever, few results are reported on finite time observer design for impulsive systems.

In this paper, we study observer design problem for nonlinear impulsive systems with persistent impact perturbation. This work is motivated by the fact that the perturbation cannot be avoided at each impact instant for some practical applications, for example, for the walking robot, the perturbation always exists when a robot's feet touches the ground. For some advanced controllers relying on the unmeasured states, a finite time observer, which can always estimate system states during a prescribed interval, is required under persistent impact perturbations.

Recently, [START_REF] Raff | An impulsive observer that estimates the exact state of a linear continuous-time system in predetermined finite time[END_REF]] proposed to simply couple two Luenberger observers to form a finite time observer for linear systems. This paper adopts this idea to synthesize a finite time observer for nonlinear impulsive systems with persistent impact perturbation, whose convergence duration is independent of the persistent perturbation at each impact instant, and can be fixed a priori.

The concept of a normal form is additionally adopted in order to study the finite time observer design problem for nonlinear impulsive systems. Roughly speaking, such a form represents a class of systems possessing the same properties of the stability, controllability and so on. A normal form used for the observability study and observer design can be found in [START_REF] Krener | Linearization by output injection and nonlinear observer[END_REF], [START_REF] Krener | Nonlinear observer with linearizable error dynamics[END_REF], [START_REF] Boutat | [END_REF]), Zheng (2007)] and in the references, quoted therein. The basic idea is to deduce a diffeomorphism, transforming a system in question into a common and simple normal form for which the observability and observer design problems have been already solved. Then, by inverting the deduced diffeomorphism, the observability and observer design problems are solved for the original system as well.

This paper is organized as follows. Notations and problem statement are given in Section 2. Section 3 presents an appropriate normal form and necessary and sufficient conditions for a transform of a nonlinear impulsive system into such a normal form to exist. The corresponding full order observer and reduced order observer are proposed in Section 4. Section 5 highlights the capability of the proposed method by means of an illustrative example.

Notations and problem statement

Consider the following nonlinear impulsive system

       ξ(t) = f (ξ(t)) + g (ξ(t)) u(t), t ∈ [t k-1 , t k ) with k ∈ Z + ξ(t + ) = γ(ξ(t -)) + Γp (t -), t = t k y(t) = h(ξ(t)) ξ(t + 0 ) = ξ 0 , t 0 = 0 (1)
with an unknown persistent impact perturbation Γp . Hereinafter, ξ ∈ R n , u ∈ R p and y ∈ R stand for the state, the admissible control input and the output, respectively; f (ξ), g(ξ), h(ξ) and γ are known smooth functions of appropriate dimensions; Z + represents the set of positive integer, and

t 0 < t 1 < • • • < t k < • • • are impact time instants such that D min = inf k∈Z + {t k -t k-1 } > 0 (2)
and

D max = sup k∈Z + {t k -t k-1 } < ∞ and lim k→∞ t k = ∞.
The left and right values at impulse time t k are defined as usual

ξ(t + k ) = lim h→0 ξ(t k + h), ξ(t - k ) = lim h→0 ξ(t k -h).
with h > 0. For certainty, it is assumed that ξ(t - k ) = ξ(t k ), thereby yielding a solution of (1) to be left continuous at t k . It is also assumed throughout that (1) is observable ( [START_REF] Hermann | Nonlinear controllability and observability[END_REF]]) in the impact-free case.

Since in what follows, the underlying system (1) is periodically affected by unknown impulsive disturbances, many conventional concepts, such as that of asymptotic observer, is to be re-worked for impulsive systems with impact uncertainties.

Given

t k-1 ≤ t < t k with k ∈ Z + , denote x(t) = x(t, t + k-1 , x(t + k-1
)) a solution of (1) initialized at t + k-1 with x(t + k-1 ).

Definition 1 The following impact system

   ż = f 1 (z, u, y), t ∈ [t k-1 , T k ) ∪ (T K , t k ) z(t + ) = g 1 (z(t -)), t = t k-1 + T k z(t + ) = g 2 (z(t -)), t = t k
with z ∈ R q where q ≥ n and some user chosen functions f 1 , g 1 and g 2 , is an exponential observer of (1) with the output η : R q → R n if the following inequality

||η (z(t)) -x(t)|| ≤ ae -b(t-t k-1 ) ||η (z(t k-1 )) -x(t k-1 )||, for t k-1 ≤ t < t k is satisfied for each interval [t k-1 , t k ), k ∈ Z + and
some positive constants a and b.

It is a finite time observer of (1) with the output η : R q → R n if there exist corresponding

σ k < t k -t k-1 , k ∈ Z + such that ||η (z(t)) -x(t)|| = 0, when t ∈ [t k-1 + σ k , t k ), k ∈ Z + .
This system is said to be a uniform finite time observer of (1) iff σ = sup k∈Z + {σ k } ∈ (0, D min ) where D min is specified in (2). Such a σ is further referred to as a settling time estimate of the observer.

To motivate the above definition it should be noted that although one can design an asymptotical observer for (1) on each interval [t k-1 , t k ), k ∈ Z + , however, the convergence of such an observer on a specific interval [t k-1 , t k ) would depend on the perturbation Γp (t - k-1 ), thereby possibly resulting in inappropriate state estimate due to accumulating an estimation error.

Transformation to normal form

By using the concept of normal form, this section is devoted to presenting a normal form for which an observer is proposed, and then deducing necessary and sufficient conditions under which (1) can be transformed into this normal form.

Throughout this article, L k f h is the k th Lie derivative of h along f and θ i = dL i-1 f h stands for the associated differentiation for 1 ≤ i ≤ n with dL 0 f h = dh. Since once not affected by impulses, system (1) is assumed to be locally observable, the 1-forms θ i are thus linearly independent, i.e. rank{θ i , 1 ≤ i ≤ n} = n. Then, one can construct the well-known [START_REF] Krener | Linearization by output injection and nonlinear observer[END_REF]

]) frame τ = (τ 1 , ..., τ n ) (3)
with the first vector field τ 1 given by the following algebraic equations

θ i (τ 1 ) = 0 for 1 ≤ i ≤ n -1 θ n (τ 1 ) = 1
whereas the rest of vector fields is obtained by iterating on i:

τ i+1 = [τ i , f ] for 1 ≤ i ≤ n -1 where [, ] denotes the conventional Lie bracket. Setting θ = (θ 1 , • • • , θ n ) T , (4)
let us introduce

Ξ = θτ =        0 0 • • • 0 1 0 0 • • • 1 l 2,n . . . . . . . . . . . . . . . 0 1 • • • l n-1,n-1 l n-1,n 1 l n,2 • • • l n,n-1 l n,n        (5) with l i,j = θ i τ j for 2 ≤ i ≤ n and n -i + 2 ≤ j ≤ n. It is clear that Ξ is invertible. With this in mind, we denote ω = Ξ -1 θ. (6)
Thus, we arrive at the following result which can be viewed as a natural extension of [START_REF] Krener | Linearization by output injection and nonlinear observer[END_REF]] to nonlinear systems with impacts.

Theorem 1 There exists a diffeomorphism x = φ(ξ) which transforms nonlinear impulsive system (1) into

       ẋ(t) = Ax(t) + β (y(t)) + ρ (y (t) , u(t)) , t ∈ [t k-1 , t k ) with k ∈ Z + x(t + ) = γ(x(t -)) + Γ p (t -), t = t k y(t) = Cx(t) x(t + 0 ) = x 0 = φ(ξ 0 ), t 0 = 0 (7)
where

A =        0 0 • • • 0 0 1 0 • • • 0 0 0 1 • • • 0 0 . . . . . . . . . . . . . . . 0 0 • • • 1 0        , C = (0, • • • , 0, 1) γ (x) = φ γ(ξ(t -)) | ξ=φ -1 (x) Γ p (t -) = φ(γ(ξ(t -)) + Γp (t -)) -φ(γ(ξ(t -))) | ξ=φ -1 (x)
if and only if one of the following conditions holds:

1.

[τ i , τ j ] = 0 for 1 ≤ i ≤ n, 1 ≤ j ≤ n and [g, τ l ] = 0 for 1 ≤ l ≤ n -1; 2. dω = 0 and [g, τ l ] = 0 for 1 ≤ l ≤ n -1.
Proof 1 This theorem is based on the result in [START_REF] Krener | Linearization by output injection and nonlinear observer[END_REF]] where nonlinear systems without impact perturbations were considered. The first part of this proof is based on [Zheng (2007)], but in order to keep the paper self contained, we still detail here the proof of the first part of this theorem.

To begin with, let us show the equivalence between the first condition and the second condition of the theorem. For this purpose, compute the differential dω of ω on two vector fields X = τ i and Y = τ j :

dω(X, Y ) = L X ω(Y ) -L Y ω(X) -ω([X, Y ]),
where L X is the Lie derivative in X direction. Since ω(τ i ) = e i and ω(τ j ) = e j are constant, then one has dω(τ i , τ j ) = ω([τ i , τ j ]). Taking into account that ω is an isomorphism, then one obtains the following equivalence

dω = 0 ⇔ [τ i , τ j ] = 0
which results in the equivalence between the afore-mentioned conditions of the theorem.

Necessity: Next let us proof that once (1) can be transformed into (7) through a diffeomorphism x = φ(ξ), Condition 1 of the theorem is satisfied. Indeed, just in case, one has φ * (τ

i ) = ∂ ∂xi for 1 ≤ i ≤ n and φ * (g) = ρ(y, u) ∂ ∂xi with y = x n , thereby ensuring that [φ * (τ i ), φ * (τ j )] = φ * ([τ i , τ j ]) = 0 for 1 ≤ i ≤ n and 1 ≤ j ≤ n and [φ * (τ l ), φ * (g)] = φ * ([τ l , g]) = 0 for 1 ≤ l ≤ n -1. Moreover, since ξ = φ(x) is a diffeomorphism, one has [τ i , τ j ] = 0 for 1 ≤ i ≤ n and 1 ≤ j ≤ n and [τ l , g] = 0 for 1 ≤ l ≤ n -1.
Sufficiency: In order to prove the sufficiency suppose that [τ i , τ j ] = 0 for 1 ≤ i ≤ n and 1 ≤ j ≤ n. Then according to Poincaré theorem, there exists

a local diffeomorphism x = φ(ξ) such that ω = dφ = φ * , φ * (τ i ) = ∂ ∂xi for 1 ≤ i ≤ n. Thus one has ∂φ * (f ) ∂x i = φ * ([τ i , f ]) = φ * (τ i+1 ) = ∂ ∂x i+1
for 1 ≤ i ≤ n -1,and by integration one arrives at

φ * (f ) = Ax + β (y) .
Moreover, one has ∂φ * (g)

∂x i = φ * ([τ i , g]) = 0
for 1 ≤ l ≤ n -1, and therefore, one obtains

ẋ = φ * (f ) + φ * (g) (8) = Ax + β (y) + ρ(y, u).
By definition ( 6) of ω, one has

ω = Ξ -1 θ =        * * • • • * 1 * * • • • 1 0 . . . . . . . . . . . . . . . * 1 • • • 0 0 1 0 • • • 0 0               dh dL f h . . . dL n-2 f h dL n-1 f h        =        * * . . . * dh        which ensures that x n = φ n (ξ) = h (ξ), thereby yielding y = Cx with C = (0, • • • , 0, 1).
Finally, with the deduced diffeomorphism x = φ(ξ), it is straightforward to check that

x(t + ) = γ(x(t -)) + Γ p (t -) at t = t k where γ and Γ p are specified in the statement of the theorem. Hence, system (1) can be transformed into (7) through the diffeomorphism x = φ (ξ) provided that Condition 1 of the theorem is satisfied. The proof of the theorem is thus completed.

Theorem 1 gives necessary and sufficient conditions which guarantee the equivalence, via a diffeomorphism x = φ(ξ), between the nonlinear impulsive systems (1) and the proposed normal form (7). Due to this equivalence, one can design an observer for the transformed form (7) to estimate x, and then by applying the inverse of the diffeomorphism one can obtain the estimation of ξ, instead of designing observers for (1) to estimate ξ directly. It should also be noted that even if the observer is global for the transformed form (7), only a local state estimation of the original system (1) is available with the present approach that can become singular, dependent of the properties of the chosen diffeomorphism [START_REF] Ciccarella | [END_REF]]). Thus motivated, the subsequent synthesis of full order and reduced order finite time observers is confined to the transformed system (7).

4 Observer design and observation error analysis

Uniform finite time observer

Provided that either condition of Theorem 1 is satisfied, the nonlinear system (1) can be represented, via a diffeomorphism, in the normal form (7) where the pair (A, C) is observable. Hence one can always find two different vectors K 1 and K 2 , such that

P i = (A -K i C) (9) is Hurwitz for 1 ≤ i ≤ 2.
For a prescribed small positive constant σ ∈ (0, D min ) where D min is defined in (2), it is obvious that one can always find two different Hurwitz matrices P 1 and P 2 in (9), such that the matrix (I n×n -e P2σ e -P1σ ) is well-posed, thus is invertible. Hence the following matrix

Λ = (Λ 1 , Λ 2 ) = (I n×n -e P2σ e -P1σ ) -1 (-e P2σ e -P1σ , I n×n ) ∈ R n×2n (10)
is well-posed. Then the following result is in order.

Theorem 2 Let either Condition 1 or Condition 2 of Theorem 1 be satisfied. Given a prescribed positive constant σ ∈ (0, D min ), let K 1 and K 2 be two different vectors such that P i , i = 1, 2 in (9) are Hurwitz, and Λ in ( 10) is well-posed.

Then the following coupled impulsive dynamics

           żi = Az i + β(y) + ρ(y, u) + K i (y -Cz i ), t = t k-1 and t = t k-1 + σ z i (t + ) = Λ z 1 (t -) z 2 (t -) , t = t k-1 + σ z i (t + ) = γ (z(t -)) , t = t k z i (t + 0 ) = z 0 , t 0 = 0 (
11) with i = 1, 2 and the output η = z 1 (alternatively, with the output η = z 2 ) form a uniform finite time observer for (7) with the prescribed settling time estimate σ.

Proof 2 Denote ǫ i = z i -x for i = 1, 2. For a prescribed σ ∈ (0, D min ), one has t k-1 + σ < t k . When k = 1, (7) and (11) yield ǫi (t) = (A -K i C)ǫ i (t) = P i ǫ i (t), for t ∈ [t 0 , t 0 + σ) which ensures that z 1 ((t 0 + σ) -) = x(t 0 + σ) + e P1σ (z 0 -x 0 ) z 2 ((t 0 + σ) -) = x(t 0 + σ) + e P2σ (z 0 -x 0 ) (12)
Multiplying ( 12) by Λ specified by (10) results in

x(t 0 + σ) = Λ z 1 ((t 0 + σ) -) z 2 ((t 0 + σ) -)
The above equality, combined with (11), yields

z i ((t 0 + σ) + ) = x(t 0 + σ),
thereby ensuring the exact estimation of x(t) when t 0 + σ ≤ t ≤ t - 1 . Since z i ((t 0 + σ) + ) = x(t 0 + σ) for i = 1, 2, thus z i of system (11) has the same dynamics as x defined in (7) when t ∈ (t 0 + σ, t - 1 ], and it implies

z i (t - 1 ) = x(t - 1 ), for i = 1, 2 and z i (t + 1 ) = γ(z i (t - 1
)), for i = 1, 2 which can be seen again that (11) has the same initial condition z i (t + 1 ) for the two coupled systems. By induction, for any k ∈ Z + and the corresponding time interval [t k-1 , t k ), one has

z 1 (t + k-1 ) = z 2 (t + k-1 ) z i ((t k-1 + σ) + ) = x(t k-1 + σ) z i (t) = x(t), for t k-1 + σ ≤ t ≤ t - k that completes the proof of Theorem 2.
Remark 1 The validity of Theorem 2 relies on the well-posedness of the matrix (I n×n -e P2σ e -P1σ ). Thus the prescribed finite time σ cannot be too close to zero, otherwise, this matrix becomes ill-posed. In this regard, given a fixed σ, the higher negative real parts of eigenvalues of both P 1 and P 2 are chosen as well as the larger distances between the real parts of eigenvalues of P 1 and P 2 are the better for inverting the matrix (I n×n -e P2σ e -P1σ ) numerically.

Observation error analysis

For the proposed observer, it is proved that it estimates the state of the original system in finite time for each interval [t k-1 , t k ). Because of the existence of periodically acting impulse perturbation, it is important to globally characterize the influence of those impulse perturbations.

Theorem 3 The observation error of the proposed finite time observer ( 11) is bounded, and the following L 2 -gain inequality is satisfied:

∞ t + 0 ǫ i (t) T ǫ i (t)dt ≤ ρ i ∞ k=1 ||Γ p (t k-1 )|| 2 2 (13)
with the attenuation level

ρ i = e 2σλ max i -1 2λ max i (14)
where Γ p (t k-1 ) = Γ p (t - k-1 ) for 1 ≤ k with Γ(t 0 ) = z 0 -x 0 being the observation error of initial conditions, and λ max i is the maximum real part of the eigenvalues of

Pi+P T i 2 defined in (9).
Proof 3 As shown in the proof of Theorem 2, for any k ∈ N and for a fixed interval [t k-1 , t k ), one always has ǫi (t) = P i ǫ(t), when t k-1 ≤ t < t k-1 + σ where ǫ(t k-1 ) = Γ p (t k-1 ) with Γ p (t k-1 ) = z 0 -x 0 , and

ǫ i (t) = 0, when t k-1 + σ ≤ t < t k
Then one has

ǫ i (t) = e Pi(t-t k-1 ) Γ p (t k-1 ), when t k-1 ≤ t < t k-1 + σ which yields ||ǫ i (t)|| 2 2 = ||e Pi(t-t k-1 ) Γ p (t k-1 )|| 2 2 ≤ e 2λ max i (t-t k-1 ) ||Γ p (t k-1 )|| 2 2
thus arriving at

t k t + k-1 ǫ i (t) T ǫ i (t)dt = t k-1 +σ t + k-1 ||e Pi(t-t k-1 ) Γ p (t k-1 )|| 2 2 dt ≤ ||Γ p (t k-1 )|| 2 2 t k-1 +σ t + k-1 e 2λ max i (t-t k-1 ) dt = e 2λ max i σ -1 2λ max i ||Γ p (t k-1 )|| 2 2 = ρ i ||Γ p (t k-1 )|| 2 2 It follows ∞ t + 0 ǫ i (t) T ǫ i (t)dt = ∞ k=1 t k t + k-1 ǫ i (t) T ǫ i (t)dt ≤ ρ i ∞ k=1 ||Γ p (t k-1 )|| 2
2 . The proof is completed.

Remark 2 The attenuation level ρ i given by ( 14) depends only on the maximum real part of the eigenvalues of P i and the prescribed finite time σ, thus it is straightforward to obtain that lim λ max i →0 -ρ i = σ, and lim

λ max i →-∞ ρ i = 0
which yields the following bounds on the observation error:

0 ≤ ∞ t + 0 ǫ i (t) T ǫ i (t)dt ≤ σ ∞ k=1 ||Γ p (t k-1 )|| 2 2
Since the pair (A, C) is observable the values of λ max i , i = 1, 2 can be prespecified to have an arbitrarily large magnitude if vectors K 1 and K 2 are properly chosen, so that ρ i can be tuned to be small enough in order to perfectly attenuate the influence of the disturbance to an arbitrarily small attenuation level, given a priori. The role of the map Γ p (t k-1 ) represents the initial error introduced by the impact disturbance at each t k-1 .

Robustness with noisy measurement

It is well known that the measurements are usually corrupted by noises, thus this section discusses the robustness issue for the proposed observer (11) with respect to noisy measurements. Under this situation, without loss of generalities, let us suppose that the output of the transformed normal form ( 7) is corrupted by a bounded noise µ(t) as follows:

y(t) = Cx(t) + µ(t) (15)
then one can state the following theorem.

Theorem 4 For the transformed canonical form (7) with the noisy output (15), the observation error for the proposed observer (11) at the prescribed time t k +σ is bounded, independent of the unknown persistent impact perturbation Γ p (t - k ), and satisfies the following inequality:

||ǫ 1 ((t k + σ) + )|| = ||ǫ 2 ((t k + σ) + )|| ≤ κµ 0 ( 16
)
where µ 0 = sup t≥0 {µ(t)} and

κ = ρ1 ||K 1 || + ρ1 ||K 1 || + ρ2 ||K 2 || |1 -||e (P2-P1)σ || 2 | with ρ1 = e λ max 1 σ -1 λ max 1 and ρ2 = e λ max 2 σ -1 λ max 2
, σ is the prescribed small positive constant, K i and P i for 1 ≤ i ≤ 2 are given in (9).

Proof 4 Firstly, let us prove the boundedness of Λ defined in (10). Since

||I n×n -e P2σ e -P1σ || 2 ≥ |||I n×n || 2 -||e P2σ e -P1σ || 2 | = |1 -||e P2σ e -P1σ || 2 | thus ||Λ 2 || 2 = ||(I n×n -e P2σ e -P1σ ) -1 || 2 ≤ 1 |1 -||e P2σ e -P1σ || 2 | and ||Λ 1 || 2 = ||I n×n -Λ 2 || 2 ≤ 1 + ||Λ 2 || 2 ≤ 1 + 1 |1 -||e P2σ e -P1σ || 2 |
Secondly, let us prove the boundedness of the observation error at the prescribed time σ. For the sake of simplicities, consider firstly t ∈ (t 0 , t 1 ). With the noisy output (15) and following the proof of Theorem 2, for 1 ≤ i ≤ 2 one has

ǫi (t) = (A -K i C)ǫ i (t) + K i µ(t) = P i ǫ i (t) + K i µ(t), for t ∈ [t 0 , t 0 + σ)
After the multiplication with Λ defined in (10), this gives the following observation error when t = t 0 + σ:

ǫ i ((t 0 +σ) + ) = Λ 1 t0+σ t0 e P1(t0+σ-τ ) K 1 µ (τ ) dτ +Λ 2 t0+σ t0 e P2(t0+σ-τ ) K 2 µ (τ ) dτ which implies that ||ǫ i ((t 0 + σ) + )|| ≤ µ 0 ||K 1 ||||Λ 1 || 2 t0+σ t0 ||e P1(t0+σ-τ ) ||dτ + µ 0 ||K 2 ||||Λ 2 || 2 t0+σ t0 ||e P2(t0+σ-τ ) ||dτ ≤ µ 0 ||K 1 ||||Λ 1 || 2 t0+σ t0 e λ max 1 (t0+σ-τ ) dτ + µ 0 ||K 2 ||||Λ 2 || 2 t0+σ t0 e λ max 2 (t0+σ-τ ) dτ ≤ µ 0 ||K 1 ||||Λ 1 || 2 e λ max 1 σ -1 λ max 1 + µ 0 ||K 2 ||||Λ 2 || 2 e λ max 2 σ -1 λ max 2 ≤ µ 0 ||K 1 ||||Λ 1 || 2 ρ1 + µ 0 ||K 2 ||||Λ 2 || 2 ρ2 = κµ 0
Since at each time t = t k , the observer states are set as z 1 (t + ) = z 2 (t + ) = γ(z(t -)), thus the above argument can be applied as well for any interval [t k , t k+1 ). Therefore, by induction, one can deduce that at each time t = t k + σ, the observation errors are always bounded as (16). Moreover, since κ in (16) does not depend on the unknown persistent impact perturbation Γ p (t - k ), one can conclude that the observation error for the proposed observer (11) at the prescribed time t k + σ is bounded, and does not depend on the unknown persistent impact perturbation.

Remark 3 Due to the existence of noise in the measurement, one can only conclude the boundedness of ǫ i ((t k + σ) + ), which is however equal to 0 for the noise free case. Thus, for t k + σ ≤ t ≤ t - k+1 , the observation errors are influenced only by the measurement noise, but not by the unknown persistent impact perturbation, and this error is governed by the following dynamics:

ǫi (t) = P i ǫ i (t) + K i µ(t)
therefore it can be easily seen that the ǫ i is ISS ([Sontag(1995)]), since the initial condition ǫ i ((t k + σ) + ) and the measurement noise µ(t) are both bounded.

Reduced order uniform finite time observer

The uniform finite time observer given in ( 11) is of dimension 2n, since two coupled Luenberger observers are involved. It is well known that, for linear system, the existence of full order observer implies the existence of reduced order observer, which is used to only estimate the states except the output. Hence we can also design a reduced order observer for the systems with impacts. For this, let us reconsider the proposed normal form (7) which can be decomposed into the following form

           ẋ1 (t) = A 11 x 1 (t) + A 12 x 2 (t) + β 1 (y) + γ 1 (y, u), t ∈ [t k-1 , t k ) ẋ2 (t) = A 21 x 1 (t) + A 22 x 2 (t) + β 2 (y) + γ 2 (y, u), t ∈ [t k-1 , t k ) x(t + ) = γ(x(t -)) + Γ p (t -), t = t k y(t) = x 2 (t) = Cx(t) x(t + 0 ) = x 0 , t 0 = 0 (17)
where x 1 ∈ R n-1 and x 2 ∈ R are the decomposed state vectors; 1) , A 12 = 0 (n-1)×1 ∈ R (n-1)×1 and A 22 = 0 1×1 ∈ R are the corresponding decomposed matrices.

A 11 =        0 0 • • • 0 0 1 0 • • • 0 0 0 1 • • • 0 0 . . . . . . . . . . . . . . . 0 0 • • • 1 0        ∈ R (n-1)×(n-1) , A 21 = [0, • • • , 0, 1] ∈ R 1×(n-
It is worth noticing that the observability of (A, C) in ( 7) ensures that (A 11 , A 21 ) in ( 17) is observable as well. Thus, similarly to full order observer design for (7), one can also find two different vectors K 1 and K 2 , such that

Q i = (A 11 -K i A 21 ) (18)
is Hurwitz for i = 1, 2. Following the same argument as that for full order observer design, for a prescribed small positive constant σ ∈ (0, D min ) where D min is defined in (2), one can also find two different Hurwitz matrices Q 1 and Q 2 defined in (18), such that the matrix (I (n-1)×(n-1) -e Q2σ e -Q1σ ) is well-posed.

Denote Ω = (I (n-1)×(n-1) -e Q2σ e -Q1σ ) -1 (-e Q2σ e -Q1σ , I (n-1)×(n-1) ) ∈ R (n-1)×2(n-1) (19) then one has the following result.

Theorem 5 For a prescribed positive constant σ ∈ (0, D min ), let K i for i = 1, 2 be two different vectors such that Q i in ( 18) is Hurwitz, and Ω in ( 19) is welldefined. Then the following coupled impulsive dynamics

               ζi = (A 11 -K i A 21 )ζ i + β 1 (y) + γ 1 (y, u) + (A 11 K i + A 12 -K i A 22 -K i A 21 K i )y -K i (β 2 (y) + γ 2 (y, u)) , t = t k-1 and t = t k-1 + σ ζ i (t + ) = Ω ζ 1 (t -) + K 1 y ζ 2 (t -) + K 2 y -K i y, t = t k-1 + σ ζ i (t + ) = γ(ζ(t -) + k i y(t -)) -k i y(t -), t = t k z i = ζ i + K i y (20) with initial conditions satisfying ζ 2 (t + 0 ) = ζ 1 (t + 0 ) + (K 1 -K 2 )y(t 0 )
and the output η = [z 1 , y] T (alternatively, with the output η = [z 2 , y] T ) form a uniform finite time observer for (7) with the prescribed settling time estimate σ.

Proof 5 Similarly to the proof of Theorem 2, denote ǫ i = z i -x 1 for i = 1, 2. According to (20), for t = t k-1 + σ where σ ∈ (0, D min ) is a prescribed constant, one obtains

żi = ζi + K i ẏ = Q i z i + A 12 y + K 1 A 21 x 1 + β 1 (y) + ρ 1 (y, u)
and this, combined with (17), yields the following dynamics

ǫi (t) = Q i ǫ i (t), for t = t k-1 + σ When k = 1, the above equality implies z 1 ((t 0 + σ) -) = x 1 (t 0 + σ) + e Q1σ (z 1 (t + 0 ) -x 0 ) z 2 ((t 0 + σ) -) = x 1 (t 0 + σ) + e Q2σ (z 2 (t + 0 ) -x 0 ) (21) Since ζ 2 (t + 0 ) = ζ 1 (t + 0 ) + (K 1 -K 2 )y(t 0 ) then one has z 1 (t + 0 ) = z 2 (t + 0 ) Consequently, multiplying (21) by Ω defined in (19) gives x 1 (t 0 + σ) = Ω z 1 ((t 0 + σ) -) z 2 ((t 0 + σ) -)
From the impulsive dynamics of ζ i in (20), one obtains

z i ((t 0 + σ) + ) = x 1 (t 0 + σ)
and the rest of the proof follows the same line of reasoning as that of the proof of Theorem 2.

Illustrative example

In order to highlight the proposed finite time observers, let us study the following system

                 ξ1 = ξ 2 2 -2ξ 2 1 ξ2-2ξ 2 1 ξ 3 2 -2ξ1ξ 2 2 1+ξ 2 2 + 1 1+ξ 2 2 u, t ∈ [t k-1 , t k ] ξ2 = ξ 1 + ξ 1 ξ 2 2 + ξ 2 ξ(t + k ) = γ(ξ 1 (t - k ) + ξ 1 (t - k )ξ 2 2 (t - k )) ξ 2 2 (t - k ) + Γ p (t - k ), t = t k y = ξ 2 ξ(t + 0 ) = ξ 0 , t 0 = 0 (22)
One has

θ 1 = dξ 2 , θ 2 = (1 + ξ 2 2 )dξ 1 + (1 + 2ξ 2 )dξ 2 thus τ 1 = 1 1+ξ 2 2 ∂ ∂ξ1 , τ 2 = -2ξ1ξ2 1+ξ 2 2 ∂ ∂ξ1 + ∂
ω = Ξ -1 θ = 1 + ξ 2 2 , 2ξ 1 ξ 2 0, 1 = d ξ 1 (1 + ξ 2 2 ) dξ 2
This gives the following diffeomorphism:

φ = ξ 1 (1 + ξ 2
2 ), ξ 2 T through which the studied system proves to be equivalent to the following form:

               ẋ(t) = 0 0 1 0 x(t) + y 2 y + 1 0 u, t ∈ [t k-1 , t k ) x(t + k ) = γ 0 0 1 x(t - k ) + Γ p (t - k ), t = t k y(t) = [0, 1] x(t) x(t + 0 ) = x 0 , t 0 = 0 ( 23
)
where Γ p represents the unknown impact uncertainties. By choosing K 1 = [4, 4] T and K 2 = [100, 20] T , one arrives at Hurwitz matrices P 1 and P 2 . Setting σ ∈ {0.2s, 0.8s}, the states of ( 23) are estimated for each [t k-1 , t k ) at the finite time σ = 0.2s or 0.8s. Thus, one can design a full order finite time observer for (23) of the form (11). In the simulation, we set t k -t k-1 = 5s for k ∈ Z + , γ = -0.8 and the unknown impact uncertainty Γ p (t - k ) belongs to the interval (-0.4, 0.4); u = ÿr -100(x 1 -ẏr ) -20(x 2 -y r ) with y r = | sin(π/5t)| for i ∈ {1, 2}. In the simulation, the Dini derivative ([Garg(1998)]) is considered for the non differential reference y r at time t = 5k for k ∈ Z + . The simulation results, supporting the theory, are depicted in Fig. 1 and2. Good performance of the proposed observer design is concluded from these figures.

Conclusion

Due to the existence of impact perturbation in nonlinear impulsive systems, the convergence of traditional observer highly depends on the initial condition and periodic impact perturbation, thus resulting in an inconsistent performance. This paper conceptually develops the finite time observer design of impact systems. First a normal form is specified for nonlinear impulsive systems and then based on such a form, two types of finite time observers of full and reduced orders are justified to estimate the state of the underlying system in a prescribed small interval in spite of persisten impact perturbations.

[ Zheng (2007)] Zheng, G., Boutat, D., and Barbot, J.P. (2007), "Signal output dependent observability normal form," SIAM Journal of Control andOptimization, 46, 2242-2255.

∂ξ2

  It is straightforward to check that [τ 1 , τ 2 ] = 0 and [τ 1 , g] = [τ 2 , g]

Figure 1 :

 1 Figure 1: Simulation of the full order finite time observer with σ = 0.8s.

Figure 2 :

 2 Figure 2: Simulation of the full order finite time observer with σ = 0.2s.


Amine Dehbaoui

Jean-Max Dutertre

Bruno Robisson

Assia Tria

Investigation of Near-Field Pulsed EMI at IC Level

This article describes the use of a near-field electromagnetic pulse EMP injection technique in order to perform a hardware cryptanalysis of the AES algorithm. This characterization technique is based on the fact that conductors, such as the rails of a Power Distribution Network PDN which is one of the primary EMI risk factors, act as antennas for the radiated EMP energy. This energy induces high electrical currents in the PDN responsible for the violation of the integrated circuit's timing constraints. This modification of the chip's behavior is then exploited in order to recover the AES key by using cryptanalysis techniques based on Differential Fault Analysis (DFA).

I. INTRODUCTION

Electromagnetic interferences (EMI) are unwanted disturbances that affect integrated circuits due to electromagnetic conduction or electromagnetic radiation emitted from an internal or external source. From an EM Compatibility (EMC) point-of-view, EMI whether intentional or not, are considered as a source of noise and interferences. Today, electromagnetic susceptibility of integrated circuits represent also vulnerabilities for hardware security modules like smart-cards. Recently, from a security point of view, researchers, industrials and governmental agencies are focusing with strong interest on these electromagnetic disturbances.

The efficiency of the EM channel is mainly due to the inner properties of EM emissions. Their ability to propagate through different materials is the most interesting one since it allows an attacker to bypass the chip package and/or some EM shields implemented as counter-measure. Moreover, the small size of EM probes permits to focus the perturbation into a small area of the targeted device. This is all the more interesting since it also allows getting around global hardware countermeasures against power glitches such as the use of detached power supplies [START_REF] Shamir | Protecting Smart Cards from Passive Power Analysis with Detached Power Supplies[END_REF] by focusing the EMP injection on reduced die areas.

Two kinds of near-field EM perturbations are usually considered: transient pulses and harmonic emissions. In [START_REF] Poucheret | Local and Direct EM Injection of Power into CMOS Integrated Circuits[END_REF] authors considered the effect of a 1 GHz electric field applied to an IC with an embedded ring oscillator (RO). The main component of that electric field was the transverse one (i.e. parallel to the surface of the chip). The perturbation impacted the output frequency of the RO. Monitoring the effect of that perturbation enabled them to draw a cartography of the sensitive areas of the chip. A cross examination between the layout of the device and the cartography demonstrates that the coupling between the injection antenna and the circuit lies mainly in the Power Distribution Network (PDN). Regarding transient EM pulses, Schmidt et al. reported the use of a spark generator to fault a CRT-based RSA algorithm running on an 8-bits micro-controller [START_REF] Schmidt | Optical and EM Fault-Attacks on CRT-based RSA: Concrete Results[END_REF]. The injected fault leads to a successful attack as it allows them to factorize the RSA modulus. Besides, their experimental setups is characterized by a very large jitter because of the use of the spark-generator. This article describes the use of the EM channel to carry out attacks against a hardware AES embedded in a FPGA with a good temporal and spacial resolution. Transient electromagnetic pulses (EMPs) are injected with a very low jitter on top of the surface of the targets using a 500 µm-diameter magnetic antenna. By doing so, we intend to analyze firstly, the effect of the EMP's polarity on the target, secondly, whether the effect of the EMP on the target is global or local, and finally, the occurrence and the behavior of the faults induced by a very short EM pulse.

The remainder of this article is organized as follows. We describe the EMP injection bench used to generate EMPs in section II. In sections III we study the effect of a localized EMP injected on top of the surface of an FPGA while executing the Advanced Encryption Standard (AES). As a conclusion, section IV summarizes our findings.

II. TEST SETUP

In this section, the near-field EMP pulse injection bench used to induce transient faults is described.

A. Pulsed EMI bench description

The near-field EMP injection bench (Figure 1) is built of a control PC, the device under test (Target), a motorized stage, a pulse generator, and a 500 µm-diameter magnetic antenna. The antenna is moved above the target by means of a high-precision mechanical positioning system (to within 0.01µm minimum). Every element of the bench is controlled by the control PC, and the communication with the target is established through a serial port or a smart card reader.

The pulse generator is capable of generating 200 V with high current (4 A) directed to a target by the antenna 50 Ω loads (with a very low jitter < 30ps), and at repetition rates up to 50 kHz. The output pulse width is variable from 10 to 200 ns. The rise and fall times are 2 ns or less, 20% -80%.

B. Test chip description

The device under test is a FPGA (Xilinx Spartan 3 family Fig. 2). Internal core logic circuits such as the configurable logic blocks CLB and programmable interconnect operate from the 1.2V VCCINT voltage supply inputs. All VCCINT inputs are connected together and to the +1.2V voltage supply, but in order to guarantee problem-free operation, a supply decoupling is present, as described in [START_REF]Spartan-3 fpga family: Complete data sheet[END_REF].

This FPGA implements a hardware 128 bits version of the AES algorithm [START_REF]Announcing the Advanced Encryption Standard (AES)[END_REF]. This algorithm is used for various security purposes. The design is written in VHDL and synthesized for the FPGA. It is built out of three main blocks: a communication and control module (FSM), a key expansion module and a cipher module (ROUND EXE). A manual Place-and-Route stage is performed in order to distinguish between the impacted logical blocks by a comparison with the Floorplan (Figure 3).

We choose to use a 128 bit-wide data path AES and to execute simultaneously on the chip the key expansion and cipher routines. As a consequence, a complete encryption round takes only one clock period, and the whole encryption process is executed in eleven clock periods.

The key expansion routine generates the round keys "on-thefly". For each clock cycle, a new round key is obtained from

III. INJECTING DELAY USING PULSED EMI

This section first reviews the principles of digital circuits synchronous operation in order to introduce DFA. From this, a susceptibility criterion is defined in order to perform a near field EMP injection cartography.

A. Timing design considerations

In this subsection we review shortly the principle of the synchronous behavior of digital ICs. In broad outline, synchronous digital circuits execute digital calculation synchronized by a common clock. They could be described as blocks of combinatorial logic separated with register banks of D flip-flop sharing the same clock as figure 4 shows. The data are generally latched by the registers on the positive edge of the clock. Between two successive clock positive edges, the computed data have to travel from one register to the next. The time needed by the data to propagate through combinatorial logic is called the propagation delay. This delay and an other

T clock > t critical + t setup (1)
Each data bit at the input of a register possesses its own logical cones from the previous register bank associated with its own propagation time. Furthermore, this propagation time is not a constant; it depends highly on the data handled across the logic and the circuit's power supply voltage.

A transient glitch in the power supply voltage, induced by an EMP, will modify the integrated circuit timings. As the master clock frequency remains the same, faults occur in the circuit. These faults may be exploited by an attacker to break a cipher like the AES [START_REF] Biham | Differential fault analysis of secret key cryptosystems[END_REF], [START_REF] Piret | A Differential Fault Attack Technique against SPN Structures, with Application to the AES and KHAZAD[END_REF], [START_REF] Giraud | DFA on AES Advanced Encryption Standard -AES[END_REF]. The attack of Piret et al. [START_REF] Piret | A Differential Fault Attack Technique against SPN Structures, with Application to the AES and KHAZAD[END_REF] is one of the most powerful. It allows to retrieve an AES key with only one pair of correct/faulty ciphertexts when a fault is induced on a single byte of the state before the penultimate MixColumn.

B. Near field susceptibility cartography

In this test, the susceptibility criterion is a timing violation of the AES design constraints. A high susceptibility value refers to a large number of violated paths, while a low susceptibility value refers to an execution without any timing violation. The pulse width value is chosen to match the clock period (TCLK = 10ns), with an amplitude of 100 Volts.

A susceptibility cartography of the design is performed during the last round of the AES. It aims at disclosing the (X,Y) coordinates where the EMP induces a timing violation (i.e faulty computation). The whole surface of the package is exposed to a localized EMP with a displacement step of 500µm (which is also the antenna diameter). The relative distance between the antenna and the surface of the package is set to 500µm. At each location, an EMP is injected during the last round of the AES and the corresponding faulted ciphertext (if Fig. 6. Faults cartography for a Negative EMP any) is retrieved. This process is done for 1,000 encryptions of the same plaintext input, and for every of the 30x30 different locations of the injection antenna on top of the FPGA package.

1) Positive EMP: In this test, and according to the polarity of the EMP, only a coupling between the antenna and the ground network of the circuit will vary the susceptibility criterion. Figure 5 reports the resulting susceptibility cartography for a positive EMP with an amplitude of +100 Volts. In this figure, the red square in the center corresponds to the FPGA die position. At each location, the number of faulted bits are reported.

Considering Figure 5, we observe that the targeted chip is very sensitive to the positive EMP. In fact, a large number of locations over the surface of the package seems more sensitive to a positive EMP than a negative one. This high susceptibility is due to a low resistivity of the ground network in comparison with a decoupled power network.

2) Negative EMP: Figure 6 reports the resulting susceptibility cartography for a negative EMP with an amplitude of -100 Volts. In this figure, the red square in the center corresponds to the FPGA die position. At each location, the number of most frequent faulted bits are reported.

Considering Figure 6, we observe that the effect of the EMP is clearly localized in space. Some locations above the surface of the circuit are more sensitive to the EMP than others. When the EMP is localized in the region near the block cipher, the number of faulted data paths increases. Moreover, we observe a good correlation between the most sensitive coordinates and the position of the ROUNDEXE in Figure 3. This logical block is the place where the critical delay path is located.

Figure 7 shows the behavior of the induced faults for a first random position (X1,Y1,Z) on top of the die's surface (7X7mm 2 ) right in the ROUNDEXE area (the cipher module). 1,000 encryptions were done with random plaintexts and a constant key while injecting EMPs during the last round of the AES calculations. The occurrence rates of both single-bit (i.e. fault affecting a single bit) and multi-bits faults are given.

The path corresponding to the 15 th byte appears to be the most sensitive to the EMP at coordinates (X1,Y1,Z). For this byte, 3% of the faults were single-bit, and 80% of the faults were multi-bits faults. It also reveals a data-dependence of the injected faults to the data handled by the target. This behavior was corroborated by an inspection of the faults.

In fact different faults values ('0' or '1') were obtained for different plaintexts with the same experimental settings. This behavior is as well a feature of faults induced by timing constraints violation (its origin lies in the data-dependence of the data propagation time through combinatorial logic). This is another sign that reinforces the assumption that the fault injection mechanism by means of EMP is related to timing constraints violation.

The same experiment was carried out for two other locations (X2,Y2,Z) and (X3,Y3,Z) on top of the die with the same 1,000 plaintexts used previously. Figures 8 and9 report the corresponding single-bit and multi-bits fault occurrence rates. These three figures (7, 8 and 9) exhibit different occurrence rates: the injection antenna location has an effect on the induced faults and on their related properties. In fact, at coordinates (X1,Y1,Z), the 15 th byte is the most sensitive to the EMP. Whereas, at coordinates (X2,Y2,Z) and (X3,Y3,Z) the most sensitive paths correspond to the 11 th byte and to the 7 th byte respectively. We observed that the faulted paths were different for different locations of the injection antenna.

The evidence of a local effect (i.e. restricted to a part of the device's area) of the EMPs, demonstrates the ability to fault sub-critical paths. In some locations, the most critical one is never faulted. This is a very interesting property (for an attacker) since it is possible to select the disturbed path without always affecting the most critical ones as it is the case for direct power injection techniques.

IV. CONCLUSION

The reported fault injection experiment reveals the ability to inject single-bit and multi-bits faults into the calculations of the AES. These faults were found data dependent. Moreover, a local effect of EMPs was underlined: the injected faults (if any) are modified when the injection antenna location is changed. According to the experiments, EMP injection mechanism may lie in a coupling between the EMP and the internal PDN of the targeted chip. This coupling induces a transient decrease of the voltage applied to the logic of the target. As a consequence, the propagation delays through the logic are increased until faults are induced by the violation of the chip's timing constraints. This property of EMP fault injection is particularly worrying. Indeed, it may allow to bypass many countermeasures intended to prevent direct power injection (e.g. power supply low-pass filtering, use of internal supply monitoring, etc.).

Fig. 1 .

 1 Fig. 1. Near-field EMP injection bench.

Fig. 2 .

 2 Fig.2. The FPGA package footprint[START_REF]Spartan-3 fpga family: Complete data sheet[END_REF]

Fig. 3 .

 3 Fig. 3. Floorplan of the device under test

Fig. 4 .Fig. 5 .

 45 Fig. 4. Synchronous Representation of Digital ICs

Fig. 7 .

 7 Fig. 7. Fault occurrence at coordinates (X 1 ,Y 1 ,Z)

Fig. 8 .Fig. 9 .

 89 Fig. 8. Fault occurrence at coordinates (X 2 ,Y 2 ,Z)

ACKNOWLEDGMENT

This work was funded by the EMAISECI project (ANR-10-SEGI-0012) .


G Perrin

C Soize

D Duhamel
email: duhamel@lami.enpc.fr

C Funfschilling
email: christine.funfschilling@sncf.fr

A POSTERIORI ERROR AND OPTIMAL REDUCED BASIS FOR STOCHASTIC PROCESSES DEFINED BY A FINITE SET OF REALIZATIONS

Keywords: Karhunen-Loève expansion, stochastic processes, a posteriori error, optimal basis AMS subject classifications. 60H35, 60H15, 60H25, 60H40, 65C50

The use of reduced basis has spread to many scientific fields for the last fifty years to condense the statistical properties of stochastic processes. Among these basis, the classical Karhunen-Loève basis corresponds to the Hilbertian basis that is constructed as the eigenfunctions of the covariance operator of the stochastic process of interest. The importance of this basis stems from its optimality in the sense that it minimizes the total mean square error. When the available information about this stochastic process is characterized by a limited set of independent realizations, the covariance operator is not perfectly known. In this case, there is no reason for the Karhunen-Loève basis associated with any estimator of the covariance that are not converged to be still optimal. This paper presents therefore an adaptation of the Karhunen-Loève expansion in order to characterize optimal basis for projection of stochastic processes that are only characterized by a relatively small set of independent realizations.

1. Introduction. For the fifty years, stochastic process analysis has been used in an increasing number of scientific fields, such as uncertainty quantification, material science, seismology, geophysics, quantitative finance, signal processing, control engineering etc. It is indeed an interesting tool for stochastic modeling, forecasting, classification, signal detection and estimation. In most of these applications, the knowledge of these stochastic processes, that we write X = {X(s), s ∈ Ω ⊂ R}, is however limited. Indeed, their statistical properties are generally known through a set of independent realizations that stem from experimental measurements. From these available measurements, several statistical techniques have therefore been developed to identify the distribution P X of stochastic process X. These methods then allow the generation of additional realizations of X that are realistic and representative of the measurements.

When Ω = R, AutoRegressive-Moving-Average (ARMA) models, that were first introduced by Whittle for time series [START_REF] Whittle | Hypothesis testing in time series[END_REF][START_REF] Whittle | Prediction and Regulation by Linear Least-Square Methods[END_REF] and popularized by Box and Jenkins [START_REF] Box | Time Series Analysis: Forecasting and Control[END_REF], allow the description of Gaussian stationary stochastic processes as a parameterized integral of a Gaussian white noise stochastic process. Based on limited knowledge of stochastic process X, these models can therefore be used to emphasize particular properties of X and to extrapolate its value. More recently, in the case Ω = [0, S], with S < +∞, methods based on a two-step approach have given promising results to identify the distribution of a priori non-Gaussian and non-stationary stochastic processes. The first step of these methods is generally the approximation of the stochastic process, X, by its projection, X B (M ) , on a set B (M) of M deterministic functions, which are supposed to be square integrable on Ω and orthonormal. The vector gathering the projection coefficients of X on each element of B (M) , denoted by C, is thus a M -dimension random vector whose components are a priori dependent. Then, the second step is the identification of the multidimensional distribution of C. There are several techniques to perform such an identification, such as the Priori Algebraic Stochastic Modeling (PASM) methods, methods based on the Information Theory and the Maximum Entropy Principle (MEP) [START_REF] Jaynes | Information theory and statistical mechanics[END_REF][START_REF] Soize | Construction of probability distributions in high dimension using the maximum entropy principle. applications to stochastic processes, random fields and random matrices[END_REF], and methods based on a direct projection of random vector C on a polynomial basis of its probability space, such as Polynomial Chaos Expansion (PCE) methods (see [START_REF] Arnst | Identification of bayesian posteriors for coefficients of chaos expansions[END_REF][START_REF] Das | Polynomial chaos representation of spatio-temporal random field from experimental measurements[END_REF][START_REF] Desceliers | Maximum likelihood estimation of stochastic chaos representations from experimental data[END_REF][START_REF] Desceliers | Identification of chaos representations of elastic properties of random media using experimental vibration tests[END_REF][START_REF] Ghanem | On the construction and analysis of stochastic models: Characterization and propagation of the errors associated with limited data[END_REF][START_REF] Marzouk | spectral methods for efficient bayesian solution of inverse problems[END_REF][START_REF] Perrin | Identification of polynomial chaos representations in high dimension from a set of realizations[END_REF][START_REF] Soize | Generalized probabilistic approach of uncertainties in computational dynamics using random matrices and polynomial chaos decompositions[END_REF][START_REF] Soize | Identification of high-dimension polynomial chaos expansions with random coefficients for non-gausian tensor-valued random fields using partial and limited experimental data[END_REF]). In such two-step approaches, the more relevant regarding X the projection basis B (M) is, the lower the dimension M has to be, to guarantee an amplitude of the residue X -X B (M ) lower than a given threshold, and so the easier and the more precise the identification of the distribution of C will be.

As an example of such a projection basis, the Karhunen-Loève (KL) basis has played for the last decades a major role and has been applied in many works [START_REF] Allery | A priori reduction method for solving the two-dimensional burgers' equations[END_REF][START_REF] Atwell | Proper orthogonal decomposition for reduced basis feedback controllers for parabolic equations[END_REF][START_REF] Berkooz | The proper orthogonal decomposition in the analysis of turbulent flows[END_REF][START_REF] Brooks | A karhunen-loève least-squares technique for optimization of geometry of a blunt body in supersonic flow[END_REF][START_REF] Christensen | Evaluation of proper orthogonal decomposition-based decomposition techniques applied to parameter dependent nonturbulent flows[END_REF][START_REF] Huang | Convergence study of the truncated karhunenloève expansion for simulation of stochastic processes[END_REF][START_REF]Track irregularities stochastic modeling[END_REF][START_REF] Le Maître | Spectral Methods for Uncertainty Quantification[END_REF][START_REF] Kunisch | Galerkin proper orthogonal decomposition methods for parabolic problems[END_REF][START_REF] Li | Comparison between karhunen-loève expansion and translation-based simulation of non-gaussian processes[END_REF][START_REF] Ma | Kernel principal component analysis for stochastic input model generation[END_REF][START_REF] Matthies | Stochastic finite elements: Computational approaches to stochastic partial differential equations[END_REF][START_REF] Nouy | Generalized spectral decomposition method for stochastic non-linear problems[END_REF][START_REF] Perrin | Identification of polynomial chaos representations in high dimension from a set of realizations[END_REF][START_REF]Karhunen-loève expansion revisited for vector-valued random fields: scaling, errors and optimal basis[END_REF][START_REF] Phoon | Implementation of karhunen-loeve expansion for simulation using a wavelet-galerkin scheme[END_REF][START_REF] Phoon | Simulation of strongly non-gaussian processes using karhunen-loeve expansion[END_REF][START_REF] Schwab | Karhunen-loeve approximation of random fields by generalized fast multipole methods[END_REF][START_REF] Soize | Identification of high-dimension polynomial chaos expansions with random coefficients for non-gausian tensor-valued random fields using partial and limited experimental data[END_REF][START_REF] Spanos | Galerkin sampling method for stochastic mechanics problems[END_REF][START_REF] Spanos | Karhunen -loève expansion of stochastic processes with a modified exponential covariance kernel[END_REF][START_REF] Wen | A multiscale approach for model reduction of random microstructures[END_REF][START_REF] Williams | The eigenfunctions of the karhunen-loeve integral equation for a spherical system[END_REF][START_REF] Wu | Statistical moving load identification including uncertainty[END_REF]. Indeed, for any integer M , it can be extracted from the KL basis associated with X the M -dimensional family that is optimal in the sense that it minimizes among all the M -dimensional families the L 2 -error associated with X. Mathematically, the KL expansion corresponds to the orthogonal projection theorem in separable Hilbert spaces, for which the Hilbertian basis, {k m , m ≥ 1}, is constructed as the eigenfunctions of the covariance operator of X, defined by the covariance function, R XX , which is assumed to be square integrable on Ω × Ω. When the available information about this stochastic process is characterized by a limited set of ν independent realizations, this covariance function is not perfectly known but can only be estimated. If we define R(ν) as the empirical estimator of R XX , there is however no reason for the eigenfunctions of R to be still optimal.

In reply to this concern, this paper presents a method to identify projection families that are as relevant as possible for X, even if the number of available realizations is relatively small. This method is based on the coupling of an alternative approximation of the covariance operator and an original use of the "leave-one-out" method to a posteriori evaluate the projection errors for X. Section 2 introduces therefore the method we propose to identify optimal projection basis from a set of independent realizations, and Section 3 then illustrates the possibilities of such a method on two examples.

2. Identification of optimal basis from a finite set of independent realizations.

2.1. Theoretical frame. Let (Θ, C, P ) be a probability space. For all M ≥ 1, let L 2 (Θ, R M ) be the space of all the second-order random vectors defined on (Θ, C, P ) with values in R M , equipped with the inner product •, • , such that for all U and V in L 2 (Θ, R M ),

U , V = Θ U T (θ)V (θ)dP (θ) = E U T V , (2.1)
where E [•] is the mathematical expectation. Let P(Ω) be the space of all the secondorder R-valued stochastic processes, indexed by the compact interval Ω = [0, S], where S < +∞. Let H = L 2 (Ω, R) be the space of square integrable functions on Ω, with values in R, equipped with the inner product (•, •), such that for all u and v in H,

(u, v) = Ω u(s)v(s)ds. (2.2)
In addition, the notation • will be used to denote the L 2 norm in P(Ω), such that for all Z in P(Ω),

Z 2 = E [(Z, Z)] . (2.3)
Let X = {X(s), s ∈ Ω} be an element of P(Ω), for which ν independent realizations, {X(θ 1 ), . . . , X(θ ν )}, are supposed to be known. For the sake of simplicity, and without any loss of generality, only centered stochastic processes X are considered in this work, which reads:

E [X(s)] = 0, s ∈ Ω.
(2.4)

It is moreover assumed that the covariance function, R XX , of X is square integrable on Ω × Ω, such that:

R XX (s, s ′ ) def = E [X(s)X(s ′ )] , Ω Ω R XX (s, s ′ ) 2 dsds ′ < +∞. (2.5)
Let B = {b m (s), s ∈ Ω} m≥1 , be a Hilbertian basis of H, such that:

X = m≥1 C m b m , (b m , b p ) = δ mp , C m = (X, b m ) , (2.6)
where the projection coefficients, {C m , m ≥ 1}, are centered random variables that are a priori correlated and δ mp is the kronecker symbol that is equal to 1 if m = p and 0 otherwise. For practical purposes, this basis has to be truncated. For all M ≥ 1, X B (M ) is thus introduced as the projection of X on any M -dimension subset B (M) of B. Finally, the relevance of B (M) to characterize X is analyzed with respect to the L 2 -error, that is denoted by ε 2 , such that:

ε 2 (B (M) ) def = X -X B (M ) 2 .
(2.7)

2.2. Difficulties concerning the identification of the Karhunen-Loève expansion from independent realizations. If S(R) is the set of all square integrable kernel functions defined on Ω × Ω, such that

S(R) = A ∈ L 2 (Ω × Ω, R), | A(s, s ′ ) = A(s ′ , s) ∈ R, (s, s ′ ) ∈ Ω × Ω , (2.8)
the solving of the eigenvalue problem associated with any function A in S(R) can be seen as a generator of a particular Hilbertian basis of P(Ω), b A m , m ≥ 1 , such that:

Ω A(s, s ′ )b A m (s ′ )ds ′ = λ A m b A m (s), λ A 1 ≥ λ A 2 ≥ . . . → 0, b A m , b A p = δ mp .
(2.9)

Due to the orthogonal projection theorem in Hilbert spaces, the KL expansion associated with X, which is denoted by K = {k m , m ≥ 1} and which can be directly computed by solving the eigenvalue problem defined by Eq. (2.8) with A = R XX , is optimal in the sense that, for all M ≥ 1, K (M) = {k m , 1 ≤ m ≤ M } minimizes error ε 2 among all the M -dimensional families of H:

K (M) = arg min B (M ) ∈H M ε 2 (B (M) ) . (2.10)
When dealing with correlated stochastic processes that are only known through a set of ν independent realizations, {X(θ 1 ), . . . , X(θ ν )}, covariance function R XX is however unknown. The eigenvalue problem defined by Eq. (2.8) cannot thus be directly solved. Indirect methods can however be proposed to identify approximations of K from this maximal available information for X. Indeed, noticing that operator R XX is also optimal for the minimization of ε 2 on S(R),

R XX = arg min A∈S(R) ε 2 (B (M) A ) , (2.11)
the best approximation for R XX from these realizations, the most relevant for X the corresponding M -dimensional family. In this prospect, the a priori best evaluation of the covariance function of X is given by the following empirical estimator, R(ν), such that:

R XX (s, s ′ ) ≈ R(ν, s, s ′ ) def = 1 ν ν n=1 X(s, θ n )X(s ′ , θ n ).
(2.12)

The solving of the eigenvalue problem defined by Eq. (2.8) with A = R(ν) allows therefore the definition of a basis of P(Ω), B R(ν) = { k m , m ≤ 1}, whose elements are likely to be well-adapted to X. However, as the rank of R(ν) is by construction lower than or equal to ν, only ν functions are needed to characterize this kernel:

R(ν, s, s ′ ) = ν m=1 λ m k m (s) k m (s ′ ), (s, s ′ ) ∈ Ω × Ω, (2.13)
such that, from Eqs. (2.12) and (2.13), we get:

X(θ n ), k m = 0, 1 ≤ n ≤ ν, m > ν. (2.14)
Hence, set {X(θ 1 ), . . . , X(θ ν )} of available realizations of X is orthogonal to the subset { k m , m > ν}, and the relevance of

B (M) R(ν)
to characterize X will be limited to values of M lower than ν.

When covariance function R XX is not perfectly known, the idea of the following section is therefore to propose an original method to identify M -dimensional projection families that will be adapted to X in the particular case M > ν.

2.3. Proposition of an original subset of S(R). From Eq. (2.11), optimal family K (M) can be seen as the solution of an optimization problem in S(R). This set being vast, this section aims at presenting a particular parametrized subspace of S(R) to identify accurate approximations of K (M) when stochastic process X is only characterized by a finite set of ν independent realizations. In this prospect, another estimator for R XX , that is denoted by R(ν), is introduced, such that for all (s, s ′ ) in Ω × Ω:

R(ν, s, s ′ ) =      1 S-(s ′ -s) S-(s ′ -s) 0 R(ν, x, x + (s ′ -s))dx if S > s ′ -s ≥ 0, 1 S-(s-s ′ ) S-(s-s ′ ) 0 R(ν, x + (s -s ′ ), x)dx if S > s -s ′ > 0, R(ν, s, s ′ ) otherwise. (2.15)
Hence, if stochastic process X is the restriction to Ω of a mean-square stationary stochastic process indexed by s in R, that is to say if R XX (s, s ′ ) only depends on the difference |s -s ′ |, function R(ν) is by construction the best stationary estimator of R XX . The interest of such an estimator is that, although it is computed from only ν realizations of X, its rank is by construction much higher than ν. As a consequence, the number of eigenfunctions associated with R(ν), which will be adapted to X, will be much higher than ν. For any error threshold T , we then define by M (T ) and M (T ) the minimal number of eigenfunctions associated with R(ν) and R(ν) respectively to guarantee a projection error lower than T :

max ε 2 (B ( M) R(ν) ), ε 2 (B ( M) R(ν) ) ≤ T ≤ min ε 2 (B ( M+1) R(ν) ), ε 2 (B ( M+1) R(ν) ) .
(2.16)

From the point of view of the minimization of Eq. (2.11), even if X is actually not mean-square stationary, R(ν) can thus be considered as a better function than R(ν) to characterize X if and only if M ≤ M . From a more general point of view, for α in [0, 1], if we denote by A(α) the function of S(R), such that:

A(α) = α R(ν) + (1 -α) R(ν) ∈ S(R), (2.17)
interesting projection basis for X can be searched as the solution, B

A(α * ) , of the following optimization problem:

     K (M) ≈ B (M) A(α * ) , α * = arg min α∈[0,1] ε 2 (B (M) A(α) ) . (2.18)
Such an approach appears indeed to be efficient when the number of available realizations, ν, is small. First, by construction of any function A(α), the relevance of optimal family B (M)

A(α * ) is at least equal to the relevance of classical families B (M) R(ν) or B (M) R(ν)
. Secondly, for α = 1, the rank of A(α) is by construction higher than ν, such that the number of eigenfunctions associated with A(α), which won't be orthogonal to each available realization of X, will be greater than ν. Finally, the optimization problem defined by Eq. (2.18) being parametrized by a single scalar α, whose definition domain is a bounded interval, the computational costs associated with its minimization can be easily controlled.

Although the convexity of such a problem has not been proved in the general case, we propose in this work to use a dichotomy-based algorithm (for which formulation is presented in detail in Appendix A) for the iterative identification of such an optimal value for α, with T a given error threshold and M ≥ 1. Under this hypothesis and using the proposed algorithm, approximations of α * with a precision of T can be computed from the solving of at most N cost = 2 + p eigenvalue problems defined by Eq. (2.9), where p is the exponent for the smallest power of two that satisfies (1/2) p ≤ T .

2.4.

A posteriori evaluation of the representativeness error. In order to solve the problem defined by Eq. (2.18), a method to a posteriori evaluate ε 2 (B (M)

A(α) ) from only ν available independent realizations, {X(θ n ), 1 ≤ n ≤ ν}, for all α in [0, 1], then is required. Indeed, from the limited set {X(θ n ), 1 ≤ n ≤ ν}, the L 2 -error, ε 2 (B (M)
), corresponding to any M -dimensional family B (M) of H M , cannot be exactly calculated, but has to be evaluated as precisely as possible. Two cases can be distinguished:

• case 1: B (M) is defined without any reference to {X(θ 1 ), . . . , X(θ ν )}.

• case 2: the knowledge of {X(θ 1 ), . . . , X(θ ν )} is used to optimize the representativeness of B (M) . In this case, B (M) depends on the available realizations of X.

2.4.1. Case 1: realizations and projection basis are independent. If B (M) has been computed without any reference to the set {X(θ 1 ), . . . , X(θ ν )}, error ε 2 (B (M) ) can be evaluated from its empirical estimation, ε 2 ν (B (M) ), such that for sufficiently high values of ν:

ε 2 (B (M) ) = E X -X B (M ) , X -X B (M ) ≈ ε 2 ν (B (M) ) def = 1 ν ν n=1 X(θ n ) -X B (M ) (θ n ), X(θ n ) -X B (M ) (θ n ) .
(2.19)

2.4.2. Case 2: the projection basis depends on the available realizations. In order to make projection family B (M) be particularly adapted to stochastic process X, it can be interesting to exploit as much as possible the information about X that is gathered in independent realizations {X(θ 1 ), . . . , X(θ ν )}. In this case, B (M) is dependent on {X(θ 1 ), . . . , X(θ ν )}, and error ε 2 ν can strongly underestimate ε 2 . This phenomenon is generally called over-fitting. For instance, if we define B (M) = {b m , 1 ≤ m ≤ M } as the Gram-Schmidt orthogonalization of the deterministic family of available independent realizations {X(θ 1 ), . . . , X(θ ν )} (see Appendix B for further details about the computation of B (M) ), the projection, X B (M ) , of X on B (M) then verifies, for 1 ≤ n ≤ ν:

X B (M ) (θ n ) = M m=1 (X(θ n ), b m ) b m = X(θ n ).
(2.20)

By construction, error ε 2 ν (B (M) ) thus is always equal to zero, whereas ε 2 (B (M) ) should be in general strictly greater than 0, as the number of available realization, ν, and the dimension of the projection basis, M, are limited. Hence, in order to correctly evaluate error ε 2 , a separation in two sets of the available realizations is generally performed:

• the first set, {X(θ 1 ), . . . , X(θ ν * )}, is a learning set, on which the definition of B (M) is based, • the second set, {X(θ ν * +1 ), . . . , X(θ ν )}, is a validation set, on which the computation of ε 2 ν-ν * is achieved to evaluate ε 2 . With such a method, it can be noticed that the higher ν * , the less precise the evaluation of ε 2 . This limits strongly the scope of such approaches when the number of available realizations ν is small, compared to the number of functions that are needed to characterize X.

In this context, the application of the jackknife theory [START_REF] Miller | The jackknife -a review[END_REF][START_REF] Quan | The prediction sum of squares as a general measure for regression diagnostics[END_REF][START_REF] Blatman | Adaptive sparse polynomial chaos expansion based on least angle regression[END_REF] to the evaluation of projection errors appears to be a good method to evaluate precise approximations for ε 2 (B (M) ) when all realizations of X have been used to compute ε 2 (B (M) ). To this end, for all 1 ≤ i ≤ ν, we denote by X (-i) = {X(θ n ), 1 ≤ n ≤ ν, n = i} the set gathering all the available realizations of X but the i th one, by B (M) (X (-i) ) the M -dimensional family that has been computed from this set X (-i) (the family that stems from the Gram-Schmidt orthogonalization defined by Eq. (B.1) is an example of such a family), and by X (-i) the projection of X on B (M) (X (-i) ). Then, the two following hypotheses are assumed:

1. First, it is supposed that error ε 2 B (M) (X (-i) ) decreases when ν increases. 2. Then, given two sets X (-i) and X (-j) with 1 ≤ i = j ≤ ν, it is assumed that:

P X (-i) ≈ P X (-j) , (2.21)
where P X (-i) and P X (-j) are the distributions of X (-i) and X (-j) respectively. In other words, the first hypothesis means that the modeling errors are expected to decrease when the information is increasing, whereas the second hypothesis asks for the application that compute the projection basis from a limited set of realizations to give robust results. It can therefore be deduced from these two hypotheses that, for ν sufficiently high:

ε 2 (B (M) ) ≈ ε 2 LOO (B (M) ), (2.22)
ε 2 LOO (B (M) ) def = 1 ν ν n=1 X(θ n ) -X (-n) (θ n ), X(θ n ) -X (-n) (θ n ) .
(2.23)

According to the central limit theorem, estimator ε 2 LOO , which is called Leaveone-out (LOO) error, converges to ε 2 at the convergence rate of 1/ √ ν. Finally, the L 2 -error, ε 2 , in the optimization problem defined by Eq. (2.18), can be replaced by its LOO estimator, such that the M -dimensional optimal projection family, K (M) , can be approximated by the following optimization problem:

     K (M) ≈ B (M) A(α * ) , α * = arg min α∈[0,1] ε 2 LOO (B (M) A(α) ) . (2.24)
As in Section 2.3, under the hypothesis that the optimization problem defined by Eq. (2.18) is convex, such that the dichotomy-based algorithm defined by Eq. (A.1) can be used, the solving of N cost (= 2 + p) × ν eigenvalue problems defined by Eq. (2.9), is required to solve the optimization problem defined by Eq. (2.24) with a precision of T for α * . However, keeping in mind that computing the LOO error associated with any projection family amounts to solving a series of slightly modified eigenvalue problems, it was noticed that iterative methods, such as the subspace iteration methods (see [START_REF] Bathe | Numerical Methods in Finite Element Analysis[END_REF] for further details about these methods), can be used to reduce drastically the computational time.

Applications.

In order to illustrate the benefits that stem from the optimization problem defined by Eq. (2.24), two applications based on simulated data are introduced in this section. The first application deals with a non-Gaussian stochastic process characterized by a rather smooth correlation function. This application aims thus at justifying the relevance of the Leave-one-out error and at emphasizing the difficulties of the classical KL expansion to identify optimal basis when the number of available realizations is low, while the generalized KL expansion, characterized by Eq. (2.24), gives promising results. The second application then shows to what extent such an approach can also be applied to processes characterized by non-smooth correlation functions.

3.1. Relevance of the generalized KL expansion to approximate stochastic processes characterized by smooth covariance functions. The objective of this section is to illustrate the results shown in Section 3. To this end, we introduce the following quantities:

• Ω = [0, 1].

• R XX is a particular element of S(R), such that for all κ in [0, 1] and all (s, s ′ ) in Ω × Ω:

R XX (s, s ′ ) = P (κ, s, s ′ ) exp (-5|s -s ′ |) cos (10π|s -s ′ |) , P (κ, s, s ′ ) = cos (2πκs) cos (2πκs ′ ) , (3.1)
• {k m , m ≥ 1} and {λ m , m ≥ 1} are the respective sets of eigenfunctions and eigenvalues associated with R XX . In this work, the eigenvalue problem defined by Eq. (2.9) with A = R XX is solved using a Galerkin-type approximation at the spatial step h = 2.5 × 10 -3 ≪ 1. For the sake of simplicity, no distinction will be made between R XX and its Galerkin approximation at the spatial step h, which amounts to assuming that for (s, s ′ ) in Ω × Ω:

R XX (s, s ′ ) = N m=1 λ m k m (s)k m (s ′ ), N = 1/h + 1 = 401. (3.2)
• C = (C 1 , . . . , C N ) is a N -dimension non-Gaussian random vector, whose components are uncorrelated but highly dependent, such that:

-0.5 0.5

κ = 0 R XX Abscissa s Abscissa s' 0 1 0 1 -0.5 0.5 κ = 0.4 R XX Abscissa s Abscissa s' 0 1 0 1 -0.5 0.5 κ = 1 R XX Abscissa s Abscissa s' 0 1 0 1 Figure 3.1: Representation of the covariance function (s, s ′ ) → R XX (s, s ′ ) for three particular values of κ in [0, 1]. C = 0≤α1+•••+α10≤4 c α1,••• ,α10 ξ α1 1 × • • • × ξ α10 10 , (3.3)
where ξ = (ξ 1 , . . . , ξ 10 ) is a 10-dimensional random vector, whose components are independent and uniformly distributed on [-1, 1], and c α1,••• ,α10 are deterministic coefficients, whose components are uniformly and independently chosen between -1 and 1 before being normalized in order to guarantee that

E [C i C j ] = δ ij .
Finally, under such a formalism, we define X as the stochastic process defined by:

X(s) = N m=1 C m k m , s ∈ Ω. (3.4)
By construction, X ∈ P(Ω) is a centered, non-Gaussian stochastic process, whose covariance function is equal to R XX . From each independent realization, ξ(θ), of ξ, we have therefore access to an independent realization, X(θ), of X. As an illustration, for κ = 1, four independent realizations of X are represented in Figure 3.2. Function P has moreover been introduced as a non-stationary perturbation such that the higher κ, the less mean-square stationary X. The influence of P can indeed be seen in Figure 3.1, where function R XX is represented for three values of κ. Two sets, X exp = {X(θ 1 ), . . . , X(θ ν )} and X valid = {XΘ 1 ), . . . , X(Θ ν valid )}, of independent realizations of X then are generated. Set X exp represents the available information for X, whereas X valid is the validation set. According to Section 2.4.1, this set will only be used to evaluate the projection error corresponding to any projection family B (M) in H M . In order to lighten the notations, no distinction will be made in the following between error ε 2 (B (M) ) and its empirical estimator, ε 2 ν valid (B (M) ), calculated from the realizations gathered in X valid .

3.1.1. Improvement of the projection basis with respect to the available information. The number of available realizations, ν, is now supposed to be in the set {0, 20, 50, 100, 400}. The case ν = 0 corresponds to a limit case when no realization of X is available. For the other cases, the empirical estimator of covariance function R XX , which is denoted by R(ν), such that:

X(θ 1 ) X(θ 2 ) X(θ 3 ) X(θ 4 )
R(ν, s, s ′ ) = 1 ν ν n=1 X(θ n , s)X(θ n , s ′ ), (s, s ′ ) ∈ Ω × Ω, (3.5)
is compared in Figure 3.3, in the case κ = 1. In these figures, it can be verified that the higher ν, the more relevant R(ν). For all ν > 0, R(ν) is in S(R) and we denote by

B (M) R(ν)
the set gathering the M eigenfunctions of R(ν) of highest eigenvalues. In particular, in the case ν = 0,

B (M) R(ν=0)
corresponds to any M -dimensional set of orthonormal functions of H. In Figure 3.4 are thus compared the evolutions of the error functions, ε 2 (B (M) R(ν) ), with respect to M , for different values of ν and κ. First, it can be noticed in these figures that ε 2 (B (M) R(ν=0) ) decreases linearly with respect to M , which means that the relevance of each element of B (M) R(ν=0) to describe X is approximatively the same. This is a direct and natural consequence of the fact that all these elements have been defined without information on X. Then, two phases can clearly be identified in the evolution of ε 2 (B are based on the available realizations of X, whereas the M -ν last elements are not.

3.1.2. Relevance of the Leave-one-out error. As presented in Section 2.4, when the assessment set, X valid , is not available, which is the general case, the Leaveone-out error allows us to evaluate the projection error from the only set X exp . For ν = 10, 20, 50, 200 and κ = 0, 0.4, 1, the relevance of the Leave-one-out error can also be noticed in is particularly relevant to characterize X when the number of available realizations, ν, is higher than M . This section aims at illustrating the benefits of the approach introduced in Section 2.3, in cases when M ≫ ν elements are needed for the projection error to be sufficiently small. To this end, the evolution of error ε 2 (B (M) A(α) ) with respect to α is therefore represented in Figure 3.5, for several values of κ and M , when only ν = 50 independent realizations of X are supposed to be available. On these figures, the optimal value of α to characterize X, α * , which is solution of the optimization problem defined by Eq. (2.18), has been added. On the first hand, if X is mean-square stationary, that is to say if κ = 0, we can thus verify in this figure that the choice α = 0 leads us particularly adapted projection basis for X. On the other hand, if X is not mean-square stationary, significant reductions for ε 2 (B (M) A(α) ) can also be obtained by searching the optimal value of α in [0, 1], especially when M = 100 > ν = 50.

More precisely, the relevance of B

(M) R(ν)
and B (M) A(α * ) is compared in Figure 3.6, for representative values of κ and ν (improvements of the same order of magnitude were obtained for the other values of (κ, ν) in {0, 0.4, 1} × {20, 50, 100, 400}). The optimal value of the projection error, ε 2 (K (M) ), has been added in these figures as a limit state. In each case, it can thus be seen that ε 2 (B the projection error associated with a particular projection family is therefore based on the computation of ν = 50 slightly modified eigenvalue problems. Moreover, as in Section 3.1.3, the dichotomy-based algorithm presented in Section 2.3, with a precision T = 5%, was once again used to identify accurate approximations of optimal value α * . Hence, for all 1 ≤ M ≤ 100, ν + 1 = 51 and N cost = 7(ν + 1) = 357 solvings of the eigenvalue problem defined by Eq. (2.9) were needed to evaluate ε 2 LOO (B

(M) R(ν)
) and ε 2 LOO (B (M)

A(α * ) ) respectively. Finally, when the available information about X is characterized by a finite set of independent realizations, the relevance of the proposed method is illustrated in Figure 3.9. Indeed, by comparing the evolutions of errors ε 2 LOO (B

(M) R(ν) ), ε 2 LOO (B (M)
A(α * ) ) and ε 2 (K (M) ) with respect to M , it appears that interesting improvements can be brought by the coupling of a Leave-one-out procedure and the solving of an original approximation of the covariance operator to identify accurate reduced basis for X.

Conclusions.

For the last fifty years, the increasing computational power has encouraged many scientific fields to take into account stochastic processes in their modeling. The development of reduced basis that could condense at best the statistical properties of these stochastic processes is therefore of great interest. In most of these applications, the knowledge of these stochastic processes is however limited to a finite set of independent realizations. In this context, using two examples, this paper emphasized the efficiency of a method based on an adaptation of the Karhunen-Loève expansion, in order to construct optimal basis from a relatively small set of independent realizations. This method defined first an original optimization problem, and secondly, required a way to a posteriori evaluate projection errors. Finally, when interested in studying complex systems that are excited by stochastic processes that are only known through a set of limited independent realizations, the proposed method opens new opportunities to optimize the projection basis with respect to the available information.

                                 
Initialize: α 1 = 0, α 2 = 1. while |α 1 -α 2 | > T :

Figure 3 . 2 :

 32 Figure 3.2: Representation of four independent realizations of X, for κ = 1.

  ) ) with respect to M , for ν = 10, 20, 50, 200: the decrease of ε 2 (B(M) R(ν)) is indeed much faster for M ≤ ν than for M > ν, where a quasi-linear decrease is found again. As explained in Section 2.2, this behavior can be justified by the fact that the ν first elements of B (M) R(ν)

Figure 3 . 3 :

 33 Figure 3.3: Empirical estimators R(ν) for four values of ν.

Figure 3 . 7 :Figure 3 . 8 :

 3738 Figure 3.7: Representation of the chosen Ornstein-Uhlenbeck covariance function, R XX , and of its empirical approximation, R(ν), for ν = 50.

ε 2 ( 6 Figure 3 . 9 :

 2639 Figure 3.9: Improvement of the projection basis for the characterization of Ornstein-Uhlenbeck process X.

Acknowledgments. This work was supported by SNCF (Innovation and Research Department).

Appendix A. Dichotomy-based identification algorithm.

ν ≪ M . Hence, for these four choices for (κ, ν), this figure underlines the great benefits brought by the formulation defined by Eq. (2.18). For instance, for ν = 20, κ = 0.4 and M = 50, we get:

Indeed, whereas the rank of R(ν = 20) is 20 < M = 50, the rank of A(α * ) is by construction much higher than ν. Therefore, more elements of B

A(α * ) are based on the knowledge of X than the elements of

, which explains such an improvement of the projection basis, even if X is not mean-square stationary.

For all 1 ≤ M ≤ 200, and all (κ, ν) in {0, 0.4, 1} × {20, 50, 100, 400}, the computation of optimal value α * was based on the dichotomy-based algorithm presented in Section 2.3, with a chosen precision T = 5%. Hence, whereas A(α) ), with respect to α, for M in {10, 15, 30, 50, 100} and κ in {0, 0.4, 1}, when the available information is characterized by a set of ν = 50 independent realizations of X.

the direct solving of the eigenvalue problem defined by Eq. (2.9), with A = R(ν), the computation of B 3.2. Relevance of the generalized KL expansion to characterize stochastic processes characterized by non-smooth correlation functions. In order to underline that the proposed approach can be applied to any type of correlation function, the Ornstein-Uhlenbeck case [START_REF] Uhlenbeck | On the theory of the brownian motion[END_REF], which is characterized by a non-smooth correlation function, is tested in this section. To this end, let X = {X(s), s ∈ [0, 1]} be the Gaussian stochastic process, whose mean function µ and covariance function R XX are given by, for (s, s ′ ) in Ω × Ω:

Abscissa s (κ = 0, ν = 50)

Abscissa s (κ = 1, ν = 100) (κ = 0.4, ν = 400)

Abscissa s

Abscissa s

where a = 1, b = 10 -2 , and c = 0.4 are deterministic parameters.

In the same manner as in Section 3.1, it is now supposed that the maximal information about X is characterized by a set of only ν = 50 independent realizations, which are once again denoted by {X(θ 1 ), . . . , X(θ ν )}. As an illustration of such a stochastic process, Figure 3.7 compares the evolutions of the covariance function of X, R XX , and of its empirical estimator, R(ν), whereas Figure 3.8 shows four particular realizations of X. Using the same notation as in Section 3.1, we define:

• K (M) and B (M) R(ν) as the sets gathering the M eigenfunctions of highest eigenvalues associated with R XX and R(ν) respectively; • α * and A(α * ) as the solution of the optimization problem defined by Eq.

(2.18) and the corresponding function in S(R) respectively. In this section, it is supposed that no validation set is available to evaluate the relevance of any projection family to characterize X. The Leave-one-out procedure, described in Section 2.4 has thus been applied. For all 1 ≤ M ≤ 100, the evaluation of