Inverse Displacement Analysis of a General 6R Manipulator Based on the Hyper-Chaotic Least Square Method

1. Introduction

The inverse displacement analysis of a robot plays an important role in robotics, which is directly related to the problems of the off-line programming, trajectory planning and real-time control. It is most difficult in relation to inverse displacement analysis on a general space 6R manipulator. This issue with the 6R mechanism belongs to the same problem, which has been called the Mount Everest of the position analysis on a space mechanism [1]. In 1973, J. Rastegar and V. Scheinman [2] pointed out that the number of the solutions was up to 32 in the inverse displacement analysis of a general 6R serial robot. In 1985, L. W. Tsai and A. P. Morgan [3] simplified the above problem in eight quadratic equations using the high dimensional approximation analysis on the inverse kinematics problem and carried out the numerical solution adopting a polynomial continuous method, and then speculated as to the conclusion that the solution number is up to 16 on a variety of different structures of 6R robot. In 1989, M. Raghavan and B. Roth [4] made use of the nature of the rationale generated by using the separation elimination method and multi-variable equations, and calculated a 16th degree polynomial from the half-angle tangent formula of the joint variables. In 1991, C. Wampler and A. P. Morgan [5] proposed that all the cases involving the inverse displacement problem were only solved using the extension method, but the algorithm based on the extension method ran very slowly. From 1990 to 1992, M. Raghavan and B. Roth [6-8] first proposed a method solving the characteristic polynomial of the general 6-DOF manipulator, including all the special circumstances of 6R, 5R1P, 4R2P, 3R3P, and with the aim of solving all of the 6-DOF robot mechanism. In 1991, C. W. Wampler and A. P. Morgan [9] completed the 6R inverse position problem using a generic-case solution methodology with 11 polynomial equations. In 1993, D. Kohli and M. Osvatic [10] used the product of the powers and the elimination method to obtain a set of equations containing 16 × 16 matrix coefficients, and to directly have 16th degree polynomials of no extraneous roots, since the solution to the problem can be simplified to the eigenvector problem. In the same year, D. Manocha and J. F. Canny [11] made use of the symbol processing and matrix computation to transform this problem into the decomposition problem of the matrix to be calculated, and solved the matrix eigenvalues and eigenvectors to obtain the robot joint variables. However, the algorithm finally computed the 24th-order square matrix to solve 24 solutions, including extraneous roots rather than directly the 16 solutions. In 1999, H. J. Su et al. [12] proposed an algorithm based on algebraic elimination to solve the robot inverse kinematics. The method can get a one-variable equation with one joint variable first of all, which greatly reduces the computation procedure of the real-time iterative. However, the algorithm needs to expand the determinant polynomial to solve 16th degree polynomials. A slower computing speed and may occur and cause problems when calculating the polynomial. Y. Q. Yu [13] solved the inverse kinematics problem of a general 6R serial robot and simulated the robot motion by using VC++ with OpenGL in 2003, but there was a poor consequence upon the real-time reverse problem. How to quickly find all the real solutions is important to the research of the mechanism, but it is also one of the basic problems for the general 6R and 7R robot mechanisms.

The research into chaos is one of the important achievements of the 21st century. How to utilize chaotic characteristics is one of the important works of the modern mechanism of research. I have solved the forward displacement problem of 6-RPS in the real area using the chaotic method in [14]. With this method, it advances that the points of a Julia set, which are generated by the Newton iteration method, will appear in the neighbourhood where the Jacobian matrix of the solving equations is equal to zero. However, this supposition has not been proven, and the process for solving the multivariable Jacobian matrix is very complex. The chaotic sequence method is new, in which the initial point of the Newton iteration is generated using the chaotic and hyper-chaotic system and all the real solutions on the mechanism synthesis can be effectively solved [15-17]. When the solutions are not in convergence using Newton method or a quasi-Newton method, the mathematical programming method can be adopted [18]. However, its computational efficiency is very low for the process of solving the inverse kinematics problem of a 6R mechanism by using the mathematics quadratic programming method based on hyper-chaos. With the Newton iteration method, the requirement for the initial value is high. However, the least square method can expand the selected scope of the initial value, and it can sometimes make one initial value convergence while Newton method cannot. For these reasons, combining the least square method with hyper-chaotic sequences and utilizing the initial value that was iterated in the least square method based on hyper-chaotic Hénon mapping, the inverse kinematics problem of a general 6R mechanism was solved in this paper. The example provided shows that the method proposed in this paper is verified as correct and effective.

2. Hénon hyper-chaotic system

The Lyapunov exponent is one of the effective methods depicting the chaos-specific property of nonlinear systems. If one of the Lyapunov exponents is positive, the system is chaotic, and if a system has two or more positive Lyapunov exponents, the system is hyper-chaotic. The greater the number of the positive Lyapunov exponents, the higher the degree of instability in the system [19]. In general, if the systematic state variable number is higher, it probably appears that the unsteady level is also higher.

A general Hénon mapping was designed [20] as follows:

{ x 1 , k + 1 x i , k + 1 } = { a − x n − 1 , k 2 − b x n , k x i − 1 , k }

(1)

where, i = 2 , 3 , ⋯ , n 1 expresses the dimension of the system, k is the discrete time, and a and b are adjustable parameters. When i=2, the above mapping is the famous Hénon mapping. When fixed parameters a=1.76 and b=0.1, and the dimensions vary from 2 to 10, after computing, it was found that by increasing n 1 , the simple relation of the number n of the positive Lyapunov exponents with the system dimension n 1 is n 1 = n + 1 in [20] - namely, when a system dimension is the larger of two, the system is hyper-chaotic. For n>10, a simulation study was conducted [18] which also obtained the same result.

3. The least squares method in nonlinear equations

The nonlinear equation is expressed as:

f ( x ) = [ f 1 ( x ) , ⋯ , f n ( x ) ] = 0

(2)

Its solution is written in x = [ x 1 , x 2 , ⋯ , x n ] T . Assuming J k = ∂ f ( x k ) ∂ x k , the least squares iterative method can be described as follows, roughly:

x k + 1 = x k − ( J k T J k ) − 1 J k T f ( x k ) ( k = 0 , 1 , 2 , ⋯ )

(3)

where f ( x k ) is the function value of f ( x ) at the point x k , and J k is the Jacobian matrix of f ( x ) at the point x k . Although J k T J k is always symmetric, positive and semi-definite, it seems that usually there exists the inverse matrix ( J k T J k ) − 1 . But the value of its determinant det ( J k T J k ) is so small that there exists a serious ill-conditioned problem. Many scholars have put forward improved algorithms, and the most well-known of these is the damped least square method, also known as the L-M method. Recently, a Chinese scholar, C. X. Zhan, proposed a new method called the C. X. Zhan method in a paper. With the new method, computational efficiency is better than with the L-M method, and the convergence rate is faster. So far, this method is the best for solving the least squares problems, and it is adopted in this paper. Its basic steps are as follows: Firstly, it is assumed that A k = J k T J k , A k is decomposed to A k = L k D k L k T (where, L k is a lower triangular matrix, and D k is a diagonal matrix). Secondly, Eq. (3) is rewritten as L k D k L κ T ( x k + 1 − x k ) = − J k T f ( x k ) .

Thirdly, after increasing the damping to D k , the iterative equation based on the C. X. Zhan method is expressed as:

L k ( D k + μ I ) L κ T ( x k + 1 − x k ) = − J k T f ( x k )

(4)

where I is the n-order unit matrix, and μ k is the damping factor where μ k > 0 . The selection of the damping factor and the convergence criteria isseen in [21].

4. The hyper-chaotic least square method in nonlinear equations

Using the least square method based on a hyper-chaotic system, the process for finding all of the solutions of nonlinear equations is as follows:

5. The mathematical model

According to the structural parameters of the robot - i.e. the D-H parameters - and the link coordinate system, the general 6R robot kinematics can be described as:

A = A 1 A 2 A 3 A 4 A 5 A 6

(5)

where A i is a 4 × 4 homogeneous transformation matrix from the i t h link coordinate system to the ( i − 1 ) t h coordinate system; A is a homogeneous transformation matrix from the coordinate system of the robot end to one of the base.

A i = [ Q i P i 0 T 1 ] , Q i = Z i X i , X i = [ 1 0 0 0 − cos α i sin α i 0 sin α i cos α i ] , Z i = [ cos θ i sin θ i 0 sin θ i − cos θ i 0 0 0 1 ] , P i = [ a i cos θ i a i sin θ i b i ] T

where, a i , b i and α i are the structural parameters of the robot, and they are determined by the specific robot structure; θ i ( i = 1 , 2 , ⋯ , 6 ) is the Joint variable. The inverse kinematics problem of the general 6R robot is that the joint variable θ i is found under the conditions of giving a i , b i , α i and A . Eq. (5) is expanded as

[ A ( 1 , 1 ) A ( 1 , 2 ) A ( 1 , 3 ) A ( 1 , 4 ) A ( 2 , 1 ) A ( 2 , 2 ) A ( 2 , 3 ) A ( 2 , 4 ) A ( 3 , 1 ) A ( 3 , 2 ) A ( 3 , 3 ) A ( 3 , 4 ) 0 0 0 1 ] = [ A 16 ( 1 , 1 ) A 16 ( 1 , 2 ) A 16 ( 1 , 3 ) A 16 ( 1 , 4 ) A 16 ( 2 , 1 ) A 16 ( 2 , 2 ) A 16 ( 2 , 3 ) A 16 ( 2 , 4 ) A 16 ( 3 , 1 ) A 16 ( 3 , 2 ) A 16 ( 3 , 3 ) A 16 ( 3 , 4 ) 0 0 0 1 ]

(6)

In Eq. (6), the corresponding parameters are equal and there are twelve scalar equations and four identities. To the mechanism, six of these scalar equations only are independent. Now, with Eq. (7) substituted into twelve scalar equations we can obtain Eq. (8):

x 1 = sin θ 1 , x 2 = cos θ 1 , x 3 = sin θ 2 , x 4 = cos θ 2 x 5 = sin θ 3 , x 6 = cos θ 3 , x 7 = sin θ 4 , x 8 = cos θ 4 x 9 = sin θ 5 , x 10 = cos θ 5 , x 11 = sin θ 6 , x 12 = cos θ 6

(7)

F = [ f 1 ( x ) , ⋯ , f 12 ( x ) ] T = 0

(8)

where:

x = [ x 1 , x 2 , ⋯ , x 12 ] T , f 1 = A ( 1 , 1 ) − A 16 ( 1 , 1 ) , f 2 = A ( 1 , 2 ) − A 16 ( 1 , 2 ) , f 3 = A ( 1 , 3 ) − A 16 ( 1 , 3 ) , f 4 = A ( 1 , 4 ) − A 16 ( 1 , 4 ) f 5 = A ( 2 , 1 ) − A 16 ( 2 , 1 ) , f 6 = A ( 2 , 1 ) − A 16 ( 2 , 2 ) f 7 = A ( 2 , 3 ) − A 16 ( 2 , 2 ) , f 8 = A ( 2 , 4 ) − A 16 ( 2 , 4 ) f 9 = A ( 3 , 1 ) − A 16 ( 3 , 1 ) , f 10 = A ( 3 , 2 ) − A 16 ( 3 , 2 ) f 11 = A ( 3 , 3 ) − A 16 ( 3 , 3 ) , f 12 = A ( 3 , 4 ) − A 16 ( 3 , 4 ) .

Without considering sin 2 θ i + cos 2 θ i = 1 , twelve scalar equations of Eq. (8) are independent. After the solution, the results should be tested as to whether they meet the condition equation sin 2 θ i + cos 2 θ i = 1 . Removing those extraneous roots which are not satisfied conditions, the solutions in Eq.(8) are found.

Eq.(5) is transformed into Eq.(9):

A 3 A 4 A 5 = A 2 − 1 A 1 − 1 A A 6 − 1

(9)

After expansion, the equation is as follows:

[ l 11 ( θ 3 , θ 4 , θ 5 ) l 12 ( θ 3 , θ 4 , θ 5 ) l 13 ( θ 3 , θ 4 , θ 5 ) l 14 ( θ 3 , θ 4 , θ 5 ) l 21 ( θ 3 , θ 4 , θ 5 ) l 22 ( θ 3 , θ 4 , θ 5 ) l 23 ( θ 3 , θ 4 , θ 5 ) l 24 ( θ 3 , θ 4 , θ 5 ) l 31 ( θ 4 , θ 5 ) l 32 ( θ 4 , θ 5 ) l 33 ( θ 4 , θ 5 ) l 34 ( θ 4 , θ 5 ) 0 0 0 1 ] = [ r 11 ( θ 1 , θ 2 , θ 6 ) r 12 ( θ 1 , θ 2 , θ 6 ) r 13 ( θ 1 , θ 2 ) r 14 ( θ 1 , θ 2 ) r 21 ( θ 1 , θ 2 , θ 6 ) r 22 ( θ 1 , θ 2 , θ 6 ) r 23 ( θ 1 , θ 2 ) r 24 ( θ 1 , θ 2 ) r 31 ( θ 1 , θ 2 , θ 6 ) r 32 ( θ 1 , θ 2 , θ 6 ) r 33 ( θ 1 , θ 2 ) r 34 ( θ 1 , θ 2 ) 0 0 0 1 ]

(10)

where l i j and r i j are the non-trivial components of the left and right in Eq.(9), respectively. Since the values of the third and fourth columns in Eq.(10) are not associated with θ 6 , utilizing the non-trivial components of the third and fourth columns obtains six equations:

F 1 = [ f 1 ( x ) , ⋯ , f 6 ( x ) ] T = 0

(11)

where:

f 1 = l 13 − r 13 , f 2 = l 14 − r 14 , f 3 = l 23 − r 23 f 4 = l 24 − r 24 , f 5 = l 33 − r 33 , f 6 = l 34 − r 34

Making use of the equation sin 2 θ i + cos 2 θ i = 1 six equations are added as:

F 2 = [ f 7 ( x ) , ⋯ , f 12 ( x ) ] T = 0

(12)

where

: f 7 = x 1 2 + x 2 2 − 1 , f 8 = x 3 2 + x 4 2 − 1 , f 9 = x 5 2 + x 6 2 − 1 , f 10 = x 7 2 + x 8 2 − 1 , f 11 = x 9 2 + x 10 2 − 1 , f 12 = x 11 2 + x 12 2 − 1

The calculation based on Eq. (11) and Eq. (12) yields the solution x = [ x 1 , x 2 , ⋯ , x 12 ] T . Because there is no θ 6 in Eq. (11), [ x 1 , x 2 , ⋯ , x 10 ] T as the known variables are substituted into Eq.(13), and then x 11 and x 12 are found.

l 11 − r 11 = 0 , x 11 2 + x 12 2 − 1 = 0

(13)

where the equation l 11 − r 11 = 0 can be replaced by the other scalar equation corresponding to the first or second column. After the solution, we can check whether the first column in Eq.(10) is equal to the second one.

6. Numerical example

Given the mechanism parameters shown in Tab.1 and the pose of the end, as follows:

Q = [ − 0.357279 − 0.85000 0.387106 0.915644 − 0.237000 0.324694 − 0.184246 0.470458 0.862973 ] P = [ 0.798811 − 0.000331 1.200658 ] T ,

where the unit of P is m, the inverse displacement of the mechanism should be found.

OPEN IN VIEWER

Table 1.

The parameters of a 6R manipulator

i	a i (m)	b i (m)	αi(°)	θ i
1	0.12	0	−57	θ 1
2	1.76	0.89	35	θ 2
3	0.07	0.25	95	θ 3
4	0.88	−0.43	79	θ 4
5	0.39	0.50	−75	θ 5
6	0.93	−1.34	−90	θ 6

OPEN IN VIEWER

Table 2.

The computing results of the variables

No.	θ 1	θ 2	θ 6
1	177.5385	−148.1792	110.9850
2	174.0832	−163.3024	141.2812
3	164.8003	−154.2911	−101.3592
4	1.2272	−7.3527	149.2076
5	−16.4788	−10.7473	5.6779
6	−22.2616	−22.4296	−17.1504
7	−22.6026	28.0947	169.8786
8	−41.6831	−29.1292	25.0910
9	−46.0139	−19.2569	−114.7688
10	−53.1779	26.1664	127.9769
11	−83.0949	57.0232	−110.9821
12	−137.1953	−156.9203	147.4464
13	−139.0594	128.1143	−119.8347
14	−148.7751	−179.7129	55.7114
15	−159.8441	−159.3355	21.6746
16	−173.9288	150.6972	−92.2836

No.	θ 3	θ 4	θ 5
1	159.4308	−148.6461	−129.2791
2	−164.7918	−107.8183	−155.7380
3	−85.3408	4.7800	−127.8081
4	142.6975	−123.8783	−29.2146
5	−58.8964	−4.1618	164.0779
6	−32.0255	−32.4075	−172.6183
7	98.6338	−176.2460	12.4561
8	52.3566	6.5588	−129.1205
9	−46.9884	−120.2180	−145.8650
10	9.1030	145.8678	136.3517
11	130.9760	67.5704	−10.8278
12	68.3066	135.6859	−51.3478
13	96.0421	25.4394	−7.3520
14	−78.5052	158.0864	148.2539
15	−111.3474	120.2692	176.5988
16	47.8112	−21.0006	−40.4385

After a random number is generated, an initial point of the hyper-chaotic set and the number of hyper-chaotic variables n 1 is selected as 13, according to the generalized Hénon hyper-chaotic system in Eq. (1), and the hyper-chaotic sequence including 12 positive Lyapunov exponents was built (the selected chaotic sequence length N was 100. If N is selected as bigger, the computing time would be longer, and if N is too small, it may be that not every real solution would be found. Therefore, in relation to these problems of the mechanism synthesis, N should first be taken as smaller when counted approximately, and then be selected as larger. In general, N is selected at (8–12) times the variable number). After an iterative calculation for 948s, Eq.(8) was solved to find x . Then, θ i ( i = 1 , 2 , ⋯ , 6 ) was converted, as shown in Table 2. These results were the same as, e.g., 8.2.3 of [22]. If solving Eq.(11), Eq.(12) and Eq.(13) from Eq. (9), 16 groups of solutions were found after running for 16.5s. These results were the same as, e.g., 8.2.3 of [22] and the computational efficiency increased greatly. Since the hyper-chaotic Newton iteration method in [15] and the hyper-chaotic sequence from Eq.(1) made the matrix singular, the result was not found. In the mathematical programming method based on hyper-chaos in [18], the fmincon function only obtained two groups of real solutions after running many times, while the fminimax function was only obtained for one or two groups. For the first kind of constraint equations, sin θ i and cos θ i ( i = 1 , ⋯ , 6 ) were not restricted, i.e., sin 2 θ i + cos 2 θ i = 1 was not used in the solution process, and the variables x 1 , x 2 , ⋯ , x 12 were dealt with directly. Thus extraneous roots were increased, i.e., there existed x i 1 > 1 ( i = 1 , 2 , ⋯ , 12 ) . The second type of constraint equation has no extraneous roots by the use of sin 2 θ i + cos 2 θ i = 1 , and so the efficiency of its solution is relatively high.

7. Conclusions

Applying the D-H method, a 4 × 4 transform matrix was obtained and the first type twelve constrained equations were established. Then, using the characteristics of the matrix after transforming, six constraint equations were established. Six supplement equations were also established by increasing six variables and the relation of the sine and cosine functions, and the second type of twelve constrained equations for the position analysis were established. Combining the least square method with hyper-chaotic sequences, a hyper-chaotic least square method based on utilizing a hyper-chaotic discrete system to obtain and locate initial points to find all the real solutions of the nonlinear questions was proposed, and the calculation steps solving the inverse displacement analysis of 6R manipulator containing two types of constrained equations were given. Two types of constraint equations can be solved, but the second type of constraint equation has a high efficiency. This method solved the problem that it is divergent when utilizing a Newton iterative method based on chaos and hyper-chaos, the quasi-Newton method, and the hyper-chaotic mathematical programming method. The numerical example shows that this new method is correct and effective. Using this method in the real area, it is more efficient and the result is the same as the homology analysis method. The new method provided in this paper is also suitable for solving the inverse displacement analysis of 7R manipulator.