Robotics: Heuristics Based Singularity Free Inverse Kinematics and Jacobian Inverse

DOI : 10.17577/IJERTV10IS120136

Text Only Version

Robotics: Heuristics Based Singularity Free Inverse Kinematics and Jacobian Inverse

Chaman Lal Sabharwal

Missouri University of Science and Technology Rolla, Missouri-65409, USA

Abstract –In Robotics, inversion is encountered at two levels: Inverse Kinematics level and Inverse Jacobian level. We present a new and intuitive approach to simplify the existing analytical complex approaches for inversion: Kinematics and Jacobian. Inverse Kinematics also employs geometric approach. We present a different style to efficient geometric approach. We show that specification of end-effector position is not necessary. Robotic Jacobian is an mÃ—n matrix with n degrees of freedom (DOF). The computation of inverse velocities is an issue when m n, for which we provide an intuitive and common sense approach to the Generalized Inverse. Thus, this paper provides intuitive exciting approaches for efficient computing of inverse kinematics and inverse Jacobian in robotics.

Keywords: Link parameters; inverse kinematics; inverse jacobian; generalized inverse; degrees of freedom (DOF)

1. INTRODUCTION

A robot is a rigid body compose of links and joints on links for motion, see Fig. 1; Fig. 2. Kinematics is the study of motion without regard to the forces that create it. The representation of the position of the robot end effecter through the Robotics Engineering (common parameters and link) forward kinematics [1], [2], [3]. Forward kinematics yields a unique orientation and position of end-effector, but the end-effector may be reached in multiple ways. There can be several paths for a robot to reach a destination. If the link parameters for joints are specified, forward kinematics uses these joint parameter values to reach a unique destination. Inverse kinematics of robotic manipulator involves obtaining the required values manipulator joint position given the desired end point and direction. There is no unique solution and close the form of direct expression of inverse kinematics mapping [2], [3], [4]. However, if the orientation and position of end-effector is specified, then Inverse Kinematics determines possible link parameter values that may result in multiple solutions, some of which may be extraneous and inconsistent. Inverse kinematics employs a bag of heuristics using one trigonometric equation and avoiding division by zero singularities. This analysis explores the multiple solutions and exploits the valid solutions. We will devise common sense approach.

For inverse kinematics in Geometric approach, to reach intermediate frames, it is computationally efficient to start from base origin [see Fig. 5, Fig. 6] instead of instead starting from end-effector. It is also intuitive to grasp. Here we provide a heuristic approach to solve the required equations. Also, Inverse kinematics employs geometric approach to break down the 6 DOF into 3 + 3 [see Fig. 5] or

4 + 2 [see Fig. 6] etc. Motion occurs about or along the axes at the joints where frames are specified or computed.

Multiple solution issue with Jacobians is resolved by using minimum norm criteria. In practice, a robot has six DOF, articulate robot may have all revolute links. For a robot with 6 DOF, normally arm has three links and wrist has three links. Some robots may have two to any number of links. In any case, in general, the number of links may not be equal to six [5]. As in forward kinematics, application of forward velocity is also straight forward. In fact, this becomes an issue only when the Jacobian is not a square matrix. It happens when the system of equations is underdetermined or overdetermined. For this Inverse Jacobian problem, there are algebraic approaches to resolve this [2],[3]. We provide more efficient ways to grasp the concepts of inverse kinematics and inverse Jacobian.

The paper is organized as follows: Section II is overview of background, Section III describes innovative techniques for computing Inverse Kinematics, Section IV is on Geometric Inverse Kinematics, efficiency, section V is intuitive Inverse Jacobian, Section VI is Conclusion followed by References.

2. OVERVIEW OF BACKGROUND

A robot is a rigid body composed of links moving around the joint axes. The linkk parameters are (k, dk, ak, k) relative to linkk-1, are used to determine a coordinate system, the link frame. The link framek is computed by starting at the origin of framek-1 in four transformations (1): rotate about zk-1 -axis by angle k , translate along zk-1-axis by amount dk , translate along xk-axis by amount k, and rotate about xk-axis by angle k . The composition of these four transformations called an A-matrix or coordinate frame. The A matrix is akin to transformation from one frame to next frame. The frame for transformation from frame k-1 to frame k becomes A- attribute matrix, synonymous with transformation matrix, T or the frame matrix, F. For notation: k-1Ak represents transformation from frame k-1 to frame k. Equivalently, it is frame k relative to frame k-1.

Figure 1. Description of joint coordinates and link parameters (k, dk, ak, k) [2].

The link parameters (k, dk, ak, k). Last two (ak, k) are constant all the time, the first two (k, dk) are variable one of which is constant at the times while other is variable for motion. The equation for composition is as follows

k-1Ak k-1Fk k-1Tk =

R(zk-1, k) T(zk-1, [0, 0, dk]), T(xk, [ak, 0, 0]) R(xk, k) (1)

with c2 (a2 + b2) . Its solution is devoid of singularities of division by zero.

For velocity,the velocity equation uses the Jacobian, J,

[ ] = J [ ] (5)

Forward Jacobian: Velocities [ ] at the joints are specified, the problem is to find end-effector velocity [ ] .

Inverse Jacobian: The end-effector velocity [ ] is given,

the problem is to find link velocities [ ] at the joints.

3. TECHNIQUES FOR COMPUTING INVERSE KINEMATICS,

The inverse kinematics involves solving trigonometric equations. If such equation is not solved properly, it may lead to divide by zero singularities in the solution [3].

For example, in the case of a six link robot with first three links create the arm, we will first solve the arm before we go into full scale 6 links. For the sake of completeness, we include the robot Fig. 2, link parameter table Table 1, and the A matrices, etc. For example, Jumbo drilling robot is displayed here in Fig. 2. The Fig. 2 has clearly marked joint axes and link parameters. One may flatten the robot, see Fig.

= [

0

0

0

0

] (2)

1

3 to set it in the rest position for ease in creation the link parameter Table 1.

This is the matrix representing frame k with respect to frame k-1. All measurements are with repect to frame k-1. For frame k-1 with respect to frame k, just invert the matrix

(2) [2],[3][6]. This matrix inverse is easy as inverse of rotation is transpose of rotation matrix. The end-effector with respect to base frame becomes

oAn = 0A1 1A2 2A3 n-1An (3)

The Linkk frame matrix with respect to universal base frame can be written as

= [ ] or [ ] (4)

= [ ] or [ ] (4)

0A

k 0 1

0 0 0 1

For notation, the symbol 0Rk is the rotational part of the matrix, onk , ook , oak , are unit vectors forming a right handed system of orthogonal vectors. This represents the orientation of framek with respect to base, and 0Pk is the position of the end-effectorwith respect to the base frame.

Inverse Kinematics and velocity are described briefly

Figure 2. Jumbo Drilling Robot [7], Joint Constraint and Home position

as:

Forward kinematics is a chain of frames or A-matrices k-1Ak

for k = 1,n, resulting in oAn .

In order to have complete and simple solution to inverse kinematics for robot manipulators, it is solved analytically to determine closed form solution. Inverse kinematics determines the chain of frames or A-matrices k-1Ak for k = 1,n from the end-effector, oAn . Inverse kinematics is based on a standard trigonometric equation a cos() + b sin() = c

Figure 3. Rest position the Jumbo Drilling Robot.

Table 1. Link parameter table for first three links Jumbo Drilling Robot Spherical Arm

The Joint mobility is constrained by mechanical limitations or physical stops. Each joint limits are tied to prismatic and revolute joint that are valid irrespective of the configuration of the remaining links.

There are three trigonometric functions: sin , cos, tan. The inverse functions sin-1(value), tan-1( value) determine

angle between –

2

and 2

whereas cos-1(value) angle is

Using the succint notation for cosines and sines, sk = sin(k), ck = cos(k), the A matrices are

c1 0 s1 a1c1

between 0 to . For robots, the rotation angle for the arm

joints is 270 degrees and the wrist joints is 180 degrees [Fig. 4], 360 degree for drilling hole [2]. The sin-1(value) , cos- 1(value) , tan-1(value) are insufficient to cover all quadrants for accuracy. Some authors use the inverse functions erroneously [9], [10]. For example, we explore solving the following equations for 1, 2, d3, from the position (6) of robot arm

Px = -c1(s2d3-a1) (9)

s

s

0A1 = 1

0 c1

a s

a s

1 1 1A2 =

Py = -s1(s2d3-a1) (10)

Pz = c2d3 (11)

0 1 0 0

0 0 0 1

Their [10] method would solve them as follows, their approach runs into undefined values for 1, 2 on using

c 0 s 0

1 0 0 0

standard inverse atan functions

2 2 From (4,1),(4,2),(4,3)

s2 0 c2 0

2A

0 0 1 0

=

Px = -c1(s2d3-a1) (12)

3

Py = -s1(s2d3-a1) (13)

0 1 0 0

0 0 0 1

0 1 0 d3

0 0 0 1

Pz = c2d3 (14)

The authors from [9], [10] simplify

= 1(231)

= 1

(15)

1(231) 1

1. Forward Kinematics

In forward kinematics, we multiply these three a- matrices, use trigonometry to combine complex expressions into simpler ones. There is a unique solution for the end-effector orientation R and position P, see equation (6).

The end-effector presentation uses the simple notation for cosines and sines, sk = sin(k), ck = cos(k), and we have

12 12 1 1(23 1)

0A3 = [12 12 1 1(23 1)] (6)

2 2 0 23

0 0 0 1

and calculate 1

1 = atan( ) (16)

There are several issues with this calculation 1. Cancelling

negative sign creates discrepancy of 180 degrees in the angle, 2. Cancelling (s2d3-a1) is not practical if it is zero. In that case division by zero is undefined. That is, if = 0, then 1 is undefined. In fact it may be well defined. We propose to use atan2, tan inverse with two parameters, if px

= 0, py 0, then

1 = atan2(py,px) = atan2(py,0) = Â± (17)

2

Secondly [9] calculates

The first three columns form a right handed system of mutually orthogonal unit vectors of the rotation matrix and the 4th column is the position vector.

2 = atan (

(111)

) (18)

2. Detour to Trigonometry

Inverse Kinematics depends on the ability to solve trigonometric equations. Basically, five trigonometric equations were required to solve the inverse kinematics problems [8]. It is not necessary to solve 5 different equations. Only one trigonometric equation is sufficient to resolve all cases. Any special case follows from the same equation. Consequently, in robotics, inverse kinematics amounts to solving the trigonometric equation

a Cos + b Sin = c. (7)

All other cases involving trigonometric functions are special case of this equation

a Cos + b Sin = c, c2 a2 + b2 (8)

which is again undefined if the expression (a1-c1Px-s1Py) =

0. It leads to both 1 and 2 as undefined as well as erroneous.

In fact, the solution can be well defined if we avert any such discrepancies involving divisions by zero that create undefined expressions. First, we use atan2 function that prevents such discrepancies of undefined terms in the calculation of angles, Second we avoid divisions by zero as follows.

In general, using atan2, the solution of

a Cos + b Sin = c, c2 a2 + b2 (19)

is

= atan2(b,a) atan2(2 + 2 2, c) (20)

There is no division by zeros, or undefined expression.

Very frequently we come across the equation

a Cos + b Sin = 0 (21)

with a2 + b2 0. Here too, the solution becomes

= atan2(b,a) atan2(2 + 2, 0) or

= atan2(b,a) atan2(1, 0) or

a1-c1Px-s1Py = s2d3 (29)

and from (4,3)

Pz = c2d3 (30)

Multiply (25) Pz = c2d3 with s2 and (29) a1-c1Px-s1Py = s2d3 with c2 and subtract, we get

s2Pz (a1-c1Px-s1Py) c2 = 0 (31)

= atan2(b,a)

2

(22)

s2Pz + (c1Px+s1Py- a1) c2 = 0 (32)

There are two solutions

2 = atan2(Pz, c1Px+s1Py- a1) Â±

2

Recall [9], [10] calculate 2 as atan( )

(111)

which is undefined when a1-c1Px-s1Py = 0.

(33)

Finally again from (4,1),(4,2),(4,3), for each value of 1 , 2 we get unique value for d3. Multiply Pz = c2d3 with c2 and a1-c1Px-s1Py = s2d3 with s2 and add, we get a unique

d3 = (a1-c1Px-s1Py)s2 + Pz c2 (34) Thus, if only arm position is specified, then there are four

possible solutions (1, 2, d3)

Fig. 4 Illustration of axis range and position MA2000 [9]

C. Inverse Kinematics analytical

With trigonometry at hand, we are in a position to solve the three link inverse kinematics for the end-effector

12 12 1 1(23 1)

 2 2 0 23 0 0 0 1
 2 2 0 23 0 0 0 1

0A3 = [12 12 1 1(23 1)] (23)

Depending on the robot, the inverse kinematics problem may be specified with just (a) Orientation of end-effector,

(2) Position of end-effector or (c) both the Position and Orientation of the end-effector.

C.1. Inverse Kinematics when only end-effector position is specified

If the position alone is specified, there can be multiple solutions. Let the postion P be provided for this example.

From (4,1),(4,2)

Px = -c1(s2d3-a1) (24)

Py = -s1(s2d3-a1) (25)

Eliminate the terms that may be cause of concern.

Multiply (24) by s1 and (25) by c1 and subtract

s1Px-c1Py = 0 (26)

C.2 Inverse kinematics when end-effector is completely specified.

If orientation and position are both specified, then there is a unique solution. That amounts to solving the following equations. Let the rotation matrix be R = [N, O,A] and position be P.

We use the same example and only atan2(y,x) function. From (1,3),(2,3)

s1 = – Ax, c1 = – Ay (35)

We find unique

1 = atan2(-Ax, -Ay) (36)

From (3,1),(3,2)

s2 = Nz, c2 = Oz (37)

We find unique

2 = atan2(Nz, Oz) (38)

From (4,1),(4,2)

Px = -c1(s2d3-a1) Py = -s1(s2d3-a1)

Now that 1 nd 2 are known, we multiply the equation (24) by -c1 and (25) by -s1, and add. We get

-c1Px-s1Py = s2d3-a1 (39)

a1-c1Px-s1Py = s2d3 (40)

From (4,3)

Pz = c2d3 (41)

Assembling these two equations (40),(41) we get

d3 = (a1-c1Px-s1Py) s2 + Pz c2 (42)

We find d3 is unique. Overall, by specifying both position and orientation, complete specification leads to a unique solution 1, 2, d3.

4. GEOMETRIC APPROACH TO INVERSE KINEMATICS

There are two solutions

1 = atan2(-Px, Py) Â±

(27)

2

Here we provided an efficient way to to solve invers kinematics and inverse Jacobian proble. We determine the

But [9], [10] calculate 1 = atan ( ) which will be

undefined when = 0

Reusing (4,1),(4,2) , multiply(24) by -c1 and (25) by -s1 and add, we get

-c1Px-s1Py = s2d3-a1 (28)

intermediate frame at the end of arm , link 3 or link 4 between the base frame ad end-effector of a 6 link robot. It makes it easy to solve 3 link or 4 link inverse kinemeic problem. For larger number of links, it becomes more challenging. To overcome this complexity, sometimes geometric approach is used. Six link robot inverse

Kinematics problem can be easilty solved by partitioning it into 3 + 3 [Fig. 5] or 4 + 2[Fig. 6] link inverse kinematics problems. The traditional geometric approaches require the specification, with respeect to base, of both Position and

Orientation of the end-effector [ ]. Now with

0 0 0 1

n,o,a,p geomery is used to calculate, with respect to base,

[ 0 0 0 1

[ 0 0 0 1

] , we calculate the parameters of first 3 or 4

links. Once it is done, we solve for parameters of the

remaing links by equating frame 3 or frame 4 (product of remaing A matrices) respectively to

1

[ 0 0 0 1] [ ]. 0 0 0 1

We propose to require only orientation [ ] of

the end-effector, that is sufficient, specification of the position of end-effector is not necessary[see Fig. 5, 6] for 3

+ 3 and 4 + 2 partitioning. In either case, we need these parameters at the end of 3 or 4 links [Fig. 5, 6]. Once, the frame at the end of arm [ ] is determined, it can be used to solve the rest of the problem, as we saw in anaylcal solutions above. The orientation and position of link three or four can be determined without the position specification of the end-effector. All other steps are similar to previous paragraph. More details on comparison of two approaches follow.

In Fig. 5 ( 3 + 3 case), visually the orientation of frame 3 is the same as frame 6, so the frame 3 orientation is the same as specified for the end-effector. Standard approach to get the postion of frame 3 with respect to base is as follows. Starting from P, we can go a4 units in direction n; a6 units along n, it gives the (a4 + a6) along n and d6 along a. This means P + (a4 + a6) n + d6 a = p, which is given. Therfore the postion of frame 3 origin is P = p – (a4 + a6) n- d6 a , completely known in terms of base coordinate system.

In figure 6 (4 + 2 case), visually the orientation of frame 4 is such that N is the same as n, A is o, O is -a. Standard approach get the postion of frame 4 with respect to base frame is as follows: starting from P, we can go a6 units along the direction n; and d6 along a. This means P + a6 n + d6 a = p, which is given. Therfore the postion of frame 4 origin is P

= p – a6 n- d6 a.

Figure 5. Partition of 6 links into 3 + 3 for inverse kinematics. Base frame expressed in terms of n, o, a.

Figure 6. Partition of 6 links into 4 + 2 for inverse kinematics. Base frame expressed in terms of n, o, a.

2. More Explantion on New Approach

We propose that specification of p is not necessary. P can be directly determined without the knowledge of p. It is observed that [N,O,A] is known, base frame orientation becomes [n,-a, o], Fig. 5, Fig. 6.

In Fig. 5, we can step a1 units along n, d3 along o resulting in P = a1 n + d3 o. In Fig. 6, we can step a1 along n, d3 along o, a4 along n resulting P = (a1 + a4) n + d3 o. In both the cases see Figure 5 and 6, we have determined that the position is P = a1 n + d3 o; orientation is [N, O, A] = [n, o, a] in Fig. 5 and the position is P = (a1 + a4) n + d3 o; orientation is [N, O, A] = [n, -a, o] in Fig. 6. This shows that arm position P can be constructed without specifying p.

5. HEURISTIC APPROACH AND JUSTIFICATION FOR INVERSE JACOBIAN,

The paper presents a heuristic method to solve inverse Jacobian for arbitrary number of degrees of freedom (DOF) robot. The Jacobian inverse computation takes long time, whereas first matrices must be multiplied, then inverse of the product is computed and then inverse is multiplied with the transpose of the original to derive the generalized inverse. We will show how this is simpler and will apply this approach to both under- and over-determined systems of velocity equations. We also provide a better method for solving inverse Jacobian problems that avoid singularity in the solutions.

The Jacobian equation

[ ] = J [ ] (43)

can be mapped into simple equation

B = AX (44)

With minimum computation effort and easily comprehensible approach, we proceed as follow. Since A is mÃ—n, there are three possibilities

1. m < n, 2. m = n, 3. m > n. (45)

For m = n, we can directly resort to inverse of A to arrive at

X = A-1 B.

For m<n or m>n we compute product of A with its transpose in such way that the product is a square matrix of smaller order. For example, if m>n , then ATA is order nÃ—n, and if m<n, then AAT is smaller order mÃ—m. Then we use the inverse of these matrices to arrive at the solution. Without going through the exercise of minimization, we can determine the solution. Vector differentiation is not trivial used for computing generalized inverses when m<n and m>n. Moreover, for mÃ—n matrix A with m>n , since ATA is a smaller square nÃ—n matrix, inverse is also square nÃ—n matrix. To derive the solution, we simply multiply nÃ—n matrix (ATA)-1 and nÃ—m matrix AT, then nÃ—m matrix (ATA)- 1 AT and mÃ—1 vector B in that order resulting in solution to the equation,

X = (ATA)-1AT B. (46)

For efficiently , we can multiply nÃ—m matrix AT and mÃ—1 vector B first, then resulting mÃ—1 vector ATB, and nÃ—n matrix (ATA)-1 yielding

X = (ATA)-1 ATB. (47)

Now we compare the computing effort in solving the equations if AT is multiplied with (ATA)-1 first or when AT is multiplied to B first in the case of m>n. Inversion is the same in both the cases, the number of multiplications is:

first case (46): m(n2 + (m + 1)n)> m(n2 + (n +

1)n)

second case(47): m(n2 + 2n)

Clearly second approach is more efficient. Since we are interested in minimum error or minimum norm solution, we will create the solution directly, heuristically minimal error is achieved by using a matrix of smaller dimension for inversion. Smaller the matrix, smaller the error in computation and representation. It also leads to conceptually comprehensible and easy to grasp the solution.

1. Underdetermined system, m<n

For underdetermined systems there can be infinitely many solutions with no error of computation. However, we are interested in solution with smallest norm, in other words, a solution that is closest to the origin or a shortest distance from the origin.

We show that for m<n,

X = AT(AAT)-1 B (48)

is a least norm solution for equation

AX = B (49)

This X is closest to the origin. Traditionally, we minimize f(X) = |X|2 subject to AX = B (50)

By method of Lagrange multipliers

f(X) = XT X + T (B-AX) (51)

Differentiate with respect to and X, and equate them to zero

() = 0 yields B – AX = 0 or B = AX (52)

() = 0 yields 2X-AT = 0 53)

Solve for

(AAT)-12B = (55)

and substitute

2X-AT(AAT)-12B = 0 (56)

We get

X = AT(AAT)-1B (57)

which is the solution to the equation AX = B.

This solution can be directly and quickly determined from the equation, AX = B.

For the purpose of creating a smaller size matrix for inversion, let us write

X = ATU for some U, (58) For this U, the equation becomes

AX = AATU or B = AATU (59)

Now A AT is mÃ—m square matrix, the inverse can be used to compute U.

U = (A AT)-1 B (60)

From X = ATU, and U = (A AT)-1 B we get (61)

X = AT(A AT)-1 B (62)

which proves the correctness of our intuitive solution shown above.

2. Overdetermined system, m>n

For overdetermined systems there may not be any solution. Any approximate solution may have an error of approximation. We can settle for an approximate solution that has minimum approximation error.

We show that for m>n,

X = (ATA)-1ATB (63)

is a least square error solution for equation

AX = B (64)

We proceed to minimize f(X) = |AX-B|2 subject to AX = B.

The function is written as

f(X) = (AX-B)T(AX-B)

= (XTAT-BT) (AX-B) (65) f(X) = XTAT AX -BT AX – XTATB + BT B (66) f(X) = XTATAX 2 XTA TB + BT B (67)

Differentiate with respect to X and set it equal to zero

() = 0 yields 2 ATAX – 2 ATB = 0 (68)

ATAX – ATB = 0 ATAX = ATB or (69)

X = (ATA)-1ATB (70)

Which is a solution to the equation AX = B. This can be directly and quickly determined from

AX = B (71)

Multiply by AT

ATAX = ATB (72)

Now ATA is nÃ—n square matrix of smaller size, the inverse can be used to compute X by multiplying both sides with (ATA)-1.

X = (ATA)-1 ATB or (73)

X = (ATA)-1AT B (74)

which proves the correctness of our simplified intuitive solution computed above.

The matrix (ATA)-1AT or AT(A AT)-1 is the pseudo inverse or generalized inverse of A. It is synonymously denoted with ay of several symbols such as: A +, A-, A #, A I etc.

Which is

2AX-AAT = 0 or

2B-AAT = 0 or 2B = AAT (54)

6. CONCLUSION

Inverse Kinematics and Inverse Jacobian are two of the problems in robot motion consideration. We have shown how to resolve these two issues with an intuitive, persuasive and provable approach. We gave closed form, efficient solution for both Inverse Kinematics and Inverse Jacobians. In fact, we have proved that our approach is accurate and efficient. We hope that the designers will find it useful in understanding, and retaining the description for reference in their future implementations.

7. REFERENCES

1. V. Grecu, N. Dumitru, and L. Grecu, Analysis of Human Arm Joints and Extension of the Study to Robot Manipulator, Proceedings of the international MultiConference of Engineers and Computer Scientists, Vol 2, March 18 – 20, 2009 pp:1348-1351.

2. Saeed B. Niku, Introduction to Robotics analysis, control, applications publisher: Wiley 2011 ISBN-978-0-470-60446-5

3. Craig, John J, Introduction to Robotics, Mechanics and Control by. 3rd Edition, Prentice Hall 2005.

4. C. Hua, C. and Wei-shan, Wavelet network solution for the inverse kinematics problem in robotic manipulator,Chen Zhejiang Univ Science , Vol. 7 ,No. 4 ,2006 PP: 525-529.

5. V.Sanchez,R.Gutierrez, G.Valdovinos and P.Ortega, 5-DOF manipulator simulation based on MATLAB-Simulink methodology , 2010 20th International Conference on Electronics, Communications and Computer (CONIELECOMP), 2010, pp: 295 300

6. Mohammed Z. Al-Faiz, Abduladhem A. Ali3,Abbas H.Miry, HUMAN ARM SIMULATION BASED ON MATLAB WITH VIRTUAL ENVIROMENT, IJCCCE, VOL.11, NO.1, 2011, pp. 86-96.

7. C.Y.Ho; Yao Jianchi, Design of Automated Jumbo Drilling Robot, Manipulator, IFAC Proceedings Volumes, Volume 20, Issue 8, August 1987, Pages 59-66

8. S. Kucuk dan Z. Bingul, Robot kinematics: Forward and inverse kinematics, in Industrial Robotics: Theory, Modelling and Control, InTech, 2006.

9. Mohammed Z. Al-Faiz ,Mohammed S.Saleh, International Journal of Robotics and Automation (IJRA), Volume (2) : Issue (4) : 2011, pp.256-264.

10. I. Agustian, N. Daratha, R. Faurina, A. Suandi, and Sulistyaningsih, Robot Manipulator Control with Inverse Kinematics PD- Pseudoinverse Jacobian and Forward Kinematics Denavit Hartenberg, Jurnal Elektronika dan Telekomunikasi, vol. 21, no. 1, pp. 818, 2021.