Da-Yan Liu
email: dayan.liu@insa-cvl.fr

Gang Zheng b Driss Boutat
email: driss.boutat@insa-cvl.fr

Hao-Ran Liu
email: liu.haoran@ysu.edu.cn

Non-asymptotic fractional order differentiator for a class of fractional order linear systems ⋆

Keywords: Fractional order linear systems, Riemann-Liouville fractional derivatives, Non-asymptotic fractional order differentiator, Algebraic method, Noise error analysis

This paper aims at designing a non-asymptotic fractional order differentiator for a class of fractional order linear systems to estimate the Riemann-Liouville fractional derivatives of the output in discrete noisy environment. The adopted method is a recent algebraic method originally introduced by Fliess and Sira-Ramirez. Firstly, the fractional derivative of the output of an arbitrary order is exactly given by a new algebraic formula in continuous noise free case without knowing the initial conditions of the considered system. Secondly, a digital fractional order differentiator is introduced in discrete noisy cases, which can provide robust estimations in finite-time. Then, some error analysis is given, where an error bound useful for the selection of the design parameter is provided. Finally, numerical examples illustrate the accuracy and the robustness of the proposed fractional order differentiator.

Introduction

Fractional calculus has a long history and has been becoming very useful in many scientific and engineering fields, including control, flow propagation, signal processing, electrical networks, etc. [START_REF] Podlubny | Fractional differential equations: an introduction to fractional derivatives, fractional differential equations[END_REF][START_REF] Miller | An Introduction to the Fractional Calculus and Fractional Differential Equations[END_REF][START_REF] Kilbas | Theory and Applications of Fractional Differential Equations[END_REF][START_REF] Diethelm | The analysis of fractional differential equations: An application-oriented exposition using differential operators of Caputo type[END_REF][START_REF] Concepción | Fractional-order Systems and Controls: Fundamentals and Applications[END_REF][START_REF] Bandyopadhyay | Stabilization and Control of Fractional Order Systems: A Sliding Mode Approach[END_REF]. For instance, fractional order systems and controllers have been applied to improve performance and robustness properties in control design [START_REF] Yin | Fractionalorder sliding mode based extremum seeking control of a class of nonlinear systems[END_REF][START_REF] Victor | Parameter and differentiation order estimation in fractional models[END_REF][START_REF] Sabatier | On observability and pseudo state estimation of fractional order systems[END_REF], where the fractional derivatives of the output usually need to be estimated from its discrete noisy observation. Consequently, an interesting research topic concerns with designing digital fractional order differentiators, which should be robust against noises.

Various robust fractional order differentiators have been proposed both in the frequency domain and in the time domain. They can be divided into two classes: model-free fractional order differentiators [START_REF] Machado | Calculation of fractional derivatives of noisy data with genetic algorithms[END_REF][START_REF] Machado | Exploiting sensor redundancy for the calculation of fractional derivatives in the presence of noise[END_REF][START_REF] Chen | Digital Fractional Order Savitzky-Golay Differentiator[END_REF][START_REF] Liu | Fractional order differentiation by integration and error analysis in noisy environment[END_REF] and modelbased ones [START_REF] Liu | An algebraic fractional order differentiator for a class of signals satisfying a linear differential equation[END_REF][START_REF] Liu | Robust fractional order differentiators using generalized modulating functions method[END_REF][START_REF] Wei | Non-Asymptotic Pseudo-State Estimation for a Class of Fractional Order Linear Systems[END_REF]. The first class is obtained by truncating an analytical expression. Hence, this generates a truncated term error which can produce an amplitude error (in the vertical sense) and a shifted error (in the horizontal sense) [START_REF] Liu | Fractional order differentiation by integration and error analysis in noisy environment[END_REF]. The second class is obtained from the differential equations of considered signals. They do not introduce any truncated term errors.

Existing fractional order differentiators are usually extensions of integer order differentiators, which generally have one of the following disadvantages:

• they are sensible to noises, such as the well-known Grünwald-Letnikov scheme [START_REF] Concepción | Fractional-order Systems and Controls: Fundamentals and Applications[END_REF], which is the extension of the finite difference scheme and only efficient in noisefree case;

• they produce a truncated term error in the estimated derivatives, such as the model-free fractional order differentiators proposed in [START_REF] Chen | Digital Fractional Order Savitzky-Golay Differentiator[END_REF][START_REF] Liu | Fractional order differentiation by integration and error analysis in noisy environment[END_REF], which are the extension of classical polynomial approximation methods; • they asymmetrically converge, for instance, the fractional order Luenberger-like observer, considered as an model-based fractional order differentiator, is devoted to estimating the pseudo-state variables which are the fractional derivatives of the output [START_REF] Sabatier | On observability and pseudo state estimation of fractional order systems[END_REF]. Moreover, a fractional order observer cannot estimate the fractional derivatives of an arbitrary order.

Among the existing methods, there is a non-asymptotic algebraic method originally introduced by Fliess and Sira-Ramirez for linear identification [START_REF] Fliess | An algebraic framework for linear identification[END_REF]. This method permits to obtain exact algebraic integral formulae for the desired estimators, which can provide estimations in finite-time. It has been shown in [START_REF] Fliess | Analyse non standard du bruit[END_REF] that, thanks to the integral formulae, these estimators exhibit good robustness properties with respect to corrupting noises. By considering these advantages, this method has been extended to design model-free integer order differentiators [START_REF] Mboup | Numerical differentiation with annihilators in noisy environment[END_REF][START_REF] Liu | Error analysis of Jacobi derivative estimators for noisy signals[END_REF] and model-based ones [START_REF] Fliess | Reconstructeurs d'états[END_REF][START_REF] Tian | Fast state estimation in linear time-varying systems: an algebraic approach[END_REF] for non-linear and linear integer order systems, respectively. More concretely, the latter ones estimate the integer order derivatives of the output of the following integer order linear system:

n i=0 a i y (i) (t) = u(t), (1)
where y and u are the output and the input, respectively. Recently, these non-asymptotic and robust differentiators have been extended to fractional order case. In [START_REF] Liu | Fractional order differentiation by integration and error analysis in noisy environment[END_REF], a model-free differentiator has been designed without considering the system model. However, it produces a time-delay in the estimation. In [START_REF] Liu | An algebraic fractional order differentiator for a class of signals satisfying a linear differential equation[END_REF], a model-based differentiator has been proposed to estimate the noninteger order derivatives of the output of the system defined in [START_REF] Podlubny | Fractional differential equations: an introduction to fractional derivatives, fractional differential equations[END_REF]. For the same purpose, another model-based differentiator has been designed by applying the modulating function method in [START_REF] Liu | Robust fractional order differentiators using generalized modulating functions method[END_REF]. The latter two differentiators do not produce any time-delay in estimations. However, they are not applicable to fractional order systems. Having these ideas in mind, the objective of this paper is to extend the algebraic method to design a model-based differentiator to estimate the fractional derivatives of the output of the following fractional order linear system:

N i=0 a i D αi t y(t) = L j=0 b j D γj t u(t), (2)
where D αi t y and D γj t u are the fractional derivatives of the output and the input, respectively.

The contributions of this paper can be outlined as follows:

• a digital model-based fractional order differentiator is introduced, which has the following advantages: (i) it can be used to estimate the fractional derivative of the output of an arbitrary order; (ii) it is given by a new algebraic formula, which does not contain any source of errors in continuous noise free case; (iii) it can provide robust estimations in finite-time in discrete noisy cases, without any truncated term error;

• an error bound is provided for the selection of the design parameter;

• there is no need on the initial conditions of the considered system. This paper is organized as follows: definitions and some useful properties of Riemann-Liouville fractional integrals and derivatives are recalled in Section 2. The main results are given in Section 3. Firstly, the algebraic method is applied to express the fractional derivatives of the output of the considered system by a new algebraic formula in continuous noise free case. Secondly, a digital fractional order differentiator is introduced in discrete noisy cases. Moreover, some error analysis is given. In Section 4, numerical results illustrate the accuracy and the robustness of the proposed fractional order differentiator. Finally, conclusions are summarized in Section 5.

Preliminaries

In this section, definitions and some useful properties of fractional integrals and derivatives are recalled. Moreover, some useful formulae related to the Laplace transform are given.

Riemann-Liouville fractional integrals and derivatives

Let I = [0, h] ⊂ R + 1 , α ∈ R +
, and l = ⌈α⌉, where ⌈α⌉ denotes the smallest integer larger than or equal to α. Then, the following definitions are given. Definition 1 ([1] p. 65) The α order Riemann-Liouville fractional integral of f is defined on ]0, h] as follows:

D -α t f (t) :=      f (t), if α = 0, t 0 κ α (t, τ )f (τ ) dτ , else, (3)
where κ α (t, •) is defined by:

κ α (t, τ ) := 1 Γ(α) (t -τ ) α-1 , (4)
and Γ(•) is the well-known Gamma function.

Definition 2 ([1] p. 68) The α order Riemann-Liouville fractional derivative of f is defined on ]0, h] as follows:

D α t f (t) := d l dt l D α-l t f (t) . ( 5
)
1 In this paper, R+ denotes the set of positive real numbers, R * + (resp. N * ) denotes the set of strictly positive real numbers (resp. strictly positive integers), and Z * -denotes the set of strictly negative integers.

Remark 1 According to (3), if 0 < α < 1, the Riemann-Liouville fractional integrals are defined by improper integrals. Thus, if l = α, the Riemann-Liouville fractional derivatives are also defined by improper integrals in [START_REF] Concepción | Fractional-order Systems and Controls: Fundamentals and Applications[END_REF].

The following formula establishes an additive index law for Riemann-Liouville fractional integrals and derivatives ( [START_REF] Podlubny | Fractional differential equations: an introduction to fractional derivatives, fractional differential equations[END_REF] pp. 71-74):

∀ n ∈ N, α ∈ R, d n dt n {D α t f (t)} = D n+α t f (t). (6)
In the following Lemma, the Leibniz formula for Riemann-Liouville integrals involving a polynomial is recalled. The general Leibniz formula for Riemann-Liouville integrals and derivatives can be found in [START_REF] Miller | An Introduction to the Fractional Calculus and Fractional Differential Equations[END_REF] (p. 75) and [START_REF] Diethelm | The analysis of fractional differential equations: An application-oriented exposition using differential operators of Caputo type[END_REF] (p. 33), respectively.

Lemma 1 [START_REF] Miller | An Introduction to the Fractional Calculus and Fractional Differential Equations[END_REF] (p. 53) Let α ∈ R * + and m ∈ N. Then, the following formula holds:

D -α t [t m f (t)] = m k=0 (-1) k m k Γ(α + k) Γ(α) t m-k D -α-k t f (t). (7)

Laplace transform formulae

Let us assume that the Laplace transform of f exists, which is denoted by f . Then, the Laplace transforms of the Riemann-Liouville integrals and derivatives of f are respectively given on C by ([3] p. 284):

∀ α ∈ R * + , L D -α t f (t) (s) = 1 s α f (s), (8)
L {D α t f (t)} (s) = s α f (s) - ⌈α⌉-1 j=0 s j D α-j-1 t f (t) t=0 , (9)
where s denotes the variable in the frequency domain.

Then, the following lemma can be obtained.

Lemma 2 [START_REF] Liu | An algebraic fractional order differentiator for a class of signals satisfying a linear differential equation[END_REF] Let α ∈ R + and n ∈ N. Then, the following formula holds:

L -1 1 s α f (n) (s) (t) = (-1) n D -α t [t n f (t)] . (10)
Finally, this section ends by giving the following formula:

∀ k ∈ N * , α ∈ R, d k ds k {s α } =        0, if α ∈ N and α < k, (-1) k (k-α-1)! (-α-1)! s α-k , if α ∈ Z * -, Γ(α+1) Γ(α+1-k) s α-k , else. ( 11
)
3 Main results

Problem formulation

Within this framework, a class of fractional order linear systems are considered, which can be described by the following equation ( [START_REF] Concepción | Fractional-order Systems and Controls: Fundamentals and Applications[END_REF], pp. 17-18):

N i=0 a i D αi t y(t) = L j=0 b j D γj t u(t) (12)
on

I = [0, h] ⊂ R + , where y ∈ R is the output, u ∈ R is the input, N ∈ N * , a N ∈ R * , a i ∈ R, for i = 0, . . . , N -1, α i = iα for i = 0, . . . , N , L ∈ N, b j ∈ R,
γ j = jα for j = 0, . . . , L, α = p q , p, q ∈ N * . In order to guarantee the stability of the considered fractional order systems, it is assumed that 0 < α < 2 [START_REF] Victor | Parameter and differentiation order estimation in fractional models[END_REF]. Without loss of generality, p is assumed to be equal to 1 in the sequel. Moreover, u is assumed to be known. Thus, D γj t u can be analytical or numerically calculated. In order to simplify the presentation, it is assumed that b 0 = 1 and b j = 0 for j = 1, . . . , L.

For the systems studied in [START_REF] Chen | Digital Fractional Order Savitzky-Golay Differentiator[END_REF], when dealing with the problem of output regulation (stabilization of the output for example), normally the controller u requires the fractional derivatives of y [START_REF] Bandyopadhyay | Stabilization and Control of Fractional Order Systems: A Sliding Mode Approach[END_REF]. Sometimes, it prefers nonasymptotic one since it might provide estimations in finite-time. This motivates the research work of this paper on non-asymptotic fractional order differentiator. It should be noticed that this paper does not adopt the conventional state-space equation, but uses the inputoutput equation [START_REF] Chen | Digital Fractional Order Savitzky-Golay Differentiator[END_REF] to describe the fractional order linear systems. This choice is naturally due to the fact that the input-output representation is much more suitable by using the algebraic method to realize the desired nonasymptotic fractional order differentiator.

Non-asymptotic fractional order differentiator

As mentioned in Remark 1, the Riemann-Liouville fractional derivatives of y and its Riemann-Liouville fractional integrals of order strictly smaller than 1 are defined by improper integrals. In order to overcome this problem, the algebraic method is applied in this subsection to express the latter by new algebraic formulae in continuous noise free case. These formulae involve proper integrals of y, which can be considered as a low-pass filter in noisy cases [START_REF] Fliess | Analyse non standard du bruit[END_REF].

Firstly, the fractional integrals of order strictly smaller than 1 are studied in the following theorem. The idea is to apply the algebraic method to transform the fractional order differential equation ( 12) into a fractional order integral equation by eliminating the unknown initial conditions. For this purpose, the following notations are introduced:

∀ β d ∈ [0, 1[, ∃ β q ∈ [0, 1 q [ and d ∈ {1, . . . , q}, such that β d = 1 -d q + β q .
Theorem 1 Let y be the output of a fractional order linear system defined by [START_REF] Chen | Digital Fractional Order Savitzky-Golay Differentiator[END_REF]. Then, the β d order Riemann-Liouville fractional integral of y can be given on ]0, h] by a recursive way as follows:

D -β d t y(t) = - N -1 i=0 a i a N D -β i,d t y(t) + m 0 t 0 P d,m (t, τ )y(τ ) dτ + t 0 Q d,m (t, τ )u(τ ) dτ, (13)
where

β i,d = α N -α i +β d = N -i q +β d , m ∈ N * , satisfying ⌈α N ⌉ ≤ m, is a design parameter, m 0 = 0 if m = 0, and m 0 = 1 else, P d,m (t, τ ) = - N β d i=0 a i a N m k=1 (-1) k m k Γ(β i,d + k) Γ(β i,d ) κ β i,d +k (t, τ ) t k - N i=0 a i a N m j=1 m j c i,j × m-j k=0 m -j k Γ(β i,d + j + k) Γ(β i,d + j) (-1) (j+k) t j+k κ β i,d +j+k (t, τ ), (14)
and

Q d,m (t, τ ) = 1 a N m k=0 (-1) k m k Γ(β 0,d + k) Γ(β 0,d ) κ β 0,d +k (t, τ ) t k , ( 15
) N β d = N -1 if β d = 0, N β d = N else, c i,j
are the coefficients of (s αi ) (j) , which are given by [START_REF] Machado | Exploiting sensor redundancy for the calculation of fractional derivatives in the presence of noise[END_REF], and κ β i,d +j+k (t, τ ) can be given by (4).

Proof. In this proof, the algebraic method will be applied, where the following steps are required.

• Step 1. Laplace transform: By applying the Laplace transform to (12) and using (9), we get: [START_REF] Wei | Non-Asymptotic Pseudo-State Estimation for a Class of Fractional Order Linear Systems[END_REF] where ŷ (resp. û) denotes the Laplace transform of ȳ (resp. ū), and ȳ (resp. ū) is defined on R in such a way that ȳ(t) = y(t) (resp. ū(t) = u(t)) if t ∈ I, and ȳ(t) = 0 (resp. ū(t) = 0), else.

N i=0 a i s αi ŷ(s)- N i=1 a i ⌈αi⌉-1 j=0 s j D αi-j-1 t y(t) t=0 = û(s),
• Step 2. Elimination of initial conditions: Since 0 ≤ α 0 < α 1 • • • < α N , if ⌈α N ⌉ ≤ m,
applying m times differentiation with respect to s allows us to eliminate the unknown initial conditions in [START_REF] Wei | Non-Asymptotic Pseudo-State Estimation for a Class of Fractional Order Linear Systems[END_REF]:

û(m) (s) = N i=0 a i d m ds m {s αi ŷ(s)} . (17)
Then, by applying the classical Leibniz formula we get:

û(m) (s) = N i=0 a i s αi ŷ(m) (s) + m 0 N i=0 a i m j=1 m j c i,j s αi-j ŷ(m-j) (s), (18)
where the coefficients c i,j of (s αi ) (j) are given by ( 11), m 0 = 0 if m = 0, and m 0 = 1 else.

• Step 3. Division by s αN +β with 0 ≤ β < 1: In order to apply Lemma 2 to obtain an integral formula when returning back into the time domain, the powers of s should become strictly negative in [START_REF] Fliess | Analyse non standard du bruit[END_REF]. Hence, dividing (18) by s αN +β gives:

û(m) (s) s αN +β = N i=0 a i ŷ(m) (s) s αN -αi+β + m 0 N i=0 a i m j=1 m j c i,j ŷ(m-j) (s) s αN -αi+β+j . (19)
• Step 4. Inverse of Laplace transform: When returning back into the time domain by applying the inverse of the Laplace transform, Lemma 2 gives:

∀ t ∈ ]0, h], N i=0 a i D αi-αN -β t [t m y(t)] = D -αN -β t [t m u(t)] -m 0 N i=0 a i m j=1 (-1) j m j c i,j D αi-αN -β-j t t m-j y(t) . (20)
Thus, the differential equation ( 12) is transformed into an integral equation in [START_REF] Liu | Error analysis of Jacobi derivative estimators for noisy signals[END_REF] without the knowledge on the initial conditions.

• Step 5. Leibniz formula for fractional integrals: By applying Lemma 1 to the two sides of ( 20), we obtain:

on the one hand, N i=0 a i D αi-αN -β t [t m y(t)] = m k=0 (-1) k m k Γ(α N + β + k) Γ(α N + β) t m-k D -αN -β-k t u(t) -m 0 N i=0 a i m j=1 m j c i,j m-j k=0 (-1) (k+j) m -j k × Γ(α N -α i + β + j + k) Γ(α N -α i + β + j) t m-j-k D αi-αN -β-j-k t y(t), (21)
on the other hand, if 0 < β < 1, we have:

N i=0 a i D αi-αN -β t [t m y(t)] = N i=0 a i × m k=0 (-1) k m k Γ(α N -α i + β + k) Γ(α N -α i + β) t m-k D αi-αN -β-k t y(t), (22)
and if β = 0, we have:

N i=0 a i D αi-αN -β t [t m y(t)] = a N t m y(t) + N -1 i=0 a i × m k=0 (-1) k m k Γ(α N -α i + k) Γ(α N -α i ) t m-k D αi-αN -k t y(t). (23)
By regrouping the terms in the double sums in the right sides of ( 22) and ( 23), respectively, the following equation is obtained:

N i=0 a i D αi-αN -β t [t m y(t)] = N i=0 a i t m D αi-αN -β t y(t) + m 0 N β i=0 a i m k=1 (-1) k m k × Γ(α N -α i + β + k) Γ(α N -α i + β) t m-k D αi-αN -β-k t y(t), (24)
where

N β = N if 0 < β < 1, and N β = N -1 if β = 0.
Hence, using ( 24), ( 21) becomes:

N i=0 a i D αi-αN -β t y(t) = m k=0 (-1) k m k Γ(α N + β + k) Γ(α N + β) 1 t k D -αN -β-k t u(t) -m 0 N β i=0 a i m k=1 m k Γ(α N -α i + β + k) Γ(α N -α i + β) (-1) k t k D αi-αN -β-k t y(t) -m 0 N i=0 a i m j=1 m j c i,j m-j k=0 (-1) (k+j) m -j k × Γ(α N -α i + β + j + k) Γ(α N -α i + β + j) 1 t j+k D αi-αN -β-j-k t y(t).
(25)

• Step 6. Recursive algorithm: Recall that α N -α i + β = N -i q +β, for i = 0, . . . , N . If β = β d = 1-d q +β q with 0 ≤ β q < 1 q , for d = 1, . . . , q, we get α N -α i + β = β i,d with β i,d = N -i-d q + 1 + β q .
Thus, according to [START_REF] Butt | Introduction to Numerical Analysis Using MATLAB[END_REF], it yields:

D -β d t y(t) = - N -1 i=0 a i a N D -β i,d t y(t) + 1 a N m k=0 (-1) k m k Γ(β 0,d + k) Γ(β 0,d ) 1 t k D -β 0,d -k t u(t) -m 0 N β d i=0 a i a N m k=1 (-1) k m k Γ(β i,d + k) Γ(β i,d ) 1 t k D -β i,d -k t y(t) -m 0 N i=0 a i a N m j=1 m j c i,j m-j k=0 (-1) (j+k) m -j k × Γ(β i,d + j + k) Γ(β i,d + j) 1 t j+k D -β i,d -j-k t y(t), (26)
where

N β d = N -1 if β d = 0, and N β d = N else.
Consequently, this proof is completed.

Remark 2 It can be seen in the proof of Theorem 1 that the m times differentiation is devoted to eliminating the unknown initial conditions with ⌈α N ⌉ ≤ m. Hence, if the initial conditions are equal to zero, m can be set to any positive integer number. In the next subsection, an error bound depending on m will be provided in discrete noisy cases, then the role of the design parameter m will be studied in the part of numerical simulations.

Thanks to Theorem 1, ∀ β d ∈ [0, 1[, D -β d t
y is exactly given by proper integrals of y and integrals of u using a recursive way. Indeed, in order to calculate D -β d t y by applying ( 13)-( 15), the required orders of the fractional integrals of u and y are given in Tab. 1 and Tab. 2, respectively, for d = 1, . . . , q, where we have:

• for d = 1, D -β1
t y is given by the fractional integrals of y (resp. of u) with orders varying from 1 + β q to α N + β 1 + m (resp. from α N + β 1 to α N + β 1 + m), which are proper integrals;

• for d = 2, . . . , q, D -β d t y is given by D -β b t y for b = 1, . . . , d -1, which are previously given by proper integrals, and by the fractional integrals of y (resp. of u) with orders varying from 1 + According to [START_REF] Concepción | Fractional-order Systems and Controls: Fundamentals and Applications[END_REF], the Riemann-Liouville fractional derivatives of y are defined by taking the integer order derivatives of the Riemann-Liouville fractional integrals of y of order strictly smaller than 1. Inspired by this idea, the following corollary can be obtained using Theorem 1. Before doing so, the following notations are introduced:

β q to α N + β d + m (resp. from α N + β d to α N + β d + m), which are proper integrals (resp. proper integrals if α N + β d ≥ 1, i.e. N ≥ d). D -β d t y D -β 0,d -k t u β1 ∈ [1 -1 q , 1[ {αN + β1, . . . , αN + β1 + m} . . . . . . β d ∈ [1 -d q , 1 -d+1 q [ {αN + β d ,
∀ ν d ∈ R * + , ∃ β q ∈ [0, 1 q [ and d ∈ {1, . . . , q}, such that ν d = l -β d with l = ⌈ν d ⌉ and β d = 1 -d q + β q ∈ [0, 1[. Hence, ∀ d ∈ {1, .
. . , q}, we have: l = ⌈ν d ⌉, i.e. l is independent of the index d.

Corollary 1 Let y be the output of a fractional order linear system defined by [START_REF] Chen | Digital Fractional Order Savitzky-Golay Differentiator[END_REF]. Then, the ν d order Riemann-Liouville derivative of y can be given on ]0, h] by a recursive way as follows:

D ν d t y(t) = d l dt l D -β d t y(t) = - N -1 i=0 a i a N D l-β i,d t y(t) + m 0 d l dt l t 0 P d,m (t, τ )y(τ ) dτ + d l dt l t 0 Q d,m (t, τ )u(τ ) dτ, (27)
where andc i,j are the coefficients of (s αi ) (j) , which are given by [START_REF] Machado | Exploiting sensor redundancy for the calculation of fractional derivatives in the presence of noise[END_REF].

β i,d = α N -α i + β d , m 0 = 0 if m = 0, and m 0 = 1 else, d l dt l t 0 P d,m (t, τ )y(τ ) dτ = - N β d i=0 a i a N m k=1 m k Γ(β i,d + k) Γ(β i,d ) × l k ′ =0 l k ′ (-1) (k ′ +k) (k ′ + k -1)! (k -1)! 1 t k+k ′ D l-k ′ -β i,d -k t y(t) - N i=0 a i a N m j=1 m j c i,j m-j k=0 m -j k Γ(β i,d + j + k) Γ(β i,d + j) × l k ′ =0 l k ′ (-1) (j+k+k ′ ) (k ′ + k + j -1)! (k + j -1)! t k+j+k ′ D l-k ′ -β i,d -j-k t y(t), (28)
d l dt l t 0 Q d,m (t, τ )u(τ ) dτ = 1 a N D l-β 0,d t u(t) + 1 a N m k=1 m k Γ(β 0,d + k) Γ(β 0,d ) × l k ′ =0 l k ′ (-1) (k ′ +k) (k ′ + k -1)! (k -1)! t k+k ′ D l-k ′ -β 0,d -k t u(t), (29)
N β d = N -1 if β d = 0, N β d = N else,
Proof. (27) can be obtained using (6) and Theorem 1. Finally, (28)-( 29) can be obtained by applying the classical Leibniz formula, where two cases of k = 0 and k = 0 are considered in (29).

Remark 3 If q = 1 in (12), the defined systems are integer order ones. On the one hand, in the case of β q = 0, Corollary 1 permits to calculate the non-integer order derivatives of the output of an integer order system, as done in [START_REF] Liu | An algebraic fractional order differentiator for a class of signals satisfying a linear differential equation[END_REF]. Thus, this paper is an extension of [START_REF] Liu | An algebraic fractional order differentiator for a class of signals satisfying a linear differential equation[END_REF]. On the other hand, in the case of β q = 0, Corollary 1 permits to calculate the integer order derivatives of the output of an integer order system, i.e. it is a model-based integer order differentiator. Thus, this paper is also an extension of [START_REF] Tian | Fast state estimation in linear time-varying systems: an algebraic approach[END_REF].

Similar to Theorem 1, ∀ ν d ∈ R * + , D ν d
t y is exactly given by proper integrals of y and a part involving u using a recursive way. Indeed, D ν d t y is obtained by taking the integer order derivative of D -β d t y and calculated by ( 27)-(29). Using Tab. 1 and Tab. 2, the required orders of the fractional integrals and derivatives of u and y are given in Tab. 3 and Tab. 4, respectively, for d = 1, . . . , q, where l = 1 in Tab. 4.

D ν d t y D l-β 0,d t u and D l-k ′ -β 0,d -k t u -ν1 = β1 -l {αN -ν1, . . . , αN + β1 + m} . . . . . . -ν d = β d -l {αN -ν d , . . . , αN + β d + m} . . . . . .
-νq = βql {αNνq, . . . , αN + βq + m} Table 3 Required orders of the fractional integrals and derivatives of u in the formula of D ν d t y. -νq = βq -1 {α1νq, . . . , -ν1, βq . . . , αN -νq}

D ν d t y D 1-β i,d t y -ν1 = βq -1 q {βq, βq + 1 q , . . . ,
D ν d t y D 1-k ′ -β i,d -k t y and D 1-k ′ -β i,d -j-k t y -ν1 = βq -1 q {α1 -ν1 + 1, . . . , αN + β1 + m} . . . . . . -ν d = βq -d q {α1 -ν d + 1, . . . , αN + β d + m} . . . . . . -νq = βq -1 {α1 -νq + 1 . . . , αN + βq + m} Table 4
Required orders of the fractional integrals and derivatives of y in the formula of D

ν d
t y in the case where l = 1.

Then, on the one hand, we have:

∀ ν d ∈ R * + , D ν d t y in- volves D l-β 0,d t u and D l-k ′ -β 0,d -k t u, whose orders vary from α N -ν d to α N + β d + m. If α N -ν d ≥ 1, all the integrals of u in D ν d t y are proper. If 0 < α N -ν d < 1, D ν d
t y contains improper integrals of u. In the two previous cases, if the input u contains noises, the noise effect can be reduced by the integrals considered as low-pass filters

. If α N -ν d = 0, D ν d t y contains u. If α N -ν d < 0, D ν d
t y contains the fractional derivatives of u.

On the other hand, we have:

• for l = 1, d = 1, D ν1 t y contains (i) D -β b t
y for b = 1, . . . , q, with 0 ≤ β b < 1, which are given by proper integrals of y by Theorem 1, (ii) the fractional integrals of y with orders varying from 1 to α N + β 1 + m, which are proper;

• for l = 1, d ∈ {2, . . . , q}, D ν d t y contains (i) the fractional derivatives D ν b t y for b = 1, . . . , d -1, which are previously given by proper integrals of y in Corollary 1, (ii) D -β b t y for b = 1, . . . , q, with 0 ≤ β b < 1, which are given by Theorem 1, (iii) the fractional integrals of y with orders varying from 1 to α N + β d + m; • similarly, for ν d > 1, D ν d t y contains (i) the fractional derivatives D α t y for α = ν d -1 q , ν d -2 q , . . . , 1 qβ q > 0, which are previously given by Corollary 1, (ii) D -β b t y for b = 1, . . . , q, with 0 ≤ β b < 1, which are given in Theorem 1, (iii) the fractional integrals of y with orders varying from 1 to α N + β d + m.

According to the previous analysis, D ν d t y is given using recursive way, which only contains proper integrals of y and a part involving u. By grouping the integrals of y, D ν d t y can be given by an algebraic formula as follows:

D ν d t y(t) = t 0 W m (τ ) y(τ ) dτ + U m (t), (30)
where U m denotes the part involving u, and W m denotes the corresponding function in the integral of y, which is the sum of the associated functions in the proper integrals of y in D α t y for α ∈ {ν d -1 q , . . . , 1 qβ q , -β q , . . . , -β 1 , -1, . . . , -α Nβ d -m}. Hence, W m is bounded on [0, t].

Consequently, the Riemann-Liouville fractional derivative of y of an arbitrary order is exactly given by a new algebraic formula in Corollary 1, i.e. it does not contain any source of errors in continuous noise-free case, and can provide non-asymptotic estimations in finite-time without knowing the initial conditions of the considered system. Moreover, the proper integrals of y can reduce the noise error contribution if the noisy observation of y is used. In the next subsection, a digital fractional order differentiator is introduced in discrete noisy cases, where some error analysis is also given based on (30).

Digital fractional order differentiator in discrete noisy cases

From now on, let y ̟ be a discrete noisy observation of

y on I = [0, h] ⊂ R + : y ̟ (t i ) = y(t i ) + ̟(t i ), (31)
where t i = iT e , for i = 0, 1, . . . , M , with an equidistant sampling period T e = h M . Moreover, the noise {̟(t), t ∈ I} is assumed to be a continuous stochastic process satisfying the following conditions: Note that a zero-mean white Gaussian noise satisfies these conditions.

(
The Riemann-Liouville fractional integrals and derivatives of the output y can be numerically approximated by the following Grünwald-Letnikov scheme [START_REF] Podlubny | Fractional differential equations: an introduction to fractional derivatives, fractional differential equations[END_REF][START_REF] Concepción | Fractional-order Systems and Controls: Fundamentals and Applications[END_REF]:

∀ α ∈ R, D α t y(t) ≈ 1 T α e ⌈ t Te ⌉ j=0 w (α) j y(t -jT e ), (32)
where the binomial coefficients can be recursively calculated as follows:

   w (α) 0 = 1, if j = 0, w (α) j = 1 -α+1 j w (α) j-1 , else. (33)
However, this scheme is not robust against noises.

Different from the Grünwald-Letnikov scheme, according to (30), the formula of D ν d t y contains two parts: the one given by a proper integral of y, and the other one involving u. In order to estimate the Riemann-Liouville fractional derivatives of y in discrete noisy cases, D ν d t y is approximated by following digital fractional order differentiator: for i = 1, . . . , M ,

D ν d m y ̟ (t i ) := T e i j=0 w j W m (t j ) y ̟ (t j ) + Ũm (t i ), ( 34
)
where w i ∈ R + are the weights of a given numerical integration method. According to previous analysis, Ũm (t i ) is calculated in the following way:

• the fractional integrals of u of order larger than 1 are calculated using a numerical integration method;

• the fractional integrals of u of order strictly smaller than 1 and the fractional derivatives of u are numerically calculated using the Grünwald-Letnikov scheme given in (32)-(33).

Consequently, the digital fractional order differentiator D ν d m y ̟ (t i ) can be used to estimate the Riemann-Liouville fractional derivative D ν d t y for any

ν d ∈ R * + ,
where y is the output of a fractional order linear system defined in [START_REF] Chen | Digital Fractional Order Savitzky-Golay Differentiator[END_REF]. It contains three sources of errors:

• the numerical error due to the numerical integration method used to approximate the proper integral involving y;

• the noise error contribution of the following form:

e ̟ m (t i ) := T e i j=0 w j W m (t j ) ̟(t j ); (35)
• the numerical errors in Ũm (t i ), which are due to the numerical integration method and the Grünwald-Letnikov scheme, respectively.

Firstly, it is well known that the numerical errors due to the numerical integration method converge to zero when T e → 0 [START_REF] Ralston | A first course in numerical analysis[END_REF]. Secondly, as shown previously, the Grünwald-Letnikov scheme is used to approximate D α t u with -1 < α in the calculation of Ũm (t i ). It is shown in [START_REF] Diethelm | The analysis of fractional differential equations: An application-oriented exposition using differential operators of Caputo type[END_REF] that in the case of -1 < α < 0 (resp. α > 0), if u is continuous on I (resp. u is ⌈α⌉-times continuously differentiable on I), the produced numerical error converges to zero when T e → 0. Thirdly, if the noise satisfies conditions (C 1 ) -(C 3 ), it can be shown that the noise error contribution e ̟ m (t i ) converges to zero in mean square when T e → 0 (see [START_REF] Liu | An algebraic fractional order differentiator for a class of signals satisfying a linear differential equation[END_REF] for more details). A similar result can also be obtained using non-standard analysis as done in [START_REF] Fliess | Analyse non standard du bruit[END_REF]. Consequently, both the numerical errors and the noise error contribution can be reduced in the proposed digital fractional order differentiator by decreasing the sampling period T e . However, if the computations are performed on a finite precision numerical machine, there are also round-off errors [START_REF] Butt | Introduction to Numerical Analysis Using MATLAB[END_REF]. Hence, the infinite reduction of T e would not lead to arbitrary reduction of estimation errors, since the round-off errors would become too large at some points.

When the sampling period is set, by applying (C 1 )-(C 3 ) and the properties of the mean value and the variance functions to (35), one obtains:

E[e ̟ m (t i )] = 0, (36)
Var

[e ̟ m (t i )] = T 2 e i j=0 w 2 j W 2 m (t j )Var [̟(t j )] ≤ δT 2 e i j=0 w 2 j W 2 m (t j ). ( 37
)
Hence, the noise error contribution e ̟ m (t i ) can be bounded using the Bienaymé-Chebyshev inequality:

∀ γ ∈ R * + , Pr |e ̟ m (t i )| < γ (Var[e ̟ m (t i )]) 1 2 > 1 - 1 γ 2 , (38)
i.e. the probability for |e ̟ m (t i )| to be smaller than γ (Var[e ̟ m (t i )])

1 2 is larger than 1 -1 γ 2 . Thus, the following error bound is deduced from (37) and (38):

|e ̟ m (t i )| pγ < γT e   δ i j=0 w 2 j W 2 m (t j )   1 2 , ( 39
)
where a pγ < b means that the probability for a real number b to be larger than another real number a is equal to

p γ with 1 -1 γ 2 < p γ < 1.
Remark that the value of p γ can be given by the probability density function of the considered noise ̟. In particular, if ̟ is a zero-mean white Gaussian noise, it can be deduced that e ̟ m (t i ) is a normally distributed random number. Thus, according to the three-sigma rule, we have: p 1 = 68.26%, p 2 = 95.44% and p 3 = 99.73%, for γ = 1, 2, 3, respectively.

Thanks to the probability properties of the mean value and the variance functions, and the Bienaymé-Chebyshev inequality, the noise error bound obtained in (39) is sharp, especially in white Gaussian noise cases. Since it depends on the design parameter m, the study of the influence of m on this error bound is useful to deduce the influence of m on the noise error contribution. Consequently, the noise error contribution can be reduced by choosing the value of m which minimizes the noise error bound.

Finally, according to the previous analysis, the following algorithm is required to realize the proposed differentiator D ν d m y ̟ (t i ), which is devoted to applying Theorem 1, Corollary 1, and contains the following steps:

Step 1. find β q ∈ [0, 1 q [ such that ν d = ⌈ν d ⌉ -1 + d q -β q with d ∈ {1, . . .

, q};

Step 2. using β q , calculate β j for j = 1, . . . , q -1;

Step 3. using ( 13)-( 15), calculate D -β1 t y, where all the integrals D -α t y and D -α t u are proper with α ≥ 1, and calculated using a numerical integration method;

Step 4. using ( 13)-( 15), calculate D -βj t y for j = 2, . . . , q in a loop. For each j, the fractional integrals D -α t y (resp. D -α t u) in the expression of D -βj t y are calculated as follows:

• if α ≥ 1, they are calculated using the numerical integration method, • otherwise they are previously calculated (resp. calculated using the Grünwald-Letnikov scheme); Step 5. using ( 27)-(29), calculate D αi t y with α i = ν di q , for i = (ν d + β q )q -1, . . . , 0 in a loop. For each i, the fractional integrals and derivatives D -α t y and D -α t u in the expression of D αi t y are calculated using a similar way as done in Step 2.

Simulation results

In order to show the accuracy and the robustness of the proposed fractional differentiator, some numerical results are presented in this section.

Example 1: In most cases, neither the fractional derivatives of a function nor the solution of a fractional order differential equation can be analytically calculated. In order to obtain the exact value of the sought fractional derivative, the following academic model is considered in this example: ∀ t ∈]0, 30],

4 i=0 a i D i 3 t y(t) = 4 i=0 a i D i 3 t u(t), (40)
where a i = 1 for i = 0, . . . , 4, u(t) = y(t) = y 1 (t)+ y 2 (t), y 1 (t) = sin(ω 1 t) and y 2 (t) = cos(ω 2 t) with ω 1 = 1.2 and ω 2 = 0.6. The analytical Riemann-Liouville fractional derivatives of y 1 and y 2 can be found in [START_REF] Miller | An Introduction to the Fractional Calculus and Fractional Differential Equations[END_REF].

D 1 2
t y will be estimated in the discrete noisy case with t i = iT e ∈ I = [0, 30], T e = 0.01, for i = 0, . . . , 3000, and y ̟ (t i ) = y(t i ) + δ̟(t i ), where the noise {̟(t i )} is simulated from a zero-mean white Gaussian iid sequence, and the value of δ is adjusted such that the signal-tonoise ratio SNR = 10 log 10

|y ̟ (ti)| 2 |δ̟(ti)| 2
is equal to SNR = 15 dB (see [START_REF] Haykin | Signals & Systems[END_REF] for this well known concept in signal processing). The original output y and the discrete noisy output y ̟ are shown in Fig. 1. Then, the fractional order differentiator obtained in (34) is applied, where the trapezoidal numerical integration method is considered. Firstly, it is shown how to choose the value of the design parameter m using the noise error bound obtained in (39) with γ = 2. Since the value of ν d is set in Corollary 1 with ν d = 1 2 and β q = 2 3 -1 2 , the expression of D

1 2
t y can be obtained. Thus, the expression of W m can be deduced in (30). Then, by taking different values of m, the values of the error bound obtained in (39) can be calculated at each t i for i = 0, . . . , 3000. Using such a way, the variation with respect to m of the noise error bound can be observed in Fig. 2. According to Fig. 2, the influence of m on the noise error contribution can be deduced. Secondly, in order to compare with the proposed method the fractional order Legendre differentiator proposed in [START_REF] Liu | Fractional order differentiation by integration and error analysis in noisy environment[END_REF] is applied. The latter is a model-free differentiator and obtained by a polynomial approximation method. Hence, it contains a truncated term error which can be significantly reduced by admitting a time-delay. Consequently, this differentiator depends on two design parameters: the order of the used polynomial N and the timedelay ϑ. According to the analysis done in [START_REF] Liu | Fractional order differentiation by integration and error analysis in noisy environment[END_REF], the design parameters are set in this example as follows: ϑ = 0. t y are shown in Fig. 3. It can be seen that different from the fractional order Legendre differentiator the proposed differentiator does not produce a time-delay. The reason is that the proposed method is based on the system model. Finally, the corresponding absolute estimation errors are shown in Fig. 4, where the one for the fractional order Legendre differentiator is obtained by shifting the estimation to avoid the time-delay. The latter principally contains the amplitude error and the noise error contribution.

Example 2: In this example, a physical model is considered, which is introduced by a fractional order differential equation [START_REF] Torvik | On the appearance of the fractional derivative in the behavior of real materials[END_REF]. Let us consider a rigid plate immersed in a Newtonian fluid of infinite extent and connected by a massless spring to a fixed point. If a force is applied to the plate, the dynamic of the plate can be represented by t y and its estimations obtained by the proposed method and the fractional order Legendre differentiator. t y produced by the proposed method and the shifted fractional order Legendre differentiator.

the following fractional order differential equation [START_REF] Torvik | On the appearance of the fractional derivative in the behavior of real materials[END_REF]:

∀ t > 0, 4 i=0 a i D i 2 t y(t) = u(t), (41)
where a 4 = M p , a 3 = 2S √ µρ, a 2 = a 1 = 0, a 0 = K, M p is the mass of the plate, S is the area of the plate, µ is the fluid viscosity, ρ is the fluid density, K is the stiffness of the spring, y is the displacement of the plate, and u is the force. Moreover, this plate-fluid system is assumed to be initially in an equilibrium state with y(0) = 0 and ẏ(0) = 0.

The Riemann-Liouville fractional derivatives of the output y will be estimated in the following situation with t i = iT e ∈ I = [0, 45], T e = 0.01, for i = 0, . . . , 4500:

• the output y bias is generated by the Grünwald-Letnikov scheme given in (32)-( 33 t y obtained by the Grünwald-Letnikov scheme in discrete noise-free case and the proposed method in the discrete noisy case, respectively.

• the parameters a j for j = 0, 3, 4, are known with errors: ãj (t i ) = a j + ̟ j (t i ), where {̟ j (t i )} are simulated from the uniform iid sequences with ̟ j (t i ) ∈]δ j , δ j [ and δ j = 1 20 a j ; • the input and output are corrupted by a zero-mean white Gaussian noise with SNR = 15 dB: u ̟ (t i ) = u(t i ) + δ u ̟ u (t i ) and y ̟ (t i ) = y bias (t i ) + δ y ̟ y (t i ). t y obtained by the Grünwald-Letnikov scheme in discrete noise-free case and the proposed method in the discrete noisy case, respectively. t y obtained by the Grünwald-Letnikov scheme in discrete noise-free case and the proposed method in the discrete noisy case, respectively.

The original output y and the discrete noisy output y ̟ are shown in Fig. 5.

Firstly, using a similar way as done in Example 1, the values of the design parameter m are set. For this purpose, D 1 2 t y is considered as an example. Then, the values of the corresponding error bound are calculated at each t i for i = 0, . . . , 4500, and the variation with respect to m is shown in Fig. 6, where the values of the parameters a j for j = 0, 3, 4, are used. Hence, according to Fig. 6, the following values of m are chosen: m = 0 for t i ∈ [T e , 2.9], m = 4 for t i ∈ [2.9 + T e , 8.7], and m = 10 for t i ∈ [8.7 + T e , 45].

Secondly, since the analytical Riemann-Liouville fractional derivatives of the output are unknown, the Grünwald-Letnikov scheme is applied in the noise-free case to verify the efficiency of the proposed fractional order differentiator. The obtained estimation of D 1 2 t y is shown in Fig. 7. It can be seen that the proposed method is robust. Hence, it is interesting to provide some robustness analysis on the proposed method in the considered example.

• By taking δ y ̟ y and the expression of W m in (35), the error due to the noisy output can be calculated at each t i using different values of m chosen previously. Moreover, the corresponding noise error bound can also be calculated using (39) with γ = 2. The obtained results are shown in Fig. 8, where it can be seen that when the time t becomes large, the noise error bound is increasing with t. Hence, the noise error can also be increasing with time. In order to solve this problem, on the one hand, according to Fig. 6, the value of m should be increased to reduce the noise error bound; on the other hand, according to the analysis done in Subsection 3.3, the value of sampling period T e should be reduced.

• According to the analysis done in Subsection 3.2, the expression of D ν d t y contains a proper integral of y and a part involving u in (30). Moreover, if ν d < α N , the part involving u is an integral of u, either proper or improper. In this example, α Nν d = 2 -1 2 = 3 2 > 1, the part involving u is a proper integral of u. This is why the noise effect in the input is reduced in the estimation of D 1 2 t y, as the one in the output.

• Assume ãj = a j + ̟ j for j = 0, . . . , N , with |̟ j | ≤ δ j < |a j |. It can be seen in ( 27)-(29) that the formula of D ν d t y is linear with respect to 1 aN and aj aN . Consequently, the error due to ̟ j in D ν d t y depends on the relative errors in 1 aN and aj aN , which can be bounded as follows: (43)

|a N | 1 ãN - 1 a N = -̟ N a N + ̟ N ≤ δ N |a N | -δ N , (42)
Hence, these relative errors are increasing with respect to δ j . In this example, δ j = 1 20 a j is chosen. This is why the obtained estimation is robust. Although there are multiplicative noises in the estimation due to the noises in the output and the input, since it is shown previously that the noise effect in the output and the input can be significantly reduced, the error due to ̟ j plays a more important role.

• If there is a disturbance in the model, it plays as a part of the input and produces a bias term to the output. Then, using (34) the error due to this bias term can be given as follows:

e bias (t i ) := T e i j=0 w j W m (t j ) (y bias (t j )y(t j )) . (44) Hence, if the fractional derivative of the bias term y biasy exists, e bias its estimation. As shown in Subsection 3.2, W m is bounded. If the bias term is also bounded, e bias is bounded. Thus, e bias depends on how the output is biased, which depends on the structure of the model and the source of the disturbance. In this example, the disturbance ̟ d (t i ) = 0.5 sin(5t i ) does not significantly the output. This is why the obtained estimation is robust.

Consequently, it can be concluded that the proposed method can well cope with the corrupting noises in the output (resp. in the input if ν d < α N ). Moreover, if the errors in the model parameters and the bias term in the output are not relatively large, robust estimations can also be obtained by the proposed method.

Finally, D 0.3 t y and D 1.3 t y are estimated by taking β q = 0.2 in Corollary 1. The same values of m are chosen as done in the estimation of D

Conclusions

In this paper, a non-asymptotic fractional order differentiator was introduced to estimate the Riemann-Liouville fractional derivatives of the output of a class of fractional order linear systems in discrete noisy environment. Firstly, by applying the algebraic method and the Leibniz formula for Riemann-Liouville fractional integrals, the fractional integral of order strictly smaller than 1 of the output was expressed by an integral formula. Then, a new algebraic formula for the fractional derivative of an arbitrary order was deduced by taking the integer order derivative of the obtained formula for fractional integrals. This formula does not contain any source of errors in continuous noise-free case. Hence, it can provide estimations in finite-time without knowing the initial conditions of the considered systems. Secondly, a digital fractional order differentiator was introduced in discrete noisy cases. Then, some error analysis was given, where an error bound useful for the selection of the design parameter was provided. Finally, numerical results were given to show the accuracy and the robustness of the proposed fractional order differentiator.

  β1, . . . , αN -ν1} . . . . . . -ν d = βq -d q {α1ν d = -ν d-1 , . . . , -ν1, . . . , αNν d } . . . . . .

  Hence, in order to reduce the noise error contribution, the following values of m are chosen: m = 2 for t i ∈ [2T e , 1.95], m = 4 for t i ∈ [1.95 + T e , 3.7], m = 8 for t i ∈ [3.7 + T e , 5.5], m = 10 for t i ∈ [5.5 + T e , 7.2], and m = 14 for t i ∈ [7.2 + T e , 30].

Fig. 1 .Fig. 2 .

 12 Fig. 1. Example 1. The output and its discrete noisy observation.

Fig. 3 .

 3 Fig. 3. Example 1. The exact value of D 1 2

Fig. 4 .

 4 Fig. 4. Example 1. The estimation errors of D 1 2

Fig. 5 .Fig. 6 .

 56 Fig. 5. Example 2. The output and its discrete noisy observation.

Fig. 7 .

 7 Fig. 7. Example 2. The numerical calculation and the estimation of D 1 2

Fig. 8 . 2 t

 82 Fig. 8. Example 2. The error due to the noisy output in the estimation of D 1 2t y obtained by the proposed method and the corresponding error bound.

Fig. 9 .

 9 Fig. 9. Example 2. The numerical calculation and the estimation of D 0.3t y obtained by the Grünwald-Letnikov scheme in discrete noise-free case and the proposed method in the discrete noisy case, respectively.

3

 3 Fig. 9. Example 2. The numerical calculation and the estimation of D 0.3t y obtained by the Grünwald-Letnikov scheme in discrete noise-free case and the proposed method in the discrete noisy case, respectively.

Fig. 10 .

 10 Fig. 10. Example 2. The numerical calculation and the estimation of D 1.3t y obtained by the Grünwald-Letnikov scheme in discrete noise-free case and the proposed method in the discrete noisy case, respectively.

3

 3 Fig. 10. Example 2. The numerical calculation and the estimation of D 1.3t y obtained by the Grünwald-Letnikov scheme in discrete noise-free case and the proposed method in the discrete noisy case, respectively.

  ja j ̟ N a j (a N + ̟ N ) ≤ |a N |δ j + |a j |δ N |a j |(|a N |δ N ) .

1 2 t

 12 y. The obtained estimations are shown in Fig.9and Fig.10, respectively.

  C 1 ) : for any t 1 , t 2 ∈ I, t 1 = t 2 , ̟(t 1 ) and ̟(t 2 ) are independent; (C 2 ) : the mean value function of {̟(t), t ∈ I} denoted by E[•] is equal to zero; (C 3 ) : the variance function of {̟(t), t ∈ I} denoted by Var[•] is bounded on I, i.e. ∃ δ ∈ R

+ , ∀ t ∈ I, Var[̟(t)] ≤ δ.

  8, N = 6 for t i ∈ [0.9, 8], N = 12 for t i ∈ [8 + T e , 18], and N = 20 for t i ∈ [18 + T e , 30].

		1
	The obtained estimations of D	2

2 ti y(t i )+a 0 y(t i )+̟ d (t i ) = u(t i ),where a

= 1, a 3 = a 0 = 0.5, u(t i ) = 1 for 0 ≤ t i ≤ 1, u(t i ) = 0 for t i > 1, and ̟ d (t i ) = 0.5 sin(5t i ) is a disturbance which biases the output;


Francisco Javier Bejarano
email: javbejarano@yahoo.com.mx

Thierry Floquet
email: thierry.floquet@ec-lille.fr

Wilfrid Perruquetti
email: wilfrid.perruquetti@inria.fr

Gang Zheng
email: gang.zheng@inria.fr

Observability and Detectability of Singular Linear Systems with Unknown Inputs $

Keywords: Singular systems, strong detectability, strong observability, algebraic observability

In this paper the strong observability and strong detectability of a general class of singular linear systems with unknown inputs are tackled. The case when the matrix pencil is non-regular is comprised (i.e., more than one solution for the differential equation is allowed). It is shown that, under suitable assumptions, the original problem can be studied by means of a regular (non-singular) linear system with unknown inputs and algebraic constraints. Thus, it is shown that for purposes of analysis, the algebraic equations can be included as part of an extended system output. Based on this analysis, we obtain necessary and sufficient conditions guaranteing the observability (or detectability) of the system in terms of the zeros of the system matrix. Corresponding algebraic conditions are given in order to test the observability and detectability. A formula is provided that expresses the state as high order derivative of a function of the output, which allows for the reconstruction of the actual state vector. It is shown that the unknown inputs may be reconstructed also.

Introduction

The problem of designing an observer for a multi-variable linear system partially driven by unknown inputs (UI) has been widely studied [START_REF] Darouach | Full-order observers for linear systems with unknown inputs[END_REF]; [START_REF] Floquet | A sliding mode approach of unknown input observers for linear systems[END_REF]; [START_REF] Floquet | On sliding mode observers for systems with unknown inputs[END_REF]; [START_REF] Guan | A novel approach to the design of unknown input observers[END_REF]). Such observers can be of important use for systems subject to disturbances or with inaccessible inputs, or when dealing with the problem of fault diagnosis. Observability and the problem of observer design have been widely studied for singular systems with perfectly known model. The problems of solvability, controllability and observability were studied in [START_REF] Yip | Solvability, controllability, and observability of continuous descriptor systems[END_REF]. There, the observability analysis is addressed and algebraic characterizations were found. By using the distributional framework, in [START_REF] Cobb | Controllability, observability, and duality in singular systems[END_REF] the algebraic duality between controllability and observability is proven. Necessary and sufficient conditions allowing for the design of a Luenbergerlike observer were found in Paraskevopoulos and Koumboulis (1992). In the three previous mentioned works, it was considered that the system has a regular matrix pencil which entails a unique state solution. Without any particular assumption over the matrix pencil of the system, the casual observability, which does not allow to use neither the derivatives of the input nor the derivatives of the output, is studied in Hou and Müller (1999a). The same authors suggest an observer design in Hou and Müller (1999b). In that work it is shown that by allowing the derivatives of the input and output to be involved in the observer (called it there as a generalized observer), detectability is enough for the convergence of the observation error. A reduced order observer is designed in [START_REF] Darouach | Design of observers for descriptor systems[END_REF]. In [START_REF] Darouach | Observers for a class of nonlinear singular systems[END_REF] an observer for nonlinear singular systems is proposed. In spite of the extended literature regarding the observability analysis and synthesis of singular systems, there exist few results dealing with such problems when the system contains UI. In Paraskevopoulos et al. (1992), the observer design problem is considered for singular linear systems with UI (SLSUI) and necessary and sufficient conditions are given for the design of a Luenberger-like observer. In [START_REF] Darouach | Reduced-order observer design for descriptor systems with unknown inputs[END_REF], a reduced order observer is proposed. Under some regularity conditions, the observer design is studied in [START_REF] Chu | Disturbance decoupled observer design for descriptor systems[END_REF]. Meanwhile, in Koenig (2005) a proportional multiple-integral observer is proposed. Using the graph-theory approach, observability conditions are found in [START_REF] Boukhobza | State and input observability for structured linear systems: A graph-theoretic approach[END_REF].

In this note, the observability problem of a general class of singular linear systems with UI is studied. The system is not required to have a regular matrix pencil. We obtain necessary and sufficient conditions for the strong observability and strong detectability. We show that the reconstruction of the state can be carried out by a formula that expresses the actual state as a high order derivative of a function of the output. The manuscript also includes Section 2, where the system is described and the strong observability and strong detectability are defined. In Section 3, we obtain the necessary and sufficient conditions which allow for the reconstruction of the state vector. Section 4 deals with an explicit formula that allows for the reconstruction of the state vector. The finite time reconstruction (observability) is considered in 4.1. The procedure for the asymptotic reconstruction (detectability) is given in 4.2. A summarized algo-rithm explaining all the estimation procedure is presented in 4.3. To reenforce the theoretical results, we present an example with simulations in Section 5. The following notation will be used throughout the paper. For a matrix X, we denote by X ⊥ a full row rank matrix such that X ⊥ X = 0, and by X ⊥⊥ a full row rank matrix such that rank X ⊥⊥ X = rank X (then the matrix X ⊥ T X ⊥⊥ T T is nonsingular). The Moore-Penrose pseudo-inverse matrix of X is denoted by X + . The rank of X is denoted by means of ρ X . By • , we mean the Euclidean norm. C -denotes de set of complex numbers with strictly negative real part. I r is the identity matrix of dimension r by r. 0 r×s is the zero matrix of dimension r by s. For the limit from above, x (0 + ) = lim t→0 + x (t).

System Description and Problem Formulation

Let us consider the SLSUI governed by the following equations

Σ : E ẋ (t) = Ax (t) + Dµ (t) y (t) = Cx + Fµ (t) , (1)
where x (t) ∈ R n is the state vector, y (t) ∈ R p is the system output, and

µ (t) ∈ R m is the unknown input vector. Matri- ces E, A ∈ R n×n , D ∈ R n×m , C ∈ R p×n , and F ∈ R p×m are all constant.
The matrix E is assumed to be singular. Given a state x 0 ∈ R n and a function µ (t), we denote by x µ (x 0 , t) the state of Σ at time t which results from taking the initial condition equal to x 0 and the input vector is equal to µ. Therefrom, in a straightforward manner we define the output y µ (x 0 , t) = Cx µ (x 0 , t) + Fµ (t). We are interested in the reconstruction of the (non-impulsive) trajectory of state vector x (t) given the output information y (τ) τ∈[0,t] . System Σ is not assumed to have a regular pencil [START_REF] Kaczorek | Polynomial and Rational Matrices Applications in Dynamical Systems Theory[END_REF]), i.e., it is allowed that det(sE -A) = 0 for all s ∈ C (then x µ (x 0 , t) may have more than a solution). Nevertheless, µ (t) must be so that x (t) be piecewise continuous for all t > 0; however, an impulse may occur at t = 0. In order to give algebraic conditions allowing the reconstruction of x (t), we consider the following definitions, which are based on classical definitions for linear time invariant systems (see, e.g. [START_REF] Trentelman | Control theory for linear systems[END_REF]).

Definition 1 (Strong observability). The system Σ is strongly observable (SO) if for all x 0 ∈ R n and for every input function µ, the following implication is satisfied y µ (x 0 , t) = 0 ∀t > 0 implies x 0 + = 0.

(2)

Definition 2 (Strong detectability). The system Σ is strongly detectable (SD) if for all x 0 ∈ R n and for every input function µ, the following implication holds

y µ (x 0 , t) = 0 ∀t > 0 implies lim t→∞ x µ (x 0 , t) = 0. ( 3
)
It is clear that strong observability is a necessary condition to reconstruct the entire trajectory of the state x (t). Indeed, let us suppose that Σ is not strongly observable, then it means that there exist μ and x0 such that y μ ( x0 , t) = 0 ∀t > 0, but x (0 + ) 0. Then, since we assume that x (t) is piecewise continuous, then x µ=0 (0, t) ≡ 0 and x μ ( x0 , t)

0 in an open interval, however, both yield a system output identically equal to zero. Thereby, it would be impossible to reconstruct the entire state trajectory. Below it would be proven that SO is a structural necessary and sufficient condition for the reconstruction in finite time of x (t). Analogously, it will be shown that SD is a structural necessary and sufficient condition for the asymptotic reconstruction of x (t).

Observability Analysis

Since E is singular, there exist non-singular matrices T ∈ R n×n and S ∈ R n×n such that E can be transformed as follows 1

T ES = I ρ E 0 0 0 (4)
Thus, let us define the vector z

:= z T 1 z T 2 T = S -1 x,
where

z 1 ∈ R ρ E and z 2 ∈ R n-ρ E .
In these new coordinates, Σ can be rewritten as follows

Ψ : T ES ż (t) = T AS z (t) + T Dµ (t) y (t) = CS z (t) + Fµ (t) . ( 5
)
In view of (4), Ψ takes the following form

ż1 (t) = T 1 AS 1 z 1 (t) + T 1 AS 2 z 2 (t) + T 1 Dµ (t) , (6a)
0 = T 2 AS 1 z 1 (t) + T 2 AS 2 z 2 (t) + T 2 Dµ (t) , (6b)
y = CS 1 z 1 (t) + CS 2 z 2 (t) + Fµ (t) . ( 6c
)
where S 1 and S 2 matrices arise from the following partition of S , S = S 1 S 2 with S 1 ∈ R ρ E ×ρ E and S 2 ∈ R ρ E ×n-ρ E . Analogously, T 1 and T 2 matrices come from the partition

T T = T T 1 T T 2 with T 1 ∈ R ρ E ×n and T 2 ∈ R n-ρ E ×n . It is clear that Σ is SO (resp. SD) 2 if
, and only if, Ψ is SO (resp. SD). Below we will see that a simple manner to study the observability of Ψ, and by extension of Σ, is by considering (6b) as part of the system output of a new pseudo system and considering z 2 as part of the vector of UI. Indeed, let us define the system Φ by means of the following equation,

Φ : ż1 (t) = Āz 1 (t) + Dv (t) ȳ (t) = Cz 1 (t) + Fv (t) , (7)
where

v (t) ∈ R n-ρ E +m , ȳ (t) ∈ R n-ρ E +p
and the matrices Ā, D, C, and F are defined as follows

Ā = T 1 AS 1 , D = T 1 AS 2 T 1 D , C = T 2 AS 1 CS 1 , F = T 2 AS 2 T 2 D CS 2 F . ( 8
)
1 We might select S = S 1 S 2 to be non-singular and so that im S 2 = ker E. Thus, ES = ES 1 0 and rank ES 1 = rank E. Then a non-singular matrix T might be selected as

T T = T T 1 T T 2 so that T 1 ES 1 = I and T 2 ES 1 = 0, one possibility is to select T 1 = (ES 1 ) + = (ES 1 ) T (ES 1 ) -1 (ES 1 ) T . 2 Respectively strongly detectable.
It is clear by ( 6) that Φ looks like system Ψ. In general, they do not represent identical systems. However, both systems are identical if these both identities hold: ȳT = 0 T y T and v T (t) = z T 2 (t) µ T (t) . In the next theorem it is claimed that the fulfillment of the SO (resp. SD) of Σ is equivalent to the fulfillment of the SO (resp. SD) of Φ (condition needed for the reconstruction of z 1 ) plus a rank condition (required for the reconstruction of z 2 ).

Theorem 1. System Σ is SO (resp. SD) if, and only if, Φ is SO (resp. SD) and the following rank condition holds

ρ B = n -ρ E + rank D F , where B := D F . (9)
Furthermore, this equivalence is independent of the choice of T and S .

Proof.

Firstly, notice that, since rank D T F T = rank (T D) T F T , the fulfillment of ( 9) is equivalent to say that B z T 2 µ T T = 0 only if z 2 = 0. Necessity. Let Σ be SO (resp. SD). Let us suppose that, for an input v and a state z 1 , ȳv (z 1 , t) = 0. For the case when z 2 (t) and µ (t) satisfy the identity z T 2 (t) µ T (t) = v T (t), we obtain that Ψ and Φ represent the same system. Thus, defining x 0 = S z (0) we conclude that y µ (x 0 , t) = 0 for all t > 0. Now, since, by assumption, (2) (resp. (3)) holds, x (0 + ) = 0 (resp. x (t) converges to zero), which in turn implies that z (0 + ) = 0 (resp. z (t) converges to zero), in particular z 1 (0 + ) = 0 (resp. lim t→∞ z 1 (t) = 0), i.e. Φ is SO (resp. SD). Now, assume that (9) does not hold. Then, there exists a vector v which can be divided as

v T = v T 1 v T 2 (v 1 ∈ R n-ρ E , v 2 ∈ R m
) so that Bv = 0 and v 1 0. By selecting z 2 (t) = v 1 and µ (t) = v 2 , and z (0 + ) = 0, eq. ( 6) is fulfilled and y (t) = 0 for all t > 0. Therefore

x (t) = S z (t) = S 0 v T 1 T = const 0.
That is, in such a case Σ is not SO (resp. Σ is not SD).

Sufficiency. Firstly, let us notice that y µ (x 0 , t) ≡ 0 implies ȳv (z 10 , t) ≡ 0 for some z 10 and v. Indeed, let us suppose that y µ (x 0 , t) = 0 for a state x 0 ∈ R n and an input function µ. Then, the function z µ (z 0 , t) = S -1 x µ (x 0 , t) satisfies ( 6), specifically the algebraic constraint in (6b) is fulfilled. Therefore, we conclude that, for system Φ, ȳv (z 10 , t) = 0, for all t > 0, with v as the extended vector of z 2 and µ. Furthermore, it is known, from linear system theory, that an input v (t) zeroing the output must have the form v (t) = K * z 1 (t) + Lw (t) (t > 0), for a particular matrix K * , a matrix L such that BL = 0, and a function w (t).

Thus, assuming that Φ is SO (resp SD), y µ (x 0 , t) = 0 for all t > 0, implies, since also ȳv (z 10 , t) = 0, for all t > 0, that z 1 (t) ≡ 0 (resp. z 1 (t) converges to zero) and Bv = 0 (resp. Bv = K * z 1 (t)). The previous identity, assuming that (9) is true, implies that z 2 (t) ≡ 0 (resp. z 2 (t) converges to zero). Therefore, in such a case, z (t) ≡ 0 and so x (t) ≡ 0. Therefore, we conclude that Σ is SO (resp. SD).

The independence from T and S is straightforward since we have proven Theorem 1 for arbitrary T and S satisfying (4).

As for Σ, we could expect that SO and SD can be completely characterized by the five-tuple (E, A, C, D, F). Indeed, let R (s) be the so-called system matrix of Σ, i.e.,

R (s) = sE -A -D C F , s ∈ C. We say that s 0 ∈ C is a zero of Σ if ρ R(s 0 ) < n + rank D F . Let
σ z (Σ) be defined as the set of zeros of Σ. Let us characterize SO and SD in terms of the zeros of Σ.

Corollary 2. System Σ is SO (resp. SD) if, and only if,

σ z (Σ) = ∅ (resp. σ z (Σ) ⊂ C -). Proof. Let us define Q (s) = sI -Ā - D C
F , which is the system matrix of Φ. In view of ( 4) and ( 8), we obtain that,

T 0 0 I R (s) S 0 0 I = Q (s) (10)
Thus, we deduce that ρ R(s) = ρ Q(s) . Then, Corollary 2 follows from the claims of Theorem 1 and the fact that Φ is SO (resp. SD) if, and only if,

σ z (Φ) = ∅ (resp. σ z (Φ) ⊂ C -) 3 .

Algebraic Observability

As we might expect SO coincides with algebraic observability (AO): we say that Σ is AO if x can be expressed as an algebraic function of y and a finite number of its derivatives (see, e.g. [START_REF] Diop | Nonlinear observability, identifiability, and persistent trajectories[END_REF]). Let M k (k ≥ 1) be the matrices obtained by the following recursive algorithm (see, [START_REF] Molinari | A strong controllability and observability in linear multivariable control[END_REF]),

M k+1 = N ⊥⊥ k+1 N k+1 , M 1 = F⊥ C ⊥⊥ F⊥ C, N k+1 = T k M k Ā C , T k = M k D F ⊥ . ( 11
)
Let us denote by l, the smallest integer such that ρ M l = ρ M l+1 . For our purposes, we point out that Φ is SO if, and only if, ρ M l = ρ E . For the case of SD we have to work a bit more with system Φ. Indeed, let us assume that ρ M l < ρ E . Let V be a full column rank matrix so that M l V = 0. There exists a pair of matrices Q and K * such that

ĀV + DK * = V Q and CV + FK * = 0. ( 12
)
From ( 12), it is clear that the im Ā + DK * V + V ⊂ im V and C + FK * V + V = 0. We can define a non-singular matrix P of dimension ρ E as P = M l V + , P -1 = M + l V , where

V + and M + l are the Moore-Penrose pseudo-inverse of V and M, respectively:

V + = V T V -1 V T and M + l = M T l M l M T l -1
. By defining the vectors w 1 = M l z 1 and w 2 = V + z 1 , we have that

z 1 = M + l w 1 + Vw 2 .
System Φ in these new coordinates can be rewritten as follows:

ẇ1 = Ā1 w 1 + D1 (v -K * w 2 ) , (13a) ẇ2 = Ā2 w 1 + Ā3 w 2 + D2 (v -K * w 2 ) , (13b) ȳ = C1 w 1 + F (v -K * w 2 ) , (13c)
where

Ā1 = M l Ā + DK * V + M + l , D1 = M l D, Ā2 = V + Ā + DK * V + M + l , D2 = V + D, Ā3 = V + Ā + DK * V + V, C1 = CM + l . (14)
Thus, as for SD, it is known that system Φ is SD if, and only if, rank D1 F = rank D F and Ā3 is a Hurwitz matrix (see, e.g. [START_REF] Bejarano | Unknown input and state estimation for unobservable systems[END_REF]). Coming back to system Φ. Define ξ 1 := F⊥ C ⊥⊥ F⊥ ȳ = M 1 z 1 , with M 1 defined as in ( 11). Let us derive the vector ξ 1 :

ξ1 (t) = M 1 Āz 1 (t) + M 1 Dv (t) . ( 15
)
Let us define a new vector ξ 2 as follows

ξ 2 := N ⊥⊥ 2 T 1 ξ1 ȳ (t) , (16)
with N ⊥⊥ 2 and T 1 defined by ( 11). Thus, taking into account ( 7), (15), and (11), we have that

d dt J 2       ȳ t t 0 ȳ (τ) dτ       = ξ 2 = M 2 z 1 (t) , t > t 0 ≥ 0, ( 17
)
where

J 2 = N ⊥⊥ 2 T 1 J 1 0 0 I p , J 1 = F⊥ C ⊥⊥ F⊥ .
In the first identity of (17), we take outside the differential operator from ( 16) and use the definition of ξ 1 . Thus, we can follow an iterative procedure to obtain the following set of equations, for k ≥ 1,

d k dt k J k+1              ȳ (t) . . . t t 0 • • • τ 2 t 0 ȳ (τ 1 ) dτ 1 • • • dτ k              = M k+1 z 1 , ( 18
)
where M k+1 is defined by ( 11), and J k+1 is defined by the following recursive algorithm, for k ≥ 1,

J 1 = F⊥ C ⊥⊥ F⊥ , J k+1 = N ⊥⊥ k+1 T k J k 0 0 I p . ( 19
)
Thus M k z 1 is expressed by a high order derivative of a function of y (t). In such a way a real-time differentiator could be used, two of them frequently used due to their finite time convergence can be found in [START_REF] Levant | High-order sliding modes: differentiation and outputfeedback control[END_REF] and Mboup et al. (2009a). For instance if rank M l = ρ E , then z 1 is algebraically observable, i.e. it could be reconstructed by using a real-time differentiator.

In order to match system Σ with system Φ, from now on, we define ȳ = 0 1×n-ρ E y T ∈ R p ( p := n -ρ E + p), and

v (t) = z T 2 (t) µ T (t) T ∈ R q (q = n -ρ E + m)
, then in view of (6), equations ( 5) and ( 7) are identical. Below, we consider two cases: when Σ is SO and when it is SD, but not SO. Of course, since Φ is a standard linear system, there might be other methods, besides the one proposed below, that might be used to carry out the algebraic reconstruction of the state.

Finite time reconstruction

Let us consider that system Σ is SO. Then, the reconstruction of entire state vector x (t) can be achieved in finite time: by means of an algebraic formula. Let us proceed in the following way. Since Φ is SO, ρ M l = ρ E . Then in this case, from (18), we obtain the equation

d l-1 dt l-1 M -1 l J l              ȳ (t) . . . t t 0 • • • τ 2 t 0 ȳ (τ 1 ) dτ 1 • • • dτ l-1              = z 1 , ( 20
)
where

M l ∈ R ρ E ×ρ E and J l ∈ R ρ E × pl . Let U be a matrix so that rank D F U = rank D F =: m, U ∈ R q× m (21)
Since ( 9) must be satisfied according to Theorem 1, we have that

z 2 (t) = I n-ρ E 0 q D F I 0 0 U + ż1 (t) ȳ -Ā C z 1 (t)
(22) where q := n -ρ E + m. Now, we are ready to give a formula to reconstruct x in finite time.

Theorem 3. If system Σ is SO, then the state x can be expressed algebraically by the following formula: l+1) and H 2 ∈ R q× p(l+1) are defined as l+1) , and M l and J defined recursively in ( 11) and ( 19), respectively, and U defined by (21).

x (t) = d l dt l S 0 n× m H 1 H 2              ȳ (t) . . . t t 0 • • • τ 1 t 0 ȳ (τ 1 ) dτ 1 • • • dτ l              , (23) where H 1 ∈ R ρ E × p(
H 1 := 0 ρ E × p M -1 l J l , H 2 := DU FU + (G 1 -G 2 ) , G 1 := M -1 l J l 0 ρ E × p 0 p× pl I p× p , G 2 := 0 ρ E × p ĀM -1 l J l 0 p× p CM -1 l J l , G 1 , G 2 ∈ R ρ E + p× p(
Proof. Let us define the extended vector Y k ∈ R p(k+1) as follows

Y k =                     ȳ (t) t t 0 ȳ (τ 1 ) dτ 1 . . . t t 0 • • • τ 1 t 0 ȳ (τ 1 ) dτ 1 • • • dτ k                     , k = 1, 2, . . . (24)
Then for a matrix X of suitable dimension the following holds:

d k dt k XY k = d k+1 dt k+1 0 X Y k+1 .
Thus, since x = S z, by manipulating ( 22) and tanking into account ( 20), ( 23) is obtained.

Remark 1. One might obtain a little more from the previous analysis, that is, one can express by an algebraic formula the part of µ that can be reconstructed (assuming Σ is SO). Indeed, let μ be implicitly defined by the equation D T F T T µ = D T F T T U μ. With the same procedure followed to obtain (23), μ can be expressed by the formula

μ (t) = d l dt l 0 m×n I m H 1 H 2              ȳ (t) . . . t t 0 • • • τ 1 t 0 ȳ (τ 1 ) dτ 1 • • • dτ l              . (25)

Asymptotic Reconstruction

Let us assume that Σ is SD, but not SO. Next we show how to carry out the estimation of x (t).

Theorem 4. Assuming that Σ is SD, but not SO, we obtain that lim t→∞ x (t)x (t) = 0, provided that x (t) is designed by following ( 26)-( 27).

x (t) = d l dt l S 0 n× m H1 H2 Y l (t) S 0 n×m V K * ŵ2 (26) ẇ2 = Ā3 ŵ2 + Ā2 d l dt l 0 ρ E -ρ M × p J l + D2 U H2 Y l (27)
where H1 ∈ R ρ E × p(l+1) , H2 ∈ R q+ p× p(l+1) , and Ḡ1 , Ḡ2 ∈ R ρ M + p× p(l+1) (ρ M = ρ M l ) satisfy the following identities,

H1 = 0 ρ E × p M + l J l , H2 := D1 U FU + Ḡ1 -Ḡ2 , Ḡ1 := J l 0 ρ M × p 0 p× pl I p , G 2 := 0 ρ M × p Ā1 J l 0 p× p C1 J l .
Proof. In this case, as Σ is not SO, by differentiation, we are able to reconstruct w 1 = M l z 1 only, where ρ M l < ρ E . Since Φ is SD and by ( 21), we have that M l D F U has full column rank. Then, from (18) and from (13a) and (13c), we have the following expression for w 1 and z 2 ,

w 1 (t) = d l-1 dt l-1 J l Y l-1 (t) , (28)
z 2 (t) = d l dt l I n-ρ E 0 n-ρ E × m H2 Y l (t)+ I n-ρ E 0 n-ρ E ×m K * w 2 . ( 29
)
Thus, since x = S z and z 1 = M + l w 1 +Vw 2 , we obtain the identity

x (t) = d l dt l S 0 n× m H1 H2 Y l (t) + S 0 n×m V K * w 2 .
(30) Therefore, by comparing ( 26) with (30), we have that

x (t) = x (t) -S 0 n× m V K * (w (t) -ŵ2 (t)) . (31)
Furthermore, in view of (13b), (28), and ( 27), we obtain that ẇ2 -ẇ2 = Ā3 (w 2 -ŵ2 ) .

Since Φ is SD, Ā3 is Hurwitz; hence ŵ2 converges to w 2 . Then by (31) the proof is finished.

Remark 2. If μ needs to be reconstructed also, then it can be done by means of μ (t), defined as follows,

μ (t) = d l dt l 0 m×n I m H1 H2 Y l (t) 0 m×n I m V K * ŵ (t) (32
) where K * is implicitly defined by the equation

D F K * = D F I n-ρ E 0 0 U K *
So, we obtain straightforwardly that μ (t) -µ (t) goes to zero.

4.3. Summarized algorithm for the reconstruction of x (t) Below, there is a step-by-step description that may be followed to carry out the estimation of x (t).

Step 1. Choose S and T to bring E into the form (4).

Step 2. Calculate matrices Ā, C, D, and F according to (8).

Step 3. Calculate matrices M l and J l following the recursive algorithms ( 11) and ( 19), respectively. Let us remind that l is such that ρ M l+1 = ρ M l > ρ M l-1 .

Step 4A. Check if Σ is SO: test both conditions i) ρ M l = ρ E and ii) condition in (9). If they both are satisfied go to

Step 5.A.

Step 4B.

Check if Σ is SD: test the conditions i) rank D1 F = n -ρ E + rank D F
and ii) Ā3 is Hurwitz ( D1 and Ā3 defined in ( 14)). If they both are satisfied, go to Step 5.B. Otherwise, the exact estimation (even asymptotic) of x is impossible.

Step 5A. x (t) is reconstructed by means of the formula (23).

Step 5B. x (t) is estimated by using x (t) described in ( 26)-( 27).

To carry out Step 5 one needs to use a real-time differentiator. In the example presented in the next section, we use algebraic numerical and sliding mode differentiators, which have been proposed in the last years and have a considerable acceptation in several observation and identification procedures thanks to their finite time convergence and robustness. A brief review of those differentiators is given in Appendix A.

Example

Let us consider that Σ has the following matrices values

E =               0 1 0 1 0 0 0 0 1 0 1 1 0 0 0 0               , A =               -1 -1 1 1 2 1 2 0 1 1 1 -1 0 0 0 0               , D =               1 -1 0 0               C = 0 1 0 1 0 -1 0 2 , F = 0 1 .
It is easy to see that, in this example, det(sE -A) = 0 for every s ∈ C. Hence, many solutions for x (t) are expected to satisfy the differential equation in (1). However, to each output y (t) corresponds only one trajectory of x (t) (a.e.). Indeed, we will see that, according to Theorem 1, Σ is SO. We describe the observer design following the step by step description.

Step 1. For this case matrices S and T are chosen as follows,

S =               0 0 -1 -1 0 0 0 -1 0 1 1 0 1 0 0 1               , T =               1 0 0 0 -1 0 1 0 0 1 0 0 0 0 0 1               ,
Step 2. Matrices Ā, C, D, and F take the following values:

Ā = 1 1 -2 0 , D = 2 3 1 -2 -6 -1 , C =               0 2 0 0 1 0 2 0               , F =               0 -3 -1 0 0 0 0 0 0 0 3 1               .
Step 3. Matrices M 2 and J 2 take the values

M 2 = 1 0 √ 2 √ 2 , J 2 = 0 0 0 0 0 0 1 0 0 0 0 0 √ 2/2 0 0 √ 2/2 .
Step 4. Here, rank M 2 = 2, rank B = 3, and n -ρ E = 2. Therefore, both conditions of Theorem 1 are satisfied.

Step 5. The reconstruction of x (t) can be done by means of the formula (23), once we define ȳT = 0 0 y T . For this example, the reconstruction of µ is also possible following formula (25). Thus, we have

x 1 = - 1 3 y 1 + 7 12 y 2 - 1 2 ẏ1 + 1 6 ẏ2 , x 2 = 2 3 y 1 - 1 6 y 2 + 1 6 ẏ2 , x 3 = - 1 4 y 2 + 1 2 ẏ1 , x 4 = 1 3 y 1 - 1 6 y 2 + 1 6 ẏ2 , µ = 1 2 y 2 + 1 2 ẏ2 .
To estimate the state x and µ, we use two different differentiators, an algebraic numerical differentiator (ALND) and a high order sliding mode differentiator (HOSMD). For simulation purposes, we have chosen µ = 2 sin (x 1x 2 + x 3 ) + cos (t) and x 1 = x 2 -2cos (3t). The state trajectories are depicted in Figure 1. The estimation error is depicted in Figures 2 and3. The corresponding estimation for µ, denoted by μ, is shown in Figure 4. We have also tested both differentiators in presence of output noise, for that we have used an approximation of white noise using the Matlab command rand with a uniform distribution on the interval (-0.1, 0.1) for the first output and (-0.05, 0.05) for the second one. The estimation error is depicted in Figure 5.

Conclusions

We have given necessary and sufficient conditions to estimate the slow (non-impulsive) trajectories, for singular systems in which more than one solution of the differential equation is allowed, i.e. the pencil of the system is not required to be regular as it is assumed in most of the previous works where the observability is studied. We have given explicit formulas to reconstruct in finite time and asymptotically the states. Nevertheless, we have to notice that when the estimate of x needs an excessive number of derivatives, the error due to noise and sampling time could increase considerably. In that case, if an asymptotic estimation is enough, we should use the on-line differentiation only the needed times allowing after to design a Luenberger-like observer. In that case, a simple but cumbersome modified procedure might be followed in order to at minimum the number of derivatives required to estimate the state, we refer the interested reader to [START_REF] Floquet | On sliding mode observers for systems with unknown inputs[END_REF] and [START_REF] Bejarano | High order sliding mode observer for linear systems with unbounded unknown inputs[END_REF].

Figure 2 :

 2 Figure 1: Trajectories of x (t).

Figure 3 :

 3 Figure 3: Estimation error e = xx using HOSMD.

Figure 4 :Figure 5 :

 45 Figure4: UI µ (solid) and its estimation μ (dashed for ALND, dotted HOSMD).

Such a statement was proven in[START_REF] Hautus | Strong detectability and observers[END_REF].

Acknowledgement

This research has been financially supported by EU Interreg IV A 2 Mers Seas Zeeën Cross-border Cooperation Programme under SYSIASS project 06-020. For more details about the SYSIASS project see the project website www.sysiass.eu. F.J. Bejarano held a posdoctoral position at INRIA -Nord Europe with the ALIEN team, from October 2010 to September 2011, period during which this work was produced.

Appendix A. Review of two differentiators

Appendix A.1. ALND This algebraic setting for numerical differentiation of noisy signals is introduced in [START_REF] Fliess | Compression différentielle de transitoires bruités[END_REF] and analyzed in Liu et al. (2011a,b); [START_REF] Mboup | A revised look at numerical differentiation with an application to nonlinear feedback control[END_REF]Mboup et al. ( , 2009b) ) (see also [START_REF] Nöthen | Beiträge zur rekonstruktion nicht direkt gemessener größen bei der silizium-einkristallzüchtung nach dem czochralski-verfahren[END_REF] for interesting discussions and comparisons). The reader may find additional theoretical foundations in [START_REF] Fliess | Analyse non standard du bruit[END_REF]; [START_REF] Fliess | An algebraic framework for linear identification[END_REF]. Consider a signal y(t) = ∞ l=0 y (l) (0) t l l! which is assumed to be analytic around t = 0 and its truncated Taylor expansion y N (t) = N l=0 y (l) (0) t l l! at order N. The usual rules of symbolic calculus in Schwartz's distribution theory [START_REF] Schwartz | Théorie des distributions[END_REF]) yield y (N+1) N (t) = y(0)δ (N) + . . . + y (N) (0)δ, where δ is the Dirac measure at zero. Multiply both sides by (-t) l : (-t) l y (N+1) N (t) = (-t) l y(0)δ (N) + . . . + y (N) (0)δ , and apply the rules tδ = 0, tδ (l) = -lδ (l-1) , l ≥ 1. We obtain a triangular system of linear equations from which the derivatives

It means that the coefficients y(0), . . . , y (N) (0) are linearly identifiable [START_REF] Fliess | An algebraic framework for linear identification[END_REF]; [START_REF] Fliess | Closed-loop parametric identification for continuous-time linear systems via new algebraic techniques[END_REF]. The time derivatives of y N (t), the Dirac measures and its derivatives are removed by integrating with respect to time both sides of Eq. (A.1) at least ν times (ν > N):

The iterated integrals may be replaced by

2) It is clear that the numerical estimation rely on lim N→+∞ [y (l) N (0)] estim (t) = y (l) (0).

Remark 3. These iterated integrals are low pass filters which attenuate the noises, which are viewed as highly fluctuating phenomena (see [START_REF] Fliess | Analyse non standard du bruit[END_REF] for more details). The above formulae may easily be extended to sliding time windows in order to obtain real time estimates (see [START_REF] Mboup | A revised look at numerical differentiation with an application to nonlinear feedback control[END_REF]Mboup et al. ( , 2009b)); Liu et al. (2011a,b) for further details). Moreover, according to the performed algebraic manipulation one can have some different formulae:

3) where x is assumed to be analytic, P µ,k n is the n th order Jacobi polynomial defined on [0, 1] with n ∈ N, k, µ ∈] -1, ∞[, β = ±1 and T > 0.

Appendix A.2. HOSMD

The HOSMD described in [START_REF] Levant | High-order sliding modes: differentiation and outputfeedback control[END_REF] can be expressed in a dynamic form by means of the equations (A.4), where the signal to be differentiated is represented by a function

Then, in the cited paper it was shown that if λ i s gains are chosen properly, the differentiator converges in a finite time T , i.e., α i (t) = f (i) (t), for all t ≥ T and i = 0, 1, . . . , r. L is a constant such that f (r+1) (t) ≤ L. For the case when r = 5, the gains could be chosen as λ 0 = 1.1, λ 1 = 1.5, λ 2 = 3, λ 3 = 5, λ 4 = 8, λ 5 = 12.