Certified Solving and Synthesis on Modeling of the Kinematics. Problems of Gough-Type Parallel Manipulators with an Exact Algebraic Method

1. Introduction

The significant advantages of parallel robots over serial manipulators are now well known. However, they still pose a serious challenge when considering their kinematics. This paper covers the state-of-the-art on modeling issues and certified solving of kinematics problems. Parallel manipulator architectures can be divided into two categories: planar and spatial. Firstly, the typical planar parallel manipulator contains three kinematics chains lying on one plane where the resulting end-effector displacements are restricted. The majority of these mechanisms fall into the category of the 3-RPR generic planar manipulator, [Gosselin 1994, Rolland 2006]. Secondly, the typical spatial parallel manipulator is an hexapod constituted by six kinematics chains and a sensor number corresponding to the actuator number, namely the 6-6 general manipulator, fig. 1.

Figure 1.

The general 6-6 hexapod manipulators

Solving the FKP of general parallel manipulators was identified as finding the real roots of a system of non-linear equations with a finite number of complex roots. For the 3-RPR, 8 assembly modes were first counted, [Primerose and Freudenstein 1969]. Hunt geometrically demonstrated that the 3-RPR could yield 6 assembly modes, [Hunt 1983]. The numeric iteration methods such as the very popular Newton one were first implemented, [Dieudonne 1972, Merlet 1987, Sugimoto 1987]. They only converge on one real root and the method can even fail to compute it. To compute all the solutions, polynomial equations were justified, [Gosselin and Angeles 1988]. Ronga, Lazard and Mourrain have established that the general 6-6 hexapod FKP has 40 complex solutions using respectively Gröbner bases, Chern classes of vector bundles and explicit elimination techniques, [Ronga and Vust 1992, Lazard 1993, Mourrain 1993a]. The continuation method was then applied to find the solutions, [Raghavan 1993], however, it will be explained why they are prone to miss some solutions, [Rolland 2003]. Computer algebra was then selected in order to manipulate exact intermediate results and solve the issue of numeric instabilities related to round-off errors so common with purely numerical methods. Using variable elimination, for the 3-RPR, 6 complex solutions were calculated [Gosselin 1994] and, for the 6-6, Husty and Wampler applied resultants to solve the FKP with success, [Husty 1996, Wampler 96]. However, resultant or dialytic elimination can add spurious solutions, [Rolland 2003] and it will be demonstrated how these can be hidden in the polynomial leading coefficients. Inasmuch, a sole univariate polynomial cannot be proven equivalent to a complete system of several polynomials. Intervals analyses were also implemented with the Newton method to certify results, [Didrit et al. 1998, Merlet 2004]. However, these methods are often plagued by the usual Jacobian inversion problems and thus cannot guarantee to find solutions in all non-singular instances. The geometric iterative method has shown promises, [Petuya et al. 2005], but, as for any other iterative methods, it needs a proper initial guess.

Hence, this justified the implementation of an exact method based on proven variable elimination leading to an equivalent system preserving original system properties. The proposed method uses Gröbner bases and the rational univariate representation, [Faugère 1999, Rouillier 1999, Rouillier and Zimmermann 2001], implementing specific techniques in the specific context of the FKP, [Rolland 2005]. Three journal articles have been covering this question for the general planar and spatial manipulators [Rolland 2005, Rolland 2006, Rolland 2007]. This algebraic method will be fully detailed in this chapter.

This document is divided into 3 main topics distributed into five sections. The first part describes the kinematics fundamentals and definitions upon which the exact models are built. The second section details the two models for the inverse kinematics problem, addresses the issue of the kinematics modeling aimed at its adequate algebraic resolution. The third section describes the ten formulations for the forward kinematics problem. They are classified into two families: the displacement based models and position based ones. The fourth section gives a brief description of the theoretical information about the selected exact algebraic method. The method implements proven variable elimination and the algorithms compute two important mathematical objects which shall be described: a Gröbner Basis and the Rational Univariate Representation including a univariate equation. In the fifth section, one FKP typical example shall be solved implementing the ten identified kinematics models. Comparing the results, three kinematics models shall be retained. The selected manipulator is a generic 6-6 in a realistic configuration, measured on a real parallel robot prototype constructed from a theoretically singularity-free design. Further computation trials shall be performed on the effective 6-6 and theoretical one to improve response times and result files sizes. Consequently, the effective configuration does not feature the geometric properties specified on the theoretical design. Hence, the FKP of theoretical designs shall be studied and their kinematics results compared and analyzed. Moreover, the posture analysis or assembly mode issue shall be covered.

2. Kinematics of parallel manipulator

2.1. Kinematics notations and hypotheses

Figure 2.

Typical kinematics chains.

The parallel Gough platform, namely 6-6, is constituted by six kinematics chains, fig. 2. It is characterized by its mechanical configuration parameters and the joint variables. The configuration parameters are thus OA_Rf as the base geometry and CB_Rm as the mobile platform geometry. The joint variables are described as ρ the joint actuator positions (angular or linear). Lets assume rigid kinematics chains, a rigid mobile platform, a rigid base and frictionless ball joints between platforms and kinematics chains.

2.2. Hexapod exact modeling

Stringent applications such as milling or surgery require kinematics models as close as possible to exactness. Realistically, any effective configuration always comprises small but significant manufacturing errors, [Vischer 1996, Patel and Ehmann 1997]. Hence, any constructed parallel manipulator never corresponds to the theoretical one where specific geometric properties may have been chosen, for example, to alleviate singularities or to simplify kinematics solving. Two prismatic actuator axes may be neither collinear nor parallel and may not even intersect. Whilst knowing joints prone to many imperfections, then rotation axes are not intersecting and the angles between them are never perpendicular. Moreover, real ball joints differ from a perfectly circular shape and friction induces unforeseeable joint shape modification, which results into unknown axis changes. However, the joint axis angles stay almost perpendicular and any rotation combination shall be feasible. In a similar fashion, the Cardan joint axes are not perpendicular and may be separated by a small offset. Finally, the articulation center is not crossed by any axis.

Identified the hexapode 138, the exact geometric model is then characterized by 138 configuration parameters. Each kinematics chain is described by 23 parameters, as shown on fig. 2 and defined hereafter:

• the 3 parameters of each base joint A_i with their error vector A_i,

• the 3 joint Ai inter-axis distances є₁ ^a, є₂ ^a and є₃ ^a

• _ii
• _r

• the angular deviation between the two prismatic actuator axes: φ,

• the 3 parameters of the platform joint B_i with their error vector B_i,

• the 3 joint B_i inter-axis distances and є₁ ^b, є₂ ^b and є₃ ^b

To solve this model includes the determination of parameters which cannot be measured neither determined. Moreover, the model includes more variables than equations and therefore, its resolution would then only be possible through optimization methods. Relying on a calibration procedure would only determine configuration parameters by specifying an error margin consisting of a radius around joint positions and would not indicate the direction of the error vector. Hence, only an error ball becomes applicable to the model. In practice, the A_i and B_i joint error vectors shall reposition the respective kinematics chains by adding an offset to the joint centers. Thus, a random function shall compute the A_i and B_i vectors with the maximum being the error ball radius. Finally, the selected model, namely the hexapod 84, is effectively based on the hexapod 42 model with errors added to the configuration data and joint variables.

2.3. Kinematics problems

Definition 2.1 The kinematics model is an implicit relation between the configuration parameters and the posture variables, F(X, , OA_|Rf,CB_|Rm)=0 where = {₁,₂,…, ₆}.

Figure 3.

Kinematics model.

Three problems can be derived from the above relation: the forward kinematics problem (FKP), the inverse kinematics problem (IKP) and the kinematics calibration problem, fig. 3. The two first problems shall be covered in this article. The inverse kinematics problem (IKP) is defined as:

Definition 2.2 Given the generalized coordinates of the manipulator end-effector, find the joint positions.

The 6-6 IKP yields explicit solutions from vector = G(X, OA_{|R f}, CB_|Rm ) and is used to prepare the FKP which is defined as:

Definition 2.3 Given the joint positions , find the generalized coordinates X of the manipulator end-effector.

The 6-6 FKP is a difficult problem, [Merlet 1994, Raghavan and Roth 1995] and explicit solutions X = G(, OA_{|R f}, CB_|Rm) have not yet been established. The difficulties in solving the FKP have hampered the application of parallel robot in the milling industry.

2.4. Vectorial formulation of the basic kinematics model

Figure 4.

The vectorial formulation.

The vectorial formulation produces an equation system which contains the same number of equations as the number of variables, fig. (4), [Dieudonne et al. 1972]. A closed vector cycle is constituted between the manipulator characteristic points: A_i and B_i, kinematics chain attachment points, O the fixed base reference frame and C the mobile platform reference frame. For each kinematics chain, a function between points A_i and B_i expresses the generalized coordinates X, such as A_iB_i = U₁(X). Inasmuch, vector A_iB_i is determined with the joint coordinates and X giving a function U₂(X, ). Finally, the following equality has to be solved: U₁(X) = U₂(X, ).

3. The inverse kinematics problem

For each kinematics chain, i = 1,..., 6, each platform point OB_i|Rf can be expressed in terms of the distance constraint, [Merlet 1997]:

$${{\displaystyle l}}_i^2=\Vert {{\displaystyle A}}_i{{{\displaystyle B}}_i\Vert }^2,i=\mathrm{1...6}$$ E1

Using the vectorial formulation, two equation families can be derived: displacement-based and position-based equations.

3.1. Displacement based equations

Any mobile platform position OB_|Rf which meets constraints 1 has a rotation matrix such that:

$${\overrightarrow{OB}}_{i|R_f}={\overrightarrow{OC}}_{i|R_f}+\Re \cdot {\overrightarrow{CB}}_{i|R_m},i=1\dots 6$$ E2

Substituting 2 in 1, we obtain:

$$l_i^2=\Vert {\overrightarrow{OC}}_{|R_f}+\Re \cdot {\overrightarrow{CB}}_{i|R_m}-{{\overrightarrow{OA}}_{i|R_f}\Vert }^2,i=1\dots 6$$ E3

This last equation system can be developed and simplified, leading to the IKP :

$$l_i^2={\left({\overrightarrow{OC}}_{|R_f}-{\overrightarrow{OA}}_{i|R_f}\right)}^2+\left({\overrightarrow{OC}}_{|R_f}-{\overrightarrow{OA}}_{i|R_f}\right)\Re \cdot {\overrightarrow{CB}}_{i|R_m}+{\overrightarrow{CB}}_i^2$$ E4

3.2. Position based equations

In 3D space, any rigid body can be positioned by 3 of its distinct non-colinear points, [Fischer and Daniel 1992, Lazard 1992b]. The 3 mobile platform distinct points are usually selected as the 3 joint centers B₁, B₂, B₃, fig. 5. The 6 variables are set as: OB_i|Rf = [x_i, y_i, z_i] for i = 1... 3. The OB_i|Rf parameters define the reference frame R_b1 relative to the mobile platform and B₁ is chosen as its center. The frame axes u₁, u₂ and u₃ are determined by the 3 platform points:

$$u_1=\frac{{\overrightarrow{B_{1}B}}_2}{\Vert \overrightarrow{B_1{B}_2}\Vert },u_2=\frac{\overrightarrow{B_1{B}_3}}{\Vert \overrightarrow{B_1{B}_3}\Vert },u_3=u_1\wedge u_2$$ E5

Any platform point M can be expressed by B₁M = a_Mu₁ + b_Mu₂ + c_Mu₃ where a_M, b_M, c_M are constants in terms of these three points. Hence, in the case of the IKP, the constants are noted a_Bi, b_Bi, c_Bi, i = i... 6 and can explicitly be deduced from CB_|Rm by solving the following linear system of equations:

$${\overrightarrow{B_{1}B}}_{i|R_{b_1}}=a_{B_i}{u}_1+b_{B_i}{u}_2+c_{B_i}{u}_3,i=1\dots 6$$ E6

Figure 5.

The platform three point coordinate system

Substituting relations 6 in the distance equations l_i ² = ║A_iB_i|Rf║, i = 1... 6, the system can be expressed with respect to the variables x_i, y_i, z_i, i = 1, 2, 3. Thus, for i = 1... 6, the IKP is obtained by isolating the ρ_i or l_i linear actuator variables in the six following equations:

$${{\displaystyle l}}_i^2={\left({{\displaystyle x}}_i-{{\displaystyle \overrightarrow{OA}}}_{ix}\right)}^2+{\left({{\displaystyle y}}_i-{{\displaystyle \overrightarrow{OA}}}_{iy}\right)}^2,i=\mathrm{1...3}$$ E7

$$l_i^2={\Vert {\overrightarrow{B}}_{k|R_{b_1}}-{\overrightarrow{OA}}_{k|R_f}\Vert }^2,i=4\dots 6$$ E8

4. The forward kinematics problem

4.1. Displacement based equations

There exist various formulations of the displacement based equation models.

4.1.1. AFD1 - formulation with the position and the trigonometry identity

The AFD1 formulation is obtained by replacing each trigonometric function of the IKP rotation matrix, 2, by one distinct variable, [Merlet 1987], for j = 1, 2, 3, then c_j = cos(Θj), s_j = sin(Θ_j). The end-effector position variables are retained. The 9 unknowns are then: {x_c, y_c, z_c, c₁, c₂, c₃, s₁, s₂, s₃}. The orientation variables can either be any Euler angles or the navigation ones (pitch, yaw and roll). The orientation variables are linked by the 3 trigonometric identities, for j = 1... 3, then c² _j + s² _j = 1 which complete the equation system:

$$F_i={\left({\overrightarrow{OC}}_{|R_f}-{\overrightarrow{OA}}_{i|R_f}\right)}^2+\left({\overrightarrow{OC}}_{|R_f}-{\overrightarrow{OA}}_{i|R_f}\right)\Re \cdot {\overrightarrow{CB}}_{i|R_m}-L_i^2,i=\mathrm{1,}\dots \mathrm{,6}$$ E9

$$F_j=c_j^2+s_j^2-\mathrm{1,}j=\mathrm{1,2,3}$$ E10

The system is constituted of 9 equations with 6 polynomials of degree 6 and 3 quadratics. The model is simply build by variable substitution without any computation. Thus, the coefficients remain unchanged. The number of variables is not minimal.

4.1.2. AFD2 - formulation with the position and the trigonometric function change

The end-effector position variables are retained. Rotation variable changes can apply the following trigonometric relations, [Griffis and Duffy 1989, Parenti­Castelli and C. Innocenti 1990, Lazard 1993]. For i = 1, 2, 3:

$$t_i=\tan \left(\frac{{\theta }_i}{2}\right),\sin \left({\theta }_i\right)=\frac{2t_i}{\left(1+t_i^2\right)},\cos \left({\theta }_i\right)=\frac{\left(1-t_i^2\right)}{\left(1+t_i^2\right)}$$ E11

The 6 variables become {x_c, y_c, z_c, t₁, t₂, t₃}. The IKP equations (2) are rewritten to obtain the 6 following equations:

$$F_i={\left({\overrightarrow{OC}}_{|R_f}-{\overrightarrow{OA}}_{i|R_f}\right)}^2+\left({\overrightarrow{OC}}_{|R_f}-{\overrightarrow{OA}}_{i|R_f}\right)\Re \cdot {\overrightarrow{CB}}_{i|R_m}-L_i^2,i=\mathrm{1,}\dots \mathrm{,6}$$ E12

The final equation system comprises 6 equations of order 8 with the high degree monomial being x_i ²x_j ²x_k ²x_n ². This model has a minimal variable number. The polynomials coefficients are expanding due to variable change computation. Moreover, this representation is not intuitive.

4.1.3. AFD3 - formulation with the translation and rotation matrix

The intuitive way to set an algebraic equation system from the IKP equations 2 is to straightforwardly use all the rotation matrix parameters and the vector OC_|Rf coordinates as unknowns, [Lazard 1993, Sreenivasan et al. 1994, Bruyninckx and DeSchutter 1996]. The variables are then {X_c, Y_c, Z_c, r_ij, j=1...3, i=1...3}. Since is a rotation matrix, the following relations hold: ^t = Id or det() = 1. These relations are redundant since ^t is symmetrical and they generate the 7 following equations:

$$\{\begin{array}{l}1=r_{11}^2+r_{12}^2+r_{13}^2\mathrm{,1}=r_{21}^2+r_{22}^2+r_{23}^2\mathrm{,1}=r_{31}^2+r_{32}^2+r_{33}^2\\ 0=r_{11}{r}_{21}+r_{12}{r}_{22}+r_{13}{r}_{23}\mathrm{,0}=r_{11}{r}_{31}+r_{12}{r}_{32}+r_{13}{r}_{33}\mathrm{,0}=r_{21}{r}_{31}+r_{22}{r}_{32}+r_{23}{r}_{33}\\ 1=r_{11}{r}_{22}{r}_{33}-r_{11}{r}_{23}{r}_{32}-r_{21}{r}_{12}{r}_{33}+r_{21}{r}_{13}{r}_{32}+r_{31}{r}_{12}{r}_{23}-r_{31}{r}_{13}{r}_{22}\end{array}$$ E13

Six rotation matrix constraints are then selected and preferably with the lowest degree polynomials. This leads to an algebraic system with 12 polynomial equations (13 and 1) in 12 unknowns.

$$F_i={\left({\overrightarrow{OC}}_{|R_f}-{\overrightarrow{OA}}_{i|R_f}\right)}^2+\left({\overrightarrow{OC}}_{|R_f}-{\overrightarrow{OA}}_{i|R_f}\right)\Re \cdot {\overrightarrow{CB}}_{i|R_m}-L_i^2,i=\mathrm{1,}\dots \mathrm{,6}$$ E14

$$F_7=r_{11}^2+r_{12}^2+r_{13}^2-1$$ E15

$$F_8=r_{21}^2+r_{22}^2+r_{23}^2-1$$ E16

$$F_9=r_{31}^2+r_{32}^2+r_{33}^2-1$$ E17

$$F_{10}=r_{11}{r}_{21}+r_{12}{r}_{22}+r_{13}{r}_{23}$$ E18

$$F_{11}=r_{11}{r}_{31}+r_{12}{r}_{32}+r_{13}{r}_{33}$$ E19

$$F_{12}=r_{21}{r}_{31}+r_{22}{r}_{32}+r_{23}{r}_{33}$$ E20

Finally, the model polynomials are quadratic and minimal. They are obtained by substitution and no computations are required. The coefficients are then unchanged. There is a very large number of variables.

4.1.4. AFD4 - formulation with the translation and Gröbner Basis on the rotation matrix

The rotation matrix constraints are not depending on the end-effector position variables. Hence, if one Gröbner Basis is computed from the rotation constraints, the Gröbner Basis is also independent of the position variables and thus constant for any FKP pose. Therefore, one preliminary Gröbner Basis can be calculated and saved into a file for later reuse.

Hence, the rotation matrix constraints in the system 20 can be replaced by their Gröbner Basis comprising 24 equations where the coefficients are only unity. Thus, the algebraic system involves 30 equations and 12 variables.

4.1.5. AFD5 - translation and quaternion algebraic model

Based on equation (2), quaternions can express mobile platform rotation, [Lazard 1993, Mourrain 1993b, Egner 1996, Murray et al. 1997]. The quaternion representation includes 4 variables q₀; q₁; q₂; q₃ where the vector q = q₁ i +q₂ j +q₃ k defines the platform specific rotation axis and q₀ = cos(/2) determines the coordinate expressing the rotation along that axis. Thus, the rotation matrix used in equations 4 may then be expressed in terms of the quaternion coordinates and with ² = q₀ ² + q₁ ² + q₂ ² + q₃ ², we can write:

$$\Re ={{\displaystyle \Delta }}^{-2}\left(\begin{array}{ccc}{q}_0^2+q_1^2-q_2^2-q_3^2& 2\left(q_1{q}_2-q_0{q}_3\right)& 2\left(q_1{q}_3+q_0{q}_2\right)\\ 2\left(q_1{q}_2+q_0{q}_3\right)& q_0^2-q_1^2+q_2^2-q_3^2& 2\left(q_2{q}_3-q_0{q}_1\right)\\ 2\left(q_1{q}_3-q_0{q}_2\right)& 2\left(q_2{q}_3+q_0{q}_1\right)& q_0^2-q_1^2-q_2^2+q_3^2\end{array}\right)$$ E21

The end-effector position variables are retained. Moreover, one may implement a unitary quaternion: ² = 1. Rewriting the IKP equations 4, we obtain 7 polynomial equations in the 7 unknowns X_c; Y_c; Z_c; q₀; q₁; q₂; q₃:

$$F_i={\left({\overrightarrow{OC}}_{|R_f}-{\overrightarrow{OA}}_{i|R_f}\right)}^2+\left({\overrightarrow{OC}}_{|R_f}-{\overrightarrow{OA}}_{i|R_f}\right)\Re \cdot {\overrightarrow{CB}}_{i|R_m}-L_i^2,i=\mathrm{1,}\dots \mathrm{,6}$$ E22

$$F_7=q_0^2+q_1^2+q_2^2+q_3^2-1$$ E23

The system contains 6 polynomials of degree 6 and 1 quadratic. The highest degree monomial is x_i ² x_j ². The quaternion has intrinsic coordinate redundancy which allows avoiding typical mathematical singularities seen in other representations. The number of variable is almost minimal. The rotation matrix system must be recomputed leading to larger resulting polynomial coefficients.

4.1.6. AFD6 - translation and dual quaternion algebraic model

Not only orientations can be formulated using quaternions, but also positions, [Husty 1996, Wampler 96]. The rotation matrix is then expressed in terms of the first quaternion = c₀; c₁; c₂; c₃. In a sense, the second = g₁, g₂, g₃, g₄ represents the end-effector position. Moreover, one relation can be written between the two quaternions: = OC . This relation unfolds in the following equations from which two constraint equations, noted FC₁ = 0 and FC₂ = 0, are selected. Lets s_i = OA_|Rf and t_i = CB_|Rm, then:

$$\begin{array}{l}{\overline{c}}^t\overline{c}=1\\ {\overline{g}}^t\overline{c}=0\\ {\overline{g}}^t\overline{g}-l_1{\overline{c}}^t\overline{c}=0\\ {\overline{c}}^t{s}_i\overline{c}+2{\overline{g}}^t{t}_i\overline{c}=0\end{array}$$ For i = 2,…,6

The dual quaternion system is thus constituted by the 8 following equations, for i = 1 … 6 :

$$F_i={\left({\overrightarrow{OC}}_{|R_f}-{\overrightarrow{OA}}_{i|R_f}\right)}^2+\left({\overrightarrow{OC}}_{|R_f}-{\overrightarrow{OA}}_{i|R_f}\right)\Re \cdot {\overrightarrow{CB}}_{i|R_m}-L_i^2$$ E24

$$F_7=F{C}_1$$ E25

$$F_8=F{C}_2$$ E26

The system comprises 6 polynomials of degree 4 and 2 quadratics. The highest degree monomials are either x_i ⁴; x_i ³ x_j or x_i ² x_j ². One more variable is added over the former quaternion model. The variable choice is not intuitive.

4.2. Position based equations

We shall examine four formulations derived from the position based equations. Every variable has the same units and their range is equivalent.

4.2.1. AFP1 - three point model with platform dimensional constraints

The 3 platform distinct points are usually selected as the three joint centers B₁, B₂ and B₃, fig. 5. The 6 variables are set as: OB_i|Rf = [x_i, y_i, z_i] for i = 1 …3.

Using the relations 6, the constraint equations L_i ² = ║A_iB_i|Rf║², i = 1, …, 6 can be expressed with respect to the variables x_i, y_i, z_i, i = 1, 2, 3. Together with equations 30, they define an algebraic system with 9 equations in 9 unknowns {x₁, y_1, z₁, x₂, y₂, z₂, x₃, y₃, z₃}. The resulting kinematics chain system becomes:

$$F_i={\left(x_i-O{A}_{ix}\right)}^2+{\left(y_i-O{A}_{iy}\right)}^2+{\left(x_i-O{A}_{ix}\right)}^2-L_i^2,i=\mathrm{1...3}$$ E27

$$F_j={\Vert {\overrightarrow{B}}_{j|R_{b_1}}-{\overrightarrow{OA}}_{j|R_f}\Vert }^2-L_j^2,j=4\dots 6$$ E28

The mobile platform geometry yields the following three distance equations:

$$\begin{array}{l}{F}_7={\Vert {\overrightarrow{B_{2}B}}_{1|R_f}\Vert }^2-{\left(x_2-x_1\right)}^2+{\left(y_2-y_1\right)}^2+{\left(z_2-z_1\right)}^2={\Vert {\overrightarrow{B_{2}B}}_{1|R_m}\Vert }^2\\ F_8={\Vert {\overrightarrow{B_{3}B}}_{1|R_f}\Vert }^2-{\left(x_3-x_1\right)}^2+{\left(y_3-y_1\right)}^2+{\left(z_3-z_1\right)}^2={\Vert {\overrightarrow{B_{3}B}}_{1|R_m}\Vert }^2\\ F_9={\Vert {\overrightarrow{B_{3}B}}_{2|R_f}\Vert }^2-{\left(x_3-x_2\right)}^2+{\left(y_3-y_2\right)}^2+{\left(z_3-z_2\right)}^2={\Vert {\overrightarrow{B_{3}B}}_{2|R_m}\Vert }^2\end{array}$$ E29

Together with equations 30, they produce an algebraic system with 9 equations with 9 unknowns {x₁, y₁, z₁, x₂, y₂, z₂, x₃, y₃, z₃}. In all instances, it can be easily proven that this 6-6 FKP formulation yields 9 quadratic polynomials.

The system variable choice is relatively intuitive. Each equation polynomial is always quadratic. However, the b₁ reference frame and the platform points B_i in the b₁ frame require computations, which usually result into coefficient size explosion. The variable number is not minimal.

4.2.2. AFP2 - the three point model with platform constraints

The former system can be slightly modified by replacing the last mobile platform constraint with a platform normal vector one. Hence, lets take the two mobile platform vectors B₁B₂ and B₁B₃, then the last constraint is calculated from these two vector multiplication:

$$\begin{array}{l}{F}_7={\Vert {\overrightarrow{B_{2}B}}_{1|R_f}\Vert }^2-{\left(x_2-x_1\right)}^2+{\left(y_2-y_1\right)}^2+{\left(z_2-z_1\right)}^2={\Vert {\overrightarrow{B_{2}B}}_{1|R_m}\Vert }^2\\ F_8={\Vert {\overrightarrow{B_{3}B}}_{1|R_f}\Vert }^2-{\left(x_3-x_1\right)}^2+{\left(y_3-y_1\right)}^2+{\left(z_3-z_1\right)}^2={\Vert {\overrightarrow{B_{3}B}}_{1|R_m}\Vert }^2\\ F_9=\left(x_3-x_1\right)\ast \left(x_2-x_1\right)+\left(y_3-y_1\right)\ast \left(y_2-y_1\right)+\left(z_3-z_1\right)\ast \left(z_2-z_1\right)-\Vert {\overrightarrow{B_{3}B}}_{2|R_m}\Vert \wedge \Vert {\overrightarrow{B_{3}B}}_{1|R_m}\Vert \end{array}$$ E30

The result is still an algebraic system with nine equations in the former nine unknowns

{x₁, y₁, z₁, x₂, y₂, z₂, x₃, y₃, z₃}. The 6-6 FKP formulation using this three point model is constituted by nine quadratic polynomials.

4.2.3. AFP3 - the three point model with constraints and function recombination

By rewriting the IKP as functions, the algebraic system comprises three equations and three functions in terms of the nine variables: x₁, y₁, z₁, x₂, y₂, z₂, x₃, y₃, z₃, equation (29).

$$F_i={\left(x_i-O{A}_{ix}\right)}^2+{\left(y_i-O{A}_{iy}\right)}^2-l_i^2,i=\mathrm{1...3}$$ E31

$$C_i={\Vert {\overrightarrow{B}}_{k|R_{b_1}}-{\overrightarrow{OA}}_{k|R_f}\Vert }^2-l_i^2,i=4\dots 6$$ E32

Hence, three constraints are derived from the following three functions, [Faugère and Lazard 1995]. Two functions can be written using two characteristic platform vector norms between the B₁,B₂ distinct points and the B₁,B₃ ones. The last function comes from these vector multiplication.

$$\begin{array}{l}{F}_7={\Vert {\overrightarrow{B_{2}B}}_{1|R_f}\Vert }^2-{\left(x_2-x_1\right)}^2+{\left(y_2-y_1\right)}^2+{\left(z_2-z_1\right)}^2={\Vert {\overrightarrow{B_{2}B}}_{1|R_m}\Vert }^2\\ F_8={\Vert {\overrightarrow{B_{3}B}}_{1|R_f}\Vert }^2-{\left(x_3-x_1\right)}^2+{\left(y_3-y_1\right)}^2+{\left(z_3-z_1\right)}^2={\Vert {\overrightarrow{B_{3}B}}_{1|R_m}\Vert }^2\\ F_9=\left(x_3-x_1\right)\ast \left(x_2-x_1\right)+\left(y_3-y_1\right)\ast \left(y_2-y_1\right)+\left(z_3-z_1\right)\ast \left(z_2-z_1\right)-\Vert {\overrightarrow{B_{3}B}}_{2|R_m}\Vert \wedge \Vert {\overrightarrow{B_{3}B}}_{1|R_m}\Vert \end{array}$$ E33

Furthermore, the three last equations (F₇, F₈, F₉) are computed by the following function sequential combinations:

$$\begin{array}{l}{F}_{\text{7}}=-C_{\text{7}}+F_{\text{1}}+F_{\text{2}}\\ F_{\text{8}}=-C_{\text{8}}+F_{\text{1}}+F_{\text{3} }\\ F_{\text{9}}=\text{ 2}*C_{\text{9}}+F_{\text{7}}+F_{\text{8}}-\text{2}*F_{\text{1}}\end{array}$$ E34

The formulation is completed with other function combinations obtained by the following algorithm leading to three middle equations (F_4, F₅, F₆). Let d₇ = ║B₂B_1|Rmj║, d₈ = ║B₃B_1|Rm║² and d₉ = ║B₃B₂|Rm║ ^ ║B₃B_1|Rm║, then for i = 4, 5, 6, we compute:

$$\begin{array}{l}{C}_i=C_i-a_{B_i}^2\ast C_7-b_{B_i}^2\ast C_8-c_{B_i}^2\ast \left(C_7\ast C_8-C_9^2\right)-a_{B_i}\ast b_{B_i}\ast \left(2\ast C_9\right)\\ C_i=2\ast C_i-a_{B_i}\ast \left(F_7-2\ast F_1\right)-b_{B_i}\ast \left(F_1+F_2-F_7\right)\\ F_i=C_i-2\ast c_{B_i}^2\ast d_7\ast \left(F_1+F_2-F_8\right)-2\ast c_{B_i}^2\ast d_8\ast \left(F_1+F_2-F_7\right)\\ -2\ast c_{B_i}^2\ast d_9\ast \left(F_7-F_9+F_8-2\ast F_1\right)+2\ast \left(a_{B_i}+b_B{}_i-1\right)\ast F_1-F_7\ast a_{B_i}-F_8\ast b_{B_i}\end{array}$$ E35

The result is an algebraic system with nine equations with the nine unknowns. The 6-6 FKP

formulation using this modified three point model includes six quadratic and three quartic polynomials. The system includes polynomials of higher degree than for the former two position based models. Computations cause to coefficient expansion.

4.2.4. AFP4 - the six point model

The six mobile platform B_i joints can be used in defining 18 variables, [Rolland 2003]. Taking the IKP equations (8), a position based variation is obtained:

$$l_i^2={\left(x_i-O{A}_{ix}\right)}^2+{\left(y_i-O{A}_{iy}\right)}^2+{\left(x_i-O{A}_{ix}\right)}^2,i=\mathrm{1...6}$$ E36

The system is completed with 12 distance constraint equations selected among the distinct B_i passive platform joints. Here are some examples:

$$\begin{array}{l}{\Vert {\overrightarrow{B_{i}B}}_{1|R_f}\Vert }^2={\left(x_i-x_1\right)}^2+{\left(y_i-y_1\right)}^2+{\left(z_i-z_1\right)}^2={\Vert {\overrightarrow{B_{i}B}}_{1|R_m}\Vert }^2,i=1\dots 6\\ {\Vert {\overrightarrow{B_{j}B}}_{2|R_f}\Vert }^2={\left(x_j-x_2\right)}^2+{\left(y_j-y_2\right)}^2+{\left(z_j-z_2\right)}^2={\Vert {\overrightarrow{B_{j}B}}_{2|R_m}\Vert }^2,j=3\dots 6\\ {\Vert {\overrightarrow{B_{k}B}}_{3|R_f}\Vert }^2={\left(x_k-x_3\right)}^2+{\left(y_k-y_3\right)}^2+{\left(z_k-z_3\right)}^2={\Vert {\overrightarrow{B_{k}B}}_{3|R_m}\Vert }^2,k=4\dots 6\end{array}$$ E37

The formulation results in 18 polynomials in the 18 unknowns:

{x₁, y₁, z₁, x₂, y₂, z₂, x₃, y₃, z₃, x₄, y₄, z₄, x₅, y₅, z₅, x₆, y₆, z₆}. The system is then constituted of

quadratic polynomials. This variable choice is intuitive and the system yields minimal degree. Finally, the number of variables is maximal.

5. Solving polynomial systems using exact computation

5.1. Mathematical system solving

Kinematics problems contain systems of several equations containing non-linear functions with various variable numbers. These systems can be difficult to solve, especially in the general 6-6 cases and response times actually makes them inappropriate for implementations in design, simulation or control. In some instances, the results may appear to be faulty bringing doubts to the reliability of the methods.

If left without any reliable and performing methods, the tendency, in engineering practice, would be to convert the difficult models into simpler linearized ones. In material handling, this proposal might suffice, however, in high speed milling where the accuracy requirements are more severe, any simplification can have a dramatic impact, whereby result certification becomes an important issue.

However, with proper polynomial formulation, algebraic methods can lead to at least certified and even exact results, whereas numeric methods, unless they implement proper interval analysis, cannot actually obtain certified results since they are prone to numeric instabilities or matrix inversion problems. Therefore, although time consuming, algebraic methods are preferred since they handle integer, rational and symbolic values as such without any truncation or approximation, even when manipulating intermediate results. Hence, there will be no loss of information.

Solving non-linear equation systems will usually result in several complex solutions, out of which a certain subset are real solutions. However, only the real solutions bear practical significance, since they correspond to effective manipulator poses.

5.2. Calculation accuracies

The calculation accuracies are depending upon the type of arithmetic, the behavior of the calculation methods and the quality of the implemented algorithms.

Definition 5.1 An exact calculation is defined as a calculation which always produces the same exact result to the same specific mathematical problem.

The result does not contain any error. Its representation is also exact.

Definition 5.2 A reliable computation is defined as one which will always give the same result from the same initial input data presented in the same format.

Definition 5.3 A certified calculation is defined as a reliable computation giving a result distant from the true solution by a certain maximum known accuracy.

Hence, such a calculation may not be exact. However, the result contains some exact digits. Hence, we shall try to apply a method that computes certified results and if possible exact ones.

For example, lets take the univariate function f₁(x) = x₂ 4/25. Computing f₁ = 0, we obtain the exact response: {2/5, 2/5}. The closed-form resolution calculates exact results with rational numbers. Therefore, the result is certified without any error.

Lets consider f₂(x) = x₂ 5. Solving f₂ = 0, the result will be two irrational numbers which can only be represented by truncation. However an interval can be certified to contain the exact result: {[2, 5/2], [5/2, 2]}. Wherefore, exact computations keep intermediate results in symbolic format whenever possible and only revert to rational or floating boundary numbers for display purposes.

Therefore, any real number can be coded by an interval which width corresponds to the required accuracy. However, the difficulty lies in insuring that the interval contains the exact result which is not known a priori.

Figure 6.

Bloc Diagram of the Continuation Method

5.3. Solving a non-linear system

Two method groups have been advocated to find all solutions of the FKP, namely: continuation methods and variable elimination ones, [Raghavan and Roth 1995]. The first approach is usually realized in a numeric environment and the later algebraic.

5.4. Continuation method with homothopy

In order to compute several solutions, the continuation method can be implemented with a homothopy process. The Continuation approach implements a numerical iterative method which is successively repeated in order to progressively transfer from an original equation system which solutions are predetermined to another system relatively close to the former, 6.

Let a system of equation be F(X) = 0 with variables X = {x₁,…,x_n}; we wish to find the solution to this equation system. Let G(X) = 0 be a similar equation system which roots are already known, namely the variety Vr (I)^G; then, we set the continuation process as H(X, ) = G(X)+ (F(X)G(X)) and commence with = 0. It provides for a mechanism to convert an original equation system into a final one through several steps. At each step, H(X, ) is successively computed with a new value of which is increased by a small arbitrary increment such that {0,…, 1}. The homothopy principle assumes connectivity between solutions of each system computed with the various . More generally, if the system H(X, ) with as a variable would be solved, it would result in paths as solutions.

Figure 7.

Examples of path following with the Continuation Method

The continuation method cannot solve any equation system as such and, at each step, when is instanciated, the H(X, ) system roots are computed by a typical iterative method, either the Newton one or the new Geometric Iterative Method which is a potentially good alternative, [Petuya et al. 2005].

This method was first applied to classical robotics kinematics, [Tsai and Morgan 1984] and then applied to solve the parallel manipulator FKP, [Kholi et al. 1992, Sreenivasan and P. Nanua 1992].

Advocating that a little change on parameters of one system shall cause only a small change on solutions, the continuation method could be used to find the 40 solutions on some 6-6 FKP cases, [Raghavan 1993].

It is feasible to construct an efficient and reliable method; however, the method is still unproven. Moreover, continuation does not alleviate the problems related to the application of a numeric iterative method.

This method can be somewhat delicate to implement. There exist several scenarios which might pose significant problems depending on how the solution paths evolve from = 0 to = 1, see fig. 7, where solutions:

1. go to or come from infinity,

2. merge or split,

3. start complex and become real,

4. start real and become complex.

Therefore, proper implementation would require a priori continuation process verification which is still an open question, since this would require solving a one-dimensional system H(X,) where is left as a variable and this is an even more difficult problem. Inasmuch, in many instances, finding a nearby equation system with known roots may not be always be feasible. Then, there is also an issue on what constitutes a sufficiently similar system.

However, it is very difficult to determine precisely what the meaning of sufficiently similar is.

5.5 Variable elimination

5.5.1. Introduction

Most algebraic methods which were implemented to solve the parallel manipulator FKP apply one form of variable elimination. Let an algebraic system F(X) = 0 be a system of polynomial functions f_i(X), i = 1,…,m with variables X = {x₁,…,x_n}, the variable elimination approach consists in the transformation of the original system F(X) = 0 into another system H(Y) = 0 with functions g_j(X), j = 1,…,p with variables Y = {y₁,…,y_r} where r < n. Ultimately, the goal is to find a method which allows to compute an equation system H(Y) in either triangular format or preferably in univariate form which would be the easiest to solve.

Most variable elimination methods are usually divided into four steps:

1. Step 1: Variable elimination.

2. Step 2: Solving the univariate equation.

3. Step 3: Return or extension to original system variables.

We will examine the variable elimination methods which were successfully applied to solve the FKP from which two can be identified:

1. method based on resultant calculation including the so-called dialytic elimination,

2. method based on Gröbner basis calculation.

5.5.2. Resultant method

Variable elimination can be implemented through a recursive method based on resultants. As input, we give a system of equations with rational coefficients. The output will be one univariate polynomial equation in terms of one of the original variables. Each elimination step involves two polynomial equations which results in one equation with the number of variable reduced by one.

Definition 5.4 [Cox et al. 1992] Let a system be F(X) = {f1,…,f_n} Q[x₁, …, x_n]; let P = f_i and

R = f_j where f_i = a_px₁ ^p +…+a₀ and f_j = b_qx₁ ^q +…+b_q with i, j {1, 2,…,n} and a_i, b_i Q[x₂,…,x_n], let p = deg(P) et r = deg(R) knowing that p, r N* ; suppose that a₀ 0 and b₀ 0 ; then Res(f, g, x₁) = det(M) which is the resultant of P and R in terms of x₁ where M is the identified as the Sylvester matrix.

Then, the Sylvester matrix can be expressed in terms of the polynomial coefficients:

$$M=\left(\begin{array}{cccccccc}{a}_0& 0& \cdots & 0& b_0& 0& \cdots & 0\\ a_1& a_0& \cdots & 0& b_1& b_0& \cdots & 0\\ ⋮& & \ddots & & ⋮& & \ddots & \\ a_l& & \cdots & a_1& b_m& & \cdots & b_1\\ 0& a_l& \cdots & a_2& 0& b_m& \cdots & b_2\\ ⋮& & \ddots & & ⋮& & \ddots & \\ 0& & \cdots & a_l& 0& & \cdots & b_1\end{array}\right)$$ E38

Inasmuch, we can write: Res(P, R, x₁) = det(Sylv(P, R, x₁). If we examine the Sylvester matrix, we can observe that part with the a_i parameters contains m columns and the one with b_i n columns.The following proposition holds and its proof is described in [Cox et al. 1992]: The resultant Res(f, g, x) is the first ideal of elimination < f, g > k[x²,…,x_n]; moreover, Res(f, g, x) = 0 iff f, g have a common factor in k[x₁,…,x_n] which has a positive degree in x. The nature of this factor has to be determined and we wish to establish if it is only one root of functions f and g. To answer that question, the following corollary will be employed: If f, g C[X] then Res(f, g, x) = 0 iff f and g contain a common root in C. This common root is determined by computing det(Sylv(P, R, x₁) = 0. The nature of this root has to be determined, notably if it is a partial one and the answer will come from the following proposition, [Cox et al. 1992] : Knowing that f, g C[X], let a₀, b₀ 0 and a₀, b₀ C[x₂,…, x_n], if Res(f, g, x) C[x₂,…,x_n] cancels at (c₂,…,c_n), then we obtain that either a₀ b₀ = 0 at (c₂,…,c_n) or either c₁ C such as f and g cancel.

In certain instances, the head terms of the polynomials can cancel which will result in the cancellation of the determinant and the process consequently adds one extraneous root.

In order to obtain the univariate equation, a recursive algorithm will be applied. Firstly, we calculate n 1 resultants h_k = Res(f_k+1, f₁, x₁) on variable x₁ for k = 1,…,n 1. Secondly, we compute n 2 resultants h⁽2)_j = Res(h_j+1, h₁, x₂) on variable x₂ for j = 1,…,n 2 and we continue in the same fashion until the univariate equation is determined. An almost triangular equation system is constructed.

$$[\begin{array}{l}{f}_1\left(X\right)=\mathrm{0,}\dots ,f_{n-2}\left(X\right)=\mathrm{0,}{f}_{n-1}\left(X\right)=\mathrm{0,}{f}_n\left(X\right)=0\\ h_1\left(x_2,\dots ,x_n\right),\dots ,h_{n-2}\left(x_2,\dots ,x_n\right),h_{n-1}\left(x_2,\dots ,x_n\right)\\ h_1^{\left(2\right)}\left(x_3,\dots ,x_n\right),\dots ,h_{n-2}^{\left(2\right)}\left(x_3,\dots ,x_n\right)\\ ⋮\\ h_1^{n-2}\left(x_{n-1},x_n\right),h_2^{\left(n-2\right)}\left(x_{n-1},x_n\right)\\ H\left(x_n\right)\end{array}$$ E39

The last H(x_n) is the targeted univariate equation. However, this equation cannot be considered equivalent to the initial algebraic system because the head terms can cancel.

The return step to original variables is performed by substituting back through the triangular system. The equation is solved H(x_n) = 0 and we obtain a series of w roots {x_n}. We take the w roots, one by one, which is introduced in one of the equations h₁ ⁽ⁿ²⁾ (x_n1) = 0 or h₂ ⁽ⁿ²⁾ (x_n1) = 0 and obtain the w roots {x_n1}. We continue until x₁ is isolated.

The 6-6 FKP has been solved applying resultants, [Husty 1996], in a computer algebra environment to avoid intermediate result truncation, since, in a sense, parameter truncation can be envisioned as changing the manipulator configuration.

A variation to the resultant method is called the dialytic elimination. Let the variable set be X = {x₁,…,x_n} of the algebraic system F(X) = 0; then select any variable x_i and set it as the hidden variable, then a monomial vector is constructed around x_i for the system F(X) which is expressed as W = (1, x_i, x_i ²,…). The FKP is rewritten as a linear system in terms of W:

$$A\overline{W}=0$$ E40

where : W 0

Being a generalization of Res(P, Q, x₁) = det(Sylv(P, Q, x₁), it is subjected to the same risks of root addition through the head term cancellation. Dialytic elimination has been implemented to solve the FKP of the 3-RS or MSSM parallel manipulators, [Griffis and Duffy 1989, Dedieu and Norton 1990, Innocenti and Parenti-Castelli 1990 ]. Satisfactory results were produced on simple parallel manipulators, [Raghavan and Roth 1995].

5.6. Gröbner Bases

Lets denote by Q[x₁,…,x_n] the ring of polynomials with rational coefficients. For any n-uple = (₁,…, _n) Nⁿ, lets denote by X the monomial X₁ ¹ ... X_n ⁿ. If < is an admissible monomial ordering and $$P={\displaystyle {\sum }_{i=0}^r{a}_i{X}^{{\mu }^{\left(i\right)}}}$$ any polynomial in Q[X₁,...,X_n], the following polynomial notations are necessary :

-LM(P,<) = max_i=0…r, <X⁽ⁱ⁾ is the leading monomial of P for the order <,

-LC(P,<) = a_i with i such that LT(P) = X⁽ⁱ⁾ is the leading coefficient of P for <,

-LT(P,<) = LC(P,<) LM(P,<) is the leading term of P for <.

Lets denote by x₁,…,x_n the unknowns and S = {P₁ =,…,P_s} any polynomial system as a subset

of Q[x₁,…,x_n]. A point Cⁿ is a zero of S if P_i() = 0 i = 1…s. Actually, any large polynomial equation system cannot be directly or explicitly solved. Thus, it is necessary to revert to mathematical objects containing sufficient information for resolution. Any polynomial system is then described by an ideal:

Definition 5.5 [Cox et al. 1992] An ideal I is defined as the set of all polynomial P(X) that can be constructed by multiplying and adding all polynomials in the ring of polynomials with the original polynomials in the set S.

A Gröbner Basis G is then as a computable polynomial generator set of a selected polynomial set S = {P₁,…,P_s} with good algorithmic properties and defined with respect to a monomial ordering. This basis is a mathematical object including the ideal I information. The lexicographic and degree reverse lexicographic (DRL) orders are usually implemented, [Cox et al. 1992, Geddes et al. 1994]. Given any admissible monomial ordering, the classical Euclidean division can be extended to reduce a polynomial by another one in Q[X₁,…,X_n]. This polynomial reduction can be generalized to the reduction by a polynomial list. The reduction output depends on the monomial ordering < and the polynomial order.

Definition 5.6 Given any admissible monomial ordering, <, a Gröbner Basis G with respect to < of an ideal I Q[X₁,…,X_n] is a finite subset of I such that: f I, g G such that LM(g,<) divides LM(f,<).

Some useful Gröbner Basis properties are described in the following theorem:

Theorem 5.1 Let G be a Gröbner Basis G of an ideal I Q[X₁,…,X_n] for any < monomial ordering, then a polynomial p Q[X₁,…,X_n] belongs to I if and only if the reduction algorithm Reduce(p, G,< ) = 0; the reduction does not depend on the order of the polynomials in the list of G; it can be used as a simplification function.

The classical method for computing a Gröbner Basis is based on Buchberger's algorithm [Buchberger et al. 1982, Buchberger 1985]. Recently, Faugère proposed more powerful algorithms, namely accel and F4 implemented in the FGb software, [Faugère 1999].

5.7. Gröbner Bases and zero dimensional systems

Definition 5.7 [Lazard 1992a] A zero-dimensional system is defined as a mathematical equation system with a variety (solution set) constituted by a finite number of complex points.

From a mathematical point of view, any variety is a valid result. From an engineering one, we seek a variety which represents an exploitable result. For any parallel manipulator FKP, the zero-dimensional systems allow pose analysis, since each real solution corresponds to each manipulator pose. Otherwise the problems fall in the category of rarely usable singular configurations. Detection of zero-dimensional systems is performed by implementing the following theorem:

Theorem 5.2 Let G = {g₁,…,g_l} be a Gröbner Basis for any ordering < of any system S = {P₁,…, P_s} Q[X₁,…,X_n]^s. The two following properties are equivalent:

-For all index i, i = 1…n, there exists a polynomial g_j G and a positive integer n_j such that X_i ^nj = LM(g_j,<);

-The system {P₁ = 0,…,P_s = 0} has a finite number of solutions in Cⁿ.

Hence, one can determine if the FKP resolution yields a finite or infinite number of solutions. This is obviously an important issue, since an algebraic system with an infinite number of solutions (variety of degree one or higher) cannot be directly exploited by a control system and failure is usually the outcome.

5.8. The lexicographic Gröbner Basis

The choice of the monomial ordering is critical for the Gröbner Basis computation efficiency which is evaluated by the computation times, the size and the shape of the result. A lexicographic Gröbner Basis with X₁< X₂< … < X_n as variable order has always the following general shape:

$$\{\begin{array}{l}{f}_1\left(X_1\right)=\mathrm{0,}{f}_2\left(X_1,X_2\right)=\mathrm{0,}\dots \\ f{}_{k_2}\left(X_1,X_2\right)=\mathrm{0,}{f}_{k_2+1}\left(X_1,X_2,X_3\right)=\mathrm{0,}\dots \\ f_{k_{n-1}+1}\left(X_1,\dots ,X_n\right)=\mathrm{0,}\dots \\ f_{k_n}\left(X_1,\dots ,X_n\right)=0\end{array}$$ E41

Whilst computing a lexicographic Gröbner Basis using Buchberger's algorithm yields doubly exponential complexity in the worst case, an order change will save computation time. Therefore, one DRL Gröbner Basis may first be computed in dⁿ polynomial time where d is the maximum total degree of the equations. Then, the DRL Gröbner Basis can be converted into a lexicographic one in O(nD³) arithmetic operations, [Faugère et al. 1991]. In the case of any zero-dimensional system, this can be done using exclusively linear algebra techniques, [Faugère et al. 1991].

Definition 5.8 Let S be a system constituted by the polynomials p1,…,p_s Q[x₁,…,x_n], without square factors. Let G = {g₁,…,g_l} be a Gröbner Basis computed for any ordering < on the system S. If the univariate polynomial f(x₁) is without any multiple factors, then the equation is square-free, the first coordinates of all solutions are distinct and the solutions yield no multiple root. Then, the system is considered in Shape Position.

If the polynomial system is in Shape Position, the Gröbner Basis has the following simplified univariate format:

$$f_{\text{1}}\left(X_{\text{1}}\right)=\mathrm{0,}{X}_{\text{2}}=f_{\text{2}}\left(X_{\text{1}}\right)=\mathrm{0,}\dots ,{\text{ X}}_{\text{n}}=f_{\text{n}}\left(X_{\text{1}}\right)$$ E42

In this format, the lexicographic Gröbner Basis can be computed directly using the F4 algorithm, [Faugère 1999].

5.9. The rational univariate representation

Another root representation is given by the Rational Univariate Representation, hereafter identified Rational Univariate Representation, [Rouillier 1999]. In the particular case of systems in Shape Position, the univariate variable is the first one in the original system ones, so that the Rational Univariate Representation can be written in the following format:

$$\{h\left(X_k\right)=\mathrm{0,}{X}_2=h_2\left(X_k\right)/h_0\left(X_k\right),\dots ,X_n=h_n\left(X_k\right)/h_0\left(X_k\right)$$ E43

Figure 8.

The correspondence between the Rational Univariate Representation and original system varieties.

h(X_k) where 1 k n denotes the univariate equation, h_i i = 0 … n the univariate return functions. These polynomials are elements of Q[X_k]. Inasmuch, h and h₀ are coprime. Moreover, if the Rational Univariate Representation system is in Shape Position and if h(X_k) is square-free, the Rational Univariate Representation can be converted into a lexicographic Gröbner Basis and the Rational Univariate Representation expression can be simply computed from this basis. This situation is generally the case in most FKP instantiations:

$$h\left(X_k\right)=f\left(X_k\right),g\left(X_k\right)=f^{\text{'}}\left(X_k\right)$$ E44

Moreover, since f and f₀ are coprimes, f₀ is invertible modulo f; this leads to: h_i = f_i (f)¹ modulo f, i = 2 … n.

Definition 5.9 Let S be a system constituted by the polynomials p1,…, p_s Q[x₁,…,x_n] and let G = {g₁,…,g_l} be a Gröbner Basis for any ordering < of the system S. Then, the RR form is defined as the system R = {p₁, p₁ ¹p₂ mod G, p₁ ¹p₃ mod G,…, p₁ ¹pk mod G} Q[x₁,…,x_n] where p₁ ¹ is calculated such that (p1 p₁ ¹ modulo G) = 1.

The rationale behind the RR form computation is to reduce the polynomial system S. The computed Gröbner Basis G on any ordering < of the system S is used for the system reduction. Then, the reduced system R resolution shall be less difficult.

Definition 5.10 Let S be a system constituted by the polynomials p₁,…,p_s Q[x₁,…, x_n], if the system S is without square factors and in Shape Position and if a lexicographic Gröbner Basis is computable, then the system RR form R can be deduced.

Since the S system reduction by a lexicographic Gröbner Basis leads to a further reduced system R, one can conjecture faster FKP resolution. However, in some instances, p₁ ¹ can hardly be computed and the RR form cannot be deduced. Hence, the FKP is solved by first computing a Gröbner Basis on a DRL order and then we rely on the aforementioned order change to obtain a lexicographic Gröbner Basis. In the rare instances where the FKP system is not in Shape Position and a lexicographic Gröbner Basis is not computable either by the RR form or by order change, the Rational Univariate Representation can be always be computed by reverting to the more robust and lengthy algorithm taking as input any monomial ordering Gröbner Basis and calculating the multiplication tables and matrices, [Rouillier 1999].

In practice, the Rational Univariate Representation coefficients are smaller than the lexicographic Gröbner Basis ones, [Alonso et al. 1996]. This is the rationale behind the preference of the Rational Univariate Representation expression which can be computed using the algorithm described in [Rouillier 1999]. Anyhow, it should be deduced from a lexicographic Gröbner Basis if the system is in Shape Position. Hence, in any case, the determination of a lexicographic Gröbner Basis is first tried and then, a Rational Univariate Representation is computed. In fact, for any parallel manipulator FKP, the computation of the RR form is preferred. Finally, the method insures that the Rational Univariate Representation system is equivalent to the former polynomial system. Each Rational Univariate Representation system solution corresponds to one and only one original system solution and the properties are preserved, fig. 8.

5.10. Real roots isolation

Once the Rational Univariate Representation is computed, then the real roots of the univariate equation h(X_k) roots can be isolated either using rapid numeric methods or exact algebraic methods. As far as algebra is concerned, such computations can be performed by a method such as the Uspensky one based on Descartes' sign rule and Sturm's theorem. An improved algorithm is introduced in [Rouillier and Zimmermann 2001] applying one Vincent theorem.

In computer algebra and in particular in the proposed method, a natural way for coding a real algebraic number is given by either a square-free polynomial f with rational coefficients such that f() = 0 or an interval [l,r] with rational bounds containing . Inasmuch, if and f() = 0 then [l, r].

Let be a root of h(T) represented by (h, [l, r]), then the resolution of h(T) = 0 produces a solution list {1,…, _r} = [(h, [l_i, r_i]), i = 1 … r] where h is the square-free part of h. Moreover, it is very easy to refine the intervals [l_i, r_i] by dichotomy since h is square-free which induce that sign(h(l_i)) sign(h(r_i)) if l_i r_i.

Since h and h₀ are coprime by definition of the Rational Univariate Representation, one can refine [l, r] so that it doesn't contain any real root of G. One can then simply use interval arithmetic to compute [l⁽ⁱ⁾, r⁽ⁱ⁾] = h_i/h₀ ([l, r]); i = 1 … n. The complete method allows one solution of the equation h(T) = 0 to correspond to each solution of the system F(X) = 0 since the bijection is guaranteed and the properties are preserved, fig. 8.

Then, using the return functions h_i/h₀, for each , we compute an interval product $${\Pi }_{i=0}^n\left[l^{\left(i\right)},r^{\left(i\right)}\right]$$ containing [ h₁/h₀ (),…, h_n/h₀()]. For refining $${\Pi }_{i=0}^n\left[l^{\left(i\right)},r^{\left(i\right)}\right]$$ , it is sufficient to refine and introduce again [l, r] in the coordinate functions hi/h0; i = 1 … n.

The Rational Univariate Representation system real root isolation requires a second step where the original system roots are determinated in terms of the original variables using the return functions provided by the Rational Univariate Representation.

By computing a Rational Univariate Representation and then isolating the real roots as described above, we obtain an exact method for isolating all the real roots of any polynomial system with a finite number of complex roots verified during the computations using theorem 5.2. The Rational Univariate Representation has rational coefficients and the isolating intervals have rational bounds. For each solution = (1,…, _n), the method produces :

- a box B $${\Pi }_{i=0}^n\left[l^{\left(i\right)},r^{\left(i\right)}\right]$$ rⁿ such that B, (li, ri) Q² with exact (rational) values and if is a root of the system, then B;

- a function for refining B up to any precision є > 0 (|ri li| < є, i = 1 … n).

5.11 The certified method description

Assembling the aforementioned algorithms, we obtain a complete certified method since it handles the intermediate results with exact representations. Taking into account the fact that most parallel manipulator FKP are in Shape Position, the method has been simplified, fig. 9.

Figure 9.

Bloc diagram of the exact method.

6. Solving the forward kinematics problem

6.1. Solving the general 6-6

The selected manipulator is a generic 6-6 in a realistic configuration, measured on a real parallel robot prototype, trying to reproduce a singularity-free theoretical SSM design. This example is a typical one among several successful trials which were performed to test the method without any failures. Even when the 6-6 FKP was in a singular pose, the method returned that the solution is a variety of degree one (infinite number of solutions), thus in singularity.

Computations were done on a personal computer equipped with a 400 MHz Pentium II microprocessor including 128 bytes of random access memory. The operating system was Red Hat LINUX. Thus, this conservative approach allows the user to expect better performance.

Figure 10.

The selected 6-6 hexapod manipulator

6.1.1. Modeling the 6-6 using the three point model

The parallel manipulator corresponding to the configuration data is shown in fig. 10. Lets take a typical 6-6 configuration example where construction realism introduce manipulator configuration deformations. It is then written in a configuration text file which includes the manipulator essential parameters: the coordinates of the joint center positions OA|Rf in the fixed base reference frame Rf and the coordinates of the joint center positions CB|Rm in the mobile platform reference frame Rm. In the computations, we use their simplified format, respectively identified as A and B. The unit is the millimeter:

B := [ [ 68410/1000, 393588/1000, 236459/1000 ],

[ 375094/1000, -137623/1000, 236456/1000 ],

[ 306664/1000, -256012/1000, 236461/1000 ],

[ -306664/1000, -255912/1000, 236342/1000 ],

[ -375057/1000, -137509/1000, 236464/1000 ],

[ -68228/1000, 393620/1000, 236400/1000 ]]:

A := [ [ 464141/1000, 389512/1000, -178804/1000 ],

[ 569471/1000, 207131/1000,-178791/1000 ],

[ 1052905/10000,-597151/1000, -178741/1000 ],

[-1052905/10000,-597200/1000, -178601/1000 ],

[-569744/1000, 206972/1000, -178460/1000 ],

[-464454/1000, 389384/1000, -178441/1000 ]]:

Li := [(1250^2)\$i=1..6]:

6.2. The ten formulation comparison

Solving the ten FKP proposed formulations has been tried on the aforementioned example. Performance results are displayed on table 1 where the total real root computation times along with model code and name are given. This table specifies if the method terminated or was interrupted by the user if computations lasted more than one hour.

Five FKP formulations led to computation termination. However, only the AFD4, AFD5, AFP2 and AFP3 models brought response times within an hour. For the AFD1 formulation, only the DRL Gröbner Basis was successfully computed. Therefore, only the three following models were retained:

1. AFD4 - Formulation with the translation and Gröbner Basis of the rotation matrix,

2. AFD5 - Formulation with Translation and Quaternion variables,

3. AFP3 - The three point model with constraints and function recombination.

Code	Method	Resp. time seconds	Ending	Comment
AFD1 AFD2 AFD3 AFD4 AFD5 AFD6 AFP1 AFP2 AFP3 AFP4	trigo. ident. trigo relations GB( ) quaternion somas 3-points+dist 3pts+dist+normal 3pts+modified 6-points	148 3600+ 3600 4,8 59,6 3600 3600+ 125,0 10,9 3600+	FGb yes no yes yes yes no no yes yes no	stopped stopped stopped stopped stopped

Table 1.

Response times for each formulations

Table 2 shows information which can be deduced from system observation. It gives the variable number, equation number and the maximum degree with respect to one variable, file size (ASCII size in bytes) and the maximum number of digits in one coefficient.

Code	Method	Variable number	Equations number	Size bytes
AFD4 AFD5 AFP3	GB on Quaternion 3 pts & recomb	12 7 9	30 7 9	7100 6100 4500

Table 2.

The selected formulation characteristics

6.3. System model analysis and selection

Solving more examples was carried using the three selected formulations in order to evaluate computations. Performance is studied on the Gröbner Basis being the largest computed mathematical object. We compare the resulting file sizes and the computation response times. The computations are performed on the preceding theoretical SSM, identified theo, and on the corresponding real general 6-6, identified eff.

As input, the manipulator configuration data include coefficients which are rational numbers with a six digit numerator and the number 1000 as denominator. Thus, the manipulator is measured at the micron accuracy. Let’s note that computations are proceeding with the RR form.

Robot + model	Strategy tgrob (s)	Accel (s)	F4 (s)	Grobner Basis Size (Kbytes)	Complex Solution	Comment
theo _ AFP3.fgb theo_AFD5.fgb theo_AFD4.fgb	29,06 124,29 28,22	6,53 36,24 3600	0,6 2,45 0,31	130 150 85	36 72 36	stopped
eff_AFP3.fgb eff_AFD5.fgb eff_AFD4.fgb	NA NA NA	1200 NA NA	33,0 31,9 38,9	800 790 768	40 80 40

Table 3.

DRL Gröbner basis computation times and memory usage

As input, the manipulator configuration data include coefficients which are rational numbers with a six digit numerator and the number 1000 as denominator. Thus, the manipulator is measured at the micron accuracy. Let’s note that computations are proceeding with the RR form.

6.4. DRL Gröbner Basis computations

Three algorithms shall be applied, namely tgrob, Accel and F4 respectively the Buchberger, the Gb and the FGb implementations, [Faugère 1999]). Table 3 presents the computations times in seconds and memory space in kilobytes for the various techniques to compute a DRL Gröbner Basis for the three selected formulations.

The response time discrepancy between the SSM and 6-6 is significant. Practically, a small even infinitesimal coefficient difference such as platform non-planarity leads to response times which can be either dramatically longer (15, 55 or 100 times). The Accel algorithm is more sensible to formulation than F4. The Accel algorithm produces lengthy computations using the AFD4 model. The quaternion formulation produces contradicting results, since it leads to significantly longer response times for the SSM, but slightly faster ones for the 6-6. For the 6-6, the response times and file sizes are almost equivalent.

Applying the fast F4 algorithm, the rotation matrix model is preferred since one part of the system equations (13) does not depend on the studied configuration.

6.5. Rational univariate representation computations

The computed DRL Gröbner Basis is provided as the input to the Rational Univariate Representation computation algorithms. Knowing that the Rational Univariate Representation computation behaviour from the three formulation Gröbner bases cannot be predicted and lead to comparable Gröbner Basis sizes, the three Rational Univariate Representation computations have been performed. The following figures come from the direct calculations of the RR form using Faugère’s F4 algorithm for the lexicographic order followed by a division leading to a Rational Univariate Representation as the end-result. This strategy is compared with the original one which computes the Rational Univariate Representation by an order change (F4+FGLM). Table 4 gives the response times in seconds, memory space usage

in kbytes, the number of complex and real roots.

Name	F4+FGLM (s)	RR (s)	File (Kbytes)	Complex #	Real #
theo _ AFP3.rur theo_AFD5.rur theo_AFD4.rur	NA NA NA	2,95 NA 3,31	210 NA 270	36 72 36	8 16 8
eff_AFP3.rur eff_AFD5.rur eff_AFD4.rur	36,0 3600 17,9	10,9 59,6 4,8	320 290 300	40 80 40	8 16 8

Table 4.

Rational Univariate Representation response times from the RR form

Problems were encountered with the quaternion based formulation to solve the SSM. This proves that even simpler hexapod configurations can be challenging. For the real 6-6 manipulator, the RR form calculation is always faster than the computation with an order change. Inasmuch, our results indicate that the RR Rational Univariate Representation calculation using the AFD4 rotation matrix formulation is the most effective one. For the quaternion formulation, two solutions express the same manipulator posture, this is confirmed by the obtained univariate polynomial which has a total degree 80 and the number of real solutions also doubled. This behavior is well known in practice since the redundant nature of quaternions leads to two results for the same pose. Finally, the direct computation of a lexicographic Gröbner Basis is the right choice when using recent algorithms such as F4. Many instantiations were tested on various manipulator configurations and all FKP systems were in Shape Position and the RR form could be computed.

6.6 FKP of typical hexapod manipulators

6.6.1. Hexapod architecture descriptions

Figure 11.

The various hexapod alternatives

Although the method has been shown to solve the 6-6 and SSM manipulator FKP, it is relevant to compare results for various hexapods which feature specific and related geometric properties.

Several well-known parallel manipulators feature theoretical architectures, fig. 11. In order to clearly examine the configuration change impact, we have selected the first five hexapods as the slightly modified versions of each other, from the theoretical MSSM until the realistic 6-6. The 6-6p design comprises the joint positions lying on the same plane for either the fixed base or mobile platform. The 6-6pp design is characterized by the two platforms being planar. The SSM design repeats the 6-6pp along with planar hexagonal platforms which are made by equally truncating an equilateral triangular platform, [Nanua et al. 1990]. The TSSM manipulators are designed by merging the ball or Cardan joints in pairs to obtain a truly triangular mobile platform. The MSSM is merging the ball or Cardan joints of the fixed base by pairs to obtain an octahedron. These theoretical design FKP are computed and their kinematics results compared with the general 6-6 case. For each manipulator, the selected configuration and joint variables of each manipulator are selected to correspond.

6.6.2. FKP solving results

On table 5, the results are including the manipulator type, the manipulator morphology, the base and platform morphology where T, S, P and NP stands respectively for triangular, truncated triangular, planar and not planar, the complete real root computation time, the number of complex roots (the univariate polynomial degree), the number of real roots and the Gröbner Basis size. Only the Gröbner Basis sizes are compared, since the Gröbner Basis is the largest mathematical object calculated by method.

Robot type	Mechanism config	Base	Platform	Time (s)	Complex sols #	Real sols #	Grobner (kbytes)
MSSM TSSM SSM 6-6pp 6-6p 6-6	octahedron hexapod (6|) (6|) (6|) (6|)	T,P S,P S,P P NP NP	T,P T,P S,P P P NP	0,07 0,08 0,67 1,1 1,8 10,4	16 16 36 40 40 40	16 16 16 16 16 16	60 76 238 390 308 402
DIET Hexa Hexaglide	(6|) 6R-6R-6R 6T-6R-6R	NP P P	NP S,P S,P	9,9 2,0 0,5	40 40 36	40 8 8	392 346 180

Table 5.

Hexapod FKP overall results and performances

6.6.3. Discussion on the results

The gradual passage to geometrically simpler parallel structures leads to significantly shorter response times. In the 6-6 case, the introduction of one planar platform reduces computation times by 20. Inasmuch, the passage from a theoretical SSM to a realistic 6-6 leads to computation times which are 40 times longer.

With model simplification, it is notable that the real solution number is maintained. One could conjecture homothopy in those cases. The number of complex roots varies when the manipulator type changes; the MSSM and TSSM only have 16 solutions whereas the SSM has 36. This last result raises an important classification issue which was not formerly identified, [Faugère and Lazard 1995, Merlet 1997]. The SSM manipulator does not go in the class of 6-6 manipulators and becomes a class by itself.

The table second section comprises FKP solving results about the hexapod obtained with the

Dietmaier's method, [Dietmaier 1998], the Hexa and the Hexaglide with a SSM type mobile platform. In all instances, the Hexa, [Pierrot et al. 1991], and Hexaglide, [Hebsacker 1998] can feature either 36 or 40 complex solutions, this only depends on the mobile platform configuration. In fact, when the truncated equilateral triangle is used, then the FKP yield 36 solutions.

Knowing that the 6-6pp and 6-6p lead to 40 complex solutions, the new SSM class is characterized not only by its complex solution number but by the mobile platform and fixed base peculiar type.

Therefore, the manipulators which fall into the SSM category are the ones which feature one specific truncated platform.

Although Dietmaier's 6-6 FKP yields the largest number of real solutions, it does not necessarily lead to the longest computation times and largest result files.

Finally, various tests were performed where leg lengths were changed such as moving the robot on a straight line or circles with the same 6-6 manipulator configuration. The real root numbers have all been an even number in the set {4,8,12,16} depending on the location inside the workspace. This could be considered a conjecture. The only case where only one real solution has been found is when a theoretical 6-6pp has its actuator values bringing the manipulator mobile platform to lie on the fixed base plane. This solution corresponds to two real coincident roots and this case can be identified as a singularity with the loss of three DOF since the manipulator can then only move in a plane.

6.7. Assembly mode analyses

The exact method allows addressing the question of assembly modes. This problem is also referred to as posture analysis. Assembly modes are defined as follows:

Figure 12.

Posture analysis for the general 6-6 manipulator

Definition 6.1 Given a manipulator configuration where the fixed base, the mobile platform and kinematics chain lengths are specified, for a set of active joint positions, determine all the possible geometric assemblies for the selected manipulator.

For the selected example, eight real solutions were obtained which geometrically represent the only possible assembly modes. These modes are then drafted using the XMuPAD environment, fig. 12. The position based models lead to root results which are directly usable to draw the effective postures. This exemplifies that posture analysis is feasible for any manipulator which can be modelled as a general 6-6 hexapod.

7. Conclusion

In this chapter, one complete exact method to solve the parallel manipulator FKP has been explained. The method was applied to the 6-6 general parallel manipulator, which is recognized to yield the most difficult problem, and also to various other manipulators such as the SSM, TSSM, MSSM, Hexa and Hexaglide.

Moreover, the modeling of the FKP was investigated. Six displacement based models and two position-based models were derived for the 6-6 general parallel manipulator.

One complete algebraic method to solve the Forward Kinematics Problem was applied. Although many methods can find solutions to some of these FKP systems, the proposed algebraic exact method insures the exactness of the real solution results, since it is based on one Gröbner Basis, which completely describes the ideal related to the original system. From this basis, one computes the Rational Univariate Representation including one univariate equation for root isolation. The selected algorithms always succeeded to solve any parallel robot FKP in all tested non-singular instances.

The selected manipulator was a typical 6-6 hexapod known as the Gough platform in a calibrated configuration, measured on a real parallel robot prototype constructed from a theoretically singularity free SSM design.

The 8 polynomial formulations were implemented and compared. We identified three models that allowed for computation termination, out of which two were retained since their computations occur with acceptable performances: a displacement based formulation with the rotation matrix Gröbner Basis with end-effector position and rotation matrix parameters as variables and a formulation with three points on a platform.

Solving typical posture examples, the Rational Univariate Representation comprised a univariate equation of degree 40 and 8 to 12 real solutions were computed depending on the position in the workspace. The total computation averaged 130 seconds on a relatively old computer. On a faster computer, the response time falls to less than one second. Hence, this method can be suitable for small-scale trajectory pursuit applications.

This result is very important since any Gröbner Basis completely describes the mathematical object related to the original system. From this basis, one can try to build an exact equivalent system with the original one including one univariate equation. Up to author’s knowledge, this is actually the only known method that is certified to establish a truly equivalent system that preserves the properties.

Further testing led to favor the first formulation, since it yields slightly faster computations and it gives directly the end-effector position. The quaternion-based models can lead to difficulties in the case of simpler configurations and thus longer computation times, since it is doubling the solution number, which is explained by the Rational Univariate Representation equations having a degree twice as large as the others. Moreover, one final advantage is that the displacement-based equations can be applied on any manipulator mobile platform.