Mikhael Tannous
email: mikhael.tannous@ec-nantes.fr

Stéphane Caro
email: stephane.caro@irccyn.ec-nantes.fr

Alexandre Goldsztejn
email: alexandre.goldsztejn@univ-nantes.fr

Sensitivity Analysis of Parallel Manipulators using an Interval Linearization Method

Keywords: parallel manipulators, sensitivity analysis, interval linearization method

The subject of this paper is about an interval linearization method for the sensitivity analysis of manipulators to variations in their geometric parameters. First, the proposed method is presented. Then, three manipulators are used as illustrative examples: The five-bar mechanism, the 3-RRR planar parallel manipulator and the Orthoglide. The benefits and restrictions of the proposed method are also discussed and appropriate indices are derived to show the efficiency of the method. The obtained results are also compared with the results obtained with frequently used methods. The proposed method is simple to implement and provides verified results in low computational time and thus can be applied to complex robots such as the Orthoglide. In particular, the standard linearization method computes unreliable results near singularities, whereas the proposed interval linearization method automatically detects such situations.

Introduction

Many robotic applications require high precision. This precision can not be achieved unless all sources of errors have been identified and integrated in dedicated computational methods. Amongst the well known sources of errors we find manufacturing errors, joint clearances, and backlashes in the actuators. The sensitivity analysis of a manipulator aims at knowing the influence of variations in its geometric parameters and/or actuators on its performance. This information is useful for the evaluation of the pose error of the end-effector and for the tolerance synthesis. The three main methods used for the sensitivity analysis to variations in geometric parameters and joint clearances of serial and parallel robots are: (i) the probabilistic methods [START_REF] Rao | Probabilistic approach to manipulator kinematics and dynamics[END_REF]; (ii) the linearization methods [START_REF] Caro | Sensitivity analysis of 3-RPR planar parallel manipulators[END_REF] and (iii) the interval analysis methods [START_REF] Merlet | Interval analysis for certified numerical solution of problems in robotics[END_REF].

As indicated in [START_REF] Muszynska | Chaotic responses of unbalanced rotor/bearing/stator systems with looseness or rubs[END_REF], science was based in the last centuries on deterministic ideas, which consist in believing that since each phenomenon is due to a determined cause, the behavior can be predicted. However, since it is usually impossible to gather enough information for the behavior to be predictable, probabilistic theories were born. The probabilistic approaches are well known in several engineering fields such as nonlinear dynamics and robotics. Recently, in [START_REF] Kaupp | Humanrobot communication for collaborative decision making -a probabilistic approach[END_REF], a probabilistic approach was used in the field of Human-Robot communication. In [START_REF] Rao | Probabilistic approach to manipulator kinematics and dynamics[END_REF], the sensitivity analysis for a two-link planar manipulator and the Stanford arm was conducted according to a probabilistic method. The authors pointed out that the probabilistic model of the kinematics and dynamics of these robots will be reliable if a sufficient number of experiments is conducted, which means that the results will not be verified. The probabilistic approaches present the advantage of having a competitive computational time while being simple and general. However, they do not provide any verified result.

Other methods for sensitivity analysis do exist such as the linearization method, which consists in analyzing the sensitivity of the linearized system. The linearization method is a numerical method that does not provide verified results [START_REF] Gregorio | Influence of the geometric parameters of the 3-UPU parallel manipulator on the singularity loci[END_REF][START_REF] Fan | Sensitivity analysis of the 3-PRS parallel kinematic spindle platform of a serial-parallel machine tool[END_REF][START_REF] Han | Kinematic sensitivity analysis of the 3-UPU parallel mechanism[END_REF]. In [START_REF] Caro | Sensitivity analysis of 3-RPR planar parallel manipulators[END_REF][START_REF] Caro | Sensitivity analysis of 3-RPR planar parallelmanipulators[END_REF][START_REF] Caro | Sensitivity analysis of planar parallel manipulators[END_REF][START_REF] Binaud | Sensitivity comparison of planar parallel manipulators[END_REF], the authors worked on the sensitivity analysis of planar parallel manipulators using the linearization method. This method is suitable for small uncertainties and as long as the robot is not close to singular configurations. The computed sensitivity may be lower than the real one. For instance, the linearization method was used to analyze the sensitivity of the Orthoglide to variations in its geometric parameters [START_REF] Caro | Sensitivity analysis of the Orthoglide : A three-dof translational parallel kinematic machine[END_REF].

The sensitivity of parallel manipulators to joint clearances is also an important issue. The great impact that joint clearances can have on the endeffector pose of the manipulator was highlighted in [START_REF] Han | Kinematic sensitivity analysis of the 3-UPU parallel mechanism[END_REF]. Furthermore, as mentioned in [START_REF] Zhuang | Self-calibration of parallel mechanisms with a case study on stewart platforms[END_REF], errors due to joint clearances can not be compensated by calibration contrary to manufacturing and assembling errors. To analyze the effect of joint clearances, a complete mathematical model of the joint should be defined. For instance, a joint can be modeled as a journal bearing [START_REF] Wang | Position errors due to clearances in journal bearings[END_REF]. The same model was used in [START_REF] Innocenti | Kinematic clearance sensitivity analysis of spatial structures with revolute joints[END_REF] in order to compute the effect of joint clearances. A general approach was proposed in [START_REF] Voglewede | Application of workspace generation techniques to determine the unconstrained motion of parallel manipulators[END_REF] but turns to be time consuming. A general error-prediction model has been recently developed in [START_REF] Binaud | The kinematic sensitivity analysis of robotic manipulators to joint clearances[END_REF]. Accordingly, two optimization problems were formulated in order to find the maximum positioning and orientation errors of the end-effector of the manipulator due to joint clearances. Again, the previous research works without minimizing their importance are based on algebraic methods that do not provide verified results.

On the contrary, interval analysis based approaches provide verified results while being simpler than their counterparts. In [START_REF] Wu | Interval approach for the modeling of tolerances and clearances in mechanism[END_REF] the authors used a Newton-Raphson based interval analysis method combined with a bisection algorithm in order to analyze the sensitivity of simple mechanisms to joint clearances. However, their approach is complicated and time consuming because it requires a bisection algorithm. Moreover, the method has not been applied to complex mechanisms, and the authors did not give enough information about the behavior of their algorithm close to singular configurations.

The standard linearization method turns to be simple to use and provides good results as long as the variations in the parameters are small enough. Therefore, the standard linearization method may provide bad results.

In this paper, we propose an interval linearization method that combines the simplicity of the standard linearization method and the verification of the results. Moreover, the proposed method is simpler than the approach described in [START_REF] Wu | Interval approach for the modeling of tolerances and clearances in mechanism[END_REF] and allows us to analyze the sensitivity of complex robots to geometric errors and joint clearances.

The paper is organized as follows. Section 2 introduces the proposed interval linearization method for the sensitivity analysis of parallel manipulators. An efficiency index of the method is also developed with regard to the amount of uncertainties. Section 3 deals with the sensitivity analysis of the five-bar mechanism and makes a comparison between the results obtained with the proposed method and the ones obtained with a standard lineariza-tion method. Section 4 deals with the sensitivity analysis of the 3-RRR planar parallel manipulator with the proposed interval linearization method. Section 5 deals with the sensitivity analysis of the Orthoglide, a three degrees of freedom translational parallel manipulator, using the proposed interval linearization method and two geometric models of the manipulator.

The Interval Linearization for Sensitivity Analysis

Interval analysis

Intervals are denoted by brackets: Intervals of reals are [x] = [x, x] with x, x ∈ R and x ≤ x. Intervals of vectors (also called boxes) are [x] = [x, x] with x, x ∈ R n and x ≤ x (the inequality being defined component wise); in this case the interval [x i ] = [x i , x i ] is the i th component of the interval vector.

Interval of matrices are defined similarly by [

A] = [A, A] with A, A ∈ R n×m and A ≤ A.
Interval analysis extends usual real operations (like scalar addition, multiplication, matrix/vector additions and multiplications) to interval objects. These interval operations are defined and implemented so as to rigorously verify the following containment property: The result of the interval operation contains all possible results of the corresponding real operation for real arguments inside the interval arguments. For example, [x][y] + [x] ⊇ {xy + x : x ∈ [x], y ∈ [y]}. Matrix and vector operations are extended in the same way, e.g.

[A][x] ⊇ {Ax : A ∈ [A], x ∈ [x]}, or [A]+[B] ⊇ {A+B : A ∈ [A], B ∈ [B]}.
See [START_REF] Moore | Introduction to Interval Analysis[END_REF] and references therein for more details. These interval operations are implemented in many environments, we used the Matlab library Intlab [START_REF] Rump | INTLAB -INTerval LABoratory[END_REF].

The interval evaluation of an expression generally provides a pessimistic enclosure of the range, which is a central issue in interval analysis. The so called mean-value extension can provide a tighter enclosure of the range, in particular when interval arguments have small widths. It uses interval evaluations of the function derivatives in order to enclose the function's graph within an interval linearization, whose range contains the range of the original function. More formally,

{f (x) : x ∈ [x]} ⊆ f (x) + i [f x i ]([x]) ([x i ] -xi ), (1)
where x ∈ [x], usually chosen as its midpoint, and [f x i ] is an interval extension of the partial derivative of f with respect to x i , which can be computed either as the interval evaluation of its derivatives or by using interval automatic differentiation.

Problem to be solved

We are given a function f : R p × R n → R n , an interval vector [a] for parameters and an approximate nominal solution (ã, x) ∈ R p × R n such that f (ã, x) ≈ 01 . Provided that the parameter domain [a] is small enough (which is verified in practice and in the experimentations presented in this paper), the solution set

Σ(f , [a]) := {x ∈ R n : ∃a ∈ [a], f (a, x) = 0} (2)
can be split in connected components that correspond to the different assembly modes of the robot. We are interested in enclosing the connected component of Σ(f , [a]), which contains the nominal solution. This connected component actually corresponds to the variations in the nominal solution when parameters vary in the domain [a] and plays a key role in the analysis of the robot sensitivity to geometric parameters. Let Σ(f , [a], x) denote this connected component. We aim at computing a box that contains Σ(f , [a], x). Such a box is of great interest for sensitivity analysis, since it quantifies rigorously the sensitivity of the pose of the manipulator with respect to variations in the parameters for a given nominal pose. The smallest box is called its interval hull, denoted by [x Hull ], and is depicted in Figure 1. However, [x Hull ] cannot be computed in general so an outer approximation [x Out ] of it will actually be computed.

The smaller the distance between [x Hull ] and [x Out ], the better the computed enclosure. To formalize this overestimation and analyze it in the experiments, we define the error for the variable x i of the pose vector by:

ǫ x i = 1 - wid[x Hull i ] wid[x Out i ] , (3)
wid denoting the width of the corresponding box. Since . Such an inner approximation leads to an upper bound of ǫ x i .

[x Out ] ⊇ [x Hull ] we have ǫ x i ≥ 0. Moreover, ǫ x i is null if and only if [x Out i ] = [x Hull i ]. Moreover,

The Krawczyk method

Our aim is to enclose the variations in the nominal solution x = x + δx when the nominal parameters change to a = ã + δa. To this end, the linearization method linearizes the function in the neighborhood of the nominal solution, which leads to the following relationship between x and a:

f (x, ã) + f x (x, ã)(x -x) + f a (x, ã)(a -ã) = 0, (4)
where f x and f a are the Jacobian of f with respect to x and a respectively. This is a linear system with respect to x -x, and the linearization method usually explicitly computes its solution as x = xf x (x, ã) -1 f a (x, ã)(aã) (with the assumption that f (x, ã) = 0, i.e., the evaluation of the function at the approximate nominal solution is usually neglected). This leads to accurate approximations of the variations in x provided that uncertainties in parameters are small enough and the matrix f x (x, ã) is well enough conditioned. These necessary conditions are difficult to quantify and verify and therefore the result of the linearization method is not reliable in general.

In the context of interval computations, the inverse of an interval matrix is usually avoided and linear systems are usually solved using iterative methods. We build such an iterative scheme following the idea of the interval Krawczyk2 method: [START_REF] Kaupp | Humanrobot communication for collaborative decision making -a probabilistic approach[END_REF] where the preconditioning matrix C is any non-singular matrix. This fixed point equation is obtained by first left-multiplying Eq. ( 4) by C and adding and subtracting (x-x) (the subtracted (x-x) being collected with its original occurrence, giving rise to the common factor (C f x (x, ã) -I)). Thus Eq. ( 4) and Eq. ( 5) have the same solutions, which can be computed by applying iteratively the fixed point equation ( 5) to an initial guessed solution. Now, we propose an interval linearization procedure by using interval evaluations of functions and derivatives. Equation 5 is reformulated as: [START_REF] Gregorio | Influence of the geometric parameters of the 3-UPU parallel manipulator on the singularity loci[END_REF] where [f ], [f x ] and [f a ] are interval extensions of f , f x and f a respectively (these expressions can be evaluated for interval arguments using interval arithmetic). We denote the right hand side of Eq. ( 6) by [K]([x]). Using the theory of [START_REF] Neumaier | Interval Methods for Systems of Equations[END_REF], it can be proved that:

x = x -C f (x, ã) -(C f x (x, ã) -I)(x -x) -C f a (x, ã)(a -ã),
[x] = x-C [f ](x, ã)-(C [f x ]([x], [a])-I)([x]-x)-C [f a ](x, [a]) ([a]-ã),
(A) If [x] satisfies [x] = [K]([x]) then ∀a ∈ [a], ∃x ∈ [x], f (x, a) = 0. (B) If [x] satisfies [x] = [K]([x]) then the square interval matrix [f x ]([x], [a])
is strongly regular, which implies according to the implicit function theorem that x changes continuously with a.

Algorithm 1: Parametric Krawczyk iteration interleaved with inflation.

Input: f : R p × R n → R n , [a] ∈ IR p , x 0 ∈ R n Output: [x Out ] ∈ IR n 1 k max ← 10; /* Maximum number of iterations */ 2 δ ← 1.01; /* Inter-step inflation ratio */ 3 success ← false; 4 [x] ← x 0 ; 5 repeat 6 [x ′ ] ← mid[x] + δ ([x] -mid[x]); 7 [x] ← [K]([x ′ ]); 8 if ( [x] ⊆ int[x ′ ] ∧ x ∈ [x] ) then success ← true; 9 k ← k + 1; 10 until ( k > k max ); 11 if ( success ) then [x Out ] ← [x] else [x Out ] ← [-∞, ∞] n ; 12 return ( [x Out ] ); Therefore, provided that x ∈ [x] we have [x] ⊇ Σ(f , [a],
x). Note that Property (B) also implies that [x] does not contain any parallel singularity. However, the exact solution of [x] = [K]([x]) cannot be computed in practice, because this iteration does not converge in general in a finite number of steps. In order to overcome this issue, Algorithm 1 interleaves the fixed point iteration

[x] ← [K]([x]
) with an inflation process at Line 6, which allows obtaining the inclusion

[K]([x ′ ]) ⊆ int[x ′ ]
after a finite (usually very small) number of steps. This latter strict inclusion also entails properties (A) and (B) (see Appendix Appendix A) and the returned box

[x Out ] con- tains Σ(f , [a], x). If the inclusion [K]([x ′ ]) ⊆ int[x ′ ] is not satisfied after k max iterations 3 , then the output [x Out ] = [-∞, ∞] n will be interpreted as a fail- ure to enclose Σ(f , [a],
x). In addition to its finite termination, Algorithm 1 presents the advantages to be easy to implement (e.g. it is straightforwardly implemented in Matlab using Intlab [START_REF] Rump | INTLAB -INTerval LABoratory[END_REF]) and to have a low computational complexity. The method we propose in Algorithm 1 can be summarized as the application of a parametric version of the Krawczyk interval operator interleaved with an inflation process, and we call it for short the Krawczyk method in the rest of the paper (and [x Out ] the Krawczyk box ).

Although the theory is correct for any non-singular preconditioning matrix C, a bad choice can prevent the fixed point iteration from converging. One typical preconditioning matrix used in the context of interval analysis is the midpoint inverse preconditioning, i.e.

C ≈ mid[f x ]([x], [a])
-1 where an approximate inverse can be used. This preconditioning matrix allows obtaining good convergence properties and is easy to compute for low dimensional systems. It is used for the illustrative examples presented in the next sections. Related Work.. Intlab [START_REF] Rump | INTLAB -INTerval LABoratory[END_REF] provides the function verifynlss, which is similar to Algorithm 1. The main difference between them consists in the evaluation of [f ](x, [a]): The proposed method uses a mean-value interval extension, while verifynlss uses a simple natural interval extension. For small uncertainties, the former is known to provide sharper results, although the former can occasionally be better. Consider for example the following simple system of two equations, two unknowns and three parameters:

(x 1 + a 1 ) 2 + (x 2 -a 2 ) 2 = a 2 3 (7) (x 1 -a 1 ) 2 + (x 2 + a 2 ) 2 = a 2 3 , (8)
where parameters domains are set to [a 1 ] = 0.5 ± 0.025, [a 2 ] = ±0.025 and [a 3 ] = 1.0 ± 0.025. Starting from the approximate solution x = [x 1 , x 2 ] T = [0.01, 0.85] T , the proposed method computes an enclosure of the solution set while verifynlss diverges. Nevertheless, since the interval natural extension can happen to compute sharper enclosures than the interval mean-value extension, it is natural to use the following hybrid iteration:

[y] = C [f ](x, ã) + (C [f a ](x, [a]))([a] -ã) ∩ C [f ](x, [a]) (9) [x] = x -[y] -(C [f x ]([x], [a]) -I)([x] -x). ( 10
)
This hybrid iteration is named the Krawczyk method and is used in the rest of the paper.

3.

Case study 1: Five-bar mechanism

Manipulator architecture

The five-bar mechanism consists of a closed kinematic chain composed of five revolute joints and five links of lengths l 0 , l 1 , l 2 , l 3 and l 4 , as shown in

Sensitivity analysis of the five-bar mechanism

Since the main sources of errors in the five-bar mechanism are the uncertainties in the leg lengths and the joint clearances, the Cartesian coordinates of point P are the components of the variable vector x of the problem at hand, while its leg lengths are its parameters, namely, the components of vector a. Here, there is no variation in l 0 as the base of the manipulator is supposed to be rigid. As a consequence, only parameters l 1 , l 2 , l 3 and l 4 are considered.

Thus, the five-bar mechanism is described with the following system of equations [START_REF] Wenger | A comparative study of serial and parallel mechanism topologies for machine tool[END_REF]:

x p + l 0 2 -l 1 cos θ 1 2 + (y p -l 1 sin θ 1 ) 2 -l 2 3 = 0 x p - l 0 2 -l 2 cos θ 2 2 + (y p -l 2 sin θ 2 ) 2 -l 2 4 = 0 (11)
that takes the form f(a, x) = 0.

We consider the following nominal values for the geometric parameters of the five-bar mechanism under study: l 0 = 3 m, l 1 = 1 m, l 2 = 1 m, l 3 = 1 m and l 4 = 1 m. Let the uncertainties in the links l i be ∆l i l i = 10 -6

for i = 1, . . . , 4. θ 1 and θ 2 are the actuated joint angles, whereas θ 3 and θ 4 are the passive joint angles. Let point P be above line CD and θ 1 and θ 2 be equal to π/6 rad and 3π/4 rad, respectively. Note that such a configuration is not singular and the mechanism is in a given assembly mode.

Figure 3 illustrates the possible positioning errors of point P in this configuration. The crosses represent the possible positions of point P of the end-effector due to geometric errors in the legs. The crosses are obtained by evaluating the position of the end-effector by means of the direct kinematic model for a finite set of parameter values inside [a] = [a, a]. The values to select in each interval [a], is the lowest value a and the highest value a. This will result in 2 n different possible positions, where n represents the number of parameters. It is obvious that this method cannot provide a verified enclosure. The set of values chosen in each interval cannot guarantee that the calculated pose of the end-effector is a minimum or a maximum. Otherwise, there is no need of proposing an interval linearization method. It is also time consuming, especially for a high number of parameters.

Another drawback of such a technique is the need of the direct kinematic model of the manipulator under study, which may be complicated for some manipulators. The box containing all of the crosses is an inner approximation [x In ] (see Sec 2.2) of the optimal box [x Out ] containing all possible positions of the end-effector due to all combinations of parameter errors. This inner approximation [x In ] and the Krawczyk box are represented with two boxes, the largest one being the Krawczyk box shown in Fig. 3. The two boxes are similar in Fig. 3.2. As a matter of fact, the overestimation of the method is ǫ x = 0.00029%, and ǫ y = 0.00029%, ǫ x and ǫ y being defined with Eq. (3).

However, the higher the uncertainties in the parameters, the less accurate the method. Table 1 represents the evolution of the overestimation of the method with respect to the uncertainties in the leg lengths. Figure 3.2 shows that the overestimation of the Krawczyk method becomes too high for relative errors larger than 10 -2 . This is the price to pay for a rigorous interval linearization. However, manufacturing errors are usually quite smaller, and the Krawczyk method overestimation is then accurate and acceptable. Relative Uncertainties ∆l i /l i , i = 1, . . . , 4 10 -6 10 -5 10 -4 10 -3 10 -2 ǫ x 0.00029 % 0.0029 % 0.0296 % 0.296 % 2.939 % ǫ y 0.00029 % 0.0029 % 0.0296 % 0.295 % 2.898 % To analyze more carefully the error due to the Krawczyk method (given in Table 1), it can be plotted in a logarithmic scale with respect to the uncertainties in the leg lengths as illustrated in Fig. 4. It turns out that the error due to the Krawczyk operator is, in the case of the five-bar mechanism, a linear function of the uncertainties in the leg lengths. The lower the uncertainties, the more accurate the Krawczyk operator.

It is obvious that the method is very accurate for low uncertainties.

Isocontours of the position error of the end-effector

The standard linearization methods are so far the most common methods used for sensitivity analysis of mechanisms. It consists in finding the relationship between the pose error of the end-effector δp and the uncertainties δa in the geometric parameters:

δp = J s δa. ( 13
)
Matrix J s characterizes the mapping between the end-effector pose error δp and the variations in the geometric parameters.

The maximum position error ∆p max is defined as follows:

∆p max = ∆x 2 p + ∆y 2 p , (14)
∆x p = x p -x p , x p and x p being the lower and upper bounds of [x p ], respectively. [x p ] is the box containing the position of the end-effector along the x-axis. Likewise, ∆y p = y p -y p , y p and y p are the lower and upper bounds of [y p ], respectively. [y p ] is the box containing the position of the end-effector along the y-axis.

Figure 5 shows the isocontours of the maximum position error ∆p max of the end-effector of the five-bar mechanism under study and for uncertainties in the parameters equal to ∆l i /li = 10 -4 for i = 1, 2, 3 and 4. Figure 3.3 represents the isocontours obtained with the Krawczyk operator, while Fig. 3.3 represents the isocontours plotted by means of the standard linearization method (all values are multiplied by 1000 for a better clarity).

However, since the chosen working mode of the mechanism contains parallel singularities, and since the sensitivity of a mechanism is highly affected by the singular configurations, it is important to detect the location of the singular configurations for a given working mode and to compare the sensitivity of the mechanism in singular and non-singular configurations using the two different methods.

Figure 6 represents the singularities of the five-bar mechanism for the chosen working mode in its Cartesian workspace.

Figure 7 illustrates a deeper comparison between the values of the sensitivity obtained with the standard linearization and the Krawczyk method by checking the sensitivity at different points along segment AF (see Figure 6) with both methods. For non-singular configurations, the values obtained with the Krawczyk operator are slightly higher than those obtained with the standard linearization method: The latter is based on a linearization of the mathematical model, while the Krawczyk method considers non linearities inside the interval linearization process and provides verified results.

However, the Krawczyk method stops in the vicinity of parallel singularities. The Krawczyk method and the standard linearization method do not diverge in the neighbourhood of serial singularities. Both the methods provide good results close to serial singularities as explained in Sec 2.3.

Nevertheless, the Krawczyk method stops working while the standard linearization method still works but provides numbers that do not make any sense in the vicinity of parallel singularities as shown in Figure 7. As the standard linearization model is badly conditioned when approaching the parallel singularities, the values provided by the standard linearization method tend toward infinity, while the Krawczyk method stops working before giving values that do not make sense. It is noteworthy that the Krawczyk method always provides verified results. Contrary to the Krawczyk method, there is no criterion indicating when the results obtained with the standard linearization method are not reliable.

Case study 2: 3-RRR planar parallel manipulator

Manipulator architecture

The 3-RRR planar parallel manipulator under study is shown in Figure 8. It is composed of three identical limbs attached to the moving platform at one end and to the base at the other end. Each limb contains two links and three revolute joints of axes normal to the moving platform. R stands for revolute joint. As a result, the moving platform can translate along the xaxis and the y-axis and rotate about an axis normal to the moving platform. The equations of motion of the manipulator can be found in [START_REF] Binaud | Sensitivity comparison of planar parallel manipulators[END_REF][START_REF] Bonev | Singularity analysis of 3-dof planar parallel mechanisms via screw theory[END_REF].

Efficiency of the Krawczyk algorithm

The Krawczyk method is used to evaluate the position error of the movingplatform for a given configuration of the manipulator. The corresponding configuration of the manipulator and the geometric parameters of the manipulator are defined thereafter: The crosses in Fig. 4.2 represent the possible positions of the geometric center P of the moving-platform for relative errors in the geometric parameters equal to 10 -6 , i.e., ∆a i /a i = 10 -6 . The Krawczyk box contains all crosses without significant overestimation. Table 2 shows the overestimations ǫ x , ǫ y and ǫ φ along the x-and y-axes and about the z-axis, respectively. Those overestimations are evaluated with Eq. ( 3) and characterize the efficiency of the Krawczyk algorithm for the 3-RRR planar parallel manipulator in a given configuration. It is noteworthy that the higher the relative errors ∆a i /a i , the 2: Overestimations ǫ x , ǫ y and ǫ φ along the x-and y-axes and about the z-axis ∆a i /a i = 10 -6 ∆a i /a i = 10 -5 ∆a i /a i = 10 -4 ∆a i /a i = 10 The crosses in Fig. 4.2 illustrate the possible positions of point P for relative errors in geometric parameters equal to 10 -3 , i.e., ∆a i /a i = 10 -3 . Figure 4.2 also highlights a large overestimation case.

ǫ x ≈ ∆a i a i 1.0034 . ( 15
)

Isocontours of the positioning and orientation errors

The maximal positioning error of the geometric center P of the moving platform is defined as in Eq. ( 14). The maximal orientation error of the moving platform is defined as:

∆φ max = φ -φ (16)
[φ] denotes the box containing the orientation of the moving platform. Figures 11 and12 show the isocontours of ∆p max and ∆φ max throughout the Cartesian workspace of the manipulator evaluated with the Krawczyk method and the standard linearization method, respectively. The nominal geometric parameters of the manipulator are the following l 11 = l 12 = l 21 = l 22 = l 31 = l 32 = r = 1 m and R = 2 m. The nominal orientation of the moving platform is equal to π 6 and the relative errors in the geometric parameters are equal to ∆a i /a i = 10 -5 , i = 1, . . . , 8.

It is apparent that the results obtained with the Krawczyk method and the standard linearization method are similar as long as the manipulator is far from singularities. In the neighborhood of singularities, which appear on the edge of the Cartesian workspace, the Krawczyk method still works while the standard linearization method provides unreliable results. Table 3 shows the comparison between the maximum and minimum values for the sensitivity obtained with the Krawczyk method and the standard linearization method. The maximum positioning and orientation errors computed by the standard linearization method do not make any sense. On the contrary, the Krawczyk method is always reliable. Therefore, we can rely on the Krawczyk method for the sensitivity analysis of 3-RRR planar parallel manipulators even in the vicinity of singular configurations. The standard linearization method cannot provide verified results. This makes the Krawczyk method more robust since it can be applied without a prior investigation of the singular configurations of the robot. The Krawczyk method is reliable near serial and parallel singularities. The Krawczyk method is used in the next section to analyze the sensitivity of the Orthoglide, which is a three degrees-of-freedom translational parallel manipulator.

5. Case study 3: Sensitivity analysis of the Orthoglide

Introduction

The Orthoglide is a three degrees-of-freedom translational parallel manipulator and is depicted in Fig. 13. It is composed of three identical limbs connected to a moving platform at one end and to the base at the other end. Each limb is a serial kinematic chain containing a prismatic joint P , a revolute joint R, a parallelogram joint P a and another revolute joint mounted in series. The frame F b is attached to the base and the frame F p is attached to the moving-platform. The dimensions of the Orthoglide are given in [START_REF] Caro | Sensitivity analysis of the Orthoglide : A three-dof translational parallel kinematic machine[END_REF][START_REF] Pashkevich | Kinematics and workspace analysis of a three-axis parallel manipulator: the Orthoglide[END_REF] and are the following:

• l i = 310.58 mm for i = 1, 2, 3.

• d i = 80 mm for i = 1, 2, 3.

• r i = 31 mm for i = 1, 2, 3.

The manipulator has a regular cubic dexterous workspace of size 200 mm × 200 mm × 200 mm. ρ 1 , ρ 2 , and ρ 3 are the actuated joint variables of the robot, namely, its command variables .

The sensitivity analysis of the Orthoglide is conducted using two parallelogram joint modelings:

1. The parallelogram joints are modeled with their global length, i.e., the variations in l i are considered only, i = 1, 2, 3. 2. The variations in the small and large links of the the parallelogram joints are considered.

First modeling of the parallelogram joints

Based on [START_REF] Caro | Sensitivity analysis of the Orthoglide : A three-dof translational parallel kinematic machine[END_REF], we can create a mathematical model of the Orthoglide. The parallelograms are modeled by their global length as shown in Fig. 14.

This model does not allow us to take into account the orientation errors of the end-effector and the variations in the links of the parallelogram joints, namely, variations in parameters d i1 , d i2 , l i1 and l i2 as they are all replaced by the global length of the parallelogram.

a 1 = -a 1 0 0 T (17a) a 2 = 0 -a 2 0 T (17b) a 3 = 0 0 -a 3 T . (17c)
b 1 = -a 1 + ρ 1 b 1y b 1z T (18a) b 2 = b 2x -a 2 + ρ 2 b 2z T (18b) b 3 = b 3x b 3y -a 3 + ρ 3 T , (18c
c 3 = 0 0 z p -r 3 T , (19c)
x p , y p and z p being the Cartesian coordinates of point P and r i the distance between point C i and point P , i = 1, 2, 3. Thus, the equation system that is suitable for the Krawczyk method takes the form: 5.2 shows the efficiency of the Krawczyk method in the (x p , y p ) plane for errors in the parameters equal to 10 -4 m. The results are similar in planes (x p , z p ) and (y p , z p ). Figure 5.2 depicts the efficiency of the method when the relative geometric errors are equal to 10 -2 m. We can notice that the overestimation is large when the geometric errors become larger than 10 -2 m.

||B 1 C 1 || 2 -l 2 1 = 0 (20) ||B 2 C 2 || 2 -l 2 2 = 0 (21) ||B 3 C 3 || 2 -l 2 3 = 0 ( 22
Table 4 shows the efficiency of the method with regard to the relative errors in the geometric parameters. Figure 16 shows the logarithm of the overestimation ǫ defined in Eq. (3) as a function of the logarithm of the relative errors in the geometric parameters, namely log(∆a/a). It appears that this overestimation is the same along the x-, y-and z-axes. From Fig. 16, the evolution of the efficiency of the Krawczyk method with respect to uncertainties in parameters a is polynomial. Indeed,

ǫ ≈ ∆a a 1.011 . (23)
The Orthoglide has a cubic dexterous workspace as shown in Fig. 17. The Cartesian coordinates of points Q 1 and Q 2 , two opposite vertices of the cubic workspace, are expressed in frame F b as follows: The sensitivity of the Orthoglide along the x-, y-and z-axes is conducted separately. The sensitivity indices are the differences between the upper and lower bounds of the end-effector position along the x-, y-and z-axes and are denoted by ∆x p , ∆y p and ∆z p . From Figs. 5.2-(c), it is apparent that the sensitivity of the Orthoglide's end-effector position is highly dependent of the robot configuration. As a matter of fact, this sensitivity is a minimum in the isotropic configuration, namely, when point P and the origin O of the base frame coincide. It also appears that the sensitivity of the Orthoglide's end-effector position is a maximum to variations in lengths l 1 , l 2 and l 3 , while it is a minimum to uncertainties in the location of points

• Q 1 = (-73.

Second modeling of the parallelogram joints

Figure 19 represents the kinematic model of the Orthoglide considering variations in all its geometric parameters. Nevertheless, links A i B i and C i P are supposed to remain orthogonal to B i1 B i2 and C i1 C i2 , respectively, i = 1, 2, 3.

As shown in Fig. 20, d i1 = d i2 = d i and l i1 = l i2 = l i in the nominal case, i.e., when the variations in the geometric parameters are null.

First, let n i be a unit vector orthogonal to the parallelogram plane, i.e, n i is orthogonal to the parallelogram B i1 B i2 C i2 C i1 , namely,

n i = B i2 C i × B i2 B 11 ||B i2 C i × B i2 B 11 || 2 , i = 1, 2, 3. (24)
Let u i and v i be the unit vectors directed along vector C i P and vector C i2 C i1 , respectively.

v i = C i P ||C i P|| 2 , i = 1, 2, 3. (25)
Vector u i is defined such that The latter amounts to

u i = n i × v i , i = 1, 2, 3. (26)
u i = C i2 C i1 ||C i2 C i1 || 2 , i = 1, 2, 3 (27)
where ||.|| 2 denotes the 2-norm. Contrary to the modeling used in Sec 5.2 that only considers the positioning errors of the Orthoglide's end-effector, the modeling introduced in this section deals with both its positioning and orientation errors as the geometric variations in the parallelogram joints are considered in the mathematical model.

In order to take into account the positioning and orientation errors of the end-effector, we should have six variables. It means that we need to write six independent expressions that connect the variables to the parameters in the form f(a, x) = 0.

Frame F p of origin O p is attached to the end-effector as shown in Fig. 19, whereas frame F b of origin O b is attached to the base.

The transformation matrix from frame F p to frame F b takes the form4 :

T b p = R b p [o p ] b 0 3 1 , (28)
with

R b p =   c φ c θ c φ s θ s ψ -s φ c ψ c φ s θ c ψ + s φ s ψ s φ c θ s φ s θ s ψ -c φ c ψ s φ s θ c ψ + c φ s ψ -s θ c θ s ψ c θ c ψ   . (29)
Euler angles φ, θ and ψ are used to parameterize the orientation of the endeffector. [o p ] b is the position vector of point O p expressed in frame F b and 0 3 denotes the three-dimensional zero vector.

In order to formulate the required expressions for the Krawczyk operator, we first express the Cartesian coordinates of points A i and B i in frame F b . Since joints are supposed to be perfect, B i1 B i2 remains perpendicular to A i B i and we can easily get the Cartesian coordinates of points B i1 and B i2 .

In frame F p , we can get the coordinates of points C i as a function of the pose of the end-effector. The Cartesian coordinates of points C i expressed in frame F b are obtained by means of matrix R b p . As a result, the six equations required for the Krawczyk operator can be formulated as:

||B i1 C i1 || 2 2 -l 2 i1 = 0 (30) ||B i2 C i2 || 2 2 -l 2 i2 = 0, i = 1, 2, 3
The relative geometric errors are supposed to be equal to 10 -4 . Moreover, the positioning error of the end-effector is expressed as: ∆p = ∆x 2 p + ∆y 2 p + ∆z 2 p (31

)
where ∆x p , ∆y p , and ∆z p denote the positioning error of the end-effector along the x-, y-and z-axes, respectively. The orientation error µ r of the end-effector is expressed as:

µ r = arccos T r(R b p ) -1 2 (32)
where T r(R b p ) is the trace of the rotation matrix R b p expressed in Eq. ( 29). The sensitivity analysis of the Orthoglide to the variations in its geometric parameters is conducted along segment Q 1 Q 2 like in Sec. 5.2,. Uncertainties are considered separately in each parameter.

Figures 21 and22 represent the sensitivity of the orientation and position of the end-effector to variations in all geometric parameters along segment Q 1 Q 2 , respectively.

From Figs. 21 and 22, links l i2 , i = 1, 2, 3, have the largest influence on the pose of the end-effector. We can also notice that the position of the end-effector is more sensitive to the variations in the actuated joints than its orientation with regard to the sensitivity to the variations in the other geometric parameters. However, it is apparent that the varitions in the actuated joints have larger influence on the end-effector's pose than variations in link lentghs d i1 and d i2 .

The comparison of these results with those obtained in [START_REF] Caro | Sensitivity analysis of the Orthoglide : A three-dof translational parallel kinematic machine[END_REF] also shows that the Krawczyk operator is efficient and relevant for the sensitivity analysis of the Orthoglide.

Conclusions

The interval linearization method proposed in this paper provides verified results for the sensitivity analysis of serial and parallel manipulators since it is based on a rigorous linearity that takes nonlinearities into consideration. This fact represents a major advantage in comparison with probabilistic and standard linearization methods. The method turns to give accurate results even when the manipulator under study is close to a singular configuration. As a matter of fact, as long as the algorithm converges, it will provide verified and relevant results. It means that we provided an auto validation algorithm: If it provides results, they will be verified even close to singular configurations. It is a major advantage in comparison with the standard linearization method where results are provided even too close to singular configurations, but are not reliable. It has also proved to be general since it works for simple mechanisms such as the five-bar mechanism, and complex ones such as the Orthoglide.

The method proposed in this paper was implemented using the Interval Laboratory Intlab [START_REF] Rump | INTLAB -INTerval LABoratory[END_REF] that works under Matlab. Implemented on a 1.86 GHz Intel Pentium Dual-core CPU, and a 3 GHz memory, thirty seconds were required to draw the isocontours for the five-bar mechanism, whereas up to one minute was required to analyze the sensitivity of the Orthoglide. However, when the full mathematical model of the Orthoglide was used, the computational time had been higher, since the corresponding model is very complex and implements hundreds of trigonometric functions in each of the six equations. We can say that our method is faster that the algorithm presented in [START_REF] Wu | Interval approach for the modeling of tolerances and clearances in mechanism[END_REF] since no bisection is needed.

As a conclusion, the proposed method for the sensitivity analysis of serial and parallel manipulators provides verified results and the computation time is acceptable. Later on, the method will be used to compute the maximum pose errors of serial and parallel manipulators due to joint clearances.

Figure 1 :

 1 Figure 1: The set Σ(f , [a], x) is depicted in light; the box in dashed line is its interval hull [x Hull ]; the box in full line [x Out ] is an outer approximation of the latter two boxes. The dotted box [x In ] is computed by solving the system for 50 randomly chosen parameter values inside [a] (solutions are illustrated with crosses) and considering the minimum and maximum values of each variable, it is included inside [x Hull ].

Figure 2 .

 2 Figure 2. The two revolute joints attached to link l 0 are actuated. Point P denotes the end of the manipulator, x p and y p are its Cartesian coordinates.

Figure 2 :

 2 Figure 2: Five-bar mechanism

  (a) ∆l i /l i = 10 -6 (b) ∆l i /l i = 10 -2

Figure 3 :

 3 Figure 3: Accuracy of the Krawczyk method

Figure 4 :

 4 Figure 4: Efficiency of the Krawczyk method

Figure 5 :

 5 Figure 5: Isocontours of the maximum position error of the end-effector throughout the manipulator Cartesian workspace

Figure 6 :

 6 Figure 6: Singularities for the selected working mode

Figure 7 :

 7 Figure 7: Comparison between the linearization and Krawczyk method along segment AF

Figure 8 : 3 - 2 .

 832 Figure 8: 3-RRR planar parallel manipulator

8 Figure 9 :

 89 Figure 9: Possible positions of point P

Figure 10 :

 10 Figure 10: Efficiency of the Krawczyk algorithm for the 3-RRR planar parallel manipulator

Figure 10

 10 Figure 10 illustrates the linear relation between the logarithm of the positioning error of the moving-platform along the x-axis and the logarithm of relative errors ∆a i /a i . From Fig 10, it turns out that:

Figure 11 :

 11 Figure 11: Isocontours of ∆p max [10 -5 m]

Figure 13 :

 13 Figure 13: CAD modeling of the Orthoglide

Figure 14 :

 14 Figure 14: Parameterization of the i th parallelogram joint

Let b 1 ,

 1 b 2 and b 3 denote the Cartesian coordinates of points B 1 , B 2 and B 3 expressed in F b :

T

  ) b 1y , b 1z , b 2x , b 2z , b 3x and b 3y represent the uncertainties in the positions of points B i for i = 1, 2, 3. Likewise, let c 1 , c 2 and c 3 denote the Cartesian coordinates of points C 1 , C 2 and C 3 expressed in F b . Therefore, c 1 = x p -r 1 0 0

)

  Let us consider the uncertainties in the following twelve parameters: b 1y , b 1z , b 2x , b 2z , b 3x , b 3y , r i , and l i for i = 1, 2, 3. The direct method requires 2 12 pose errors to be computed. It gives us the boundaries of the possible end-effector positions represented with crosses in Figs. 5.2 and 5.2.

Figure

  Figure 5.2 shows the efficiency of the Krawczyk method in the (x p , y p ) plane for errors in the parameters equal to 10 -4 m. The results are similar in planes (x p , z p ) and (y p , z p ). Figure5.2 depicts the efficiency of the method when the relative geometric errors are equal to 10 -2 m. We can notice that the overestimation is large when the geometric errors become larger than 10 -2 m.Table4shows the efficiency of the method with regard to the relative errors in the geometric parameters.

  21, -73.21, -73.21) [mm] • Q 2 = (126.79, 126.79, 126.79) [mm]Here the sensitivity of the Orthoglide to variations in its geometric parameters is analyzed along segment Q 1 Q 2 .

Figure 15 :

 15 Figure 15: Orthoglide's end-effector positions in the (x p , y p ) plane for relative geometric errors equal to (a) 10 -4 and (b) 10 -2

Figure 16 :

 16 Figure 16: Efficiency of the Krawczyk method for the sensitivity analysis of the Orthoglide

Figure 17 :

 17 Figure 17: The Orthoglide's workspace

Figures 5 . 2 -

 52 (c) represent the sensitivity of the Orthoglide's end-effector pose to variations in the twelve geometric parameters that are considered along the x-, y-and z-axes. The uncertainties in parameters l i and r i , i = 1, 2, 3, are chosen arbitrarily to be equal to 1/1000 times the leg length and the uncertainties in parameters b 1y , b 1z , b 2x , b 2z , b 3x , and b 3y are supposed to be equal to 10 -3 mm.

  B i , namely to variations in geometric parameters b 1y , b 1z , b 2x , b 2z , b 3x and b 3y .

Figure 19 :Figure 20 :

 1920 Figure 19: Schematic representation of the Orthoglide with variations in all its geometric parameters

Figure 21 :

 21 Figure 21: Sensitivity of the orientation of the end-effector to variations in geometric parameters along segment Q 1 Q 2 .

Figure 22 :

 22 Figure 22: Sensitivity of the orientation of the end-effector to variations in geometric parameters along segment Q 1 Q 2 .

Table 1 :

 1 Efficiency of the Krawczyk operator

Table 3 :

 3 Comparison between the minimum and maximum values of ∆p max and ∆φ max obtained with the Krawczyk and the standard linearization methods

		Krawczyk Method Standard linearization method
	Min ∆p max [m]	4.4899 10 -5	4.4878 10 -5
	Max ∆p max [m]	43.4090 10 -5	2.4889 10 -2
	Min ∆φ max [rad]	3.6326 10 -5	3.6314 10 -5
	Max ∆φ max [rad]	26.4128 10 -5	2.1387 10 -2

Table 4 :

 4 Efficiency of the Krawczyk method for the sensitivity analysis of the Orthoglide

	Relative Uncertainties
	10 -4	10 -3	10 -2
	ǫ x 0.071632 0.719432 7.559653
	ǫ y 0.071632 0.719432 7.559653
	ǫ z 0.071628 0.719072 7.525974

A manipulator is usually modeled by f (a, u, x) = 0 where a are geometric parameters, u are commands and x are positions. Since u is fixed in the scope of our study, we include it into the geometric parameters in order to also deal with uncertainties in the command. Note that x and u play symmetric roles, therefore the dual problem consisting in enclosing the set of commands necessary to keep a pose constant in spite of variations in geometric parameters could be solved similarly.

Other operators such as the Hansen-Sengupta operator exist and would lead to similar results (see e.g.[START_REF] Goldzstejn | Sensitivity analysis using a fixed point interval iteration[END_REF] for a comparison between the Krawczyk and Hansen-Sengupta operators in the context of sensitivity analysis). We have chosen the Krawczyk operator because it is straightforward to implement using e.g. the Matlab library Intlab[START_REF] Rump | INTLAB -INTerval LABoratory[END_REF].

The fixed point iteration may converge slowly or diverge because of large uncertainties or because [f x ]([x], [a]) contains some singularities.

(a) Krawczyk method (b) standard linearization method Figure 12: The isocontours of ∆φ max [10 -5 rad]

c φ , c θ , c ψ , s φ , s θ and s ψ stand for cos φ, cos θ, cos ψ, sin φ, sin θ and sin ψ, respectively.

Appendix A. Proof of the main result

We prove here that [K]([x]) ⊆ int[x] implies properties (A) and (B) of Section 2.3. First of all, by Corollary 4.4.7 of [START_REF] Neumaier | Interval Methods for Systems of Equations[END_REF]

) is an H-matrix, and thus by definition of strong regularity that [f x ]([x], [a]) is strongly regular, which proves Property (B).

For an arbitrary a ∈ [a], let us define g(x) := f (x, a). Using the well known properties of interval centered extensions,

which is Cg(x). Furthermore, we have g x (x) = f x (x, a) and therefore

) is an interval extension of the derivatives of g. So far we have proved that

We can finally apply Theorem 5.1.8 of [START_REF] Neumaier | Interval Methods for Systems of Equations[END_REF] which proves that

implies g has an unique solution inside [x]. Since this holds for an arbitrary a ∈ [a], we have proved Property (A).