Romain Delpoux

Laurentiu Hetel

Alexandre Kruszewski

Permanent Magnet Synchronous Motor control via Parameter Dependent Relay Control

Keywords: PMSM, relay feedback control, linear matrix inequalities, switched systems

The article presents a novel control strategy for the control of Permanent Magnet Synchronous Motor (PMSM). The approach is motivated by the fact that PMSM are usually controlled by relays and thus only a finite set of control inputs is available. However in classical control design the use of Pulse-Width Modulation (PWM) ignores the relay nature of the actuators. Here we propose a direct relay control. As PMSM may be modeled as Linear Parameter Varying (LPV) systems, we propose a Parameter Dependent Relay (PDR) control. A design based on Linear Matrix Inequalities (LMI) allows to derive the switching surfaces, which depend on the motor position. The theory described is illustrated by simulations results.

I. INTRODUCTION

Permanent Magnet Synchronous Motors (PMSM's) are widely used in the industry. They are more reliable than brush Direct Current (DC) motors and produce a higher torque per volume. From a control engineer perspective, using the field oriented reference frame, referred to as d -q frame, offers several advantages since it provides a simplified structure for the control, by avoiding sinusoidal functions. The control laws in the d -q frame are referred as Field Oriented Control (FOC). They represent the most used control schemes in the field. Several control methodologies have been proposed in the literature : input-output linearization [START_REF] Bodson | High performance nonlinear feedback control of a permanent magnet stepper motor[END_REF], flatness based approach [START_REF] Sira-Ramírez | A Passivity plus Flatness Controller for Permanent Magnet Stepper Motor[END_REF] or sliding mode control [START_REF] Zribi | Static and dynamic sliding mode control schemes for a permanent magnet stepper motor[END_REF], [START_REF] Laghrouche | Robust second order sliding mode control for a permanent magnet synchronous motor[END_REF], [START_REF] Defoort | A Third-Order Sliding-Mode Controller for a Stepper Motor[END_REF]. However most of the existing control laws are usually implemented via Pulse-Width Modulation (PWM) hardware and ignore the relay nature of the actuators. Here we propose a direct relay control which may use the advantages of the switching actuators in power electronics.

Relay feedback was proposed in variable structure systems [START_REF] Tsypkin | Relay Control Systems[END_REF], [START_REF] Edwards | Sliding mode control -Theory and applications[END_REF], [START_REF] Emel'yanov | Variable Structure Automatic Control Systems[END_REF], [START_REF] Utkin | Sliding modes in control optimization[END_REF] and has very interesting robustness properties faced to matched perturbations. In this field, several problems remain unsolved, such as choosing the switching surfaces to optimize system performances, the robustness properties or the size of the domain of attraction. Recently, Linear Matrix Inequalities (LMI) conditions have been proposed for designing a relay control in the case of Linear Time Invariant (LTI) systems with time varying sampling [START_REF] Hetel | Sampled-data control of LTI systems with relay: a convex optimization approach[END_REF]. The method considers LMI techniques used for systems with saturation [START_REF] Boyd | Linear Matrix Inequalities in System and Control Theory[END_REF], [START_REF] Hindi | Analysis of Linear Systems with Saturation using Convex Optimization[END_REF], [START_REF] Hu | Control Systems With Actuator Saturation: Analysis and Design[END_REF], [START_REF] Hu | Analysis and design for discretetime linear systems subject to actuator saturation[END_REF] and convex optimization techniques [START_REF] Liberzon | Switching in Systems and Control[END_REF]. Up to now, it exists few switching surfaces optimization tools applied to nonlinear systems. Here we extend the results from [START_REF] Hetel | Sampled-data control of LTI systems with relay: a convex optimization approach[END_REF] to the case of Linear Time-Varying systems [START_REF] Biannic | Parameter varying control of a high performance aircraft[END_REF], [START_REF] Bara | Parameter-dependent state observer design for affine LPV systems[END_REF], [START_REF] Geromel | Robust stability of time varying polytopic systems[END_REF] which are interesting since they may be used for approximating locally nonlinear systems. The main contribution is to propose an LMI based approach for the design of switching surfaces which are parameterized by time varying parameters. The objective is to find a Parameter Dependent Relay (PDR) controlling directly the commutations of the actuators. When applied to the case of nonlinear systems, the method leads to designing nonlinear switching surfaces guaranteeing the stability of the system under actuator constraints. An illustration of the theory through simulation results using real motor parameters is proposed.

The paper is structured as follow: Section II describes the LPV model of the PMSM for the purpose of relay control. Section III is dedicated to the generic design of the PDR control for LPV systems. The simulation results on the PMSM are proposed in Section IV.

Notations: The identity (or null) matrix with appropriate dimension is denoted by I (or 0). For a square symmetric matrix, M > 0 (M < 0) indicates that M is positive (negative) definite. For a full rank square symmetric matrix M , M -1 denotes the inverse of M . For a symmetric matrix,

M = A B * C . ( 1
)
where * denotes an element that may be inferred by symmetry.

We denote by I N the set {1, 2, . . . , N }.

The symbol ∆ denotes the unit simplex in R N defined by:

∆ = µ ∈ R N : µ i ≥ 0, i ∈ I N , N i=1 µ i (t) = 1 . (2)
For a given set S ∈ R n the symbol conv{S} denotes the closed convex hull of the set. For a symmetric positive definite matrix P ∈ R n×n and a positive scalar c we denote by E(P, c) the ellipsoid:

E(P, c) = {x ∈ R n : x T P x < c}. (3)
Let B(x, c) denotes the open ball centered on x ∈ R n with radius c > 0:

B(x, c) = {y ∈ R n : |x -y| < c}. (4)
Given a constant scalar U > 0, we define the finite set:

Ψ m (U ) = {-U, U } m .
Given a continuous function f : R n → R and a compact set S we denote by:

argmin s∈S f (s) = {s ∈ S : f (s) ≤ f (r), ∀r ∈ S}.
For a set S ⊂ R m , K ∈ R m×n , we denote:

C S (K) = {x ∈ R n : Kx ∈ conv{S}}.
(5)

II. MOTIVATION

The equations [START_REF] Defoort | A Third-Order Sliding-Mode Controller for a Stepper Motor[END_REF] give the standard PMSM model in the phase (or winding) variables:

               L di α dt = v α -Ri α + KΩ sin(n p θ), L di β dt = v β -Ri β -KΩ cos(n p θ), J dΩ dt = K (i β cos(n p θ) -i α sin(n p θ)) -f v Ω -τ, (6)
where v α and v β are the voltages applied to the two phases of the PMSM, i α and i β are the two phase currents, L is the inductance of a phase winding, R is the resistance of a phase winding, K is the back-EMF constant (and also the torque constant), θ is the angular position of the rotor, Ω = dθ/dt is the angular velocity of the rotor, n p is the number of pole pairs (or rotor teeth), J is the moment of inertia of the rotor (including the load), f v is the coefficient of viscous friction and τ represents the load torque. The non-linear state space representation of the system of equations ( 6) is given by:

ẋαβ (t) = f (x αβ , t) + Bv αβ (t) + D̟(t), (7)
where

x T αβ = i α i β Ω , v T αβ = v α v β and ̟ = τ . The function f (x αβ , t) is defined by: f (x αβ ) =           - R L i α + K L Ω sin(n p θ) - R L i β - K L Ω cos(n p θ) K J (i β cos(n p θ) -i α sin(n p θ)) - f v J Ω           , B =     1 L 0 0 1 L 0 0     and D =    0 0 1 J    .
Considering that each motor phase is actuated via commutation, the control vector v αβ may take values in a finite set defined by:

Ψ 2 (V ) = {v αβ ∈ R 2 : v α , v β ∈ {-V, V }},
where V represents the maximal voltage that the control can deliver.

Model in the rotating frame (d -q)

In the phases frame the signals i α and i β vary at n p times the frequency of rotation. This high frequency problem is alleviated by the use of the direct quadrature (d -q) transformation. This transformation changes the frame of reference from the fixed phase axes to axes moving with the rotor. Equation [START_REF] Edwards | Sliding mode control -Theory and applications[END_REF] gives the transformation performed to obtain the rotating frame:

R(θ) = cos(n p θ) sin(n p θ) -sin(n p θ) cos(n p θ) . (8)
The variables in the rotating frame are expressed by:

i d i q = R(θ) i α i β , (9)
and

v d v q = R(θ) v α v β . (10)
The state space representation is then given by:

ẋdq (t) = A dq (Ω(t))x dq (t) + Bv dq (t) + D̟(t), (11)
where x T dq = i d i q Ω , v T dq = v d v q , and:

A dq (Ω(t)) =           - R L n p Ω(t) 0 -n p Ω(t) - R L - K L 0 K J - f v J           .
The matrices B and D remain unchanged. Consider that Ω(t) ranges between known extremal values

Ω(t) ∈ [Ω, Ω].
In this frame the PMSM can be described using a LPV state space representation. The state space representation of the system depends linearly on a vector of time varying parameters: Ω(t), sin(n p θ(t)), cos(n p θ(t)). The model may be represented as follows:

             ẋdq (t) = A(µ(t))x dq (t) + Bv dq (t) + D̟(t), A(µ(t)) = N A i=1 µi(t)Ai, ∀i, µi(t) ≥ 0, i µi(t) = 1, (12)
where N A = 2, with A ∈ [A dq (Ω), A dq (Ω)], and v d,q (t) are defined such that ∀θ ∈ [0, 2π]. The control is defined by:

v dq (t) = k(x dq (t), θ(t)), (13)
where:

k : R 3 × R → v dq (θ). ( 14
)
v dq,4 (θ)

v dq,3 (θ) v dq,2 (θ) v dq,1 (θ) v αβ,1 v αβ,2 v αβ,3 v αβ,4 R(θ) 0 V V -V -V v α v β v d v q
Fig. 1. Finite set of control in the fixe frame and in the rotating frame.

The control vector v dq (t) is a PDR control which takes values in a finite set v dq (θ) defined by:

v dq (θ) = {v dq ∈ R 2 : v dq = R(θ)v αβ , v αβ ∈ Ψ 2 (V )}.
(15) For a given V , the objective is to determine the switching surfaces in the state space which ensure the closed loop stability of the system [START_REF] Hetel | Sampled-data control of LTI systems with relay: a convex optimization approach[END_REF] with the control law [START_REF] Hindi | Analysis of Linear Systems with Saturation using Convex Optimization[END_REF]. The control inputs in the different frames are represented Fig. 1.

III. PARAMETER DEPENDENT RELAY (PDR) CONTROL

We consider the class of Linear Parameter-Varying (LPV) systems with the state-space realization:

ẋ = A(µ(t))x + B(µ(t))u, (16)
where x ∈ R n is the state vector and u ∈ R m is the control vector, the matrices A ∈ R n×n , B ∈ R n×m , are polytopic matrices with the following form:

A(µ(t)) = N i=1 µ i (t)A i , B(µ(t)) = N i=1 µ i (t)B i , (17)
with A 1 , . . . , A N , B 1 , . . . , B N being known constant matrices. The vector µ(t) = µ 1 (t) . . . µ N (t)

T is a vector of real and known parameters which evolves in the unit simplex ∆.

We assume that for each µ ∈ ∆ the control u may only takes values in a finite set which depends on the parameter µ. We define this set of finite values V µ by:

V µ = {v i (µ(t)), i ∈ I k }, v i : ∆ → R m , ∀i ∈ I k . ( 18
)
We assume that conv{V µ } is a non empty set containing the origin in its interior for any µ ∈ ∆.

The objective is to find a PDR control u(x, µ) which locally stabilize the system (16):

u : R m × ∆ → V µ . (19)
The vector set V µ is a known set corresponding to the actuators constraints.

Given that u is discontinuous, the closed loop vector field is discontinuous. Therefore system solutions are considered in the sense of Filippov [START_REF] Filippov | Differential Equations with Discontinuous Righthand Sides: Control Systems[END_REF]. We recall the definition of locally asymptotically stability. One poses:

ẋ = f (t, x), (20)
where

f (t, x) = A(µ(t))x + B(µ(t))u. Definition 1: [1]
The equilibrium point x = 0 of ( 20) is

• stable if, for each ǫ > 0, there is δ = δ(ǫ, t 0 ) > 0 such that:

||x(t 0 )|| < δ ⇒ ||x(t)|| < ǫ, ∀t ≥ t 0 ≥ 0, (21)
• locally asymptotically stable at the origin if it is stable and there is c = c(t 0 ) > 0 such that x(t) → 0 as t → ∞, for all ||x(t 0 )|| < c.

A. Stabilization by PDR control

In this paragraph we show how classical LMIs [START_REF] Boyd | Linear Matrix Inequalities in System and Control Theory[END_REF] may be used to find a control law under the form [START_REF] Ogata | Modern Control Engineering. Instrumentation and controls series[END_REF], that is a PDR control.

Proposition 1: Consider system (16) with the description [START_REF] Laghrouche | Robust second order sliding mode control for a permanent magnet synchronous motor[END_REF]. Assume that the set of LMI conditions:

QA T i + A i Q + B i Y + Y T B T i < -2δQ, i = i, . . . , N, (22)
hold with Q = Q T > 0, Y ∈ R m×n and δ > 0. Then there exists a positive γ such that the function

V (x) = x T Q -1 x, satisfies: ∂V ∂x (A(µ)x + B(µ)u(x, µ)) < -2δV (x), ∀x ∈ Γ \ {0}, ∀µ ∈ ∆, (23)
where:

u(x, µ) ∈ argmin v∈Vµ x T Q -1 B(µ)v, x ∈ Γ \ {0}, (24)
with: Γ = E(Q -1 , γ). Proof: It is desired to design the PDR control (24) that ensures locally the stabilization. The condition [START_REF] Utkin | Sliding modes in control optimization[END_REF] means that the Lyapunov function V (x) = x T Q -1 x satisfies:

2x T Q -1 A(µ) + B(µ)Y Q -1 x < -2δV (x), ∀x = 0, ∀µ ∈ ∆, (25)
i.e. it is a Lyapunov function for the system ( 16) with the state-feedback control law Y Q -1 x.

Since for all µ ∈ ∆, conv{V µ } is non empty and contains the origin in its interior, there exists a positive scalar γ such that:

∀x ∈ E(Q -1 , γ), µ ∈ ∆ ⇒ Y Q -1 x ∈ conv{V µ }, (26) one has: C Vµ (Y Q -1 ) : {x ∈ R n : Y Q -1 x ∈ conv{V µ }. (27)
Let us denote Γ := E(Q -1 , γ). Then for any x ∈ Γ and µ ∈ ∆ there exist k scalars α j (x, µ) ≥ 0, j ∈ I k , with k j=1 α j (x, µ) = 1 such that:

Y Q -1 x = k j=1 α j (x, µ)v j (µ). (28)
From ( 25) and (28) we have:

k j=1 α j (x, µ)2x T Q -1 (A(µ)x + B(µ)v j (µ)) < -2δV (x), ∀µ ∈ ∆ and ∀x ∈ Γ \ {0}. (29)
Considering α j (x, µ) > 0, µ ∈ ∆ and j ∈ I k , there must be at least one j ∈ I k such that:

2x T Q -1 (A(µ)x + B(µ)v j (µ)) < -2δV (x), ( 30
) ∀x ∈ Γ \ {0}, ∀µ ∈ ∆.
Since Γ represents a sub-level set of V (x), local stabilization in Γ with a finite control set is ensured by choosing the control u(x, µ) with the steepest descend of the Lyapunov function:

u(x, µ) ∈ argmin v∈Vµ x T Q -1 B(µ)v. ( 31
)
Since all minimizers ensure the decay of V and the switching surfaces are continuous the the stability (in the sense of Filipov solutions) is insured.

B. Characterization of the domain of attraction

In the previous section we have shown how to guarantee the local stability of the system [START_REF] Labit | SeDuMi Interface 1.02: a Tool for Solving LMI Problems with SeDuMi[END_REF] with the PDR control law (24) in an ellipsoid Γ. In practice it is often useful to provide an estimation of this domain of attraction.

Since µ ∈ V µ , conv{V µ } is non empty and contains the origin, remark that there exists a polytopic region:

Q = conv{q 1 , q 2 , . . . , q p } = {y ∈ R m : h i y ≤ 1, i ∈ I N h }, (32)
such that:

Q ∈ conv{V µ }, ∀µ ∈ ∆ and 0 ∈ Int{Q}.
Using the polytope Q one can adjust the design condition ( 22) so to include an optimization of the domain of attraction.

Proposition 2: Consider system ( 16) and assume that there exists Q = Q T > 0, Y and a positive scalar δ such that:

QA T i + A i Q + B i Y + Y T B T i < -2δQ, i ∈ I N , (33a)   1 h i Y * Q   > 0, i ∈ I N h , (33b)
  eI I * Q   > 0. ( 33c
)
Then the equilibrium point x = 0 of the closed-loop system ( 16)-( 24) is locally asymptotically stable. An estimation of the domain of attraction is provided by the ball B(0, √ ǫ)

with ǫ = 1 e . Proof:
The objective is to find Q and Y such that:

∀x ∈ E(Q -1 , 1) ⇒ Y Q -1 x ∈ Q. ( 34
)
Using Lemma 1 ( [START_REF] Hindi | Analysis of Linear Systems with Saturation using Convex Optimization[END_REF], see appendix), (34) leads to:

1 ≤ min hiKx=1 x T Q -1 x = min i∈IN h h i Y Q -1 Q(Y Q -1 ) T (h i ) T -1 , (35)
which may be characterized by the sufficient set of LMI

1 h i Y * Q > 0, i ∈ I N h . (36)
To guarantee that the ball B(0, √ ǫ) is included in the domain of decay E(Q -1 , 1) of the Lyapunov function we add the constraints [START_REF] Hu | Control Systems With Actuator Saturation: Analysis and Design[END_REF], [START_REF] Hindi | Analysis of Linear Systems with Saturation using Convex Optimization[END_REF]:

eI I * Q > 0. (37)
with e = 1 ǫ . Remark 1: The feasibility of the LMI optimization problem (33) with some matrix Q, guarantees that any LPV system [START_REF] Labit | SeDuMi Interface 1.02: a Tool for Solving LMI Problems with SeDuMi[END_REF] with the control (24) originating from the invariant ellipsoid E(Q -1 , 1) is converging to the origin with a decay rate δ. By minimizing e, the size of the invariant ellipsoid is maximized, in the sense that it contains the ball B(0, √ ǫ) with the maximum radius √ ǫ. This allows to provide an inner ellipsoid approximation of the domain of attraction.

IV. SIMULATION RESULTS

In this section we propose simulation results to illustrate the PDR control, applied to a PMSM. For this purpose, we design a velocity tracking control where we consider that only four control inputs are available. In this article, the control design is considered with the assumption that there is no external torque (i.e. τ = 0).

The objective is to follow a velocity reference y ref . Given that the plant has no natural integrator, an integral action is implemented [START_REF] Ogata | Modern Control Engineering. Instrumentation and controls series[END_REF], to ensure tracking performance. The integral action is given by:

ζ = y -y ref = Cx -y ref , ( 38
)
where ζ is the output of the integrator (ζ(0) = 0). A combination of the state space representation [START_REF] Geromel | Robust stability of time varying polytopic systems[END_REF] and the integral action (38) without torque can be re-written as:

ẋdq ζ ż = A dq (Ω) 0 C 0 A(Ω) x dq ζ z + B 0 B u - 0 I R y ref . ( 39
)
The objective is to design an asymptotically stable system such that at steady state one has ζ = 0, and get the output at steady state y = y ∞ . Moreover, at steady state, one has:

ż∞ = A(Ω)z ∞ + Bu ∞ -Ry ref,∞ , (40)
Considering a constant reference,

y ref = y ref,∞ for t > 0. Define e = z -z ∞ , u c = u -u ∞ .
The dynamic error equation ė = ż -ż∞ is defined by:

ė = A(Ω)e + Bu c , (41)
where u c is a state feedback control law defined by:

u c = Y Q -1 e, (42)
with Q = Q T > 0 and Y ∈ R 2×4 . The parameter Ω takes values in [0, 30] (in rad.s -1 ). For each θ, the control u is constrained to switch among four different values in the set {R(θ)ρ, ρ ∈ Ψ 2 (V )}. The matrix R(θ) is defined by equation [START_REF] Edwards | Sliding mode control -Theory and applications[END_REF]. For the simulations we consider the real parameters of the stepper motor bench developed in LAGIS at École Centrale de Lille [START_REF] Delpoux | Parameter estimation of permanent magnet stepper motors without mechanical sensors[END_REF]. The parameters are L = 9mH, R = 3.01Ω, K = 0.27N.m.A -1 and J = 3.18.10 -4 kg.m 2 .

The number of pole pairs is n p = 50. The power supply provides a maximum voltage v max = 20V . The sampling period for this simulation is constant and equal to 10 -4 s.

The relay value V is equal to v max = 20V leading to a set of two input each taking values in {-V, V }.

Considering as bounded time varying parameter sin(θ(t)), cos(θ(t)) and Ω(t), the system (41) may be rewritten as an LPV system of the form (16) defined by:

ẋ = A(µ)x + Bu, (43)
with 8 vertices.

The control u takes values in the finite set (18) defined by:

V µ = {v i (µ), i ∈ 1, . . . , 4} = {R(θ)ρ, ρ ∈ Ψ 2 (V )}. (44)
In order to optimize the domain of attraction, it is necessary to find a polytopic region Q ∈ V µ . Note that the input u takes values in the vertices of a square of size 2×V submitted to a rotation defined by R(θ). The region containing u is the inscribed disc of radius V . This disc can be approximated by the polytope Q represented for which the vertices q i are given by:

q i+1 = V       cos 2iπ p sin 2iπ p       , i = 0, . . . , p -1. (45)
Each face of the polytope can be characterized by its normal:

h i+1 = q i + q i+1 1 + cos 2π p , i = 0, . . . , p -1. (46)
The simulation are realized using Matlab/Simulink and the LMI are solved using the SeDuMi solver [START_REF] Labit | SeDuMi Interface 1.02: a Tool for Solving LMI Problems with SeDuMi[END_REF]. In order to approximate the inscribed disc by the polytope Q we take p = 30. Choosing a decay rate δ = 1, the LMI solver computes the following Q and Y matrices:

Q =    
The Q matrix defines the PDR control (24) and thus the switching regions. To illustrate the theoretical results, one proposes to compare the continuous state feedback (CSF) control law u c = Y Q -1 e with the PDR control.

Fig. 2 represents the simulation results when no external torque is added. The velocity tracking is accurate in both cases: it shows that at steady state the desired trajectory is tracked with a precision of 0.1rad.s -1 for the PDR control. It must be noted that chattering phenomena appears in the PDR case leading to a tracking error slightly higher. However in this case only four control inputs are used for control. In order to show the robustness of the PDR control, Fig. 3 shows the same tracking but, in this case, an external torque is added from t = 7s until the end of the simulations. Due to the control saturation, the CSF is not able to track the desired velocity while the PDR shows good robustness properties.

V. CONCLUSION

This article has presented a novel control strategy for the control of PMSM. While classical control design using PWM ignores the relay nature of the actuators, in this paper a direct relay control has been considered. Based on an LPV modeling of the PMSM we have proposed a PDR control. LMI have allowed to derive the switching surfaces, which depend on the motor position. Simulation results have been presented based on real PMSM parameters and have shown that, it was possible to control a PMSM using simply the available inputs of the hardware. Moreover the simulations have pointed out the interesting robustness properties of relay control. As future work it would be interesting to consider the perturbations for the control design. A generalization to a larger class of non-linear system is under consideration. Experimental results are also in progress.

APPENDIX

Lemma 1: ( [START_REF] Hindi | Analysis of Linear Systems with Saturation using Convex Optimization[END_REF]). Let V (x) = x T Q -1 x, where Q = Q T > 0, C be a row vector in R n and r be a nonzero scalar. Then the minimum of V along the hyperplane {x|Cx = r} is given by: α r = r 2 CQC T .

(48)

R.

  Delpoux is with Laboratoire des Sciences de l'Information et des Sytèmes (CNRS UMR 7296) INSM team, 8, Bd Louis XIV, 59046 Lille Cedex and Non-A Team, Inria Lille Nord Europe, France romain.delpoux@ensam.eu L. Hetel and A. Kruszewski are with Laboratoire d'Automatique, Génie Informatique et Signal (CNRS UMR 8219), Ecole Centrale de Lille, 59651 Villeneuve d'Ascq cedex, France laurentiu.hetel@ec-lille.fr and alexandre.kruszewski@ec-lille.fr The research leading to these results has received funding from the European Community's 7th Framework Program (FP7/2007-2013) (grant agreement No 257462) HYCON2 Network of Execellence.

Fig. 2 .

 2 Fig. 2. CSF and PDR velocity tracking simulation without perturbation (in rad.s -1 ).

Fig. 3 .

 3 Fig. 3. CSF and PDR velocity tracking simulation with additive external torque (in rad.s -1 ).


Romain Delpoux
email: romain.delpoux@ensam.eu

Richard Béarée
email: richard.bearee@ensam.eu

Adel Olabi
email: adel.olabi@ensam.eu

Olivier Gibaru
email: olivier.gibaru@ensam.eu

Julien Bancalin

Toward on-line robot vibratory modes estimation

Keywords: Manipulator robots, dominating oscillatory mode, parameter estimation, algebraic approach, sliding modes

published or not. The documents may come

I. INTRODUCTION

Manipulator robots are widely used in many fields of industry. Such processes can be used to carry out repetitive tasks, for example, pick and place or assembly tasks. However to improve the performance in terms of speed, such robots are becoming more lightweight and thus more flexible. Speed and accuracy require consideration of vibration of the end effector [START_REF] Olabi | Feedrate planning for machining with industrial six-axis robots[END_REF].

In the literature, solutions are proposed to guarantee trajectories which do not excite the vibration modes of the systems. Among these techniques it can be mentioned the Input Shaping (IS) [START_REF] Singer | Preshaping Command Inputs to Reduce System Vibration[END_REF]. IS methodology consists in the convolution of impulse sequences with a desired system command to produce a shaped input that is used to drive the system [START_REF] Singhose | A comparison of input shaping and time-optimal flexible-body control[END_REF]. However, IS are designed for a given frequency. In industrial applications, where uncertain or time-varying parameters are considered, IS can lose efficiency. IS with parameter adaptations have been proposed, known as Adaptive Input Shaper (AIS). AIS solutions can be designed based on frequency domain [START_REF] Tzes | An adaptive input shaping control scheme for vibration suppression in slewing flexible structures[END_REF] or time domain [START_REF] Bodson | An Adaptive Algorithm for the Tuning of Two Input Shaping Methods[END_REF], [START_REF] Pereira | Adaptive input shaping for manoeuvring flexible structures using an algebraic identification technique[END_REF]. The development of such algorithms has motivated the comparison of two methods for on-line parameter. identification.

Many different methods for the parameter identification exist in the literature. One of the most popular concept is the regression (linear or nonlinear) [START_REF] Söderström | System Identification[END_REF]. Observer based approaches can also be found, such as asymptotic observers using the extended Kalman filter [START_REF] Bittanti | On the parametrization and design of an extended Kalman filter frequency tracker[END_REF] or finite time ones like sliding modes observers. Another approach for the parameter identification is based on an algebraic method. In this paper we propose to compare the algebraic method and the sliding modes for the parameter identification. The objective to characterize the oscillatory behaviour of manipulator robots, i.e. the natural frequency and damping ratio in order to compensate the vibrations. Both methods lead to non-asymptotic convergence estimation procedure.

The algebraic approach was introduced by M. Fliess and H. Sira-Ramírez in [START_REF] Fliess | An algebraic framework for linear identification[END_REF], [START_REF] Fliess | Closed-loop parametric identification for continuous-time linear systems via new algebraic techniques[END_REF]. The method is based on differential algebra and operational calculus. The desired parameters are expressed as a function of integrals of the measured outputs and inputs of the system. It does not need any statistical knowledge of the noise (for instance the assumption that the noise is Gaussian is not required). This method has already been successfully applied to parameter estimation [START_REF] Delpoux | On-line parameters estimation via algebraic method : An experimental illustration[END_REF], [START_REF] Mboup | Parameter estimation for signals described by differential equations[END_REF], [START_REF] Morales | Nonlinear Control for Magnetic Levitation Systems Based on Fast Online Algebraic Identification of the Input Gain[END_REF], [START_REF] Perruquetti | An algebraic approach for human posture estimation in the sagittal plane using accelerometer noisy signal[END_REF], to abrupt change detections and the efficient identification of time delays [START_REF] Belkoura | Parameters estimation of systems with delayed and structured entries[END_REF], [START_REF] Fliess | Algebraic change-point detection[END_REF]. Numerical differentiation of noisy signals may also benefit from this approach, as demonstrated in [START_REF] Mboup | Numerical differentiation with annihilators in noisy environment[END_REF], [START_REF] Riachy | Différentiation numérique multivariable I : estimateurs algébriques et structure[END_REF].

Sliding modes have been popularized by the precursor article of V.-I. Utkin [START_REF] Utkin | Variable structure systems with sliding modes[END_REF]. Their popularity is due to the robustness properties with respect to perturbations and uncertainties [START_REF] Edwards | Sliding mode control -Theory and applications[END_REF], [START_REF] Perruquetti | Sliding mode control in engineering[END_REF]. Chattering phenomenon was a main drawback of the method, however the introduction of highorder sliding modes has overcome this problem. In this paper second order sliding modes observers are presented [START_REF] Davila | Second-order sliding-mode observer for mechanical systems[END_REF]. These observers ensure the finite time convergence to the observed variables, providing equivalent output injection (EOI). The EOI is exploited to obtain the desired parameter estimations [START_REF] Davila | Observation and Identification of Mechanical Systems via Second Order Sliding Modes[END_REF]. Sliding modes have been used in a wide range of application for the control, the observation and the identification [START_REF] Bartolini | Global output-feedback tracking and load disturbance rejection for electrically-driven robotic manipulators with uncertain dynamics[END_REF], [START_REF] Drakunov | Stabilization and tracking control for an extended Heisenberg system with a drift[END_REF], [START_REF] Floquet | Higher-order sliding mode stabilization for a class of nonholonomic perturbed systems[END_REF], [START_REF] Fridman | High-order sliding-mode observation and fault detection[END_REF], [START_REF] Pisano | Cascade control of PM-DC drives via second-order sliding mode technique[END_REF], [START_REF] Riachy | Second order sliding mode control of underactuated Mechanical systems I: Local stabilization with application to an inverted pendulum[END_REF].

The paper is organized as follows: Section II describes the problem statement, the robotic system and a description of vibratory phenomena. The algebraic and sliding mode identifications algorithms are presented III. Finally, the last section presents experimental results on the manipulator robot.

II. PROBLEM STATEMENT

A. Robotic system description

The vibratory study presented in this article is realized on a Stäubli RX-170. This manipulator, depicted in Fig. 1, is a 6 degrees of freedom (DOF) industrial robot with revolute joints. Industrial robots are known to have a good repeatability. The static and dynamic accuracy is far beyond, typically around some millimeters. Dynamic accuracy is mainly deteriorated by joint deformations, which induced low-damped vibrations of the tip of the robot. Classically, each robot axis are submitted to a dominating flexible vibratory mode [START_REF] Olabi | Feedrate planning for machining with industrial six-axis robots[END_REF]. This paper focused on the fast estimation of this vibration for one axis of the robot.

The Stäubli CS8 controller provides access to the actual joints positions and velocities measured through encoders on the motor shafts. The controller is sampled at a sampling frequency of 4ms (250Hz). In order to measure the absolute position of the robot, a Laser Tracker from API inc. is used, where a retroreflector target is mounted on the tool-tip of the robot. The system resolution is 10µ.m -1 and the sampling frequency is 3ms (333Hz).

B. Vibratory dynamics

Considering the first dominating flexible mode, robot axis can be considered as a two mass coupled system, having the rigid link driven by electrical motor through a rotational joint transmission undergoing a viscoelastic joint deformation with a constant stiffness K and a viscosity D (see Fig. 2).

The objective of the paper is to identify the flexible vibratory mode of this axis, i.e. the equivalent harmonic oscillator parameters between the gearbox output and the end of the axis. This behavior can be represented in continuous time domain by the second order system:

θ l (s) θ r (s) = Y (s) U (s) = Kω 2 n s 2 + 2ξω n s + ω 2 n , (1)
where Y (s) is the output (the angular position of the axis), U (s) is the input (angular position of the gearbox output shaft), ξ is the damping ratio, ω n is the natural frequency and K is the gain of the system. This model is considered with the assumption that the system is governed by one vibration mode and that the others have negligible contributions.

Consider the equation ( 1), expressed as a second order differential equation:

ÿ(t) + 2ξω n ẏ(t) + ω 2 n y(t) = Kω 2 n u(t). (2)
In order to simplify the following developments, equation ( 2) is expressed as:

ÿ(t) + α 2 ẏ(t) + α 1 y(t) = α 3 u(t), (3)
where:

ξ = α 2 2 √ α 1 , ω n = √ α 1 , K = α 3 α 1 . (4)
Introducing the variables x 1 = y, x 2 = ẏ, the model ( 3) can be rewritten under the state-space form:

ẋ1 = x 2 , ẋ2 = -α 1 x 1 -α 2 x 2 + α 3 u. (5)
The different representations introduced in this section will be thereafter used to develop the identification algorithms.

III. PARAMETERS ESTIMATION

A. Algebraic Approach

The algebraic estimator presented in this article is based on the basic approach introduced by M. Fliess and H. Sira-Ramírez and can be found in [START_REF] Pereira | Adaptive input shaping for manoeuvring flexible structures using an algebraic identification technique[END_REF] for its application on a second order system. In this article a theoretical development was proposed with the objective to tune an Adaptive Input Shaping. Modifications are proposed to estimate the system's gain. Consider the differential equation (3). Its Laplace Transform is given by:

s 2 Y (s) -sy(0) -ẏ(0) + α 1 (sY (s) -y(0)) +α 2 Y (s) -α 3 U (s) = 0. ( 6
)
The initial conditions which appear in the equation ( 6) are annihilated by taking two derivatives with respect to the complex variable s. One obtains

s 2 d 2 Y ds 2 + 4s dY ds + 2Y + α 1 s d 2 Y ds 2 + 2 dY ds +α 2 d 2 Y ds 2 -α 3 d 2 U ds 2 = 0. (7)
Recall that derivation w.r.t. s in the operational domain translates into multiplication by -t in the time domain.

Multiplication by s in the operational domain corresponds to derivation in the time domain. Applying the linear estimator ( 7) is not appropriate. Derivation amplifies the high frequency components and consequently, the noise contribution.

A simple solution is to make the estimator proper. It is enough to multiply both sides of ( 7) by s -2 , to eliminate the derivation terms and obtain a relationship in function of integral operators.

After algebraic manipulations, one has:

d 2 Y ds 2 + 4s -1 dY ds + 2s -2 Y +α 1 s -1 d 2 Y ds 2 + 2s -2 dY ds +α 2 s -2 d 2 Y ds 2 -α 3 s -2 d 2 U ds 2 = 0. (8)
By application of the Laplace inverse, the equivalent time domain expression is:

η 1 (t) + α 1 η 2 (t) + α 2 η 3 (t) -α 3 η 4 (t) = 0, (9)
in which:

η1(t) = t 2 y(t) -4 t 0 σy(σ)dσ + 2 t 0 σ 0 y(λ)dλdσ, η2(t) = t 0 σ 2 y(σ)dσ -2 t 0 σ 0 λy(λ)dλdσ, η3(t) = t 0 σ 0 λ 2 y(λ)dλdσ, η4(t) = t 0 σ 0 λ 2 u(λ)dλdσ. (10)
As mentioned in [START_REF] Pereira | Adaptive input shaping for manoeuvring flexible structures using an algebraic identification technique[END_REF], the set of equations can be implemented by means of time varying linear (unstable) filters.

From equation ( 9) we have one equation for three unknown parameters. A solution would consist in integrating [START_REF] Davila | Observation and Identification of Mechanical Systems via Second Order Sliding Modes[END_REF] successively twice to obtain a set of three independent equations linear with respect to the parameter to be identified.

The resulting equation [START_REF] Davila | Observation and Identification of Mechanical Systems via Second Order Sliding Modes[END_REF], is linear in the unknown parameters. We rewrite it as:

p 1 (t)θ = q 1 (t), (11)
where

p 1 (t) = -η 1 (t) -η 2 (t) η 3 (t) , q 1 (t) = η 1 (t) and θ = α 1 α 2 α 3 .

B. Sliding Modes Approach

Consider a second order system written under the state space form:

ẋ1 (t) = x 2 (t), ẋ2 (t) = f (t, x 1 (t), x 2 (t), u(t)) + ζ(t, x 1 (t), x 2 (t), u(t)), y(t) = x 1 (t), (12)
where f (t, x 1 (t), x 2 (t), u(t)) is a known function while the uncertainties are concerned in the term ζ(t, x 1 (t), x 2 (t), u(t)).

1) Observer design: The proposed Super-Twisting observer has the form:

ẋ1 (t) = x2 (t) + z 1 (t), ẋ2 (t) = f (t, x 1 (t), x2 (t), u) + z 2 (t), (13)
where x1 (t) and x2 (t) are the state estimations, and the correction variables z 1 (t) and z 2 (t) are the output injections of the form:

z 1 (t) = λ|x 1 (t) -x1 (t)| 1/2 sign(x 1 (t) -x1 (t)), z 2 (t) = αsign(x 1 (t) -x1 (t)). (14)
where z 1 (t) and z 2 (t) are the variables of the Super Twisting Algorithm proposed in [START_REF] Davila | Second-order sliding-mode observer for mechanical systems[END_REF].

At the initial moment, x1 (0) = x 1 (0) and x2 (0) = 0.

Taking e 1 (t) = x 1 (t) -x1 (t) and e 2 (t) = x 2 (t) -x2 (t) the error equations are given by:

ė1 (t) = e 2 (t) -λ|x 1 -x1 (t)| 1/2 sign(x 1 -x1 (t)), ė2 (t) = F (t, x 1 (t), x2 (t), u(t)) -αsign(x 1 (t) -x1 (t)), (15)
where

F (t, x 1 (t), x2 (t), u(t)) = f (t, x 1 (t), x 2 (t), u(t)) - f (t, x 1 (t), x2 (t), u(t))+ζ(t, x 1 (t), x 2 (t), u(t))
. Suppose that the system states can be assumed bounded then the existence is ensured for a constant f + , such that the inequality:

|F (t, x 1 (t), x2 (t), u(t))| < f + , (16)
holds for any possible t, x 1 (t), x 2 (t) and x2 (t) < 2sup|x 2 (t)|.

Let α and λ satisfy the inequalities:

α > f + , β > 2 α -f + (α + f + )(1 + p) (1 -p) , (17)
where p is some chosen constant, 0 < p < 1. Theorem 3.1: Suppose that the parameters of the observer (13), ( 14) are selected according to [START_REF] Fridman | High-order sliding-mode observation and fault detection[END_REF] and condition (16) holds for system [START_REF] Edwards | Sliding mode control -Theory and applications[END_REF]. Then, the variables of the observer converge in finite time to the states of the system, i.e. (x 1 (t), x2 (t)) → (x 1 (t), x 2 (t)). Proof, see [START_REF] Davila | Second-order sliding-mode observer for mechanical systems[END_REF].

2) Parameter Identification Formulation: The parameter identification developed in this section comes from [START_REF] Davila | Observation and Identification of Mechanical Systems via Second Order Sliding Modes[END_REF]. The finite time convergence to the second order sliding mode set ensures that there exists the time constant t 0 > 0 such that for all t ≥ t 0 , from (15) the following identity holds:

0 ≡ ė2 (t) 0 ≡ F (t, x 1 (t), x2 (t), u(t)) -αsign(x 1 (t) -x1 (t)), (18)
notice that F (t, x 1 (t), x2 (t), u(t)) = ζ(t, x 1 (t), x 2 (t), u(t)) because x2 (t) = x 2 (t). Then the equivalent output injection z eq is given by:

z eq ≡ α 1 sign(e 1 (t)) ≡ ζ(t, x 1 (t), x2 (t), u(t)). (19)
Consider that ζ(t, x 1 (t), x2 (t), u(t)) can be decomposed using the regressor notation [START_REF] Söderström | System Identification[END_REF] as:

ζ(t, x 1 (t), x2 (t), u(t)) = θ(t)ϕ(t, x 1 (t), x2 (t), u(t)), (20)
where θ(t) ∈ R n×l is a matrix composed by the value of the uncertain parameters and ϕ(t, x 1 (t), x 2 (t), u(t)) ∈ R l is a known nonlinear function vector.

For the case where the system parameters are time invariant, i.e. θ(t) = θ, the equivalent output injection can be represented in the form:

ζ(t, x 1 (t), x 2 (t), u(t)) = θϕ(t, x 1 (t), x 2 (t), u(t)). (21)
Applied to the article configuration, none of the parameters are known. The equation ( 12) is expressed with:

f (t, x 1 (t), x 2 (t), u(t)) = 0, ζ(t, x 1 (t), x 2 (t), u(t)) = -α 1 x 1 (t) -α 2 x 2 (t) + α 3 u(t). ( 22
)
Equation ( 21) can then be expressed by

p 2 (t)θ = q 2 (t), (23)
where

p 2 (t) = -x 1 (t) -x 2 (t) u(t) , q 2 (t) = ζ(t, x 1 (t), x 2 (t), u(t)) and θ = α 1 α 2 α 3 .
Remark 1: For the purpose of this article, we have considered that none of the parameters were known. Another configuration could have consider that we have nominal parameters expressed by f (t, x 1 (t), x 2 (t), u(t)) and parameter variations to be identified expressed by ζ(t, x 1 (t), x 2 (t), u(t)).

C. Parameters Identification

The proposed approaches for the parameter identification based on an algebraic approach and on sliding modes led to the two similar linear expressions [START_REF] Drakunov | Stabilization and tracking control for an extended Heisenberg system with a drift[END_REF] and [START_REF] Perruquetti | Sliding mode control in engineering[END_REF]. The solution for θ is obtained as a classical solution given by the Least Squares method [START_REF] Cortes-Romero | Algebraic parameter identification for induction motors[END_REF]

θi = t 0 p T i (σ)p i (σ)dσ -1 t 0 p T i (σ)q i (σ)dσ . ( 24
)
for i ∈ {1, 2}. The algebraic parameter estimation is given for i = 1 while the sliding modes parameters identification is given for i = 2.

In the next section, the algorithms for the parameter identification are applied experimentally for the robot system.

IV. EXPERIMENTAL RESULTS

A. Experimental setup

The experimental results are carried out on the manipulator robot described in the Section II. In order to exhibit the oscillatory behavior, the desired trajectory was planned as an angular motion of the first joint (30 • ), represented Fig. 1, while the five others are fixed. The positions measured by the encoders are collected during the displacement. The position of the robot in the cartesian frame is obtained using the kinematic model of the robot, which was previously identified. At the same time, a laser tracker measures the position of the tool-tip. The frames of the robot and the laser tracker have been matched using an Iterated Closest Points (ICP) algorithm [START_REF] Besl | A Method for Registration of 3-D Shapes[END_REF]. The two trajectories are represented Fig.

B. Off-line deformation analysis

In this section is presented the identification of the parameters ω n , ξ and K along the X position. Fig. 4 represents the temporal evolution of the X position. On top, the blue curve represents the encoder measure, the red one, the laser. In order to exhibit the oscillatory behavior, the second subplot represents the difference between the two measures whether the deformation noted ∆θ. This figure shows that indeed the assumption of a second order system for the modelling of the flexible mode makes sense (see section II-B for more details). Using classical results on temporal response of second order system one can easily define the parameters to be compared with the ones identified using the on-line approaches. Indeed the damping ratio can be defined using the formula:

ξ = ln ∆1 ∆2 (2π) 2 + ln ∆1 ∆2 2 , (25)
From Fig. 4 one has d 1 = 0.74mm and d 2 = 0.43mm thus from this response, we have a damping ratio ξ = 8.6%. The figure also show a static gain close to zero. A measure of the oscillation period could give the natural frequency, however to be more precise we propose to compute the Fourier transform of the deformation. The analysis of the Fourier transform highlight the different frequencies. The Fourier transform of this signal is represented Fig. 5 and shows that dominating mode as a pure natural frequency of 8.13Hz.

The off-line analysis previously presented gives an idea of the parameters to be estimated. Note that these parameters have been identified experimentally and cannot be considered as reference parameters. These values are used to give an order of magnitude to be compared with the online estimation presented in the next section.

C. On-line parameter estimation

Before comparing the estimation results, it is important to show the convergence of the sliding mode observer. Indeed as mentioned in Section III, the estimation via sliding modes relies on the design of an observer. The finite time convergence of the observer is based on the assumption of bounded system states. Without loss of generalities, one can assume that the modal deformation of the robot axis is bounded (its derivative equally). The observer gains are chosen in accordance with equation [START_REF] Fridman | High-order sliding-mode observation and fault detection[END_REF]. In Fig. 6 the axis deformation ∆θ is represented with its estimation ∆ θ. The estimation error plotted in the second subplot show the good behavior of the proposed observer. Note that the observer tracking error represents an interesting criterion for the parameter estimation convergence. The experimental comparison of the identification methods is represented Fig. 7. The figure shows that after convergence, the estimations give the same results, whether ξ = 8.4%, ω n = 8.12Hz and K = 0.02 moreover the results are close to the ones obtained off-line. However, for the sliding mode approach, the convergence is faster (around 0.07s for the Sliding modes and 0.12s for the algebraic method). The algebraic approach does not require gain tuning, which can be a complicated task. Note that the parameter convergence time can be obtained using a method based on the standard deviation of the estimated parameters [START_REF] Pereira | Adaptive input shaping for manoeuvring flexible structures using an algebraic identification technique[END_REF].

V. CONCLUSION

In this paper was proposed an experimental comparison between two on-line parameter estimation methods. The algebraic approach was compared to the sliding modes ap- proach. The algorithms have been evaluated experimentally on an industrial robot axis. The objective was to identify the first modal deformation of the robot axis. This vibratory dynamics were modelled by a second order system, where the natural frequency and the modal damping were the parameters to be identified. The comparison between both algorithms has shown similar results in terms of estimated parameter results although the time of convergence is faster for the sliding modes approach. The main difference between these algorithms concern the gains tuning and the algorithmic complexity. The sliding modes structure is simpler to implement, however observers are based on gains which depend on the perturbation amplitudes. The algebraic method, on the other hand, does not depend on parameter tuning but are more complex to implement.

Regarding the convergence time of both the algorithms, experiments have pointed that around one period of the vibration signal was necessary for the estimators to converge. This result can be considered insufficient for most of the vibration shaping methods require half a period to be efficient. The observation can be relativized considering the sampling frequency of the sensor. Higher sampling frequency could lead to faster estimation convergence. One can note that the sensing device used for these estimations, a 3d absolute position system, works at a maximum sampling frequency of 333Hz. A simpler and lower-cost vibration sensor, such as accelerometer, can be used with higher sampling frequency.

Future works on vibration control concern the on-line implementation of these estimation approach for real-time adaptation of input-shaping parameters.

Fig. 1 .

 1 Fig. 1. Stäubli RX-170B robot arm.

Fig. 2 .

 2 Fig. 2. Flexible mode interpretation.

3 .

 3 Although the trajectories are close the figure exhibits the oscillatory behavior at the end effector.

Fig. 3 .

 3 Fig. 3. Experimental measurement of the robot trajectory tracking in the (X, Y, Z) Cartesian frame.

Fig. 4 .

 4 Fig. 4. Temporal evolution of the X position (top) and its deformation (bottom).

Fig. 5 .

 5 Fig. 5. Fourier transform of the X deformation.

Fig. 6 .

 6 Fig.6. Sliding mode deformation observation and its error.

Fig. 7 .

 7 Fig.7. On-line parameter estimation, (blue) sliding modes approach, (red) algebraic approach.