Analysis of computational efficiency for the solution of inverse kinematics problem of anthropomorphic robots using Gröbner bases theory

Abstract

This article presents an analysis of computational efficiency to solve the inverse kinematics problem of anthropomorphic robots. Two approaches are investigated: the first approach uses Paul’s method applied to the matrix obtained by the Denavit–Hartenberg algorithm and the second approach uses Gröbner bases theory. With each approach, the problem of inverse kinematics for an anthropomorphic robot will be solved. When comparing each method, this article will demonstrate that the method using Gröbner bases theory is more computationally efficient.

Introduction

Two classic problems in the kinematics of robotic manipulators are forward kinematics (FK) and inverse kinematics (IK). In FK, the position and orientation of the end-effector are calculated from the joint variables. In the IK, the movements of individual joints are calculated from the position and orientation of the end-effector of a robotic manipulator.

According to Daya et al., ¹ one of the biggest problems in kinematics and robot control lies in determining the solution of IK. Still, according to the same authors, as the complexity of the robotic manipulator increases, obtaining the IK is difficult and computationally expensive. Traditional methods, such as geometric, iterative, and algebraic, are inadequate as the geometric structure of the manipulator becomes more complicated.

Paul et al. ² used matrix algebra to obtain a solution to the IK problem of a specific manipulator robot. Lee ³ and Yih and Youm ⁴ used a geometric approach to solve this problem. Subsequently, Huang and Milenkovic ⁵ and Manseur and Doty ⁶ developed iterative procedures to solve the IK problem.

Denavit and Hartenberg ⁷ introduced a convention to standardize the reference coordinate systems for spatial links; after a decade, the same authors presented an algorithm for solving the kinematics problems of articulated systems. ⁸ The importance of this algorithm for the kinematics analysis of robotic systems was studied by Paul ⁹ 16 years ago.

An alternative method for the solution of IK will be used in this article, namely the Gröbner bases theory. The computational performance data obtained in this article will prove that this choice will reduce the computational cost for the solution of the IK problem of serial robotic manipulators compared to Paul’s method.

The Gröbner bases theory, developed by Bruno Buchberger in his PhD thesis, consists of an algebraic algorithm that can be applied to a given set of nonzero polynomials, producing a finite set of generators, so that it is possible to identify from these generators when a specific polynomial belongs to the original given set. ¹⁰ When finding solutions for this set of generators, those are also solutions for the whole set. With the aid of MAPLE, an algebraic computing software, all these solutions can be calculated, including cases not so easily solved using matrix algebra or geometric applications due to the complexity of the geometry of some manipulator robots.

Kendricks, ¹¹ in her PhD thesis, compared the algorithm of Denavit and Hartenberg with Gröbner bases theory for the solution of the IK problem. Her analysis focused on the robot Fanuc A-510, a SCARA type robot, with three degrees of freedom.

To familiarize the reader with each method, the anthropomorphic robot Unimation PUMA 560, a manipulator with six rotary joints, will be used as a case study applying Paul’s method from the Denavit–Hartenberg matrix and the method proposed in this work to solve the IK problem. The main goal of this article is to analyze the computational efficiency in solving the problem of the IK of anthropomorphic robots using both methods.

Guzmán-Giménez et al. ¹² presented a systematic procedure for the solution of IK problem of nonredundant open-chain robotic systems using Gröbner bases theory after solving the FK from the Denavit–Hartenberg algorithm; however, it is worth mentioning that the authors restricted themselves to the first three joints of the anthropomorphic robot Unimation PUMA 560.

Ni and Wu ¹³ presented an algorithm to solve the IK of a generic 6R robot based on matrix employing quaternions as spatial operators and the Gröbner bases theory. Through the algorithm presented by the authors, 16 different solutions are determined, but, among these, only four solutions are real. As this method is not compared to any other, its computational efficiency could not be confirmed.

Unlike Wang et al., ¹⁴ who also presented a solution for the IK of a 6R serial manipulator using Gröbner bases theory, a base with a reduced number of polynomials will be determined in this article, since these authors, using the solution method of the kinematics of mechanisms proposed by Duffy and Crane, ¹⁵ produced a Gröbner basis with 72 polynomials. With a reduced number of polynomials, which will be presented in this article, this base will increase the computational efficiency of the calculations to obtain all the joint configurations.

El-Sherbiny et al. ¹⁶ presented a comparative study between different soft-computing-based methods (artificial neural network, adaptive neuro-fuzzy inference system, and genetic algorithms) applied to the problem of IK. Their analysis focused on a serial manipulator with five degrees of freedom, and the results of the analysis demonstrate that the methods are not accurate and computationally efficient at the same time. The choice of method depends on the application, being necessary to decide which one has more priority, the calculation time or the accuracy.

However, the present work intends to go further. It brings an analysis of the computational efficiency in the solution of the IK problem of a serial manipulator with six degrees of freedom using Gröbner bases theory. The obtained results will be compared with those obtained by Paul’s method from the matrix obtained by the Denavit–Hartenberg algorithm, demonstrating that the Gröbner bases theory method is more computationally efficient and that the solutions obtained from this method are accurate.

Gröbner bases theory

The core of Gröbner bases theory is an algorithm used to solve systems of polynomial equations. It is known as Buchberger’s algorithm and applies computational algebra techniques to polynomial ideals, producing a so-called Gröbner basis, a new polynomial set.

Briefly, a Gröbner basis for a polynomial equation system is a new collection of polynomials in which the null set coincides with the solution set of the original system; however, this new system is simpler to solve. This definition justifies the choice of Gröbner bases theory to solve the IK of serial manipulators with the reduction of computational cost.

Gröbner bases theory is widely applied in other areas besides mathematics. For example, Cox et al. ¹⁷ have dedicated their work to the applications of Gröbner bases theory in geometry, graph theory, and robotics.

According to Sturmfels, ¹⁸ a Gröbner basis can be computed in order that the polynomials have a property known as the elimination property. This elimination property makes the computation of solutions of a polynomial system, which is a Gröbner basis to be mathematically efficient because, through retroactive substitutions, a solution can be found for the remaining set of variables, with a reduction in the number of floating-point operations.

For polynomials with several variables, the complexity for the determination of a Gröbner basis depends not only on the number of polynomials in the set but also on how their terms are ordered. The lexicographic order is particularly suitable for solving systems of polynomial equations. In a lexicographic order, a variable dominates any monomial involving only minor variables, regardless of the total degree. For example, for the lexicographic order $x\succ y\succ z$ , thus $x^2{\succ }_{lex}{y}^{12}{z}^8$ .

It can be seen in Figure 1 that all variables can be solved consecutively because the obtained Gröbner basis is a triangular equation system.

Figure 1.

Input and output of Buchberger’s algorithm.

In some cases, one chooses to consider the total degree of a monomial and to order the monomials of larger degrees first. In this specific case, the graduated lexicographical order (grlex) is used. Another order is given by the reverse graduated lexicographic order (grevlex), which is commonly used because it generally provides the basis more quickly due to a reduction in computational cost.

Buchberger’s algorithm constructs a Gröbner basis starting from a finite set of polynomials. According to Cox et al., ¹⁷ to develop Buchberger’s algorithm, it is necessary to understand the nature of S-polynomials.

Let f and g be two nonzero polynomials and consider $L=\text{lcm}\left(\text{lp}\left(f\right),\text{lp}\left(g\right)\right)$ , where lcm is the least common multiple among the leading power products $\left(\text{lp}\right)$ of polynomials. The combination $S\left(f,g\right)=\frac{L}{\text{lt}\left(f\right)}f-\frac{L}{\text{lt}\left(g\right)}g$ is called the S-polynomial of f and g, where (lt) stands for the leading term of each polynomial.

As an example, $f=x^3{y}^2-x^2{y}^3+x$ , $g=3x^{4}y+y^2\in ℝ\left[x,y\right]$ , with lexicographical order $x\succ y$ , it is intended to determine: $S\left(f,g\right)$ . To solve this problem, it is necessary to determine $L=\text{lcm}\left(\text{lp}\left(f\right),\text{lp}\left(g\right)\right)=\text{lcm}\left(x^3{y}^2,x^4{y}^2\right)$ . Thus

$S\left(f,g\right)=\frac{x^4{y}^2}{x^3{y}^2}f-\frac{x^4{y}^2}{3x^{4}y}g$

1

$S\left(f,g\right)=x\left(x^3{y}^2-x^2{y}^3+x\right)-\frac{y}{3}\left(3x^{4}y+y^2\right)$

2

$S\left(f,g\right)=-x^3{y}^3+x^2-\frac{y^3}{3}$

3

Note that the S-polynomial is produced to cancel the leading terms of the polynomials involved. The Buchberger’s algorithm receives as input a finite set $F=\left\{f_1,f_2,\cdots ,f_n\right\}$ of polynomials and provides a set G, which is a Gröbner basis of $F.$ According to Stifter, ¹⁹ this algorithm can be described as follows.

Algorithm 1.

Algorithm for computing Gröbner bases.

Faugère ²⁰ presented an efficient algorithm for computing Gröbner bases: the F4 algorithm. This algorithm computes successive truncated Gröbner bases, and it replaces the classical polynomial reduction found in the Buchberger’s algorithm by the simultaneous reduction of several polynomials. This powerful reduction mechanism utilizes symbolic precomputation and sparse linear algebra methods.

With the help of algebraic computing software, a Gröbner basis can be calculated more quickly. Many algebraic computing systems include packages to deal with Gröbner bases, for example, MAPLE, Singular, muMath, AXIOM, Mathematica, Macaulay, CoCoA, and Reduce.

Inverse kinematics using Paul’s method

The PUMA 560 is a robotic manipulator with six degrees of freedom and has only rotary joints, that is, it is a 6R mechanism. According to Ghosal, ²¹ PUMA 560 is a robotic manipulator extensively seen in textbooks, industry, and academic and research institutions.

The PUMA 560 is shown in Figure 2, with its six frames attached to the links.

Figure 2.

Establishing coordinate frames for the Puma 560. Source: Adapted from Chakravarti and Sivakumar. ²²

It can be seen that the reference system $\left\{0\right\}$ , not shown in Figure 2, coincides with the reference system $\left\{1\right\}$ when ${\theta }_1$ is null. It is also noted that for the Puma 560 manipulator, as for other industrial robots, the axes of joints 4, 5, and 6 intersect at a common point, and that point of intersection coincides with the origin of frames $\left\{4\right\}$ , $\left\{5\right\}$ , and $\left\{6\right\}$ , the axes of these three joints being mutually orthogonal.

According to Craig, ²³ to determine a general transformation matrix T for the solution of the FK of the Puma 560 manipulator, it is necessary to construct a transformation that defines the system $\left\{i\right\}$ in relation to the system $\left\{i-1\right\},$ that is, it is intended to determine ${}^{i-1}{A}_i$ . By defining a reference system for each link, the kinematic problem is divided into n subproblems, and to solve them, each subproblem is further broken up into four others.

Initially, three intermediate frames $\left\{P\right\}$ , $\left\{Q\right\}$ , and $\left\{R\right\}$ are defined for each link of the manipulator. It can be seen in Figure 3 that the reference system $\left\{R\right\}$ differs from the reference system $\left\{i-1\right\}$ only by a rotation of ${\alpha }_{i-1}$ , the reference system $\left\{Q\right\}$ differs from $\left\{R\right\}$ by a translation $a_{i-1}$ , the reference system $\left\{P\right\}$ differs from $\left\{Q\right\}$ by a rotation ${\theta }_i$ , and the reference system $\left\{i\right\}$ differs from $\left\{P\right\}$ by a translation d_i.

Figure 3.

Location of intermediate frames $\left\{P\right\}$ , $\left\{Q\right\}$ , and $\left\{R\right\}$ . Source: Adapted from Craig. ²³

Thus

$$\begin{gathered} {}^{i-1}A{}_i=R_{x,{\alpha }_{i-1}}\cdot T{r}_{x,a_{i-1}}\cdot R_{z,{\theta }_i}\cdot T{r}_{z,d_i} \\ =\left[\begin{array}{cccc}{c}_{{\theta }_i}& -s_{{\theta }_i}& 0& a_{i-1}\\ s_{{\theta }_i}{c}_{{\alpha }_{i-1}}& c_{{\theta }_i}{c}_{{\alpha }_{i-1}}& -s_{{\alpha }_{i-1}}& -s_{{\alpha }_{i-1}}{d}_i\\ s_{{\theta }_i}{s}_{{\alpha }_{i-1}}& c_{{\theta }_i}{s}_{{\alpha }_{i-1}}& c_{{\alpha }_{i-1}}& c_{{\alpha }_{i-1}}{d}_i\\ 0& 0& 0& 1\end{array}\right] \end{gathered}$$

4

It is used $s_{\theta }$ for $\text{sin}\left(\theta \right)$ and $c_{\theta }$ for $\text{cos}\left(\theta \right)$ .

All Denavit–Hartenberg parameters of this manipulator are listed in Table 1, where i represents the joint number.

Table 1.

Denavit–Hartenberg parameters of Puma 560.

i	${\alpha }_{i-1}(\text{rad})$	$a_{i-1}(\text{mm})$	$d_i(\text{mm})$	${\theta }_i(\text{rad})$
1	0	0	0	${\theta }_1$
2	$-\frac{\pi }{2}$	0	0	${\theta }_2$
3	0	a 2	d 3	${\theta }_3$
4	$-\frac{\pi }{2}$	a 3	d 4	${\theta }_4$
5	$\frac{\pi }{2}$	0	0	${\theta }_5$
6	$-\frac{\pi }{2}$	0	0	${\theta }_6$

After constructing all six transformation matrices using the Denavit–Hartenberg parameters and multiplying them (5), the homogeneous transformation matrix (6) is obtained

$$\begin{gathered} \begin{array}{c}T={}^{0}A{}_1\cdot {}^{1}A{}_2\cdot {}^{2}A{}_3\cdot {}^{3}A{}_4\cdot {}^{4}A{}_5\cdot {}^{5}A{}_6=\\ \left[\begin{array}{cccc}{c}_1& -s_1& 0& 0\\ s_1& c_1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]\cdot \left[\begin{array}{cccc}{c}_2& -s_2& 0& 0\\ 0& 0& 1& 0\\ -s_2& -c_2& 0& 0\\ 0& 0& 0& 1\end{array}\right]\cdot \left[\begin{array}{cccc}{c}_3& -s_3& 0& a_2\\ s_3& c_3& 0& 0\\ 0& 0& 1& d_3\\ 0& 0& 0& 1\end{array}\right] \\ \cdot \\ \left[\begin{array}{cccc}{c}_4& -s_4& 0& a_3\\ 0& 0& 1& d_4\\ -s_4& -c_4& 0& 0\\ 0& 0& 0& 1\end{array}\right]\cdot \left[\begin{array}{cccc}{c}_5& -s_5& 0& 0\\ 0& 0& -1& 0\\ s_5& c_5& 0& 0\\ 0& 0& 0& 1\end{array}\right]\cdot \left[\begin{array}{cccc}{c}_6& -s_6& 0& 0\\ 0& 0& 1& 0\\ -s_6& -c_6& 0& 0\\ 0& 0& 0& 1\end{array}\right]\end{array} \end{gathered}$$

5

$T=\left[\begin{array}{cccc}{n}_x& o_x& a_x& x\\ n_y& o_y& a_y& y\\ n_z& o_z& a_z& z\\ 0& 0& 0& 1\end{array}\right]$

6

where the entries of T are given by

$n_x=c_1\left(c_{23}\left(c_4{c}_5{c}_6-s_4{s}_5\right)-s_{23}{s}_5{c}_5\right)+s_1\left(s_4{c}_5{c}_6+c_4{s}_6\right)$

7

$n_y=s_1\left(c_{23}\left(c_4{c}_5{c}_6-s_4{s}_6\right)-s_{23}{s}_5{c}_6\right)+c_1\left(s_4{c}_5{c}_6+c_4{s}_6\right)$

8

$n_z=-s_{23}\left(c_4{c}_5{c}_6-s_4{s}_6\right)-c_{23}{s}_5{c}_6$

9

$o_x=c_1\left(c_{23}\left(-c_4{c}_5{s}_6-s_4{c}_6\right)+s_{23}{s}_5{s}_6\right)+s_1\left(c_4{c}_6-s_4{c}_5{s}_6\right)$

10

$o_y=s_1\left(c_{23}\left(-c_4{c}_5{s}_6-s_4{c}_6\right)+s_{23}{s}_5{s}_6\right)-c_1\left(c_4{c}_6-s_4{c}_5{s}_6\right)$

11

$o_z=-s_{23}\left(-c_4{c}_5{s}_6-s_4{c}_6\right)+c_{23}{s}_5{c}_6$

12

$a_x=-c_1\left(c_{23}{c}_4{s}_5+s_{23}{c}_5\right)-s_1{s}_4{s}_5$

13

$a_y=-s_1\left(c_{23}{c}_4{s}_5+s_{23}{c}_5\right)+c_1{s}_4{s}_5$

14

$a_z=s_{23}{c}_4{s}_5-c_{23}{c}_5$

15

$x=c_1\left(a_2{c}_2+a_3{c}_{23}-d_4{s}_{23}\right)-d_3{s}_1$

16

$y=s_1\left(a_2{c}_2+a_3{c}_{23}-d_4{s}_{23}\right)+d_3{c}_1$

17

$z=-a_3{s}_{23}-a_2{s}_2-d_4{c}_{23}.$

18

It is used $s_{23}$ for $\text{sin}\left({\theta }_2+{\theta }_3\right)$ and $c_{23}$ for $\text{cos}\left({\theta }_2+{\theta }_3\right)$ .

According to Craig, ²³ for the solution of IK of a robotic manipulator with six degrees of freedom, there are 12 equations and 6 unknowns. However, among the nine equations that arise from the portion of the rotational matrix of $T,$ only three are independent. These, added to the equations of the portion of the position vector of T, result in six equations with six unknowns. Such equations are nonlinear, transcendental, which can be quite difficult to solve.

The systematic approach to solve kinematic equations, proposed by Paul, ⁹ rewrites FK equations in several different ways. Thus, both sides of the FK matrix equation are left and right multiplied by inverse transformation matrices. In doing so, one gets equivalent equations having the same solutions. Each matrix equation provides 12 nonlinear equations from the three main matrix lines. The fourth line of the matrix is always trivial. The main problem is to find, in this universe, suitable equations that can be solved analytically.

The solution of IK of this robotic manipulator, as determined by Craig, ²³ is given by equations to

${\theta }_1=\mathrm{Atan}2\left(y,x\right)-\mathrm{Atan}2(d_3,\pm \sqrt{x^2+y^2-d_3^2})$

19

${\theta }_3=\mathrm{Atan}2\left(a_3,d_4\right)-\mathrm{Atan}2\left(\begin{array}{c}\left(\frac{x^2+y^2+z^2-a_2{}^2-a_3{}^2-d_3{}^2-d_4{}^2}{2a_2}\right),\\ \pm \sqrt{a_3{}^2+d_4{}^2-{\left(\frac{x^2+y^2+z^2-a_2{}^2-a_3{}^2-d_3{}^2-d_4{}^2}{2a_2}\right)}^2}\end{array}\right)$

20

${\theta }_2=\mathrm{Atan}2\left(\left(-a_3-a_2{c}_3\right)z-\left(c_{1}x+s_{1}y\right)\left(d_4-a_2{s}_3\right),\left(a_2{s}_3-d_4\right)z-\left(a_3+a_2{c}_3\right)\left(c_{1}x+s_{1}y\right)\right)-{\theta }_3$

21

${\theta }_4=\mathrm{Atan}2\left(\left(-a_x{s}_1+a_y{c}_1\right),\left(-a_x{c}_1{c}_{23}-a_y{s}_1{c}_{23}+a_z{s}_{23}\right)\right)$

22

${\theta }_4\text{'}={\theta }_4+\pi$

23

${\theta }_5=\pm \mathrm{Atan}2\left(-\left(a_x\left(c_1{c}_{23}{c}_4+s_1{s}_4\right)+a_y\left(s_1{c}_{23}{c}_4-c_1{s}_4\right)-a_z\left(s_{23}{c}_4\right)\right),\left(a_x\left(-c_1{s}_{23}\right)+a_y\left(-s_1{s}_{23}\right)+a_z\left(-c_{23}\right)\right)\right)$

24

${\theta }_6=\mathrm{Atan}2\left(\begin{array}{c}\left(-n_x\left(c_1{c}_{23}{s}_4-s_1{c}_4\right)-n_y\left(s_1{c}_{23}{s}_4+c_1{c}_4\right)+n_z\left(s_{23}{s}_4\right)\right),\\ \left(n_x\left(\left(c_1{c}_{23}{c}_4+s_1{s}_4\right)c_5-c_1{s}_{23}{s}_5\right)+n_y\left(\left(\left(s_1{c}_{23}{c}_4-c_1{s}_4\right)c_5\right)-s_1{s}_{23}{s}_5\right)-n_z\left(c_4{s}_{23}{c}_5+s_1{c}_{23}{c}_4{s}_5\right)\right)\end{array}\right)$

25

${\theta }_6\text{'}={\theta }_6+\pi$

26

All the equations of the FK and IK of the PUMA 560 manipulator were deduced in detail by Craig ²³ and Ghosal. ²¹

The tree structure in Figure 4 represents the order of computations to determine all eight joint configurations of the PUMA 560.

Figure 4.

Echelon form of the eight inverse kinematics configurations of PUMA 560.

To measure the computational cost to determine all joint configurations to position and orient the end-effector arbitrarily in the space, the following joint values: $\theta =\left[-\frac{5\pi }{6},\frac{\pi }{2},0,\pi ,\frac{\pi }{2},\frac{\pi }{6}\right]$ rad and the nonzero parameters: $a_2=432$ , $a_3=20$ , $d_3=125$ , and $d_4=430$ were replaced in equation (5), thus obtaining equation (27)

$T=\left[\begin{array}{cccc}1& 0& 0& \frac{125}{2}+215\sqrt{3}\\ 0& -1& 0& 215-\frac{125\sqrt{3}}{2}\\ 0& 0& -1& -452\\ 0& 0& 0& 1\end{array}\right]$

27

From this point, equations (7) to (18) have been rewritten by substituting n_x, n_y, n_z, o_x, o_y, o_z, a_x, a_y, a_z, x, y, and z by the respective values of the obtained numerical matrix.

The Solve command in MAPLE is used to find all the values of $\theta =\left[{\theta }_1,{\theta }_2,{\theta }_3,{\theta }_4,{\theta }_5,{\theta }_6\right]$ , which allow the robot to assume the configuration given by equation (27).

Performing 34,514 floating-point operations (2409 additions, 4261 multiplications, 175 divisions, and 3393 functions), eight different joint configurations are obtained, as given in Table 2.

Table 2.

Joint configurations of the Puma 560, determined through Paul’s method.

${\theta }_1(\text{rad})$	${\theta }_2(\text{rad})$	${\theta }_3(\text{rad})$	${\theta }_4(\text{rad})$	${\theta }_5(\text{rad})$	${\theta }_6(\text{rad})$
−2.617993878	3.091716375	−3.048636391	3.141592654	0.04307998349	0.523598777
−2.617993878	3.091716375	−3.048636391	0	−0.04307998349	−2.617993878
−0.04220293154	1.570796328	−3.048636391	0	1.477840064	−0.04220293154
−0.04220293154	1.570796328	−3.048636391	3.141592654	−1.477840064	3.099389722
−2.617993878	1.570796328	0.00000000008	0	−1.570796328	−2.617993878
−2.617993878	1.570796328	0.00000000008	3.141592654	1.570796328	0.523598777
−0.04220293154	0.04987627939	0.00000000008	0	−0.04987627939	−0.04220293154
−0.04220293154	0.04987627939	0.00000000008	3.141592654	0.04987627939	3.099389722

This article will prove that using Gröbner bases theory, the computational cost to determine the solution of the IK of serial manipulators is considerably reduced.

Inverse kinematics using Gröbner bases theory

The IK problem consists of finding all combinations of joint configurations that will position the end-effector of a robot at a given point in space. ²³ To solve this problem, one must determine polynomial equations that model the movement of the robotic arm in each joint configuration. With all these elements, one must seek a set of solutions from the possible movements of each joint.

The joint values are set to $\theta =\left[-\frac{5\pi }{6},\frac{\pi }{2},0,\pi ,\frac{\pi }{2},\frac{\pi }{6}\right]$ rad, and the parameters are set to $a_2=432$ , $a_3=20$ , $d_3=125$ , and $d_4=430$ in equation (5), exactly as in the third section, to determine all joint configurations to position and orient the end-effector of PUMA 560 robotic manipulator in space.

From this point, equations (7) to (18) are rewritten by substituting n_x, n_y, n_z, o_x, o_y, o_z, a_x, a_y, a_z, x, y, and z by the respective values of the matrix (27). These 12 polynomials, added to the six polynomials obtained from the fundamental trigonometric relation for each joint angle, are used to generate a Gröbner basis.

The selected lexicographic order for this robot is the one shown in equation (28). This order represents the same sequence in which the variables were solved while applying Paul’s method presented in the third section

$c_6\succ s_6\succ c_5\succ s_5\succ c_4\succ s_4\succ c_2\succ s_2\succ c_3\succ s_3\succ c_1\succ s_1.$

28

The same values assigned in the third section are used, intending to compare the computational efficiency to solve the IK problem using Paul’s method from the Denavit–Hartenberg matrix and Gröbner bases theory of the same 6R manipulator.

Performing 466 floating-point operations (46 additions, 100 multiplications, and 40 functions), the following Gröbner basis consisting of 12 polynomials is obtained using Faugère’s F4 algorithm ²⁰

$g_1=8600\sqrt{3}{s}_1+32084s_1{}^2+4300\sqrt{3}+2500s_1-6771$

29

$g_2=-8021\sqrt{3}{s}_1-1250\sqrt{3}+5521c_1-8600s_1-4300$

30

$g_3=1853s_3{}^2+172s_3$

31

$g_4=2c_3-43s_3-2$

32

$\begin{array}{c}{g}_5=18578641250s_1{s}_3\sqrt{3}+862257500\sqrt{3}{s}_1\\ +9289320625\sqrt{3}{s}_3+63910525900s_1{s}_3\\ +431128750\sqrt{3}+2966165800s_1\\ +2148795284s_2+9959875000s_3-665712384\end{array}$

33

$$\begin{gathered} g_6=-226593250\sqrt{3}{s}_1-113296625\sqrt{3}+537198821c_2 \\ -779480780s_1-5780183345s_3-389740390 \end{gathered}$$

34

$g_7=s_4$

35

$g_8=c_4{}^2-1$

36

$\begin{array}{c}{g}_9=433134000c_4{s}_3{s}_1\sqrt{3}-431128750s_1{c}_4\sqrt{3}\\ +216567000c_4{s}_3\sqrt{3}+1489980960c_4{s}_3{s}_1\\ -215564375\sqrt{3}{c}_4-1483082900s_1{c}_4\\ +11821264848c_4{s}_3+332856192c_4+1074397642s_5\end{array}$

37

$\begin{array}{c}{g}_{10}=-4656190500s_1{s}_3\sqrt{3}-226593250\sqrt{3}{s}_1\\ -2328095250\sqrt{3}{s}_3-16017295320s_1{s}_3\\ -113296625\sqrt{3}+537198821c_5\\ -779480780s_1-2226636864s_3-389740390\end{array}$

38

$g_{11}=-c_4{s}_1+s_6$

39

$$\begin{gathered} g_{12}=-8021s_1{c}_4\sqrt{3}-1250\sqrt{3}{c}_4-8600s_1{c}_4-4300c_4 \\ +5521c_6 \end{gathered}$$

40

It is worth noting that this reduced number of polynomials, compared to the 72 polynomials found by Wang et al., ¹⁴ will positively reflect on the computational efficiency for the IK solution of a 6R serial manipulator.

To determine all the values of $\theta =\left[{\theta }_1,{\theta }_2,{\theta }_3,{\theta }_4,{\theta }_5,{\theta }_6\right]$ , the Solve command is used in MAPLE. Eight different joint configurations are obtained, performing 954 floating-point operations (46 additions, 100 multiplications, and 101 functions), as can be seen in Table 3. These are approximately the same configurations obtained in Table 2.

Table 3.

Joint configurations of the Puma 560, obtained through Gröbner bases theory.

${\theta }_1(\text{rad})$	${\theta }_2(\text{rad})$	${\theta }_3(\text{rad})$	${\theta }_4(\text{rad})$	${\theta }_5(\text{rad})$	${\theta }_6(\text{rad})$
−2.617993878	3.091716375	−3.048636391	3.141592654	0.04307998349	0.5235987758
−2.617993878	3.091716375	−3.048636391	0	−0.04307998349	−2.617993878
−0.04220293154	1.570796327	−3.048636391	0	1.477840064	−0.04220293154
−0.04220293154	1.570796327	−3.048636391	3.141592654	−1.477840064	3.099389722
−2.617993878	1.570796327	0	0	−1.570796327	−2.617993878
−2.617993878	1.570796327	0	3.141592654	1.570796327	0.5235987758
−0.04220293154	0.04987627939	0	0	−0.04987627939	−0.04220293154
−0.04220293154	0.04987627939	0	3.141592654	0.04987627939	3.099389722

It is worth noting that all the equations presented in the third section, using Paul’s method to solve the IK of the PUMA 560 manipulator, were calculated using the computer algebra system MAPLE performing 34,514 floating-point operations to determine the eight different configurations, as given in Table 2.

Adding the 466 floating-point operations to obtain the Gröbner basis to the number of floating-point operations to obtain the eight different joint configurations through the generated Gröbner basis (954 flop), 1420 floating-point operations are obtained. There was a reduction from 34,514 to 1420 in the number of floating-point operations performed to solve the IK problem for this set of configurations.

Performance analysis

To confirm the computational efficiency for the solution of the IK of anthropomorphic robotic manipulators using Gröbner bases theory, 19 new joint configurations will be analyzed. The following joint values are used: $\left[0,0,0,0,0,0\right]\text{,}$ $\left[\frac{\pi }{6},\frac{\pi }{3},-\frac{\pi }{6},\frac{2\pi }{3},\frac{\pi }{3},-\frac{\pi }{3}\right]$ , $\left[0,-\frac{\pi }{6},0,\frac{\pi }{3},-\frac{2\pi }{3},0\right]$ , $\left[\frac{\pi }{6},-\frac{\pi }{6},\frac{\pi }{6},-\frac{\pi }{6},\frac{\pi }{6},-\frac{\pi }{6}\right]$ , $\left[\frac{\pi }{6},-\frac{\pi }{6},0,0,0,-\frac{\pi }{6}\right]$ , $\left[0,0,\frac{\pi }{2},-\frac{\pi }{2},0,\frac{\pi }{2}\right]$ , $\left[-\frac{\pi }{2},\frac{\pi }{2},0,-\pi ,-\frac{\pi }{2},\frac{\pi }{6}\right]$ , $\left[-\frac{\pi }{3},\frac{\pi }{6},\frac{\pi }{3},\frac{2\pi }{3},\frac{\pi }{2},-\frac{\pi }{2}\right]$ , $\left[0,\frac{\pi }{3},-\frac{\pi }{6},\frac{2\pi }{3},-\frac{\pi }{2},\frac{\pi }{2}\right]$ , $\left[0,\frac{\pi }{3},-\frac{\pi }{2},\frac{\pi }{2},0,\frac{\pi }{2}\right]$ , $\left[\frac{\pi }{2},0,0,-\frac{\pi }{2},0,\frac{\pi }{2}\right]$ , $\left[\frac{\pi }{6},-\frac{\pi }{6},\frac{\pi }{6},-\frac{\pi }{6},\frac{\pi }{6},-\frac{\pi }{6}\right]$ , $\left[\frac{\pi }{3},-\frac{\pi }{6},\frac{\pi }{2},-\frac{\pi }{2},\frac{2\pi }{3},-\frac{\pi }{2}\right]$ , $\left[\frac{\pi }{6},-\frac{\pi }{3},\frac{\pi }{2},\frac{\pi }{2},-\frac{\pi }{3},0\right]$ , $\left[\frac{\pi }{6},-\frac{\pi }{3},0,0,-\frac{\pi }{3},0\right]$ , $\left[\frac{\pi }{2},\frac{5\pi }{6},0,\frac{\pi }{3},-\frac{\pi }{3},\frac{\pi }{2}\right]$ , $\left[\frac{\pi }{3},\frac{5\pi }{6},-\frac{\pi }{6},\frac{\pi }{3},-\frac{\pi }{3},-\frac{\pi }{2}\right]$ , $\left[0,\frac{5\pi }{6},0,\frac{\pi }{3},-\frac{\pi }{3},-\frac{\pi }{2}\right]$ , and $\left[\frac{\pi }{6},\frac{\pi }{2},-\frac{\pi }{6},\frac{\pi }{3},-\frac{\pi }{2},\frac{\pi }{6}\right]$ .

Unlike the results presented by Guzmán-Giménez et al., ¹² all joint configurations provide the same amount of solutions determined by Paul’s method. The two methods were simulated in MAPLE, testing their ability to solve the IK problem, and all the 100 solutions obtained for both methods were compared, testing their efficiency in solving.

Comparisons of the number of mathematical operations performed using the two methods are shown in Figure 5. The x-axis presents a partition of 20 sets used to determine 100 different configurations for the PUMA 560 using Paul’s method and Gröbner bases theory, while the y-axis shows the total number of operations performed on a logarithmic scale.

Figure 5.

(a–d) Amount of operations performed using the two methods.

The performance comparisons of the two methods are shown in Figure 6. The x-axis presents a partition of 20 sets used to determine 100 different configurations for the PUMA 560 using both methods, while the y-axis shows the total number of floating-point operations on a logarithmic scale.

Figure 6.

Performance comparisons of the two methods.

It was possible to observe in Figure 5 that there was a significant reduction in the number of mathematical operations to obtain all solutions using Gröbner bases theory. As expected, the number of floating-point operations is reduced, as shown in Figure 6.

The histogram of the RMS errors in the calculations of all six joint angles of the PUMA 560 manipulator in all analyzed configurations is shown in Figure 7. The x-axis presents a partition of 100 sets of the obtained root mean square (RMS) errors in ascending order, while the y-axis shows the frequency of occurrence of each of these error sets. The RMS error is approximately $8.14\times {10}^{-10}$ rad.

Figure 7.

Histogram of the RMS errors in the computations of all six joint angles of the PUMA 560. RMS: root mean square.

It can be verified in the histogram presented that most of the obtained RMS errors occur within the first sets, the ones that represent less error. It can also be seen that the higher sets still represent very insignificant errors. Figure 8 shows the simulation results of joint variable errors in all 100 configurations analyzed.

Figure 8.

Simulation results of joint variables errors. (a) θ₁, (b) θ₂, (c) θ₃, (d) θ₄, (e) θ₅, and (f) θ₆.

The comparative analysis of IK solutions is shown in Figure 9. The red dots represent the angles obtained through Paul’s method for each joint, while the continuous blue lines represent the angles obtained through Gröbner bases theory. The data show that the outputs follow those obtained by Paul’s method with high accuracy and a negligible error.

Figure 9.

Comparative analysis of inverse kinematics solutions ( ${\theta }_1$ to ${\theta }_6$ ): (a) θ₁, (b) θ₂, (c) θ₃, (d) θ₄, (e) θ₅, and (f) θ₆.

The histograms of the RMS errors in the calculations of all positions and rotation angles of the end-effector of the PUMA 560 are shown in Figure 10. The x-axis presents a partition of 100 sets of the obtained RMS errors, while the y-axis shows the frequency of occurrence of each of these error sets. The RMS error obtained is $4.34\times {10}^{-7}$ mm for the position and $1.27\times {10}^{-9}$ rad for the orientation of the robot’s end-effector.

Figure 10.

Histograms of the RMS errors in the computations of all positions and rotation angles of the end-effector. RMS: root mean square.

All the histograms presented (Figures 7 and 10) show that most of the obtained RMS errors occur within the first sets, the ones that represent less error. It can also be seen that the higher sets still represent a completely insignificant error.

Figure 11 shows the simulation results of positioning and orientation errors in all 100 configurations analyzed.

Figure 11.

(a–f) Results of simulation of positioning and orientation errors.

The comparative analysis of FK solutions is shown in Figure 12. The green dots represent the x, y, and z coordinates of the points obtained from the angles determined through Paul’s method, while the solid orange lines represent the coordinates of the points obtained from the angles determined through Gröbner bases theory. The blue dots represent the rotation angles obtained from the angles determined through Paul’s method, while the solid red lines represent the rotation angles obtained from the angles determined through Gröbner bases theory. The data show that the results continue to maintain high accuracy and negligible error.

Figure 12.

(a–f) Comparative analysis of forward kinematics solutions.

El-Sherbiny et al. ¹⁶ presented a comparative study between three different soft-computing-based methods (artificial neural network, adaptive neuro-fuzzy inference system, and genetic algorithms) applied to the problem of IK. The analysis results demonstrate that the methods are not accurate and computationally efficient at the same time. Using the artificial neural network method with a minimized error function and the adaptive neuro-fuzzy inference system method with the same function, the calculation time is quite similar. The difference is in the accuracy of the end-effector position and the angle of rotation. When the genetic algorithm method is used, the position error is very low, which is less than ${10}^{-3}$ mm, while the calculation time is high compared to the other methods.

It is worth mentioning that in this work, the RMS error obtained is $4.34\times {10}^{-7}$ mm for the position of the end-effector and $1.27\times {10}^{-9}$ rad for the orientation of the end-effector. As shown in Figure 12, the solutions obtained through Gröbner bases theory produce negligible positioning and orientation errors in all configurations.

Open problem: Analysis of singularities

Stifter ¹⁹ presented an algorithm to solve the IK of a 2R robot using Gröbner bases theory. The author presented a brief list of subjects that were not yet fully resolved but that could play an essential role in applying Gröbner bases theory for solving the problem of IK. The question about the possibility of determining singular configurations using a Gröbner basis is highlighted in this list.

Kendricks ¹¹ used a Gröbner basis to determine the occurrence of a singularity of the workspace boundary, configuration when the manipulator is fully extended to the outer boundary or fully retracted to the inner boundary of its workspace. Her analysis focused on the robot Fanuc A-510, a SCARA type robot with three degrees of freedom, analyzing the discriminant of quadratic polynomials on the Gröbner basis generated to solve the IK of this serial robot. In this case, the author determined intervals of angles that produce a single real solution, whose discriminants in the quadratic equations are equal to zero.

The anthropomorphic robot Unimation PUMA 560, used as a case study in this work, has the last three axes intersecting at a point. As a consequence, the position of the end-effector is only a function of ${\theta }_1$ , ${\theta }_2$ , and ${\theta }_3$ . In this way, the following three nonlinear equations in three unknowns are used

$x=c_1\left(a_2{c}_2+a_3{c}_{23}-d_4{s}_{23}\right)-d_3{s}_1$

41

$y=s_1\left(a_2{c}_2+a_3{c}_{23}-d_4{s}_{23}\right)+d_3{c}_1$

42

$z=-a_3{s}_{23}-a_2{s}_2-d_4{c}_{23}$

43

To check whether it is possible to determine from a Gröbner base whether there is a singularity, equations (41) to (43) are rewritten by substituting x, y, and z by the respective values of the matrix (27). These three polynomials, added to the three polynomials obtained from the fundamental trigonometric relation for each joint variable, are used to generate a Gröbner basis.

The selected lexicographic order for this robot is the one shown in equation (44). This order represents the same sequence in which the first three joint variables were solved while applying Paul’s method presented in the third section

$c_2\succ s_2\succ c_3\succ s_3\succ c_1\succ s_1$

44

The following Gröbner basis consisting of six polynomials is obtained using Faugère’s F4 algorithm ²⁰

$g_1=(x^2+y^2)s_1{}^2+250x{s}_1-y^2+15625$

45

$g_2=c_{1}y-s_{1}x-125$

46

$\begin{array}{c}{g}_3=138325708800s_3{}^2+(743040x^2+743040y^2\\ +743040z^2-287964408960)s_3+x^4+2x^2{y}^2\\ +2z^2{x}^2+y^4+2y^2{z}^2+z^4-775098x^2\\ -775098y^2-775098z^2+149895629001\end{array}$

47

$g_4=-x^2-y^2-z^2+17280c_3-371520s_3+387549$

48

$\begin{array}{c}{g}_5=(16009920x^2+16009920y^2)s_1{s}_3\\ +(43x^4+86x^2{y}^2+43x^2{z}^2+43y^4+43y^2{z}^2\\ -16664607x^2-16664607y^2)s_1\\ +(1728x^{2}y+1728y^3+1728y{z}^2-27000000y)s_2\\ +2001240000x{s}_3+2x^{2}yz+2y^{3}z+2y{z}^3+5375x^3\\ +5375x{y}^2+5375z^{2}x-28602yz-2083075875x\end{array}$

49

$\begin{array}{c}{g}_6=(1728x^{2}y+1728y^3+1728y{z}^2-27000000y)c_2\\ +(-2x^4-4x^2{y}^2-2x^2{z}^2-2y^4-2y^2{z}^2+28602x^2\\ +28602y^2)s_1+16009920s_{3}yz+43x^{2}yz+43y^{3}z\\ +43y{z}^3-250x^3-250x{y}^2-250z^{2}x-16664607yz\\ +3575250x\end{array}$

50

Note that the first element in the basis is a quadratic polynomial in terms of s₁. By examining the discriminant of this equation, it can be observed that there will be no real solutions for s₁, when

$4y^2\left(x^2+y^2-15625\right)<0$

51

Moreover, by analyzing the discriminant of the third element in the basis, there will be no real solutions to s₃ when

$$\begin{gathered} -x^4-2x^2{y}^2-2x^2{z}^2+775098x^2-y^4-2y^2{z}^2 \\ +775098y^2-z^4+775098z^2-11868518601<0 \end{gathered}$$

52

To obtain a quadratic polynomial in terms of s₂, a new lexicographic order for this robot is used: $c_3\succ s_3\succ c_2\succ s_2\succ c_1\succ s_1$ . The following Gröbner basis consisting of six polynomials is obtained using Faugère’s F4 algorithm ²⁰

$g_1=(x^2+y^2)s_1{}^2+250x{s}_1-y^2+15625$

53

$g_2=c_{1}y-s_{1}x-125$

54

$\begin{array}{c}{g}_3=(746496x^2+746496y^2+746496z^2-11664000000)s_2{}^2\\ +(1728x^{2}z+1728y^{2}z+1728z^3-24712128z)s_2\\ +x^4+2x^2{y}^2+2z^2{x}^2+y^4+2y^2{z}^2+z^4-775098x^2\\ -775098y^2-28602z^2+11868518601\end{array}$

55

$\begin{array}{c}{g}_4=(864x^{2}y+864y^3-13500000y)c_2\\ +(-864x^{2}z-864y^{2}z)s_1{s}_2+(-x^4-2x^2{y}^2\\ -x^2{z}^2-y^4-y^2{z}^2+14301x^2+14301y^2)s_1\\ -108000xz{s}_2-125x^3-125x{y}^2-125z^{2}x+1787625x\end{array}$

56

$\begin{array}{c}{g}_5=(1728x^4+3456x^2{y}^2+1728x^2{z}^2+1728y^4\\ +1728y^2{z}^2-27000000x^2-27000000y^2)s_1{s}_2\\ +(2x^{4}z+4x^2{y}^{2}z+2x^2{z}^3+2y^{4}z+2y^2{z}^3-28602x^{2}z\\ -28602y^{2}z)s_1+(216000x^3+216000x{y}^2+216000x{z}^2\\ -3375000000x)s_2+(16009920x^{2}y+16009920y^3\\ -250155000000y)s_3+43x^{4}y+86x^2{y}^3+43x^{2}y{z}^2\\ +43y^5+43y^3{z}^2+250x^{3}z+250x{y}^{2}z+250x{z}^3\\ -17336482x^{2}y-17336482y^3-671875y{z}^2\\ -3575250xz+260384484375y\end{array}$

57

$\begin{array}{c}{g}_6=(16009920x^{2}y+16009920y^3-250155000000y)c_3\\ +(37152x^4+74304x^2{y}^2+37152x^2{z}^2+37152y^4\\ +37152y^2{z}^2-580500000x^2-580500000y^2)s_1{s}_2\\ +(43x^{4}z+86x^2{y}^{2}z+43x^2{z}^3+43y^{4}z+43y^2{z}^3\\ -614943x^{2}z-614943y^{2}z)s_1\\ +(4644000x^3+4644000x{y}^2+4644000x{z}^2\\ -72562500000x)s_2-2x^{4}y-4x^2{y}^3-2x^{2}y{z}^2-2y^5\\ -2y^3{z}^2+5375x^{3}z+5375x{y}^{2}z+5375x{z}^3+806348x^{2}y\\ +806348y^3+31250y{z}^2-76867875xz-12110906250y\end{array}$

58

Note that the third element in the basis is a quadratic polynomial in terms of s₂. By examining the discriminant of this equation, it can be observed that there will be no real solutions for s₂, when

$\begin{array}{c}-(2985984(x^2+y^2-15625))(x^4+(2y^2+2z^2-775098)x^2\\ +y^4+(2z^2-775098)y^2+z^4-775098z^2+11868518601)<0\end{array}$

59

When it is possible to determine an analytical solution for the IK of a manipulator robot using Gröbner bases theory, given a point $\left(x;y;z\right)$ , it can determine the maximum number of geometric configurations to position the end-effector of the robot, and, more importantly, it is possible to discover all the real solutions that produce these settings.

However, it is essential to highlight that when all six joint variables are considered in the model, the 6R manipulator presents more than four solutions; however, the Gröbner bases method does not produce analytical solutions, thus requiring numerical methods to find the final solution. In many industry fields, the determination of all rotation angles of the end-effector is indispensable; thus, it justifies the importance of determining all joint angles of a robotic manipulator.

Furthermore, although it is possible to determine angular intervals that produce points that are external to the manipulator’s workspace, singularities such as boundary singularity or inner singularity remain to be further studied.

Conclusions

When using the Gröbner bases method for the solution of IK of anthropomorphic robots, one notes in the tests performed during this research that it is computationally more efficient than Paul’s method since the equations produced for determination of six variables are mathematically simpler to be solved. Since through retroactive substitutions, all joint solutions can be determined to precisely position and orient the robotic manipulator end-effector in space.

According to the performance data obtained (number of floating-point operations) in all the configurations analyzed, it can be stated that the method proposed in this work for the solution of IK of anthropomorphic robots is more efficient computationally than Paul’s method.

In future work, it is intended to demonstrate that this method can also be computationally efficient compared to the other methods for solving the IK problem of redundant manipulators.