domain: spi phys math info

text:
Alexandre Janot
email: alexandre.janot@onera.fr

Maxime Gautier
email: maxime.gautier@ls2n.fr

Mathieu Brunot
email: mathieu.brunot@onera.fr

Mathieu Brunot Data

Data Set and Reference Models of EMPS

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

EMPS OVERVIEW, MODELLING AND CONTROL

EMPS Overview

The EMPS illustrated Fig. 1 is a standard configuration of a drive system for prismatic joint of robots or machine tools. It is connected to a dSPACE digital control system for easy control and data acquisition using Matlab and Simulink software. Its main components are:

• A Maxon DC motor equipped with an incremental encoder. As we shall see later, this DC motor is position-controlled with a PD controller. • A Star high-precision low-friction ball screw drive positioning unit. Its extremity is equipped with an incremental encoder that is not used in this benchmark. • A load in translation. An accelerometer is mounted on it and is not used in this benchmark.

The measurements provided by the encoder mounted at the extremity of the ball screw and the accelerometer are not used because such measurements are not available in most industrial robots, see e.g. [START_REF] Khalil | Modeling, identification and control of robots[END_REF].

EMPS Modelling

The Inverse Dynamic Model (IDM) of a robot expresses the joint torque/force as a function of the joint position, velocity and acceleration, [START_REF] Khalil | Modeling, identification and control of robots[END_REF]. From Newton's law, the IDM of EMPS is given by τ idm (t) = M q(t) + F v q(t) + F c sign ( q(t)) + of f set, (1) where q, q and q are respectively the joint position, velocity and acceleration; τ idm is the joint torque/force; M is the inertia of the arm; F v and F c are respectively the viscous and Coulomb frictions; of f set is an offset of measurements of τ idm . M , F v and F c are referred to as the dynamic parameters.

Equation ( 1) is linear in relation to the dynamic parameters,

τ idm (t) = [q(t) q(t) sign ( q(t)) 1] θ = φ(t)θ , (2)
where φ(t) is the (1 × 4) observation matrix of basis functions of the IDM and θ = [M F v F c of f set] T is the (4 × 1) vector of the 3 dynamic parameters plus the offset.

The Direct Dynamic Model (DDM) of a robot expresses the joint acceleration as a function of the joint torque/force, position and velocity, [START_REF] Khalil | Modeling, identification and control of robots[END_REF]. In the case of the EMPS, we have

q(t) = τ idm (t) M - F v M q(t) - F c M sign ( q(t)) - of f set M . (3)
Again, the DDM is linear in relation to a set of parameters that is a nonlinear combination of the dynamic parameters i.e.

q(t) = [τ idm (t) -q(t) -sign ( q(t)) -1] θ ddm (4) = φ ddm (t)θ ddm ,
where φ ddm (t) is the (1 × 4) observation matrix of basis functions of the DDM and

θ ddm = [1/M F v /M F c /M of f set/M ] T is the (4 × 1)
vector of the nonlinear combination of the 3 dynamic parameters plus the offset. In the case of EMPS, the 3 dynamic parameters and the offset are divided by M .

In [START_REF] Brunot | Physical parameter identification of a one-degree-of-freedom electromechanical system operating in closed loop[END_REF], a comparison between the IDM and the DDM is made for the identification of the EMPS. However, because the DDM is usually nonlinear with respect to the dynamic parameters for multi degrees-offreedom robots, it is rarely used [START_REF] Swevers | Dynamic model identification for industrial robots[END_REF] and [START_REF] Gautier | A new closed-loop output error method for parameter identification of robot dynamics[END_REF]. Finally, the DDM can be written as a transfer-function (TF) like expression as usually done in robotics, [START_REF] Gautier | A new closed-loop output error method for parameter identification of robot dynamics[END_REF]

q(t) = b 0 s (s + a 0 ) (τ idm (t) -d(t)) , (5)
where s = jω is the Laplace's variable; b 0 = 1/M ; a 0 = F v /M ; d(t) = F c sign( q(t)) + of f set a perturbation encompassing the Coulomb and offset effects.

EMPS Control

As well enlightened by relation ( 5), the EMPS model (1) has a pure integrator and cannot, therefore, be identified in open loop. The system is thus driven by a Proportional-Derivative (PD) controller consisting of an outer loop controlling the position and an inner loop controlling the velocity. The control signal ν is given by ν(t) = k p k v (q r (t) -q(t)) -k v q(t), (6) where k p is the proportional gain and k v is the derivative gain. With the bandwidth of the position loop equal to 20 Hz, the gains haven chosen such as k p = 160.18 1/s and k v = 243.45 V /(m.s -1 ). In [START_REF] Gautier | A new closed-loop output error method for parameter identification of robot dynamics[END_REF], the authors have shown that a PD control is enough to identify the dynamic parameters of robots because an excellent tracking is not needed.

The control signal, coming from the control law, is linked to τ idm by the following relation τ idm (t) = g τ ν(t), (7) where g τ is the drive gain of the EMPS. Although the drive gain is usually given by the manufacturer, it can be identified with special tests [START_REF] Gautier | Global identification of joint drive gains and dynamic parameters of robots[END_REF]. In this case, it has been estimated to g τ = 35.15 N/V .

DESCRIPTION OF DATA SET, EXCITATION OF EMPS, IDENTIFICATION GOALS AND CHALLENGES

Excitation of EMPS

According to [START_REF] Gautier | Numerical calculation of the base inertial parameters of robots[END_REF], M , F v , F c and of f set constitute the set of base parameters i.e. the set of structurally identifiable. Hence, if the EMPS is appropriately excited, then these parameters must be well identified. Of course, if M , F v , F c and of f set are well excited, so are 1/M , F v /M , F c /M and of f set/M .

To excite the dynamic parameters, bang-bang accelerations are applied to the EMPS. Such trajectories are common in robotics, see [START_REF] Khalil | Modeling, identification and control of robots[END_REF], because they allow the user identifying the inertia (or mass), gravity and friction parameters. Indeed, inertia (or mass) is identified while the acceleration is varying (or constant) over a period and gravity and friction parameters are identified while the velocity is constant over a period with a maximum angular range. Finally, to asses the excitation of the system, the condition number of the observation is considered, [START_REF] Gautier | Dynamic identification of robots with power model[END_REF]. In this case, the condition number being close to 30, the dynamic parameters can be considered as well excited.

Description of Data Set

In most applications, the available information is the measurements of the joint position denoted as q m and the control signal denoted as ν m . As previously mentioned, all data are measured by a dSPACE digital control system in order to ensure quality and periodicity.

The sampling frequency is 1kHz while the duration is 25s approximately. q m is measured by an incremental encoder whose the resolution is 12500 counts per revolution. Since this encoder is working in quadrature count mode, its resolution is finally 12500•4 = 50000 counts per revolution. Such a resolution ensures an extremely good quality of q m . Concerning ν m , it is considered as the output of the controller and is directly measured by the dSPACE card.

To compute the velocity involved in the control law (6), the motor position is first filtered by a FIR filter. Hence, at time k, on has

qm (k) = q f m (k) -q f m (k -1) dt with q f m (k) = 0.5 • q m (k) + 0.5 • q m (k -1)
where dt is the sampling time (1ms here); and q f m (k) the filtered motor position at time k.

The EMPS DATA mat file contains all the data required to identify the EMPS. The names of variables are the following:

• qm: the vector of motor position expressed in the load side i.e. in m; • qg: the vector of reference of motor position expressed in the load side i.e. in m; • vir: the vector of motor force expressed in the load side i.e. in N ; • t: the vector of time expressed in seconds;

• kp: the gain of the outer loop of position in 1/s; • kv: the gain of the inner loop of velocity in V /(m.s -1 ); • gtau: the drive gain in N/V .

It is important to stress that all the data were not treated and they can be, therefore, considered as raw. Furthermore, the input of the EMPS is the force i.e. vir MATLAB variable while the output is the position i.e. qm MATLAB variable.

To validate the estimates, cross-test validations must be executed. To run cross-test validations, the EMPS is still excited with bang-bang accelerations but pulses are also The complete diagram of the EMPS is illustrated in Fig. 2 where the variables of signals stored in the mat file are framed in green while the variables of gains are framed in orange. Finally w q is the noise caused by the encoder.

Identification goals

In robotics, the main goal is to identify the physical parameters involved in the IDM i.e. M , F v and F c in order to develop accurate model-based controllers such as the computed torques, see [START_REF] Khalil | Modeling, identification and control of robots[END_REF]. In addition, identifying these physical parameters allows the user having a direct physical assessment of the results.

Always in robotics, another goal is the enhancement of the mechanical or electrical design, e.g. to enhance the haptic rendering and/or backdrivability, to increase the stiffness of the structure, to lower the effect of frictions... It is interesting to notice that identifying the EMPS (or robots) to develop a simple controller (e.g. PID controller) is not usual because the Computer-Aided Design (CAD) values of dynamical parameters provided by manufacturers are enough to design such controllers in order to get rough performances.

Challenges

Despite its simplicity, it is really difficult to use common toolboxes to identify the EMPS: First, like any electro-mechanical systems, the EMPS exhibits a pure integrator which means that it must be identified with a closed-loop process.

Second, because the controller consists of two nested loops, it cannot be written with the following form commonly adopted in toolboxes ν(t) = C(s)(q ref (t) -q(t)) , where C(s) is the transfer function of the controller; and q ref (t) the position reference.

Third, the IDM of EMPS includes a nonlinear effect: the friction. As we shall see later, the friction model which is first assumed to be symmetric is, in fact, asymmetric. Finally, these tree features make the usual toolboxes rather ineffective if it is aimed to obtain accurate dynamic models. This will be emphasized in a later section.

THE STANDARD METHOD FOR ROBOT IDENTIFICATION: IDIM-LS

This section presents the standard approach employed to identify the dynamic parameters of robots termed the Inverse Dynamic Identification Model with Least-Squares estimation (IDIM-LS). Such an approach has been validated on several industrial robots and prototypes see e.g. [START_REF] Gautier | Dynamic identification of robots with power model[END_REF], [START_REF] Khalil | Modeling, identification and control of robots[END_REF], [START_REF] Swevers | Dynamic model identification for industrial robots[END_REF], [START_REF] Gautier | A new closed-loop output error method for parameter identification of robot dynamics[END_REF], (Janot et al., 2014a) and the references therein.

Joint velocity and acceleration estimation

The joint velocities and accelerations have to be calculated from q m in order to build the observation matrix φ as described in [START_REF] Gautier | Dynamic identification of robots with power model[END_REF]. q m is firstly filtered to obtain q. From this filtered position, the derivatives can be calculated with finite differences. The filter type and the cut-off frequency, ω fq , are selected such as q, q, q ≈ (q, q, q) in the range [0, ω fq ]. The filter, which is usually a Butterworth one, is applied in both forward and reverse directions to avoid lag introduction. The rule of thumb for the cut-off frequency is ω fq ≥ 5ω dyn , where ω dyn is the natural frequency of the highest mode of the closed-loop system. The combination of the Butterworth filter and the central differentiation is referred to as the bandpass filtering process.

High-frequency ripples rejection

The IDIM differs from the IDM by an error term v(t), resulting of perturbations coming from measurement noise and modelling errors. The IDIM is then given by:

τ (t) = τ idm (t) + v(t) = φ(t)θ + v(t). (8)
In practice, the torque is perturbed by high-frequency ripples which are rejected by the controller. Those ripples are removed with a parallel lowpass filtering of each basis function at the cut-off frequency ω Fp ≥ 2ω dyn . The choice of ω Fp is involved to keep enough information while avoiding the high frequency noise. Since there is no more useful information beyond the cut-off frequency, the data are also re-sampled by keeping one sample over n d , i.e. n d is the decimation factor. This combination of parallel filtering and re-sampling is referred to as the decimate process. After data acquisition and parallel filtering, we obtain

τ Fp (t) = F p (z -1 )τ (t) (9) = φ Fp q(t), q(t), q(t) θ + v Fp (t),
where z -1 is the backward shift operator and F p is the decimate filter applied to each element of the observation matrix, the torque signal and the error term.

IDIM-LS estimates

The IDIM-LS estimates are given by

θ LS = Φ T Fp Φ Fp -1 Φ T Fp y Fp , (10)
where Φ Fp is the (N ×4) observation matrix resulting from the sampling of φ Fp q(t), q(t), q(t) ; y Fp is the (N × 1) vector of measurements resulting from the sampling of τ Fp ; and N is the number of samples obtained after sampling and decimation.

If the two following conditions hold

Finally, a direct comparison between the simulated and measured data is plotted, see Fig. 6. In this case, only the relative error between the measured and simulated acceleration calculated with 11 has significantly increased whereas the others do not excess 10%. This rise may be explained by a flexibility that was accidentally excited by the pulses whereas the model does not account for it.

Comparison with toolboxes

We tried to identify a dynamic model with the CAPTAIN and the CONTSID toolboxes by assuming that it can be expressed with a simple transfer function, poor results were obtained:

• if ν(t) (i.e. vir Matlab variable) is considered as the input and q m (t) as the output (i.e. qm Matlab variable), a rel err of 74.6 % (!!) with a R 2 T of 0.5 are obtained;

• if ν(t) (i.e. vir Matlab variable) is considered as the input and qm (t) = d dt (q m (t)) (calculated with the diff Matlab function) as the output, a rel err of 47.3 % (!) with a R 2 T of 0.8 is obtained. For instance, a direct comparison between the measured velocity and the one estimated with the rivcbjid function of the CAPTAIN toolbox by considering a linear model is illustrated in Fig. 7. It is clear that the matching is not really good. One notable exception is the identification of the whole closed loop i.e. by taking q r (t) as the input (i.e. qg Matlab variable), q m (t) as the output (i.e. qm Matlab variable) and a second-order structure. In this case, a rel err of 0.005 % (!!) with a R 2 T of 1.0 are obtained. This result is not surprising. Indeed, it simply indicates that d(t) has been well rejected by the controller which implies that the tracking is excellent. According to (Brunot et al., 2018a), this leads to a second-order system for the whole closed loop.

IMPROVING THE IDIM-LS METHOD

Considering other approaches

It is known that cautions must be taken when the LS estimation is used to identify a system operating in closed loop, see e.g. [START_REF] Söderström | Instrumental variable methods for system identification[END_REF], [START_REF] Ljung | System identification: theory for the user[END_REF], [START_REF] Gilson | Optimal instrumental variable method for closed-loop identification[END_REF], [START_REF] Young | Recursive estimation and timeseries analysis: an introduction for the student and practitioner[END_REF] among others. This explains why a tailor-made data filtering is employed to identify electro-mechanical systems.

Recently, a method based on Instrumental Variable approach was proposed in (Janot et al., 2014a), [START_REF] Janot | An instrumental variable approach for rigid industrial robots identification[END_REF] and then further improved in [START_REF] Brunot | An instrumental variable method for robot identification based on time variable parameter estimation[END_REF] and (Brunot et al., 2018a). This approach combines the IDM and DDM and has a good robustness against noises and/or data filtering and converges quickly. In [START_REF] Gautier | A new closed-loop output error method for parameter identification of robot dynamics[END_REF], the authors have proposed and validated an original Input-Error identification method. Indeed, because the system is identified in closed loop, it has been shown that the simulated position which is the usual output in classical Output-Error methods is not really sensitive to parameters' variations. Hence, the authors proposed to circumvent this problem by considering the simulated joint torques instead of the simulated joint positions. Like IV approaches, this input-error method has a good robustness against noises and data filtering and converges quickly. It is worth to note that it is possible to implement the IV and input-error methods mentioned above by using the Simulink file provided by the authors.

In [START_REF] Olsen | Maximum Likelihood Identification of a Dynamic Robot Model: Implementation Issues[END_REF], the authors present an interesting Maximum Likelihood approach that is suitable for identification of industrial robots while in [START_REF] Ramdani | Robust dynamic experimental identification of robots with set membership uncertainty[END_REF] the authors have tested a set membership uncertainty method. A multi-objective genetic algorithm was used to identify the dynamic parameters of a hydraulic robot manipulator in [START_REF] Montazeri | Dynamic modelling and parameter estimation of a hydraulic robot manipulator using a multi-objective genetic algorithm[END_REF]. The authors suggested to employ the Linear Matrix Inequalities to identify the dynamics of a SCARA robot in [START_REF] Calafiore | Robust calibration and control of robotic manipulators[END_REF].

Finally, in [START_REF] Janot | Identification and control of electro-mechanical systems using state-dependent parameter estimation[END_REF], the authors propose an original method by combining the IDM and a black-box identification method termed the State-Dependent Parameter identification method developed by P.C. Young, [START_REF] Young | Recursive estimation and timeseries analysis: an introduction for the student and practitioner[END_REF] chapter 11. In so doing, he authors have shown that the friction model is asymmetric i.e. described by τ f ric = F + v 0 + ( q) + F + c sign(0 + ( q))+ (12) F - v 0 -( q) + F - c sign(0 -( q)) , where 0 + ( q) = q • (1 + sign( q))/2 (resp. 0 -( q) = q • (1sign( q))/2) is a mathematical operator which returns q when q > 0 (resp. q < 0) and 0 otherwise; and F + v and F + c (resp. F - v and F - c ) are the viscous coefficient and Coulomb coefficients for positive (resp. negative) velocities respectively. It is interesting to note that this asymmetry is due to a mechanical fatigue of the screw which is a phenomenon that can be hardly a priori accounted for.

Fig

  Fig. 1. EMPS prototype

Fig. 2 .

 2 Fig. 2. Diagram of EMPS and signals directly applied. The EMPS DATA PULSE mat file contains all the data required to execute the cross-test validations. The pulses applied to the EMPS are stored in the variable pulses N given in N .

Fig. 5 .

 5 Fig. 5. Direct comparison between the simulated and measured data

Fig. 7 .

 7 Fig. 7. Direct comparison between the measured velocity and estimated one obtained with the CAPTAIN toolbox

where E [•] denotes the mathematical expectation, then the IDIM-LS (10) are unique and asymptotically unbiased.

A source code to run the IDIM-LS method

The Script IDIM LS m file presents a simple source code allowing the user running the IDIM-LS method. The code is commented and divided into 8 main parts:

• data loading;

• parameters of filters, both Butterworth and decimate;

• construction of the vector measurements;

• data filtering;

• construction of the observation matrix;

• decimation process;

• calculation of the IDIM-LS estimates;

• plotting and displaying the results.

Once the code executed, something like that must be displayed on the MATLAB command window: Then, a direct comparison between the reconstructed and measured force is plotted, see Fig. 3, and a cross test validation is also performed by using a data set which was not used to run the IDIM-LS method, see Fig. 3. In robotics, it is common to assess the quality of the identification results by computing the relative error given by

where y is the (N × 1) vector of measurements; and y is the (N × 1) vector of estimation. For instance, in the case of the IDIM-LS method, one has y = Φ Fp θLS and y = y Fp when performing direct validations while one has y = Φ CT V Fp θLS and y = y Fp when performing cross test validations; Φ CT V Fp is the (N × 4) observation matrix constructed with the set of data which was not used to run the identification process.

To execute the cross-test validation, the Script IDIM LS CTV m file is provided. It loads the EMPS DATA PULSE mat file and reconstructs the force with the new set of data and the IDIM-LS estimates. Once the code executed, something like that must be displayed on the MATLAB command window:

Validating the IDIM-LS estimates with cross-test validation Absolute value of the diagonal element of R Finally, the whole system is simulated with the estimated values of the dynamic parameters (i.e. with θLS ) and the simulated force, position, velocity and acceleration are compared with the measured ones, see Fig. 5. Again, all the relative errors between the simulated data and the measured ones are calculated with (11), e.g. when comparing the measured and the simulated positions, one has y = q s and y = q m where q s is the simulated position. Once the code executed, something like that must be displayed on the MATLAB command window: Finally, a direct comparison between the simulated and measured data is plotted, see Fig. 5. According to the relative errors calculated with 11, the IDIM-LS estimates can be considered as good enough since these relative errors do not excess 10%.

The second one is the Simulation EMPS Validation m file which calls the Simulink EMPS Rigide Pulse Simulink file generated with the MATLAB 2016b version. This m file is used to perform cross-test validations by comparing the simulated data with the measured one with a trajectory that was not used to run the identification process. This m file loads the EMPS DATA PULSE mat file described in section 2.2.

Once the code executed, something like that must be displayed on the MATLAB command window: By considering the relative error, we can notice that the IDM was slightly improved.