Ngoc Thach
email: thach.dinh@lss.supelec.fr

Dinh

Vincent Andrieu

Madiha Nadri
email: nadri@lagep.univ-lyon1.fr

Ulysse Serres
email: ulysse.serres@univ-lyon1.fr

Continuous-discrete time observer design for Lipschitz systems with sampled measurements (long version)

Keywords: Continuous discrete-time observers, reachable sets, Pontryagin Maximum Principle, LMI. I

This paper concerns observers design for Lipschitz nonlinear systems with sampled output. Using reachability analysis, an upper approximation of the attainable set is given. When this approximation is formulated in terms of a convex combination of linear mappings, a sufficient condition is given in terms of linear matrix inequalities which can be solved employing a linear matrix inequalities solver. This novel approach seems to be an efficient tool to solve the problem of observer synthesis for a class of Lipschitz systems of small dimensions.

of nonlinear systems has already been investigated in the literature. It can be traced back to Jazwinski who introduced the continuous-discrete Kalman filter to solve a filtering problem for stochastic continuous-discrete time systems (see [START_REF] Jazwinski | Stochastic processes and filtering theory[END_REF]). Inspired by this approach, the popular highgain observer introduced in [START_REF] Gauthier | A simple observer for nonlinear systems applications to bioreactors[END_REF] has been adapted to the continuous-discrete context in [START_REF] Deza | High gain estimation for nonlinear systems[END_REF]. Since then, different approaches, have been investigated by different authors. In [START_REF] Arcak | A framework for nonlinear sampled-data observer design via approximate discrete-time models and emulation[END_REF], the robustness of an observer with respect to time discretization is studied. In [START_REF] Karafyllis | From continuous-time design to sampled-data design of observers[END_REF][START_REF] Ahmed-Ali | Continuous-discrete observer for state affine systems with sampled and delayed measurements[END_REF][START_REF] Fall | State observers for a class of continuous nonlinear systems with discrete measurements[END_REF] observers are designed from an output predictor. Some other approach based on time delayed techniques have also been considered in [START_REF] Raff | Observer with sample-and-hold updating for Lipschitz nonlinear systems with nonuniformly sampled measurements[END_REF].

In our study we focus on the preliminary work of [START_REF] Deza | High gain estimation for nonlinear systems[END_REF], and consider the case in which the continuous-discrete observer is obtained in two steps: i) when no measurement is available, the state estimate is computed by integrating the model; ii) when a measurement occurs, the observer makes an impulsive correction on the state estimate.

Note that in [START_REF] Deza | High gain estimation for nonlinear systems[END_REF] the correction gain of this impulsive correction is obtained by integrating a continuous-discrete time Riccati equation. However, in the following, inspired by [START_REF] Nadri | Observer design for uniformly observable systems with sampled measurements[END_REF], we will consider a constant correction term.

In most of the works cited above, the asymptotic convergence of the estimate to the state is obtained by dominating the Lipschitz nonlinearities with high-gain techniques. However, as this is now well understood there is a trade-off between the high-gain parameter and the measurement step size. This can lead to restrictive design conditions on the sampling measurement time (see for instance [START_REF] Andrieu | Extended-Kalman-Filter-like observers for continuous time systems with discrete time measurements[END_REF]) which may prevent the use of such technique in practice.

Recently, a new observer design methodology for Lipschitz nonlinear systems with continuous time measurements has been introduced in [START_REF] Zemouche | Observers for a class of Lipschitz systems with[END_REF]. In this approach, it is shown that the differential equation satisfied by the estimation error can be rewritten in the form of a linear parameter varying system (LPV). Hence, the convergence to zero of the estimation error can be obtained by solving some specific linear matrix inequalities (see Section II for a brief summary of this approach). The aim of our paper is to extend the approach presented in [START_REF] Zemouche | Observers for a class of Lipschitz systems with[END_REF] to the discrete time measurement case. In the adopted strategy, the main problem is decomposed into two subproblems: i) Computation of an upper approximation of the reachable set for a bilinear system. This set characterizes the possible expansion of the estimation error between two measurements February 10, 2014 DRAFT when the estimate is given by integrating the model.

ii) Construction of a correction term K ensuring the convergence to zero of a quadratic error Lyapunov function. As in [START_REF] Zemouche | Observers for a class of Lipschitz systems with[END_REF], this step is performed through linear matrix inequalities (LMIs) techniques.

In this paper, we address the first problem by considering systems with specific structure: upper triangular nonlinearities (as in the preliminary version of this work in [START_REF] Andrieu | Observer design for lipschitz systems with discrete-time measurements[END_REF]), and in observability canonical form. This allows us to obtain a constructive approach to the synthesis of observers for a wide class of Lipschitz nonlinear systems while avoiding standard high-gain approach.

Consequently, we hope to build observers with larger sampling period than those obtained by usual techniques. However, the large size of the linear matrix inequality restricts the use of this approach only to low dimensional systems.

The paper is organized as follows. In Section II, after having defined the considered class of systems (in Subsection II-A), some preliminary results and in particular the approach of [START_REF] Zemouche | Observers for a class of Lipschitz systems with[END_REF] for continuous time measurements are recalled (in Subsection II-B). An approach of observer design based on reachability set analysis is presented in Subsection II-C. Section III-A, is devoted to the study of a first particular class of systems: feedforward systems. Section III-B concerns uniformly observable systems. Simulations of two academic examples are given in Section IV to illustrate the methodology proposed in this paper. Some notations:

• M ∈ R m×n means that M is a m × n matrix.

• If M ∈ R m×n , M ∈ R n×m denotes the transpose matrix.

• A symmetric matrix M ∈ R n×n is positive definite (resp. negative definite) if for all vectors v ∈ R n , v M v > 0 (resp. v M v < 0).
We then write M > 0 (resp. M < 0);

• The dot ˙denotes the derivative with respect to time.

• |v| denotes the usual Euclidean norm of x ∈ R n , and v, w the inner product of two such vectors.

• If N is a positive integer, we denote by

C 1 0 (R + , R N ) the set of mappings from R + to R N which are C 1 at zero.

II. PRELIMINARIES

A. Problem statement

The class of nonlinear systems under consideration is described by the following differential equation

ẋ(t) = Ax(t) + φ(x(t), u(t)), (1)
where x ∈ R n is the state variable, u : R → R p is a known input, A is a matrix in R n×n and φ : R n × R p → R n is a C 1 globally Lipschitz (uniformly in the input) function. In other words, the following assumption is made.

Assumption 1. For every i, j in {1, . . . n}, there exists a real number c ij 0 such that

∂φ i ∂x j (x, u) c ij , ∀ (x, u) ∈ R n × R p .
The state x of system (1) is accessible via discrete time measurements given as a sequence of m dimensional real vectors (y k ) k∈N of the form

y k = Cx(t k ), (2)
where C is a real matrix in R m×n and (t k ) k∈N is a sequence of positive real numbers defined by t k+1 = t k + δ, δ > 0 representing the sampling measurement time.

The main objective of this work is to design a global observer for system (1) which gives an estimate x that converges asymptotically to x from the knowledge of the output y k given in [START_REF] Andrieu | Extended-Kalman-Filter-like observers for continuous time systems with discrete time measurements[END_REF].

Inspired by [START_REF] Deza | High gain estimation for nonlinear systems[END_REF] and [START_REF] Nadri | Observer design for uniformly observable systems with sampled measurements[END_REF], the analysis is restricted to a specific class of continuous-discrete time observers defined by the hybrid system

   ẋ(t) = Ax(t) + φ(x(t), u(t)), t ∈ [t k , t k+1 ), x(t k ) = x(t - k ) + K(y k -C x(t - k )), (3)
where

x(t - k ) = lim t→t k , t<t k x(t). (4)
The estimation problem consists in determining a gain K such that the estimation error e(t) =

x(t) -x(t) converges asymptotically toward zero. The proposed approach is based on a result obtained in [START_REF] Zemouche | Observers for a class of Lipschitz systems with[END_REF] which is recalled in the following section.

B. An approach for continuous time measurements

In [START_REF] Zemouche | Observers for a class of Lipschitz systems with[END_REF], the authors addressed the problem of observer design for system [START_REF] Ahmed-Ali | Continuous-discrete observer for state affine systems with sampled and delayed measurements[END_REF] where the output is considered as a continuous time function and given by

y(t) = Cx(t), ∀ t > 0. (5)
In the rest of the paper, we extend this procedure when the measurement is discrete in time.

February 10, 2014 DRAFT Let R ⊂ R n×n be the set of matrices defined as

R = R ∈ R n×n | R ij = A ij ± c ij , ∀ i, j = 1, . . . , n .
Note that R is composed of 2 ρ elements where ρ is the number of c ij = 0.

One of the results obtained in [START_REF] Zemouche | Observers for a class of Lipschitz systems with[END_REF] can be summarized by the following theorem.

Theorem II. 1 ([21]). Assume that Assumption 1 is satisfied for system [START_REF] Ahmed-Ali | Continuous-discrete observer for state affine systems with sampled and delayed measurements[END_REF]. If there exist a symmetric positive definite (SPD) matrix P in R n×n and a matrix L in R n×m such that the following matrix inequalities hold:

R P + P R -C L -LC < 0, ∀ R ∈ R, (6)
then the system ẋ(t) = Ax(t) + φ(x(t), u(t))

+ P -1 L(y(t) -C x(t), (7)
is an asymptotic observer for system (1) where y is the continuous time measurement given by (5), i.e. lim t→+∞ |x(t) -x(t)| = 0.

C. An approach based on reachability analysis

The approach of observer design proposed in this work, is based on a reachable set computation. To develop this approach, let us consider system (1) and observer [START_REF] Andrieu | Observer design for lipschitz systems with discrete-time measurements[END_REF]. The estimation error

e = x -x is solution to    ė(t) = Ae(t) + ∆φ(x(t), u(t), e(t)), t ∈ [kδ, (k + 1)δ), e(kδ) = (Id -KC)e(kδ -), (8)
where the notation (4) is used and ∆φ is the continuous function defined as ∆φ(x, u, e) = φ(x, u) -φ(x -e, u).

Employing the mean value theorem, it yields the existence of n functions z i : R n → R n , i = 1, . . . , n such that the components of the function ∆φ satisfy

∆φ i (x, u, e) = ∂φ i ∂x (z i (x, e), u)e, i = 1, . . . , n. (9)
Hence, the error e(t) satisfies the following equation

ė(t) = Ae(t) + V (t)e(t), (10)
where t ∈ [kδ, (k + 1)δ) and the elements V ij (t) = ∂φ i ∂x j (z i (x(t), e(t)), u(t)) of the matrix V (t) satisfy

|V ij (t)| c ij , ∀ i, j ∈ {1, . . . , n}. (11)
In our approach, we consider the error equation [START_REF] Karafyllis | From continuous-time design to sampled-data design of observers[END_REF] as a bilinear control system where the elements V ij are bounded control inputs. The proposed observer should converge for all possible values of V ij satisfying [START_REF] Fall | State observers for a class of continuous nonlinear systems with discrete measurements[END_REF]. One way to formalize this problem is to introduce the notion of attainable set in finite time.

Let A δ (e 0 ) denote the attainable set from e 0 at time δ 0 of system [START_REF] Karafyllis | From continuous-time design to sampled-data design of observers[END_REF], i.e.,

A δ (e 0 ) = {e(δ) | δ 0, e(•) is solution to [START_REF] Karafyllis | From continuous-time design to sampled-data design of observers[END_REF] with e(0) = e 0 } .

The following theorem uses the set A δ (e) to give a condition (formulated as LMIs) guaranteeing the convergence of a continuous-discrete observer (3) for system (1)-( 2).

Theorem II.2 (Sufficient condition for observer design). Let Assumption 1 hold for system [START_REF] Ahmed-Ali | Continuous-discrete observer for state affine systems with sampled and delayed measurements[END_REF] and let δ be the sampling measurement time. Assume that there exist a finite set M of matrix functions (mappings from R + into R n×n ), a positive definite matrix P in R n×n and a vector W in R n×m such that 1

A δ (e) ⊂ Conv{M (δ)e, M ∈ M }, ∀ e ∈ R n , (12) and 
 P M (δ) (P -C W ) (P -W C)M (δ) P   > 0, ∀ M ∈ M . (13)
Then, for K = P -1 W , the estimation error given by the observer (3) converges asymptotically to zero.

Proof. For M ∈ M define N M = (Id -KC)M . On the one hand, equation [START_REF] Lohmiller | Contraction analysis of non-linear distributed systems[END_REF] and the Schur complement leads to

N M P N M -P < 0, ∀ M ∈ M .
Hence, for every M ∈ M we have

N M P N M (1 -κ M )P,
1 Conv denotes the convex closure.

February 10, 2014 DRAFT where 0 < κ M < 1 is defined as

κ M = λ min (P -N M P N M ) λ max (P ) ,
with λ max (•) and λ min (•) being the largest and the smallest eigenvalue respectively. Set κ = min M ∈M {κ M }. Since M is finite, κ ∈ (0, 1). Consequently,

N M P N M (1 -κ)P, ∀ M ∈ M .
On the other hand, take k ∈ N. The error e(•) being solution of [START_REF] Karafyllis | From continuous-time design to sampled-data design of observers[END_REF], it follows from Assumption 1 and the inclusion ( 12) that

e(t - k+1 ) ∈ Conv{M (δ)e(t k ), M ∈ M }.
With the discrete dynamics of the error in [START_REF] Gauthier | A simple observer for nonlinear systems applications to bioreactors[END_REF], it yields

e(t k+1 ) ∈ Conv{N M (δ)e(t k ), M ∈ M }.
In other words, there exist M 1 , . . . , M in M and α 1 , . . . , α in R + , with i=1 α i = 1 such that

e(t k+1 ) = i=1 α i N M i e(t k ).
We have e(t k+1 ) P e(t k+1 ) = i=1 j=1

α i α j e(t k ) N M i P N M j e(t k ) 1 2 i=1 j=1 α i α j e(t k ) N M i P N M i +N M j P N M j e(t k ) = i=1 j=1 α i α j e(t k ) N M i P N M i e(t k ) (1 -κ)e(t k ) P e(t k ) (1 -κ) k+1 e(0) P e(0). (14)
According to [START_REF] Karafyllis | From continuous-time design to sampled-data design of observers[END_REF], we have for all t ∈ [t k+1 , t k+2 ), ˙ e(t) P e(t) = e(t) (P (A + V (t)) + (A + V (t)) P )e(t).

Assumption 1 implies that for all t, n i=1

(A ij + V ij (t))e(t) ∈ Conv{Re(t), R ∈ R}.
Thus, ˙ e(t) P e(t) ηe(t) P e(t), for all t ∈ [t k+1 , t k+2 ), where η is the positive real number defined as η = max R∈R {λ max (P R + R P ), 0} λ min (P ) .

Using ( 14), the following holds e(t) P e(t) (1 -κ) k+1 exp(ηδ)e(0) P e(0),

∀ t ∈ [t k+1 , t k+2 ).
Since 0 < κ < 1, the function t → e(t) P e(t) goes to zero as t goes to infinity. The function e → e P e being proper and positive definite, it follows that e(t) converges to zero which ends the proof.

D. Some remarks on the approach of Theorem II.2

The first step of the proposed approach is the computation of an approximation of the attainable set A δ (e) for a bilinear control system. Reachability analysis has received numerous attentions in the literature; for instance, in [START_REF] Brockett | On the reachable set for bilinear system[END_REF], the author analyzes the geometry of the reachable set of bilinear systems. In [START_REF] Topunov | The convexity of the reachable set for a bilinear controllable system[END_REF], the author gave sufficient conditions guaranteeing that the reachable set of a bilinear controllable system is convex. If some results on the characterization of the reachable set are now available for low dimension systems (see for instance the recent result in [START_REF] Margaliot | Nice reachability for planar bilinear control systems with applications to planar linear switched systems[END_REF]) its characterization is still an open problem in general.

However, the novelty of the studied problematic is that the exact computation of this set is not needed. As a matter of fact, only an upper approximation in terms of the matrix functions M ∈ M as expressed in [START_REF] Löfberg | Yalmip : A toolbox for modeling and optimization in MATLAB[END_REF] is required. As it will be seen in the remaining part of the paper, for uniformly observable systems, an upper approximation can be explicitly given. Hence, for these two classes of systems our observer design strategy can be performed.

Given the matrix function set M , the second step of the design is to solve the linear matrix inequality [START_REF] Lohmiller | Contraction analysis of non-linear distributed systems[END_REF]. In fact, the usual detectability property is embedded in this inequality. For instance, it is a necessary condition that the couple (exp(Aδ), C) is detectable for this inequality February 10, 2014 DRAFT to have a solution. Note however that inequality [START_REF] Lohmiller | Contraction analysis of non-linear distributed systems[END_REF] is much stronger than this detectability condition since all Lipschitz nonlinearities have to be taken into account.

Moreover, note that this condition is not necessary for the existence of a continuous discrete observer. Indeed, as in all high-gain based designs, the nonlinearities are considered in the design as disturbances. For instance, the trivial Lipschitz system with no measurement defined on R as

ẋ = - x 1 + |x| , y = 0, (15)
defines a contraction (see [START_REF] Lohmiller | Contraction analysis of non-linear distributed systems[END_REF]) and consequently a simple copy of the system

ẋ = - x 1 + |x| ,
defines an asymptotic state observer and fits in the class of observers considered in this note (the impulsive correction is simply made with K = 0). However, our approach cannot apply to system (15) since the associated linear part (which is null) is not detectable.

In fact, it can be shown that if Theorem II.1 applies for the continuous time measurement case, then for a small sampling measurement time, the proposed approach can be applied provided that the elements of the matrix functions set M satisfy some local properties. Indeed, the link between the two matrix inequalities ( 6) and ( 13) can be expressed as follows.

Proposition II.1 (Local properties of matrices in M ). Assume that there exist P and L such that the matrix inequality (6) holds for a given set of matrices R. If the set M of matrix functions is such that every function M ∈ M is C 1 at time t = 0 and satisfies

M (0) = Id, Ṁ (0) ∈ R, (16)
then, for all δ small enough, the matrix inequality ( 13) is satisfied with the same matrix P and for W = δL.

Proof of Proposition II.1. Let M be in M and let S be the matrix defined by,

S = M (Id -P -1 W C) P (Id -P -1 W C)M -P.
We have

S = (M -Id) P M + P (M -Id) +M (-W C -C W + C W P -1 W C)M. Taking W = δL yields S(δ) δ = M (δ) -Id δ P M (δ)+ P M (δ) -Id δ + M (δ) (-LC -C L + δC L P -1 LC)M (δ),
which, using the fact that M is a C 1 matrix function and that M (0) = Id, implies

lim δ→0 S(δ) δ = P Ṁ (0) + Ṁ (0) P -LC -C L .
Note that with ( 16) and the fact that ( 6) is satisfied for all matrices in R, it yields that

lim δ→0 S(δ)/δ < 0 for all M ∈ M .
This implies that for all sufficiently small positive δ the inequality S(δ) < 0 holds. With the Schur complement, it shows that the matrix inequality ( 13) is satisfied for sufficiently small δ.

III. SOME SPECIFIC STRUCTURES ON THE SYSTEM

From Proposition II.1, it follows that a good upper approximation of the reachable set A δ (e 0 ) in terms of the matrix functions M ∈ M should be those which satisfy [START_REF] Nadri | Observer design for uniformly observable systems with sampled measurements[END_REF]. As it will be seen in this section, this is indeed the case when considering some particular class of Lipschitz nonlinear systems: The class of feedforward systems is considered in Subsection III-A. Uniformly observable systems are considered in Subsection III-B.

A. Case of feedforward systems

The approach presented in the previous section can be applied when considering a specific class of Lipschitz nonlinear systems [START_REF] Ahmed-Ali | Continuous-discrete observer for state affine systems with sampled and delayed measurements[END_REF]. Indeed, when the matrix A and the function φ have an upper triangular structure, the computation of the compact set M of matrix functions involved in the procedure, can be made explicitly.

In this subsection we consider the case in which the matrix A and the function φ satisfies the following assumption.

Assumption 2 (Feedforward systems). The matrices A and V in [START_REF] Karafyllis | From continuous-time design to sampled-data design of observers[END_REF] are upper triangular. In other words, A ij = c ij = 0 for all i, j ∈ {1, . . . , n} such that j < i.

Under Assumption 2, the following result can be stated.

Theorem III.1 (Case of feedforward systems). Consider system [START_REF] Karafyllis | From continuous-time design to sampled-data design of observers[END_REF] with control constraint [START_REF] Fall | State observers for a class of continuous nonlinear systems with discrete measurements[END_REF] where the matrix A and the positive real numbers c ij satisfy Assumption 2. Then, there exists a set M ⊂ C 1 0 (R + , R × R) of matrix functions such that for all δ 0, the inclusion ( 12) is satisfied. Moreover, the conditions ( 16) are satisfied for all M in M .

Theorem III.1 allows to construct an observer which estimates asymptotically the system state if the matrix inequality ( 13) is satisfied. Since conditions [START_REF] Nadri | Observer design for uniformly observable systems with sampled measurements[END_REF] are satisfied for all M in M , we know from Proposition II.1 that this matrix inequality has a solution provided that the LMI [START_REF] Clarke | Optimization and nonsmooth analysis[END_REF] has one and δ is sufficiently small. Note however that as it will be seen in the proof of Theorem III.1 the set M may be composed of 2

n(n+1) 2
matrices. This implies that the associted linear matrix inequality may be of very high dimension. Hence, this approach may fail to be applied for high dimensional systems. which can be described by the following lemma the proof of which is given in Appendix A at the end of the paper.

Lemma III.1 (Iterative design of M ). Let M ny be a subset of C 1 0 (R + , R ny×ny ) such that M * (0) = Id for all M * in M ny . Consider the following control system with (x, y) in R × R ny defined for all positive time t by:

ẋ(t) = u(t)x(t) + v(t), y(t) , (u, v) ∈ U × Ω, (17)
y(t) ∈ Conv{M * (t)y 0 , M * ∈ M ny }, (18)
x(0) = x 0 (19)
where U = [u min , u max ] ⊂ R and Ω = ny j=1 [v j min , v j max ] ⊂ R ny are the control sets. Then, there exists a set of C 1 matrix functions M 1+ny ⊂ C 1 (R + , R (1+ny)×(1+ny) ) such that, 1) for all positive time t,  

x(t) y(t)   ∈ Conv    M (t)   x 0 y 0   , M ∈ M 1+ny    .
2) for every M in M 1+ny , M (0) = Id.

3) for every M in M 1+ny , there exist M * in M ny and (u

* , v * ) ∈ ∂(U × Ω) 2 such that Ṁ (0) =   u * v * 0 Ṁ * (0)   . ( 20
)
Using this lemma, we give below the iterative procedure for the construction of M and prove Theorem III.1.

Proof of Theorem III.1. We consider the error system (10) with Assumption 2.

Step 0: initialization step. At the first step of the iterative procedure, the dynamics of the last component e n of the error in system [START_REF] Karafyllis | From continuous-time design to sampled-data design of observers[END_REF] Hence,

e n (t) ∈ Conv M (t)e n (0), M ∈ M 1 ,
where

M 1 = {exp ([A nn -c nn ]t) , exp ([A nn + c nn ]t)} . Note that, for M ∈ M 1 , M (0) = Id, Ṁ (0) ∈ {A nn -c nn , A nn + c nn } .
Step ( n): Apply Lemma III.1 with

x = e , y = (e ( +1) , . . . , e n ),

u = A + V , v = (A ( +1) + V ( +1) , . . . , A n + V n ), U = [A -c , A + c ], Ω = n- j=1 [A ( +j) -c ( +j) , A ( +j) + c ( +j) ].
It is also assumed that the previous step yield the existence of a set of matrix functions M n- such that ( 18) is satisfied (n y = n -). According to Assumption 2 equation ( 17) is satisfied.

Consequently, Lemma III.1 and the structure of the matrix functions in (20) yield the existence a set of matrix functions M 1+(n-) which satisfies items 1, 2 and 3 (of Lemma III.1).

Finally, at

Step n of this iterative design, the result is obtained.

B. Case of system in observability canonical form

Another context of interest is the one in which the nonlinear model ( 1) is given in observability canonical form. In this section, we assume the following assumption.

Assumption 3 (Obsevability canonical systems). The matrix A and the positive real number c ij are such that

A =           0 1 0 . . . 0 . . . . . . . . . . . . . . . 0 • • • 0 1 0 0 • • • • • • 0 1 0 • • • • • • • • • 0           , c ij = 0, ∀ (i, j) ∈ {1, . . . , n -1} × {1, . . . , n}.
In this particular context, we have V (t) = (0, . . . , 0, v(t)), where v = (v 1 , . . . , v n ). Consequently, we may rewrite system [START_REF] Karafyllis | From continuous-time design to sampled-data design of observers[END_REF] with c j := c nj as

       ėj = e j+1 , j = 1, . . . , n -1 ėn = v, e , v ∈ n j=1 [-c j , c j ], (21)
1) Statement of the main result: In this framework, the following result is established.

Theorem III.2 (Case of observability canonical systems). Let Assumptions 1 and 3 hold for system [START_REF] Ahmed-Ali | Continuous-discrete observer for state affine systems with sampled and delayed measurements[END_REF]. Then, there exist δ * > 0 and a set M ⊂ C 1 0 (R + , R n×n ) of matrix functions such that for all δ ∈ [0, δ * ], the inclusion ( 12) is satisfied. Moreover, the conditions ( 16) are satisfied for all M in M .

Before proving Theorem III.2 let us make a few remark.

As Theorem III.1, Theorem III.2 allows to construct an observer which estimates asymptotically the system state if the matrix inequality (13) is satisfied. Since conditions [START_REF] Nadri | Observer design for uniformly observable systems with sampled measurements[END_REF] are satisfied for all M in M we know from Proposition II.1 that this matrix inequality has a solution provided that the LMI (6) has one and δ is sufficiently small. Note that when C = [1, . . . , 0], (6) has a solution which is the well-known high-gain observer (see [START_REF] Gauthier | A simple observer for nonlinear systems applications to bioreactors[END_REF]). Hence, Theorem III.2 result gives a continuous discrete version of the high-gain observer. This result is not new since the system is in observable form and based on the results presented in [START_REF] Nadri | Observer design for uniformly observable systems with sampled measurements[END_REF] we know that there exists an observer for sampling-time δ small enough.

Moreover, not all systems considered in [START_REF] Nadri | Observer design for uniformly observable systems with sampled measurements[END_REF] can be addressed since general lower triangular nonlinearities are not allowed.

The interest of Theorem III.2 lies in its constructive proof and in the fact that the approach is not based on high-gain technics. Hence, we expect that δ may be chosen larger than the one allowed employing the high-gain approach of we expect to choose δ larger than the one allowed employing the high-gain approach of [START_REF] Nadri | Observer design for uniformly observable systems with sampled measurements[END_REF] (see [START_REF] Andrieu | Extended-Kalman-Filter-like observers for continuous time systems with discrete time measurements[END_REF] for a study of the limitation of the usual high-gain approach).

However, this approach may be difficult to apply due to the fact that for systems of high dimension (n large), the set M may have 2 n 2 elements.

2) Preliminaries for the proof of Theorem III.2: The idea of the proof of Theorem III.2 is to approximate the attainable set A δ (e) by an n-dimensional rectangle (a direct product of n intervals). This approximation is obtained solving some optimal control problems. Indeed, the next proposition which is proved employing Pontryagin's maximum principle gives a way to build such a rectangle from the solutions of the two following Lipschitz dynamical systems: [e - j (e 0 , δ), e + j (e 0 , δ)] ⊂ R n .

Proposition III.2. There exists δ * > 0 such that for every e 0 ∈ R n and every δ ∈ [0, δ * ], A δ (e 0 )

is contained in T (e 0 , δ).

Proposition III.2 establishes that the solutions to systems (22) and (23) are bounding state trajectories for system [START_REF] Zemouche | Observers for a class of Lipschitz systems with[END_REF]. Note that this result can be related to Müller's theorem (see [START_REF] Müller | Über das fundamentaltheorem in der theorie der gewöhnlichen differentialgleichungen[END_REF]) but is however different and has the advantage to be more realizable in practice since the bounding systems are autonomous.

The proof of Proposition III.2 is obtained solving 2n Mayer problems of the form:

           Minimize a, x(T ) Subject to ẋ = f (x, v), v ∈ Ω x(0) = x 0 (24)
where a ∈ R n is a constant vector, f is a smooth functions, and Ω ⊂ R n is compact.

The main tool to solve the Mayer problem (24) is the Pontryagin's Maximum Principle (PMP)

arising in optimal control. We refer the reader to [START_REF] Vinter | Optimal control. Systems & Control: Foundations & Applications[END_REF] for a general version of PMP adapted to the resolution of Mayer problems. The next theorem is a version of PMP that we state in our own context only.

Theorem III.3 (PMP for problem (24)). Consider the control system [START_REF] Zemouche | Observers for a class of Lipschitz systems with[END_REF]. Let H : R n ×R n ×Ω → R denote the control dependent Hamiltonian function defined by a pair (x(•), p(•))) satisfying the conditions given by the PMP is said to be an extremal (resp. an extremal pair). An extremal corresponding to λ = 0 is said to be an abnormal extremal, otherwise we call it a normal extremal.

H(x, p, v) = p, f (x, v) . If the pair (x(•), v(•)) : [0, T ] → R n × Ω is a (local)
• Notice that a Mayer problem of the form (24) do not admit abnormal extremals. Indeed, otherwise, due to the transversality condition iv and equation ii (which is linear), the pair (p(t), λ) would be trivial (everywhere) which contradicts the PMP. Consequently, p(t) is never zero as soon as the constant vector a = 0.

We are now ready to prove Proposition III.2.

Proof of Proposition III.2. First of all, notice that the result is obvious for e 0 = 0 since for every δ the attainable set of system ( 21) at e 0 = 0 is A δ (0) = {0}. Suppose from now that e 0 = 0.

Let a 1 , . . . , a n denote the canonical basis vectors of R n , and let Ω = n j=1 [-c j , c j ]. The proof is based upon the resolution of the 2n Mayer's problems (P k, ,δ ) with k ∈ {1, . . . , n}, ∈ {-, +}, δ > 0 defined as follows:

(P k, ,δ )            Minimize a k , e(δ)
Subject to [START_REF] Zemouche | Observers for a class of Lipschitz systems with[END_REF] e(0) = e 0

We note that the solution of each of these problems allow us to find the minimum and maximum of each component of the solutions of the error system [START_REF] Zemouche | Observers for a class of Lipschitz systems with[END_REF]. Consequently, to show Proposition III.2, it suffices to prove the existence of a δ * > 0 such that, for every δ < δ * , and every k = 1, . . . , n we have

min v(•)∈L ∞ (R + ,Ω) a k , e(δ) = e k (δ), (25)
where e k (δ) is the solution of ( 22 Moreover, the following transversality condition holds:

   ṗ1 = -p n v 1 ṗj = -p j-1 -p n v j , j = 2, . . . , n, (26)
p(δ) = -a k . ( 28
)
Notice that, according to Remark III.3, the adjoint vector p(•) is nontrivial. Hence, p n (•) cannot vanish on an interval. Indeed, otherwise, since (26), we get also p n-1 (t) = -ṗn (t) = 0 for every t in this interval. Then iteratively p(•) = 0 in this interval.

Similarly, for j = 1, . . . , n, e j (•) cannot vanish on an interval. Indeed, otherwise, from [START_REF] Zemouche | Observers for a class of Lipschitz systems with[END_REF] we get also e j+1 (t) = ėj (t) = 0 for t in this interval and iteratively, we obtain e (•) = 0 for = j, . . . , n. If j 2, we have e

(j)
1 (•) = 0 and therefore e 1 , . . . , e j-1 are polynomials of degree j -2, . . . , 0 respectively. On the other hand, in this interval we have ėn (t) = j-1 =1 v (t)e (t) = 0. Thus, with (27) we get that, for almost all t in this interval

p n (t) j-1 =1 e (t)v (t) = |p n (t)| j-1 =1 c |e (t)| = 0.
Now, since the e j (•)s are polynomials of different degrees and p n (•) cannot vanish on an interval, we conclude that e 1 (t) = • • • = e j-1 (t) = 0 for a.e. t. Consequently, e(t) = 0, which is a contradiction with e 0 = 0 since the system ( 21) is linear.

Summing up, we obtain that p n (•)e j (•), j = 1, . . . , n cannot vanish on an interval, which implies that for almost all t in R + v j (t) = |c j | sign(p n (t)e j (t)).

Note that if we show that there exists a time δ * > 0 such that p n does not change its sign if t < δ * , then with the transversality condition we get that

v j (t) = |c j | sign(p n (δ * )) sign(e j (t)) = -|c j | sign(e j (t)) ∀ t ∈ [0, δ * ]. (29)
In this case equation ( 25) is thus obtained and the proof of Proposition (III.2) is complete.

Consequently, to finish the proof it remains to show that p n (•) does not change its sign if t is small enough. This property is obtained from the following technical lemma, the proof of which is given in Appendix B.

Lemma III.4. There exist a δ 1 positive and a positive coefficient d such that for every t in [0, δ 1 ] and for every δ > 0, we have

|p n (δ -t)| < d t n-k .
From Lemma III.4, (26) and using the fact that |v 1 (t)| c 1 we get,

| ṗ1 (δ -t)| c 1 d t n-k , ∀ t ∈ [0, δ 1 ].
Integrating this inequality, we obtain

|p 1 (δ -t)| c 1 d n -k + 1 t n-k+1 + |p 1 (δ)|, ∀ t ∈ [0, δ 1 ].
By integrating successively the previous inequality, and using the transversality conditions (28), it yields that there exists a positive number d k-1 such that:

|p k-1 (δ -t)| d k-1 t n-k+1 , ∀ t ∈ [0, δ 1 ]. (30)
On the other hand, we have the following inequality for all t in [0, δ 1 ]:

p k (δ -t) p k (δ) - t 0 |p k-1 (δ -s)| + c k |p n (δ -s)|ds.
From ( 28) and (30), it follows that for t in [0, δ 1 ]:

p k (δ -t) 1 -d k-1 (δ -t) n-k+2 -c k d t n-k 1 -d k t n-k ,
where d k is positive. Now, by integrating the previous inequality, we obtain for all t in [0, δ 1 ]:

p k+1 (δ -t) t -d k t n-k+1 -c k+1 d t n-k t -d k+1 t n-k ,
where d k+1 is positive. Proceeding in the same manner successively, we obtain a positive constant

d n-1 such that p n-1 (δ -t) t n-k-1 -d n-1 t n-k .
Comparing the two degrees of monomials of the last equation, we obtain the existence of a time

δ 2 in [0, δ 1 ] such that p n-1 (δ -t) > 0, ∀ t ∈ [0, δ 2 ].
Finally, we have ṗn v n p n + p n-1 and p n (δ) = 0. Thus, we get for every t in [0, δ 2 ],

p n (δ -t) = δ-t δ exp δ-t s v n (r)dr p n-1 (s)ds > 0.
Proceeding in the same way, one infers that there exists a sufficiently small time δ * , such that for every t in [0, δ * ] and every k ∈ {1, . . . , n}, p n (δ -t) is positive (resp. negative) if = -(resp.

= +).

From the structure of minimizing control, given in (29), it is concluded that δ < δ * .

Hence, v k (t) = -|e k (t)|.

February 10, 2014 DRAFT

3) Proof of Theorem III.2: With the approximation of the attainable set A δ (e 0 ) by the rectangle T (e 0 , δ) (see Proposition III.2), we can now give the proof of Theorem III.2. In order to do that, let us define the Clarke's gradient in the direction of x j (denoted by ∂ j C f (x)) of a vector function f to be the generalized gradient 3 of the function x j → f (x 1 , . . . , x n ).

Since ( 22) and ( 23) are globally Lipschitz, e + (e 0 , δ) and e -(e 0 , δ) are globally Lipschitz 4 .

Consequently, the Clarke gradients ∂ j C e + i (e, δ) and ∂ j C e - i (e, δ) exist for every i = 1, . . . , n. Let S denote the unit sphere of R n . Introduce the two functions

m - ij (δ) = min ν∈S min ∂ j C e + i (ν, δ), ∂ j C e - i (ν, δ) , m + ij (δ) = max ν∈S max ∂ j C e + i (ν, δ), ∂ j C e - i (ν, δ) .
Let M be the set of 2 n 2 matrices M given by

M ij (δ) ∈ m - ij (δ), m + ij (δ) . (31)
To complete the proof of Theorem III.2 we have to show the following two properties:

1) the conditions ( 16) are satisfied for all M in M ;

2) there exists δ * > 0 such that for all 0 δ δ * , the inclusion ( 12) is satisfied.

About the conditions [START_REF] Nadri | Observer design for uniformly observable systems with sampled measurements[END_REF]. First of all, note that since e - i (ν, 0) = e + i (ν, 0) = ν, it yields that M (0) = Id. For the local property of the time derivative, let M + ij (δ) be the function defined by

M + ij (δ) = max ν∈S ∂ j C e + i (ν, δ).
Note that by definition of the generalized gradient, for each δ there exists a sequence of (ν (δ)) ∈N converging toward ν * (δ) in S such that for all in N, ∂e + i ∂e j (ν (δ), δ) is well defined and

M + ij (δ) = lim →+∞ ∂e + i ∂e j (ν (δ), δ).
3 Following F. Clarke [START_REF] Clarke | Optimization and nonsmooth analysis[END_REF], we define the generalized gradient of a scalar function f at x0 as the convex envelop of all possible limits of derivatives of f at points xn ∈ R n , xn → x0. Note that, in general, ∂ C f (x0) is a set. 4 To see this, note that for (ea, e b ) in R n , we have On another hand, with v j the unit vector with 1 at the jst component, we have for ( , δ) in

N × R + ∂e + i ∂e j (ν (δ), δ) = lim r→0 e + i (ν (δ) + rv j , δ) -e + i (ν (δ), δ) r .
However, we have for (r, , δ

) in R + × N × R + e + i (ν (δ) + rv j , δ) -e + i (ν (δ), δ) = r Id ij + δ F + i (ν (δ) + rv j ) -F + i (ν (δ)) + G + i (ν (δ), r, δ),
where

G + i (ν (δ), r, δ) = δ 0 F + i (e + i (ν (δ) + rv j , s)) -F + i (ν (δ) + rv j ) -F + i (e + i (ν (δ), s)) + F + i (ν (δ))ds. ( 32
)
Note that F + i being a globally Lipschitz vector field with Lipschitz constant denoted c max , it yields,

G + i (ν (δ), r, δ) δ 0 c max e + i (ν (δ) + rv j , s) -ν (δ)
-rv j -e + i (ν (δ), s) + ν (δ) ds, (33) which gives,

G + i (ν (δ), r, δ) δ 0 c max s 0 F + i (e + i (ν (δ) + rv j , τ )) -F + i (e + i (ν (δ), τ ))dτ ds c max δ 0 s 0 F + i (e + i (ν (δ) + rv j , τ )) -F + i (e + i (ν (δ), τ )) dτ ds c 2 max δ 0 s 0 e + i (ν (δ) + rv j , τ ) -e + i (ν (δ), τ ) dτ ds.
Moreover, note that for all and δ, we have

lim r→0 F + i (ν (δ) + rv j ) -F + i (ν (δ), δ) r max ν∈S ∂ j C F + i (ν).
This implies,

M + ij (δ) Id ij +δ max ν∈S ∂ j C F + i (ν) + c 2 L δ 2 m+ ij (δ),
where

m+ ij (δ) = max s∈ ν∈S ∂ j C e + i (ν,δ)
|s|.

This gives finally

lim δ→0 M + ij (δ) -M + ij (0) δ max ν∈S ∂ j C F + i (ν) A ij + c ij .
A similar property may be obtain for the minimum. In other words, we get

lim δ→0 M + ij (δ) -M + ij (0) δ A ij -c ij .
This implies that

lim δ→0 M + ij (δ) -M + ij (0) δ ∈ Conv {R ij , R ∈ R} ,
where R is the set of matrices defined in Section II-B. From this we conclude that the set M of matrix functions defined in (31) satisfies

Ṁij (0) ⊂ Conv{R ij , R ∈ R}.
About the inclusion [START_REF] Löfberg | Yalmip : A toolbox for modeling and optimization in MATLAB[END_REF]. To finish the proof of the Theorem, it remains to show that ( 12) is also satisfied. The systems ( 22) and (23) being homogeneous of degree zero, their solutions are homogeneous of degree one. We deduce that for all λ > 0, all e 0 ∈ R n , e + (λe 0 , δ) = λe + (e 0 , δ) , e -(λe 0 , δ) = λe -(e 0 , δ) .

Consider a vertice S(e 0 , δ) of T (e 0 , δ) and assume without loss of generality that S i (e 0 , δ) = e + i (e 0 , δ). Given e 0 = 0 in R n , let ν 0 = e 0 |e 0 | in S. We can also introduce ν a sequence in R n such that lim →+∞ ν = ν 0 and e + (ν , δ) is C 1 in its first argument. We have By the homogeneity property, we have,

e + i (ν , δ)|e 0 | = |e 0 | n j=1 ∂e + i ∂e j (ν , δ)(ν ) j .
Hence, it yields

e + i (e 0 , δ) = lim →+∞ n j=1 ∂e + i ∂e j (ν , δ)(ν ) j |e 0 | (34) = lim →+∞ n j=1 ∂e + i ∂e j (ν , δ)(e 0 ) j . (35)
Consequently, from the definition of the set of matrices M we obtain

S i (e 0 , δ) = e + i (e 0 , δ) = lim →+∞ n j=1 ∂e + i ∂e j (ν , δ)(ν ) j |e 0 | ∈ Conv{M i (δ)e 0 , M ∈ M }.
By iterating this procedure on each component of the vertice S, it yields that S(e 0 , δ) ∈ Conv{M i (δ)e 0 , M ∈ M }.

By iterating this procedure for each vertice of the multidmensional rectangle which is a convex set, finally we obtain T (e 0 , δ) ⊂ Conv{M (δ)e 0 , M ∈ M }.

The result is then obtained using Proposition III.2.

IV. ILLUSTRATIVE EXAMPLES

We consider the following simple model of a pendulum ẋ1 = x 2 , ẋ2 = sin x 1 , where

x 1 , x 2 denote the angle between the pendulum and the vertical axis and the pendulum speed, respectively, and

y k = x 1 (t k ), t k = t k-1 + δ. The associate estimation error equation is given by ė1 = e 2 , ė2 = ve 1 , v ∈ [-1, 1].
Following the proposed approach, we get that the approximation by the rectangle denoted by .

Note that by considering the minimum and maximum of each element, we can reduce the number of matrices to 16.

We also note that the obtained set of matrices satisfies the local properties ( 16) of Proposition II.1. The system is in uniformly observable form, so we know that for small values of δ, the assumptions of Theorem II.2 are satisfied.

Employing the Yalmip package ( [START_REF] Löfberg | Yalmip : A toolbox for modeling and optimization in MATLAB[END_REF]) in Matlab in combination with the solver Sedumi ( [START_REF] Sturm | Using SeDuMi 1.02, a MATLAB toolbox for optimization over symmetric cones[END_REF]), it can be checked that the LMI ( 13) is satisfied5 for δ 0.668. The observer gain for δ = 0.668

is K = [-1, -1.8361] .
It is interesting to notice that the bound obtained is much larger then the one obtained employing the usual high-gain approach as exposed in [START_REF] Andrieu | Extended-Kalman-Filter-like observers for continuous time systems with discrete time measurements[END_REF]. The figure 1 represents a simulation of the observer using an integration algorithm of the model with a semi-implicit integration step of 0.01.

V. CONCLUSION

In this paper, the problem of designing an observer for nonlinear systems with discrete time measurements and globally Lipschitz nonlinearities is addressed. A solution based on the synthesis of an upper approximation of a reachable set have been presented. When this approximation is given in terms of a convex combination of linear mappings, a sufficient condition of the global convergence of the proposed observer is obtained in terms of a linear matrix inequality. The good performances obtained on an illustrative example demonstrate that the proposed approach is an efficient tool.

This implies that, for all y 0 in R ny there exist γ 1 and γ 2 in Γ * such that for all M ∈ M ny and all (s, t) With the definition of Γ * , (20) is satisfied. This concludes the proof.

B. Proof of Lemma III.4:

At first, we note that because of the transversality condition, we have p n (δ) = 0. The application t → p n (t) being continuous, for every > 0 there exists a time δ 1 such that |p n (δ -t)| < for every t in [0, δ 1 ]. Now, assume that there exist a C ∞ function α : R + → R + and an integer q α such that for every t in [0, δ 1 ]

|p n (δ -t)| α(t), α ( ) (0) = 0, = 0, . . . , q α -1, α (qα) (0) > 0.

Thus, using (26) and the fact that Using the transversality conditions (28), this inequality becomes

|p n-1 (δ -t)| n-1 =1 c α n-(t) + t n-k-1 (n -k -1)! , ∀ t ∈ [0, δ 1 ],
with k being the one considered in the proof of Proposition III.2. Therefore, using Gronwall's lemma, this gives that the last component and for every t in [0, δ 1 ]: Note that β satisfies β ( ) (δ) = 0, = 0, . . . , q β -1, β (q β ) (δ) > 0, with q β = min{q α + 2, n -j}. Reiterating the procedure with the function β instead of the function α, we obtain the existence of a C ∞ function α such that for every t in [0, δ 1 ]

|p n (δ -t)|
|p n (δ -t)| α(t), α ( ) (0) = 0, = 0, . . . , n -k -1, α (n-k) (0) > 0.

This completes the proof of the lemma.

  ė = F + (e), F + (e) = e 2 , . . . , e n , n j=1 c j |e j | , (22) ė = F -(e), F -(e) = e 2 , . . . , e n , -n j=1 c j |e j | . (23) Let e + (e 0 , •) (resp. e -(e 0 , •)) denote the solution of (22) (resp. (23)) emanating from a point e 0 at time 0 and let T (e 0 , δ) be the n-dimensional rectangle defined by T (e 0 , δ) = n j=1

  ) or (23) depending on the value of . The rest of the proof is devoted to the resolution of the 2n optimal control problems (P k, ,δ ) to prove the existence of such a δ * . Let us apply Theorem III.3 to problem (P k, ,δ ). The control dependent Hamiltonian associated to (21) reads H(e, p, v) = n-1 j=1 p j e j+1 + p n n j=1 e j v j , p ∈ R n . Assume that (e(•), v(•)) is an extremal pair associated with the minimization problem (P k, ,δ ), then, according to Theorem III.3, for a.e. t ∈ [0, δ]: the adjoint system reads

e 0 F

 0 + (ea, δ) -e + (e b , δ) δ + (e + (ea, s)) -F + (e + (e b , s)) ds δ 0 cmax e + (ea, s) -e + (e b , s) ds exp(cmax δ) |ea -e b | .

  e + (e 0 , δ) = e + (ν 0 , δ)|e 0 | = lim →+∞ e + (ν , δ)|e 0 |. February 10, 2014 DRAFT

Figure 1 .

 1 Figure 1. Evolution with time of the state componant x 2 and its estimation x2 .

γ 1 ( 0 γ 2 MṀM

 102 s, t), y 0 expt s u(r)dr (v M )(s)y 0 γ 2 (s, t), y 0 . dr (v M )(s)y 0 ds t (s, t), y 0 ds, t), y 0 , γ ∈ Γ * .Finally, M 1+ny is the set of 2 × 4 ny × #M ny (#M ny denotes the cardinal of M ny ) * ∈ M ny , γ ∈ Γ * , u * ∈ {u min , u max } .Moreover, notice that for every M ∈ M 1+ny * ∈ M ny , γ ∈ Γ * , u * ∈ {u min , u max } .

  |v 1 (t)| c 1 we have | ṗ1 | c 1 |p n | and it is obtained that for every t in [0, δ 1 ] |p 1 (δ -t)| c 1 t 0 |p n (δ -s)|ds + |p 1 (δ)| c 1 α 1 (t) + |p 1 (δ)|, where α 1 (t) = t 0 α(s)ds. Similarly, we have | ṗ2 | |p 1 | + c 2 |p n | and it yields,

|p 2 ( 0 α - 1

 201 δ -t)| c 1 α 2 (t) + |p 1 (δ)|t + c 2 α 1 (t) + |p 2 (δ)|, ∀ t ∈ [0, δ 1 ],where α 2 (t) = t 0 α 1 (s)ds. Then, by integrating successively along the trajectories of the system we obtain the inequality for every t in [0, δ 1 ]|p n-1 (δ -t)| n-1 =1 c α n-(t) + |p (δ)| t n--1 (n --1)! ,where the functions α , = 1, . . . , n -1 are defined by the following iterative procedure:α (t) =t (s)ds, = 1, . . . , n.

e

  cn(s+t-δ) |p n-1 (s)|ds,β(t),February 10, 2014 DRAFTwhere β : R + → R + is the C ∞ function defined by:β(t) = e cnδ 1 +1 (t) + t n-k c n (n -k)! .

  is considered. According to Assumption 2, this component satisfies ėn = [A nn + V nn ]e n , with the control constraint |V nn | c nn . The solutions of this system are given by e

n (t) = exp t 0 (A nn + V nn (s))ds e n (0).

  By integrating solutions of systems (22) and (23), we obtain a set M of 64 matrix functions M satisfying the assumptions of Theorem III.2, where M 11 ∈ {cos δ, cosh δ} , M 21 ∈ {-sin δ, sinh δ} , M 12 ∈ {sin δ, sinh δ, cosh δ tan δ, tanh δ cos δ} , M 22 ∈ cos δ, cosh δ, 1+sinh δ sin δ

	cos δ	, 1+sin δ sinh δ cosh δ

T (e 0 , δ) is possible for δ * = √ 3 -1.

February 10, 2014 DRAFT

∂(U × Ω) denotes the boundary of U × Ω February 10, 2014 DRAFT

The Matlab files can be downloaded from https://sites.google.com/site/vincentandrieu/ February 10, 2014 DRAFT

This work was supported by ANR LIMICOS contract

A. Proof of Lemma III.1

With the variation of constant methode, the solution of the system ( 17)-( 18) is given by:

Let Ω * = ny j=1 {v j min , v j max }. Let Γ * denote the set of 4 ny vector functions (s, t) → γ(s, t) ∈ R ny defined as follows. For γ ∈ Γ * , the j th component of γ(s, t) (denoted by γ j (s, t)) satisfies,

, where u 1 , u 2 ∈ {u min , u max }. Note that for every j in {1, . . . , n y }, it holds:

and,


Philippe Marchal
email: marchal@math.univ-paris13.fr

Permutations with a prescribed descent set

We give a formula to compute the number of permutations with a prescribed descent set in quadratic time. We give the generating function of the number of permutations with a periodic descent set. We introduce an algorithm generating uniformly distributed random permutations with a prescribed descent set.

Statement of the results

Let σ be permutation of {1, 2, . . . N }. The descent set of σ, D(σ), is the set of integers i ∈ [1, N -1] such that σ(i) > σ(i + 1). For instance, D(Id) = ∅. If σ is an alternating permutation, that is, if σ(1) > σ(2) < σ(3) > σ(4) . . .

then D(σ) = [1, N -1] -2N.
The aim of this paper is to study permutations of {1, 2, . . . N } with a prescribed descent set A ⊂ [1, N -1]. We first give a method to compute this number efficiently, namely, using O(N 2 ) elementary computations, while the classical formula based on an inclusion-exclusion argument uses an exponential number of elementary computations. Next, given a periodic pattern of ascents and descents, we give a complete expression of the exponential generating function of permutations whose descent set follows this pattern. Finally, we introduce an algorithm generating at random, with uniform probability, permutations with a given descent set. Our first theorem is the following:

Theorem 1 Fix an integer N ≥ 2 and a set A ⊂ [1, N -1]. Define by descending induction a sequence of polynomials (f i , 1 ≤ i ≤ N ) as follows:

• f N = 1 • If i ∈ A, f i (x) = x 0 f i+1 (y)dy • If i / ∈ A, f i (x) = 1 x f i+1 (y)dy
Let Q N (A) be the number of permutations of {1, 2, . . . N } with descent set A. Then

Q N (A) = N ! 1 0 f 1 (y)dy
The reason why we use descending induction, rather than regular induction, will appear in Theorem 3. An easy inclusion-exclusion argument (see for instance [B], Chapter 1) gives

Q N (A) = B⊂A (-1) |A|-|B| g N (B) (1)
where the function g N is defined as follows. If B = ∅, then g N (B) = 1. Otherwise, there exists k ∈ [1, N -1] such that B has the form B = {i 1 < i 2 . . . < i k } and in that case

g N (B) = N ! (N -i k )!(i k -i k-1 )! . . . (i 2 -i 1 )!i 1 !
However, in order to use (1), one needs to compute 2 N -|A| terms. By comparison, Theorem 1 only requires the computations of the coefficients of N polynomials of degree 0 to N -1, which means computing O(N 2 ) coefficients. Moreover, each coefficient can be computed using only one elementary computation, except for the constant coefficients. It is easily seen that calculating the constant coefficient of f d requires Nd elementary computations. Hence the number of elementary computations in Theorem 1 is O(N 2 ). Another efficient method to compute Q N is given by Viennot [V].

Theorem 1 also gives access to comparison results. For a permutation σ of

[1, N ], say that m ∈ [2, N -1] is a local extremum if either σ(m -1) < σ(m) > σ(m + 1) or σ(m -1) > σ(m) < σ(m + 1). Then we have Corollary 1 For B ⊂ [2, N -1], let F (B) be the number of permutations of [1, N ] with set of local extrema B. Then F : P({2, 3, . . . n -1}) → N is an increasing function.
Of course, such results are more difficult to derive using alternating sums such as (1). The fact that the maximum of F is achieved for alternating permutations was first observed by Niven [N], see also [BR]. Our next result deals with the case of a periodic descent set.

Theorem 2 Let p ≥ 2 be an integer, fix a set A ⊂ {1, 2, . . . p} and put

A ∞ = ∪ ∞ n=0 A + np
For each integer n ≥ 1, let d n be the number of permutations of {1, . . . , n} with descent set

A ∞ ∩ [1, n -1] Let ω 1 , . . . ω p be the p-th roots of 1 (resp -1) if p -|A| is even (resp. odd). For k, l ∈ [1, p], put g k,l (t) = ∞ n=0 (ω k t) np+l (np + l)!
Then there exists a rational function R A ∈ Z(X 1 , . . . , X p(p+1) ), which can be computed explicitely using equations ( 7) through (17), such that

∞ n=1 d n t n n! = R A (ω 1 , . . . , ω p , g 1,1 (t), . . . , g p,p (t))
Equations ( 7) to ( 17) are given in Section 3. We shall see that these expressions involve matrices of size µp × µp, where µ = p -|A| is the number of ascents in a period. Theorem 2 generalizes the classical theorem by André [A] for alternating permutations:

∞ n=0 d n t n n! = 1 + sin(t) cos(t)
where we agree that d 0 = 1. More recent results were established by Carlitz (see [C] for instance) and Mendes-Remmel-Riehl [MRR] in the case µ = 1. We shall see how to recover these results from our method in Section 3.

For the first-order asymptotics, see [LM] and [BHR]. For every set A, one can compute R A using (7) to (17) and derive the asymptotics from the zeros of the denominator, using the classical tools of analytic combinatorics. See for instance the analysis of alternating permutations as Example IV.35 in [FS]. However, our approach does not provide the general form of the denominator.

Finally, we introduce an algorithm generating uniform random permutations with a given descent set. We begin by a simple remark, which was already used in [ELR]

among others. Let Y = (Y 1 , . . . Y N ) be a sequence of distinct reals in [0, 1]. We can construct from Y a permutation σ Y as follows. Let k 1 ∈ {1, 2, . . . N } be the integer such that Y k1 is minimal in {Y 1 , . . . Y N } and put σ Y (k 1 ) = 1. Then, let k 2 ∈ {1, 2, . . . N } -{k 1 } be the integer such that Y k2 is minimal in {Y 1 , . . . Y N } -{Y k1 }, put σ Y (k 2 ) = 2
and so on. To recover σ Y from Y , one can use a sorting algorithm.

One can define the descent set D(Y ) of a sequence of reals Y as for a permutation and clearly, D(Y ) = D(σ Y ). Moreover, if Y is chosen according to the Lebesgue measure on [0, 1] N , then σ Y is uniformly distributed over all permutations of {1, 2, . . . N }. As a consequence, if Y is chosen uniformly at random among all sequences with descent set A, that is, if its density with respect to the Lebesgue measure on [0, 1] N is

C1 {D(Y )=A} (2)
for some constant C > 0, then σ Y is uniformly distributed over all permutations with descent set A. Therefore, we want to find an algorithm constructing a random sequence with descent set A.

Algorithm Using iid random variables (U 1 , . . . U N ), uniform on [0, 1], and the functions f i defined in Theorem 1, construct a sequence Y = (Y 1 , . . . Y N ) as follows:

• Y 1 is the real in [0, 1] such that Y1 0 f 1 (y)dy = U 1 1 0 f 1 (y)dy • For i ∈ [1, N -1], Y i+1 is the only solution in [0, 1] of the equation f i (Y i+1 ) = U i+1 f i (Y i ) (3)
Finally, recover the permutation σ Y from Y by sorting.

Theorem 3 The algorithm described above yields a random permutation of {1, 2, . . . N }, uniformly distributed over all permutations with descent set A.

To see that the algorithm is well-defined, suppose for instance that i ∈ A.

Then (3) becomes

Yi+1 0 f i+1 (y)dy = U i+1 f i (Y i ) (4)
Since the function f i is nonnegative on the interval [0, 1], (4) clearly has a unique solution on [0, 1].

We shall not study in detail the complexity of the algorithm. Nevertheless, let us make a few comments.

First, remark that a simple rejection algorithm would have exponential complexity. Indeed, if we draw a permutation at random , the most likely profile is the alternating case, and the probability for a random permutation of length N to be alternating is ∼ 4 π (2/π) N (see Example IV.35 in [FS]). Therefore, the average number of times one has to draw a permutation before finding one with the prescribed descent set is at least c(π/2) N .

Using our algorithm, we need to compute the functions f i , which can be done using O(N 2 ) elementary computations, as already seen. The last step, namely sorting the sequence Y to deduce σ Y , has an average complexity O(N log N ) if one uses randomized quicksort. What is more intricate is to evaluate the time needed to generate the variables Y i , which amounts to solving N equations of the form (4). We shall not discuss this from a theoretical point of view. However, some experimental results are given at the end of the paper. Typically, our algorithm can generate permutations of length a few thousands.

In the particular case of alternating permutations, other algorithms with a better complexity can be found in [BRS] and [M].

We now proceed to the proof of the results. We first prove Theorem 3 and deduce Theorem 1 and Corollary 1 in Section 2. Section 3 is devoted to the proof of Theorem 2. Finally, a possible implementation of our algorithm in Maple is given in Section 4.

Proof of Theorems 1 and 3

First, observe that there are obvious symmetries in the problem. If one replaces

σ(i) with N + 1 -σ(i) for each i, then A is replaced with its complement. If σ(i) is replaced with σ(N + 1 -i) for each i, then A is replaced with N -A.
Let us prove Theorem 3. One can re-formulate the algorithm, using the notion of conditional density of a random variable, as follows:

• Y 1 has density f 1 / 1 0 f 1 (y)dy. • If i ∈ A, conditionally on Y i , Y i+1 has density 1 [0,Yi] f i+1 /f i (Y i ). • If i / ∈ A, conditionally on Y i , Y i+1 has density 1 [Yi,1] f i+1 /f i (Y i ).
As a consequence, the density of the sequence Y = (Y 1 , . . . , Y N ) with respect to the Lebesgue measure on [0, 1] N is equal to the telescopic product

f 1 (Y 1 ) 1 0 f 1 (y)dy × f 2 (Y 2 ) f 1 (Y 1 ) × . . . × f N (Y N ) f N -1 (Y N -1 ) 1 {D(Y )=A} = 1 1 0 f 1 (y)dy 1 {D(Y )=A}
which is exactly the form required in (2). Therefore, σ Y is uniformly distributed over all permutations with descent set A. This proves Theorem 3. Let us deduce Theorem 1. If a permutation of {1, . . . n} is drawn uniformly at random, the probability that its descent set is A is

1 0 dy 1 . . . 1 0 dy N 1 {D(Y )=A}
But this probability is also equal to

Q N (A) N !
Using the formula for the density of Y , we find that

1 1 0 f 1 (y)dy
h i , 1 ≤ i ≤ N by h N = 1 and h i (x) = f i (x) if i ∈ A while h i (x) = f i (1 -x) if i / ∈ A.
Remark that for i ≤ N -1, the h i are increasing functions on [0, 1] such that h i (0) = 0. Besides, one can directly define the h i by descending induction as follows. First, h N = 1 and then, if i / ∈ B,

h i (x) = x 0 h i+1 (y)dy (5) while if i ∈ B, h i (x) = x 0 h i+1 (1 -y)dy (6)
Likewise, define the polynomials hi , 1

≤ i ≤ N by hN = 1, hi (x) = fi (x) if i ∈ Ã and hi (x) = fi (1 -x) if i /
∈ Ã. Since B = B ∪ {m}, we have hi = h i for i > m. Moreover, using ( 5) and (6), we get that hm (x) > h m (x) for every x ∈ (0, 1] and by induction, for every i < m and every x ∈ (0, 1], hi (x) > h i (x). Integrating this inequality for i = 1, we get the result.

Remark

The coefficients of the polynomials f i can be computed by induction. Put

f N -i (x) = i k=0 a (i) k x k First, a (0) 0 = 1. Next, if N -i ∈ A, then a (i+1) 0 = 0 and for k ≥ 0, a (i+1) k+1 = a (i) k k + 1 On the other hand, if N -i / ∈ A, then for k ≥ 0, a (i+1) k+1 = - a (i) k k + 1 and a (i+1) 0 = i k=0 a (i) k k + 1 Therefore, if N -i / ∈ A, a (i+1) 0 = i j=0 a (i-j) 0 (j + 1)! (-1) |A c ∩[N -i+1,N -1]|
By recursion, this leads to (1), which provides an alternative proof of Theorem 1.

Proof of Theorem 2

Define by induction the sequence of polynomials (f n ) by

• f 0 = 1 • If i ∈ A ∞ , i ≥ 1 f i (x) = x 0 f i-1 (y)dy • If i / ∈ A ∞ , i ≥ 1 f i (x) = 1 x f i-1 (y)dy
We claim that

d N = N ! 1 0 f N -1 (x)dx
Indeed, Theorem 1 tells us that the right-hand side is the number of permutations with length N and descent set N -D N . Using the symmetries of the problem, we see that this is equal to the number of permutations with length N and descent set D N .

Consider the bivariate generating function

F (x, y) = ∞ n=0 y n f n (x)
Then we have

∞ n=1 d n n! y n = 1 0 yF (x, y)dx
So it is sufficient to compute F . Let µ be the number of ascents in a period, µ = p -|A|. By differentiating,

∂ p F (x, y) ∂x p = (-1) µ y p F (x, y)
Solutions of this differential equation have the general form

F (x, y) = p k=1 a k (y)e ω k yx
where the ω k are the p-th roots of 1 if µ is even and the p-th roots of -1 if µ is odd, and where the a k are power series:

a k (y) = n≥0 a k,n y n Since f 0 = 1, p k=1 a k,0 = 1
Next, we can re-express

f n (x) = n l=0 p k=1 ω l k x l l! a k,n-l
and so

f ′ n+1 (x) = n l=0 p k=1 ω l+1 k x l l! a k,n-l If n ≥ 1 is a descent, f n (0) = 0 whence p k=1 a k,n = 0 If n ≥ 1 is an ascent, f n (0) = 1 whence p k=1 a k,n = - n j=1 p k=1 a k,n-j ω j k j!
Moreover, f ′ n+1 = f n or -f n according as whether n is a descent or an ascent. Therefore, for all n, j ≥ 0, identifying the coefficient of x j in the polynomials f n+j and f ′ n+j+1 , we find that if n + j + 1 is a descent,

p k=1 a k,n ω j k = p k=1 a k,n ω j+1 k and if n + j + 1 is an ascent, p k=1 a k,n ω j k = - p k=1 a k,n ω j+1 k
This leads to the following matrix expression. Put

B n =    a 1,n . . . a p,n   
Let V be the Vandermonde matrix

V =      1 . . . 1 ω 1 . . . ω p . . . . . . ω p-1 1 . . .ω p-1 p      (7)
For each n ≥ 0 put

U n =    u n,1 0. . .0 . . . . . . . . . u n,p 0. . .0    where u n,1 = 1 (8)
and for every j ≥ 1, if n + j is a descent, u n,j+1 = u n,j while if n + j is an ascent, u n,j+1 = -u n,j . In other words,

u n,j+1 = u n,j × (-1) 1 {n+j / ∈A∞ } (9)
Next, put ρ 0 = 1 and for n ≥ 1, if n is a descent, put ρ n = 0 and if n is an ascent,

ρ n = - n j=1 p k=1 a k,n-j ω j k j!
Then we have the matrix equality

V B n = U n      ρ n 0 . . . 0      That is, B n = V -1 U n      ρ n 0 . . . 0     
Let us now express the generating series of ρ n . For every

i ∈ [1, p], ∞ n=1 ρ n t n u n,i = - p k=1 ∞ n=1 n-1 j=0 a k,j ω n-j k (n -j)! u n,i t n 1 {n / ∈A∞} = - p k=1 ∞ j=0 ∞ n=j+1 a k,j ω n-j k (n -j)! u n,i t n 1 {n / ∈A∞} = - p k=1 ∞ j=0 ∞ m=1 a k,j ω m k m! u m+j,i t m+j 1 {m+j / ∈A∞} = - p k=1 p-1 l=0 ∞ n=0 a k,np+l t np+l ∞ m=1 (ω k t) m m! u m+l,i 1 {m+l / ∈A∞} Thus let us put, for k ∈ [1, p] h k,l (t) = ∞ n=0 a k,np+l t np+l and for k ∈ [1, p], l ∈ [0, p -1] and i ∈ [1, p], f k,l,i (t) = ∞ n=1 (ω k t) n n! u n+l,i 1 {n+l /

∈A∞}

Then we can re-express, for every

i ∈ [1, p], ∞ n=1 ρ n t n u n,i = - p k=1 p-1 l=0 h k,l (t)f k,l,i (t)
In matrix terms, this gives

∞ n=1 t n U n      ρ n 0 . . . 0      =    p k=1 p-1 l=0 h k,l (t)f k,l,1 (t) . . . p k=1 p-1 l=0 h k,l (t)f k,l,p (t)    whence ∞ n=1 t n B n =    p k=1 h 1,l (t) -a 1,0 . . . p k=1 h p,l (t) -a p,0    = -V -1    p k=1 p-1 l=0 h k,l (t)f k,l,1 (t) . . . p k=1 p-1 l=0 h k,l (t)f k,l,p (t)   
We want to deduce an expression of the functions h k,l from the equality above.

To do so, put, for k ∈

[1, p], l ∈ [0, p -1], m ∈ [1, p] and i ∈ [1, p], g k,l,m,i (t) = ∞ n=0 (ω k t) np+m (np + m)! u m+l,i 1 {m+l / ∈A∞} Then for l ∈ [0, p -1],    h 1,l (t) -a 1,0 1 {l=0} . . . h p,l (t) -a p,0 1 {l=0}    = -V -1    p k ′ =1 p-1 l ′ =0 h k ′ ,l ′ (t)g k ′ ,l ′ ,l-l ′ ,1 (t) . . . p k ′ =1 p-1 l ′ =0 h k ′ ,l ′ (t)g k ′ ,l ′ ,l-l ′ ,p (t)   
where the index ll ′ in the right-hand side is to be taken modulo p. So let H(t) be the p 2 -dimensional vector with entries in C [[t]] defined by

H(t) =                     h 1,0 (t) . . . h p,0 (t) h 1,1 (t) . . . h p,1 (t) . . . h 1,p-1 (t) . . . h p,p-1 (t)                     (10)
Likewise, let B be the p 2 -dimensional vector with entries in C defined by

B =           a 1,0 . . . a p,0 0 . . . 0           (11)
where the first p entries are given by

   a 1,0 . . . a p,0    = V -1    u 0,1 . . . u 0,p    (12) For l ∈ [0, p -1], let G l be the p × p 2 matrix with entries in C[[t]] defined as follows. For k ∈ [1, p], l ′ ∈ [0, p -1], k ′ ∈ [1, p], put (G l ) k,l ′ p+k ′ = g k ′ ,l ′ ,l-l ′ ,k (t) Then for l ∈ [0, p -1],    h 1,l (t) . . . h p,l (t)    = -V -1 G l H(t)
Finally, let G be the p 2 × p 2 matrix with entries in C [[t]] defined by

G =    G 0 . . . G p    In other words, if l > l ′ , G lp+k,l ′ p+k ′ = ∞ n=0 (ω k ′ t) np+l-l ′ (np + l -l ′ )! u l,k 1 {l / ∈A∞} (13) whereas, if l ≤ l ′ , G lp+k,l ′ p+k ′ = ∞ n=1 (ω k ′ t) np+l-l ′ (np + l -l ′ )! u l,k 1 {l / ∈A∞} (14)
Remark that G has valuation 1. Let V be the p 2 × p 2 matrix, V being diagonal by blocks and given by

V = -      V -1 0 . . . 0 0 V -1 . . . 0 . . . . . . . . . 0 . . . 0 V -1      (15)
Then

H(t) -B = VGH(t) whence H(t) = (I -VG) -1 B ( 16
)
Remark that each entry of G is either 0 or has the form ±g k,l (t) with the functions g k,l (t) defined in Theorem 2. Moreover, using the formula for the inversion of matrices, we check that the entries of V have a rational expression, with integer coefficients, in terms of the roots ω 1 , . . . , ω p . Therefore, H itself has a rational expression, with integer coefficients, in terms of the ω k and the g k,l (t), and so have all the functions h k,l (t). Finally,

∞ n=1 d n n! t n = 1 0 tF (x, t)dx = p k=1 p l=1 h k,l (t) e ω k t -1 ω k (17)
and each term e ω k t -1 ω k also has a rational expression, with integer coefficients, in terms of the ω k and the g k,l (t). This establishes the existence of the rational function R A as stated in Theorem 2.

Note that, according to ( 13) and ( 14), if l is a descent, all the lines of the matrix G between l + 1 and l + p are 0. So we may suppress these lines and we are lead to a system of size µp × µp instead of p 2 × p 2 .

Example. Let us examine the case when when µ = 1. Suppose that the only ascent is at position j ∈ [1, p]. First, the ω k are roots of -1: for 1 ≤ k ≤ p, ω k = exp (2k + 1)iπ p

Next, it turns out that h k,l (t) = 0 whenever l = j, whence

∞ n=1 d n n! t n = p k=1 h k,j (t) e ω k t -1 ω k
In that case, the vector H from (10) becomes p-dimensional:

H(t) =    h 1,j (t) . . . h p,j (t)   
To compute H, define the p-dimensional vector U by

U =           1 . . . 1 -1 . . . -1          
where the first j entries are 1 and the last pj entries are -1. The matrix G(t) becomes a p × p matrix, since we can restrict equations ( 13) and ( 14) to the case l = l ′ = j. Using the definition of the u l,m , we find that G(t) is a matrix where all the lines are equal, namely

G k,k ′ (t) = ∞ n=1 (ω k ′ t) np (np)!
and ( 16) becomes

H(t) = (I + V -1 G(t)) -1 V -1 U
In the particular case of alternating permutations,

ω 1 = i, ω 2 = -i, G(t) = [cos(t) -1] 1 1 1 1 V -1 = 1 2 1 -i 1 i and U = 1 1
Diagonalizing the matrix V -1 G(t), we find

-V -1 G(t) = 1 -cos(t) 2 1 -i 1 -i 1 + i 1 + i = 1 -cos(t) 2
P -1 2 0 0 0 P for some invertible matrix P . By geometric summation,

(I + V -1 G(t)) -1 = ∞ n=0 (-V -1 G(t)) n = P -1 1/ cos(t) 0 0 0 P and therefore, (I + V -1 G(t)) -1 = 1 2 cos(t) 1 -i 1 -i 1 + i 1 + i
The generating series can then be expressed as

∞ n=1 d n n! t n = (1 -i)(e it -1) 2i cos(t) + (1 + i)(e -it -1) -2i cos(t)
Adding the term d 0 = 1, we recover André's result.

Implementation

We have implemented our algorithm in Maple 1 . The length of the permutation is nmax. The random variables U [n], which are drawn independently at random, represent the descent profile: U [n] = 0 or 1 according as whether nmaxn is a descent or an ascent. The reals V [n] in the code correspond to the random variables Y i in the algorithm. As pointed out in the introduction, generating the variables Y i amounts to solving N equations of the form (4). Of course, we can only get approximate solutions and since the Y i are computed recursively, there is a propagation of errors. We shall not discuss this point from a theoretical point of view here. Observe, however, that these errors do not affect the permutation we generate as long as, for each i, the difference between Y i and its approximation is smaller than the minimal value of 2|Y i -Y j | over all i, j ∈ [1, nmax].

Since equations of the form (4) are polynomial and since f ′ i = ±f i-1 , we can directly implement Newton's method in the algorithm. We stop the iterations in Newton's method as soon as two successive approximations of the solution are at distance less that 10 -p , where p is a parameter that can be tuned. We give the time needed to generate permutations of variable length and for various values of p. For a permutation of length 2000, we have run the algorithm for the values p = 4 and p = 5 and we observe empirically that this yields the same permutation.

Here are some tables giving the run time of the algorithm for permutations of various lengths. In the first table, we have also let the parameter p vary. The first line is the length of the permutation, the second line the run time (in seconds) for p=4 and the third line the run time for p=5. Empirically, the run time is more than quadratic. This is due to the fact that we have to perform computations with polynomials of large degree. Also, the quantity of memory is quite large if we want to record all the coefficients of all the polynomials. Of course, only keeping in memory the polynomial of highest degree and the descent profile would be sufficient, since the other polynomials can be obtained by differentiation. However, this would mean computing twice these polynomials.

  us prove Corollary 1. It suffices to prove that if B, B are two subsets of [2, N -1] with B = B ∪ {m}, m / ∈ B, then the number of permutations with set of extrema B is smaller than the number of permutations with set of extrema B. Using the symmetries of the problem, it suffices to count permutations with set of extrema B, B ending with a descent. If we impose this condition, let A, Ã be the corresponding descent sets. Using A, Ã we define by descending induction the polynomials f i , fi as in Theorem 1. Now define the polynomials

  .936 2.212 4.468 6.757 11.360 16.145 23.366 31.734 37.930 0.196 0.972 2.521 4.656 7.360 12.293 17.705 24.662 34.594 41.151 In the next simulations, we have taken p=4 and we only give the integer part of the run time (in seconds).

	100	200	300	400	500	600		700	800	900	1000
	0.184 01200 1400 1600 1800 2000 2200 2400 2600 2800 3000
	68	90	119	174	213	269	351	627	678	766
	3300 3600 3900 4200 4500 4800 5100 5400 5700 6000
	840	1210 1410 1822 2385 3283 5395 5877 6380 9121

The code is available at http://www.math.univ-paris13.fr/∼marchal/algopermut.m

Acknowledgments I am very grateful to Cyril Banderier for his help in the implementation of the algorithm. I thank Olivier Bodini and various participants of the conference "Permutation patterns 2013" for useful comments and references.