Experimental Identification of the Inverse Dynamic Model: Minimal Encoder Resolution Needed Application to an Industrial Robot Arm and a Haptic Interface

To develop appropriate control laws and use fully the capacities of robots, a precise modelization is needed. Classic models such as ARX or ARMAX can be used but in the robotic field the Inverse Dynamic Model (IDM) gives far better results. In this model, the motor torques depend on the acceleration, speed and position of each joint, and of the physical parameters of the link of the robots (inertia, mass gravity, stiffness and friction). The parametric identification estimates the values of these last parameters. These estimations can also help to improve the mechanical conception during retro-engineering steps... It comes that the identification process must be as accurate and reliable as possible. The most popular identification methods are based on the least-squares (LS) regression “because of its simplicity” (Atkeson et al., 1986), (Swevers et al., 1997), (Ha et al., 1989), (Kawasaki & Nishimura 1988), (Khosla & Kanade 1985), (Kozlowski 1998), (Prufer et al., 1994) and (Raucent et al., 1992). In the last two decades, the IRCCyN robotic team has designed an identification process using IDM of robots and LS regressions which will be developed in the second part of this chapter. This technique was applied and improved on several robots and prototypes (see Gautier et al., 1995 – Gautier & Poignet 2002 for example). More recently, this method was also successfully applied on haptic devices (Janot et al., 2007). However, it is very difficult to know how much these methods are dependent on the measurement accuracy, especially when the identification process takes place when the system is controlled by feedback. So, we ignore the necessary resolution they require to produce good quality and reliable results. Some identification techniques seem robust with respect to measurement noises. They are called “robust identification methods”. But even if they give reliable results, they are only applied on linear systems and, overall, they are very time consuming as can be seen in (Hampel, 1971) and (Hubert, 1981). Finally, it seems difficult to apply them on robots and we do not know how much they are robust with respect to these noises.


Introduction
To develop appropriate control laws and use fully the capacities of robots, a precise modelization is needed. Classic models such as ARX or ARMAX can be used but in the robotic field the Inverse Dynamic Model (IDM) gives far better results. In this model, the motor torques depend on the acceleration, speed and position of each joint, and of the physical parameters of the link of the robots (inertia, mass gravity, stiffness and friction). The parametric identification estimates the values of these last parameters. These estimations can also help to improve the mechanical conception during retro-engineering steps… It comes that the identification process must be as accurate and reliable as possible. The most popular identification methods are based on the least-squares (LS) regression "because of its simplicity" (Atkeson et al., 1986), (Swevers et al., 1997), (Ha et al., 1989), (Kawasaki & Nishimura 1988), (Khosla & Kanade 1985), (Kozlowski 1998), (Prüfer et al., 1994) and (Raucent et al., 1992). In the last two decades, the IRCCyN robotic team has designed an identification process using IDM of robots and LS regressions which will be developed in the second part of this chapter. This technique was applied and improved on several robots and prototypes (see Gautier et al., 1995-Gautier & Poignet 2002 for example). More recently, this method was also successfully applied on haptic devices . However, it is very difficult to know how much these methods are dependent on the measurement accuracy, especially when the identification process takes place when the system is controlled by feedback. So, we ignore the necessary resolution they require to produce good quality and reliable results. Some identification techniques seem robust with respect to measurement noises. They are called "robust identification methods". But even if they give reliable results, they are only applied on linear systems and, overall, they are very time consuming as can be seen in (Hampel, 1971) and (Hubert, 1981). Finally, it seems difficult to apply them on robots and we do not know how much they are robust with respect to these noises.
Another simple and adequate way consists in derivating the CESTAC method (Contrôle et Estimation Stochastique des Arrondis de Calculs developed in Vignes & La Porte, 1974) which is based on a probabilistic approach of round-off errors using a random rounding mode. The third part of this chapter introduces the design and the application of a derivate of the CESTAC method enabling us to estimate the minimal resolution needed for an accurate parametric identification. This theoretical technique was successfully applied on a 6 degrees of freedom (DOF) industrial arm (Marcassus et al., 2007) and a 3 DOF haptic device , the major results obtained will be used to illustrate the use of this new tool of reliability.

General Inverse Dynamic Model
The IDM calculates the joint torques as a function of joint positions, velocities and accelerations. It is usually represented by the following equation: vs =A(q)q+H(q,q)+F q+F sign(q) (1) where is the torques vector of the actuators, q, and are respectively the joint positions, velocities and accelerations vector of each links, () coefficients and an offset of current measurement noted OFFSET j . The kinetic and potential energies being linear with respect to the inertial parameters, so is the dynamic model (Gautier & Khalil, 1990). Equation (1) can thus be rewritten as: where s D (q,q,q) is a linear regressor and χ is a vector composed of the inertial parameters, it is written: with j χ the dynamic parameters of link j and its actuator written: j χ = [ XX j , XY j , XZ j , YY j , YZ j , ZZ j , MX j , MY j , MZ j , M j , Ia j , FV j , FS j , OFFSET j ] T To calculate the dynamic model we do not need all these parameters but only a set of base parameters which are the ones necessary for this computation. They can be deduced from the classical parameters by eliminating those which have no effect on the dynamic model and by regrouping some others. Actually, they represent the only identifiable parameters. Two main methods have been designed for calculating them: a direct and recursive method based on calculation of the energy (Gautier & Khalil, 1990) and a method based on QR numerical decomposition (Gautier, 1991). The numerical method is particularly useful for robots consisting of closed loops. By considering only the b base parameters, equation (2) has to be rewritten as follows: where ( ) D q,q,q is the linear regressor and b χ is the vector composed of the base parameters.

General theory
Generally, ordinary LS technique is used to estimate the base parameters by solving an overdetermined linear system obtained from the sampling of the dynamic model, along a specifically dedicated trajectory (q, q , q ), (Gautier et al., 1995) or (Khalil et al. 2007).
X being the b minimum parameters vector to be identified, Y the torques measurements vector, W the observation matrix and the vector of errors, the system is described as follows: X being the solution of the LS regression, it minimizes the 2-norm of the errors vector . W is a r×b full rank and well conditioned matrix, obtained by tracking exciting trajectories and by considering the base parameters, r being the number of samplings along a given trajectory, r>>b.
Hence, there is only one solution X , (Gautier, 1997) : with W + the pseudo-invert matrix of W.
Standard deviations of the identified parameters, σ î X , are estimated using classical and simple results from statistics considering that the matrix W is deterministic and is a zeromean additive independent noise with a standard deviation such as: where E is the expectation operator and I r the r×r identity matrix. An unbiased estimation of σ is: The variance-covariance matrix of the standard deviation is calculated as follows: However, in practice, W is not deterministic. This problem can be solved by filtering the measurement vector Y and the columns of the observation matrix W.

Data Filtering
Vectors Y and q are measures sampled during an experimental test. We know that the LS method is sensitive to outliers and leverage points so a median filter is applied to eliminate them.
The derivatives q and q are obtained without phase shift using a centered finite difference of the vector q. Knowing that q is perturbed by high frequency noises, which will be amplified by the numeric derivations, a lowpass filter, whose order is greater than 2, is applied on q and q .
The choice of the cut-off frequency ω f is very sensitive because the filtered data fff (q ,q ,q ) must be equal to the vector (q,q,q) in the range [0, ω f ] in order to avoid distortion of the dynamic regressor. So the filter must have a flat amplitude characteristic without phase shift in the range [0, ω f ], with the rule of thumb ω f >10*ω dyn , where ω dyn represents the dynamic frequency of the system. Considering an off-line identification, it is easily achieved with a noncausal zero-phase digital filter by processing the input data through an IIR lowpass Butterworth filter in both the forward and reverse directions. Since the measurement vector Y and matrix W have been filtered, a new filtered linear system is defined: Because there is no more signal in the range [ω f , ω s /2], where ω s is the sampling frequency, the vector Y f and the columns of W f are resampled at a lower rate after a lowpass filtering, keeping one sample over nd samples, in order to obtain the final system to be identified. This process called parallel filtering is done thanks to the "decimate" function available in the Signal Processing Toolbox of Matlab. To have the same cut-off frequency ωf for the lowpass filter, we choose:

Exciting Trajectories
Knowing the base parameters, exciting trajectories must be designed. In fact, they represent the trajectories which excite well these parameters. Design and calculations of these trajectories can be found in (Gautier & Khalil, 1991). If the trajectories are enough exciting, the conditioning number of W, (denoted cond(W)) is close to 1. However, in practice, this conditioning number can reach few hundreds depending of the high number of base parameters. If the trajectories are not enough exciting, the system is ill conditioned, undesirable gatherings occur between inertial parameters and finally the results have no sense whatever the encoder resolutions.

Theory of the CESTAC method
From a theoretical point of view, the LS assumptions are violated in practical applications. In equation (6), the observation matrix W is built from the joint positions q which are measured and from q , q which are often computed numerically from q. Therefore the observation matrix is noisy. Moreover identification process takes place when the robot is controlled by feedback. It is well known that these violations of assumption imply that the LS estimator is biased. An adequate and simple way to evaluate the robustness of the LS estimator with respect to the quantization errors (which contribute significantly to the bias of the estimator) consists in deriving the CESTAC method which is based on a probabilistic approach of round-off errors using a random rounding mode defined below: Definition: each real number x, which is not a floating-point number, is bounded by two consecutive floating-point numbers denoted respectively − X (rounded down) and + X (rounded up). The random rounding mode defines the floating-point number X representing x as being one of two values − X or + X with probability 1/2. Thus, with this random rounding mode, the same program run several times provides different results, due to different round-off errors. Under some assumptions, X can be considered as a quasi-Gaussian distribution (Jezequel, 2004). In our case, the position is perturbed by the encoder resolution. This measurement can be counted up (q + ) as it can be counted down (q -) at each sample. So, we can derive the CESTAC method by building a new position: Equation (15) defines our rounding mode. Then, q CESTAC is filtered thanks to the data filtering previously described. Finally, we build a new linear regression called W CESTAC : Each column of W CESTAC is filtered by using the decimation filter as explain in section xx. Hence, we obtain a new estimation of base parameters denoted X CESTAC : We run this rounding mode N times and because of equation (15), we obtain different results. Hence, one obtains N estimation vectors denoted X CESTAC/k . The mean value is computed thanks to (18): The standard deviations are given by (19): Then, the relative standard deviations are calculated: We calculate the relative estimation error of the main parameters given by the following equation: where: • j CESTAC X is the j th main parameter identified by means of the CESTAC method, • j LS X is the initial value of the j th main parameter identified through LS technique.
Finally, we calculate the relative variation of the norm of the residual torque: where: If the contribution of the noise from the observation matrix W is negligible compared to modeling errors and current measurements noises, then j %e will be practically negligible (less than 1%), as %e(|ρ|) (less than 1%). In this case, the bias of the LS estimator proves to be negligible.
Otherwise, the relative error j %e can not be neglected (greater than 5%), as %e(|ρ|) (greater than 5%). This means that the LS estimator is biased and the results are controversial.
The previous reasoning can be justified by considering the following equation: where: • W ∂ represents the variation of W due to joint positions, velocities and accelerations errors, • Y ∂ represents modeling errors and current measurements noises, • X is the "true" solution. • Y is the "perfect" measurements vector Hence, one obtains: This gives: While a rounding mode is defined for q and if no variations are observed, then, it comes that the solution X is weakly sensitive to quantizing noises. It comes also that the norm of the residual vector the is poorly correlated with W ∂ . We can write that: And, this implies that W is practically uncorrelated with ρ. Finally, the LS estimator is practically unbiased.

6 DOF industrial robot arm
The robot to be identified is presented Figure 1, it is a Stäubli TX40. Its structure is a classical anthropomorphic arm with a 6 DOF serial architecture and its characteristics can be found at the Stäubli web site. The initial encoder resolution is less than 2.10-4 degree per count so the original resolution is more than enough to provide valuable measures. In order to establish the IDM, firstly we define the Modified Dennavit and Hartenberg (DHM) geometric parameters (Khalil & Kleinfinger, 1986), Table 1, based on the schematic of Figure 2. Next, the linear regressor W and the base parameters are computed thanks to the software SYMORO+ (Khalil & Creusot 1997). We notice that this robot has 60 base parameters, some inertial parameters are gathered with others, the letter "R" is added at the end of the regrouped parameters.

Identification Results using the LS Technique
The exciting trajectories, illustrated on Figure 3, consist of polynomial positions which are designed to have constant velocities over a period (the gravity and friction parameters are thus well excited) and to reach maximum admissible accelerations (the inertia parameters are also well excited).

Figure 3. Typical Exciting Trajectory
In our case, W is a 42620x60 matrix and its conditioning number is close to 200, so the trajectories are enough exciting. The cut-off frequencies of the Butterworth filter and the decimate filter are close to 50Hz. This value was found thanks to a spectral analysis. Finally, only 28 essential parameters are enough to characterize the dynamic model of the TX40. Our computations give us the whole 60 parameters, but at the end of our algorithm the parameter with the higher relative standard deviation is removed. Then the LS Method is applied on the new dynamic model until the relative standard deviation of the error vector is above a threshold: σ e ≥1.04σ , σ being the initial relative standard deviation.  Figure 4 shows the measurement vector Y fd , the results of W fd *X fd and the error vector ρ fd . It is obvious that the measured torques vector is matched by the estimated torques vector which validates the identification. The results of the identification are summed up in Table 2, only the most significant parameters for the purpose of this paper are written.

Identification with various lower resolution encoders
The identification protocol of the derivate CESTAC method is designed as explained in section 3. Considering the initial encoder resolution, for the first identification with a reduced resolution, we decide to use 1=2 /10000. Then the resolution is decreased by 1000 down to 2 /1000, 2 /500 is also used. Only 2=2 /5000, 3=2 /2000 and 4=2 /1000, the more relevant, are presented in Table 3.
The exciting trajectories are the identical to those used for the LS technique. The cut-off frequencies of the Butterworth filter and of the decimate filter are closed to 50Hz and the parameters previously identified as essential are the same from those exposed in Table 2.
All tests are carried out 5 times, and each result of these tests is used to compute a set of dynamic parameters. The results summed up in Table 3 Table 3 give that 4 is the threshold beyond which the identification becomes irrelevant.

Identification with lower resolution torque sensors
The protocol is the same as in the previous sections. The ratio between the torque measurement resolution Φ i is of 1/5.  Table 4. Identified values of the mechanical parameters with various torque sensors resolution Table 4 shows that there is no "a priori" lower boundary for the torque measurement resolution.

Development
It appears that the results exposed in Table 2 and those exposed in the columns 1,2 and 3 of Table 3 are close to each others, values in Table 4 are also generally close to the ones in Table  2. All the relative estimation errors, given by (21), have been calculated. Of the various statistic indicators which can be calculated from these results, the relative error of , the error vector, is the best one because it indicates the general accuracy of the identification from a global point of view.
The last line of Table 3 indicates that while the resolution is higher than 1000 counts per revolution the LS identification is consistent and excepted for few parameters, the identified values can be considered as correct.
Overall, it comes that the variations of parameters for the resolutions 1, 2 and 3 are not very important, excepted for MX3, OFFSET4 and OFFSET6 whom respective %ej exceed 20%. The observed major variations for MX3, OFFSET4 and OFFSET6 can be explained by the fact that these parameters do not contribute significantly to the dynamic of the system, and also that they are very sensitive to noises. On the other hand the last line of Table 4 indicates that despite the variations of resolution of torque measurement, the identification still gives accurate results. Indeed, the variations of main inertia parameters (i.e. ZZ j and XX j included in the essential parameters) are less than 5%, while those of the main gravity parameters are less than 3% (excepted the variation of MX3) and the variations of friction parameters are less than 5%.

Presentation and modeling
The CEA LIST has recently developed a 6 DOF high fidelity haptic device for telesurgery (Gosselin et al., 2005). As serial robots are quite complex to actuate while fully parallel robots exhibit a limited workspace, this device makes use of a redundant hybrid architecture composed of two 3 DOF branches connected via a platform supporting a motorized handle, having thus a total of 7 motors. Each branch is composed of a shoulder (link 1), an arm (link 2) and a forearm (link 3) actuated by a parallelogram loop (link 5 and 6). To provide a constant orientation between the support of the handle (link 4) and the shoulder, a double parallelogram loop is used (see Figure ). Figure 6 presents the DHM frames and Table 5 the DHM parameters. Our purpose is to apply the CESTAC method to the serial upper branch of the interface (the handle is disconnected). We note that q 1 , q 2 and q 5 are the active joint positions. The complete modeling of the branches can be found in (Janot et al., 2007a). z 0 ,z 1 z 2 , z 5 z 3 z 4 z 7 , z 8 z 6 x 0 , x 1 x 5 x 2 x 3 ,x 7 , x 8 ,u 3  The joint torque is calculated through current measurement. We have a =NK T I, where a is the joint torque, N the transmission gear ratio, K T the motor torque constant and I the current of the motor. The results are summed up in Table 6. The relative standard deviations are not given, although they have also been calculated, because they don't give important information.
The physical meaning of these parameters is explained in section 2.1. The subscript R means that this parameter is regrouped with others. The symbol * means that these values have been identified through a specific tests.
The parameters having small influence have been removed. We checked that when identified they have a large relative deviation, and that when removed from the identification model, the estimation of the other parameters is not perturbed. Moreover, the norm of the residual torque does not vary (its variation is less than 1%). Finally, only 14 parameters are enough for characterizing the dynamic model of the 3 DOF branch. They are the main parameters of the interface. Direct comparisons have been performed. These tests consist in comparing the measured and the estimated torques just after the identification procedure. An example is given Fig. 4.

www.intechopen.com
Now, we identify the base parameters thanks to the derivate CESTAC method as explained in section 3. For each motor, the encoder resolution equals 4000 counts per revolution. Thus, we have = 2 /4000. The identification process is running 5 times. The results are also summed up in Table . It appears that the results are close to each others. The calculation of the relative estimation error given by (21) shows that the relative errors of the main parameters are less than 1%. Moreover, the variation of the residual torque proves to be negligible. In this case, the LS estimator is practically unbiased and its use is justified. So, the values given in Table can be considered as the "real" values. Now, the limit of the encoder resolution beyond which the bias of the LS estimator can not be neglected is estimated. To do so, we carry out several identification tests using the derivate CESTAC method by increasing in equation (15). Finally, %e j is calculated with the LS estimations given in Table 7. In other words, one supposes that position measurement is increasingly corrupted and we analyze the LS estimator behavior. In our case, small variations appear when equals 2 /100: %e j reaches 1% for main inertia parameters while it reaches 3% for the main gravity parameters. Strong and unacceptable variations occur if is less than 2 /80. Indeed, %e j reaches 5% for main inertia parameters while it reaches 6% for the main gravity parameters. In addition, the variation of the residual torque is greater than 10%. In this case, the bias of the LS can not be neglected and the results are controversial. The results are summed Table . These results prove that the estimation of the base parameters is reliable and practically unbiased. Thus, by writing the dynamic model in the operational space, the apparent mass, stiffness and friction felt by the user can be calculated with a good accuracy. One can assess the quality of the interface as made in (Janot et al., 2007b)

Discussion
The derivate CESTAC method enables us to evaluate the bias of the LS estimator. When the relative errors are very small, its use is fully justified. Otherwise, the LS estimator is biased and the results are controversial.
In the case of the haptic device, small variations appear when equals 2 /100. This value is 40 times as great as the initial value. From (Diolaiti et al., 2005), we know that the stability of the haptic rendering depends explicitly on the encoder resolution. Hence, one can not have a poor encoder resolution. It is one of particularities of haptic devices compared to classical industrial robots. So, due to the high resolution of the encoders, there is a high probability that the LS estimator is practically unbiased. However, an external verification using the derivate CESTAC method is useful.
In addition, we have also tried different probability distributions in equation (15): a uniform distribution and a Gaussian distribution. For the haptic device, one obtains the following results: small variations appear when equals 2 /100 for a uniform distribution while they appear when equals 2 /200 for a Gaussian distribution. The use of a Gaussian distribution means that the probability such that the position error varies between -and + is close to 67%. Compared with a uniform distribution and the rounding mode defined in this paper, this distribution is pessimistic. It can be considered as the worst of cases. Hence, a classical Gaussian can be used to define the rounding mode. Equation (15)  where: • q is the measured motor position, • the encoder resolution, • N(0, 2 ) a Gaussian distribution.

Conclusion
Considering the importance for the robotic applications of a precise parametric modelization, especially to optimize the performances, it is clearly important that the method used for this modelization have too be accurate and unbiased. The CESTAC method let us decrease virtually the resolution of the robot sensors in order to evaluate the bias of the LS estimation. The different results presented throughout the last pages show that the identification process give reliable results and that this technique does not require unusually accurate encoder resolution according to the values found for the lower resolution limit. It is also important to recall that this technique makes it possible to assess the qualities and drawbacks of haptic interfaces during reengineering process.
In the future the CESTAC method will be mixed up with the Direct and Inverse Dynamic Model method (DIDIM) to build up an identification technique requiring less measurements and hence being less dependent of experimental conditions.