Validation of Formal Specification: the Case for Animation

Atif Mashkoor, Faqing Yang, and Jean-Pierre Jacquot LORIA -Nancy Université Vandoeuvre-lès-Nancy, France {firstname.lastname}@loria.fr Formal methods such as B [START_REF] Abrial | The B Book[END_REF] or Event-B [START_REF] Abrial | Modeling in Event-B: System and Software Engineering[END_REF] are designed around the idea that a piece of code can be "correct per construction." They use the usual notion of correctness: the program is a mathematically proven implementation of the specification. They are good candidates for industrial use for two major reasons: they embody a development process, refinement, which breaks the notoriously difficult correctness proof into many small and manageable proof obligations, and they have effective tool support.

Through the realization and the analysis of two important case-studies using Event-B in the area of transportation [START_REF] Mashkoor | B événementiel pour la modélisation du domaine: application au transport[END_REF][START_REF] Mashkoor | Domain Engineering with Event-B: Some Lessons We Learned[END_REF][START_REF] Yang | Prouvé ? Et après ? In: 10th Approches Formelles dans l'Assistance au Développement de Logiciels (AFADL'10)[END_REF][START_REF] Yang | Scaling up with Event-B: A Case Study[END_REF], we have discovered that proofs alone are not sufficient to produce a "good" software: we need also to execute, i.e., to test, the software. Two main reasons justify this proposition. Some properties, notably temporal, are virtually impossible to model with constructs such as states, invariants, or events which Event-B provides us. We need to ensure that the specification is an adequate model of the problem we want to solve. Software must be verified and validated.

Waiting to have an executable program to begin the validation leads to the same difficulties as proving a program against its specification: costly, very complex, soon unmanageable. The strategy which works with the verification of specifications built by a stepwise refinement process could also work for validation. We should be able to enrich the development process by adding validation activities to each refinement step. Such an enhancement is discussed in [START_REF] Yang | Prouvé ? Et après ? In: 10th Approches Formelles dans l'Assistance au Développement de Logiciels (AFADL'10)[END_REF].

The idea of executing specifications is not new [START_REF] Balzer | Operational specification as the basis for rapid prototyping[END_REF]. Some tools have already been developed and integrated into Rodin, ProB [START_REF] Leuschel | ProB: An Automated Analysis Toolset for the B Method[END_REF] or Brama1 , for instance. However, these tools are often unable to animate well written abstract specifications. In fact, qualities sought after for a well written specification, such as abstractness, non determinism, non constructive definitions, are contradictory with what is required for effective computation.

To use validation early in the development cycle, which in our opinion is highly desirable, we need to find techniques to adapt the existing tools. The cost of the technique should be minimal so that it can be frequently used during refinement process.

The first approach to the problem is based on the observation that being provable and being animatable are distinct qualities for specifications: they can be one, the other, or both. It is then possible, in many cases, to "downgrade" a proven, non animatable, specification into a "behaviorally equivalent" animatable, but non provable, specification. In [START_REF] Mashkoor | Transformation Heuristics for Formal Requirements Validation by Animation[END_REF], we have proposed several transformations to realize this idea. The validity of such a technique depends on the semantics of the transformations. Since ours are not semantic preserving in a strict sense, we had to develop an ad-hoc semantics based on the behavior of the model. The preservation property we want to show can be summed up as: "whatever we observe on the execution of the transformed specification would have been observed on the proven specification." We have obtained some results in this direction.

Despite having transformation rules at our disposal, sometimes animators still fail to execute a specification. For such cases, we expect that translation of Event-B text into an executable programming language can be a reasonable fall-back strategy. We are currently studying this on a 2D platooning model [START_REF] Yang | Scaling up with Event-B: A Case Study[END_REF] which is out of reach of both ProB and Brama. Current experiments with hand-written systematic translation are promising.

http://www.brama.fr


Jérémy Malaizé
email: jeremy.malaize@newport-fr.com

Jean Lévine
email: jean.levine@ensmp.fr

Active Estimation of the Initial Phase for Brushless Synchronous Motors

This paper addresses the initial phase estimation problem for brushless synchronous motors. Only displacement measurements are used (no current) and friction, load and motor parameters are supposed to be unknown. Because of friction, the system is modelled by a differential equation with discontinuous right-hand side. Specific open-loop inputs are designed (active method) to get the initial phase as a function of the magnitude of the displacements along the corresponding trajectories. The estimation relies on a complete classification of the possible dynamical behaviours of the considered discontinuous right-hand side system with periodic input, whatever values the unknown parameters may take. We propose an approximated formula of the initial phase. Some experimental results are given, together with a comparison of our method to a classical procedure.

I. INTRODUCTION

We consider both rotary and linear synchronous motors used in high precision positioning applications, such as the semiconductor wafer industry. Brushless motors are preferred to DC motors with brushs for their longer life time, improved cleanliness of operation and better resistance to wear and tear. For DC motors (with brushs), current commutations, as well as initialization, are performed mechanically while, for brushless motors (without mechanical contact), additionnal measurements and an algorithm are required for the currents and the magnetic field to be in phase. The initialization procedure aims at getting the currents in phase with the position of the motor coils inside the magnetic field. In this paper, we propose a method using only displacement measurements to determine the initial rotor position (modulo the magnetic pitch) for brushless synchronous motors 1 .

When it comes to reducing the number of sensors for economical reasons, current sensors are usually preferred to position encoders. Several papers address sensorless control of synchronous motors, that is to say the determination of the initial rotor position from current measurements. Some existing techniques consist in using an observer to estimate the back-electromotive force induced in the coils, see [START_REF] Dhaouadi | Design And Implementation Of An Extended Kalman Filter For The State Estimation Of A Permanent Magnet Synchronous Motor[END_REF] and [START_REF] Lim | Observers For Sensorless Control Of Permanent Magnet Synchronous Motor Drive[END_REF] for example. Other approaches determine the motor windings inductance by injecting either specific carrier-frequency signals ( [START_REF] Ostlund | Sensorless Rotor-Position Detection From Zero To Rated Speed For An Integrated PM Synchronous Motor Drive[END_REF], [START_REF] Schmidt | Initial Rotor Angle Detection Of A Non-Salient Pole Permanent Magnet Synchronous Machine[END_REF] or [START_REF] Noguchi | Initial Rotor Position Estimation Method Of Sensorless PM Motor With No Sensitivity To Armature Resistance[END_REF]) or no signal at all ( [START_REF] Kulkarni | A Novel Position Sensor Elimination Technique For The Interior Permanent Magnet Synchronous Motor Drive[END_REF] or [START_REF] Ogasawara | An Approach To Real Time Position Estimation At Zero And Low Speed For A PM Motor Based On Saliency[END_REF]). Both these techniques rely on current measurements to determine an electrical variable in phase with the magnetic field. 1 US patent pending.

For the aforementionned applications featuring nanometric positioning precision, high resolution position sensors turn out to be necessary. Optical resolvers, or absolute position sensors, are too expensive and, in this paper, only optical incremental position sensors measuring the relative displacement around the initial position are considered. We propose a method that provides an autonomous current controller (whose current measurements are not made available to the initialization algorithm) with specific inputs computed offline and simply gets the corresponding displacement measurements back to estimate the initial rotor position. This estimation scheme, illustrated figure 1, spares the use of additional and unnecessary Hall effect sensors. An existing solution 2 , falling under the scope of the previously defined scheme, consists in maintaining a constant current in an electrical phase until the mechanical equilibrium. The corresponding position shift is then used to compute an estimate of the initial phase. Unfortunately, before stabilizing, erratic oscillations around the equilibrium position are witnessed. These movements are as large as the magnetic pitch (typically a couple of millimeters) and cannot be avoided. This behaviour might not be suitable for the aforementionned high precision applications.

In this paper, the method can be tuned to generate arbitrarily small magnitude displacements (a couple of microns). According to scheme 1, neither magnetic field nor currents have to be fed to the algorithm and the motor parameters are unknown (gain, load and friction). In these conditions, our method features an accurate initial phase estimation with little computation, which makes real-time implementation easier.

This paper is organized as follows. We first model the dynamics of the brushless synchronous motors and highlight the need for a precise initial position estimation : guaranteeing stable and accurate positioning together with an optimal use of the motor, from an energetic point of view (section II). We then derive the relation between the initial phase, the other unknown parameters and the magnitude of the displacements when the system is driven by a periodic and open-loop input (section III). This analysis is based on a complete classification of the periodic orbits of the considered system under forced oscillations. In particular, it is shown how the friction, modelled by a function of the sign of the speed (thus giving rise to a discontinuous right-hand side system in the sense of Filippov (see [START_REF] Filippov | Differential Equations with Discontinuous Righthand Sides[END_REF])), affects the trajectories followed by the system. Then, we compare the measured displacements to those predicted by the classification, and, thanks to an approximation, the initial position is obtained. In section IV, experimental results are provided to consolidate the different assumptions and show that our initialization procedure outperforms the classical method.

II. INITIAL PHASE ESTIMATION

A. Modelling

For simplicity's sake, only two-phase linear motors (depicted on figure 2) are allowed for since the generalization is rather obvious. Referring to figure 2, a track made up of permanent magnets creates a sinusoidal magnetic field whose magnitude and spatial period are respectively noted B 0 (unknown) and P (given). Let l be the length of the active windings so that a force F 1 is created by the electrical phase #1 :

F 1 = i 1 l B 0 sin 2π P x + (-i 1 ) l B 0 sin 2π P (x -P 2 ) = K m i 1 sin 2π P x = K m i 1 sin 2π P (x 0 + d) , (1)
with K m = 2 l B 0 the motor gain, x 0 the initial position and d the measured displacement. By calculations similar to (1), the force created by the second electrical phase is derived :

F 2 = K m i 2 cos 2 π P (x 0 + d) . (2)
Let's assume a controller stabilizes the currents i 1 and i 2 so that their dynamics turn out to be fast and negligible compared to the mechanical behaviour of the system. The whole system is thus modelled by :

ẍ = K m m i 1 sin 2π P (x 0 + d) + K m m i 2 cos 2π P (x 0 + d) - f m sign( ẋ). (3)
The friction force f , the gain K m and the load m are unknown.

To drive (3) when x 0 is unknown, let's assume x 0 , an estimate of x 0 , is available (keep in mind the determination of x 0 motivates this paper) and make the following choice for the currents i 1 and i 2 :

i 1 = I sin 2π P ( x 0 + d) + ϕ i 2 = I cos 2π P ( x 0 + d) + ϕ , (4)
with I a current reference and ϕ a supplementary degree of freedom. Both I and ϕ may be time-varying functions or not. Equation ( 3) together with the choice (4) finally writes :

ẍ = K m m cos(ϕ 0 -ϕ 0 -ϕ) I - f m sign( ẋ),
with ϕ 0 = 2 π P x 0 the initial phase and ϕ 0 = 2 π P x 0 . Let K m and m be a priori estimates of K m and m, respectively. Since Km m I is homogenous to an acceleration, from now on, let's note :

K m m I = ẍref , (5)
which, together with α = Km m × m

Km

, leads to the following model of the system :

ẍ = α cos(ϕ 0 -ϕ 0 -ϕ) ẍref - f m sign( ẋ). (6)

B. Role of the initial phase

From now on, our concern is to get an estimate of the initial phase from (6) by appropriately choosing ẍref and ϕ and despite the lack of knowledge on the other parameters.

Let ϕ = 0, then the force created by the motor is :

F = F 1 + F 2 = K m cos(ϕ 0 -ϕ 0 ) I,
and, consequently :

|F | ≤ |K m I| = m α |ẍ ref | .
Thus, for a given acceleration reference ẍref , as soon as the relation ϕ 0 = ϕ 0 holds, the norm of the force is lower than the desired value. The energetic transfer between the current and the force is optimal when ϕ 0 is accurately determined. Moreover, the gain of the open loop from x to ẍref depends on cos(ϕ 0 -ϕ 0 ). Therefore, the accuracy of the estimation is also material to guarantee both the precision and the stability of the positioning.

III. ESTIMATION PROCEDURE A. Choice of the trajectories

We now propose a specific input ẍref that both generates displacements as small as possible and differentiates between the contributions of the different unknown parameters of the model [START_REF] Guckenheimer | Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields[END_REF].

Note that, under the following assumptions :

α = 1, ϕ 0 = 0, f = 0,
if ẍref is the second derivative of an oscillating trajectory of magnitude ξ 1 -ξ 0 , the magnitude of the measured displacements is immediately given by cos (ϕ 0 -ϕ) (ξ 1 -ξ 0 ). Choosing ϕ = 0 and ϕ 0 = π 2 successively gives an estimation of cos (ϕ 0 ) and sin (ϕ 0 ). This way, the small magnitude constraint is fulfilled by appropriately tuning ξ 1 -ξ 0 and the initial phase is easily obtained.

Let's design such an oscillating trajectory by first defining an elementary trajectory from an initial point ξ 0 , at rest, to a final point ξ 1 , also at rest, in a time T > 0 :

1) Let's first impose the following initial conditions :

x elem (0) = ξ 0 , ẋelem (0) = 0, ẍelem (0) = 0.

2) Similarly, the final conditions are :

x elem (T ) = ξ 1 , ẋelem (T ) = 0, ẍelem (T ) = 0.
3) A possible choice for x elem (t) is the polynomial interpolation of degree 5 given by :

x elem (t) = ξ 0 + (ξ 1 -ξ 0 ) 5 i=1 a i t T i . (7)
4) The coefficients {a i } 1≤i≤5 are derived by identifying [START_REF] Kulkarni | A Novel Position Sensor Elimination Technique For The Interior Permanent Magnet Synchronous Motor Drive[END_REF] at t = 0 and t = T with the initial and final conditions, to eventually get :

a 1 = a 2 = 0, a 3 = 10, a 4 = -15, a 5 = 6.
The reverse trajectory from ξ 1 , at equilibrium, to ξ 0 , at equilibirum, in a time T > 0 is similarly derived. Finally, let x M be the trajectory made up of M round trips from ξ 0 to ξ 1 and back to ξ 0 , using interpolation [START_REF] Kulkarni | A Novel Position Sensor Elimination Technique For The Interior Permanent Magnet Synchronous Motor Drive[END_REF].

Because of the unknown parameters in [START_REF] Guckenheimer | Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields[END_REF], when ẍref = ẍM , the system does not follow the trajectory x M (keep in mind this is an open-loop identification procedure) unless the following conditions are fullfilled :

α = 1, ϕ = ϕ 0 , ϕ = 0, f = 0,
that is to say when all the parameters are perfectly known. As long as f = 0, only the magnitude of the displacements is affected by the unknown parameters and the measured trajectory is proportional to x M . On the other hand, when f = 0, the friction introduces a delay.

This brief study suggests that, when ẍref = ẍM , it is possible to ditsinguish between the different unknown parameters as they delay or rescale the reference trajectory x M . We now propose a complete analysis for in-depth explanation of how to extract the initial phase from the measured displacements despite unknown parameters.

B. Dynamical behaviours classification

The idea underlying the proposed method is to integrate equation ( 6) twice, with ẍref = ẍM , to explicit the relation between the measured displacements around the initial position and the different parameters of the model featuring a discontinuous right-hand side. As addressed in [START_REF] Filippov | Differential Equations with Discontinuous Righthand Sides[END_REF], in the vicinity of the discontinuity surface, whether the system crosses the surface or slides on it directly depends on the model parameters. This link is pointed out hereafter by introducing the reduced parameter µ, whose definition is given by :

µ = α m| cos(ϕ 0 -ϕ 0 -ϕ)|ẍ max /f, (8)
with ẍmax the maximal value of the reference acceleration ẍM :

ẍmax = max 0≤t≤T |ẍ elem (t)| . (9)
As stated in theorem 1, the value µ leads to a complete classification of the dynamical behaviours of the system (6) driven by ẍM . This theorem makes up an essential prerequisite to integrate (6) twice.

Theorem 1: For each µ > 1, there exists only one attractive periodic orbit for the system (6) initially at rest with ẍref = ẍM .

There also exist two real numbers µ 1 and µ 2 , 1 < µ 1 ≤ µ 2 , depending only on ẍM , that define three and only three possible behaviours :

(i) 1 < µ ≤ µ 1 ≈ 1.4 : the system instantaneously reaches a periodic orbit along which sliding on the surface of discontinuity ẋ = 0 occurs (see simulation results 3 figure 3). (ii) µ 1 < µ ≤ µ 2 ≈ 1.7 : the system reaches a periodic orbit in finite time (lower than 2 T ) along which sliding on the surface of discontinuity ẋ = 0 occurs (see simulation results figure 4). (iii) µ > µ 2 : the system asymptotically reaches a periodic orbit and no sliding occurs along this orbit (see simulation results figure 5). Moreover, if µ ≤ 1, the system remains at rest. x(t) ×10 -3 Fig. 3. State space (left) and temporal (right) trajectories for 1 < µ ≤ µ 1 . 3 The following reduced system is simulated : ẍ = ẍM ẍmax -µ -1 sign( ẋ).

C. Integration of (6)

According to theorem 1, whatever value µ may take, the system reaches a periodic orbit for sure and it is possible to distinguish between orbits sliding on the surface of discontinuity (µ ≤ µ 2 ) or not (µ > µ 2 ). This classification makes the integration of ( 6) tremendously easier and let δ be the magnitude of the displacements on the orbit :

δ = lim k→∞    max 0≤τ ≤T x(k T + τ ) -x(k T )    . ( 10
)
δ is naturally related to µ and thus to the initial phase. The theorem 2, that, for simplicity reasons, we give without proving it, states this relation :

Theorem 2: For the system (6) with ẍref = ẍM , there exists a function ∆ defined on R + and taking values in [0 , 1[ so that :

δ = α |cos(ϕ 0 -ϕ 0 -ϕ)| (x 1 -x 0 )∆(µ). (11)
Moreover, ∆ (illustrated figure 6(a)) depends only on ẍM .

D. Estimation

From now on we impose ϕ 0 = 0 in ( 6) and ( 11), and we propose a way of extracting ϕ 0 from the relation [START_REF] Ostlund | Sensorless Rotor-Position Detection From Zero To Rated Speed For An Integrated PM Synchronous Motor Drive[END_REF].

Let {ϕ i } 1≤i≤N be a set of N real numbers so that, according to [START_REF] Ostlund | Sensorless Rotor-Position Detection From Zero To Rated Speed For An Integrated PM Synchronous Motor Drive[END_REF], N values δ i are defined :

∀i , δ i = ε i cos(ϕ 0 -ϕ i )α (x 1 -x 0 )∆(µ i ), (12)
with

µ 0 = α ẍmax m/f ε i = sign (cos(ϕ 0 -ϕ i )) µ i = ε i µ 0 cos(ϕ 0 -ϕ i ). (13)
Let (i , j) ∈ {1 , . . . , N } 2 and i = j, according to the previous notations, the following expression is derived :

δ i δ j = µ i ∆ (µ i ) µ j ∆ (µ j ) .
Let's note J ij (µ 0 , ϕ 0 ) the function defined by :

J ij (ϕ 0 , µ 0 ) =   δ i µ j ∆(µ j ) -δ j µ i ∆(µ i )   2 . ( 14
)
Keep in mind that µ i , whose definition is given by ( 13), is a function of both ϕ 0 and µ 0 .

Estimating the initial phase comes down to solving the following optimization problem :

  ϕ 0 , µ 0   = arg min ϕ0∈R µ0∈R +            i ,j j>i J ij (ϕ 0 , µ 0 )            (15)

E. Approximated solution

Even if different iterative methods could be used to solve (15), thanks to an approximation, a less computation demanding solution is derived. The relevancy of the approach is illustrated by experimental results in the next section.

The basic idea developped in this section is to obtain an optimization problem close to (15) but whose solution does not rely on a numerical algorithm. To this end, suppose there exists a constant γ ∈ R + so that µ ∆(µ) is linearized :

∀µ ≥ 1, µ ∆(µ) ≈ γ (µ -1). ( 16
)
Keep in mind the approximation ( 16), illustrated figure 6(b) in dotted line, is only valid for µ ≥ 1.

The elementary criterion J ij (ϕ 0 , µ 0 ), given by ( 14), together with (16), is rewritten in J * ij (ϕ 0 , µ 0 ) defined by :

J * ij (ϕ 0 , µ 0 ) = γ 2   δ i (µ j -1) -δ j (µ i -1)   2 = γ 2    (µi µ j ) A ij µ i µ j + b T ij µ i µ j + c 2 ij     (17)
with the following notations :

-

A ij = δ 2 j -δ i δ j -δ i δ j δ 2 i , -b T ij = -2 (δ i -δ j ) (-δ j δ i ), -c ij = (δ i -δ j )
2 . The criterion J(ϕ 0 , µ 0 ) to be minimized in (15) is approximated by J * (ϕ 0 , µ 0 ), a quadratic function of the vector µ = (µ 1 . . . µ N )

T :

J * (ϕ 0 , µ 0 ) = γ 2   µ T A µ + b T µ + c   . (18)
A is a N × N matrix, b an N -dimensional vector and c a scalar. They are defined by :

-

∀ (i , j) , i = j , A(i , i) = j =i δ 2 j j > i , A(i , j) = A(j , i) = -δ i δ j -∀ i , b(i) = 2 j =i (δ j (δ i -δ j )) -c = i ,j c ij
Note that the value of the unknown parameter γ can be ommitted to compute the minimum of J * .

To make the most of J * being a linear function of µ, a coordinates change is necessary. According to the definitions (13), for all i, µ i is rewritten in cartesian coordinates :

µ i = µ 0 ε i cos(ϕ i ) ε i sin(ϕ i ) cos (ϕ 0 ) sin (ϕ 0 ) .
Then, we define θ ∈ R 2 and M an N × 2 matrix :

θ = µ 0 cos (ϕ 0 ) sin (ϕ 0 ) M =    ε 1 cos(ϕ 1 ) ε 1 sin(ϕ 1 ) . . . . . . ε N cos(ϕ N ) ε N sin(ϕ N )   
so that µ = M θ is a linear function of the new optimization variables. Note that ϕ 0 is the phase of the vector θ and µ 0 its norm. J * is a function of the optimization variable θ with :

J * (ϕ 0 (θ) , µ 0 (θ)) = J(θ) = θ T M T AM θ + b T M θ + c.
Since approximation ( 16) is valid for µ ≥ 1, the minimization of J(θ) is subjected to N constraints :

ε i cos(ϕ i ) sin(ϕ i ) θ ≥ 1.
Let θ be the solution of this quadratic optimization problem under constraints :

       θ = arg min θ∈R 2 J(θ) Subject to : ∀i, ε i cos(ϕ i ) sin(ϕ i ) θ ≥ 1 (19)
θ is easily obtained by writting the optimality conditions of Kuhn and Tucker (see [START_REF] Culioli | Introduction à l'Optimisation[END_REF] for details) and the phase of the vector θ is ϕ 0 , the estimate of the initial phase.

IV. EXPERIMENTAL RESULTS AND COMPARISON

A. Implementation

To solve the approximated optimization problem (19), for each ϕ i , δ i and ε i need to be determined from the available measurements.

M is chosen large enough for the periodic orbit to be reached, and, to get rid of transient phenomena, an integer n 0 is picked up so that δ, defined by [START_REF] Ostlund | Sensorless Rotor-Position Detection From Zero To Rated Speed For An Integrated PM Synchronous Motor Drive[END_REF], is experimentally determined by :

δ ≈ 1 2 M -n 0 + 1 2 M k=n0 max 0≤τ ≤T |x(k T + τ ) -x(k T )|. (20)
Applying relation (20) for each ϕ i , we build up the set {δ i } 1≤i≤N .

Let's now consider the system (6) (with ϕ 0 = 0 and ϕ = ϕ i ) and recall this system is initially at rest. If ξ 1 > ξ 0 (resp. ξ 1 < ξ 0 ), then x(T ) is larger (resp. smaller) than x(0) if cos(ϕ 0 -ϕ i ) > 0 and, conversely, x(T ) is smaller (resp. larger) than x(0) if cos(ϕ 0 -ϕ i ) < 0. Thus, ε i is easily derived from the following relation :

ε i = sign ξ 1 -ξ 0 x(T ) -x(0) . (21)

B. Classical method

We now give a brief description of another initialization method only requiring displacement measurements. For this method, in phase #1, the current is constant and non-zero and, in phase #2, the current is zero :

i 1 = m K m ẍ0 , i 2 = 0.
ẍ0 is constant and homogenous to an acceleration (similarly to ( 5)). The evolution of the system is then given by :

ẍ = α sin 2 π P d + ϕ 0 ẍ0 - f m sign( ẋ).
The corresponding force of the motor generates damped oscillations (because of friction) before an equilibrium position, noted x ∞ , is reached. The corresponding displacement is d ∞ :

ϕ 0 = 2 π P x ∞ - 2 π P d ∞ mod 2 π.
If friction is not too significant, the following relation holds :

∃k ∞ ∈ Z, x ∞ ≈ x ∞ = P 2 + k ∞ P.
With no additionnal information, a possible estimate for ϕ 0 is :

ϕ 0 = 2 π P x ∞ - 2 π P d ∞ mod 2 π = π - 2 π P d ∞ mod 2 π. (22)
Actually, when the friction f is significant, the previous estimate is highly biased and it can be shown that the estimation error is related to the parameter µ ′ = α m ẍ0 /f . The smaller µ ′ , the worse the estimation.

C. Results

We illustrate the results of our method in comparison with the classical procedure. To this end, let's consider two linear motors whose parameters are unknown but the friction is known to be more important for the second motor, even though this information is not used by the initialization algorithms.

According to the previous notations, the first motor is defined by α 1 , m 1 , f 1 and the second by α 2 , m 2 , f 2 .

To compare the methods, we impose µ 0 = µ ′ (µ 0 is defined by ( 13)) which turns out to be equivalent to ẍ0 = ẍmax (ẍ max is given by ( 9)). We estimate the average estimation error ϕ 0 -ϕ 0 when ϕ 0 is varied from 0 to 2 π.

Let's consider the first linear motor to experimentally check out that the precision of the classical method worsens as µ ′ = α 1 m 1 ẍ0 /f 1 decreases while our solution does not depend upon µ 0 = α 1 m 1 ẍmax /f 1 . To this end, we consider the following experimental conditions :

-ẍ0 = ẍmax = 1000 mm/s 2 (figure 7 As expected, the method of the paper offers the same accuracy when µ 0 is varied. For this motor, the precision on the initial phase is better than 10 • , which confers an efficiency of 98% to the motor. For high accelerations (figure 7(a)), the classical method gives slightly less accurate results but generates erratic and large magnitude displacements. As expected, when µ ′ decreases (figure 7(b)), this method exhibits a precision of 30 • and a corresponding efficiency of 87% while our procedure still offers the same 98% efficiency. For figure 7 and the motor with little friction, the precision of our method is constant and independent of the motor parameters (µ 0 ) and gives a substantially better energetic efficiency to the motor.

Let's now move on to the second motor which exhibits important friction so that the reference acceleration has to be larger (ẍ max = ẍ0 = 4000 mm/s 2 ). The average estimation error is given for both the methods on figure 8. Our method is very robust to dry friction since it still achieves the same precision (approximately 10 • ). On the contrary, the classical does not perform that well.

The method of the paper turns out to be independent of µ 0 which means insensibility to friction, load and motor gain. Whatever values these parameters may take, for the considered motors, our method outperforms the classical procedure and a precision of 10 • is achieved.

V. CONCLUSION

A new method to determine the initial phase for synchronous motors has been proposed. No information about the motor internal functioning is required and the algorithm only has to be fed with relative position measurements from incremental coders. For high precision applications, no additionnal hardware is needed which allows to get rid of Hall effect sensors. Our solution exhibits better results than the classical method though implemented in the same conditions. Despite friction, poorly identified motor gain or unknown motor load, our method features precise initial phase estimation which guarantees optimal energetic use of synchronous motors.

Fig. 1 .

 1 Fig. 1. Estimation Scheme.

Fig. 2 .

 2 Fig.2. Two-phase ironless linear motor.

  x

Fig. 4 .Fig. 5 .

 45 Fig. 4. State space (left) and temporal (right) trajectories for µ 1 < µ ≤ µ 2 .

Fig. 6 .

 6 Fig. 6. Functions ∆(µ) on the left and µ ∆(µ) on the right (solid line).

  (a)), -ẍ0 = ẍmax = 500 mm/s 2 (figure 7(b)). ẍ0 = ẍmax = 500 mm/s 2

Fig. 7 .

 7 Fig.7. Average estimation error ϕ 0 -ϕ 0 versus initial phase ϕ 0 . Our method (19), solid line, and the classical method (22), dotted line.

Fig. 8 .

 8 Fig.8. Average estimation error ϕ 0 -ϕ 0 versus initial phase ϕ 0 . Our method (19), line, and the classical method (22), dotted line.


Thomas Devos
email: thomas.devos@ensmp.fr

Jean Lévine
email: jean.levine@ensmp.fr

A FLATNESS-BASED NONLINEAR PREDICTIVE APPROACH FOR CRANE CONTROL

Keywords: nonlinear predictive control, differential flatness, motion planning, crane control

. The controller is composed of two parts: the first one is a traditional PD output feedback to track the reference trajectory and reject small perturbations, the second one consists of updating the reference trajectory from the current estimated state of the crane to the desired equilibrium point on a receding horizon each time the pursuit error exceeds a given threshold. Simulations are presented to illustrate its performances.

During operations, the crane is supposed to move as fast as possible to carry the load from an initial position to a desired final destination, avoiding obstacles and sway. Controlling the position of the load requires both motion planning and feedback control to attenuate the perturbations. To this aim, the mechanical degrees of freedom of the crane are actuated, all but those corresponding to the load. In addition, sensors measuring the respective actuators' positions together with two synchronized digital cameras giving the three coordinates of the load via image processing, at a a slower but sufficiently high rate, are available. Note that the whole state of the crane model is not directly available through these measurements and must be estimated.

According to the nature of this weight handling system, the presence of large perturbations such as those created by frictions of the motors or belt elasticity (a belt being used to transmit the motor rotation to the platform) or unknown load mass, and limitations on the motors, make the tuning tradeoff between performance and robustness delicate: to stabilize the load at a rest point, small gains are preferred to keep a good observability of the load oscillations through the motor sensors, whereas to follow fast trajectories, large gains are required. Moreover, if large errors appear, the actuators are soon saturated.

To circumvent this difficulty, we propose a predictive control approach (for recent surveys on this topic see e.g. [START_REF] Morari | Model predictive control: past, present and future[END_REF][START_REF] Qin | An overview of nonlinear model predictive control applications[END_REF]): in addition to a first PD output feedback loop for reference trajectory tracking, each time the pursuit error exceeds a given threshold, the reference trajectory is updated from the (estimated) present state to the desired rest point on a receding horizon. Thus, since we start anew with a reference trajectory at the present estimated state, the tracking error is instantaneously reduced to the estimation error. With this approach, we expect that the rate of decrease of the tracking error is much faster than with pure feedback and certainty equivalence. This nonlinear predictive control method makes an extensive use of the differential flatness property of the system (see other related linear and nonlinear flatness-based approaches in [START_REF] Fliess | Continuoustime linear predictive control and flatness : a module-theoretic setting with examples[END_REF][START_REF] Delaleau | Commande prédictive non linéaire fondée sur la platitude différentielle[END_REF]) to generate such reference trajectories and stabilize the load along them. It also depends on the observer design used to estimate the tracking error.

The remaining part of the paper is organized as follows. The next section gives a description of our crane reduced size model. In section 3, we present the crane model in three dimensions. The differential flatness of the model and its consequences on motion planning are studied in section 4. We describe in section 5 the basic PD controller used as a first step towards stabilization and perturbation rejection. In section 6, we design an observer estimating the successive derivatives of the flat output and the performances of the observer are studied in simulations. Then the overall predictive control algorithm is presented and simulations showing its efficiency to reduce the oscillations of the load are given in section 7.

DESCRIPTION OF THE REDUCED SIZE MODEL OF THE CRANE

The considered model is the reproduction of a US Navy crane at the scale 1/80 (Fig. 1 and2). It has been created by a designer1 in 1998 in a brass structure. It is composed of 4 DC motors for the 4 following displacements.

• A motor for the platform rotation (motor 4).

• A motor for the vertical motion of the load (motor 2). • A motor for the vertical motion of the free mobile pulley (motor 3). • And next, a motor for the horizontal motion of the load via the translation of the free mobile pulley (motor 1).

L 2 T 1 ,m 1 T 2 ,m 2 L 1 ( x A C B P O 0 ,y 0 , z 0 ) ( x ,y , z) C 4 ,M L 3 ( ( x x 1 2 , , y y 1 2 , , z z 1 2 ) ) m r Fig. 2.

US Navy crane in 3 dimensions

Sensors of the four motor positions allow to determine the length of the ropes via a correct initialization. Some important coordinates of the crane are used, such as the positions of the fixed pulleys (the coordinates of the pulleys located at A and P are (x 1 , y 1 , z 1 ) and (x 2 , y 2 , z 2 )), the coordinates of the load is (x, y, z) (point C) and the position of the mobile pulley is (x 0 , y 0 , z 0 ) (point B). Moreover, the three fixed pulleys are lined up and we have

-→ OA = α 1 --→ OP .

CRANE MODEL

The model of the crane is obtained by Newton's second principle applied to the various rigid bodies constituting the crane (load, motors, free pulley ...). In the model used in this article, the mass of the free pulley is neglected (no dynamics associated to the free pulley). Under this assumption, the load, the free pulley and the ropes BA, BP and CB are in the same plane (Fig. 3(b)). This particular plane's dynamics are described by the angles ξ and ϕ (Fig. 3(a)), and then, the positions of the ropes and of the load in this plane by the angles γ and β and by the L 3 . More details are given in [START_REF] Kiss | Planification de trajectoires et commande d'une classe de systèmes mécaniques plats et liouvilliens[END_REF]).

First, we introduce three frames. The origin of the first one, denoted by K b , is fixed at O and its zaxis coincides with the vertical rotation axis of the crane. The two other frames have also their origin at O. The second one, denoted by K g is chosen such that the z b -axis of the frame K b and the point P determine its x g z g -plane. These frames are transformed into one another by the rotation of angle ξ around z b . The last frame, denoted by K is chosen such that the points P , A, and B determine its xz-plane. The transformation between K g and K is a combination of three rotations of angle α, ϕ and -α (α is the angle between the jib of the crane and the vertical rotation axis of the crane).

The coordinates of the point C in the frame K are thus obtained as functions of the angles γ and β and the length L 3 (see Fig. 3(b)). Denoting by x K , y K and z K the coordinates of the load in the frame K, we have the following (non differential) equations:

       x K = k sin α + L 1 sin(γ + (α -β)) +L 3 sin(2γ + (α -β)) y K = 0 z K = k cos α + L 1 cos(γ + (α -β)) +L 3 cos(2γ + (α -β)) (1)
with L 1 function of β and γ:

L 1 = l sin β sin γ . ( 2
)
Note that we also have

L 2 = l sin(γ -β) sin γ . ( 3
)
Denoting by x K b , y K b and z K b the coordinates of the point C in the frame K b , the second principle applied to the load reads:

m ẍK b ÿK b zK b + g = Ω K b K g • Ω K g K -T 3 sin θ 0 T 3 cos θ (4)
with θ the angle beetween --→ BC and the z-axis of the frame K, and T 3 the rope tension modulus at C. Finally, the torque balance at A, P and O is given by

J 1 ρ 1 L1 = T 1 ρ 1 -u 1 J 2 ρ 2 L = T 2 ρ 2 -u 2 J pf ξ = proj z b T 2 --→ OP × --→ P B --→ P B + T 1 -→ OA × --→ AB --→ AB +u 4 (5)
with L = L 2 + L 3 , T 1 and T 2 the respective rope tension moduli at the points A and P , u 1 and u 2 the torques delivered by the motors related to the pulleys at A and P , u 4 the torque applied to the platform and J 1 , J 2 and J pf the respective inertias of the motors reported to the points A, P and 0.

We can eliminate all the variables but ξ, ϕ, γ, β, L 3 and their derivatives by combining these equations , to obtain a model of the form A(X, u) Ẋ = b(X, u) (see [START_REF] Kiss | Planification de trajectoires et commande d'une classe de systèmes mécaniques plats et liouvilliens[END_REF] for more details). We then invert A(X, u), to obtain the (nonlinear) explicit form Ẋ = f (X, u). The measured signal is the vector Y = (ξ, L 1 , L, x K b , y K b , z K b ). The state vector X of dimension 10, is made up with the angles ξ, ϕ, γ and β, the length L 3 and their first derivatives.

DIFFERENTIAL FLATNESS

The US Navy crane, as for more general cranes, is a differentially flat system. A flat output is remarkably the position of the load ( [START_REF] Kiss | Planification de trajectoires et commande d'une classe de systèmes mécaniques plats et liouvilliens[END_REF]). So, if we know the position of the load and its derivatives up to the fourth order, we are able to compute the open loop controls that generate the desired trajectory. The proof can be found in [START_REF] Kiss | Planification de trajectoires et commande d'une classe de systèmes mécaniques plats et liouvilliens[END_REF].

First, the flatness property is used to generate a rest-to-rest trajectory corresponding an idle-to-idle displacement of the load. The trajectory, allowing to avoid obstacles, is obtained using polynomial interpolation. In the present case, since we have 10 conditions (5 at the initial point: x i , ẋi , ẍi , x

(3) i , x (4) i and 5 at the end:

x f , ẋf , ẍf , x (3) f , x (4)
f ), we construct a 9-th degree polynomial for the three coordinates of the load, which reads, for the x coordinate of the load:

x(t) = x i +(x f -x i ) t -t i t f -t i 5 4 j=0 a j t -t i t f -t i j
(6) with a 0 = 126, a 1 = -420, a 2 = 540, a 3 = -315 and a 4 = 70.

The same formula is obtained for y(t) and z(t).

From these trajectories, all the system variables can be deduced. In particular, the reference controls of u 1 , . . . u 4 that exactly generate the above trajectories, assuming the exactness of the model and without perturbations, can be deduced ( [START_REF] Fliess | Flatness and defect of nonlinear systems: introductory theory and examples[END_REF][START_REF] Fliess | A Lie-Bäcklund approach to equivalence and flatness of nonlinear systems[END_REF][START_REF] Kiss | Modelling, flatness and simulation of a class of cranes[END_REF]Kiss et al., 2000a)). The corresponding lengthy formulas are omitted.

PD FEEDBACK

We know from [START_REF] Kiss | Planification de trajectoires et commande d'une classe de systèmes mécaniques plats et liouvilliens[END_REF]) that around any rest-to-rest reference trajectory (L 1ref , L ref , ξ ref ) and (u 1ref , u 2ref , u 4ref ) of for (L 1 , L, ξ) and (u 1 , u 2 , u 4 ) respectively, obtained from the previous section, the system is locally stabilized by the following PD controller:

u 1 = u 1ref + K p1 (L 1 -L 1ref ) + K d1 ( L1 -L1ref ) u 2 = u 2ref + K p2 (L -L ref ) + K d2 ( L -Lref ) u 4 = u 4ref + K p4 (ξ -ξ ref ) + K d4 ( ξ -ξref ) (7)
with the gains K p1 , K d1 , K p2 , K d2 , K p4 and K d4 all positive.

Moreover, this controller makes the end point a globally stable equilibrium.

As announced in the introduction, the tuning of these gains may be quite delicate, depending on the perturbations (frictions, belt elasticity, . . . ), the parametric errors (on masses, inertias, . . . ) and the time constants of the displacements. Therefore, to alleviate the duty of this controller, we propose to update the reference trajectory each time the pursuit error exceeds a given bound, the gains of the controller (7) being fixed once for all. However, since the current state from which the new reference trajectory must start is not measured, we need to construct an observer to deduce it from the available measurements. This is the purpose of the next section.

OBSERVER

The rope lengths L 1 and L = L 2 + L 3 , the angle ξ and the coordinates of the load (x, y and z) are measured by two synchronized digital cameras but their successive derivatives have to be estimated.

We propose an extended Kalman filter to estimate the vector (x, ẋ, ẍ, x (3) , x (4) , y, . . . , y (4) , z, . . . , z (4) ), which, according to section 4, is necessary to compute all the system state.

Recall that we have denoted by X = (ξ, ϕ, γ, β, L 3 , ξ, φ, γ, β, L3 ) the state vector of the system and that the crane is represented by:

Ẋ = f (X, u) Y = h(X) (8)
with the mesurements Y = (ξ, L 1 , L, x, y, z).

Its linear tangent approximation along the above trajectory is easily deduced and gives the matrices

A lin (t) = ∂f ∂X (X ref (t)), u ref (t)), B lin (t) = ∂f ∂u (X ref (t)), u ref (t)) and C lin (t) = ∂h ∂X (X ref (t)). The observer has the form δ X = A lin (t)δ X + B lin (t)δu + K(Y -Ŷ ) Ŷ = h( X) (9) with δ X = X -X ref .
We then calculate the asymptotic gain K of the Kalman filter as a function of the dynamical and observation noise covariance matrices, that are tuned to obtain a fast enough convergence.

Once the state vector X of the system is estimated, we obtain the successive derivatives of the flat output by the change of coordinates that maps X = ( ξ, φ, γ, β, L3 , ξ, φ, γ, β, L3 ) to (x, . . . , x(4) , ŷ, . . . , ŷ(4) , ẑ, . . . , ẑ(4) ).

The formulas of this change of coordinates are omitted for simplicity's sake.

To illustrate the results obtained with the above extended Kalman filter, we introduce an output white noise of variance 10 -6 for the rope lengths, the angle ξ and the coordinates of the load, which corresponds to an average linear deviation of 1 mm and angular deviation of approx. 0.03 • , the order of magnitude of the noises on our reduced size setup.

The results are presented in Fig. 4. The observer shows a satisfactory convergence up the 3rd order derivative, and the convergence of the 4th order is still acceptable.

PREDICTIVE CONTROL

Define the pursuit error by: and let σ be a given positive real number corresponding to the threshold on the error e above which the reference trajectory must be updated.

e = (x -x ref ) 2 + (ŷ -y ref ) 2 (10)
The predictive control algorithm is the following: let t be the first time after the initial time t i such that e(t) = σ. From X(t), computed online, we deduce x(t), . . . , x(4) (t), ŷ(t), . . . , ŷ(4) (t) (z is not considered here since the perturbations acting on it are negligible). Therefore, we can update the reference trajectory starting from the current estimated state at t and finishing at rest at a given time t + τ , with τ , the receding horizon, to be precised later. Note that the formula ( 6) is no longer valid since we don't start from a rest point anymore. The new formula is given by:

x ref (t) = x(t) + 9 j=1 b j t -t τ j (11)
with the b j 's appropriately computed to arrive at rest at t + τ .

The receding horizon τ is chosen according to the following heuristic rule (which can be verified on simulations): the closed-loop input maximal amplitudes remain comparable to the ones of the open-loop references if the horizon τ is of the same order of magnitude as the inverse of the lowest closed-loop eigenfrequency. According to simulations, the frequencies of the response of the closed-loop system with (7) to a perturbation are shown on Fig. 5. Therefore, we choose τ corre-sponding to an average period. We can verify that it allows to attain the rest point of the load as fast as possible thanks to the reference update (11), without saturating the different motors. Note that this problem is all the more crucial that we want to generate fast motions.

Once the updated trajectory is computed, we check if some inequality constraints on the position, velocity and acceleration of the load are satisfied all along. In case of an affirmative answer, the trajectory y is declared admissible and the corresponding updated reference control is generated. Oherwise, τ is increased.

The same steps apply each time the error happens to increase up to σ.

If the initial state of the trajectory is well estimated, the real motion of the load in response to the updated control fits with the desired trajectory and no oscillations or only small ones remain after t + τ . An external perturbation is applied to the load at the beginning, at t = 4 s and t = 8 s, inducing oscillations of an amplitude of about 4 cm on the x and y axis. Piecewise constant voltage perturbations are also applied to the motors. Finally, a real load mass 25 % higher than the one considered in the model is used. The receding horizon of the updated trajectory is τ = 0.8 s. The coordinates of the load and the corresponding motor voltages are presented in the subfigures (d), (e) and (f).

It may be seen that the load oscillations and the perturbations on the motors are satisfactorily attenuated thanks to the updated trajectories and the PD corrector.

CONCLUSION

A flatness-based predictive control to reduce the load oscillations of a reduced size crane has been presented. Simulation results show that this method works as expected. This work is a first step towards its implementation on the crane and further studies, in particular concerning the robustness of this approach, will be done in the near future.

Fig. 1 .

 1 Fig. 1. Reduced size model of the US Navy crane Cranes are used in various industries such as construction and naval transport. Any automatic control improvement of the assistance to crane

Fig. 3 .

 3 Fig. 3. Geometry of the crane

  Fig. 4. simulations: Estimated derivatives of the coordinate x with an output noise; (a): x and x; (b): ẋ and ẋ; (c): ẍ and ẍ; (d): x (3) and x(3) ; (e): x (4) and x(4) .

Fig. 5 .

 5 Fig. 5. simulations: Fourier transform of the signals x(t) and y(t) in the frame K g for different z: z = 0, 0.25, 0.5, 0.75, 0.1; (a): oscillation frequencies of x K g ; (b): oscillation frequencies of y K g ; Fig. 6 illustrates our predictive control performances. In this simulation, the coordinates of the load (x, y, z) are sampled at the frequency of 10 Hz, whereas the frequency of the sampled measurements sent by the motor sensors (L, L 1 and the angle ξ) and is about 200 Hz. Moreover, the measurements are corrupted by a bounded noise in the interval [-10 -3 , 10 -3 ].

  Fig. 6. simulations: Predictive control; (a) Coordinate x real and reference ; (b) Coordinate y real and reference ; (c) Coordinate z real and reference ; (d) Voltage corresponding to the controlled torque u 1 ; (e) Voltage corresponding to the controlled torque u 2 ; (f ) Voltage corresponding to the controlled torque u 4 ;

Walter Rumsey : Artenciel,

24, rue le Regrattier, 75004 Paris, FRANCE. Tel : 01 43 25 73 40, Fax : 01 43 25 73 40


Felix Antritter
email: felix.antritter@unibw.de

Regelungstechnik

Jean Lévine
email: jean.levine@ensmp.fr

Towards a Computer Algebraic Algorithm for Flat Output Determination

Keywords: I, 1, 4 [Computing Methodologies]: Symbolic and Algebraic Manipulation-Applications; G, 4 [Mathematics of Computing]: Mathematical Software

), that yields an abstract algorithm for its computation. In this paper it is discussed how these conditions can be checked using computer algebra. All steps of the algorithm are discussed for the simple (but rich enough) example of a non holonomic car.

INTRODUCTION

We consider a nonlinear control system ẋ = f (x, u)

where x = (x1, . . . , xn) is the state vector, u = (u1, . . . , um) the control vector, m ≤ n, and f is a meromorphic function of its arguments. We say that this system is differentially flat, or shortly flat ( [START_REF] Martin | Contribution à l'étude des systèmes diffèrentiellement plats[END_REF][START_REF] Fliess | Sur les systèmes non linéaires différentiellement plats[END_REF]), if and only if there exists a vector y = (y1, . . . , ym) such that: (i) y and its successive time derivatives ẏ, ÿ, . . . are functionnally independent, (ii) y is a function of x, u and a finite number of time derivatives of the components of u, (iii) x and u can be expressed as functions of the components of y and a finite number of their successive time derivatives:

x = ϕ(y, ẏ, . . . , y (α) ), u = ψ(y, ẏ, . . . , y (α+1) ) for some multi-integer α = (α1, . . . , αm), and with the notation y (α) = ( d α 1 y 1 dt α 1 , . . . , d αm ym dt αm ). A vector y having these properties is called a flat output. This concept has inspired an important literature and a large number of practical and industrial applications (see e.g. [START_REF] Martin | Flat systems[END_REF] for a survey). Its main advantages rely on the simplicity to solve the motion planning and stable tracking problems. Various formalisms have been introduced to study this remarkable class of systems: finite dimensional differential geometric approaches ([2, 10, 28, 29]), differential algebra and related approaches ( [START_REF] Fliess | Flatness and defect of nonlinear systems: introductory theory and examples[END_REF][START_REF] Aranda-Bricaire | A linear algebraic framework for dynamic feedback linearization[END_REF][START_REF] Jakubczyk | Invariants of dynamic feedback and free systems[END_REF]), infinite dimensional differential geometry of jets and prolongations ([9, 30, 22, 20, 24]). Among these contributions, the characterization of differential flatness takes a large part ( [START_REF] Aranda-Bricaire | A linear algebraic framework for dynamic feedback linearization[END_REF][START_REF] Charlet | Sufficient conditions for dynamic state feedback linearization[END_REF][START_REF] Chetverikov | New flatness conditions for control systems[END_REF][START_REF] Franch | Flatness, Tangent Systems and Flat Outputs[END_REF][START_REF] Jakubczyk | Invariants of dynamic feedback and free systems[END_REF][START_REF] Martin | Systems without drift and flatness[END_REF][START_REF] Pereira Da Silva | Flatness of nonlinear control systems : a Cartan-Kähler approach[END_REF][START_REF] Pomet | On dynamic feedback linearization of four-dimensional affine control systems with two inputs[END_REF][START_REF] Rathinam | Configuration flatness of Lagrangian systems underactuated by one control[END_REF][START_REF] Rouchon | Necessary condition and genericity of dynamic feedback linearization[END_REF][START_REF] Shadwick | Absolute equivalence and dynamic feedback linearization[END_REF][START_REF] Sluis | A necessary condition for dynamic feedback linearization[END_REF][START_REF] Lévine | On flatness necessary and sufficient conditions[END_REF][START_REF] Lévine | On necessary and sufficient conditions for differential flatness[END_REF]). We follow here the results of [START_REF] Lévine | On flatness necessary and sufficient conditions[END_REF][START_REF] Lévine | On necessary and sufficient conditions for differential flatness[END_REF] in the formalism of manifolds of jets of infinite order ( [START_REF] Fliess | A Lie-Bäcklund approach to equivalence and flatness of nonlinear systems[END_REF][START_REF] Krasil'shchik | Geometry of Jet Spaces and Nonlinear Partial Differential Equations[END_REF][START_REF] Pomet | A differential geometric setting for dynamic equivalence and dynamic linearization[END_REF][START_REF] Zharinov | Geometrical Aspect of Partial Differential Equations[END_REF]). For the stated flatness conditions implicit systems are considered which are obtained from (1) by eliminating the input vector u. We recall the notions of Lie-Bäcklund equivalence and Lie-Bäcklund isomorphism in this context and the flatness necessary and sufficient conditions in terms of polynomial matrices and differential forms. Note that this approach may be seen as an extension to nonlinear systems of [START_REF] Lévine | Flat output characterization for linear systems using polynomial matrices[END_REF] and provide flatness conditions that are invariant by endogeneous dynamic feedback extension. The derived conditions use differential operators which combine differential geometric concepts like exterior derivative and wedge product as well as algebraic concepts as operations on skew polynomials with coefficients that are meromorphic functions of the coordinates. Existing computer algebra systems offer lots of functionalities for each of the mentioned fields but their combination is not considered. In this paper we show how to implement such operators in Maple 11, and include them in an algorithm to check the flatness necessary and sufficient conditions. Note that this algorithm doesn't necessarily finish in a finite number of steps. The paper is organized as follows: Section 2 is devoted to the basic description of implicit control systems on manifolds of jets of infinite order. The notions of Lie-Bäcklund equivalence and Lie-Bäcklund isomorphism in this framework are recalled and the flatness property is defined precisely. In Section 3 the necessary and sufficient conditions for differential flatness derived in [START_REF] Lévine | On flatness necessary and sufficient conditions[END_REF][START_REF] Lévine | On necessary and sufficient conditions for differential flatness[END_REF] are stated and Section 4 discusses in more detail the introduced operators and conditions. General formulae for the operators are derived from its defining relations and it is explained how the operators and the algorithm developed in Section 3 can be implemented using a computer algebra system like Maple 11. Finally, in Section 5 the application of the implemented method is illustrated for the well-known nonholonomic car example. This example, though simple, may however call for all the resources of the algorithm. A first draft of this algorithm can be obtained by sending an e-mail to the first author.

IMPLICIT CONTROL SYSTEMS ON MAN-IFOLDS OF JETS OF INFINITE ORDER

Given an infinitely differentiable manifold X of dimension n, we denote its tangent space at x ∈ X by TxX, and its tangent bundle by TX. Let F be a meromorphic function from TX to R n-m . We consider an underdetermined implicit system of the form F (x, ẋ) = 0 (2) regular in the sense that rk ∂F ∂ ẋ = n -m in a suitable dense open subset of TX. According to the implicit function theorem, any explicit system (1) with x ∈ X, (x, f (x, u)) ∈ TxX for every u in an open subset U of R m , and rk ∂f ∂u = m in a suitable open subset of X × U , can be locally transformed into (2), and conversely. A vector field f that depends, for every x ∈ X, on m independent variables u ∈ R m in a meromorphic way with rk ∂f ∂u = m in a suitable open subset of X × R m , satisfying F (x, f (x, u)) = 0 for every u ∈ U , is called compatible with [START_REF] Charlet | Sufficient conditions for dynamic state feedback linearization[END_REF]. Note that this elimination step, though easy for some classes of systems, e.g., affine with respect to u, may be non trivial in general. Remark also that the implicit representation (2), as opposed to [START_REF] Aranda-Bricaire | A linear algebraic framework for dynamic feedback linearization[END_REF], is invariant by endogeneous dynamic extension (see [START_REF] Fliess | A Lie-Bäcklund approach to equivalence and flatness of nonlinear systems[END_REF] for a precise definition). In [START_REF] Fliess | A Lie-Bäcklund approach to equivalence and flatness of nonlinear systems[END_REF] (see also [START_REF] Pomet | A differential geometric setting for dynamic equivalence and dynamic linearization[END_REF] where a similar approach has been developed independently), infinite systems of coordinates (x, u) = (x, u, u, . . .) have been introduced to deal with prolonged vector fields

f (x, u) = n i=1 fi(x, u) ∂ ∂xi + m j=1 k≥0 u (k+1) j ∂ ∂u (k) j
the original system being in explicit form [START_REF] Aranda-Bricaire | A linear algebraic framework for dynamic feedback linearization[END_REF]. Following [START_REF] Lévine | On flatness necessary and sufficient conditions[END_REF][START_REF] Lévine | On necessary and sufficient conditions for differential flatness[END_REF], we adopt an external description of the prolonged manifold containing the solutions of (2): we consider the infinite dimensional manifold X defined by

X def = X × R n ∞ def = X × R n × R n × .
. ., made of an infinite (but countable) number of copies of R n , endowed with the product topology, and we assume that we are given the global infinite set of coordinates of X:

x = x, ẋ, . . . , x (k) , . . . , .

Recall that, in this topology, a function ϕ from X to R is continuous (resp. differentiable) if ϕ depends only on a finite (but otherwise arbitrary) number of variables and is continuous (resp. differentiable) with respect to these variables. C ∞ or analytic or meromorphic functions from X to R are then defined as in the usual finite dimensional case since they only depend on a finite number of variables. We endow X with the so-called trivial Cartan field ( [START_REF] Krasil'shchik | Geometry of Jet Spaces and Nonlinear Partial Differential Equations[END_REF][START_REF] Zharinov | Geometrical Aspect of Partial Differential Equations[END_REF])

τ X = n i=1 j≥0 x (j+1) i ∂ ∂x (j) i . (4)
We also denote by

Lτ X ϕ = n i=1 j≥0 x (j+1) i ∂ϕ ∂x (j) i = dϕ dt
the Lie derivative of a differentiable function ϕ along τ X and L k τ X ϕ its kth iterate. Thus x

(k) i = d k x i dt k = L k τ X xi for every i = 1, . . . , n and k ≥ 1, with the convention x (0) i = xi. Since d dt x (j) i def = ẋ(j) i = x (j+1) i
, the Cartan field acts on coordinates as a shift to the right. X is thus called manifold of jets of infinite order. From now on, x y, . . . stand for the sequences of jets of infinite order of x, y,. . .

A regular implicit control system is defined as a triple (X, τ X , F ) with X = X × R n ∞ , τ X its associated trivial Cartan field, and

F meromorphic from TX to R n-m ) satisfying rk ∂F ∂ ẋ = n -m in a suitable open subset of TX.

Lie-Bäcklund equivalence for implicit systems

We recall from [START_REF] Lévine | On flatness necessary and sufficient conditions[END_REF][START_REF] Lévine | On necessary and sufficient conditions for differential flatness[END_REF] the following definition: Let us consider two regular implicit control systems (X, τ X , F ), with

X = X × R n ∞ , dim X = n and rk ∂F ∂ ẋ = n -m, and (Y, τ Y , G), with Y = Y × R p ∞ , dim Y = p, τ Y its trivial Cartan field, and rk ∂G ∂ ẏ = p -q. Set X0 = {x ∈ X|L k τ X F (x) = 0, ∀k ≥ 0} and Y0 = {y ∈ Y|L k τ Y G(y) = 0, ∀k ≥ 0}.
They are endowed with the topologies and differentiable structures induced by X and Y respectively.

Definition 1. We say that the regular implicit control systems (X, τ X , F ) and (Y, τ Y , G) are Lie-Bäcklund equivalent (or shortly L-B equivalent) at the pair of points (x0, y 0 ) ∈ X0 × Y0 if and only if (i) there exist neighborhoods X0 and Y0 of x0 in X0 and y 0 in Y0 respectively and a one-to-one meromorphic mapping Φ = (ϕ0, ϕ1, . . .) from Y0 to X0 satisfying Φ(y 0 ) = x0 and such that the trivial Cartan fields are Φ-related, namely Φ * τ Y = τ X ;

(ii) there exists Ψ one-to-one and meromorphic from X0 to Y0, with Ψ = (ψ0, ψ1, . . .), such that Ψ(x0) = y 0 and Ψ * τ X = τ Y . The mappings Φ and Ψ are called mutually inverse Lie-Bäcklund isomorphisms at (x0, y 0 ). The two systems (X, τ X , F ) and (Y, τ Y , G) are called locally L-B equivalent if they are L-B equivalent at every pair (x, Ψ(x)) = (Φ(y), y) of an open dense subset Z of X0 × Y0, with Φ and Ψ mutually inverse Lie-Bäcklund isomorphisms on Z. As a result, local L-B equivalence preserves equilibrium points, namely points ỹ (resp. x) such that G(ỹ, 0) = 0 (resp. F (x, 0) = 0), and coranks (m = q).

Differential Forms

Let us introduce a basis of the tangent space TxX of X at a point x ∈ X consisting of the set of vectors

{ ∂ ∂x (j) i |i = 1, . . . , n, j ≥ 0}.
A basis of the cotangent space T * x X at x is given by {dx

(j) i |i = 1, . . . , n, j ≥ 0} with < dx (j) i , ∂ ∂x (l) k >= δ i,k δ j,l , δ i,k being the Kronecker symbol.
The differential of F is thus given, in matrix notations, by

dF = ∂F ∂x dx + ∂F ∂ ẋ d ẋ. (5)
Note that the shift property of d dt on coordinates extends to differentials:

d dt dx = d ẋ = d d dt x, i.e. d
dt commutes with d. Since a smooth function depends on a finite number of variables, its differential contains only a finite number of non zero terms. Accordingly, we define a 1-form on X as a finite linear combination of the dx (j) i 's, with coefficients meromorphic from X to R or, equivalently as a local meromorphic section of T * X. The set of 1-forms is noted Λ 1 (X). We also denote by Λ p (X) the module of all the p-forms on X, by (Λ p (X)) m the space of all the m-dimensional vector p-forms on X, by (Λ(X)) m the space of all the m-dimensional vector forms of arbitrary degree on X, and by Lq

((Λ(X)) m ) = L (Λ p (X)) m , Λ p+q (X)
m , p ≥ 1 , the space of all linear operators from (Λ p (X)) m to Λ p+q (X) m for all p ≥ 1, where L (P, Q) denotes the set of linear mappings from a given space P to a given space Q.

Note that if Φ is a meromorphic mapping from Y to X, the definition of the (backward) image by Φ of a 1-form is the same as in the finite dimensional context.

Flatness

First recall from [START_REF] Fliess | A Lie-Bäcklund approach to equivalence and flatness of nonlinear systems[END_REF] that a system in explicit form is flat if and only if it is L-B equivalent to a trivial system. The reader may easily check that this definition is just a concise restatement of the definiton given in the Introduction. In our implicit context, it reads:

Definition 2. The implicit system (X, τ X , F ) is flat at (x0, y 0 ) ∈ X0 × R m ∞ if and only if it is L-B equivalent at (x0, y 0 ) to the trivial implicit system (R m ∞ , τ R m ∞ , 0).
In this case, the mutually inverse L-B isomorphisms Φ and Ψ are called inverse trivializations. The next result is proven in [START_REF] Lévine | On necessary and sufficient conditions for differential flatness[END_REF].

Theorem 1. The system (X, τ X , F ) is flat at (x0, y 0 ) ∈ X0 × R m
∞ if and only if there exists a local meromorphic invertible mapping Φ from R m ∞ to X0, with meromorphic inverse, satisfying Φ(y 0 ) = x0, and such that Φ * dF = 0. (6)

NECESSARY AND SUFFICIENT CONDI-TIONS FOR DIFFERENTIAL FLATNESS

We now analyze condition [START_REF] Fliess | A remark on Willems' trajectory characterization of linear controllability[END_REF] in more details: it characterizes the linear tangent mapping of Φ whose image entirely lies in the kernel of dF . The set of such mappings may be obtained in a systematic way in the framework of polynomial matrices by considering the following matrices polynomial with respect to the differential operator d dt (we use indifferently d dt for Lτ X or Lτ R m ∞ , the context being unambiguous):

P (F ) = ∂F ∂x + ∂F ∂ ẋ d dt , P (ϕ0) = j≥0 ∂ϕ0 ∂y (j) d j dt j (7)
with P (F ) (resp. P (ϕ0)) of size (n -m) × n (resp. n × m. Equation ( 6) reads: Φ * dF = P (F )P (ϕ0)dy = 0.

Clearly, the entries of the matrices in [START_REF] Fliess | Sur les systèmes non linéaires différentiellement plats[END_REF] are polynomials of the differential operator d dt with meromorphic coefficients from X to R. We denote by K the field of meromorphic functions from X to R and by K[ d dt ] the principal ideal ring of polynomials of d dt with coefficients in K.

Note that K[ d dt ] is non commutative, even if n = 1: for every a ∈ K, a = 0, we have d dt • x -x • d dt (a) = ẋa + x ȧ -x ȧ = ẋa = 0, or d dt • x -x • d dt = ẋ. For r, s ∈ N,
(F ) ∈ Mn-m,n[ d dt ]: V P (F )U = (∆, 0n-m,m) (9) with 0n-m,m the (n -m) × m matrix whose entries are all zeros, V ∈ Un-m[ d dt ], U ∈ Un[ d dt ] and ∆ ∈ Mn-m,n-m[ d dt
] a diagonal matrix whose entries di,i divide dj,j for all 0 ≤ i ≤ j ≤ n -m. Moreover, the degrees of the di,i's are uniquely defined (see [START_REF] Cohn | Free Rings and Their Relations[END_REF]).

Definition According to the equivalence between flatness and controllability of the tangent linear system (see [START_REF] Fliess | A Lie-Bäcklund approach to equivalence and flatness of nonlinear systems[END_REF]) and controllability and freeness of the module associated to the tangent linear system (see [START_REF] Fliess | A remark on Willems' trajectory characterization of linear controllability[END_REF]), it is proven in [START_REF] Lévine | On necessary and sufficient conditions for differential flatness[END_REF] that P (F ) is hyperregular around every integral curve of the system (X, τ X , F ) if and only if its corresponding tangent module is free.

Algebraic characterization of the differential of a trivialization

From now on, we assume that P (F ) is hyper-regular in a neighborhood of x0. In other words, there exist V and U such that V P (F )U = (In-m, 0n-m,m) .

U and V satisfying [START_REF] Franch | Flatness, Tangent Systems and Flat Outputs[END_REF] are indeed non unique. We say that

U ∈ R -Smith (P (F )) and V ∈ L -Smith (P (F )) if they are such that V P (F )U = (Im, 0). Accordingly, if M ∈ Mn,m[ d dt ] is hyper-regular with m ≤ n, we say that V ∈ L -Smith (M ) and W ∈ R -Smith (M ) if V ∈ Un[ d dt ] and W ∈ Um[ d dt ] satisfy V M W = Im 0 .
In place of (8), we first solve the matrix equation: P (F )Θ = 0 (11) where the entries of Θ ∈ Mn,m[ d dt ] are not supposed to be gradients of some function ϕ0.

Lemma 1. The set of hyper-regular matrices Θ ∈ Mn,m[ d dt ] satisfying [START_REF] Gohberg | Matrix Polynomials[END_REF] is nonempty and given by

Θ = U 0n-m,m Im W ( 12
)
with U ∈ R -Smith (P (F )) and W ∈ Um[ d dt ] arbitrary.

Lemma 2. For every Q ∈ L -Smith Û , with Û given by

Û = U 0n-m,m Im ( 13
)
there exists Z ∈ Um[ d dt ] such that QΘ = Im 0n-m,m Z. (14)
Moreover, for every

Q ∈ L -Smith Û , the submatrix Q = (0n-m,m, In-m) Q is equivalent to P (F ) (∃L ∈ Un-m[ d dt ] such that P (F ) = L Q).

Integrability

Let us denote by Qi,j

= k≥0 Q k i,j d k dt k the (i, j)-th poly- nomial entry of Q ∈ L -Smith Û obtained from Lemma 2.
We also denote by ω the m-dimensional vector 1-form defined by

ω(x) =    ω1(x) . . . ωm(x)    = (Im, 0m,n-m) Q(x)dx X 0 =      n j=1 k≥0 Q k 1,j (x)dx (k) j X 0 . . . n j=1 k≥0 Q k m,j (x)dx (k) j X 0      (15)
the restriction to X0 meaning that x ∈ X0 satisfies L k τ X F = 0 for all k and that the dx (k) j are such that dL k τ X F = 0 in X0 for all k. Since Q is hyper-regular, the forms ω1, . . . , ωm are independent by construction. Let us also recall that, if τ1, . . . , τm are given independent 1forms in Λ1 (X0), the K[ d dt ]-ideal T generated by τ1, . . . , τm is the set of all combinations with coefficients in K[ d dt ] of forms η ∧ τi with η an arbitrary form on X0 of arbitrary degree and i = 1, . . . , m.

Definition 4. We say that the K[ d dt ]-ideal T generated by τ1, . . . , τm is strongly closed if and only of there exists a matrix M ∈ Um[ d dt ] such that d(M τ ) = 0. This definition is indeed independent of the choice of generators.

Theorem 2. A necessary and sufficient condition for system (2) to be flat at the point (x0, y 0 ) is that there exist

U ∈ R -Smith (P (F )) and Q ∈ L -Smith Û , with Û
given by ( 13), such that the K[ d dt ]-ideal Ω generated by the 1-forms ω1, . . . , ωm defined by ( 15) is strongly closed in X0.

In order to develop the expression d(M τ ) for M polynomial matrix, we define the operator d by:

d (H) κ = d(Hκ) -Hdκ (16)
for all m-dimensional vector p-form κ in (Λ p (X)) m and all p ≥ 1. Note that (16) uniquely defines d (H) as an element of L1 ((Λ(X)) m ). We can prolong d for all µ ∈ Lq ((Λ(X)) m ) and for all κ ∈ (Λ p (X)) m and all p ≥ 1 by the formula:

d (µ) κ = d(µ κ) -(-1) q µ dκ. ( 17
)
Theorem 3. The K[ d dt ]-ideal Ω generated by the 1-forms ω1, . . . , ωm defined by ( 15) is strongly closed in X0 (or, equivalently, the system (X, τ X , F ) is flat) if and only if there exists µ ∈ L1 ((Λ(X)) m ), and a matrix

M ∈ Um[ d dt ] such that dω = µ ω, d (µ) = µ 2 , d (M ) = -M µ. ( 18
)
where we have noted µ 2 = µµ.

In addition, if [START_REF] Martin | Flat systems[END_REF] holds true, a flat output y is obtained by integration of dy = M ω. Note that conditions (18) may be seen as a generalization in the framework of manifolds of jets of infinite order of Cartan's well-known moving frame structure equations (see e.g. [START_REF] Chern | Lectures on Differential Geometry[END_REF]).

A Theoretical Algorithm

From the necessary and sufficient conditions ( 18), we derive the following abstract algorithm:

1. We first compute a Smith decomposition 1 of P (F ) and then Û as decribed in Lemmas 1 and 2. If P (F ) is not hyperregular, the system is non flat. Otherwise compute the vector 1-form ω defined by [START_REF] Lévine | On necessary and sufficient conditions for differential flatness[END_REF].

2. We compute the operator µ such that dω = µω by componentwise identification. It is easy to prove that such µ always exists.

3. Among the possible µ's, only those satisfying d (µ) = µ 2 are kept. If no µ satisfy this relation, the system is non flat.

4.

We then compute M such that d (M ) = -M µ, still by componentwise identification.

5. Finally, only those matrices M which are unimodular are kept. If there are no such M , the system is non flat. In the opposite case, a flat output is obtained by integration of dy = M ω, which is possible since d(M ω) = 0.

COMPUTER ALGEBRA IMPLEMENTA-TION

Before discussing the computer algebra implementation, useful general structures and formulae are derived. Note that computational cost aspects are not considered here. We only aim at showing that Algorithm 3.3 can be implemented using a standard computer algebra system.

The structure of elements of M r,s [ d dt ]

The elements of matrices A ∈ Mr,s[ d dt ] have the structure

[A] ij = k≥0 a ijk d k dt k , i = 1, 2, . . . r; j = 1, 2, . . . , s (19)
where the a ijk 's are smooth functions on X. Thus, for ω ∈ (Λ p (X)) s , Aω ∈ (Λ p (X)) r is the vector whose ith component is

[Aω]i = s j=1 k≥0 a ijk L k τ X ωj, i = 1, 2, . . . , r (20)
The multiplication of two matrices

A ∈ Mr 1 ,r 2 [ d dt ] and B ∈ Mr 2 ,r 3 [ d dt ] given by [A]i,j = k≥0 a i,j,k d k dt k , and [B]i,j = k≥0 b i,j,k d k dt k , is [AB]i,j = r 2 l=1 k 1 ,k 2 ≥0 a i,l,k 1 d k 1 dt k 1 b l,j,k 2 d k 2 dt k 2 = r 2 l=1 k 1 ,k 2 ≥0 k 1 k 3 =0 k1 k3 a i,l,k 1 L k 1 -k 3 τ X b l,j,k 2 d k 2 +k 3 dt k 2 +k 3 (21)
with

k1 k3 = k 1 ! k 3 !(k 1 -k 3 )! .
By explicitely indicating the dependence of the coordinates on the independent variable t it becomes possible to use the mult command (with DDt being used in our implementation to symbolize d dt ) of the DETools package of Maple 11. The function mskew (Input:

A ∈ Mr 1 ,r 2 [ d dt ], B ∈ Mr 2 ,r 3 [ d dt ]; Output: C = AB ∈ Mr 1 ,r 3 [ d dt ]
) provides a multiplication which corresponds to [START_REF] Pereira Da Silva | Flatness of nonlinear control systems : a Cartan-Kähler approach[END_REF]. Note furthermore that specifying the associated Cartan field to deal with functions of t only is not necessary (a discussion of this is done Section 4.4).

Smith-Decomposition of elements of M r,s [ d dt ]

With the matrix multiplication over K[ d dt ] using mskew, a Smith decomposition of matrices A ∈ Mr,s[ d dt ] can be implemented by adapting, e.g., the algorithm given in [START_REF] Gohberg | Matrix Polynomials[END_REF] for polynomial matrices with constant coefficients to the noncommutative case by constructing suitable unimodular matrices for left and right actions (see [START_REF] Cohn | Free Rings and Their Relations[END_REF] for more details). The resulting Maple procedure has been called Smith sa

(Input: A ∈ Mr,s[ d dt ]; Output: U ∈ R -Smith (A) and V ∈ L -Smith (A))
. The DETools package of Maple 11 provides all necessary operations.

4.3

The structure of elements of L q ((Λ(X)) m )

The elements of an operator µ ∈ Lq ((Λ(X)) m ) have the structure

[µ] ij = k≥0 µ ijk ∧ d k dt k , i, j = 1, 2, . . . , m (22)
where µ ijk is an arbitrary q-form, which means that, for every ω ∈ (Λ p (X)) m , µω ∈ Λ p+q (X) m is given by

[µω]i = m j=1 k≥0 µ ijk ∧ L k τ X ωj, i = 1, 2, . . . , m (23)
To implement operators µ ∈ Lq ((Λ(X)) m ) in Maple 11 we define the operator itself by specifying it of (Maple-) type Matrix, the components being differential forms on a previously specified Frame. A Frame fixes the coordinates of a jet space in the DifferentialGeometry package. The coefficients of the differential forms are then specified as polynomials in d dt . Thus the evaluation of µ on a p-form ω, according to [START_REF] Pomet | On dynamic feedback linearization of four-dimensional affine control systems with two inputs[END_REF], uses both ∧ and d dt . Therefore, a special function, called Dtwedge (Input:

ω ∈ (Λ p (X)) m , ω ∈ (Λ p (X)) m ; Out- put: µω ∈ Λ p+q (X)
m ) is introduced. The core algorithm of this function, for fixed i and j, results to (assuming that µ has some finite degree l w.r.t d dt )

• set k=0

• set γ = 0 while k ≤ l do • extract the coefficient µ ijk of µ ij (using the commmand DGmap of the DifferentialGeometry package together with the standard Maple command coeff )

• compute d k dt k ω • set γ = γ + µ ijk ∧ d k dt k ω • set k=k+1 end while • set µω = γ
Note that, in this case, the time derivative of ω in step 2 has to be implemented using the Lie-Derivative with respect to a cartan field on the previously defined jet space of a suitable prolongation order (a more detailed discussion of the choice of the truncation order is given in Section 4.4). It becomes clear that, though the evaluation of ( 23) is a combination of standard operations, its implementation is not obvious in standard computer algebra systems.

Multiplication of two elements

µ ∈ Lq 1 ((Λ(X)) m )
and κ ∈ Lq 2 ((Λ(X)) m )

We consider µ = k≥0 µ k ∧ d k dt k where the µ k 's are matrices whose entries are in Λ q 1 (X), and

κ = k≥0 κ k ∧ d k dt k
where the κ k 's are matrices whose entries are in Λ q 2 (X). The product µκ is evaluated by considering the product µκω for every ω ∈ (Λ(X)) m :

µκω = µ(κω) = k 1 ≥0 µ k 1 ∧ L k 1 τ X   k 2 ≥0 κ k 2 ∧ L k 2 τ X ω   = k 1 ,k 2 ≥0 k 1 k 3 =0 k1 k3 µ k 1 ∧ L k 1 -k 3 τ X κ k 2 ∧ L k 2 +k 3 τ X ω Therefore: µκ = k 1 ,k 2 ≥0 k 1 k 3 =0 k1 k3 µ k 1 ∧ L k 1 -k 3 τ X κ k 2 ∧ d k 2 +k 3 dt k 2 +k 3 (24)
In partcular, setting µ = d dt , we immediately deduce that:

d dt κ = k Lτ X κ k ∧ d k dt k + κ k ∧ d k+1 dt k+1 (25)
The multiplication of two such operators is done with the function Dtwedgeop. Its input is: µ ∈ Lq 1 ((Λ(X)) m ), κ ∈ Lq 2 ((Λ(X)) m ) and output: µκ ∈ Lq 1 +q 2 ((Λ(X)) m ). Its core algorithm using Maple code results to (assuming again that µ has a finite degree l w.r.t. d dt )

• set j = 0 • compute γ = µ ∧ κ while k ≤ l do • compute κ (j+1) = d dt κ (j) according to (25) • set γ = γ + µj+1κ (j+1) end while • set µκ = γ

The operator d

At this point we investigate the definition (17) of the operator d ∈ L1 ((Λ(X)) m ).

We first remark that, for m = 1, if µ is a 0th-order polynomial w.r.t. d dt , i.e. µ = µ0∧ with µ0 ∈ Λ q (X), ( 17) boils down to the usual anti-derivation property of the exterior derivative, i.e. for every ω ∈ Λ p (X)

d(µω) = dµ0∧ d(µ) ω + (-1) q µ0 ∧ dω (26)
Then, going back to the general case, with µ ∈ Lq ((Λ(X)) m )

and ω ∈ (Λ p (X)) m , we have, for i = 1, 2, . . . , m,

[d(µω)] i = d m j=1 k≥0 µ ijk ∧ L k τ X ωj = m j=1 k≥0 dµ ijk ∧ L k τ X ωj +(-1) q m j=1 k≥0 µ ijk ∧ d L k τ X ωj = m j=1 k≥0 dµ ijk ∧ L k τ X ωj +(-1) q m j=1 k≥0 µ ijk ∧ L k τ X (dωj) (27)
We then calculate d(µ)ω by combining ( 17) with ( 27):

[d (µ) ω] i = [d(µ ω) -(-1) q µ dω] i = m j=1 k≥0 dµ ijk ∧ L k τ X ωj (28)
As a consequence, the entries of d(µ) are simply given by

[d(µ)] ij = k≥0 dµ ijk ∧ d k dt k , i, j = 1, 2, . . . , m (29)
Note that operators µ ∈ Lq ((Λ(X)) m ) are specified for Maple as polynomials in DDt, the latter operator being a constant with respect to the Exterior Derivative on the Jet-Manifold X with coordinates (x, ẋ, . . .). Thus, for the chosen implementation in Maple, d boils down to applying the Exteri-orDerivative command of the DifferentialGeometry package to the chosen representation of µ, which readily gives the procedure gdmu (Input: µ0 ∈ Λ q (X), Output: d (µ0) ∈ Λ q+1 (X)).

Iterative increase of truncation order and degree

The Algorithm 3.3 includes a large number of degrees of freedom. There are many choices in the Smith decompositions yielding the vector 1-form ω, but any other choice of a basis of the ideal Ω is, at least algebraically, equivalent. Nevertheless, the set of operators µ such that dω = µω, which is always non empty, generally contains infinitely many elements. A lower bound of its degree w.r.t. d dt is easily computed but there is no a priori upper bound. The only restriction is that the equation d (µ) = µ2 must be satisfied. If we note µ = k≥0 µ k ∧ d k dt k , as a consequence of ( 29) and ( 24), the matrices µ k must satisfy the infinite sequence of differential equations:

dµ k = k k 2 =0 k 1 ≥k-k 2 k1 k -k2 µ k 1 ∧ L k 1 +k 2 -k τ X µ k 2 (30)
for all k ≥ 0. However, since the degree of µ w.r.t. d dt is finite, i.e. only finitely many µ k 's are non zero, on the one hand, and since µ depends only on a finite number of coordinates of X on the other hand, the number of non trivial equations in [START_REF] Van Nieuwstadt | Differential flatness and absolute equivalence of nonlinear control systems[END_REF] is indeed finite. Moreover, [START_REF] Van Nieuwstadt | Differential flatness and absolute equivalence of nonlinear control systems[END_REF] establishes a link between the number of coordinates that are active via the expression of dµ k , and the polynomial degree of µ. Therefore, for a given truncation order 2 and a given degree which are compatible relatively to [START_REF] Van Nieuwstadt | Differential flatness and absolute equivalence of nonlinear control systems[END_REF], solutions µ and M , if they exist, may be found using the above computer algebra program. If not, the truncation order and/or the degree can be increased. Unfortunately, there is no simple answer to the question "is this process ending?". Once µ is determined, according to step 4, a suitable unimodular matrix M ∈ Up d dt satisfying d (M ) = -M µ, has to be found. As a first guess an initial M can always be chosen as an upper triangular Matrix with ones on the main diagonal, the north-eastern entries of M being polynomials in d dt of suitable degree with coefficients depending on the truncated coordinates. If the provided free parameters in M are not sufficient to find a solution, more complicated unimodular matrices can be constructed by left and right multiplication with elementary left and right actions (see e.g. [START_REF] Cohn | Free Rings and Their Relations[END_REF]). If in step 4 no suitable matrix M can be found then it may be necessary to go back to step 3 and to increase the truncation order and/or the degree of µ to introduce additional degrees of freedom. Using Maple 11 and the chosen implementation of the operators, construction methods of general operators µ ∈ and matrices M ∈ for given truncation orders and degrees with respect to d dt could be implemented. This easily allows to iteratively increase the used truncation order and degrees which is essential for the application of this approach.

NON HOLONOMIC CAR

Consider the 3 dimensional system in the x -y plane, representing a vehicle of length l, whose orientation is given by the angle θ, the coordinates (x, y) standing for the position of the middle of the rear axle, and controlled by the velocity modulus u and the angular position of the front wheels ϕ.

ẋ = u cos θ ẏ = u sin θ θ = u l tan ϕ (31)
Since n = 3 and m = 2, n -m = 1 and ( 31) is equivalent to the single implicit equation obtained by eliminating the inputs u and ϕ: F (x, y, θ, ẋ, ẏ, θ) = ẋ sin θ -ẏ cos θ = 0 (32)

We immediately have:

P (F ) = ∂F ∂x + ∂F ∂ ẋ d dt ∂F ∂y + ∂F ∂ ẏ d dt ∂F ∂θ + ∂F ∂ θ d dt = sin θ d dt -cos θ d dt ẋ cos θ + ẏ sin θ .
(33) In the following, all steps of Algorithm 3.3 are performed using the above presented Maple functions. Note that this well-known example has been chosen as it is at the same time challenging enough to illustrate most of the properties of the proposed approach and can at the same time be discussed in detail.

Step 1:

Setting E = ẋ cos θ + ẏ sin θ, we apply the Smith decomposition algorithm (i.e. we apply the function Smith sa) and we get U ∈ R -Smith (P (F )) with

U =   0 0 1 0 1 0 1 E cos θ E d dt -sin θ E d dt   . Thus Û = U 01,2 I2 =   0 1 1 0 cos θ E d dt -sin θ E d dt   with I2 the identity matrix of R 2 . Again, computing Q ∈ L -Smith Û yields Q =   0 1 0 1 0 0 sin θ E d dt -cos θ E d dt 1   .
Multiplying Q by the vector (dx, dy, dθ) T , the last line reads is trivially strongly closed with M = I2, which finally gives the flat output (which we denote here by y f to avoid confusion with the coordinate y) y f = (y, x) T . We have thus recovered the flat output originally obtained in [START_REF] Rouchon | Flatness, motion planning and trailer systems[END_REF][START_REF] Rouchon | Flatness and motion planning: the car with n-trailers[END_REF], up to a permutation of the components of y.

Step 1 b:

Other decompositions of P (F ), given by (33), may indeed be obtained. They are all equivalent in the sense that one decomposition may be deduced from another one by multiplication by a unimodular matrix. However, the resulting vector 1-form ω, contrarily to what happens in the previous example, may not be integrable. Our aim is here to show how the generalized moving frame structure equations [START_REF] Martin | Flat systems[END_REF]

U =   cos θ -1 θ cos 2 θ d dt 1 θ E cos θ sin θ 1 -1 θ sin θ cos θ d dt 1 θ E sin θ 0 0 1   , i.e. Û =   -1 θ cos 2 θ d dt 1 θ E cos θ 1 -1 θ sin θ cos θ d dt 1 θ E sin θ 0 1   The Smith decomposition of Û then yields Q ∈ L -Smith Û with Q =   -tan θ 1 0 0 0 1 -1 θ sin θ cos θ d dt 1 θ cos 2 θ d dt -1 θ E cos θ   .
Step 2:

In this case ω is obtained as Step 4/5:

The simplest unimodular matrix has truncation order 0 and degree 0: Thus, we get as differential of a flat output, which we denote here as above by y f dy f = M ω = -tan(θ)dx + dy + (-x cos 2 (θ) + C2(θ))dθ dθ This one-form is closed and, using pdsolve we obtain the flat output y f = (y -x tan(θ) + C3(θ), θ) T .

CONCLUSIONS

In this contribution the necessary and sufficient conditions for differential flatness of nonlinear control systems have been discussed with the aim of arriving at a computer algebra implementation. To this end general formulae for the used operators have been deduced from its defining relations. It could be shown that all used opeartors can be implemented using, e.g., the computer algebra package Maple 11. However, as operations from differential geometry as well as from algebra were needed, the operators could not be directly implemented but special functions had to be created to implement the action of the results on differential forms or other operators. We want to emphasize that this paper is a first step towards the formal computation of flat output where computational costs are voluntarily ignored. Obtaining more efficient algorithms will be the subject of future work.

1 E

 1 (sin θd ẋ -cos θd ẏ + ( ẋ cos θ + ẏ sin θ)dθ) = 1 E d( ẋ sin θẏ cos θ) and, by (32), identically vanishes on X0. The remaining part of the system,

  may be used to obtain an integrable M ω. Such an example is provided by restarting the right-Smith decomposition of P (F ) by right-mutiplying it by formula sin θ d dt (cos θ) -cos θ d dt (sin θ) = -θ. The Smith decomposition algorithm yields

ω= 0 1 cos 2 θ 0 ∧ 3 : 1 cos 2 θ 0 ∧ 1 cos 2 θ

 120312012 = (ω1, ω2) T = Q (dx, dy, dθ) T = (-tan θdx + dy, dθ) T Its exterior derivative is non zero:dω = (dω1, dω2) T = (-1 cos 2 θ dθ ∧ dx, 0) Tshowing that ω is not closed. The simplest possible operator µ is of truncated order 0 and degree 0:with µ120 = µ1210(x, y, θ)dx+µ1220(x, y, θ)dy+µ1230(x, y, θ)dθ. Evaluation of µω using the Dtwedge function yields µω = µ120 ∧ dθ 0 with µ120 ∧dθ = µ1210dx∧dθ+µ1220dy ∧dθ. The comparison of dω and µω yields µ1210 = 1 cos 2 θ and µ1220 = 0, or µ dx + µ1230(x, y, θ)dθ 0Step By direct computation (using Dtwedgeop): dx + µ1230(x, y, θ)dθ 0From µ 2 = d (µ) we obtain the system of P.D.E.'s ∂ ∂x µ1230(x, y, θ) = 2 sin(θ) cos 3 (θ) ∂ ∂y µ1230(x, y, θ) = 0 whose solution, using pdsolve, is µ1230(x, y, θ) = 2 sin(θ)x cos 3 (θ) + C1(θ). dx + 2 sin(θ)x cos 3 (θ) + C1(θ) dθ 0 0 .

  x, y, θ) = -1 cos 2 θ ∂ ∂y m120(x, y, θ) = 0 ∂ ∂θ m120(x, y, θ) = -2 sin(θ)x cos 3 (θ) + C1(θ)whose solution is (determined again with pdsolve)m120 = -x cos 2 (θ) + C2(θ) (35) with C1(θ) in (34) given by C1(θ) = -d dθ C2(θ). It results that M = 1 -x cos 2 (θ) + C2(θ) 0 1 .

Clearly, in this context, Smith-Decomposition is far from being optimal with respect to computational cost, but is only used to prove the existence of the resulting polynomial matrices. Improvements will be investigated in future work

by truncation, we refer to the elimination of the infinitely many coordinates of which µ doesn't depend