C Casenave
email: casenave@laas.fr

G Montseny
email: montseny@laas.fr

Identification of Nonlinear Volterra Models by means of Diffusive Representation

Keywords: Volterra model, non rational operator, state realization, diffusive representation

We introduce a new identification method for nonlinear Volterra models of the form Hx = f (u, x) with H a causal convolution operator. It is mainly based on a suitable parameterization of H deduced from the so-called diffusive representation, devoted to state representations of integral operators. Following this approach, the complex dynamic nature of H can be summarized by a few numerical parameters on which the identification of the dynamic part of the model will focus. For illustration, we implement this method on a concrete numerical example.

INTRODUCTION

We introduce an identification method for nonlinear Volterra models of the form:

H(∂ t )x = f (u, x),
(1) where u(t), x(t) ∈ R, t > 0, f is a sufficiently regular function defined on R 2 and H(∂ t ) is a causal convolution operator defined on a suitable space of continuous functions with support in R + t . We suppose that problem (1) is well-posed in the sense of existence, uniqueness and continuous dependency on the input u of the solution x. The operator H(∂ t ) is also supposed to be invertible, so (1) can be rewritten under the standard form: x(t) = t 0 k(ts) f (u(s), x(s)) ds ∀t > 0, where k is the impulse response of operator H(∂ t ) -1 .

Such models are frequently encountered in various domains:

• thermic phenomena (we will consider in section 4 the combustion model elaborated in [START_REF] Joulin | Point-source initiation of lean spherical flames of light reactants: an asymptotic theory[END_REF]), • electrical engineering [START_REF] Rumeau | Behavior modeling of a CaCu 3 Ti 4 O 12 ceramic for capacitor applications[END_REF]),

• linear SISO differential systems with nonlinear feedback:

Ẋ = AX + Bf (u, x) x = CX, with H(p) -1 = C(pI -A) -1 B,
• SISO partial differential systems on R + * t ×Ω z ⊂ R n+1 of the form:

∂ t ϕ = A(z, ∇)ϕ + f (u, T ϕ dz), ϕ 0 = 0, x = T ϕ dz, which can be rewritten under the synthetic equation t 0 k(ts) f (u(s), x(s)) ds = x(t) where the function k is the impulse response associated with the operator f → x, that is the solution of: φ = A(z, ∇)φ + δ, φ 0 = 0 k = T φ dz. • etc.

The identification problem under consideration in the sequel is to build estimations (if possible optimal) of H(∂ t ) and/or f from (noised) data x * = x+ v where v designates some additive measurement noise and x is generated by an experimental process driven by a known input u. The main difficulty when identifying such models results from the coupling between the dynamic operator H(∂ t ) and the (static) function f via equation ( 1). This difficulty is strengthened when operator H(∂ t ) is non rational (which is in general the case), when f is singular and when there exists some dynamic bifurcations (see section 4).

The proposed identification method is based on a suitable parameterization of H(∂ t ) deduced from the so-called diffusive representation [START_REF] Montseny | Représentation Diffusive[END_REF]), devoted to state realizations of integral causal operators. Following this approach, the complex dynamic nature of H(∂ t ) is in some sense summarized by a few numerical parameters on which the identification problem will focus.

The paper is organized as follows. In section 2, we briefly present a simplified version of the diffusive representation. In section 3, we describe the identification method in a general framework and we give some indications for numerical implementation. In section 4, we finally implement this method on a concrete example, which allows us to comment the quantitative relevance of the method.

DIFFUSIVE FORMULATION OF CAUSAL CONVOLUTION OPERATORS

A complete statement of diffusive representation will be found in [START_REF] Montseny | Représentation Diffusive[END_REF]. Various applications and questions relating to this approach will be found for example in [START_REF] Audounet | An asymptotic fractional differential model of spherical flame[END_REF], [START_REF] Ph | Fractional Brownian Motion and the Markov Property[END_REF], [START_REF] Casenave | Optimal identification of delay-diffusive operators and application to the acoustic impedance of absorbent materials[END_REF],

Casenave [2008], [START_REF] Degerli | Bandlimited 1/f α noise source[END_REF], [START_REF] Garcia | Identification of the dynamics of a lead acid battery by a diffusive model[END_REF], [START_REF] Lenczner | Diffusive realization of operators solutions of certain operational partial differential equations[END_REF], [START_REF] Levadoux | Diffusive formulation of the impedance operator on circular boundary for 2D wave equation[END_REF], [START_REF] Montseny | Diffusive representation for operators involving delays[END_REF], Mouyon [2002], [START_REF] Rumeau | Behavior modeling of a CaCu 3 Ti 4 O 12 ceramic for capacitor applications[END_REF].

Mathematical framework

We consider a causal convolution operator defined, on any continuous function u : R + → R, by

u → t 0 k(t -s) u(s) ds.
(2)

We denote K the Laplace transform of k and K(∂ t ) the convolution operator defined by (2).

Let u t (s) = 1 ]-∞,t] (s) u(s) the restriction of u to its past and u t (s) = u t (ts) the so-called "history" of u. From causality of K(∂ t ), we deduce: [K(∂ t )(uu t )](t) = 0 for all t; then, we have for any continuous function u:

[K(∂ t ) u](t)= L -1 (K Lu) (t)= L -1 K Lu t (t
∂ t Ψ(t, p) = p Ψ(t, p) + u, t > 0, Ψ(0, p) = 0. 2. For any b 0, [K(∂ t ) u] (t) = 1 2iπ b+i∞ b-i∞ K(p) Ψ u (t, p) dp.
(5)

We denote Ω the holomorphic domain of K (after analytic continuation). Let γ a closed 1 simple arc in C -; we denote Ω + γ the exterior domain defined by γ, and Ω - γ the complementary of Ω + γ . By use of standard techniques (Cauchy theorem, Jordan lemma), it can be shown:

Lemma 2. For γ ⊂ Ω such that K is holomorphic in Ω + γ , if K(p) → 0 when p → ∞ in Ω + γ , then: [K(∂ t ) u] (t) = 1 2iπ γ K(p) Ψ u (t, p) dp, ( 6
)
where γ is any closed simple arc in Ω + γ such that γ ⊂ Ω - γ .

We now suppose that γ, γ are defined by functions of the Sobolev space 2 W 1,∞ loc (R; C), also denoted γ, γ and such that:

γ(0) = 0. (7) We use the convenient notation μ, ψ = μ ψ dξ; in particular, when μ is atomic that is μ = k a k δ ξ k , we have: μ, ψ = k a k ψ(ξ k ).
Under hypothesis of lemma 2, we have [START_REF] Montseny | Représentation Diffusive[END_REF]): Theorem 3. If the possible singularities of K on γ are simple poles or branching points such that |K • γ| is locally integrable in their neighborhood, then:

1. with μ = γ 2iπ K • γ and ψ(t, .) = Ψ u (t, .) • γ: [K(∂ t ) u] (t) = μ, ψ(t, .) ; ( 8 ) 1 Possibly at infinity 2 W 1,∞ loc (R; C) is the topological space of measurable functions f : R → C such that f, f ∈ L ∞
loc (that is f and f are locally bounded in the almost everywhere sense).

2. with3 γn → γ in W 1,∞ loc and μ = γ 2iπ lim K • γn in the sense of measures:

[K(∂ t ) u] (t) = μ, ψ(t, .) , (9) where ψ(t, ξ) is solution of the following evolution problem on (t, ξ) ∈ R * + ×R (of diffusive type):

∂ t ψ(t, ξ) = γ(ξ) ψ(t, ξ) + u(t), ψ(0, ξ) = 0.
(10) Definition 4. The measure μ defined in theorem 3 is called γ-symbol of operator K(∂ t ). The function ψ solution of ( 10) is called the γ-representation of u.

Note in particular that thanks to (9), the Dirac measure δ is clearly a γ-symbol of the operator u

→ t 0 u(s) ds, denoted ∂ -1 t . We indeed have (∂ -1 t u)(t) = δ, ψ(t, .) = ψ(t, 0), with ∂ t ψ(t, 0) = u, ψ(0, 0) = 0.
Beyond the measure framework, the general space of γsymbols is a quotient space of distributions, denoted Δ γ ; it is the topological dual of the space Δ γ ψ(t, .) [START_REF] Montseny | Représentation Diffusive[END_REF]). The composition product of operators has an equivalent in Δ γ , denoted γ or simply : if μ and ν are respective γ-symbols of H(∂ t ) and

K(∂ t ), then μ γ ν is a γ-symbol of H(∂ t ) • K(∂ t ).
Note that the product γ is inner, commutative and continuous 4 in Δ γ and so (Δ γ , γ ) is an algebra (of γ-symbols) isomorphic to a commutative algebra of causal convolution operators. Formulation (10,9) can be extended to operators of the form

K(∂ t ) • ∂ n t where K(∂ t ) admits a γ-symbol in Δ γ . We have (formally): [K(∂ t ) • ∂ n t u](t) = μ, ∂ n t ψ(t, .
) , (11) with ψ(t, ξ) solution of ( 10) and μ the γ-symbol of K(∂ t ). In the particular case where n = 1, (11) becomes:

[K(∂ t ) • ∂ t u](t) = μ, γ ψ(t, .) + u(t) .
(12) In the same way, Δ γ can be extended to an algebra denoted Σ γ whose elements are the γ-symbols of operators of the form

K(∂ t ) • ∂ n t where K(∂ t ) is associated to a γ-symbol μ ∈ Δ γ and n ∈ N. γ-symbols ν of K(∂ t ) • ∂ n t are characterized by the relation: μ = ν δ n where δ n = δ δ ... δ n times ∈ Δ γ is a γ-symbol of ∂ -n t .
The inversion of γ-symbols cannot be defined in Δ γ because this algebra is not unitary; this operation is nevertheless well-defined in Σ

γ . If μ ∈ Σ γ is a γ-symbol of K(∂ t ) such that K(∂ t ) -1 • ∂ -n t has a γ-symbol ν ∈ Δ γ , then ν = μ -1 δ n and we have: [K(∂ t ) -1 u](t) = μ -1 δ n , ∂ n t ψ(t, .
) , (13) with ψ solution of (10). Note in particular that operator ∂ -1 t defined above is the unique inverse of the derivative operator ∂ t .

About numerical approximations

We only give a few indications. More details can be found in the references cited above.

The state equation ( 10) is infinitedimensional. To get numerical approximations, we consider M L a sequence of Ldimensional spaces of atomic measures on suitable meshes {ξ L l } l=1:L on the variable ξ; L-dimensional approximations μ L of the γ-symbol μ ∈ Δ γ are then defined in the sense of atomic measures, that is:

μ L = L l=1 μ L l δ ξ L l , μ L l ∈ C. (14) If ∪ L M L is dense in the topological space Δ γ (that is, concretely, if ∪ L {ξ L l } is dense in R)
, then we can have [START_REF] Montseny | Représentation Diffusive[END_REF]):

μ L , ψ -→ L→+∞ μ, ψ ∀ψ ∈ Δ γ ;
(15) so, we have the following L-dimensional approximate state formulation of K(∂ t ) (with γ-symbol μ):

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ ∂ t ψ(t, ξ L l ) = γ(ξ L l ) ψ(t, ξ L l ) + u(t), l = 1 : L, ψ(0, ξ L l ) = 0 [K(∂ t ) u](t) L l=1 μ L l ψ(t, ξ L l ).
Note that an approximate state formulation of operator ∂ t • K(∂ t ) can be easily deduced under the form:

∂ t • K(∂ t ) u l γ(ξ L l ) μ L l ψ(., ξ L l ) + l μ L l u.
One of the properties of the approach presented above is that most of non rational operators encountered in practice can be closely approximate with small L (see for example [START_REF] Montseny | Simple approach to approximation and dynamical realization of pseudodifferential time-operators such as fractional ones[END_REF]). In the context of identification of Volterra models, this will be a great advantage because only a few numerical parameters μ L l will have to be identified from experimental data, while the property (15) will ensure the well-posedness and the robustness of the problem as soon as the operator to be identified admits a γ-symbol in Σ γ .

IDENTIFICATION OF VOLTERRA MODELS

In this section, we focus on the problem of identification of Volterra models of the form (1). We simply introduce the principle of the approach and give some indications about the numerical aspects of the problem. The statement is essentially formal. More details will be given in a further work.

Principle

In the sequel, we suppose for simplicity that the dynamic operators

H(∂ t ) -1 and H(∂ t )•∂ -1 t both admit a γ-symbol in Δ γ ; we denote μ any γ-symbol of H(∂ t ) • ∂ -1
t . For any fixed u, equation (1) expresses the balance between two trajectories obtained from x, relating respectively to the linear dynamic operator H(∂ t ) and the (nonlinear) static operator defined from function f by5 : [f (u, x)](t) = f (u(t), x(t)) ∀t. So, from the point of view of trajectories (and if the numerical evaluation of H(∂ t )x is of reasonable cost), we can remark that the operators H(∂ t ) and f (u, .) play a comparable role in expression (1). The proposed identification method then consists in parameterizing both operator H(∂ t ) by means of its γ-symbol and function f by means of a suitable functions basis. So, we will get an equivalent problem in which the unknown parameters are linearly dependent on the data and of reasonable dimension under numerical approximation (thanks to the properties of diffusive representation mentioned in the previous section).

Given a suitable γ and denoting ψ x the γ-representation of x, according to (12), we consider the γ-realization of H(∂ t ) defined by: H(∂ t ) x = μ, γ ψ x + μ, 1 x. By denoting A x the linear operator defined on any γsymbol μ by: ∀t > 0, [A x μ](t) = μ, γ ψ x (t, .) + μ, 1 x(t), (16) we then get:

H(∂ t ) x = A x μ.
Now suppose that: f = f + f where f is known a priori and f is unknown (to be identified from experimental data). We consider a topological basis {g i ⊗ k j } i,j=1:+∞ of a tensorial product of Hilbert spaces to which belongs the function f ; we then have:

f = f + i,j a ij g i ⊗ k j ,
where the unknown is now the set of real parameters a := (a ij ). If the trajectories u and x are known (from measurements), so is it for the trajectories g i (u), k j (x) and f (u, x); equation (1) can then be equivalently expressed under the linear form:

A x μ - i,j g i (u) k j (x) a ij = f (u, x). ( 17
)
By denoting

G u,x : (μ, a) → A x μ - i,j g i (u) k j (x) a ij , (18)
b = f (u, x), Ξ = (μ, a) and given two suitable Hilbert spaces E and F , the identification problem of model (1) from data (u, x * ), x * = x + v where v designates some additive measurement noise, can be expressed as: min

Ξ∈E G u,x * Ξ -b * 2 F . ( 19
)
Thanks to the linearity of operator G u,x * , the solution of ( 19) is obtained by orthogonal projection:

Ξ * = G † u,x * b * , ( 20
)
where G † u,x * denotes the pseudo-inverse of G u,x * (Ben-Israel [2003]). In the sense of the hilbertian norm of F , the estimation Ξ * of Ξ is then optimal.

The physical system under consideration can then be described by the model H * (∂ t )x = f * (u, x) where H * and f * are deduced from the identified parameters μ * , a * . Furthermore, we deduce from the equivalent expression x = H * (∂ t ) -1 f * (u, x) that the following input-output state realization is available (up to the computation of δ (μ * ) -1 which can be numerically performed) for simulation or control purposes:

∂ t ψ = γ ψ + f * u, δ (μ * ) -1 , ψ , ψ(0, .) = 0 x = δ (μ * ) -1 , ψ .
(21)

Numerical formulation

We consider now a finitedimensional subspace E k of E.

After suitable time discretization and approximation of μ such as described in section 2, we deduce from (16,18):

∀n = 1 : N, [G u,x (μ L , a)](t n ) = L l=1 μ L l γ(ξ L l ) ψ x (t n , ξ L l ) + x(t n ) - I i=1 J j=1 g i (u(t n )) k j (x(t n )) a ij ;
then the operator G u,x can be expressed by means of a matrix G u,x ∈ M N,L+I×J whose pseudo-inverse is classically given by 6 (with N L + I × J):

G † u,x = (G * u,x G u,x + I) -1 G * u,
x . Remark 5. Note that thanks to the state realization (10) of ψ x and the static nature of operators g i , k j , recursive formulations of (20) can be established, which allows to treat real-time identification or even pursuit problems.

APPLICATION TO A MODEL OF FLAME

In this section, we illustrate the identification method introduced above by implementing it on data (u, x * ) elaborated from numerical simulations of a complex dynamic phenomenon studied in [START_REF] Joulin | Point-source initiation of lean spherical flames of light reactants: an asymptotic theory[END_REF], [START_REF] Audounet | An asymptotic fractional differential model of spherical flame[END_REF].

The model under consideration

In [START_REF] Joulin | Point-source initiation of lean spherical flames of light reactants: an asymptotic theory[END_REF], Joulin elaborated a Volterra model to describe, in suitable thermodynamic conditions, the evolution of a spherical flame initiated by a source at point 0 in a mixture of reactive species. Under some reasonable physical hypothesis, such a phenomenon can be described by a system of two partial differential equations relating to the temperature and the mass density of the mixture. By considering the reactive zone as a thin sheet located on a sphere with radius x(t), Joulin has established that when the flame is developing in free space, x is solution of the following nonlinear singular Abel-Volterra equation 7 (u(t) designates the source strength at time t):

x(t) t 0 ẋ(s) ds π(t -s) = 2 x(t) ln x(t) + 2 u(t) ∀t > 0, ( 22
)
with the additional conditions: x(0) = 0, u 0, x 0 (whose physical interpretation is obvious). By denoting

H(∂ t ) the convolution operator 8 x → t 0 ẋ(s) ds √ π(t-s) and f (u, x) := 2 ln x + 2 u
x , ( 22) can be formally rewritten under the form (1). It has been shown in [START_REF] Audounet | An asymptotic fractional differential model of spherical flame[END_REF] that the evolution problem ( 22) is well-posed, that is the solution x exists, is unique and depends continuously on u.

In real conditions however, various perturbations are involved in the evolution of x (due for example to the loss of spatial symmetries), and both the convolution operator 6 G * is the dual of G with respect to the scalar product under consideration. As usual, 0 is devoted to numerical reconditioning. 7 Here adimensional for simplicity.

8 Note that in this ideal case,

H(∂t) = ∂ 1 2 t .
H(∂ t ) and the function f will be more or less far from the ideal ones. So, an identification process can be justified if accuracy of the model is required.

We will consider in the sequel the problem of identification of (1) from data (u, x * ) obtained from highly accurate numerical simulations of (22) (not described here).

On the dynamic behavior of the flame radius

It has been shown in [START_REF] Audounet | An asymptotic fractional differential model of spherical flame[END_REF] that there exists a threshold relating to the power of the source u, beyond which the flame is developing whereas a quenching occurs below. In that sense, this evolution phenomenon is essentially unstable with two qualitatively different behaviors.

Because of the hereditary nature of the problem, it is difficult to know the value of this threshold which must be evaluated on the basis of numerical simulations. In figures 1 and 2, the flame is either quenching or developing, the source function being given, as in [START_REF] Audounet | An asymptotic fractional differential model of spherical flame[END_REF], by: Such dynamic behaviors generate ill-conditioned identification problems due to the sensitivity of the solutions x of (1) with respect to the source power. As a consequence, this model can be viewed as a significant test for the proposed identification method.

u(t) = E t 0.3 (1 -t) 1 [0,1] (t). ( 23
) 0 0.5 1 1.5 2 2.5 3

Parameterization of the problem

The flame model under consideration can be written under the form (1) with, on the one hand:

f (u, x) = 2 k(x) + 2 u x
where the function k is to be identified, and on the other hand, the operator H(∂ t ) which is also to be identified with the assumption that H(∂ t ) • ∂ -1 t admits a γ-symbol μ ∈ Δ γ . As described in section 3, we search an estimation of f of the (simpler) form:

f (u, x) = 2 u x + 2 j a j k j (x),
where k j are (known) basis functions by means of which the function f (u, x) -2 u x can be approximate, for example in the sense of L 2 (x min , x max ). We simply consider power functions:

k j (x) = x j-1 .

According to section 3, the vector b * is then defined by:

b * n = 2 u(t n ) x * (t n ) ,
and the matrix G u,x * by:

G u,x * = [Ψ | K], Ψ nl = γ(ξ L l ) ψ x * (t n , ξ L l ) + x * (t n ), K nj = 2 (x * (t n )) j-1 .

Numerical identification results and comments

Recall that the problem consists in identifying the model (1) with the following functions to be estimated from numerical experimental data via the parameterizations μ, a defined above:

H(p) = p 1 2 , f(u, x) = 2 ln(x) + 2 u x .
In a first time, only the operator H(∂ t ) is identified, while the function f is suppose to be known. Then, H(∂ t ) and f are identified together.

The data are composed of 4 trajectories (u, x * ) (see figure 3) associated with 4 different sources of the form9 (23) with E = 1.5, 1.738, 1.8 and 5.0 respectively. The time step Δt has been taken equal to 5.10 -5 and the maximal final time is equal to 10. The measured output is x * = x + ε v where v is a unity gaussian white noise and ε = 10 -3 .

(ξ) = ξ e i 3π 4 if ξ > 0 -ξ e -i 3π 4 if ξ < 0.
First, we suppose that f is known and we identify H(∂ t ) via the γ-symbol μ. We use L = 20 values ξ L l geometrically spaced to cover 4 decades from 10 -1 to 10 3 . The frequency response of the so-identified operator H(∂ t ) is given in figure 4. In the accessible frequency band, the identification is good. To complete the validation, we can compare in figures 5 and 6 the respective responses of the identified model and of the exact one for inputs u different from those used for identification. In conclusion, when f is known, the identified model closely behaves like the exact one.

(t) = 4.0 t 0.3 (1 -t) 1 [0,1] (t).
We can add that results with comparable quality are obtained when f is identified only (i.e. H is supposed to be known).

We now consider the problem of identification of operator H(∂ t ) and function f together (both are supposed unknown). We use L = 20 values ξ L l geometrically spaced to cover 4 decades from 10 0 to 10 4 and the estimation of f is searched as: The so-identified operator H(∂ t ) and function f are given in figures 7 and 8. The identified H and f remain close to the exact ones, which agrees with the theoretical analysis. However, due to the increasing number of parameters to be identified, the ill-conditioning of the problem and the dynamic poverty of the data (u, x * ), the estimation error on H and f is greater than previously. It follows that the validation by simulation of the identified model reveals itself more delicate, in particular in the case of inputs with power close to the bifurcation threshold, which leads to the conclusion that more data will be necessary to improve the identification quality.

f (u, x) = 2 u x + 2 9 j=0 a j x j .

CONCLUSION

Relating to the first attempt to identify nonlinear Volterra models by use of the method presented in section 3, the numerical results obtained in section 4 can be considered as positive, from both points of view of implementation simplicity and accuracy. Several questions must be studied in order to improve such results. For example, among the most significant, the involved hilbertian norms should be judiciously chosen and adapted to the specific properties of the class of models under consideration; indeed, this choice is crucial in terms of sensitivity with respect to perturbations of any nature. It can also be shown that the measurement noise induces some estimation bias which should be significantly reduced by appropriate treatments. All these questions are currently under study and will be discussed in a deepened paper.

Fig. 2 .

 2 Fig. 1. Source (---) and radius of the flame (-) with E = 1.7390.

Fig. 3 .

 3 Fig. 3. The trajectories x used for identification. The function γ is chosen as γ(ξ) = ξ e i 3π 4 if ξ > 0 -ξ e -i 3π 4 if ξ < 0.

Fig. 5 .

 5 Fig. 4. Exact (--) and identified (-) frequency responses H(iω) when f is known.

Fig. 6 .

 6 Fig.6. Evolution of x for the exact (--) and identified (-) models, with u(t) = 4.0t 0.3 (1t) 1 [0,1] (t).

  Fig. 7. Exact (--) and identified (-) frequency responses H(iω).

  Fig. 8. Exact (--) and identified (-) function f .

This convergence mode means that on any bounded set P , γn |Pγ |P → 0 and γ n |P -γ |P → 0 uniformly.

We distinguish carefully trajectories (u, x,...), which are functions of the time, and values taken by trajectories at time t (i.e. u(t), x(t),...).

Note that the source function is physically realistic but rather poor from the point of view of information, which strengthen the difficulty of the problem.


C Casenave
email: casenave@laas.fr

E Montseny
email: emontseny@laas.fr

H Camon
email: camon@laas.fr

Identification of Electrostatically Actuated MEMS Models from Real Measurement Data

Keywords: Dynamic models, parameters identification, nonlinear models, least-squares identification

This paper focuses on the identification of nonlinear dynamic models for physical systems such as Micro-Electro-Mechanical Systems (MEMS) from measurement data. The proposed approach consists in transforming, by means of suitable global operations, the specific input-output differential physical model of the system elaborated from physical analysis, in such a way that we get a new equivalent model formulation specifically adapted to the identification problem. Thanks to the equivalence of the dynamic model and the derived identification problem, the so-identified model remains of continuous-time type, with a clear physical meaning of any of its components, which is not the case when using, for example, black-boxes approaches. The method is implemented on real measurement data from a physical system.

INTRODUCTION

This paper focuses on the identification of nonlinear dynamic models for physical systems such as MEMS from measurement data. Due to their very small size, many parameters of those systems cannot be accurately known; in this case, identification process can be the only way to get reliable models, useful for example for control [START_REF] Zhu | Modeling and control of electrostatically actuated Mems in the presence of parasitics and parametric uncertainties[END_REF]; [START_REF] Liao | Closedloop adaptive control for torsional micromirrors[END_REF]; [START_REF] Sane | Robust control of electrostatic torsional micromirrors using adaptive sliding-mode control[END_REF]; [START_REF] Daqaq | Input-Shaping Control of Nonlinear MEMS Proceedings of the IMAC-XXII Conference and Exposition on Structural Dynamics[END_REF]; [START_REF] Bryzek | Control issues for MEMS Proceedings of the 42nd IEEE Conference on Decision and Control[END_REF]) or even dimensioning purposes. Several informations and techniques about identification of dynamic systems will be found in [START_REF] Ljung | System identification -Theory for the user[END_REF] and [START_REF] Garnier | Time-domain approaches to continuous-time model identification of dynamical systems from sampled data[END_REF].

In this paper, we identify parameters and functional components of a nonlinear dynamic model presenting a bifurcation point. The proposed approach consists in transforming, by means of suitable global operations, the specific input-output differential physical model of the system elaborated from physical analysis, in such a way that we get a new equivalent model formulation well adapted to the identification problem. By using spline interpolations, we establish a time-continuous identification problem which will be solved by classical least-square methods, with the only shortcoming of a small and easy to suppress identification bias. Moreover, thanks to the equivalence of the dynamic model and the derived identification problem, the so-identified model remains of continuous-time type, with a clear physical meaning of any of its components, which is not the case when using, for example, black-boxes approaches.

The paper is organized as follows. We present the physical system and its dynamic model, and then state the associated identification problem. After that we express the solution of this problem in terms of operatorial pseudoinversion and describe how to solve the problem numeri-cally from a discrete set of measurement data. Finally, we implement the method on real measurement data.

THE PHYSICAL SYSTEM UNDER CONSIDERATION

The system under consideration is an electrostatically actuated micromirror1 (cf. Fig. 1). It is composed in two assembled parts. The upper one is a thin plate, the mirror, linked to a thick external rigid frame by two thin and narrow arms, the springs. This part is tailored in the same microcrystalline silicon layer of a SOI (Silicon On Insulator) wafer. The lower one comprises a balance-knifeedge with two electrodes distributed on both sides of it.

The two parts are assembled in such a manner that the axis of the springs and balance-knife-edge are identical. So the electrodes are located underneath the mirror inducing its rotation (left or right) when a voltage V is applied.

The development of those micro-mirrors has been conducted with Tronics Microsystems (France), a manufacturer of custom MEMS components. Several configurations of the electrode can be envisaged, namely the case where electrodes are flat or inclined (cf. Fig. 2 and 3 for flat electrodes configuration). Due to the fabrication process, essentially the gluing process of the two parts, it is realistic to consider that the physical structure is significantly different from the ideal one, and so do the associated dynamic model.

In the sequel the physical limit angle the mirror can reach is denoted α, whereas θ(t) denotes the angle of the mirror at a given time (see Fig. 3). The dynamic model of the system is obtained by application of the fundamental principle of dynamics:

I θ + (μ 0 + v(θ)) θ + K θ = V 2 k(θ), (1)
with the constraint :

|θ| |α| , ( 2
)
where I is the inertia moment of the system, V 2 k(θ) is the electrostatic moment, K θ the spring moment, and (μ 0 + v(θ)) θ the viscous friction moment. The system (1) of input V and output θ, is completed by the initial conditions: θ(0) = 0 and θ(0) = 0.

(3) Remark 1. In the model, θ is a negative angle.

The expression of the terms v(θ) and k(θ) involved in model (1) depends on the configuration of the physical system (flat or inclined electrodes), its geometry, the materials used etc. The analysis of this dynamic system (Casenave (2008)) reveals the existence of a threshold voltage denoted V dyn pullin , below which the mirror stabilizes to an angle θ stab (V ) (see Fig. 4a), whereas beyond this voltage the mirror switches (see Fig. 4b).

We focus in this paper on the problem of identification of model (1) (namely identification of inertia moment I, electrostatic moment function k, etc.) from real measurement data. This problem is tricky for several reasons: the model is nonlinear, and presents a bifurcation point that separate two different dynamic behaviors of the system. Moreover, inertia identification is difficult because its contribution is often (but not always) dominated by viscosity terms.

FORMULATION AND RESOLUTION OF THE IDENTIFICATION PROBLEM

In this section, we state and study the identification problem under consideration, whose aim is to estimate the parameters of model (1,3) from measurements of trajectories θ.

Preliminaries for the identification problem

As said previously, the expression of electrostatic and viscous friction moments functions k and v depends on the physical system, and in general it may not be accurately known. We identify the two functions under the form:

k(θ) = l c l k l (θ) ( 4 ) v(θ) = q μ q v q (θ), (5)
where the parameters to be identified are the real coefficients c l and μ q . From the theoretical point of view, these sums can be finite or not and we implicitly refer to functional Hilbert spaces generated by k l and v q , to which belong k and v. In the sequel, we will use the notations:

c := (c 1 , c 2 , ...) T ∈ 2 , k(θ) := (k 1 (θ), k 2 (θ), ...), (6) μ := (μ 1 , μ 2 , ...) T ∈ 2 , v(θ) := (v 1 (θ), v 2 (θ), ...). (7)
So we can write in a condensed form:

k(θ) = k(θ) c , v(θ) = v(θ) μ. (8)
Then, we consider the possibility of transforming the model ( 1) in an equivalent manner by composition with any causal convolution operator

H := H(∂ t ) with impulse response h = L -1 H: (H • ∂ 2 t ) θ I + (H • ∂ t ) θ μ 0 + H( v(θ)∂ t θ) μ + (Hθ) K = V 2 H( k(θ)) c.
(9) To avoid the amplification of noise measurement by operators ∂ t and ∂ 2 t , the operator H(∂ t ) has to be chosen in such a way that high frequencies are attenuated, that is

|H(iω)| 2 ∼ H.F 1 ω 2 , without amplifying low ones H(iω) 2 ∼ L.F
1: basically, it behaves like a second order low-pass filter.

Finally, it is obvious that, in model ( 9), the parameters can only be identified up to a multiplicative constant. Then, we consider that μ 0 = 1 and denote:

λ := (I, μ, K, c) T ∈ E := R× 2 × R× 2 (10)
the vector of parameters to be identified.The exact (but unknown) value of λ is denoted λ 0 . We can then rewrite (9) under the synthetic form, used for identification purposes:

A θ λ = b θ , (11) where H 1 := H • ∂ t , H 2 := H • ∂ 2 t , b θ := -H 1 θ ∈ L 2 (0, T ; R J ) and A θ : E -→ L 2 (0, T ; R J ) is the matrix operator defined by: A θ = H 2 θ H( v(θ)∂ t θ) Hθ -V 2 H (k(θ)) . (12) Remark 1.
In practice, we will manipulate multiple trajectories θ j , j = 1 : J obtained with different input voltages V j (the associated initial conditions are all taken equal to 0). Then, model ( 9) can be extended to the general case of multiple trajectories by considering:

θ = (θ 1 , ..., θ J ) T , ∂ t θ = (∂ t θ 1 , ..., ∂ t θ J ) T v(θ) = ⎡ ⎢ ⎣ v 1 (θ 1 ) v 2 (θ 1 ) • • • . . . . . . v 1 (θ J ) v 2 (θ J ) • • • ⎤ ⎥ ⎦ , k(θ) = ⎡ ⎢ ⎣ k 1 (θ 1 ) k 2 (θ 1 ) • • • . . . . . . k 1 (θ J ) k 2 (θ J ) • • • ⎤ ⎥ ⎦ , and V = ⎡ ⎢ ⎣ V 1 0 . . . 0 V J ⎤ ⎥ ⎦ .
(13)

Formulation and resolution of the identification problem

The aim is to estimate λ from experimental data of the form: θ m = θ + η, with η some additive measurement noise such that, for simplicity, E(η) = 0. We then write the identification problem under the following form, taking advantage from the linearity of (11) with respect to parameter λ:

min λ∈E A θ m λ -b θ m 2 L 2 (0,T ;R J ) . ( 14
)
Its solution λ * is formally obtained by:

λ * = A † θm b θ m , ( 15
)
where A † θm designates the pseudo-inverse of operator A θ m , defined by [START_REF] Ben-Israel | Generalized inverses: theory and applications[END_REF]):

A † θm = (A * θm A θ m ) -1 A * θm n (16) with A *
θm the adjoint operator of A θ m . In particular, we have from (11): A θ λ 0 = b θ ; so, the estimator λ * is exact when the data are unoised:

A † θ b θ = λ 0 . ( 17
)
Because operator A θ m depends on the measurement noise, the solution λ * of ( 14) will be a biased estimator of λ 0 . However, an efficient iterative procedure can be used to reduce the bias and thus increase the identification accuracy (Casenave (2008)).

NUMERICAL FORMULATION

In this section, we show how to numerically solve the continuous time (infinite dimensional) identification problem introduced above.

Discrete data interpolation

In practice we only have access to discrete trajectories, that is several sets of data {θ j,k m } k=1:K , j = 1 : J, where θ j,k m is a measurement of θ(t k , V j ) with t k+1 = t k + Δt. To go back to a continuous-time formulation, we use cubic spline interpolations which allow the analytical computation of the involved quantities (such as A m and b m introduced later).

In the sequel we will denote θ j m : t → θ j m (t) the cubic spline interpolation of the set {(t k , θ j,k m )} k . We have:

θ j m = θ j + η,
(18) where θ j is the exact (unknown) trajectory solution of (9) (with V = V j ), and η is an additional noise resulting both from measurement errors and interpolation ones. In the same way, cubic spline interpolation will be made on the sets

{(t k , v q (θ j,k m )∂ t θ j,k m )} k and {(t k , k l (θ j,k m ))
} k to get a continuous-time identification problem of the form ( 14), whose numerical solution is obtained as described here after.

Numerical resolution of (14)

To numerically solve the minimisation problem ( 14) which is a priori of infinite dimension, we first have to replace the infinite sums ( 8) by truncations at finite orders L and Q:

k(θ j ) c L l=1 c l k l (θ j ) and v(θ j ) μ Q q=1 μ q v q (θ j ). (19)
In the sequel, k(θ) and v(θ) will now designate the respective finitedimensional matrices:

⎡ ⎢ ⎣ v 1 (θ 1 ) • • • v Q (θ 1 ) . . . . . . v 1 (θ J ) • • • v Q (θ J ) ⎤ ⎥ ⎦ and ⎡ ⎢ ⎣ k 1 (θ 1 ) • • • k L (θ 1 ) . . . . . . k 1 (θ J ) • • • k L (θ J ) ⎤ ⎥ ⎦ . ( 20
)
We then discretize the problem (14). Let {t n } n=1:N a time discretization. For any t n we have:

(A θ λ)(t n ) = b θ (t n ). ( 21
) Let A m : R × R Q ×R × R L -→ R J N the block-matrix: A m : = ⎡ ⎢ ⎣ A θ m (t 1 ) . . . A θ m (t N ) ⎤ ⎥ ⎦ (22)
where

A θ m (t n ) : R × R Q ×R × R L -→ R J is the matrix defined by: A θ m (t n ):= [(H 2 θ m ) (t n )(H(v(θ m )∂ t θ m )) (t n )(Hθ m ) (t n ) -V 2 (t n ) (H(k(θ m ))) (t n ) . (23) Let also define: b m := ⎡ ⎢ ⎣ -(H 1 θ m ) (t 1 ) . . . -(H 1 θ m ) (t N ) ⎤ ⎥ ⎦ ∈ R J N . ( 24
)
Now, we can write the numerical formulation of the identification problem ( 14) (with possible reconditioning small parameter ε 0):

min λ ||A m λ -b m || 2 + ε λ 2 , ( 25
)
whose solution is classically given by:

λ * = (A T m A m + εI) -1 A T m b m .
(26) Remark 2. The dimension of the square matrix A T m A m +εI is rather small (equals to the number of coefficients to be identified).

APPLICATION TO REAL MEASUREMENT DATA

In Casenave (2008) the method was first validated on simulated data. The results obtained were good: all parameters of the model (inertia, stiffness, electrostatic moment function etc.) were accurately identified (error less than 1%), and identified curves perfectly fit the data. We implement in this section the identification method on real data, measured on a MEMS like described in section 2.

Physical system and measurement data

The MEMS under consideration has flat electrodes as described in Fig. 3, with α = -0.0215 rad, a 1 = 150 μm, a 2 = 280 μm and a 3 = 300 μm. The width of the mirror is equal to 600 μm and its thickness is equal to 10 μm [START_REF] Camon | Solving functional reliability issue for an optical electrostatic switch[END_REF]).

Data available for identification are composed of J = 36 sets of data {θ j,k m } k=1:K with K = 501, obtained by sampling at frequency 200 kHz during 10 ms of 36 measured trajectories of θ. The associated inputs voltages are described here-after. During the first 5 ms, the voltage between electrodes is fixed at a constant value V j ; thus, the system either switches or stabilizes at θ j ∞ < α. During the next 5 ms, the two electrodes are at the same potential: the system then asymptotically returns to its rest position (say θ = 0, θ = 0). The 36 voltage values V j are distributed from 6.085 V to 97.36 V, with

V j < V j+1 .
Two examples of data trajectories are given in Fig. 4a and4b. Note, that only the parts of data which are compatible with the model are used for identification. For example, in Fig. 4b, the saturation zone [0.0045 s, 0.005 s] where θ = α is not relevant (saturation).

Identification model

The identification model considered in the sequel is of the form (9) with: k(θ) identified under the form ( 19) with L = 6 and: k l (θ) = θ l-1 , l = 1 : 6. ( 27) Indeed, we suppose the electrostatic moment can be accurately described (in the domain [0, α] in which θ evolves) and identified by means of a fifth order polynomial.

Remark 3. In fact, a simplified physical analysis leads to the expression: M e (θ, V ) ∝ V 2 g(θ) with [START_REF] Camon | Fabrication, simulation and experiment of a rotating electrostatic silicon mirror with large angular deflection, 13th[END_REF]):

g(θ) = 1 θ 2 1 1 -a2 a3 θ α - 1 1 -a1 a3 θ α + ln 1 -a2 a3 θ α 1 -a1 a3 θ α
.

(28) If this expression would be reliable, then we could identify the electrostatic moment under this form by simply considering ( 19) with L = 1 and k 1 = g. However, due to several approximations in the physical analysis, the function (28) in fact reveals itself ill adapted for large values of θ.

(a) V = 18.25 V < V dyn pullin (b) V = 28.72 V > V dyn pullin Fig. 4. Data θ j,k
m for different values of V j v(θ) identified under the form ( 19) with Q = 1 and:

v 1 (θ) = θ 3 . ( 29
)
Indeed, a linear viscosity term μ 0 θ is not sufficient to correctly describe the viscous moment in its whole, namely when both |θ| and | θ| become large.

for the prefiltering operator H, a simple second order low-pass filter with transfer function:

H(p) = σ 2 (p + σ) 2 , ( 30
)
with σ = 9.0 × 10 3 > 0 (the cutoff frequency) chosen in such a way that A θ m λb θ m 2 be "small". We so have:

(H(∂ t ) u)(t) = σ 2 t 0 (t -s) e -σ(t-s) u(s) ds, ( 31
)

Identification results

The algorithm described in section 4 has been implemented on the sets of trajectories {θ j,k m }, j = 8 : 20, associated with V j ∈ [18.25 V, 31.64 V]: only significant data of the first 5 ms are considered. The system being initially at rest, the initial conditions are equal to 0. We obtain the following identified parameters: As said previously, the parameters are identified up to the multiplicative constant μ 0 . To deduce the physical parameters from the identified quantities, we must know (or compute) the value of one parameter, from which we can deduce all the others. For example, we compute here the physical value Ĩ of inertia moment from the dimensions of the MEMS (see Fig. 3) and the density of silicon (2.33 × 10 3 kg/ m 3 ), and then deduce other physical quantities:

I * = 1.944 × 10 -5 , K * = 3.874 × 10 3 , (32) μ * 1 = -5.061 × 10 5 , c * = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 0.
Ĩ = 2.693 × 10 -16 N m s 2 / rad, (computed) (34) μ0 = Ĩ I * μ 0 = 1.385 × 10 -11 N m s/ rad, ( 35
) K = Ĩ I * K * = 5.366 × 10 -8 N m/ rad, (36)
and the same for μ * 1 and c * . In order to globally estimate the identification quality, we consider the quantity:

E = j,k θ j * (kΔt) -θ j,k m j,k θ j,k m j KΔt 0 θ j * (t) -θ j m (t) j KΔt 0 θ j m (t)
, where θ j * is the trajectory of θ obtained by simulation of the identified model for an input voltage equal to V j . This quantity represents the cumulated relative error on all the trajectories used for the identification. In the same way, we will denote E g the cumulated relative error extended to all the available measurement trajectories (thus including prediction errors). We obtain: E = 0.0071 and E g = 0.0095, (37) which highlight the good precision of the identification process. In Fig. 5 we can see some examples of trajectories obtained with the identified model for V 10 = 20.69 V < V dyn pullin < V 18 = 29.82 V. Those trajectories fit well the measured data. As a validation test, we show here after the results obtained in predictive situation. We give in Fig. 6 some data which has not been used for identification and their associated trajectories predicted by the identified model. Here again, all the trajectories fit measured data with a high accuracy. Note that the voltage value V max = V 36 = 97.36 V we used for prediction results is far from the ones used for identification (whose larger value is 31.64 V): that means there is a ratio almost equal to 10 between the input V 2 of predicted curves and this of curves used for identification. In such a situation, significant identification errors on the system parameters would generate large prediction errors; in our case, the predictive results remain good and confirm the identification accuracy.

Finally, the identified electrostatic function k is given in Fig. 7 and compared with physical approximation of the form (28), which clearly reveals itself inadequate for deviations of θ close to α. As an illustration, we can see in Fig. 8 an example of simulated trajectory when the

CONCLUSION

In this paper, we proposed an identification method for a nonlinear dynamic model of electrostatically actuated MEMS, using operatorial prefiltering to take advantage from the linearity with respect to the parameters to be identified. This method is applied to real measurement data, and the results are quite good: the identified model is relevant and accurate, even in prediction situations.

The accuracy of the so-identified models suggests several applications in the future. Namely, a high precision tension reference can be built by controlling the micromirror around its semi-stable bifurcation point.

Fig. 1 .

 1 Fig. 1. View of the physical system

Fig. 5 .

 5 Fig. 5. Measured trajectories and the associated θ simulated from the identified model for two values of V function k is under the form (28): the model does not fit with the data, contrary to the the case when k is fully identified (see Fig. 5(b)).

Fig. 7 .

 7 Fig. 7. The identified k(θ) and its theoretical estimation (28).

Fig. 6 .Fig. 8 .

 68 Fig. 6. Measured trajectories and the associated θ predicted by the identified model for different values of V

The picture is published with courtesy of Tronics Microsystems (France).


C Casenave
email: casenave@laas.fr

G Montseny
email: montseny@laas.fr

Diffusive Identification of Volterra Models by Cancellation of the Nonlinear Term

Keywords: Volterra Model, Nonlinear Model, Identification, Pseudo-inverse, Diffusive Representation, Operatorial Cancellation

We present a new identification method for nonlinear Volterra models of the form HX = F (u, X) + v with H a causal convolution operator. It is based first on a suitable parameterization of H deduced from the so-called diffusive representation, devoted to state representations of integral operators, and second on a suitable operatorial transformation of the problem with the property that the nonlinear term F (u, X) is cancelled, allowing to identify the operator H(∂ t ) alone. Then, the nonlinear term can be identified from standard regression methods. For illustration, we implement this method on a concrete numerical example.

INTRODUCTION

We present an identification method for nonlinear Volterra models of the form:

H(∂ t )X = F (u, X) + v, (1)
where u(t), v(t), X(t) ∈ R, t > 0, F is a nonlinear function defined on R 2 and H(∂ t ) is a causal convolution operator defined on a suitable space of continuous functions with support in R + t . We suppose that problem (1) is well-posed in the sense of existence, uniqueness and continuous dependency on the inputs u, v of the solution X. The operator H(∂ t ) is also supposed invertible, so (1) can be rewritten under the form:

X(t) = t 0 k(t -s) [F (u(s), X(s)) + v(s)] ds ∀t > 0,
where k is the impulse response of operator H(∂ t ) -1 . Such models are frequently encountered in various domains:

• thermic phenomena (we will consider in section 4 the combustion model elaborated in [START_REF] Joulin | Point-source initiation of lean spherical flames of light reactants: an asymptotic theory[END_REF]), • electrical engineering [START_REF] Rumeau | Behavior modeling of a CaCu 3 Ti 4 O 12 ceramic for capacitor applications[END_REF]),

• linear SISO differential systems with nonlinear feedback:

Ż = AZ + BF (u, X) X = CZ, with H(p) -1 = C(pI -A) -1 B, • SISO partial differential systems on R + * t ×Ω z ⊂ R n+1 of the form: ∂ t ϕ = A(z, ∇)ϕ + F (u, T ϕ dz), ϕ 0 = 0, X = T ϕ dz,
which can be rewritten under the synthetic equation t 0 k(t-s) F (u(s), X(s)) ds = X(t) where the function k is the impulse response solution of k = T φ dz, φ = A(z, ∇)φ + δ, φ 0 = 0. • etc.

The identification problem under consideration in the sequel is to build estimations of H(∂ t ) and F from (noised) data X = X + η where η designates some additive measurement noise and X is generated by an experimental process driven by known inputs u, v. The main difficulty when identifying such models results from the coupling between the dynamic operator H(∂ t ) and the (static) function F via equation (1). This difficulty is strengthened when operator H(∂ t ) is non rational (which is in general the case), when F is singular, when the system is unstable and/or when there exists some dynamic bifurcations (see section 4).

The proposed identification method is based first on a suitable parameterization of H(∂ t ) deduced from the socalled diffusive representation [START_REF] Montseny | Représentation Diffusive[END_REF]), devoted to state realizations of general integral causal operators. Following this approach, the complex dynamic nature of H(∂ t ) is in some sense summarized by a few numerical parameters on which the identification problem will focus. Second, we introduce a suitable operatorial transformation of the problem with the property that the nonlinear term F (u, X) is cancelled, allowing to identify the operator H(∂ t ) alone from the available data. The nonlinear term can then be identified from the data completed by the knowledge of the previously identified H(∂ t ). This paper must be viewed as an extension of a previous one (Casenave (2008-2)) to which we broadly refer. We only present here the basic principles of the approach; more details, namely numerical questions, will be given in a future publication.

The paper is organized as follows. In section 2, we briefly present a simplified version of the diffusive representation. In section 3, we describe the proposed identification method in a general framework. In section 4, we finally implement this method on a concrete example, in order to highlight the relevance of the method.

DIFFUSIVE FORMULATION OF CAUSAL CONVOLUTION OPERATORS

We only give a few basic notions; more details will be found in Casenave (2008-2). A complete statement of diffusive representation will be found in [START_REF] Montseny | Représentation Diffusive[END_REF]. Various applications and questions relating to this approach will be found for example in [START_REF] Audounet | An asymptotic fractional differential model of spherical flame[END_REF], Carmona (1998), [START_REF] Casenave | Optimal identification of delay-diffusive operators and application to the acoustic impedance of absorbent materials[END_REF], Casenave (2008-1), [START_REF] Degerli | Bandlimited 1/f α noise source[END_REF], [START_REF] Garcia | Identification of the dynamics of a lead acid battery by a diffusive model[END_REF], [START_REF] Lenczner | Diffusive realization of operators solutions of certain operational partial differential equations[END_REF], [START_REF] Levadoux | Diffusive formulation of the impedance operator on circular boundary for 2D wave equation[END_REF], [START_REF] Montseny | Diffusive representation for operators involving delays[END_REF][START_REF] Ph | Identification of a 2D turbulent wind spectrum[END_REF], [START_REF] Rumeau | Behavior modeling of a CaCu 3 Ti 4 O 12 ceramic for capacitor applications[END_REF].

We consider a causal convolution operator defined, on any continuous function u : R + → R, by

u → t 0 k(t -s) u(s) ds. ( 2
)
We denote K the Laplace transform of the (locally integrable) function k, and K(∂ t ) the convolution operator defined by (2).

Let u t (s) = 1 ]-∞,t] (s) u(s) the restriction of u to its past and u t (s) = u t (ts) the so-called "history" of u. From causality of K(∂ t ), we deduce: [K(∂ t )(uu t )](t) = 0 for all t; then, we have for any continuous function u:

[K(∂ t ) u](t)= L -1 (K Lu) (t)= L -1 K Lu t (t). (3)
We define:

Ψ u (t, p) := e p t Lu t (p) = (Lu t ) (-p); (4)
by computing ∂ t Lu t , Laplace inversion and use of (3), we have: Lemma 1. 1. The function Ψ u is solution of the differential equation:

∂ t Ψ(t, p) = p Ψ(t, p) + u, t > 0, Ψ(0, p) = 0. 2. For any b 0, [K(∂ t ) u] (t) = 1 2iπ b+i∞ b-i∞ K(p) Ψ u (t, p) dp.
(5)

We denote Ω the holomorphic domain of K (after analytic continuation). Let γ be a closed 1 simple arc in C -; we denote Ω + γ the exterior domain defined by γ, and Ω - γ the complementary of Ω + γ . By use of standard techniques (Cauchy theorem, Jordan lemma), it can be shown:

Lemma 2. For any γ ⊂ Ω such that K is holomorphic in Ω + γ , if K(p) → 0 when p → ∞ in Ω + γ , then: [K(∂ t ) u] (t) = 1 2iπ γ K(p) Ψ u (t, p) dp. ( 6
)
In the sequel, we use the convenient notation μ, ψ = μ ψ dξ.

Under hypothesis of lemma 2 (except γ ⊂ Ω), we have [START_REF] Montseny | Représentation Diffusive[END_REF] Formulation (8,7) can easily be extended to operators of the form

K(∂ t ) • ∂ t where K(∂ t ) admits a γ-symbol. We have (formally): [K(∂ t ) • ∂ t u](t) = μ, ∂ t ψ u (t, .) = μ, γ ψ u (t, .) + u(t) , (9) with μ the γ-symbol of K(∂ t ).

IDENTIFICATION VIA OPERATORIAL CANCELLATION OF THE NONLINEAR TERM

The basic principle

Given E, F, G separable Hilbert spaces, f : E → G a nonlinear function and A : F → G linear, consider the abstract model:

f (x) + Ay + z = 0, (10)
where x, y, z are trajectories defined on t ∈ [0, T ] with values respectively in E, F, G. The problem under consideration is to estimate both A and f from (possibly) noised data trajectories ỹ, x, z, that is to solve (in some sense):

min A,f f (x) + Aỹ + z 2 , ( 11
)
with . an hilbertian norm, for example 10) and ε i are some (sufficiently small) noises. Such a problem often presents some serious difficulties, in particular when f is highly nonlinear (for example non differentiable, discontinuous or even weakly singular) and no a priori information about it is available. In such cases indeed, linear parametrizations2 of f of the form f = i a i f i (with (f i ) a topological basis of E chosen a priori) cannot be accurate enough except if the number of a i is quite large; but on the other hand, when the number of parameters to be estimated is excessive, estimation is generally bad due to excessive sensitivity to the measurement noises. Now remark that thanks to linearity, operator A could be easily estimated if the term f (x) was known, the problem (11) then reducing to a standard least-square one:

T 0 . 2 G dt, x = x + ε 1 , ỹ = y + ε 2 , z = z + ε 3 where (x, y, z) is solution of (
min A Y • A -b 2 ,
with solution formally given by the pseudo-inverse of operator Y (Ben-Israel ( 2003)):

A * = Y † • b.
(12) Following this remark, the method presented here-after allows to estimate A alone from a suitable transformation of model ( 10) with the special property that the term f (x) is cancelled 3 in a particular subset of equations specifically defined from the data x. Under some hypothesis frequently satisfied in practice, this "submodel" then reveals itself sufficient for a good estimation of A under the form (12). Once operator A is correctly estimated, it will then be easy to deduce estimations of f from standard methods.

The f -cancellation operator D x,ε

We briefly introduce, in an abstract framework, the "difference" operator D x,ε on which the proposed method is based.

Given X a Banach space and Ω a metric space, let denote by C 0 (Ω; X) the Banach space of continuous 4 functions defined on Ω with values in X, with norm

x = sup ξ∈Ω x(ξ) X . Given X 0 ⊂ X, let C 0 (Ω; X 0 ) := {x ∈ C 0 (Ω; X); ∀ξ ∈ Ω, x(ξ) ∈ X 0 }. Given metric spaces 5 Ω i , i = 1 : n, the norm on the product space n i=1 C 0 (Ω i ; X)
is defined by x = sup i sup ξ∈Ω i x i (ξ) X . Finally, we denote Ω i,j := Ω i × Ω j and the following operator D is defined:

D : n i=1 C 0 (Ω i ; X) → n i,j=1 C 0 (Ω i,j ; X) x = (x i ) i=1:n → Dx = (D i,j x) i,j=1:n , ∀i, j ∈ {1 : n}, ∀(t, τ ) ∈ Ω i,j , (D i,j x)(t, τ ) = x i (t)-x j (τ ).
(13) The operator D is clearly linear. Thanks to the continuity of the + operation in Banach spaces, it is continuous; note that this last property is essential, namely for numerical treatments.

Let Y be another Banach space, X 0 ⊂ X and f : X 0 → Y continuous. For x ∈ n i=1 C 0 (Ω i ; X 0 ), we denote by convention:

f • x := f (x i ) i=1:n ∈ n i=1 C 0 (Ω i ; Y). (14)
Note that if f is linear, we have:

D(f • x) = f • (Dx), this
property being false in general.

Now consider the sets 6 (see Figure 1):

Ω i,j x,ε := (t, τ ) ∈ Ω i,j ; (D i,j x)(t, τ ) ε ; (15) for any (i, j) s.t. Ω i,j x,ε = ∅, we denote D i,j
x,ε and D x,ε the operators defined by ∀y, D i,j x,ε y :

= D i,j y| Ω i,j x,ε and D x,ε := (D i,j x,ε ) i,j ; (16)
we have: Proposition 1. For any x ∈ n i=1 C 0 (Ω i ; X 0 ) and any continuous f : X 0 → Y, the following properties hold:

(1) D x,0 x = 0, (2) D x,0 (f • x) = 0, (3) D x,ε x ε,
3 Up to a constant term: f 0 := f (x(t 0 )) ∈ G which will be identified simultaneously. 4 This hypothesis can be weakened by means of suitable technical adaptations. 5 In practice, Ω i will be time intervals

[t i 0 , t i f ]. 6 The sets Ω i,i
x,ε are never empty because (D i,i x)(t, t) = 0. Furthermore, thanks to the continuity of Dx, we have ε>0 Ω i,j

x,ε = Ω i,j x,0 .

(4)

∀(t ε , τ ε ) ∈ Ω i,j x,ε s.t. (t ε , τ ε ) → ε→0 (t, τ ) ∈ Ω i,j x,0 , lim ε→0 D i,j x,ε (f • x) (t ε , τ ε ) = 0, (5) If f is k-Lipschitz, then D x,ε (f • x) kε.
Proof.

(1) and ( 2) ∀i, j ∈ {1 : n}, ∀(t, τ ) ∈ Ω i,j x,0 , we have (D i,j x)(t, τ ) = x i (t)x j (τ ) = 0, that is x i (t) = x j (τ ), and so f (x i (t)) = f (x j (τ )); then:

(D i,j (f •x))(t, τ ) = f (x i (t))- f (x j (τ )) = 0.
(3) is a direct consequence of the definition of Ω i,j

x,ε and of the norm on i,j C 0 (Ω i,j

x,ε ; X).

(4) First remark that (t, τ ) ∈ Ω i,j x,0 ⇒ x i (t) = x j (τ ); thanks to the continuity of f and x, we then have: lim

ε→0 D i,j x,ε (f • x) (t ε , τ ε ) = lim ε→0 f (x i (t ε )) -f (x j (τ ε )) = f (x i (t)) -f (x j (τ )) = 0. (5) D x,ε (f • x) = sup i,j sup Ω i,j x,ε f (x i (t)) -f (x j (τ )) sup i,j sup Ω i,j x,ε k x i (t) -x j (τ ) kε.
Remark 1. The estimate ( 5) is uniform, while (4) relates to the weaker topology of simple convergence. For this reason, it can be expected that when f is not Lipschitz-continuous on X 0 , numerical treatments can be somewhat more delicate. However, in many concrete situations, function f is Lipschitz-continuous on subsets X η ⊂ X 0 , in such a way that the property ||D x,ε ( fη • x)|| k η ε holds for any k η -Lipschitz-continuous function fη (defined on X 0 ) such that fη | Xη = f | Xη . We will see on a concrete example (with X = Y = R, X 0 = R + * , f = ln, X η = [η, +∞[) that such a property (in addition to (4)) can be sufficient to get good identification results.

For convenience, we will denote in the sequel:

Ω x,ε := i,j {(i, j)} × Ω i,j x,ε ; (17)
this set will be called the ε-kernel of Dx. Roughly speaking, it can be viewed as the set of all the (i, j, t, τ ) such that x i (t)x j (τ ) ε. So, for any y ∈ n i=1 C 0 (Ω i ; Y), the function D x,ε y is simply defined by: D

x,ε y : Ω x,ε → Y (i, j, t, τ ) → y i (t) -y j (τ ). ( 18
)
Note that for numerical purposes, we will mainly have to consider discrete subsets of Ω x,ε , with elements (i, j, t k , τ l ) ∼ (i, j, k, l) ∈ N 4 . Numerical questions (in particular the choice of ε and the construction of a discrete set Ω x,ε from a discrete-time trajectory x) will be studied in a further publication.

The D x,ε -transformed identification problem

Let us consider the model (10): f • x + Ay + z = 0, now with (x, y, z) a possibly vector trajectory with values in

E n × F n × G n , with x i , y i , z i defined on Ω i = [t i 0 , t i f ]
. By applying operator D x,0 to both members of this equation, we get D x,0 (f • x) + A D x,0 y + D x,0 z = 0 and, thanks to property (2) of proposition 1, A D x,0 y + D x,0 z = 0; note that this last equation is linear with respect to the (unknown) operator A. However, due to the fact that Fig. 1. Examples of (t, τ ) in Ω i,i

x,ε (top) and in Ω i,j

x,ε (bottom) operator D x,0 is based on differences, equivalence is lost.

To restore equivalence, it is necessary 7 to consider an "initial" condition (at least one!) on (10):

f (x i (t i * )) + Ay i (t i * ) + z i (t i * ) = 0 for some t i * ∈ [t i 0 , t i f ]; so, with n * n and by denoting f 0 := f (x i (t i * )) ∈ G n * , y 0 := (y i (t i * )) ∈ F n * , z 0 := (z i (t i * )) ∈ G n * , the model: A D x,0 y = -D x,0 z f 0 + A y 0 = -z 0 ( 19
)
can be equivalent to (10) (that is (x, y, z) is solution of (10) iff it is solution of ( 19)) if the set Ω x,0 is large enough 8 . The interest of this new formulation is that both the unknowns A and f 0 are linearly involved; formally, these unknowns can then be determined by means of the pseudoinverse of operator Y 0 : (A, f 0 ) → (A D x,0 y, f 0 + A y 0 ).

In general however, on the one hand the set Ω x,0 may be too poor to get the strict equivalence of ( 19) and ( 10) and, on the other hand, the available data x, ỹ, z are noised. But in fact, it is only expected that some "good" identification of (A, f 0 ) can be computed from the available data; so, thanks to properties ( 4), ( 5) of proposition 1 which allow to neglect (in some sense) the term D x,ε (f • x), problem ( 19) is replaced by the weakened one, based on operator

Y ε : (A, f 0 ) → (A D x,ε ỹ, f 0 + A ỹ0 ): min (A,f0) Y ε • (A, f 0 ) + (D x,ε z, z0 ) 2 , ( 20
) with formal solution (A * , f * 0 ) = -Y † ε • (D x,ε z, z0 ). If data x, ỹ, z allow to get (A * , f * 0 ) (A, f 0 ) 9
, then we can expect from continuity properties: A * (ỹ i (t)ỹj (τ )) + zi (t)zj (τ ) 0, ∀(i, j, t, τ ) ∈ Ω x,ε , and so, although f is yet entirely unknown 10 at this stage:

f (x i ) + A * ỹi + zi 0 in the sense f (x i (t)) + A * ỹi (t) + zi (t) 0, ∀t ∈ [t i 0 , t i f ]
. Now, thanks to the property 7 But maybe not sufficient if data are not rich enough. 8 This should be clarified in a further work. 9 In particular if n is sufficiently big and noises ε i are small enough. This will be studied in a further work.

10 Except its value at point

x 0 = (x i (t i * )) ∈ E n * .
f (x i (t)) -Aỹ i (t)zi (t), identification of the nonlinear function f can be achieved by means of classical regression methods implemented on the following "pseudo graph" of f deduced from the available data:

G f = i={1:n},t∈[t i 0 ,t i f ] { xi (t), -A * ỹi (t) -zi (t) } ⊂ E × G.
(21)

Application to Volterra models under diffusive formulation

We start from the diffusive formulation of model (1). Given the following diffusive realization of operator H(∂ t ): H(∂ t ) X = μ, γ ψ X + X , with ψ X the diffusive representation of X, and μ the γsymbol associated with operator

H(∂ t ) • ∂ -1 t , the model (1) is equivalently rewritten: μ, γ ψ X + X -F • (u, X) -v = 0. ( 22
)
Thanks to the static nature of the linear operator ϕ → μ, ϕ , that is ∀t, ( μ, ϕ )(t) = μ, ϕ(t) , ( 22) can be rewritten (in opposite to (1) which is a t-nonlocal model):

∀i, t, μ, γ ψ X i (t) + X i (t) -F (u i (t), X i (t)) -v i (t) = 0. ( 23
)
Then, to apply the identification method presented above to model ( 22), we must take:

Ω i := [0, t i f ], x := (u, X), y := γ ψ X + X, z := -v, f := -F , A : y → μ, y . Indeed,
it can be verified that model ( 22) is then written under the form (10). In particular, the unknowns to be identified are (μ, f 0 ) and the operator Y ε defined in the previous section is expressed by

Y ε : (μ, f 0 ) → ( μ, D x,ε ỹ , f 0 + μ, ỹ0 ).
Remark 2. Thanks to proposition 1, some prefiltering Q(∂ t ) can be envisaged on the data11 X and v. Indeed, by applying operator Q(∂ t ) to (1), we get the equivalent model

H(∂ t ) • Q(∂ t ) X = Q(∂ t )F (u, X) + Q(∂ t )v; if, in a suitable sense, Q(∂ t )F (u, X) F (u, X), then D x,0 Q(∂ t )F (u, X) D x,0 F (u, X) = 0.

A CONCRETE EXAMPLE

In this section, we illustrate the identification method introduced above by implementing it on data elaborated from numerical simulations of a complex dynamic phenomenon studied in [START_REF] Joulin | Point-source initiation of lean spherical flames of light reactants: an asymptotic theory[END_REF], [START_REF] Audounet | An asymptotic fractional differential model of spherical flame[END_REF].

The model under consideration

In [START_REF] Joulin | Point-source initiation of lean spherical flames of light reactants: an asymptotic theory[END_REF], Joulin elaborated a Volterra model to describe, in suitable thermodynamic conditions, the evolution of a spherical flame initiated by a source at point 0 in a mixture of reactive species. Under some reasonable physical hypothesis, such a phenomenon can be described by a system of two partial differential equations relating to the temperature and the mass density of the mixture. By considering the reactive zone as a thin sheet located on a sphere with radius x(t), Joulin has established that when the flame is developing in free space, x is solution of the following nonlinear singular Abel-Volterra equation 12 (u(t) designates the source strength at time t):

x(t) t 0 ẋ(s) ds π(t -s) = x(t) ln x(t) + u(t) ∀t > 0, ( 24
)
with the additional conditions: x(0 + ) = 0, u 0, x 0 (whose physical interpretation is obvious). By denoting H(∂ t ) the convolution operator 13 x → t 0 ẋ(s) ds √ π(t-s) and F (u, x) + v := ln x + u x , (24) can be formally rewritten under the form (1).

It has been shown in [START_REF] Audounet | An asymptotic fractional differential model of spherical flame[END_REF] that the evolution problem ( 24) is well-posed, that is the solution x exists, is unique and depends continuously on u. In real conditions however, various perturbations are involved in the evolution of x (due for example to the loss of spatial symmetries), and both the convolution operator H(∂ t ) and the function F will be more or less far from the ideal ones. So, an identification process can be justified if accuracy of the model is required.

We consider the problem of identification of ( 10) from (noised) data (u, x) obtained from highly accurate numerical simulations of (24) (not described here, see Casenave (2008-1) for numerical method for such simulations).

On the dynamic behavior of the flame radius

It has been shown in [START_REF] Audounet | An asymptotic fractional differential model of spherical flame[END_REF] that there exists a threshold relating to the power of the source u, beyond which the flame is developing whereas a quenching occurs below. In that sense, this evolution phenomenon is essentially unstable with two qualitatively different behaviors.

Because of the hereditary nature of the problem, it is difficult to know the value of this threshold which must be evaluated on the basis of numerical simulations. In Casenave (2008-2), examples are given of flame either quenching or developing, when the source function is defined, as in [START_REF] Audounet | An asymptotic fractional differential model of spherical flame[END_REF], by:

u(t) = E t 0.3 (1 -t) 1 [0,1] (t). ( 25
)
The shortcomings of such sensitive dynamic behaviors are reinforced by the fact that the input u is dynamically poor and f is singular, which generates ill-conditioned identification problems. As a consequence, this model can be viewed as a significant test for the proposed identification method.

Formulation of the identification problem

For identification purposes, the flame model under consideration can be written under the form (1); by denoting:

f (x) : = F (u, x), z := u x , y := γ ψ x + x, Ay = μ, y , X = Y = R, X 0 = R + * ,
we can then apply the method described in section 3. A detailed description of the practical implementation of this method will be given in a future paper.

Numerical results and comments

Recall that the problem mainly consists in identifying the diffusive symbol μ of operator

H(∂ t ) • ∂ -1 t 13 Note that in this ideal case, H(∂t) = ∂ 1 2
t . Under weak hypothesis on γ, this operator admits a γ-symbol μ (see (Casenave (2008-2))).

from numerical experimental data. These data are composed of 9 trajectories (u i , xi ) i=1:9 (see figure 2) associated with 9 sources u i of the form ( 25) with E = 1.5, 1.6, 1.7, 1.73, 1.74, 1.75, 2.0, 3.0 and 5.0. The time step Δt has been taken equal to 10 -5 and the maximal final time is t f = 6. The measured trajectories are xi = x i + η v i where v i is a unity gaussian white noise and η = 10 -2 .

The function γ has been chosen as γ(ξ) = |ξ| e iα sign(ξ) with α = 160 180 π. The ξ-discretization has been defined from 130 values ξ l geometrically spaced to cover 5 frequency decades from 10 -2 to 10 5 (see Casenave (2008-2)).

In figure 3, we can see the so-identified diffusive symbol μ * associated with H * (∂ t )• ∂ -1 t . The associated frequency response H * (iω) is given in figure 4. In the covered frequency band, this frequency response is correctly identified. The results are close to the ones obtained in Casenave (2008-2) in the case where function f was known.

As an example, two of the resulting pseudo graphs of f , namely (see section 3.3):

G i f = t k { xi (t k ), -A * ỹi (t k )zi (t k ) } are visible in figure 5. The function f = ln is correctly reconstructed. Note that although the function f is not Lipschitz-continuous around x = 0, property (5) of proposition 1 remains satisfied around this point.

Finally, identification of function f has been achieved by means of a standard least-square method implemented on the whole available pseudo graph, given by ( 21); the result is shown in figure 6. In spite of the poverty of the data, in particular at high frequencies, the identification remains good.

Fig. 2 .

 2 Fig. 2. The trajectories x used for identification.

Fig. 3 .

 3 Fig. 3. The identified diffusive symbol μ * (ξ l ).

Fig. 4 .

 4 Fig. 4. Frequency response of the identified H(iω).

Fig. 5 .

 5 Fig. 5. Two of the pseudo grapfs G i f . C. Casenave, E. Montseny, Time-local dissipative formulation and stable numerical schemes for a class of integrodifferential wave equations, SIAM Journal on Ap-plied Mathematics, Vol.68, N6, pp.1763-1783[START_REF] Casenave | Identification of Nonlinear Volterra Models by means of Diffusive Representation[END_REF]

Fig. 6 .

 6 Fig. 6. Identified f * .

The advantage of which is to transform (10) into a linear problem.

The new prefiltered data are v := Q(∂t)v, X = Q(∂t)X.