Next Article in Journal
Special Issue “Legged Robots into the Real World”
Next Article in Special Issue
CAD-Based Robot Programming Solution for Wire Harness Manufacturing in Aeronautic Sector
Previous Article in Journal
Recent Advances and Perspectives in Deep Learning Techniques for 3D Point Cloud Data Processing
Previous Article in Special Issue
Simulated and Real Robotic Reach, Grasp, and Pick-and-Place Using Combined Reinforcement Learning and Traditional Controls
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Inverse Kinematics of an Anthropomorphic 6R Robot Manipulator Based on a Simple Geometric Approach for Embedded Systems

by
Michael Anschober
1,*,
Raimund Edlinger
1,2,*,
Roman Froschauer
1 and
Andreas Nüchter
2
1
Department of Smart Automation and Robotics, University of Applied Sciences Upper Austria Campus Wels, Stelzhamerstrasse 23, 4600 Wels, Austria
2
Department of Informatics XVII: Robotics at the Julius-Maximilians Universität Würzburg, Am Hubland, 97074 Würzburg, Germany
*
Authors to whom correspondence should be addressed.
Robotics 2023, 12(4), 101; https://doi.org/10.3390/robotics12040101
Submission received: 31 May 2023 / Revised: 29 June 2023 / Accepted: 7 July 2023 / Published: 12 July 2023
(This article belongs to the Special Issue Advanced Grasping and Motion Control Solutions)

Abstract

:
This manuscript presents an efficient algorithm for solving the inverse kinematics problem of a 6R robot manipulator to be deployed on embedded control hardware. The proposed method utilizes the geometric relationship between the end-effector and the base of the manipulator, resulting in a computationally efficient solution. The approach aims to minimize computational complexity and memory consumption while maintaining the accuracy and real-time performance demonstrated by simulations and verified by experimental results on an embedded system. Furthermore, the manipulator is analyzed in terms of singularities, limits, the workspace, and general solvability. Due to the simplicity of the algorithm, a platform-independent implementation is possible. As a result, the average calculation time is reduced by a factor of five to eight and the average error is decreased by a factor of fifty compared to a powerful analytical solver.

1. Introduction

The capability to perform fine and intricate movements with its end-effector is essential for mobile manipulators, as it enables the robot to interact with its environment in a more effective and precise manner. A modular adaptable search and rescue robot [1] with dexterous manipulation capabilities can perform a variety of tasks, such as picking up objects, opening doors, and manipulating objects in tight spaces, see Figure 1. The development of such complex systems requires the integration of several different technologies, including robot kinematics, control, sensing, and actuation. The robot’s end-effector is designed to have a high degree of dexterity, with the ability to move in many different directions and orientations. An intuitive Human–Robot Interaction (HRI) approach that allows operators to precisely and accurately control end-effector movements is essential [2].
One approach to achieve dexterous manipulation is to use robotic arms or hands with multiple degrees of freedom (DOF). These types of manipulators can perform a wide range of movements, allowing the robot to interact with its environment in a more flexible manner. The kinematics of the manipulator must be carefully designed to ensure that it has the necessary range of motion and precision for the task at hand. For this purpose, manipulators with six DOF are typically used, e.g., a manipulator consisting of serial links and an end-effector connected with six revolute joints (6R). This is due, among other reasons, to the fact that a minimum of six DOF are required to reach any pose of the tool center point (TCP) within the workspace, which provides a high degree of usability and flexibility. Due to the high nonlinearity of the system, resource-intensive computational algorithms are typically executed on powerful computers to control six DOF manipulators.
For mobile robot platforms, the challenge is that either the data have to be transmitted to the robot or the control unit has to be moved to the platform. For smooth and accurate motions, a large amount of data is generated. Consequently, real-time capable transmission systems with high data transfer rates and a stable transmission link are needed. Since this condition generally cannot be guaranteed in rescue robotics, it is reasonable to perform the calculation directly on the mobile platform. In the case of lightweight rescue robots with a compact design (cf. Figure 1), it is impractical or impossible to move large and heavy computers to the platform, therefore small and lightweight embedded systems are used. However, due to the limited computing power, developing an efficient algorithm that consumes as few computing resources as possible is necessary. Further difficulties arise from the fact that the communication with the robot may be interrupted at any time. In this case, the stable behavior of the manipulator needs to be ensured. As a consequence, differential kinematic approaches to end-effector control are generally excluded, since a stop of the manipulator cannot be ensured.
For the reasons stated, an efficient inverse kinematics (IK) solver to control a custom-built manipulator attached to a mobile robot platform is presented. Based on the requirements for the solver and the geometry of the manipulator, a geometric approach for solving the IK problem is chosen. The developed algorithm is supported by experimental results, which were investigated for the accuracy and real-time capability of the solver on embedded hardware and compared to a powerful analytical kinematics solver. The modularization of the control unit in combination with the manipulator allows the manipulator to be addressed via a predefined interface. As a result, the manipulator can be controlled independently of higher-level control systems.
The presented manuscript contains the following key contributions:
  • The fundamental derivation of the forward and inverse kinematics of a 6R robot manipulator, as well as the development of the IK solver in order to resolve the joint angles based on arbitrary TCP poses. The geometric approach is an extension of the inverse kinematic solutions based on industrial robots with respect to the deployment of embedded hardware.
  • The analysis of important characteristics, which are essential for the practical application of the solver.
  • Simulation results that validate the presented solution to the IK problem and the solver accuracy with high computational power.
  • Experimental results that examine the solver accuracy and real-time capability on an embedded open-source 32-bit ARM Cortex®-M7 board with an FPU and a clock frequency of 260 MHz (https://emanual.robotis.com/docs/en/parts/controller/opencr10/ accessed on 6 July 2023).
When implementing IK on an embedded platform, there are several innovations and considerations that can improve its performance and efficiency. Compared to related work, there are some notable innovations here:
  • Efficient IK Solvers: Developing lightweight and efficient IK solvers is crucial for embedded platforms with limited computational resources.
  • Real-Time Optimization: Embedded platforms often require real-time performance for controlling robots in dynamic environments. The developed algorithm can efficiently handle constraints and optimize IK solutions within the given time constraints.
  • Memory Optimization: Embedded platforms often have limited memory resources. Optimizing data structures, algorithms, and parameter representations can help reduce memory usage.
  • Energy Efficiency: Embedded systems typically operate on limited power sources. Designing energy-efficient IK algorithms by minimizing unnecessary computations, optimizing data flow, or leveraging low-power modes when idle can help prolong the system’s battery life for mobile manipulators.
  • Platform-Independent: The platform-independent IK solver presented in this paper offers a versatile and adaptable solution for solving the inverse kinematics problem of robotic manipulators. Decoupling the solver from platform-specific considerations provides a unified framework that can be easily integrated into diverse robotic systems, regardless of the underlying computing platform.
It is important to note that the specific innovations for an IK implementation on an embedded platform may vary depending on the system’s requirements, constraints, and available hardware resources.

2. Related Work

Due to the fact that IK of robots is of central importance in robotics, it has been extensively covered. Fundamental studies of robot manipulators are described by Pieper [3]. Furthermore, it is shown that the IK problem for manipulators with six DOF and spherical wrists can be separated and thus solved geometrically. Lee and Ziegler [4] describe a geometric approach to obtain a closed-form solution for the IK of a robot with a geometric configuration equivalent to that of the PUMA robot. The IK problem is solved in two stages. First, the solution of the first three joints is derived geometrically. Second, the solution of the last three joints is calculated based on the solution of the first three joints to result in the desired orientation of the TCP. The selection of the most suitable solution among several possible solutions is realized via indicators. Lloyd and Hayward [5] provide forward kinematics, analytical IK, and Jacobian solutions of common serial manipulators.
Raghavan and Roth [6] present a solution of the IK problem for robots with general geometry. Through extensive mathematical transformations and simplifications, a univariate 16th order polynomial can be formulated, whose roots can be solved. Subsequently, all remaining unknowns can be solved. The description of the robot kinematics is compatible with the notation according to Denavit and Hartenberg [7], which is commonly used to relate coordinate systems in robotics. Manocha and Canny [8] develop an efficient algorithm to solve the system of equations stated in [6] using eigendecomposition. As a consequence, the calculation time can be significantly reduced.
IK solutions for sixteen different types of manipulators with spherical wrists are provided by Kucuk and Bingul [9]. Thus, closed-form solutions exist for all common manipulators where the last three axes intersect at one common point.
Husty et al. [10] describe an algorithm with a new approach to solving the IK problem compared to previous publications. The geometric structure of the manipulator is considered in the analysis using the study model of Euclidean displacements. Further, techniques of multidimensional geometry are applied to solve the problem. As a result, all sixteen solutions to the IK problem of serial 6R robots with general geometry can be calculated analytically and more efficiently than in [6].
An alternative and more generic description of robot geometry compared to the Denavit–Hartenberg notation [7] is described by Brandstötter et al. [11]. A corresponding forward and closed-form IK solution is stated. As a result, serial 6R robots with ortho-parallel basis and spherical wrists can be defined by seven characteristic parameters, which improves the practicability and applicability of common robots.
Khatamian [12] provides a forward and geometric IK solution to robots with an ortho-parallel basis and a spherical wrist, based on the industrial robot KUKA KR60. An alternative solution to [12], based on the industrial robot COMAU NM45, is proposed by Asif and Webb [13]. Solutions to similar approaches are presented in [14,15,16]. Krishnan and Ashok [17] describe a forward and closed-form analytical IK solution to robots with an ortho-parallel basis and a spherical wrist, based on the industrial robot ABB IRB 1200. Recent publications concentrate on solving the IK problem based on screw theory, e.g., [18,19,20]. In summary, a variety of solutions exist for calculating the IK of 6-DOF robotic manipulators, among other reasons due to their importance in the industry.
The open source framework MoveIt [21], which is integrated into ROS [22], is a widely used software for manipulation. MoveIt provides core functionalities such as motion planning, 3D perception, navigation, and control. Despite the advantages mentioned, MoveIt cannot be executed on embedded systems. The kinematics solver IKFast (http://openrave.org/docs/latest_stable/openravepy/ikfast/ accessed on 6 July 2023) [23] can be used as a plugin for MoveIt, and the source code of the closed-form analytical solution can be exported for independent use. This enables the development of the IK solver in MoveIt and a platform-independent application of the algorithm. Since the source code generated by IKFast has to be extensively adapted to be used on embedded systems and the geometry of the manipulator is not considered in the analytical solution (and thus is more extensive), a geometric IK solver with an approach similar to that proposed in the publications [12,13,15,16] is presented.

3. Fundamental Geometry

Robot kinematics is based on the geometry of the manipulator and specifies the procedure of the solver development. For this reason, the geometry is analyzed in advance. The kinematic relationships are described by homogeneous transformation matrices (cf. Equation (1)). The homogeneous transformation matrix T j i consists of the translational vector p i j i and the rotational matrix R j i (cf. Equations (2) and (3), respectively).
T j i = R j i p i j i 0 T 1
p i j i = p i j , x i p i j , y i p i j , z i
R j i = m j i n j i o j i = m j , x i n j , x i o j , x i m j , y i n j , y i o j , y i m j , z i n j , z i o j , z i
The following notation commonly used in robotics is introduced:
C i = cos φ i C i j = cos ( φ i φ j ) S i = sin φ i S i j = sin ( φ i φ j )
A common method to describe the manipulator geometry is the parameterization according to Denavit and Hartenberg [7]. The kinematics of the manipulator can be completely described by the Denavit–Hartenberg parameters a i , α i , d i , and θ i of each link i and the corresponding coordinate systems, illustrated in Table 1 and in the schematic sketch in Figure 2. The corresponding CAD drawing of the initial pose of the manipulator is depicted in Figure 3. The joint angles φ 1 φ 5 are limited by the mechanical structure of the manipulator, and the joint { 6 } is designed as a continuous joint, thus allowing infinite rotation. The transformation matrix A i i 1 is used to represent the coordinate system { i } in the coordinate system { i 1 } based on Denavit–Hartenberg parameters and allows the complete description of the kinematic chain of a rigid body system. The Denavit–Hartenberg transformation matrix is given in Equation (5).
A i i 1 = R o t ( z i 1 , θ i ) · T r a n s ( z i 1 , d i ) · T r a n s ( x i , a i ) · R o t ( x i , α i ) = cos θ i sin θ i · cos α i sin θ i · sin α i a i · cos θ i sin θ i cos θ i · cos α i cos θ i · sin α i a i · sin θ i 0 sin α i cos α i d i 0 0 0 1
In the context of the 6R manipulator, each revolute joint of the robotic arm is actuated by a Dynamixel Pro (https://emanual.robotis.com/docs/en/dxl/pro/ accessed on 6 July 2023) motor. These motors provide precise control over the joints, enabling accurate positioning and movement of the arm. For communication, the RS485 interface is used, which allows seamless integration of the arm with the embedded board. Table 2 shows the motor configuration. Each joint is connected by rigid links made of lightweight and durable materials such as aluminum or carbon fiber, ensuring strength and flexibility.
Figure 4 illustrates the workspace of a 6R manipulator, which refers to the region in 3D space where the end-effector (the tool or gripper attached to the robot) can reach and operate. The left view orientation shows the vertical and horizontal range of the workspace, which depends on the height of the arm and the orientation of the joints. The arm can reach 3D poses within a sphere with a radius of 1090 mm from the center of joint { 1 } of the manipulator. The working area from the top view of the manipulator is shown as a two-dimensional projection of the plane. The shape of the working area is influenced by the lengths of the arm segments, the range of the joints, and all physical constraints that limit the movement of the arm with a radius of 1090 mm.

4. Forward Kinematics

Forward kinematics is the determination of the position and orientation of the TCP in relation to the base depending on the joint variables. The forward transformation of serial robots, also referred to as direct transformation, results from a chain transformation of the individual transformation matrices stated in Equation (6) and provides a unique solution for serial robots. The pose of the TCP can thus be expressed by the transformation matrix T 6 0 .
T 6 0 = T 1 0 · T 2 1 · T 3 2 · T 4 3 · T 5 4 · T 6 5
Hence, the calculation of the individual transformation matrices of each link based on the joint variables and manipulator geometry is required. The transformation matrices can be obtained according to the Denavit–Hartenberg convention (cf. Equation (5)). It should be noted that the assignment of the coordinate systems according to the Denavit–Hartenberg convention does not necessarily correspond to the desired assignment of the joint coordinate systems due to the fact that they have to be assigned according to certain criteria in order to comply with the Denavit–Hartenberg convention. Due to this reason, additional transformations may be necessary to achieve the desired representation.
For manipulators with simple geometries, the transformation matrices can alternatively be obtained by translations and elementary rotations in relation to arbitrarily assigned joint coordinate systems. For this robot arm, the forward and inverse kinematics are developed according to the sketch illustrated in Figure 5. The procedure is generally the same, however, the description of the transformation matrices is simplified because they consist of only one translation and rotation instead of two.

5. Inverse Kinematics

Inverse kinematics or backward transformation describes the inverse operation of forward kinematics, i.e., the determination of the joint variables at a given pose of the TCP. In contrast to forward kinematics, inverse kinematics generally does not provide a unique solution. Thus, the robot arm is analyzed to develop an inverse solver in order to resolve the joint angles ( φ 1 φ 6 ) based on the pose of the TCP.

5.1. Geometric Derivation

Given that the axes of joints { 4 } , { 5 } , and { 6 } intersect at a single point, the IK can be solved geometrically (described by Pieper [3], cf. e.g., [24] (p. 29)). The position of the intersection point of the axes is independent of the joint angles φ 4 φ 6 and results directly from the given goal pose of the TCP. Consequently, the joint angles φ 1 φ 3 of joints { 1 } , { 2 } , and { 3 } can be calculated by geometric considerations based on the known position of the wrist point (WP, intersection point of the axes) illustrated in Figure 6. Due to the given goal orientation R 6 0 of the TCP and the calculated angles φ 1 φ 3 , the angles φ 4 φ 6 can then be resolved analytically. Mathematically, this means a separation of the problem in two sub-problems. The basic procedure is described by the following Equations (7) and (8):
p 05 0 = p 04 0 = p 02 0 ( φ 1 ) + p 25 0 ( φ 2 , φ 3 )
R 6 3 ( φ 1 , φ 2 , φ 3 ) = R 6 3 ( φ 4 , φ 5 , φ 6 ) ( R 1 0 ( φ 1 ) · R 2 1 ( φ 2 ) · R 3 2 ( φ 3 ) ) 1 · R 6 0 = R 4 3 ( φ 4 ) · R 5 4 ( φ 5 ) · R 6 5 ( φ 6 )
Hence, the calculated solutions for the joint angles φ 1 φ 6 result in the exact position p 06 0 and orientation R 6 0 of the TCP (specified by the goal pose T 6 0 ).
Due to the kinematics of the manipulator, there are eight different solutions of combinations of joint angles to obtain the desired pose of the end-effector. This results because
  • the link between the coordinate systems { 1 } and { 2 } can be oriented towards or opposite from the WP,
  • the elbow defined by the coordinate system { 3 } can be oriented upwards or downwards, and
  • the orientation of the wrist is identical every half turn of the joints assigned to the coordinate systems { 4 } and { 6 } and the corresponding angle of the joint assigned to the coordinate system { 5 } ,
hence 2 3 = 8 solutions. Based on the combinations of joint angles, the joints of the manipulator can be divided into the following three subgroups consisting of
  • joint { 1 } (shoulder);
  • joints { 2 } and { 3 } (elbow);
  • joints { 4 } , { 5 } , and { 6 } (wrist).
As a result, the IK problem can be further divided into three problems corresponding to each joint subgroup, whose solutions are described in Section 5.2.

5.2. Solution

5.2.1. Joint 1

The given goal pose of the TCP represented by the transformation matrix T 6 0 is defined in Equation (9):
T 6 0 = m n o p 0 0 0 1 = m x n x o x p x m y n y o y p y m z n z o z p z 0 0 0 1
In this specific case, superscripts and subscripts of the coordinate systems are omitted due to better readability and because of the significance of the matrix T 6 0 and its entries in the further course. Using the goal pose of the TCP  T 6 0 , the link length d 6 and the joint angle φ 6 , the transformation between base { 0 } , WP { 5 } , and TCP { 6 } can be described (cf. Equation (10)).
T 6 0 = T 5 0 · T 6 5 = T 5 0 · T r a n s ( d 6 , 0 , 0 ) · R o t ( x , φ 6 ) = T 5 0 · 1 0 d 6 0 C 6 S 6 0 0 S 6 C 6 0 0 0 0 1
Using Equation (10) and the inverse transformation matrix T 6 5 1 (cf. Equation (11)),
T 5 6 = T 6 5 1 = 1 0 d 6 0 C 6 S 6 0 0 S 6 C 6 0 0 0 0 1
the position p 05 0 of the WP with reference to the base coordinate system { 0 } can be calculated (cf. Equations (12) and (13)).
T 5 0 = T 6 0 · T 6 5 1 = T 6 0 · T 5 6 = p x d 6 · m x R 5 0 p y d 6 · m y p z d 6 · m z 0 0 0 1
p 05 0 = p 06 0 d 6 · m = p x d 6 · m x p y d 6 · m y p z d 6 · m z
Solution for Joint 1: The plane spanned by the rotation of joint angle φ 1 is orthogonal to the plane spanned by the rotation of joint angles φ 2 and φ 3 , because the axis { 1 } is perpendicular to the axes { 2 } and { 3 } . Consequently, the joint angle φ 1 can be calculated using the x and y components of the vector p 05 0 from the base to the WP (cf. Equation (14)).
φ 1 = atan2 ( p 05 , y 0 , p 05 , x 0 ) towards WP atan2 ( p 05 , y 0 , p 05 , x 0 ) away from WP

5.2.2. Joints 2 and 3

Using the calculated joint angle φ 1 , the transformation matrix T 2 0 from the base { 0 } to coordinate system { 2 } can be calculated via forward transformation (cf. Equations (15)–(17)).
T 1 0 = T r a n s ( 0 , 0 , d 1 ) · R o t ( z , φ 1 ) = C 1 S 1 0 0 S 1 C 1 0 0 0 0 1 d 1 0 0 0 1
T 2 1 = T r a n s ( a 1 , 0 , 0 ) · R o t ( y , φ 2 ) = C 2 0 S 2 a 1 0 1 0 0 S 2 0 C 2 0 0 0 0 1
T 2 0 = T 1 0 · T 2 1 = a 1 · C 1 R 2 0 a 1 · S 1 d 1 0 0 0 1
The translation vector p 02 0 (cf. Equation (18)) results directly from the transformation matrix T 2 0 .
p 02 0 = a 1 · C 1 a 1 · S 1 d 1
From the vectors p 02 0 and p 05 0 (cf. Equations (18) and (13), respectively), the vector p 25 0 can be calculated (cf. Equation (19)), which represents the translation between the coordinate system { 2 } and the coordinate system { 5 } of the WP.
p 25 0 = p 05 0 p 02 0 = p x d 6 · m x a 1 · C 1 p y d 6 · m y a 1 · S 1 p z d 6 · m z d 1
The vectors p 23 2 and p 35 3 can be set up by geometrical considerations (cf. Equations (20) and (21), respectively).
p 23 2 = a 2 0 0
p 35 3 = d 4 0 a 3
The vectors p 25 , p 23 , and p 35 form a triangle whose vertices are located at the origins of the coordinate systems { 2 } , { 3 } , and { 5 } (WP). The square of the vector length is obtained via the vector norm (cf. Equations (22)–(24)). Subsequently, the angles φ 2 and φ 3 can be calculated using the law of cosines.
p 25 2 = p 25 0 2 = p 25 , x 0 2 + p 25 , y 0 2 + p 25 , z 0 2 = ( p x d 6 · m x a 1 · C 1 ) 2 + ( p y d 6 · m y a 1 · S 1 ) 2 + ( p z d 6 · m z d 1 ) 2
p 23 2 = p 23 2 2 = a 2 2
p 35 2 = p 35 3 2 = d 4 2 + a 3 2
To calculate the elevation angle β 2 , the vector p 25 0 is transformed into the coordinate system { 1 } using the inverse of the rotation matrix R 1 0 (cf. Equation (25)).
p 25 1 = R 1 0 1 · p 25 0 = R 0 1 · p 25 0
The calculation of the elevation angle β 2 can thus be done via the x and z components of the vector p 25 1 , with the z component remaining unchanged by the transformation from the base coordinate system { 0 } to coordinate system { 1 } (cf. Equations (26) and (27)).
p 25 , x 1 = p x · C 1 + p y · S 1 a 1 d 6 · ( m x · C 1 + m y · S 1 )
p 25 , z 1 = p 25 , z 0 = p z d 6 · m z d 1
Solution for Joint 2: The solution for φ 2 can be calculated based on the elevation angle β 2 and the triangular angle α 2 (cf. Equations (28) and (29), respectively).
β 2 = atan2 ( p 25 , z 1 , p 25 , x 1 )
α 2 = arccos p 23 2 + p 25 2 p 35 2 2 · p 23 · p 25
Due to geometrical considerations, the following two solutions for the joint angle φ 2 result (cf. Equation (30)):
φ 2 = π α 2 β 2 elbow up π + α 2 β 2 elbow down
Solution for Joint 3: Likewise, the offset angle β 3 and the triangular angle α 3 can be calculated (cf. Equations (31) and (32), respectively).
β 3 = arctan p 35 , z 3 p 35 , x 3 = arctan ( a 3 d 4 )
α 3 = arccos p 23 2 + p 35 2 p 25 2 2 · p 23 · p 35
The angles α 3 and β 3 result in the joint angle φ 3 (cf. Equation (33)):
φ 3 = α 3 β 3 elbow up 2 π α 3 β 3 elbow down

5.2.3. Joints 4, 5, and 6

Based on the joint angles φ 1 φ 3 , the rotation matrix R 3 0 can be calculated by forward transformation (cf. Equation (34)).
R 3 0 = R 1 0 · R 2 1 · R 3 2 = R o t ( z , φ 1 ) · R o t ( y , φ 2 ) · R o t ( y , φ 3 )
Subsequently, the rotation matrix R 6 3 can be determined using the given orientation matrix of the TCP  R 6 0 (cf. Equation (35))
R 6 3 = R 3 0 1 · R 6 0 = R 3 0 T · R 6 0 = R 0 3 · R 6 0
to obtain the corresponding matrix entries (cf. Equation (36)).
m 6 , x 3 = ( m x · C 1 + m y · S 1 ) · C 23 m z · S 23 m 6 , y 3 = m y · C 1 m x · S 1 m 6 , z 3 = ( m x · C 1 + m y · S 1 ) · S 23 + m z · C 23 n 6 , x 3 = ( n x · C 1 + n y · S 1 ) · C 23 n z · S 23 o 6 , x 3 = ( o x · C 1 + o y · S 1 ) · C 23 o z · S 23
In addition, the rotation matrix R 6 3 can also be derived via forward transformation using the joint angles φ 4 φ 6 (cf. Equation (37)).
R 6 3 = R o t ( x , φ 4 ) · R o t ( y , φ 5 ) · R o t ( x , φ 6 ) = C 5 S 5 · S 6 S 5 · C 6 S 4 · S 5 C 4 · C 6 S 4 · C 5 · S 6 C 4 · S 6 S 4 · C 5 · C 6 C 4 · S 5 S 4 · C 6 + C 4 · C 5 · S 6 S 4 · S 6 + C 4 · C 5 · C 6
Equations (38)–(41) can be set up by transforming the entries of the rotation matrix R 6 3 stated in Equation (37).
m 6 , y 3 m 6 , z 3 = S 4 · S 5 C 4 · S 5 = tan φ 4
m 6 , y 3 2 + m 6 , z 3 2 = n 6 , x 3 2 + o 6 , x 3 2 = sin φ 5
m 6 , x 3 = cos φ 5
n 6 , x 3 o 6 , x 3 = S 5 · S 6 S 5 · C 6 = tan φ 6
Solution for Joints 4, 5, and 6: As a result, the joint angles φ 4 φ 6 can be calculated by substituting the matrix entries (cf. Equation (36)) into Equations (38)–(41) to obtain the solutions stated in Equations (42)–(44). The solution is derived from the calculation of Euler angles from rotation matrices (cf. e.g., [24] (p. 13)).
φ 5 = atan2 m 6 , y 3 2 + m 6 , z 3 2 , m 6 , x 3 wrist orientation 1 atan2 m 6 , y 3 2 + m 6 , z 3 2 , m 6 , x 3 wrist orientation 2
φ 4 = atan2 m 6 , y 3 sin φ 5 , m 6 , z 3 sin φ 5
φ 6 = atan2 n 6 , x 3 sin φ 5 , o 6 , x 3 sin φ 5
As a consequence of the Euclidean norm in the function atan2() in Equation (42), the joint angle φ 5 is always positive within the range [ 0 , π ] for wrist orientation 1 and always negative within the range [ π , 0 ] for wrist orientation 2. Hence, the calculation of both joint angle combinations can be separated, resulting in alternative Equations (45)–(47). Due to the independence of the equations, the joint angles φ 4 φ 6 can be calculated independently.
φ 4 = atan2 ( m 6 , y 3 , m 6 , z 3 ) wrist orientation 1 atan2 ( m 6 , y 3 , m 6 , z 3 ) wrist orientation 2
φ 5 = atan2 m 6 , y 3 2 + m 6 , z 3 2 , m 6 , x 3 wrist orientation 1 atan2 m 6 , y 3 2 + m 6 , z 3 2 , m 6 , x 3 wrist orientation 2
φ 6 = atan2 ( n 6 , x 3 , o 6 , x 3 ) wrist orientation 1 atan2 ( n 6 , x 3 , o 6 , x 3 ) wrist orientation 2
An equivalent solution for the joint angle φ 5 , stated in Equation (48), can alternatively be derived from Equation (40), however, the result is generally slightly less accurate than the solution in Equation (46) (verified in Section 6.6).
φ 5 = arccos ( m 6 , x 3 ) wrist orientation 1 arccos ( m 6 , x 3 ) wrist orientation 2

6. Experiments and Results

In this section, the procedure to validate the IK solution, presented in Section 5, is described and further improvements are stated.

6.1. Evaluation and Test Procedure

To verify the solution of the inverse solver, a test procedure is defined and implemented in MATLAB®. The test procedure (depicted in Figure 7) is as follows:
  • Import the manipulator geometry and set the parameters described in Table 1.
  • Random generation of a combination of joint angles φ 1 φ 6 ( φ c u r r ) within the limits.
  • Forward transformation of the generated joint angle combination to obtain the pose of the TCP represented by the transformation matrix T 6 0 .
  • Inverse transformation of the transformation matrix T 6 0 to resolve for eight different combinations of joint angles φ 1 φ 6 ( φ inv ).
  • Selection of the combination of joint angles which causes the smallest angular deviation as the solution ( φ i n v ).
  • Analysis of the IK solution in comparison with the generated joint angles and the remaining joint angle combinations.

6.2. Selection of the Solution

As described in Section 5.1, a maximum of eight different solutions to the IK problem exist, each leading to the desired pose of the TCP. However, when using the solver in practice, only one of the eight solutions is important. It is reasonable to select the solution which results in the smallest angular difference with respect to the current joint angles. As a consequence, the joints of the manipulator move at a minimum to reach the goal pose.
The minimum joint motion can be determined by several approaches. The most straightforward approach is the one shown in Figure 7, where all eight possible combinations of joint angles are calculated and then the most suited solution is selected. Alternatively, the minimum joint motion can be obtained by considering the three subgroups of joints (described in Section 5.1). It is assumed that a minimum motion of each subgroup leads to a minimum total motion of the manipulator, which means that all subgroups can be considered separately. Since the smallest angle difference can be determined after each solved subgroup of joint angles, the desired combination of joint angles can be narrowed down after each subgroup. This approach corresponds to a fusion of the components Inverse Kinematics and Joint Angle Selection (cf. Figure 7). As a result, the desired combination of joint angles remains and moreover, the calculation effort can be considerably reduced.
The angular difference Δ φ j for each possible joint angle combination j of a subgroup can be calculated using Equation (49), where i n v , i represents the IK solutions of joint { i } and φ c u r r , i represents the current angle of joint { i } . The joint indices k and l define the subgroup of joints, which is a subset of { 1 , , 6 } .
φ j = i = k l | φ i n v , i , j φ c u r r , i |
Subsequently, the joint angle combination j with the smallest angular difference φ j is selected as the solution.
In some cases, joint angle combinations can be excluded in advance, because, e.g., the joint limits are exceeded. For instance, in the case of the manipulator analyzed here, the elbow-down configurations are excluded due to the limits of joint { 3 } , i.e., only the elbow-up configurations of joint angles are relevant. Moreover, for combinations where the joint angles can be calculated independently, such as in Equations (45)–(47), a single joint angle may be sufficient to select the solution, which eliminates the need to calculate all joint angles of the subgroup. Hence, the calculation time may be further reduced.

6.3. Conversion between Singleturn and Multiturn Range

For the calculation with joint angles, it is essential that angles can be transformed between multiturn range [ , ] and singleturn range (e.g., [ π , π ] ) while maintaining the position of the joint. These transformations are particularly necessary when using continuous rotational joints with angles exceeding the singleturn range. Moreover, the conversion to a specific range may be required after the addition or subtraction of joint angles, e.g., when calculating the joint angles of the elbow due to the angular offset β 3 or when calculating the difference of two joint angles in order to select a solution. Hence, the functions wrap() and unwrap() are defined. The function wrap() is used to convert the multiturn angle φ to the singleturn angle φ w r a p p e d (cf. Equation (50))
φ w r a p p e d = wrap ( φ , φ i n , φ e x ) : = φ floor φ φ i n φ e x φ i n · ( φ e x φ i n )
and the function unwrap() is used to transform the singleturn angle φ w r a p p e d to the multiturn angle φ using a reference angle φ r e f (cf. Equation (51)).
φ = unwrap ( φ w r a p p e d , φ r e f , φ i n , φ e x ) : = φ r e f + wrap ( φ w r a p p e d φ r e f , φ i n , φ e x )
In order to convert the angles, the range is set by means of the included angle φ i n and the excluded angle φ e x . For example, by setting φ i n to π and φ e x to π , the angles are converted to the range [ π , π ] . The reference angle φ r e f limits the possible multiturn range with the result that there remains only one valid solution. Hence, the maximum angular difference between the multiturn angle φ and the reference angle φ r e f is π . As a reference angle, e.g., the current joint angle φ c u r r can be used. The functional principle of the conversion is depicted in Figure 8.

6.4. Singularities and Limits

Joint angles exist where the manipulator loses DOF due to its geometry, and these are commonly referred to as singularities of the manipulator. For singular configurations, there are often an infinite number of joint positions that lead to the same TCP pose. Furthermore, jump discontinuities of joint positions may occur for adjacent TCP poses near singularities. Due to the stated reasons, motions in the vicinity of singular configurations can cause unintended behavior, such as sudden joint rotations, and thus have to be considered separately.
In order to prevent critical behavior, the singularities of the manipulator have to be identified. Despite the fact that solutions for IK solvers based on geometric methods can be identified analytically, the calculation is performed numerically using computer-based systems. This causes the numerical errors to increase as the singular configurations are approached and can further lead to inaccurate results or, at worst, invalid floating point numbers (NaN) in the immediate vicinity of the singularities. If the required accuracy of the solver is known in advance, the numerical range for the identification of the singularities can be set more generously and thus the calculation of invalid solutions can be prevented.
The singularities and limits of the manipulator are considered separately for each of the three subgroups described in Section 5.1.

6.4.1. Joint 1

The singularity of joint { 1 } , also referred to as shoulder singularity, occurs when the WP lies on the rotation axis of joint { 1 } . A special case arises if, in addition, the rotation axes of joints { 1 } and { 6 } coincide. In the case of a shoulder singularity, the x and y components of the vector p 05 0 , stated in Equation (13), are zero and hence the joint angle φ 1 is undefined and cannot be resolved (cf. Equation (14)). In order to avoid an undefined angle, a numerically small threshold is defined which is used to specify a small area on the x y -plane around the rotation axis of joint { 1 } . If this critical area is reached, i.e., the length of the vector p 05 falls below the threshold, e.g., the current joint angle φ c u r r can be selected as the solution or, alternatively, the angle φ 1 can be calculated based on the goal orientation of the TCP and the current joint angles φ c u r r of the wrist. Subsequently, the remaining joint angles can be calculated.

6.4.2. Joints 2 and 3

The elbow singularity occurs when the WP lies on the straight line through joints { 2 } and { 3 } , i.e., the vectors p 23 and p 35 are collinear. Thereby, the elbow-up and elbow-down configurations degenerate to one singular configuration. This case can only occur when the manipulator is either fully unfolded or fully folded, although it is typically not physically possible to fully fold the manipulator due to collisions. The elbow singularity represents the exact limits of the domain of the functions stated in Equations (29) and (32). In the case of a singularity, the argument of the arc cosine function is either 1 or 1 , which corresponds to an angle of 0 or π , respectively. Hence, the angles φ 2 and φ 3 can be resolved using Equations (30) and (33), with the numerical accuracy decreasing as the singularity is approached. Further research in this context is discussed in Section 6.5.

6.4.3. Joints 4, 5, and 6

The coincidence of the rotation axes of joints { 4 } and { 6 } is referred to as wrist singularity and occurs if the joint angle φ 5 is equal to k · π ( k Z ). As a consequence, the angles φ 4 and φ 6 can no longer be resolved via the matrix entries stated in Equation (36) at singular configurations of the rotation matrix R 6 3 (cf. Equations (52) and (53)).
R 6 3 = R o t ( x , φ 4 ) · R o t ( y , 0 ) · R o t ( x , φ 6 ) = R o t ( x , φ 4 ) · I · R o t ( x , φ 6 ) = R o t ( x , φ 4 + φ 6 ) = 1 0 0 0 cos ( φ 4 + φ 6 ) sin ( φ 4 + φ 6 ) 0 sin ( φ 4 + φ 6 ) cos ( φ 4 + φ 6 )
R 6 3 = R o t ( x , φ 4 ) · R o t ( y , π ) · R o t ( x , φ 6 ) = R o t ( x , φ 4 ) · R o t ( y , π ) · R o t ( x , φ 6 ) = 1 0 0 0 cos ( φ 4 φ 6 ) sin ( φ 4 φ 6 ) 0 sin ( φ 4 φ 6 ) cos ( φ 4 φ 6 )
Singularities of the wrist can be detected via the matrix entry m 6 , x 3 , and subsequently, the corresponding joint angle φ 5 can be set. For this purpose, a numerical threshold is defined, since it is very unlikely that the matrix entry m 6 , x 3 is exactly 1 or 1 . For the matrix entry m 6 , x 3 reaching approximately 1, Equation (54) applies. For the matrix entry m 6 , x 3 reaching approximately 1 , an analogous procedure can be followed but shall not be further elaborated here since the findings are equivalent.
cos ( φ 5 ) = 1 t h r e s h o l d
Rearranging Equation (54) yields Equation (55) for calculating the numerical threshold in order to achieve a certain maximum angular deviation φ t h r e s h o l d of the joint angle φ 5 .
t h r e s h o l d = 1 cos ( φ t h r e s h o l d )
The relationship of the threshold and the maximum angular deviation φ t h r e s h o l d is shown in the double logarithmic graph depicted in Figure 9.
For example, if a wrist singularity is to be detected at an angular deviation φ t h r e s h o l d of joint angle φ 5 smaller than one-tenth of a degree, the numerical threshold is to be set to
t h r e s h o l d = 1 cos ( 0.1 ° ) 1.5 × 10 6 .
In the case of a singularity, the joint angles φ 4 and φ 6 are calculated via the matrix entries n 6 , y 3 and o 6 , y 3 (cf. Equations (57)) due to the fact that the calculation is more efficient than using the matrix entries n 6 , z 3 and o 6 , z 3 .
n 6 , y 3 = n y · C 1 n x · S 1 o 6 , y 3 = o y · C 1 o x · S 1
Subsequently, Equation (58) is used to calculate the combination of the joint angles φ 4 and φ 6 based on the rotation matrix R 6 3 at singular configurations (cf. Equations (52) and (53)).
φ 4 + φ 6 = atan2 ( o 6 , y 3 , n 6 , y 3 ) φ 4 φ 6 = atan2 ( o 6 , y 3 , n 6 , y 3 ) = atan2 ( o 6 , y 3 , n 6 , y 3 )
By rearranging Equation (58), the joint angle φ 6 can be expressed as a function of the rotation matrix R 6 3 and the angle φ 4 (cf. Equation (59)).
φ 6 = atan2 ( o 6 , y 3 , n 6 , y 3 ) φ 4 if m 6 , x 3 > ( 1 t h r e s h o l d ) atan2 ( o 6 , y 3 , n 6 , y 3 ) + φ 4 if m 6 , x 3 < ( 1 + t h r e s h o l d )
Hence, the joint angle φ 6 can be calculated if the joint angle φ 4 is set to an arbitrary angle, e.g., the current joint angle φ 4 , c u r r . The joint angle φ 5 can be set according to the matrix entry m 6 , x 3 .

6.5. Workspace and Solvability

Due to the fact that the solver is based on an analytical algorithm, a (non-imaginary) solution exists only if the pose of the end-effector is within the workspace of the manipulator. This condition can be checked for any TCP pose. The WP, the angle φ 1 , and consequently the vector p 25 can be calculated by using the goal pose of the TCP and the geometry of the manipulator (cf. Equation (19)). Furthermore, the vectors p 23 and p 35 are known (cf. Equations (20) and (21), respectively). Since these three vectors span a triangle, with only the length of the vector p 25 being variable, the solvability constraints
p 25 p 23 + p 35 p 25 a 2 + a 3 2 + d 4 2
and
p 25 | p 23 p 35 | p 25 | a 2 a 3 2 + d 4 2 |
apply (cf. triangle inequality). Geometrically, this means that the distance between the coordinate systems { 2 } and { 5 } (WP) cannot be greater than the distance resulting from a maximum deflection of the manipulator and cannot be smaller than the difference of the links. In addition, it must be considered that although the manipulator can reach the maximum deflection, the folded position cannot be reached due to the collision of the links and the limitations of the joints. For this reason, the solvability constraint is adapted to
p 25 ( a 2 d 4 ) 2 + a 3 2 .
Consequently, if the stated constraints are not fulfilled, there is no solution to the problem.
Although no solution can be found when the vector length is not within the limits and thus the position of the WP cannot be reached, the joint angles φ 2 and φ 3 can be calculated so that the vector p 25 is oriented towards the WP. Subsequently, the remaining joint angles φ 4 φ 6 can be calculated and the resulting pose approximates the desired pose. At maximum deflection, the joint angles result in
φ 2 = π β 2 = π atan2 ( p 25 , z 1 , p 25 , x 1 ) φ 3 = π β 3 = π arctan ( a 3 d 4 )
and the folded position exhibits the joint angles
φ 2 = π β 2 arctan ( a 3 | a 2 d 4 | ) = π atan2 ( p 25 , z 1 , p 25 , x 1 ) arctan ( a 3 | a 2 d 4 | ) φ 3 = 0
(see Figure 6). This allows the calculation of the joint angles for TCP poses outside the workspace of the manipulator.

6.6. Solver Accuracy and Execution Time

The inverse solver consisting of the equations presented in Section 5 is simulated according to the test procedure described in Section 6.1. The TCP poses T 6 0 are generated by forward transformation of randomly generated joint angles φ c u r r within the range [ π , π ] . Physical constraints of the manipulator such as collisions and joint limits are ignored, as only the mathematical capabilities of the solver are examined, thus all TCP poses within the workspace of the manipulator are evaluated. In addition, all eight combinations of joint angles φ inv are calculated (cf. Figure A1 and Figure A2 in Appendix A). However, only the solution φ i n v with the smallest angle difference is considered in the evaluation, which is determined on the basis of the joint angles φ c u r r (cf. Equation (49)). In order to obtain representative results, a data set of m = 10,000 randomly generated joint angle combinations is created. Consequently, valid solutions exist for each combination. The generated data set consisting of the joint angles φ 1 φ 6 is approximately uniformly distributed (illustrated in the histograms in Figure 10).
To evaluate the capability and accuracy of the inverse solver, the following solver variants based on
are examined. The detection of singularities, described in Section 6.4, is disabled, i.e., the threshold for the detection of shoulder and wrist singularities is set to zero. The solver variants are both simulated on PC and tested with the embedded board using the implicit and explicit equations. For the explicit calculation variants, the explicit equations derived in Section 5.1 are used, whereas for the implicit calculation variants, the stated matrix multiplications are calculated and then the resulting matrix entries are substituted into the equations. To compare the presented solver with the analytical solver IKFast, the source code of IKFast is generated based on the geometry of the manipulator, exported, and adapted for use on the embedded board. In order to obtain comparable results, the identical, randomly generated data set depicted in Figure 10 is both used for the simulation and the evaluation of the inverse solver. The data are transmitted without loss of numerical precision. After executing steps 3 to 5 of the test procedure (cf. Section 6.1) on the embedded board, the results are transferred to the PC for further analysis.
The capability and accuracy of the solver are evaluated based on the maximum joint angle error φ e r r , m a x for each combination j of all n = 6 joints (cf. Equation (65)). The desired result is to exactly resolve the joint angles φ c u r r , i.e., the angular difference to the resolved joint angles φ i n v expressed with the joint angle errors φ e r r is zero and thus the desired maximum joint angle error φ e r r , m a x is zero.
φ e r r , m a x = max ( φ e r r ) = max | φ c u r r φ i n v | = max i = 1 n | φ c u r r , i φ i n v , i |
The histogram and the associated box plot of the maximum joint angle error φ e r r , m a x using the solver variant 2 with explicit equations are depicted in Figure 11 (calculated with MATLAB® on the PC) and Figure 12 (calculated with the embedded board). The histogram and the associated box plot of the maximum joint angle error φ e r r , m a x of IKFast executed on the embedded board are depicted in Figure 13 (note the different scaling of the x-axis).
For further examination, all solver variants are tested and examined. The histograms of the remaining variants 1 and 3 present the same course as the histograms shown in Figure 11 and Figure 12. Due to the fact that the histograms represent an irregular distribution with a wide value range, the percentiles of the distribution of the maximum joint angle errors φ e r r , shown in Table 3, are used for the characterization.
For further investigation, the mean calculation time t ¯ m e a n of the presented solver variants and IKFast are compared to the average joint angle error φ ¯ e r r . For this purpose, the calculation on the embedded board is repeated k = 100 times for each joint angle combination j. The required calculation time t l is measured for each repetition l with a resolution of 1   μ s . Subsequently, the average, minimum, and maximum calculation times t m e a n , t m i n , and t m a x of each joint angle combination are determined. The calculation times for variant 2 are depicted in Figure 14 and in Figure A4 and Figure A5 in Appendix B. Variant 2 is shown as a representative of all solver variants, which show similar normally distributed courses. The calculation times for IKFast are depicted in Figure 15, and the interesting area is enlarged and depicted in Figure 16.
An overall arithmetic mean value t ¯ m e a n of all m = 10,000 combinations is calculated from the arithmetic mean values t m e a n in order to obtain a reasonable parameter for the comparison of the solver variants (cf. Equation (66)).
t ¯ m e a n = 1 m j = 1 m t m e a n = 1 m · k j = 1 m l = 1 k t l
The average joint angle error is obtained via the arithmetic mean of all joint angle errors of all combinations and is calculated according to Equation (67), where the joint index i runs from one to the maximum number of joints n = 6 for each combination j running from one to the number of m = 10,000 combinations.
φ ¯ e r r = 1 m · n j = 1 m i = 1 n φ e r r , i , j = 1 m · n j = 1 m i = 1 n | φ c u r r , i , j φ i n v , i , j |
The results of the average calculation time stated in Equation (66) and the average joint angle error stated in Equation (67) are compared in Table 4.

7. Discussion

Based on the histograms and box plots depicted in Figure 11 and Figure 12, which cover the majority of the data, no difference in accuracy between the simulation and the execution on the embedded board can be observed for solver variant 2 (the histograms and box plots of the solver variants 1 and 3 show an analogous course). Moreover, it can be observed that the majority of the joint angle combinations exhibit a minor joint angle error. On the basis of the percentiles stated in Table 3, it can be determined that 99% of the joint angle combinations exhibit a maximum error below 10 13   rad (cf. 99th percentile in Table 3). In the worst case, the maximum error is less than 10 10   rad (cf. 100th percentile in Table 3). Hence, it can be concluded that the presented solution is valid and moreover that the inverse solver is sufficiently accurate.
Furthermore, it can be seen that there is no significant numerical degradation due to the different variants since all variants exhibit similar results at the examined percentiles of the maximum joint angle errors. Additionally, it can be shown that the solver provides similarly accurate results both when simulated on the PC and when deployed on the embedded board.
Compared to IKFast, the presented solver variants show significantly lower error (cf. Table 3). The error of IKFast is more than ten times as large as that of the geometric solvers over the entire range of values, and the maximum error is at least more than eighty times greater (depending on the solver variant). This trend can be confirmed when comparing the histograms and box plots depicted in Figure 11, Figure 12 and Figure 13 (note the different scaling of the x-axis).
In Table 4, it can be observed that the average joint angle error of the calculation on the embedded board is slightly greater than the average joint angle error compared to the simulation on the PC. Hence, the solutions are slightly less accurate when calculated on the embedded board than the simulation results calculated on the PC. The largest joint angle errors are caused by poses in the vicinity of singularities (cf. Figure A3 in Appendix A). Moreover, the explicitly calculated variants exhibit a higher average error than the implicitly calculated variants. In addition, the error of variant 3 is greater compared to the other variants, with variants 1 and 2 having similar joint angle errors.
Furthermore, it can be observed that the average calculation times of the explicitly calculated variants are considerably lower than those of the implicitly calculated variants (cf. Table 4). In addition, a significant reduction of the average calculation time can be achieved with variants 2 and 3 compared to variant 1, whereas the average calculation times of variants 2 and 3 are approximately the same.
The average joint angle error φ ¯ e r r shows that the presented algorithm achieves significantly more accurate results than IKFast (cf. Table 4). However, the greatest achievement becomes clear when comparing the average calculation times t ¯ m e a n . It can be shown that the presented algorithm provides a considerably faster calculation of the solution in every case. Moreover, the presented algorithm provides more predictable calculation times due to the fact that the values are located in a more condensed range (cf. Figure 14, Figure 15 and Figure 16).
In summary, the solver variant 2 exhibits a low average calculation time in combination with a low average joint angle error and therefore combines the advantages of the different solver variants. The average joint angle error can be minimized by using the implicit calculation variants at the expense of higher calculation times, otherwise, the explicit calculation variants are sufficiently accurate. The proposed algorithm is suitable for use on embedded systems due to the real-time capability (cf. Figure A5 in Appendix B) and furthermore requires significantly less computation time than alternative algorithms such as IKFast. In addition, the accuracy of the solver is noticeably better than using IKFast.

8. Conclusions and Future Work

IK is the process of determining the joint angles of a robot manipulator that correspond to a desired end-effector position and orientation. This is a fundamental problem in robotic control, as it allows the manipulator to move to a desired location in the workspace. However, solving the IK problem can be computationally expensive and time-consuming, especially for complex robot manipulators with a high number of DOF.
One solution to this problem is to use a geometric approach for solving the IK of a 6R robot manipulator. In this approach, a geometric relationship between the end-effector and the joints of the manipulator is utilized to solve the IK problem. Thus, the algorithm can be simplified compared to analytical approaches. As a result, this method is computationally efficient and can be implemented on embedded boards in real-time, which makes it suitable for a wide range of mobile robotic applications. The proposed method is based on the following steps:
  • The geometric relationship between the end-effector and the joints of the manipulator is established.
  • The end-effector pose is expressed in terms of the joint angles of the manipulator.
  • The IK problem is solved using this geometric relationship.
The effectiveness of the proposed method is demonstrated through simulation and experimental results. The simulation results show that the proposed method can accurately determine the joint angles of the manipulator for a wide range of end-effector positions and orientations. The experimental results show that the proposed method can be implemented in real-time and can be used to control the 6R robot manipulator to move to a desired location in the workspace. Compared to universal out-of-the-box solutions such as IKFast, the computation time can be significantly reduced and the accuracy of the solver can be noticeably improved. In summary, the proposed method provides an efficient and accurate solution to the IK problem for an anthropomorphic 6R robot manipulator with a spherical wrist.
The developed manipulator software framework will be crucial for developing safer and more powerful mobile manipulation tasks that can operate effectively in human environments in the future. Further algorithms which provide distance-checking capabilities that should generate collision detection and avoidance are a major goal. In addition, another significant milestone is to develop a differential kinematic algorithm in order to establish the relationship between the velocity of the TCP and the joint velocities. Furthermore, investigations concerning path and trajectory planning are to be conducted in order to achieve a smooth movement of the TCP. These measures are primarily intended to intuitively enable teleoperated control of the manipulator to perform a variety of handling and inspection tasks to be accomplished in the field of rescue robotics. As a result, safe, controllable, and efficient motions of the manipulator in space should be achieved.
In conclusion, the integration of the platform-independent IK solver presents a promising avenue for achieving efficient and optimized motion trajectories in robotic manipulators. The proposed framework provides a foundation for future research and development in the field of motion planning for robotic systems on embedded platforms.

Author Contributions

M.A. is a team member of robotics research group “Team Dynamics” at the University of Applied Sciences Upper Austria. He developed the algorithm and analyzed the data and wrote the manuscript. R.E. conducted research and experiments as an active member and team leader of the mobile robotics research group “Team Dynamics” at the University of Applied Sciences Upper Austria (https://sar.fh-ooe.at/index.php/en/robo-racing/team-dynamics accessed on 6 July 2023). He also analyzed the data, developed the key studies, and wrote the manuscript. Roman Froschauer and Andreas Nüchter advised the research and experiments of this project. All authors read and approved the final manuscript.

Funding

This research was funded by the University of Applied Sciences Upper Austria via internal grants to assist the project at the Robotics Lab.

Data Availability Statement

Not applicable.

Acknowledgments

Open Access Funding by the University for Continuing Education Krems, the University of Applied Sciences BFI Vienna and the University of Applied Sciences Upper Austria. The authors would like to thank the members of the Team Dynamics team who have enthusiastically and collaboratively supported the robotics projects of the SAR— Smart Automation Robotics research group at the Upper Austria University of Applied Sciences. This publication was supported by the Open Access Publication Fund of the University of Würzburg.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:
DOFDegrees of Freedom
IKInverse Kinematics
TCPTool Center Point
WPWrist Point

Appendix A. Simulation with MATLAB®

Figure A1. Solutions 1…4 for a randomly generated goal pose in MATLAB®.
Figure A1. Solutions 1…4 for a randomly generated goal pose in MATLAB®.
Robotics 12 00101 g0a1
Figure A2. Solutions 5…8 for a randomly generated goal pose in MATLAB®.
Figure A2. Solutions 5…8 for a randomly generated goal pose in MATLAB®.
Robotics 12 00101 g0a2
Figure A3. Solutions with the largest joint angle errors simulated in MATLAB®. Pose 1—elbow singularity; Pose 2—shoulder and wrist singularity; Pose 3—elbow singularity; Pose 4—elbow singularity; Pose 5—elbow singularity; Pose 6—shoulder singularity; Pose 7—wrist singularity; Pose 8—wrist singularity.
Figure A3. Solutions with the largest joint angle errors simulated in MATLAB®. Pose 1—elbow singularity; Pose 2—shoulder and wrist singularity; Pose 3—elbow singularity; Pose 4—elbow singularity; Pose 5—elbow singularity; Pose 6—shoulder singularity; Pose 7—wrist singularity; Pose 8—wrist singularity.
Robotics 12 00101 g0a3

Appendix B. Calculation Times of the IK Solver Executed on the Embedded Board

Figure A4. Histogram and box plot of the minimum calculation times of the solver variant 2 executed on the embedded board.
Figure A4. Histogram and box plot of the minimum calculation times of the solver variant 2 executed on the embedded board.
Robotics 12 00101 g0a4
Figure A5. Histogram and box plot of the maximum calculation times of the solver variant 2 executed on the embedded board.
Figure A5. Histogram and box plot of the maximum calculation times of the solver variant 2 executed on the embedded board.
Robotics 12 00101 g0a5

References

  1. Edlinger, R.; Nuechter, A. Marc-modular autonomous adaptable robot concept. In Proceedings of the 2019 IEEE International Symposium on Safety, Security, and Rescue Robotics (SSRR), Würzburg, Germany,, 2–4 September 2019; pp. 1–7. [Google Scholar]
  2. Edlinger, R.; Anschober, M.; Froschauer, R.; Nüchter, A. Intuitive HRI Approach with Reliable and Resilient Wireless Communication for Rescue Robots and First Responders. In Proceedings of the 2022 31st IEEE International Conference on Robot and Human Interactive Communication (RO-MAN), Napoli, Italy, 29 August–2 September 2022; pp. 75–82. [Google Scholar]
  3. Pieper, D.L. The Kinematics of Manipulators under Computer Control; Stanford University: Stanford, CA, USA, 1969. [Google Scholar]
  4. Lee, C.S.G.; Ziegler, M. Geometric Approach in Solving Inverse Kinematics of PUMA Robots. IEEE Trans. Aerosp. Electron. Syst. 1984, AES-20, 695–706. [Google Scholar] [CrossRef]
  5. Lloyd, J.; Hayward, V. Kinematics of common industrial robots. Robot. Auton. Syst. 1988, 4, 169–191. [Google Scholar] [CrossRef]
  6. Raghavan, M.; Roth, B. Kinematic analysis of the 6R manipulator of general geometry. In Proceedings of the Fifth International Symposium on Robotics Research, Tokyo, Japan, 28–31 August 1989; pp. 314–320. [Google Scholar]
  7. Denavit, J.; Hartenberg, R.S. A kinematic notation for lower-pair mechanisms based on matrices. ASME J. Appl. Mech. 1955, 22, 215–221. [Google Scholar] [CrossRef]
  8. Manocha, D.; Canny, J.F. Efficient inverse kinematics for general 6R manipulators. IEEE Trans. Robot. Autom. 1994, 10, 648–657. [Google Scholar] [CrossRef] [Green Version]
  9. Kucuk, S.; Bingul, Z. The inverse kinematics solutions of industrial robot manipulators. In Proceedings of the IEEE International Conference on Mechatronics, ICM‘04, Istanbul, Turkey, 5 June 2004; pp. 274–279. [Google Scholar] [CrossRef]
  10. Husty, M.L.; Pfurner, M.; Schröcker, H.P. A new and efficient algorithm for the inverse kinematics of a general serial 6R manipulator. Mech. Mach. Theory 2007, 42, 66–81. [Google Scholar] [CrossRef]
  11. Brandstötter, M.; Angerer, A.; Hofbaur, M. An Analytical Solution of the Inverse Kinematics Problem of Industrial Serial Manipulators with an Ortho-parallel Basis and a Spherical Wrist. In Proceedings of the Austrian Robotics Workshop, Linz, Austria, 22–23 May 2014. [Google Scholar]
  12. Khatamian, A. Solving Kinematics Problems of a 6-DOF Robot Manipulator. Int. Conf. Sci. Comput. 2015, 2, 228. [Google Scholar]
  13. Asif, S.; Webb, P. Kinematics analysis of 6-DoF articulated robot with spherical wrist. Math. Probl. Eng. 2021, 2021, 6647035. [Google Scholar] [CrossRef]
  14. Pratheep, V.; Chinnathambi, M.; Priyanka, E.; Ponmurugan, P.; Thiagarajan, P. Design and Analysis of six DOF robotic manipulator. In Proceedings of the IOP Conference Series: Materials Science and Engineering; IOP Publishing: Erode, India, 2021; Volume 1055, p. 012005. [Google Scholar]
  15. Zhang, H.; Xia, Q.; Sun, J.; Zhao, Q. A Fully Geometric Approach for Inverse Kinematics of a Six-Degree-of-Freedom Robot Arm. J. Phys. Conf. Ser. 2022, 2338, 012089. [Google Scholar] [CrossRef]
  16. Dikmenli, S. Forward & Inverse Kinematics Solution of 6-DOF Robots Those Have Offset & Spherical Wrists. Eurasian J. Sci. Eng. Technol. 2022, 3, 14–28. [Google Scholar]
  17. Krishnan, M.G.; Ashok, S. Kinematic Analysis and Validation of an Industrial Robot Manipulator. In Proceedings of the TENCON 2019—2019 IEEE Region 10 Conference (TENCON), Kochi, India, 17–20 October 2019; pp. 1393–1399. [Google Scholar] [CrossRef]
  18. Chen, I.M.; Gao, Y. Closed-form inverse kinematics solver for reconfigurable robots. In Proceedings of the Proceedings 2001 ICRA. IEEE International Conference on Robotics and Automation (Cat. No.01CH37164), Seoul, Republic of Korea, 21–26 May 2001; Volume 3, pp. 2395–2400. [Google Scholar] [CrossRef]
  19. Tab, Y.; Xiao, A. Extension of the Second Paden-Kahan Sub-problem and its’ Application in the Inverse Kinematics of a Manipulator. In Proceedings of the 2008 IEEE Conference on Robotics, Automation and Mechatronics, Chengdu, China, 21–24 September 2008; pp. 379–381. [Google Scholar] [CrossRef]
  20. Chen, Q.; Zhu, S.; Zhang, X. Improved Inverse Kinematics Algorithm Using Screw Theory for a Six-DOF Robot Manipulator. Int. J. Adv. Robot. Syst. 2015, 12, 140. [Google Scholar] [CrossRef]
  21. Chitta, S. MoveIt!: An introduction. In Robot Operating System (ROS) The Complete Reference (Volume 1); Springer: Cham, Switzerland, 2016; pp. 3–27. [Google Scholar]
  22. Quigley, M.; Conley, K.; Gerkey, B.; Faust, J.; Foote, T.; Leibs, J.; Wheeler, R.; Ng, A.Y. ROS: An open-source Robot Operating System. In Proceedings of the ICRA Workshop on Open Source Software, Kobe, Japan, 12–17 May 2009; Volume 3, p. 5. [Google Scholar]
  23. Diankov, R. Automated Construction of Robotic Manipulation Programs. Ph.D. Thesis, Carnegie Mellon University, Robotics Institute, Pittsburgh, PA, USA, 2010. [Google Scholar]
  24. Siciliano, B.; Khatib, O. (Eds.) Springer Handbook of Robotics; Springer Handbooks; Springer: Berlin/Heidelberg, Germany, 2016. [Google Scholar] [CrossRef]
Figure 1. A chain-driven rescue robot equipped with a 6-DOF manipulator for dexterous manipulation tasks, e.g., turning valves or reading pressure gauges.
Figure 1. A chain-driven rescue robot equipped with a 6-DOF manipulator for dexterous manipulation tasks, e.g., turning valves or reading pressure gauges.
Robotics 12 00101 g001
Figure 2. Schematic sketch of the robot arm with coordinate systems assigned according to the Denavit–Hartenberg convention, φ 1 φ 6 = 0 rad .
Figure 2. Schematic sketch of the robot arm with coordinate systems assigned according to the Denavit–Hartenberg convention, φ 1 φ 6 = 0 rad .
Robotics 12 00101 g002
Figure 3. CAD drawing of the initial pose of the robot arm, φ 1 φ 6 = 0 rad .
Figure 3. CAD drawing of the initial pose of the robot arm, φ 1 φ 6 = 0 rad .
Robotics 12 00101 g003
Figure 4. The workspace of the 6R manipulator from left and top view orientation (all dimensions are in millimeters).
Figure 4. The workspace of the 6R manipulator from left and top view orientation (all dimensions are in millimeters).
Robotics 12 00101 g004
Figure 5. Schematic sketch of the robot arm with assigned coordinate systems used for forward and inverse kinematics, φ 1 φ 6 = 0 rad .
Figure 5. Schematic sketch of the robot arm with assigned coordinate systems used for forward and inverse kinematics, φ 1 φ 6 = 0 rad .
Robotics 12 00101 g005
Figure 6. Schematic sketch of the robot arm (coordinate systems {0}…{5}) representing the inverse solver functionality, φ 1 φ 3 arbitrary.
Figure 6. Schematic sketch of the robot arm (coordinate systems {0}…{5}) representing the inverse solver functionality, φ 1 φ 3 arbitrary.
Robotics 12 00101 g006
Figure 7. Overview of the evaluation and its main components for the IK solver.
Figure 7. Overview of the evaluation and its main components for the IK solver.
Robotics 12 00101 g007
Figure 8. Conversion between wrapped and unwrapped joint angles, φ i n = π , φ e x = π .
Figure 8. Conversion between wrapped and unwrapped joint angles, φ i n = π , φ e x = π .
Robotics 12 00101 g008
Figure 9. Relationship of the threshold angle φ t h r e s h o l d in and the numerical threshold.
Figure 9. Relationship of the threshold angle φ t h r e s h o l d in and the numerical threshold.
Robotics 12 00101 g009
Figure 10. Histograms of the randomly generated joint angles φ 1 φ 6 representing the frequency per bin of the joint angles in radian (360 bins in total).
Figure 10. Histograms of the randomly generated joint angles φ 1 φ 6 representing the frequency per bin of the joint angles in radian (360 bins in total).
Robotics 12 00101 g010
Figure 11. Histogram and box plot of the maximum joint angle errors φ e r r , m a x of the solver variant 2 with explicit equations simulated with MATLAB® on the PC.
Figure 11. Histogram and box plot of the maximum joint angle errors φ e r r , m a x of the solver variant 2 with explicit equations simulated with MATLAB® on the PC.
Robotics 12 00101 g011
Figure 12. Histogram and box plot of the maximum joint angle errors φ e r r , m a x of the solver variant 2 with explicit equations executed on the embedded board.
Figure 12. Histogram and box plot of the maximum joint angle errors φ e r r , m a x of the solver variant 2 with explicit equations executed on the embedded board.
Robotics 12 00101 g012
Figure 13. Histogram and box plot of the maximum joint angle errors φ e r r , m a x of IKFast executed on the embedded board.
Figure 13. Histogram and box plot of the maximum joint angle errors φ e r r , m a x of IKFast executed on the embedded board.
Robotics 12 00101 g013
Figure 14. Histogram and box plot of the mean calculation times of the solver variant 2 executed on the embedded board.
Figure 14. Histogram and box plot of the mean calculation times of the solver variant 2 executed on the embedded board.
Robotics 12 00101 g014
Figure 15. Histogram and box plot of the mean calculation times of IKFast executed on the embedded board.
Figure 15. Histogram and box plot of the mean calculation times of IKFast executed on the embedded board.
Robotics 12 00101 g015
Figure 16. Enlarged interesting area of the histogram and box plot of the mean calculation times of IKFast executed on the embedded board from Figure 15.
Figure 16. Enlarged interesting area of the histogram and box plot of the mean calculation times of IKFast executed on the embedded board from Figure 15.
Robotics 12 00101 g016
Table 1. Denavit–Hartenberg parameters of the first six DOF from the base to the TCP of the gripper.
Table 1. Denavit–Hartenberg parameters of the first six DOF from the base to the TCP of the gripper.
LinkJoint Variable θ i [ rad ] d i [ m ] a i [ m ] α i [ rad ]Range
1 φ 1 φ 1 0.0810.077 π 2 ±90
2 φ 2 φ 2 π 00.520 π 0 – 180
3 φ 3 φ 3 π 2 00.066 π 2 0 – 180
4 φ 4 φ 4 0.4090 π 2 ±180
5 φ 5 φ 5 00 π 2 ±90
6 φ 6 φ 6 0.18000 ±
Table 2. Motor configurations of the 6-DOF manipulator.
Table 2. Motor configurations of the 6-DOF manipulator.
JointMotorMotor Resolution [pulse/rev]
1H54-100-S500-R501,923
2H54-200-S500-R501,923
3H54-200-S500-R501,923
4H54-100-S500-R501,923
5H42-20-S300-R303,751
6H42-20-S300-R303,751
Table 3. Percentiles of the distribution of the maximum joint angle errors φ e r r .
Table 3. Percentiles of the distribution of the maximum joint angle errors φ e r r .
Solver VariantPercentile [ 10 16   rad ]
0th (Min.)5th25th50th (Med.)75th95th99th100th (Max.)
PCimplicit10.02.24.48.917.8125.4770.5107,642.8
20.02.24.48.917.8125.5770.5107,642.8
30.02.24.48.917.8127.4783.3107,642.8
explicit10.02.24.47.813.3102.1532.8309,219.3
20.02.24.48.313.3102.1532.8309,219.3
30.02.24.48.915.5106.6537.3309,219.3
Embedded Boardimplicit10.02.24.48.914.4102.1562.3267,226.2
20.02.24.48.914.4102.1562.3267,226.2
30.02.24.48.915.5106.6579.5267,226.2
explicit10.02.24.48.914.4102.1532.8309,219.3
20.02.24.48.914.4102.1532.8309,219.3
30.02.24.48.915.5106.6532.8309,219.3
IKFast13.340.071.1113.2340.83466.120,643.524,375,716.2
Table 4. Average joint angle error φ ¯ e r r and average calculation times t ¯ m e a n on the embedded board.
Table 4. Average joint angle error φ ¯ e r r and average calculation times t ¯ m e a n on the embedded board.
Solver VariantPCEmbedded Board
φ ¯ err [ 10 15   rad ] φ ¯ err [ 10 15   rad ] t ¯ mean [ μ s ]
implicit12.863.05378.3
22.863.05313.7
32.933.13310.1
explicit13.163.19316.1
23.163.20252.0
33.243.27252.9
IKFast-153.991856.2
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Anschober, M.; Edlinger, R.; Froschauer, R.; Nüchter, A. Inverse Kinematics of an Anthropomorphic 6R Robot Manipulator Based on a Simple Geometric Approach for Embedded Systems. Robotics 2023, 12, 101. https://doi.org/10.3390/robotics12040101

AMA Style

Anschober M, Edlinger R, Froschauer R, Nüchter A. Inverse Kinematics of an Anthropomorphic 6R Robot Manipulator Based on a Simple Geometric Approach for Embedded Systems. Robotics. 2023; 12(4):101. https://doi.org/10.3390/robotics12040101

Chicago/Turabian Style

Anschober, Michael, Raimund Edlinger, Roman Froschauer, and Andreas Nüchter. 2023. "Inverse Kinematics of an Anthropomorphic 6R Robot Manipulator Based on a Simple Geometric Approach for Embedded Systems" Robotics 12, no. 4: 101. https://doi.org/10.3390/robotics12040101

APA Style

Anschober, M., Edlinger, R., Froschauer, R., & Nüchter, A. (2023). Inverse Kinematics of an Anthropomorphic 6R Robot Manipulator Based on a Simple Geometric Approach for Embedded Systems. Robotics, 12(4), 101. https://doi.org/10.3390/robotics12040101

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop