Nathalie Verdière

Lilianne Denis-Vidal

Ghislaine Joly-Blanchard

A new method for estimating derivatives based on a distribution approach

Keywords: Derivatives estimation, Distributions, Test Functions

In many applications, the estimation of derivatives has to be done from noisy measured signal. In this paper, an original method based on a distribution approach is presented. Its interest is to report the derivatives on infinitely differentiable functions. Thus, the estimation of the derivatives is done only from the signal. Besides, this method gives some explicit formulae leading to fast calculus. For all these reasons, it is an efficient method in the case of noisy signals as it will be confirmed in several examples.

Introduction

In many estimation and observation problems [START_REF] Diop | Interpolation and numerical differentiation for observer design[END_REF], [START_REF] Gauthier | A simple observer for nonlinear systems: Application to bioreactors[END_REF], [START_REF] Verdière | Identifiability and estimation of pharmacokinetic parameters of ligands of macrophage mannose receptor[END_REF], derivatives, sometimes of high order, have to be evaluated from a noisy signal given by discrete samples of measurements. It is wellknown that it is an ill-posed problem in numerical analysis and in signal processing and control. Several approaches for estimating the derivatives of a noisy signal have been proposed in the literature. We can name, for example, methods which consist in approximating the signal by polynomials using least-squares and adding a regularization criteria [START_REF] Boor | A practical Guide to Splines[END_REF], [START_REF] Ibrir | A numerical procedure for filtering and efficient high-order signal differentiation[END_REF]. A recent one proposed in [START_REF] Ditkowski | Computing derivatives of noisy signals using orthogonal functions expansions[END_REF] consists in projecting the signal onto a set of orthogonal functions or modes and approximating the projection by a linear combination of a few of these modes. Fliess et al. have introduced a method based on the local approximation of the signal by a truncature Taylor expansion and the properties of the Mikusinski field [START_REF] Fliess | Questioning some paradigms of signal processing via concret examples[END_REF], [START_REF] Mboup | Numerical differentiation with annihilator in noisy environment[END_REF]. A study of errors due to the method and the noise on the signal has been treated in [START_REF] Liu | Error analysis of a new class of derivative estimators for noisy signals[END_REF]. A new method using the distribution theory and the test functions is proposed in this paper. It is based too on the idea of approximating locally the signal by a truncature Taylor expansion but then it uses the properties of the distribution theory. Its first advantage comes from the fact that it allows one to report the derivatives on the test functions which are explicitly known. Thus, the estimation of the derivatives is done only from the signal. Then, this theory introduces naturally integrals which enable us to annihilate a part of the noise on the signal in our estimations. In [START_REF] Mboup | Numerical differentiation with annihilator in noisy environment[END_REF], their formulae are based on the estimation of integrals too which are obtained by integral iterations. An other advantage of our method is that it does not require any statistic knowledge about the noise and gives some explicit expressions leading to fast calculus. Compared to the one proposed by [START_REF] Mboup | Numerical differentiation with annihilator in noisy environment[END_REF], no delay has to be taken into account for improving the estimation of the derivatives. Thus, the derivatives of the signal can be estimated from the transmission of the signal.

This paper is organized as follows. In section 2, the derivative estimation method based on the distribution theory is presented. From some theoretical recall done in section 2.1, the derivative estimation method is deduced and presented at section 2.2 as well as the explicit formulae. In section 3, two types of errors are analyzed: the first one is the error due to the method and which does not consider the noise, the second one is the error due to the noise on the signal. Two inequalities will be given and will serve in the numerical applications, first in giving a criteria for choosing the best test function among others, then for reducing the search of test functions support length. In section 4, the method is applied on different examples. In order to see how the method works, a simple model based on a sinusoid is used. Then, a pharmacokinetic model is taken again from [START_REF] Verdière | Identifiability and estimation of pharmacokinetic parameters of ligands of macrophage mannose receptor[END_REF]. In this paper, the estimation of the derivatives specially at 0 + , the right limit of 0, was important for a parameters estimation problem. Thus, the method is particularly well adapted since it gives explicit formulae at 0 + . Finally, in 4.4, our method is compared with two others. The first one is the method proposed by S. Ibrir and S. Diop in [START_REF] Ibrir | A numerical procedure for filtering and efficient high-order signal differentiation[END_REF]. In this article, they proposed a regularization method based on the B-splines. The second one comes from the article [START_REF] Mboup | Numerical differentiation with annihilator in noisy environment[END_REF] where Fliess et al. analyzed and improved their method by considering barycentric coordinates and by authorizing a delay.

A derivative estimation method based on the distribution theory

Theoretical results

The signal x is supposed to be an analytic function on ]0, a[, equal to zero on ] -b, 0] where a > 0 or a = +∞ and b > 0 or -b = -∞. Besides, x is supposed to coincide on ]0, a[ with an analytic function x defined on ] -b, a[. Afterwards, we will say that x verifies the assumption H.

The signal x can also be associated to a distribution which, in this case, is regular and which operates on the test functions space noticed D(] -b, a[). These functions are infinitely differentiable and their supports are a compact in ] -b, a[. The associated distribution is noticed T x and is defined by:

< T x , ϕ >= ]-b,a[ x(t)ϕ(t)dt where ϕ ∈ D(] -b, a[).
One of the motivations of distributions' use is that they are infinitely differentiable and that the derivation can be reported on the test function. Besides they can take into account the discontinuities, in our case at 0 as it is explicited in the following theorem [START_REF] Rudin | Functional Analysis[END_REF]:

Theorem 1. The Dirac distribution δ 0 verifies for all integer i: < δ (i) 0 , ϕ >= (-1) i ϕ (i) (0).

Moreover, if x verifies the assumption H, one gets for all

k ∈ [0, N ]: < T (k) x , ϕ >= (-1) k < T x , ϕ (k) >, (1)
T (k) x = T x (k) + k-1 i=0 σ k-1-i δ (i) 0 , k > 0 ( 2
)
with σ j = x (j) (0 + ) -x (j) (0 -) = x (j) (0 + ) for 0 ≤ j ≤ N .
The notation 0 + (respectively 0 -) designates the right limit at 0 (respectively left limit) of the corresponding functions.

The use of the truncated Taylor expansion for the derivatives estimation is classical (finite difference, etc...). The difference is the way to use it. Let us consider x a function verifying H and t l ∈]0, a[ or t l = 0 + for the right limit at 0. Let D(t l , R) defined for all R by {t ∈] -b, a[; | t -t l |< R} Then x is a real analytic function in ] -b, a[ and it can be expanded in a power series about the point t l which is its Taylor expansion at t l . There exists a number h, h > 0, such that the series converges in D(t l , h) and converges uniformly in the compact set D(t l , ǫ) for all ǫ such that 0 < ǫ < h. Let ǫ be given such that 0 < ǫ < h.

If 0 ∈]t l -ǫ, t l + ǫ[ (respectively 0 ∈ [t l -ǫ, t l + ǫ]),
the function x is defined by:

x(t) =          +∞ n=0 x (n) (t l ) (t -t l ) n n! , t ∈]0, t l + ǫ] (resp. t ∈ [t l -ǫ, t l + ǫ]) 0, t ≤ 0, (3)
and the function x N which corresponds to the truncated Taylor expansion

x N (t) =          N n=0 x (n) (t l ) (t -t l ) n n! , t ∈]0, t l + ǫ] (resp. t ∈ [t l -ǫ, t l + ǫ]) 0, t ≤ 0, (4)
is an approximation of x(t). Now consider some test functions ϕ ∈ D(] -b, a[) whose supports are in [t l -ǫ, t l + ǫ]. Afterwards, ϕ is supposed to be positive without any restriction. According to theorem 1, three cases appear. Before giving them, let us

introduce α i = 1 i! Vt l (t -t l ) i ϕ(t)dt
, where V t l is specified in each cases. The function test are chosen such that α 0 = 0. According to (1), ( 2) and ( 4), one gets the three following cases:

-If 0 ∈ [t l -ǫ, t l + ǫ], V t l = [t l -ǫ, t l + ǫ] and for k = 1, . . . , N , N i=k x (i) (t l )α i-k = (-1) k t l +ǫ t l -ǫ x N (t).ϕ (k) (t)dt. (5)
-If

t l = 0 + , V t l = [0, ǫ] and for k = 1, . . . , N , N i=k x (i) (0 + )α i-k + k-1 i=0 (-1) k-1-i x (i) (0 + )ϕ (k-1-i) (0) = (-1) k ǫ l x N (t).ϕ (k) (t)dt, (6)
-If 0 ∈]t l -ǫ, t l + ǫ[, V t l = [0, t l + ǫ] and for k = 1, . . . , N ,

k-1 i=0 (-1) k-1-i x (i) N (0 + )ϕ (k-1-i) (0) + N i=k x (i) (t l )α i-k = (-1) k t l +ǫ l x N (t)ϕ (k) (t)dt (7)
And in each previous case, one gets:

N i=0 x (i) (t l )α i = Vt l x N (t).ϕ(t)dt, (8)
In systems ( 5), ( 6), ( 7) and ( 8), x N (t) which is an approximation of x(t) is generally an unkown value. In concrete applications, only measured values of x(t) corrupted by noise are available. Afterwards, y(t) = x(t)+ η(t) where η will denote the noise is the measured signal. The aim of the following section will be to give systems depending uniquely of the measured value y(t) and allowing the estimation of the derivatives.

Estimation of the derivatives

In this section, formulae deduced from ( 5), ( 6), ( 7) and ( 8) for evaluating derivatives are explicited in substituting x N (t) by y(t). The deduced systems provide us an estimator C (i) ǫ,N (t l ) of x (i) (t l ) when 0 ∈ [t l -ǫ, t l + ǫ] (respectively t l = 0 + and 0 ∈]t l -ǫ, t l + ǫ[). As previously, the test function

ϕ ∈ D(] -b, a[) has a support in [t l -ǫ, t l + ǫ]. Consider, A ǫ,ϕ,N =      α 0 α 1 . . . α N 0 α 0 . . . α N -1 . . . . . . . . . . . . 0 . . . . . . α 0      , (9)
B ǫ,ϕ,N =        0 0 . . . . . . 0 ϕ(0) 0 . . . . . . 0 -ϕ ′ (0) ϕ(0) 0 . . . 0 . . . . . . . . . . . . . . . (-1) N-1 ϕ (N-1) (0) . . . . . . ϕ(0) 0        , (10)
C ǫ,N (t l ) = C (i) ǫ,N (t l ) 0≤i≤N , X(t l ) = (x (i) (t l )) 0≤i≤N , (11)
and finally,

I ǫ,ϕ,x (t l ) = (-1) i Vt l x(t)ϕ (i) (t)dt 0≤i≤N . (12)
Clearly, matrices A ǫ,ϕ,N and B ǫ,ϕ,N are Triangular Toeplitz matrix.

System (5) completed by [START_REF] Liu | Error analysis of a new class of derivative estimators for noisy signals[END_REF] with V t l = [t l -ǫ, t l +ǫ] and ϕ 1 the corresponding test function, becomes

A ǫ,ϕ1,N X(t l ) = I ǫ,ϕ1,xN (t l ) (13)
and the estimate C ǫ,N (t l ) verifies:

A ǫ,ϕ1,N C ǫ,N (t l ) = I ǫ,ϕ1,y (t l ). ( 14
)
Matrix A ǫ,ϕ1,N is triangular and det(A ǫ,ϕ1,N ) = α N 0 = 0. System (6) completed by [START_REF] Liu | Error analysis of a new class of derivative estimators for noisy signals[END_REF] with V t l = [0, ǫ] and ϕ 2 the corresponding test function becomes

(A ǫ,ϕ2,N + B ǫ,ϕ2,N )X(0 + ) = I ǫ,ϕ2,xN (0 + ) ( 15
)
and the estimate C ǫ,N (0 + ) verifies

(A ǫ,ϕ2,N + B ǫ,ϕ2,N )C ǫ,N (0 + ) = I ǫ,ϕ2,y (0 + ). ( 16
) Matrix D ǫ,ϕ2,N = A ǫ,ϕ2,N + B ǫ,ϕ2
,N is a Toeplitz matrix and has a non-zero determinant with a good choice of test function ϕ 2 .

System [START_REF] Ibrir | A numerical procedure for filtering and efficient high-order signal differentiation[END_REF] completed by ( 8) with V t l = [0, t l + ǫ] and ϕ 3 the corresponding test function becomes

A ǫ,ϕ3,N X(t l ) = I ǫ,ϕ3,xN (t l ) -B ǫ,ϕ3,N X(0 + ) (17)
and the estimate C ǫ,N (t l ) verifies:

A ǫ,ϕ3,N C ǫ,N (t l ) = I ǫ,ϕ3,y (t l ) -B ǫ,ϕ3,N C ǫ,N (0 + ) ( 18
)
where C ǫ,N (0 + ) is given by ( 16). Matrix A ǫ,ϕ3,N is triangular and

det(A ǫ,ϕ3,N ) = α N 0 = 0.
In all the cases, the estimation

C (i)
ǫ,N (t l ) of the derivatives x (i) (t l ) will be done only from the observation y.

Error analysis

Two types of error have to be analyzed. The first one which does not consider the noise is the error due to the method. It permits to know if the problem is well-posed. The second one is the error due to the noise on the signal. A compromise has to be down between these two errors. Indeed, in order to minimize the first one, the test functions support length must be as small as possible. However, it is in opposition for reducing the second one since the noise is better filtered if the integration interval is greater. Afterwards, the errors will be evaluated in the more usual case that is, the first one when the support of the test function, denoted ϕ 1 , is included in [t l -ǫ, t l +ǫ] and [t l -ǫ, t l + ǫ] ⊂]0, a[. Thus, only the matrix A ǫ,ϕ1,N is considered. From the studie of these errors, a criteria will be given for choosing the best test function among others. Then a method for reducing the search of the test functions support length will be deduced. Finally, since our estimates are obtained as output of FIR filters with a noisy input, the output signal to noise ratio can be estimated in the case of a zero-mean white Gaussian and will be given at the end of this section.

Error due to the method

In this section, the notations given in sections 2.1 and 2.2 are taken again and the test functions are supposed to be positive. In order to ease the presentation, the proofs of this section are given at the end of this article.

According to the following theorem, the error due to the method can be controlled owing to the choice of the test function and its support length.

Theorem 2. For all N ∈ IN by using previous notations, there exist real numbers S, R and r such that for ǫ ≤ R:

C ǫ,N (t l ) -X(t l ) 1 ≤ S ǫ r N +1 2ǫ φ(N ) (19) with φ(N ) = N i=0 φ i and φ i = sup t∈[t l -ǫ,t l +ǫ] (| ϕ (i) 1 (t) |).
φ(N ) gives a criteria for choosing a test function from a family of them, a similar criteria will be given in section 3.2.

The following proposition gives the way that the condition number evolves according to the support length of the test function.

Proposition 1. For all N ∈ IN, 1 ≤ χ(A ǫ,ϕ1,N ) ≤ e ǫ (1 + h(ǫ)) (20) with lim ǫ→0 h(ǫ) = 0
The latter proposition ensures that for an interval length sufficiently small, the condition number of the matrix defined by system ( 14) is near to one and this for any test functions. Thus, the resolution of this system is a stable problem according to the perturbations of the second member. However, since in the numerical applications the integrals will be estimated by a discrete method and the noise has to be considered, the support length can not be taken as small as possible.

Error due to the noise on the signal

In [START_REF] Fliess | Questioning some paradigms of signal processing via concret examples[END_REF], the authors bring out two kinds of perturbations. The first ones, like the constant perturbations whose amplitude is unknown are said structured and are solutions of a linear homogeneous differential equation. They can be annihilated by linear differential operators and it is, for example, the case of a random noise whose mean is constant and unknown. The second ones are said unstructured and correspond to a high frequency perturbation which can be reduced by low pass filters like iterative integrals. Thus, any noise η(t) can be written as η(t) = η 0 (t) + γ where η 0 (t) is a noise whose mean is equal to zero and γ is a constant representing the mean of η(t).

The method based on the distribution theory permits to annihilate a part of the structured noises. When the estimation of the derivatives are done with the formulae ( 14), ( 16), (18), the measured signal is only in the second member of the equations. According to the integral linearity, noise appears as

t l +ǫ t l -ǫ η(t)ϕ (k) (t). Since η(t) = η 0 (t) + γ, one gets t l +ǫ t l -ǫ η(t)ϕ (k) (t) = t l +ǫ t l -ǫ η 0 (t)ϕ (k) (t)dt + γ t l +ǫ t l -ǫ ϕ (k) (t)dt. If k ≥ 1, t l +ǫ t l -ǫ ϕ (k) (t) = ϕ (k-1) (t l + ǫ) -ϕ (k-1) (t l -ǫ) is equal to zero as soon as 0 + / ∈ [t l -ǫ, t l + ǫ]
which will be generally our case in the numerical applications. Thus, constant perturbations are annihilated.

On the other hand, if 0 + ∈ [t l -ǫ, t l + ǫ] then

t l +ǫ t l -ǫ ϕ (k) (t) = -ϕ (k-1) (0 + ).
Thus, constant perturbations can be controlled owing to the test function. Furthermore, according to the following inequality, the choice of the test function can permit to attenuate the zero-mean noise effect. Indeed, one gets:

t l +ǫ t l -ǫ η 0 (t)ϕ (i) (t)dt 0≤i≤N 2 1 ≤ t l +ǫ t l -ǫ (η 0 (t)) 2 dt × ψ ǫ,ϕ,N (21)
and

E( t l +ǫ t l -ǫ (η 0 (t)) 2 dt × ψ ǫ,ϕ,N ) ≤ t l +ǫ t l -ǫ var(η 0 (t))dt × ψ ǫ,ϕ,N (22)
where var denotes the variance, E(X) the expected value of a random variable

X and ψ ǫ,ϕ,N = N i=0 t l +ǫ t l -ǫ |ϕ (i) (t)| 2 dt. ψ ǫ,N,ϕ
gives a criteria for choosing a test function among others. Besides, it permits to determine a lower bound for the test functions support length in order to attenuate a part of the noise. This criteria will be tested on the first example of section 4.2.

Let us return to the first case that is when ϕ 1 is considered. In the case of a zero-mean white Gaussian, the bias and variance can be evaluated in order to validate the results. Indeed, if y(t) = x(t) + η 0 (t) where x(t) is the non noisy signal and η 0 is a zero-mean white gaussian noise then the signal derivatives estimates are obtained as output of FIR filters with a noisy input. Suppose that L is the number of samples in [a, b] regularly spaced 1 and the test function used to estimate the derivatives of x at the sampling time t l (0 ≤ l ≤ L) has its support contained in an interval of the form [t l-M , t l+M ]. Each estimated derivative can be expressed as:

C (n) ǫ,N (t l ) = tM -tM β(s)y(t l -s) ds,
where β is a linear combination of the test function ϕ 1 (t l -.) and its derivatives. If W m are the weights associated with a given numerical integration method and β m = β(t m ), for each sampling time t l , the following expression is computed:

C (n) ǫ,N (t l ) ≈ C(n) ǫ,N (l) = M m=-M W m β m y l-m .
1 Notice that the following reasoning can be done with a non-uniform sampling.

Denote c m = W m β m and x

(n) l = x (n) (t l ). Thus, the output signal to noise ratio (SNRout) is given by:

SN Rout = 1 L-2M L-M l=M |x (n) l | 2 σ 2 M m=-M |c m | 2 .
Concerning the bias, the relative sample mean squared error due to the bias term and to the numerical approximation of the integrals is computed as:

SBR = L-M l=M |x (n) l | 2 L-M l=M |x (n) l -E s (C (n) ǫ,N (l))| 2 , where E s (C (n) ǫ,N (l)) is the sample means.
In the numerical applications, the SNRout and SBR will be given in the decadic logarithm.

Reduction of the support length search

The estimation of a upper bound of the relative error between C ǫ,N (t l ) and X(t l ) will allow one to determine the maximal value of the test functions support length. Without noise, one gets:

C ǫ,N (t l ) -X(t l ) 2 2 X(t l ) 2 2 ≤ χ(A ǫ,ϕ1,N ) 2 I ǫ,ϕ1,x (t l ) -I ǫ,ϕ1,xN (t l ) 2 2 I ǫ,ϕ1,x (t l ) 2 2 (23)
where χ(A ǫ,ϕ1,N ) is the condition number of the matrix A ǫ,ϕ1,N . Thus,

X(t l ) 2 2 C ǫ,N (t l ) -X(t l ) 2 2 ≥ Ẽǫ,ϕ1,N M 2 χ(A ǫ,ϕ1,N ) 2 (24)
where M = sup

t∈[t l -ǫ,t l +ǫ] |x (N +1) (t)| and Ẽǫ,ϕ,N = N i=0 t l +ǫ t l -ǫ x(s)ϕ (i) (s)ds 2 N i=0 t l +ǫ t l -ǫ (s -t l ) N +1 (N + 1)! ϕ (i) (s)ds 2 (25)
In order to obtain good numerical results, Ẽǫ,ϕ,N should be maximum. Ẽǫ,ϕ,N can be put in correspondance with the SBR. Thus, in the same way as previously and since in the numerical applications only the measured signal y(t) will be accessible at discrete samples, E ǫ,ϕ,N defined by

E ǫ,ϕ,N = 10 log 10 N i=0 l+M k=l-M y(t k )ϕ (i) (t k ) 2 N i=0 l+M k=l-M (t k -t l ) N +1 (N +1)! ϕ (i) (t k ) 2 (26)
will be used for reducing the search of the support length of the test functions.

Numerical applications

In the previous section, the method for estimating the derivatives from the distribution theory have been established as well as two inequalities. The first one depends on ψ ǫ,ϕ,N (22) and minimizes the error due to the noise. The second one depends on E ǫ,ϕ,N and has to be maximized. Since ψ ǫ,ϕ,N does not depend on the signal, it can be used as a criteria for choosing the best test function and as it will be seen in our examples, it can give the minimal value of the test functions support length. Thus, ψ ǫ,ϕ,N will be analyzed with two test functions in a first subsection. Then, a simple example is treated for testing the formulae presented at section 2. In this example, (26) will be analyzed and a method for reducing the search of the interval length will be deduced from ψ ǫ,ϕ,N and E ǫ,ϕ,N . An application in the pharmacokinetic domain will be then treated. And finally, two examples will be used in order to compare this derivative method with two classical ones. Thereafter, all the integrals are evaluated with the trapezoidal method.

Study of ψ ǫ,ϕ,N

ψ ǫ,ϕ,N is an expression depending uniquely on the chosen test function, its derivatives and N and it minimizes the error due to the noise. Thus, in choosing correctly ǫ, ϕ, N for minimizing ψ ǫ,ϕ,N , this choice permits to attenuate the effect of zero-mean noise on the signal. In this paper, the measured data are supposed to come from an analytic function. When the derivatives are estimated at t l , bell-shaped functions centered at t l have been chosen in order to take into account all the informations of the measured signal in the neighborhood of t l . ψ ǫ,ϕ,N will be evaluated with the two following test functions:

φ(t) = e ǫ 2 t 2 -ǫ 2 if | t |< ǫ 0 if | t |≥ ǫ, (27) and φ
(t) = cos 2 (t)e ǫ 2 t 2 -ǫ 2 if | t |< ǫ 0 if | t |≥ ǫ, (28)
where [-ǫ, ǫ] is the support of the test function. A translation allows to obtain an infinitely differentiable function on any interval. Afterwards, the solid line curve corresponds to φ and the dashed line curve to φ. Figure 1 (resp. 2) represents Ψ ǫ,ϕ,1 (resp. Ψ ǫ,ϕ,2 ). First, one can notice that for attenuating the effect of the noise with the distribution derivative method, the support length of the chosen test functions must be greater than one. Furthermore, the test function φ seems more efficient than φ for attenuating the noise on the signal. These remarks will be confirmed in the following example.

Example 1: Study of E ǫ,ϕ,N

In this subsection, a very simple model is chosen:

y(t) = sin(̟ 0 + ̟ 1 t) + n(t), t ∈ [0, 10]
where ̟ 1 = 0.8 is the angular frequency of the signal, ̟ 0 = 1 the phase and n(t) is a perturbation. The latter is supposed to follow the normal law with mean-zero and the measures are supposed to be done at each step h = 0.05, thus there are 200 measured points on the interval [0, 10]. The signal y is drawn at figure 3. For estimating the first and second order derivatives, the systems The derivative estimation method presented at section 2 will be tested with the two infinitely differentiable test functions φ (27) and φ (28). E ǫ,ϕ,N given by (26) provides in our examples a criteria for estimating the maximal value of the test functions support length. Figure (4) (resp. ( 5)) represents E ǫ,ϕ,1 (resp. E ǫ,ϕ,2 ). According to the figures, results on the derivative estimates are expected to degrade rapidly as soon as the support length is greater than 3. Furthermore, in this example, for having equivalent results with φ and φ, the support length of φ must be smaller than φ when it is greater than 1.5. In return, the effect of the noise will be less attenuate. In taken again the remark done at section 4.1, that is the test functions support length must be greater than one, the support length is chosen between 1 and 3 according to the value of N . Indeed, when seeking to estimate a highter order derivative, the support length must be greater for better filter the noise.

For N = 1, the interval length of φ and φ is equal to 1.5 and their curves coincide. Thus, only the first derivative evaluated with φ is represented at For N = 2, the interval length for the two test functions is equal to 3. According to figure 5, the results for φ should not vary much when the support length vary near 3 and this is confirmed by numerical tests. The second derivative evaluated with φ is represented at figure 7. Some irregularities can be noticed for the second derivative unlike the one evaluated with φ (figure 8).

For avoiding this phenomenon, the support length of φ must be chosen greater.

The results on the SNRout and the SBR are given at table I One can notice that numerically, the best results are obtained when the order of the truncature Taylor expansion corresponds to the order of the derivative. Furthermore, when the order of the derivative increases, the length of the test function support has to be taken greater. Afterwards, the function test φ is chosen in the numerical estimation of the derivatives. Besides, the derivatives will be estimated with ( 14), ( 16), ( 18) and the truncature Taylor expansion will correspond to the desired derivative order.

Example 2

The following example ( 29) is now taken from [START_REF] Verdière | Identifiability and estimation of pharmacokinetic parameters of ligands of macrophage mannose receptor[END_REF]. Since in this article the estimation of the derivatives near 0 was important for finding the system parameters, the distribution approach is appropriate. This example has been chosen owing to the difficulties of estimating the derivatives, specially the second order derivative whose values are of order 10 -3 . Even if they seem negligible, they were important for estimating the system parameters.

     ẋ1 = k 12 (x 2 -x 1 ) -k v x 1 1 + x 1 , x 1 (0) = x 10 , ẋ2 = k 21 (x 1 -x 2 ), x 2 (0) = 0 y = x 1 (29)
where x ∈ IR 2 and y ∈ IR denote the state variables and the output respectively.

The simulated signal y represented at figure 9 is calculated from the true signal ȳ computed with θ = (0.011, 0.02, 0.1) and with the initial condition x 10 = 0.625. The output y is supposed to follow a random law with ȳ mean and (σ ȳ) 2 variance. The coefficient σ is computed so that the relative error has a maximum value of 0.05 with an error probability less than 0.003. The observations are supposed to be done at discrete time (t i ) i=1,...,N on the interval [0, 117] with a step equal to 0.5.

With the same principle as previously, one can deduce from E ǫ, φ,1 and E ǫ, φ,2 that the intervals length of φ must be less than 30. The best results are given for an interval length equal to 24 for the first derivative and 29 for the second derivative. The first and second derivatives curves are represented at figure 10 and 11 respectively. As it can be noticed, the results are very satisfactory.

Comparison with other methods

In this section, two other methods for estimating derivatives are considered and compared with the distribution method. The first one is presented in [START_REF] Ibrir | A numerical procedure for filtering and efficient high-order signal differentiation[END_REF] and is based on the B-splines. The second one is proposed by M. Fliess et al. ( [START_REF] Fliess | Questioning some paradigms of signal processing via concret examples[END_REF], [START_REF] Mboup | Numerical differentiation with annihilator in noisy environment[END_REF]). As for the distribution method, these two methods do not require any knowledge of the statistics of measurement uncertainties.

Example 3

This example taken from [START_REF] Ibrir | A numerical procedure for filtering and efficient high-order signal differentiation[END_REF] is based on the following noisy signal:

y(t) = cos(30t) sin(t) + ǫ(t)
and is supposed to be measured for each step δ = 0.01 s. The function ǫ follows the normal law whose the variance is unknown. The noisy signal is represented at figure [START_REF] Verdière | Identifiability and estimation of pharmacokinetic parameters of ligands of macrophage mannose receptor[END_REF]. The method proposed in [START_REF] Ibrir | A numerical procedure for filtering and efficient high-order signal differentiation[END_REF] is a regularization method based on the De Boor criteria, 1978, [START_REF] Boor | A practical Guide to Splines[END_REF] :

J = 1 n n i=1 (y i -ŷi ) 2 + λ t l ŷ(m) (s)ds, (30)
which contains a least-squares term symbolizing the connection between the solution of the optimization problem and the measures and a second smoothing term measuring the degree of the solution filtering. The equilibrium between the two distances is controlled by the choice of a parameter λ which will be obtained from the optimality condition of GCV criterion (Generalized Cross-Validation criterion) introduced by Craven and Wahba in 1979. De Boor has shown in 1978 that the best results were obtained with a spline function of order 2m. In [START_REF] Ibrir | A numerical procedure for filtering and efficient high-order signal differentiation[END_REF], the 3-B-splines have been chosen and an explicit algorithm is deduced.

In drawing E ǫ, φ,1 and E ǫ, φ,2 , the support length of φ must be less than 0.3. The interval on which the integrals are evaluated have a length equal to 0.09 (resp. 0.11) for the first derivative (resp. second derivative) which corresponds to 9 points of integration (resp. 11).

The results on the first and second derivatives are represented at figures 13 and 14 respectively. which is the third Jacobi orthogonal polynomial. The parameters κ and µ will be given afterwards.

Let y(t) = x(t) + ̟(t), 0 ≤ t ≤ 5 (without unit) be the noisy measurement of the signal x(t) = tanh(t -1) + e -t/1.2 sin(6t + π) with t ∈ [0, 5] and which is measured for each δ = 0.025. It corresponds to 1000 samples in the interval [0, 5]. The signal is supposed to be noisy by a zero-mean white gaussian iid sequence. The variance is adjusted in such a way the signal to noise ratio in dB, i.e, SN R = 10 log 10 |yi| 2 |̟i| 2 , corresponds to SN R = 25 dB (see figure 15).

In [START_REF] Mboup | Numerical differentiation with annihilator in noisy environment[END_REF], the authors test their method with different values of parameters. Moreover, in order to validate their results, they investigate the bias and variance in each simulation.

For N = 1, the integrals involved in the formulae are estimated with 60 For the second derivative, the integrals are calculated with 110 points for the two methods except in the neighborhood of 0 + for the distribution method. Indeed, they are evaluated with 80 points. However, the derivatives in the neighborhood of 0 + are not taken into account for calculating the SBR and the SNRout. The parameters for the Mboup et al. method are N = 2, κ = 0, µ = 1 and the SBR and SNRout are equal respectively to 18.1dB and 22.6dB. For the distribution method, we obtain SBR=16.2dB and SNRout=27.2dB. The obtained figures are drawn at figure 17. In the case of 70 points in the test functions support lenght, the results are: SBR=20.3 and SNRout=20.2.

Conclusion

In this paper, a new derivative method based on the distribution theory and the test functions is presented. Its first interest comes from the fact that it allows one to report the derivatives on the test functions. Thus, the estimation of the derivatives is done only from the signal. Then, some explicit expressions requiring any knowledge of the statistics of measurement uncertainty are deduced and they lead to fast calculus. Finally, the fact that a regular distribution is a linear integration functional and the use of test functions whose support is chosen adequately permit to annihilate a part of the structured noises. This method was applied and compared with two other methods in order to estimate the derivatives of noisy signals and as it was seen, it gives satisfactory results.

Proofs

Recall that the inverse of an upper triangular Toeplitz invertible matrix is also an upper triangular Toeplitz matrix.

Proof (Proof of theorem 2) From ( 13) and ( 14) the following inequality is obtained:

C ǫ,N (t l ) -X(t l ) 1 ≤ A -1 ǫ,ϕ1,N 1 I ǫ,ϕ,x (t l ) -I ǫ,ϕ,xN (t l ) 1 . (31)
The two following lemmae will give a priori bounds on A -1 ǫ,ϕ1,N 1 and I ǫ,ϕ,x (t l ) -I ǫ,ϕ,xN (t l ) 1 .

Lemma 1. There exist real numbers M and r such that:

I ǫ,ϕ,x (t l ) -I ǫ,ϕ,xN (t l ) 1 ≤ M ǫ r N +1 2ǫ φ(N ) (32) with φ(N ) = N i=0 φ i and φ i = sup t∈[t l -ǫ,t l +ǫ] (| ϕ (i) 1 (t) |).
Proof (Proof of Lemma 1) Since x is analytic, there exists θ

i,s ∈]t l -ǫ, t l + ǫ[ such that I ǫ,ϕ,x (t l ) -I ǫ,ϕ,xN (t l ) = t l +ǫ t l -ǫ (s -t l ) N +1 (N + 1)! x (N +1) (θ i,s )ϕ (i) 1 (s)ds 0≤i≤N . (33)
Furthermore, the analycity of x on ]0, a[ implies the existence of real numbers M and r ( [START_REF] Valiron | Sur les fonctions analytiques d'une variable réelle[END_REF]) such that for all t ∈ [t l -h, t l + h]

| x (N +1) (t) | (N + 1)! ≤ M r N +1 , (34)
where [t l -ǫ, t l + ǫ] ⊂ [t l -h, t l + h] ⊂]0, a[. Thus lemma 1 is proved.

Lemma 2. There exist real numbers A 0 , ..., A N depending only on the test function ϕ 1 such that:

A -1 ǫ,ϕ1,N 1 ≤ i=N i=0 A i ǫ i-1 (35)
with for all k ∈ {1, ..., N }:

-A 0 = 1 2ϕ 1 (c) -0 < A 0 = A 1 < ... < A N and -A k = A 0 k! + A 1 (k -1)! + .. + A k-1 1! . Proof (Proof of Lemma 2)
The knowledge of the first line or the last column of the matrix A -1 ǫ,ϕ1,N gives all the entries of this matrix. Putting,

A -1 ǫ,
The last column of this matrix is solution of the system: This leads to: Given any k ∈ {1, ...N -1} and suppose that for all i ∈ {1, .., k -1},

       β 0 = 1 α0 β 1 = -α1
| α i α 0 |≤ ǫ i i! (38)
| β i |≤ A i ǫ i-1 with A i such that A i A 0
is a rational number. According to (38), one gets:

| β k | ≤ ǫ k k! A 0 ǫ + ... + ǫ i i! A k-i ǫ k-i-1 + ... + ǫ 1! A k-1 ǫ k-2 , ( 41
) that is, | β k |≤ A k ǫ k-1 with A k = A 0 k! + A 1 (k -1)! + .. + A k-1 1! .
Notice that A i A 0 is a rational number because A 1 A 0 , ..., A k-1 A 0 are rational numbers. Thus, the proof of the lemma 2 is done. Finally, one can deduce the following lemma.

Lemma 3.

C ǫ,N (t l ) -X(t l ) 1 ≤ 2 i=N i=0 A i ǫ i M ǫ r N +1 φ(N ). ( 42
)
For proving the theorem 2, it is sufficient to prove that the series i=+∞ i=0

A i z i is convergent when ǫ is near to zero. This is done by the following lemma. Proof (Proof of Lemma 4) By using the lemma 2, one gets A i z i has a radius of convergence R ≥ 1 e -1 . Finally by using lemmae 3 and 4, the proof of theorem 2 is done.

A k+1 A k = A k A k 1 1! + ... + A 1 A k 1 (k -1)! + A 0 A k 1 k
Proof (Proof of Proposition 1) By using (38), one gets

A ǫ,ϕ1,N 1 ≤ i=N i=0 | α i |≤ α 0 i=N i=0 ǫ i i! ≤ α 0 e ǫ . ( 43
)
And since there exists c ∈]t l -ǫ, t l + ǫ[ such that α 0 = 2ǫϕ 1 (c), the following inequalities are obtained: A i ǫ i = 0

α 0 ≤ A ǫ,

Figure 2 :

 2 Figure 1: Ψ ǫ,ϕ,1

Figure 3 : Example 1 (

 31 Figure 3: Example 1 (14), (16), (18) are used with N = 1 (resp. N = 2). The derivative estimation method presented at section 2 will be tested with the two infinitely differentiable test functions φ (27) and φ (28).E ǫ,ϕ,N given by (26) provides in our examples a criteria for estimating the maximal value of the test functions support length. Figure (4) (resp. (5)) represents E ǫ,ϕ,1 (resp. E ǫ,ϕ,2 ). According to the figures, results on the derivative estimates are expected to degrade rapidly as soon as the support length is greater than 3. Furthermore, in this example, for having equivalent results with φ and φ, the support length of φ must be smaller than φ when it is greater than 1.5. In return, the effect of the noise will be less attenuate. In taken again the remark done at section 4.1, that is the test functions support length must be greater than one, the support length is chosen between 1 and 3 according to the value of N . Indeed, when seeking to estimate a highter order derivative, the support length must be greater for better filter the noise.For N = 1, the interval length of φ and φ is equal to 1.5 and their curves coincide. Thus, only the first derivative evaluated with φ is represented at

Figure 6 :Figure 7 :Figure 8 :

 678 Figure 4: E ǫ,ϕ,1

Figure 9 :

 9 Figure 9: Example 2

Figure 11 :

 11 Figure 10: First derivative

Figure 12 :

 12 Figure 12: Example 3

Figure 14 :

 14 Figure 13: Estimation of the first derivative

Figure 15 :

 15 Figure 15: Example 4

Figure 17 :

 17 Figure 16: Estimation of the first derivative

| β 1

 1 |≤ ǫ | β 0 | (39) Putting A 0 = 1 2ϕ1(c) and A 1 = A0 1! , this gives | β 0 |= A 0 ǫ | β 1 |≤ A 1 ǫ 0

Furthermore, A 0 = 1 2ϕ 1

 1 (c) > 0 and by simple recurrenceA k > 0 for all k ∈ IN. Since A k = A kA k > A k-1 for all k ∈ IN * -{1}.

Lemma 4 .

 4 The power seriesi=+∞ i=0 A i z i converges for | z |< R with R > 1 e -1 .

ϕ1,N 1 ≤

 1 2ϕ 1 (c)ǫ e ǫ = e ǫ α 0 (44) and (45), one gets1 ≤ χ(A ǫ,ϕ1,N ) ≤ e ǫ + 2ϕ 1 (c)e ǫ +∞ i=1 A i ǫ i (46)By using lemma 4, the series+∞ i=0 A i ǫ i is convergentwhen ǫ is near zero, then lim ǫ→0 +∞ i=1 A i ǫ i = 0 and finally lim ǫ→0 2ϕ 1 (c)e ǫ +∞ i=1

  . The results confirm that, in this example, the second test function gives better results when N = 2.

		SNRout SBR
	N = 1 φ′	36.6	24
	φ′	34.1	24.01
	N = 2 φ′	36.4	19.58
	φ′′	23.5	18.38
	φ′	33.9	23.82
	φ′′	24.3	21.56

  α0 β 0 ..........β N = -αN α0 β 0 -αN-1 α0 β 1 -... -α1 α0 β N -1 (t)dtand ϕ 1 is a continuous function there exists c ∈ ]t l -ǫ, t l + ǫ[ such that α 0 = 2ǫϕ 1 (c) and because ϕ 1 is a positive function 1 i! |

				(37)
		t l +ǫ
	Because α 0 =	t l -ǫ	ϕ 1 ǫ i i!	α 0

t l +ǫ t l -ǫ (t -t l ) i ϕ 1 (t)dt |≤


Nathalie Verdiere

Ghislaine Joly-Blanchard

Lilianne Denis-Vidal

Identifiability

Identifiability and identification of a pollution source in a river by using a semi-discretized model

teaching and research institutions in France or abroad, or from public or private research centers.


Carlos Alberto

Cavichioli Gonzaga

Marc Jungers
email: marc.jungers@cran.uhp-nancy.fr

Carlos A C Gonzaga

Jamal Daafouz

Stability analysis and stabilization of switched nonlinear systems: a switched nonlinearity-dependent Lyapunov function

Keywords: Lyapunov function, cone bounded nonlinearity, local stability and stabilization, input saturation, basin of attraction estimate, LMIs

come

Stability analysis and stabilisation of switched nonlinear systems

Introduction

Switched systems are a subclass of hybrid systems characterized by a number of subsystems and a switching strategy imposing the active mode at each time [START_REF] Liberzon | Switching in Systems and Control[END_REF]. In a large number of contributions, the switching modes are formulated as linear ones. This is, in general, only an approximation for physical systems models to be valid in a limited domain. Nonetheless, because of the actuators cannot provide unbounded magnitude signals, the resulting dynamics could be nonlinear. Thus, it is reasonable to consider some nonlinearities dependent on the state and the control input. These nonlinearities are naturally mode-dependent and should be taken into account to refine the modelling step.

The saturation nonlinearity draw attention in the classical literature of control design. Two representations of the saturation are mostly considered [START_REF] Tarbouriech | Stability and Stabilization of Linear Systems with Saturating Actuators[END_REF]): as a polytope and as a cone bounded nonlinearity. These models were also used to study switched systems subject to saturation, [START_REF] Benzaouia | Stability and control synthesis of switched systems subject to actuator saturation[END_REF][START_REF] Lu | Design of switched linear systems in the presence of actuator saturation[END_REF][START_REF] Jungers | Stabilization of discretetime swithing systems including modal nonlinearities and saturation actuators[END_REF][START_REF] Jungers | Finite L 2 -induced gain and λ-contractivity of discrete-time switching systems including modal nonlinearities and actuator saturations[END_REF]. References [START_REF] Jungers | Stabilization of discretetime swithing systems including modal nonlinearities and saturation actuators[END_REF][START_REF] Jungers | Finite L 2 -induced gain and λ-contractivity of discrete-time switching systems including modal nonlinearities and actuator saturations[END_REF] deal with more general systems including, in addition, mode-dependent cone-bounded nonlinearities, as investigated in this paper.

Stability of linear systems connected to a cone bounded nonlinearity has been widely studied, so-called Lur'e problem [START_REF] Khalil | Nonlinear Systems -Third Edition[END_REF]. In continuous-time, a dedicated Lyapunov function has been proposed, named Lur'e-type function, which is composed of a quadratic term with respect to the state and an integral term of the nonlinearity. The integral is justified because it induces the presence of the nonlinearity in the Lyapunov time-derivative. Stability conditions are, then, achieved only by managing with the sector condition. In discrete-time, the idea of using the same function has also been considered [START_REF] Haddad | Parameter-dependent Lyapunov Functions and the Discrete-time Popov Criterion for Robust Analysis[END_REF][START_REF] Kapila | A Multivariable Extension of the Tsypkin Criterion Using a Lyapunov-Function Approach[END_REF][START_REF] Gomes Da | Stability and Stabilization of Linear Discrete-Time subject to Control Saturation[END_REF], Jury and Lee 1964a,b, Sharma and Singh 1981). Nevertheless, additional assumptions about the nonlinearity slope are required to upper-bound the remaining integral term in the Lyapunov difference. At last but not least, the Lur'e function requires timeinvariance of the nonlinearity. Extending such a tool to switched systems with mode-dependent nonlinearities is difficult due to the time dependency.

In this paper, we tackle the problems of stability analysis and stabilization of discrete-time switched systems including switched nonlinearities and input saturation. The same issue has been investigated in [START_REF] Gonzaga | Stability analysis of discrete-time Lur'e systems[END_REF] in the simple case where the systems do not switch. The main contribution here is to consider a switched Lyapunov function which depends on the switched cone bounded nonlinearities without additional assumptions on their slope. Sufficient conditions for stability analysis are provided as linear matrix inequalities (LMI). On the other hand, the stabilization purpose is directly given by bilinear matrix inequalities (BMI) which are relaxed through the Finsler's Lemma obtaining sufficient conditions as LMIs.

The paper is organized as follows: Section 2 is devoted to the system description, some preliminaries and a switched nonlinearity-dependent Lyapunov function is presented. The main results related to the stability analysis problem are provided in Section 3 and a numerical example is given. Section 4 is dedicated to the stabilization issue and the same example of Section 3 is considered. Concluding remarks are given in Section 5.

Notation: Components of vector x ∈ R n are denoted x (ℓ) , ∀ ℓ = 1, ..., n. Vectors inequalities are component-wise: x ≥ 0 means x (ℓ) ≥ 0 and x ≥ y means x (ℓ)y (ℓ) ≥ 0. A (ℓ) (resp. A i,(ℓ) ) is the ℓ-th row of matrix A (resp. A i ). For two symmetric matrices A, B, A > B means that A -B is positive definite. A ′ is the transpose of A. I n (resp. 0 n ) and 0 m×n are the n-order identity (null) matrix and the m × n-order null matrix, respectively. Operator diag[x] means a diagonal matrix obtained from vector x. Symbol ⋆ stands for symmetric block in matrices. The ellipsoidal set E(M, γ) associated with M > 0 is given by {x ∈ R n ; x ′ M x ≤ γ} and the shortcut E(M ) = E(M, 1) is used. For a given matrix X ∈ R n×n , He(X) = X + X ′ .

System description and preliminaries

Consider the discrete-time switching system composed of N nonlinear modes (N ∈ N * )

x k+1 = A σ(k) x k + F σ(k) ϕ σ(k) (y k ) + B σ(k) sat(u k ), (1)
y k = C σ(k) x k , (2)
where

x k ∈ R n , u k ∈ R m , y k ∈ R p and ϕ σ(k) (•) ∈ R p
are respectively the state, input, output and the active nonlinearity vector of the system (1)-(2) at the time instant k ∈ N. This system will satisfy the following assumptions.

Assumption 2.1 The switching rule σ : N → N * takes its values in the finite set

I N = {1, • • • , N }
and is assumed to be not known a priori, but its current value σ(k) to be available in real-time.

The notation M σ(k) means that at each time k, M σ(k) takes its values in the set

{M 1 , • • • , M N } indexed by σ(k). The matrices A i , B i , C i and F i , ∀i ∈ I N , have appropriate dimensions.
Assumption 2.2 The N nonlinearities ϕ i (•) : R p → R p associated with each mode i ∈ I N are assumed to be decentralized and to satisfy their own cone bounded sector condition [START_REF] Khalil | Nonlinear Systems -Third Edition[END_REF]. These sector conditions will be mode-dependent.

This assumption, which writes ϕ i (•) ∈ [0 p , Ω i ] with a conventional abuse of notations [START_REF] Khalil | Nonlinear Systems -Third Edition[END_REF], means that ϕ i (0) = 0 and there exist N diagonal positive definite matrices

Ω i ∈ R p×p such that, ∀y ∈ R p , ∀ i ∈ I N , ϕ i,(ℓ) (y) [ϕ i (y) -Ω i y] (ℓ) ≤ 0, ∀ ℓ = 1, • • • , p. Hence, we have the following inequality: SC(ϕ i (•), y, Λ i ) = ϕ ′ i (y)Λ i [ϕ i (y) -Ω i y] ≤ 0, (3)
where Λ i △ = diag{λ q,i } q=1;••• ;p ∈ R p×p are any positive definite diagonal matrices. Note that Ω i is given by the designer and assumed to be known hereafter for each mode i ∈ I N .

It is simple to show that relation (3

) is equivalent to [Ω i y] ′ (ℓ) [ϕ i (y) -Ω i y] (ℓ) ≤ 0, ∀ℓ = 1, • • • , p, ∀y ∈ R p ; ∀i ∈ I N , which implies, with Λ i diagonal positive definite, that 0 ≤ ϕ ′ i (y)Λ i ϕ i (y) ≤ ϕ ′ i (y)Λ i Ω i y ≤ y ′ Ω ′ i Λ i Ω i y, ∀y ∈ R p . (4)
Assumption 2.3 The control input u k is subject to actuation limits, and the standard saturation function is considered:

sat(u k ) (ℓ) = sign((u k ) (ℓ) ) min ρ (ℓ) , (u k ) (ℓ) , ∀ℓ = 1, ..., m. The vector 0 m < ρ ∈ R m is assumed to be given.
Throughout this paper, the following class of control law is considered [START_REF] Jungers | Stabilization of discretetime swithing systems including modal nonlinearities and saturation actuators[END_REF]:

u k = K σ(k) x k + Γ σ(k) ϕ σ(k) (y k ), (5)
where the m × n matrix K σ(k) is a switching state feedback gain and the m × p matrix Γ σ(k) is a switching feedback gain associated with the active nonlinearity ϕ σ(k) (•). If Γ σ(k) = 0 m×p , the feedback control law requires the following assumption Assumption 2.4 The value of the active nonlinearity ϕ σ(k) (y k ) is assumed to be known as a signal [START_REF] Arcak | Circle and Popov criteria as tools for nonlinear feedback designs[END_REF] at time k, either by model estimate or measuring.

The saturation is described as a dead-zone nonlinearity Ψ(u k ) = u ksat(u k ). By replacing u k defined in (5) and using Ψ(u k ) into (1), the closed-loop model is given by

x k+1 = A cl,σ(k) x k + F cl,σ(k) ϕ σ(k) (y k ) -B σ(k) Ψ(u k ), (6)
where

A cl,i = A i + B i K i and F cl,i = F i + B i Γ i , ∀i ∈ I N .
The following set is used to characterize the dead-zone as belonging to a generalized sector condition. For given matrices

H i ∈ R m×(n+p) , i ∈ I N , consider S({H i } i∈IN , ρ) = θ ∈ R n+p ; -ρ ≤ H i θ ≤ ρ, ∀i ∈ I N . (7) Lemma 2.5 Let m × (n + p)-matrices K σ = [K σ Γ σ ] and J σ = [J 1,σ J 2,σ ], σ ∈ I N . If the vector x k = [x ′ k ϕ ′ σ (y k )] ′ is an element of S({ K σ -J σ } σ∈IN , ρ), the nonlinearity Ψ(u k ) satisfies the following sector condition SC uk = Ψ ′ (u k )U -1 σ(k) Ψ(u k ) -J 1,σ(k) x k -J 2,σ(k) ϕ σ(k) (y k ) ≤ 0, ∀σ ∈ I N , (8)
with u k defined in (5), for any diagonal positive definite matrix U σ ∈ R m×m .

Proof The proof follows the same lines of [START_REF] Tarbouriech | Stability analysis and stabilization of systems presenting nested saturations[END_REF], Lemma 1).

After presenting the studied class of systems and the necessary background, the following problems will be addressed in this paper.

Problem 2.6 (Local Stability Analysis) Given the switched gains K i and Γ i (i ∈ I N ) of the control law (5), and under Assumptions 2.1-2.4, determine a region in the state space, as large as possible, included in the basin of attraction B 0 of the system (6), allowing to conclude the local asymptotic stability for any switching rule.

Problem 2.7 (Closed-loop Stabilization) Determine the control law ( 5) such that the closed-loop system ( 6) is locally asymptotically stable for any switching rule, in a region as large as possible, included in the basin of attraction B 0 .

Instead of using a switched quadratic Lyapunov function, as considered in [START_REF] Jungers | Stabilization of discretetime swithing systems including modal nonlinearities and saturation actuators[END_REF], to study both problems, we want to investigate a function which takes the cone bounded nonlinearity ϕ i (•) into account. However, the wide-spread Lur'e-type function [START_REF] Khalil | Nonlinear Systems -Third Edition[END_REF]) cannot be employed due to the mode-dependent cone bounded nonlinearities. By taking advantage of the class of Lyapunov functions presented in our recent work [START_REF] Gonzaga | Stability analysis of discrete-time Lur'e systems[END_REF], an extension to cover the time-varying nature of these nonlinearities is possible. Moreover it will allow to improve the results of the literature considering switched quadratic Lyapunov functions and to relaxe the bilinear nature of control design sufficient conditions.

Switched nonlinearity-dependent Lyapunov functions

The following class of non-quadratic switched Lyapunov candidate functions will be considered as main tool in this paper. Let:

V : I N × R n × R p -→ R, (i, x, ϕ(C i x)) -→ x ′ P i x + 2ϕ ′ (C i x)∆ i Ω i C i x, (9)
where matrices P i ∈ R n×n are symmetric positive definite and ∆ i ∈ R p×p are diagonal positive semi-definite (i ∈ I N ). It is noteworthy that, by setting ∆ i = 0 p (∀i ∈ I N ), the switched quadratic Lyapunov function is recovered and, thus, it may be considered as a particular case of (9). From Inequalities (4), it is possible to define a lower and an upper bounds given by two switched quadratic functions

x ′ P i x ≤ V (i; x; ϕ i (C i x)) ≤ x ′ (P i + 2C ′ i Ω ′ i ∆ i Ω i C i )x, ∀i ∈ I N . (10)
The function V can be considered as a candidate, because it verifies the following properties:

• V (i; x; ϕ i (C i x)) ≥ 0, ∀x ∈ R n , ∀i ∈ I N , due to the left inequality in (10); • V (i; x; ϕ i (C i x)) = 0, if and only if x = 0, because

P i > 0 n , ∀i ∈ I N ; • V (i; x; ϕ i (C i x)) is unbounded, i.e., V (i; x; ϕ i (C i x)) → ∞ whenever x → ∞, ∀i ∈ I N .
The Lyapunov difference is defined as

δ k V = V (σ(k + 1); x k+1 ; ϕ σ(k+1) (C σ(k+1) x k+1 )) -V (σ(k); x k ; ϕ σ(k) (C σ(k) x k )). (11)
The Level Set L V (γ) associated with V (i; x; ϕ i (C i x)) and a γ > 0 is given by

L V (γ) = {x ∈ R n ; V (i; x; ϕ i (C i x)) ≤ γ, ∀i ∈ I N } = i∈IN {x ∈ R n ; V (i; x; ϕ i (C i x)) ≤ γ} , (12)
and is naturally related to the two ellipsoids intersections associated with the upper and lowerbounding switched quadratic functions i∈IN

E(P i , γ) ⊇ L V (γ) ⊇ i∈IN E(P i + 2C ′ i Ω ′ i ∆ i Ω i C i , γ). (13)
Inclusion ( 13) emphasizes an important property of function V (i; x; ϕ i (C i x)): because of taking the nonlinearities ϕ i (•) (i ∈ I N ) into account, the set L V (γ) may be non-convex and disconnected. Hence, the level set L V (1) of our function V (i; x; ϕ i (C i x)) will be considered to study Problems 2.6 and 2.7.

Remark 1 : The non-convexity and the disconnection of the level set are important properties, justified by the fact that the transition between x k and x k+1 in the discrete-time domain is not continuously done.

Local Stability Analysis

This section is dedicated to solve Problem 2.6. Consider the auxiliary matrices

M 0,i,j =       P j 0 n 0 n×p 0 n×m 0 n×p ⋆ -P i C ′ i Ωi [T i -∆ i ] J ′ 1,i 0 n×p ⋆ ⋆ -2T i J ′ 2∆ j       ; M 1,i,j =       I n -(A i + B i K i ) ′ -(F i + B i Γ i ) ′ U ′ i B ′ i 0 p×n       ; (14) M ′ 2,i,j = 0 p×n -Ω j C j (A i + B i K i ) -Ω j C j (F i + B i Γ i ) Ω j C j B i U i I p , (15)
N Gj = -G j ; 0 n×(m+2p+n) ; N Θj = 0 p×(m+2n+p) ; -Θ j , where Θ j = W j + ∆ j , ∀j ∈ I N , which will be used in the next theorem Theorem 3.1 For given matrices K i ∈ R m×n , Γ i ∈ R m×p and fixed positive diagonal matrices

U i ∈ R m×m (i ∈ I N ), consider matrices G i ∈ R n×n , J 1,i ∈ R m×n , J 2,i ∈ R m×p , symmetric positive definite matrices P i ∈ R n×n and positive semi-definite diagonal matrices ∆ i ∈ R p×p and positive definite diagonal matrices R i , Q i , T i , Θ i ∈ R p×p
, and a scalar µ in the optimization problem min µ, Gi, Pi, J1,i, J2,i, Qi, Ri, Ti, Θi, ∆i µ subject to LMIs:

µI n -P i -C ′ i Ωi [R i + ∆ i ] ⋆ 2R i > 0 (n+p) , ∀i ∈ I N , (16)
   P i C ′ i Ωi [∆ i -Q i ] (K i -J 1,i ) ′ (ℓ) ⋆ 2Q i (Γ i -J 2,i ) ′ (ℓ) ⋆ ⋆ ρ 2 (ℓ)    > 0 (n+p+1) , ∀i ∈ I N ; ∀ℓ = 1, • • • , m; (17) M 0,i,j + He M 1,i,j N Gi + M 2,i,j N Θj < 0 (2n+2p+m) , ∀(i, j) ∈ I 2 N . ( 18
)
Then, Problem 2.6 is solved with an estimate of B 0 of system (6) given by the set L V (1), induced by the Lyapunov function (9).

Proof By noting the block matrices as defined in ( 14)-( 15), Inequality ( 18) is the same as the Taylor & Francis and I.T. Consultant following one

      P j -G j -G ′ j G ′ j (A i + B i K i ) G ′ j (F i + B i Γ i ) -G ′ j B i U i 0 n×p ⋆ -P i C ′ i Ωi [T i -∆ i ] J ′ 1,i (A i + B i K i ) ′ C ′ j Ωj Θ j ⋆ ⋆ -2T i J ′ ⋆ -2 [Θ j -∆ j ]       < 0 (2n+2p+m) .
(19) Hence, if Inequality (18) holds, one has Θ j -∆ j > 0 p , which combined with ∆ j ≥ 0 p , implies Θ j > 0 p , and so W j > 0 p (∀j ∈ I N ). In addition, it means P j -G ′ j -G j < 0 n and P i > 0 n . Thus, G j is of full rank and -G ′ j P -1 j G j ≤ P j -G ′ j -G j (see [START_REF] Daafouz | Parameter dependent Lyapunov functions for discrete time systems with time varying parametric uncertainties[END_REF]). This implies, combined with the change of basis diag[G -1 j ; I n+p ; U -1 i , I p ], and a Schur complement [START_REF] Boyd | Linear Matrix Inequalities in System and Control Theory[END_REF])

    A ′ cl,i F ′ cl,i -B ′ i 0 p×n     P j     A ′ cl,i F ′ cl,i -B ′ i 0 p×n     ′ +     -P i C ′ i Ω i [T i -∆ i ] J ′ 1,i U -1 i A ′ cl,i C ′ j Ω j [W j + ∆ j ] ⋆ -2T i J ′ 2,i U -1 ⋆ -2W j     < 0 (2n+2p+m) . (20)
By multiplying Inequality (20) on the right by

[x ′ k ϕ ′ i (y k ) Ψ ′ (u k ) ϕ ′ j (y k+1 )]
′ and on the left by its transpose, with i = σ(k) and j = σ(k + 1), it leads to

δ k V -2SC uk -2SC(ϕ σ(k+1) (•), y k+1 , W σ(k+1) ) -2SC(ϕ σ(k) (•), y k , T σ(k) ) ≤ 0. (21)
By applying a Schur complement on Inequality (17), one has

P i C ′ i Ω i [∆ i -Q i ] ⋆ 2Q i - 1 ρ 2 (ℓ) ( K i -J i ) ′ (ℓ) ( K i -J i ) (ℓ) > 0 (n+p) , (22)
with K i and J i as defined in Lemma 2.5. By multiplying on the right by

x k = [ x ′ k ϕ ′ i (y k ) ]
′ and on the left by its transpose, with i = σ(k), it leads to

V (σ(k); x k ; ϕ σ(k) (y k )) + 2SC(ϕ σ(k) (•), y k , Q σ(k) ) ≥ 1 ρ 2 (ℓ) ( K σ(k) -J σ(k) ) (ℓ) x k 2 . ( 23
)
Because the ϕ σ(k) (•) verifies the sector condition, one has

V (σ(k); x k ; ϕ σ(k) (y k )) ≥ 1 ρ 2 (ℓ) ( K σ(k) -J σ(k) ) (ℓ) x k 2 , (24)
which induces the inclusion

L V (1) ⊂ S ( K σ -J σ ) σ∈IN , ρ , (25)
that is, in the set L V (1), the dead-zone sector condition (8) is verified. By multiplying Inequality (16) on the right by

x 0 = [ x ′ 0 ϕ ′ i (C i x 0 )
] ′ and on the left by its transpose, with i = σ(0), one has subject to LMIs: ( 16), ( 17), and

µx ′ 0 x 0 + 2SC(ϕ σ(0) (•), C σ(0) x 0 , R σ(0) ) ≥ V (σ; x 0 ; ϕ σ(0) (C σ(0) x 0 )). ( 26
)
  M 0,i,j -He diag[X G,j ; 0 (n+p+m) ; X Θ,j ] ⋆ ⋆ M ′ 1,i,j -(X G,j + F ′ G,j ) I n 0 n×(n+2p+m) -He (F G,j ) ⋆ M ′ 2,i,j -(X Θ,j + F ′ Θ,j ) 0 (2n+2p+m)×p I p 0 p×n -He (F Θ,j )   < 0 (3n+3p+m) , (28)
∀(i, j) ∈ I 2 N . Then, with

G i = X G,i F -1 G,i ; Θ i = X Θ,i F -1 Θ,i , ∀i ∈ I N , (29)
Problem 2.7 is fulfilled with an estimate of B 0 given by the set L V (1) .

Proof The Inequalities ( 16)-( 17) remain as stated in Theorem 3.1, implying the Inclusions ( 25) and ( 27). Moreover, if Inequality (28) holds, it means that there exists a multiplier

F i,j =     F G,j 0 (2n+p+m)×p 0 (n+2p+m)×n F Θ,j F G,j 0 n×p 0 p×n F Θ,j     , (30)
with F G,j ∈ R n×n of full rank and F Θ,j ∈ R p×p diagonal positive definite, verifying

  M 0,i,j M 1,i,j M 2,i,j M ′ 1,i,j 0 n 0 n×p M ′ 2,i,j 0 p×n 0 p   + He F i,j N Gj -I n 0 n×p N Θj 0 p×n -I p < 0 (3n+3p+m) , (31)
where N Gj and N Θi are as previously defined. By the Finsler's Lemma, this is equivalent to have

  I 2n+2p+m N Gj N Θj   ′   M 0,i,j M 1,i,j M 2,i,j M ′ 1,i,j 0 n 0 n×p M ′ 2,i,j 0 p×n 0 p     I 2n+2p+m N Gj N Θj   < 0 (2n+2p+m) (32)
verified, and which is, in fact, the Inequality ( 18). Hence, one can conclude that Theorem 4.1 allows to exhibit a solution for Problem 2.7 by considering the proof of Theorem 3.1.

Let us reconsider the Example 1 to perform the control synthesis by applying Theorem 4.1. It will illustrate some special features of this novel Lyapunov function.

Control synthesis for Example 1: The optimization problem of Theorem 4.1 provides µ = 2.5563 with the control gains given by K 1 = -0.7168 -1.0136 ; Γ 1 = -1.2923; K 2 = -1.2581 -0.7326 ; Γ 2 = -1.4650. The estimate L V (1) is depicted (solid line) in Fig. 2 with largest sphere E(µI n ) ⊂ L V (1) (dashed line) and the estimate set given by Theorem 3.1 (dotted line).

It is possible to see that our estimate presents disconnected sets and also contains the ellipsoids intersection of [START_REF] Jungers | Stabilization of discretetime swithing systems including modal nonlinearities and saturation actuators[END_REF]) (dot-dashed line). By calculating the area, one has A LV (1) = 3.2471, corresponding an improvement of 61% of with respect to the one of [START_REF] Jungers | Stabilization of discretetime swithing systems including modal nonlinearities and saturation actuators[END_REF]) and 24% regarding the area related to the estimate of Theorem 3.1 given in Section 3. Two trajectories of the closed-loop system (6) with the control gains calculated by our framework and two different realizations of arbitrary switching rules are also depicted in Fig. 2. The initial conditions are settled in different disconnected sets of our estimate. It is noteworthy that, the trajectories converge asymptotically to the origin.

Because Theorem 4.1 provides sufficient stability conditions, a question about the gap between the basin of attraction B 0 and the disconnected estimate set L V (1) may arise. Hence, we have analysed the trajectories for initial conditions located in the following region in the phase portrait: -2 ≤ x 0,(1) ≤ 2 and -2 ≤ x 0,(2) ≤ 2 (including the set L V (1)). Four switching rules were considered: It can be seen in Fig. 3 numerous points (shadow region) x 0 / ∈ L V (1) for which system ( 6) is unstable under one of the considered switching rules. In fact, these several points filling up the region between the disconnected sets, in Fig. 4, point out the suitability of the set L V (1) in estimating the basin of attraction. Thus, one has highlighted a remarkable characteristic of our Lyapunov function which is a less conservative estimate of B 0 of the system (6).

• σ a : σ a (k) = 1; ∀k ∈ N; • σ b : σ b (k) = 2; ∀k ∈ N; • σ c : σ c (2k) = 1; σ c (2k + 1) = 2; ∀k ∈ N; • σ d : σ d (2k) = 2; σ d (2k + 1) = 1; ∀k ∈ N.

Conclusion

The class of discrete-time switched systems with mode-dependent cone bounded nonlinearities subject to input saturation was addressed in this paper. As main contribution, a new switched Lyapunov function is considered in order to take the switched nonlinearities into account. The great feature of this function, in respect to the Lur'e-type one, is the possibility of covering modedependent nonlinearities and requiring only the mode-dependent sector condition assumption.

Local stability analysis and control design were tackled by optimization problems aiming to enlarge the size of the basin of attraction estimate. This estimate is induced by the Lyapunov function level set and it may be composed of disconnected sets. Stability analysis is given as LMIs conditions which were used in convex programming. For the stabilization purposes, the Finsler's Lemma was applied with a particular structure of multipliers such that sufficient LMI conditions were achieved. Numerical examples illustrate the relevance of our Lyapunov function regarding the switched quadratic function.

  Figure1. Level set L V (1) of Theorem 4.1 (solid line). Ellipsoids intersection of[START_REF] Jungers | Stabilization of discretetime swithing systems including modal nonlinearities and saturation actuators[END_REF]) (dot-dashed line) and largest sphere E(µIn) ⊂ L V (1) given by Theorem 4.1 (dashed line).

Figure 2 .

 2 Figure 2. Set L V (1) obtained by Theorem 4.1 (solid line) and Theorem 3.1 (dotted line). Ellipsoids intersection of (Jungers et al. 2009) (dot-dashed line). Two stable trajectories under two realizations of arbitrary switching rules.

Figure 3 .

 3 Figure3. Set L V (1) (solid line) and initial conditions x 0 = [x 0,(1) ; x 0,(2) ] ′ leading to unstable trajectories. For the switching rules σa (in blue), σ b (in green), σc and σ d (in red).

Figure 4 .

 4 Figure 4. Zoom of Fig. 3: Set L V (1) (solid line) obtained by Theorem 4.1. Ellipsoids intersection of (Jungers et al. 2009) (dot-dashed line). Unstable trajectories between the disconnected L V (1).

Taylor & Francis and I.T. Consultant

2,i 0 p ⋆ ⋆ ⋆ -2U i 0 m×p ⋆ ⋆ ⋆ ⋆

2,i(F i + B i Γ i ) ′ C ′ j Ωj Θ j ⋆ ⋆ ⋆ -2U i -U i B ′ i C ′ j Ωj Θ j ⋆ ⋆ ⋆

i F ′ cl,i C ′ j Ω j [W j + ∆ j ] ⋆ ⋆ -2U -1 i -B ′ i C ′ j Ω j [W j + ∆ j ] ⋆ ⋆

Acknowledgement

This work was partially supported by ANR project ArHyCo, Programme "Systèmes Embarqués et Grandes Infrastructures" -ARPEGE, contract number ANR-2008 SEGI 004 01-30011459 and by the European Community's Seventh Framework Programme (FP7/2007-2013) under grant agreement n • 257462: HYCON2 Network of Excellence "Highly-Complex and Networked Control Systems".

The nonlinearity ϕ σ(0) (•) verifying the sector condition (3), ∀σ ∈ I N , we have

By having the sector condition SC uk ≤ 0 verified inside L V (1), Inequality (21) implies δ k V < 0 (∀x = 0), proving asymptotic stability in the set L V (1). Finally, by minimizing µ, one has the maximization of the radius of the ball E(µI n ) included in L V (1).

Remark 2 : Due to the fact that the control gains K i , Γ i and the weighting matrices U i in SC uk are fixed, ∀i ∈ I N , the Inequality ( 18) is a LMI, ∀(i, j) ∈ I 2 N . Matrices U i can be set to the values obtained by the control design algorithm, like in [START_REF] Jungers | Stabilization of discretetime swithing systems including modal nonlinearities and saturation actuators[END_REF]. The optimization problem of Theorem 3.1 can be considered to solve Problem 2.7, though the control gains must be also determined and (18) becomes a BMI, which means loss of convexity.

An example is presented to show that Theorem 3.1 is able to estimate the basin of attraction, for a fixed control law, solving the Problem 2.6. The control gains are the ones associated with the optimal solution of the framework proposed in [START_REF] Jungers | Stabilization of discretetime swithing systems including modal nonlinearities and saturation actuators[END_REF].

Example 1: Consider a switching system with N = n = 2; p = m = 1; ρ = 1.5, Ω 1 = 0.7; Ω 2 = 1.3; C 1 = 0.9 0.5 ; C 2 = 1 -0.7 ;

The modal nonlinearities are ϕ 1 (y) = 0.5Ω 1 y(1exp(-0.2y 2 ); ϕ 2 (y) = 0.5Ω 2 y(1 + cos( 100y 3 )). The control synthesis algorithm in [START_REF] Jungers | Stabilization of discretetime swithing systems including modal nonlinearities and saturation actuators[END_REF]), based on a switched quadratic Lyapunov function, provides µ = 2.7140 with the control gains Γ 1 = -1.1411; K 1 = -0.9266 -0.9892 ; K 2 = -1.7245 -0.5038 ; Γ 2 = -1.1890. We wish to estimate the basin of attraction of this system. Hence, by applying Theorem 3.1, with matrices U i ∈ R m×m set to the values given by in [START_REF] Jungers | Stabilization of discretetime swithing systems including modal nonlinearities and saturation actuators[END_REF], ∀i ∈ I N , one gets µ = 2.3789.

Both estimates are shown in Fig. 1. The largest ball E(µI n ) solution of our optimization problem is also depicted (dashed line). It is possible to see that, for this example, the intersection of ellipsoids given by [START_REF] Jungers | Stabilization of discretetime swithing systems including modal nonlinearities and saturation actuators[END_REF]) (dot-dashed line) is included in the L V (1) (solid line). By calculating the area of both regions A LV (1) = 2.6135 and A E = 1.8957 of [START_REF] Jungers | Stabilization of discretetime swithing systems including modal nonlinearities and saturation actuators[END_REF]. One can conclude that the area has been increased by 37.87%.

Closed-loop Stabilization

This section is devoted to tackle Problem 2.7 based on the function (9). According to Remark 2, this task is naturally described as BMIs, which could be then sub-optimally solved by considering an iterative algorithm [START_REF] Tarbouriech | Stability and Stabilization of Linear Systems with Saturating Actuators[END_REF]). However, it is known that a solution provided by this approach is highly sensitive to the initial values of the algorithm. Hence, a relaxation is proposed by means of the Finsler's Lemma (de Oliveira and Skelton 2001) leading to sufficient conditions given by LMIs which are considered in the following statement.

symmetric positive definite matrices P i ∈ R n×n and positive semi-definite diagonal matrices ∆ i ∈ R p×p and positive definite diagonal matrices X Θ,i , F Θ,i , R i , Q i , T i ∈ R p×p , U i ∈ R m×m and a scalar µ solutions of the convex optimization problem min µ, FG,i, XG,i, Pi, Ki, J1,i,Γi, J2,i, Qi, Ri, Ti, Ui, FΘ,i, XΘ,i, ∆i µ