A method for the automatic generation of inverse kinematic maps in modular robotic systems

Abstract

Flexible manufacturing based on rapidly reconfigurable robotic systems will enable factories to meet time-sensitive and fast-changing industrial demands. However, with the rise of modular systems there is also the need to quickly and easily determine which configuration is optimal for performing a certain task. In this article, we present a path-based ad hoc technique for determining the inverse and forward kinematics map based on relative joint space variable to reduce the computational complexity. The proposed technique is nonsingular and suits kinematic analysis and optimization of robotic systems with undetermined configuration, and it can be extended to solve generalized inverse kinematic of robotics system involving large number of joint variable.

Introduction

In the ever-changing world of technology, electronic products are continually becoming smaller and smaller, implying the need for the micromanufacturing of system components. The microfactory is an industrial robotic setup that allows for the quick changing of a factory’s layout to meet the changing needs of industry. ^{1 –3} It is popular to use modular components for the construction of robots, as modular components tend to be independent agents, making them amenable to relocation within the robotic infrastructure. ⁴ While downsizing for the purpose of micromanufacturing and modularizing joints for quick reconfigurability, designers must keep in mind the main concerns of the micromanufacturing industry, that is, the research of minimal actuations factory parts and fabrication accuracy.

The fabrication process in an automated manufacturing setting is normally a set of commands given by a central processing unit to a robotic manipulator to follow a certain trajectory. Therefore, the path that a manipulator must take in Cartesian space is often known before a specific robotic manipulator is designed. At this stage, the problem of fabrication comes down to selecting a configuration of prismatic and rotary robotic links that will minimize path error and also minimize the amount of space and time needed to complete the task. The most common way to obtain the joint angles needed for a robotic system is to solve kinematic problem. The approach can be applied when both the path and the robotic manipulator are known.

The determination of configuration using inverse kinematic is a method widely used in biomechanics for the analysis of human motion. For example, inverse kinematic analysis may be used to reconstruct anatomical parameters during some motor tasks such as handwriting or upper limb motor strategies. Vimercati et al. ⁵ observed that participants with Down syndrome tend to draw faster but with less accuracy, or that the motor sequence adopted in a pointing and tapping task may be also affected by the pathology. ⁶ Sale et al. ⁷ adopted an end effector-based robot to control the kinematics of the lower limb in patients with Parkinson’s disease in order to improve their mobility. Another example is the work of Ancillao et al. ⁸ who developed a protocol to reconstruct the kinematics of the upper limb and upper body during a drawing/handwriting task.

At the qualitative level, previous work in determining configurations for robotic and kinematic tasks has taken an agent-based approach, where each link is an agent and the objective is to coordinate the behaviors of all agents. For example, in 1996, Carnegie Mellon University’s Robotics Institute began a program called the Agile Assembly Architecture, whose purpose was to create an Internet interface through which users could input which robotic components they desired to use; a decentralized cooperative framework enabled each component to be modular, and thus allowed the system to be easily reconfigurable, while maintaining high operability. ^{9 –12} In similar work by Chen, ¹³ the tools for developing a reconfigurable robotic workcell were considered (with a workcell being a specific system component or robot). A model for automatic robot generation was presented, which was modular and geometry independent. Optimization of the robotic workcell was specified, and a communication and control framework presented for coordination among the system workcells, with the end goal being a CAD-like program that would enable the appropriate robotic configuration to be determined for a specific task. ¹³

At the quantitative level, Barraquand et al. ¹⁴ presented a robot path-planning method to handle high degrees of freedom (DOFs). In their method, a graph of the local minima of a potential function was constructed incrementally during the search for the goal configuration. Their method relied on numerical methods for a solution. Optimization from a slightly different approach was carried out by Snyman and Berner ¹⁵ who optimized over robotic parameters (link lengths, masses, etc.) rather than positions.

Other work includes the development of inverse kinematic mappings for “hyper-redundant” robots. In the study by Chirikjian, ¹⁶ the author presents a “backbone curve” technique to compare manipulator configurations. Goldenberg et al. ¹⁷ present a technique based on the Newton–Raphson numerical method to iteratively find an inverse kinematic mapping of a kinematically redundant robot. On the other hand, one work that closely addresses accuracy concerns is that by Park and Bobrow, ¹⁸ in which the authors develop methods for optimizing robotic configurations via the “rigid body guidance problem.” Fault tolerance, that is, trajectory planning that takes into account the failure of a joint, was proposed by Paredis and Khosla. ¹⁹ The authors discussed that added redundancy in the form of multiple actuators and joints can aid the robot in completing a task. Their most important point was that certain joint configurations make it difficult to use the manipulator.

In recent years, neural network-based technique has found various applications for inverse kinematics. ²⁰ The work by Li et al. ²¹ proposes a strategy with hierarchical tree generation on topology and nonlinear neural dynamics to control multiple manipulator with only partial command. The method is globally stable and the optimality is guaranteed, which makes it scalable to large number and types of manipulators including redundant DOF system. To solve redundant manipulator under physical constraints, Zhang et al. ²² proposed an online approach based on quadratic programming problem formulation for manipulator torque optimization using various levels of acceleration and velocity schemes.

Another approach to solving the inverse kinematics problem for redundant robots was given by Tabandeh et al., ²³ who used a genetic algorithm to solve the inverse kinematics problem for a redundant robot by way of a niching method. They also used a numerical inverse kinematics solver to attain arbitrary accuracies of joint angles. Work by Chen² and Chen and Yang, ²⁴ addresses and solves the inverse kinematics problem for redundant serial and tree-type modular robots. Assuming that each link has only 1-DOF, is symmetric to all other links, the Jacobian is used to solve differential kinematics equations for the joint angles using the Newton–Raphson method. The work by Pac et al. ²⁵ uses interval analysis in inverse kinematics problem to provide a range of joint space variables for a desired end effector position. Wu et al. ²⁶ describe a rooted tree in module-level approach to generate kinematics of modular robots based on task-oriented concepts.

Meredith and Maddock ²⁷ describe four different types of inverse kinematic mappings, one of which is called the analytical method; analytical methods are described as those that “reduce the inverse kinematics problem to a mathematical equation that need only be evaluated in a single step to produce a result.” ²⁷ One of the drawbacks of analytical methods is that they are often not scalable enough to handle large chains of integrators. ²⁷ In many robotics applications, particularly in the case of a modular robot, it is economically efficient to determine the robotic configuration that will be the best for a certain task from an existing pool of components. For a system with many links, determining an inverse kinematic mapping for each possible configuration is time consuming, and limiting to the time-sensitive efficiency of many factories. Thus, an inverse kinematic mapping that can be applied to many different configurations of robotic components is needed. Given a path, a global definition of the translation and rotation necessary to move from one point to another must be reiterated in terms of the joint variables and joint limits.

Many of the concerns for modular robotic systems have been addressed in the literature, particularly in developing globalized inverse kinematic maps whose coordinates can be mapped to various robotic structures. In addition, many continuous inverse kinematic problems based on linearization method have been addressed in literature to compute a configuration of a manipulator capable of producing a prescribed trajectory in the task space. ^{28 –30} However, much of this work relies on the Homotopy technique or Newton-Raphson method, whose success depends on how close the initial guess of the solution is to the actual solution. ³¹ For initial guesses that are close to the actual solution, the run time of the algorithm is low, and Newton–Raphson techniques are ideal. For large joint displacements over small periods of time, in which smaller sampling intervals must be considered to maintain reasonable initial guesses for Newton–Raphson, computational complexity will increase.

Research aim

The first goal of this work is to present an algorithm that estimates the amount of rotation and translation needed to traverse a planar path using the path’s geometry. The second goal of this work is to outline a mapping scheme from a simplified robot to an actual 6-DOF manipulator, and then compare the manipulator’s performance for different configurations.

The contribution of this article is to present a framework in which the joint angles may be retrieved for motion along planar paths in Cartesian coordinates. This method is different from previously derived inverse kinematics methods in that no knowledge of the robot’s structure is needed to determine the joint angles. The method used is a non-gradient-based optimization problem that focuses on changes in the joint angles rather than the actual joint angles of the robots at each point during the task. Since these coordinates may be mapped to the workspace of any robotic configuration, this framework is ideal for reconfigurable robots. The contribution is also to compare the performance of a more complicated robotic arm to a simplified system, enabling the fast determination of manipulator efficiency along a path.

This article is organized as follows: we present the materials and methods in the “Materials and methods” section. We introduce inverse kinematic mapping to find joint angles along a path. Then, we introduce a simple robot with four links, which is used to implement an existing inverse kinematics algorithm. Then, we map the joint angles from the inverse kinematics map to a physical robot by introducing constraint equations. We present the results and discussion in the “Results and discussion” section. In this section, we discuss simulation results highlighting the use of the inverse kinematic map and the results of mapping joint angles to a physical robot. Finally, in the “Conclusion” section, we offer concluding remarks and discuss future directions for research.

Materials and methods

Path-based determination of the inverse kinematic map

The inverse kinematic mapping of a robot—serial, parallel, and mixed—is based on a method that involves approximating the curve by sampling points along the curve, and using the joint center of gravity and the vectors between two consecutive path points to obtain translation and rotation of the joint. For consecutive movement increments, the amount of rotation and translation is divided among the joints in proportion with the total amount of movement the joint is capable of. Suppose we have a path $P\in {ℝ}^3$ parameterized by time, so that $P=P\left(t\right)$ . At each time step t, the path is specified by a point $p\left(t\right)\in {ℝ}^3$ , where $t=1,2,\mathrm{...}\dots ,n$ and n is the total number of time-parameterized data points that comprise path $P$ as illustrated in Figure 1. Let two consecutive points be given by p and q, where $p=p\left(t\right)$ and $q=p\left(t+1\right)$ . The homogeneous transformation of the end effector from point p to point q is given by

$T\overline{p}=\overline{q}$

1

where $\overline{p}={[p, 1]}^T$ , $\overline{q}={[q, 1]}^T$ , and where $T\in SE\left(3\right)$ is the homogeneous transformation

$T=\left[\begin{array}{ll}R\hfill & \overline{r}\hfill \\ 0_{1\times 3}\hfill & 1\hfill \end{array}\right]$

2

$R=R\left(\Theta \right)\in SO\left(3\right)$ is a rotational transformation and $\overline{r}=\left[r_x,r_y,r_z\right]$ is the translation of the center of mass.

Figure 1.

Projections of the center of mass on a hemisphere, with $\Delta t=0.1$ The black squares represent the paths points p and q , and the lines between points p and q represent the line segments $l_i$ . The dotted lines are the end effector vectors $e_i$ , and the black circles are the center of mass.

Notice that to solve this equation, two quantities must be specified: the translational joint angles, given by $\overline{r}$ , and the rotational joint angles, given by Θ. With only knowledge of the path, finding this information reduces to a geometric problem, detailed below.

Translation of the end effector

Suppose that between every consecutive pair of points p and q there is a straight line $l_t$ connecting point p to point q. This line segment has a length $\left|l_t\right|$ , and an orientation in the base frame $O=O\left({\theta }_t,{\rho }_t,{\phi }_t\right)$ , where θ is the angle of rotation about x, ρ is the angle of rotation about y, and φ is the rotation about z. We assume that the end effector of the manipulator is defined in coordinate frame with a vector in ${ℝ}^3$ denoted by $e_t$ , is always orthogonal to the line segment $l_t$ and axis of rotation for the path $\widehat{u}$ , that is, it lies in the plane of rotation. Therefore, the end effector $e_t$ originates at the center of mass of the final link and ends at the path point q. To find the line segment $e_t$ , we assume that the axis of rotation for the path, is known; in fact, it is the axis orthogonal to the plane of the path. With knowledge of the path’s axis of rotation (assumed to be parallel to one of the global axes $\widehat{x}$ , $\widehat{y}$ , or $\widehat{z}$ ), we can form a 3-axis body coordinate system in which the axis of rotation is the first axis, and $l_t$ is the second axis. We would like for the final axis to point inward toward the center of rotation, thus serving as the end effector $e_t$ and completing the 3-axis coordinate system. Therefore, the body frame’s coordinate system is given by $b=b\left(\widehat{u},l_t,e_t\right)$ . We have that

$e_t=\widehat{u}\times l_t$

3

If we let d to be the distance from the path to the center of mass of the final link, then $e_t$ is given by

$e_t=d\frac{e_t}{\left|e_t\right|}$

4

We assume that $e_t$ acts through point q; thus, defining the center of mass to be the variable g, the location of the center of mass in ${ℝ}^3$ at point p is given by $g_t=l_t+e_t$ . The change in location of the center of mass is what gives us the incremental translation of the manipulator. Hence, we have the translation as

$\delta g_t=g_{t+1}-g_t$

5

Theorem 1

Recall that the translation of the mth link’s center of mass is given by $\overline{r}$ , and let us denote the change in the center of mass’s location as $\delta g\left(t\right)=g\left(t+\delta t\right)-g\left(t\right)$ . Then

$\underset{\delta t\to 0}{\text{lim}} \overline{r}-\delta g\left(t\right)=0$

6

Proof. In appendix 1

The results of this proof indicate that using a large number of data points enhances the accuracy of the inverse kinematic mapping. On the other hand, the number of data points is proportional to the amount of time taken for the robot to complete the path. If we let the time needed to compute a single data point be $t_c$ , the desired time limit be t, and the time increment be $\delta t$ , then the proportionality relationship is given by

$t_{path}=\frac{T}{\delta t}{t}_c$

7

Thus, although smaller time increments increase accuracy, they also increase the run time of the robot along the path. The engineer must balance throughput with accuracy by tweaking such parameters.□

Rotation of the end effector

For any given vector in ${ℝ}^3$ , there are 3 rotational DOFs that may be used to characterize its orientation; they are rotation about the x-, y-, and z-axes. The overall pose of the end effector is often a coupling of all 3 rotational DOFs, making the individual determination of all 3 rotational DOFs elusive at best. In this section, we present a procedural and algorithmic approach to determine the rotation of the end effector in the plane.

Since we have approximated the path with line segments between path-points, we may estimate the amount of rotation that the manipulator has undergone using the dot product between two consecutive end effector orientations $e_{t+1}$ and $e_t$ . Letting α represent the rotational angle needed to define the orientation of the end effector, we have

$\alpha =\text{arccos}\left(\frac{〈e_{t+1},e_t〉}{\left|e_{t+1}\right|\left|e_t\right|}\right)$

8

This gives us the amount of rotation in the plane formed by the two vectors $e_{t+1}$ and $e_t$ .

We measure the accuracy of the robot along the path by comparing the “fit” of the path generated from the inverse kinematics mapping—which is defined by a polygon that is passing through a set of vertices—to the original path, that is, we measure the Euclidean distance between points on the generated path and the corresponding points on the original path. Since error is cumulative, we compare performance at different path parcings by examining the average error per time increment. The error at time t is given by

$\epsilon \left(t\right)=\overline{p}-T\left(0\right)\overline{p}\left(0\right)-T\left(t\right)\overline{p}\left(t\right)$

9

where $T\left(0\right)\overline{p}\left(0\right)$ is the initial location of the end effector. Therefore, for a total of n data points, the average error per unit time along the path is

${\epsilon }_{avg}\left(t\right)=\frac{1}{n}{\displaystyle \sum }_{t=1}^n\left|\epsilon \left(t\right)\right|$

10

Global joint angles from a simplified manipulator

In this section, we demonstrate the use of joint angles from an inverse kinematic map using a simple robot. We consider a robotic manipulator that consists of three linear joints: one that traverses the x direction, one that traverses the y direction, and one that traverses the z direction. It also consists of one rotary joint that is assumed to rotate about the local x- or z-axis. Such kind of robots essentially boil down to the inverse kinematics of a spherical dyad. The four links are shown in Figure 2. There are 64 configurations possible from these components, that is

${\displaystyle \sum }_{k=1}^n\frac{n!}{\left(n-k\right)!}$

11

configurations possible from the four links, and we compare them based on average error and computational cost, where n is the original set of joints. For a given path we measure the error in the manipulator’s performance under the inverse kinematic map based on the equation (10).

Figure 2.

The simple robot in the (1234) configuration shown in the top view and side view. Links 1, 2, and 3 are the x-, y-, and z-directional linear joints, respectively. Link 4 is the 2-DOF rotary joint. The white dot on link 4 is representative of an end effector.

We obtain the joint angles from the paths using the well-known Paden–Kahan subproblems produced by Murray et al. ³² Since there are essentially only two revolute joints, we determine joint angles using Paden–Kahan subproblem 2, featured. The joint angle trajectories are determined from the path by applying subproblem 2 repeatedly to consecutive points along the path. Again, let the first point sampled be denoted by $\overline{p}$ , and the second point by $\overline{q}$ . We wish to find the amount of rotation and translation that has occurred between both points during the time increment $\delta t$ . Starting with rotation, we determine the axis of revolution for the first rotation in sequence to be about the z-axis, whose unit vector is given by ${\widehat{a}}_1={[0,0,1]}^T$ ; we denote this angle by ${\theta }_z$ . The second rotation in sequence, which we can denote ${\theta }_x$ , is about the local x-axis ${\widehat{a}}_2={[1,0,0]}^T$ . Subproblem 2 requires finding an intermediate point $\overline{h}$ between $\overline{p}$ and $\overline{q}$ that can be reached by either of the transformations ³²

$\begin{array}{ccc}{T}_2\left({\theta }_x\right)\overline{p}& =& \overline{h}\\ T_1\left(-{\theta }_z\right)\overline{q}& =& \overline{h}\end{array}$

12

We find the vector, $\overline{z}$ , that connects $\overline{h}$ to a predetermined point on the axis of rotation of the first rotation; here, we’ve chosen that point to be $\overline{r}={[0,0,0,1]}^T$ . Thus, $\overline{z}=\overline{h}-\overline{r}$ . We also can define vectors between point $\overline{p}$ and $\overline{r}$ and point $\overline{q}$ and $\overline{r}$ as ${\overline{w}}_0=\overline{p}-\overline{r}$ , and ${\overline{v}}_0=\overline{q}-\overline{r}$ . Letting ${\widehat{a}}_3={\widehat{a}}_1\times {\widehat{a}}_2$ , and then using a linear combination of the triad ${\widehat{a}}_1$ , ${\widehat{a}}_2$ , and ${\widehat{a}}_3$ gives z as

$z=\alpha {\widehat{a}}_1+\beta {\widehat{a}}_2+\gamma {\widehat{a}}_3$

13

The coefficients α, β, and γ are found by

$\alpha =\frac{〈{\widehat{a}}_1,{\widehat{a}}_2〉{\widehat{a}}_2^T{\overline{w}}_0-{\widehat{a}}_1{\overline{v}}_0}{{(〈{\widehat{a}}_1,{\widehat{a}}_2〉)}^2-1}$

14

$\beta =\frac{〈{\widehat{a}}_1,{\widehat{a}}_2〉{\widehat{a}}_1^T{\overline{v}}_0-{\widehat{a}}_2{\overline{w}}_0}{{(〈{\widehat{a}}_1,{\widehat{a}}_2〉)}^2-1}$

15

and

${\gamma }^2=\frac{|{\overline{w}}_0{|}^2-{\alpha }^2-{\beta }^2-2\alpha \beta 〈{\widehat{a}}_1,{\widehat{a}}_2〉}{|{\widehat{a}}_1\times {\widehat{a}}_2{|}^2}$

16

Using $\overline{z}=\overline{h}-\overline{r}$ , we can solve for the point $\overline{h}$ . To solve for the joint angles, we follow the steps outlined in Murray ³² by solving subproblem 1 for ${\theta }_x$ and ${\theta }_z$ . To find ${\theta }_z$ , if we define

$\overline{w}={\overline{w}}_0-{\widehat{a}}_1{\widehat{a}}_1^T{\overline{w}}_0$

17

and

$\overline{v}={\overline{v}}_0-{\widehat{a}}_1{\widehat{a}}_1^T{\overline{v}}_0$

18

then the rotation may be found as

$\delta {\theta }_z={\text{tan}}^{-1}\left(\frac{{\widehat{a}}_1^T\left(\overline{w}\times \overline{v}\right)}{{\overline{w}}^T\overline{v}}\right)$

19

An identical procedure may be used to find ${\theta }_x$ . To find the amount of translation that has occurred between two points, we first use homogeneous transformations using the two rotational angles on the point p, and then subtract the result from the point $\overline{q}$ . Thus,

$R_1=\left[\begin{array}{lll}\text{cos}\left({\theta }_z\right)\hfill & -\text{sin}\left({\theta }_z\right)\hfill & 0\hfill \\ \text{sin}\left({\theta }_z\right)\hfill & \text{cos}\left({\theta }_z\right)\hfill & 0\hfill \\ 0\hfill & 0\hfill & 1\hfill \end{array}\right]$

20

and

$R_2=\left[\begin{array}{lll}1\hfill & 0\hfill & 0\hfill \\ 0\hfill & \text{cos}\left({\theta }_x\right)\hfill & -\text{sin}\left({\theta }_x\right)\hfill \\ 0\hfill & \text{sin}\left({\theta }_x\right)\hfill & \text{cos}\left({\theta }_x\right)\hfill \end{array}\right]$

21

so that

$T_1=\left[\begin{array}{ll}{R}_1\hfill & 0_{3\times 1}\hfill \\ 0_{1\times 3}\hfill & 1\hfill \end{array}\right]\hspace{1em}\mathrm{and}\hspace{1em}{T}_2=\left[\begin{array}{ll}{R}_2\hfill & 0_{3\times 1}\hfill \\ 0_{1\times 3}\hfill & 1\hfill \end{array}\right]$

22

we have that the translation is given by $\overline{d}=\overline{q}-T_1{T}_2\overline{p}$ . An example using this method is given in the “Results” section.

Mapping joint angles to a physical system

The variables that we wish to find at each time step are the changes in the robot’s joint angles. The robot’s movement in 3-space can then be defined as $f_R=f_R\left(\delta \overline{x},\delta \overline{y},\delta \overline{z},\delta \overline{\theta }\right)$ , where $\delta \overline{x}$ , $\delta \overline{y}$ , $\delta \overline{z}$ , and $\delta \overline{\theta }$ are the joint variables associated with each of the manipulator’s links; each variable can be a multidimensional vector depending on the redundancy of links exhibiting such motion. We define the global position variables obtained from the inverse kinematic mapping to be $\delta x_{\text{global}}$ , $\delta y_{\text{global}}$ , and $\delta z_{\text{global}}$ , while the global rotational displacement is given by $\delta {\gamma }_{\text{global}}$ . Constraints between the physical robot and the three components of Cartesian space are given by equating the global position variables to the projections of the robot’s motion onto the three ordinal directions. That is, for the inertial unit vectors $\widehat{x}$ , $\widehat{y}$ , and $\widehat{z}$

$\begin{array}{ccc}\delta x_{\text{global}}& =& f_R\cdot \widehat{x}\\ \delta y_{\text{global}}& =& f_R\cdot \widehat{y}\\ \delta z_{\text{global}}& =& f_R\cdot \widehat{z}\end{array}$

23

We also represent constraints on the revolute motion of the robot as

$\Delta {\gamma }_{\text{global}}={\displaystyle \sum }_{i=1}^m\hspace{0.17em}{\zeta }_i\delta {\theta }_i,$

24

where m is the number of revolute variables affected by the constraint, ${\gamma }_{\text{global}}$ is the constraint condition, $\delta {\theta }_i$ is the i th revolute variable, and ${\zeta }_i$ is a weight.

To determine the optimal solution angles based on minimum path error between time t and $t+1$ , we propose constrained optimization, with the Lagrangian ³³

$L=f_R\left(\delta \overline{x},\delta \overline{y},\delta \overline{z},\delta \overline{\theta }\right)+{\lambda }_1\left(\delta x_{\text{global}}-f_R\cdot \widehat{x}\right)+{\lambda }_2\left(\delta y_{\text{global}}-f_R\cdot \widehat{y}\right)\hspace{0.17em}+{\lambda }_3\left(\delta z_{\text{global}}-f_R\cdot \widehat{z}\right)+{\lambda }_4\left(\delta {\gamma }_{\text{global}}-{\displaystyle \sum }_{i=1}^m{\zeta }_i\delta {\theta }_i\right)$

25

By taking partial derivatives of the Lagrangian with respect to each variable and setting the resulting equations equal to zero, we have

$\begin{array}{c}\frac{\partial L}{\partial \delta \overline{x}}=0,\hspace{1em}\frac{\partial L}{\partial \delta \overline{y}}=0,\hspace{1em}\frac{\partial L}{\partial \delta \overline{z}}=0,\hspace{1em}\frac{\partial L}{\partial \delta \overline{\theta }}=0\\ \frac{\partial L}{\partial {\lambda }_1}=0,\hspace{1em}\frac{\partial L}{\partial {\lambda }_2}=0,\hspace{1em}\frac{\partial L}{\partial {\lambda }_3}=0,\hspace{1em}\frac{\partial L}{\partial {\lambda }_4}=0\end{array}$

26

With this framework, we are able to solve for each of the desired joint angles of the target robot.

Example using a 6-DOF manipulator

Inverse kinematics using the simplified robotic structure gives us global definitions of the amount of rotation about the x- and y-axes, and translation in all three directions. The next step is to map these globally referenced joint variables to the 6-DOF robotic manipulator. Unlike the simplified representation, only 6-DOF robot has rotary links.

To accomplish the mapping, we note that we have two rotational DOFs centered at the same center of rotation; during our comparison process, we sampled different configurations of the simplified robot, and moved the center of rotation to various positions within the robot’s structure. We were able to compare which rotation center produced the most accurate following of the three-dimensional (3-D) path. As can be seen from Figure 3, the manipulator has six joints like the center of rotation for the simplified robot, each successive pair of joints can rotate in about two directions. Thus, in the mapping process, we can choose which two joints we would like to be the center of rotation and assign the two rotational joint angles to that pair. For the remaining joint angles, a combination of the unused rotational DOFs would need to be used to produce the needed translational motion. Let $\delta {\theta }_i$ ( $i=1,\dots \dots 6$ ) represent the i th joint angle of the 6-DOF manipulator, and let the optimal joint angles be given by $\delta {\overline{\theta }}^{*}={\left[\delta {\theta }_1^{*},\delta {\theta }_2^{*},\delta {\theta }_3^{*},\delta {\theta }_4^{*},\delta {\theta }_5^{*},\delta {\theta }_6^{*}\right]}^T$ . Choosing $\delta {\theta }_1^{*}$ to be fixed and the manipulator’s 5th and 6th rotational joints, $\delta {\theta }_5^{*}$ and $\delta {\theta }_6^{*}$ , respectively, to be the center of rotation, we may remove these angles from the vector $\delta \overline{\theta }$ . For the remaining rotational DOFs, we have the following constraint relationships

$\begin{array}{l}\delta x_{\text{global}}=l_2\text{cos}\left(\delta {\theta }_2\right)+l_3\text{cos}\left(\delta {\theta }_3\right)+l_5\text{cos}\left(\delta {\theta }_5\right)\\ \delta y_{\text{global}}=l_2\text{sin}\left(\delta {\theta }_2\right)+l_3\text{sin}\left(\delta {\theta }_3\right)+l_5\text{sin}\left(\delta {\theta }_5\right)\end{array}$

27

Figure 3.

6-DOF robot. Joints that enable rotation about the local z-axis are denoted by the cylindrical prisms, while joints that enable rotation about the local x-axis are denoted by circles. 6-DOF: 6 degrees of freedom.

We now use constrained optimization to determine the remaining joint. The Lagrangian is given by

$L=\delta \overline{\theta }+{\lambda }_1\left(\delta x_{\text{global}}-l_2\text{cos}\left(\delta {\theta }_2\right)-l_3\text{cos}\left(\delta {\theta }_3\right)-l_5\text{cos}\left(\delta {\theta }_5\right)\right)\hspace{0.17em}+{\lambda }_2\left(\delta y_{\text{global}}-l_2\text{sin}\left(\delta {\theta }_2\right)-l_3\text{sin}\left(\delta {\theta }_3\right)-l_5\text{sin}\left(\delta {\theta }_5\right)\right)\hspace{0.17em}+{\lambda }_3\left({\theta }_B-\delta {\theta }_6-\delta {\theta }_4-\delta {\theta }_1\right)+{\lambda }_4\left({\theta }_A-\delta {\theta }_2-\delta {\theta }_3-\delta {\theta }_5\right)$

28

where ${\theta }_A$ and ${\theta }_B$ are the two rotational DOFs retrieved from the all DOF robot, and ${\lambda }_1$ , ${\lambda }_2$ , ${\lambda }_3$ , and ${\lambda }_4$ are Lagrange multipliers. The partial derivatives of the Lagrangian are given by

$\begin{array}{l}\frac{\partial L}{\partial \delta {\theta }_1}=\frac{\partial L}{\partial \delta {\theta }_6}=\frac{\partial L}{\partial \delta {\theta }_4}=1-{\lambda }_3=0\\ \frac{\partial L}{\partial \delta {\theta }_2}=1+{\lambda }_1{l}_2\text{sin}\left(\delta {\theta }_2\right)-{\lambda }_2{l}_2\text{cos}\left(\delta {\theta }_2\right)-{\lambda }_4=0\\ \frac{\partial L}{\partial \delta {\theta }_3}=1+{\lambda }_1{l}_3\text{sin}\left(\delta {\theta }_3\right)-{\lambda }_2{l}_3\text{cos}\left(\delta {\theta }_3\right)-{\lambda }_4=0\\ \frac{\partial L}{\partial \delta {\theta }_5}=1+{\lambda }_1{l}_5\text{sin}\left(\delta {\theta }_5\right)-{\lambda }_2{ł}_5\text{cos}\left(\delta {\theta }_5\right)-{\lambda }_4=0\\ \frac{\partial L}{\partial \delta {\lambda }_1}=\delta x-l_2\text{cos}\left(\delta {\theta }_2\right)-l_3\text{cos}\left(\delta {\theta }_3\right)-l_5\text{cos}\left(\delta {\theta }_5\right)=0\\ \frac{\partial L}{\partial \delta {\lambda }_2}=\delta y-l_2\text{sin}\left(\delta {\theta }_2\right)-l_3\text{sin}\left(\delta {\theta }_3\right)-l_5\text{sin}\left(\delta {\theta }_5\right)=0\\ \frac{\partial L}{\partial \delta {\lambda }_3}={\theta }_B-\delta {\theta }_6-\delta {\theta }_4-\delta {\theta }_1=0\\ \frac{\partial L}{\partial \delta {\lambda }_4}={\theta }_A-\delta {\theta }_2-\delta {\theta }_3-\delta {\theta }_5=0\end{array}$

29

Suppose we choose joints 5 and 6 to be the center of rotation for the manipulator, and since joint 1 is at the base of the manipulator, we fix it so that it has zero rotation at all times. We can then reduce the size of the problem to one in which we have only three unknown joint angles to solve for. Thus, along with the Lagrange multipliers, the optimization variables we want to find are given by $\delta \overline{\theta }={\left[\delta {\theta }_2,\delta {\theta }_3,\delta {\theta }_4,{\lambda }_1,{\lambda }_2,{\lambda }_3,{\lambda }_4\right]}^T$ , and the new reduced system to be optimized is

$\begin{array}{l}{f}_1=1+{\lambda }_1{l}_2\text{sin}\left(\delta {\theta }_2\right)-{\lambda }_2{l}_2\text{cos}\left(\delta {\theta }_2\right)-{\lambda }_4\\ f_2=1+{\lambda }_1{l}_3\text{sin}\left(\delta {\theta }_3\right)-{\lambda }_2{l}_3\text{cos}\left(\delta {\theta }_3\right)-{\lambda }_4\\ f_3=1-{\lambda }_3\\ f_4=\delta x_{\text{global}}-l_2\text{cos}\left(\delta {\theta }_2\right)-l_3\text{cos}\left(\delta {\theta }_3\right)-l_5\text{cos}\left(\delta {\theta }_5\right)\\ f_5=\delta y_{\text{global}}-l_2\text{sin}\left(\delta {\theta }_2\right)-l_3\text{sin}\left(\delta {\theta }_3\right)-l_5\text{sin}\left(\delta {\theta }_5\right)\\ f_6={\theta }_B-\delta {\theta }_6-\delta {\theta }_4-\delta {\theta }_1\\ f_7={\theta }_A-\delta {\theta }_2-\delta {\theta }_3-\delta {\theta }_5\end{array}$

30

The procedure for finding the core set of partial derivatives when the center of rotation is either links 3 and 4 or 1 and 2 is a similar process of removing those equations for the joint variables that are already known.

Equation (29) must be solved at every time step, and they have no closed form solution. Performing a Newton–Raphson technique at every time step could be computationally time consuming when path is discretized into very small step size. Thus, since we are taking small increments of movement and time, we linearize the equation (29) about the manipulator’s current pose, and solve the linearized system. Let “·” denote the dot product, $\overline{\lambda }=\left[{\lambda }_1,{\lambda }_2\right]$ , and $K\left(\delta {\theta }_i\right)=\left[\text{cos}\left(\delta {\theta }_i\right),\text{ sin}\left(\delta {\theta }_i\right)\right]$ . With the Jacobian matrix

$J=\left[\begin{array}{ccccccc}{l}_2\overline{\lambda }\cdot K\left(\delta {\theta }_2\right)& 0& 0& l_2\text{sin}\left(\delta {\theta }_2\right)& -l_2\text{cos}\left(\delta {\theta }_2\right)& 0& -1\\ 0& l_3\overline{\lambda }\cdot K\left(\delta {\theta }_3\right)& 0& l_3\text{sin}\left(\delta {\theta }_3\right)& -l_3\text{cos}\left(\delta {\theta }_3\right)& 0& -1\\ 0& 0& 0& 0& 0& -1& 0\\ l_2\text{sin}\left(\delta {\theta }_2\right)& l_3\text{sin}\left(\delta {\theta }_3\right)& 0& l_5\text{sin}\left(\delta {\theta }_5\right)& 0& 0& 0\\ -l_2\text{cos}\left(\delta {\theta }_2\right)& -l_3\text{cos}\left(\delta {\theta }_3\right)& 0& -l_5\text{cos}\left(\delta {\theta }_5\right)& 0& 0& 0\\ 0& 0& -1& 0& 0& 0& 0\\ -1& -1& 0& 0& 0& 0& 0\end{array}\right]$

31

The linearized system is given by $J\delta {\overline{\theta }}^{*}=0$ . Thus, we survey the null space of equation (30) to find our optimal solution $\delta {\overline{\theta }}^{*}={[\delta {\theta }_2^{*},\delta {\theta }_3^{*},\delta {\theta }_4^{*},{\lambda }_1,{\lambda }_2,{\lambda }_3,{\lambda }_4]}^T$ at each time step. To show that the joint angles obtained from the null space of the linearized system are suitable for equation (31), we have the following theorem.

Theorem 2

As the size of the state $\delta {\overline{\theta }}^{*}\to 0$ , the linear approximation of the partial derivatives given in equation (31) approaches the nonlinear partial derivatives of equation (30).

Proof. In Appendix 2

For sufficiently small δθ, we can use the null space of the Jacobian matrix J to find the joint angles $\delta {\theta }^{*}$ . This raises the question of what is “sufficiently small” for the constructed framework, and how do we measure whether or not the joint angles commanded by the optimization fall within limits for which the linearization is valid? These questions lead to the next theorem.□

Theorem 3

The Jacobian matrix J given by equation (31) is evaluated at the current value of the state δθ. That is

$J=\frac{\partial f_{i,L}}{\partial \delta \overline{\theta }}{|}_{\delta \overline{\theta }=\delta {\overline{\theta }}^{*}}\to J\delta \overline{\theta }=J\left(\delta {\overline{\theta }}^{*}\right)\delta {\overline{\theta }}^{*}$

32

Therefore, the linearization is valid for arbitrarily large $\delta {\theta }^{*}$ .

Proof. In appendix 3

Upon finding the optimal joint angles, we substitute them into the forward kinematic map for the 6-DOF robot. If we let the matrix describing local rotations about the y-axis be given by

$R_y\left(\delta {\theta }_i\right)=\left[\begin{array}{lll}\text{cos}\left(\delta {\theta }_i\right)\hfill & 0\hfill & \text{sin}\left(\delta {\theta }_i\right)\hfill \\ 0\hfill & 1\hfill & 0\hfill \\ -\text{sin}\left(\delta {\theta }_i\right)\hfill & 0\hfill & \text{cos}\left(\delta {\theta }_i\right)\hfill \end{array}\right]$

33

and the matrix describing local rotations about the z-axis be given by

$R_z\left(\delta {\theta }_i\right)=\left[\begin{array}{lll}\text{cos}\left(\delta {\theta }_i\right)\hfill & -\text{sin}\left(\delta {\theta }_i\right)\hfill & 0\hfill \\ \text{sin}\left(\delta {\theta }_i\right)\hfill & \text{cos}\left(\delta {\theta }_i\right)\hfill & 0\hfill \\ 0\hfill & 0\hfill & 1\hfill \end{array}\right]$

34

and furthermore, let the center of mass of the i th link be given by $d_i$ , the respective homogeneous transformations for rotary links about the y-axis and about the z-axis are given by

$T_y\left(\delta {\theta }_i\right)=\left[\begin{array}{ll}{R}_y\hfill & d_i\hfill \\ 0\hfill & 1\hfill \end{array}\right]\hspace{1em}and\hspace{1em}{T}_z\left(\delta {\theta }_i\right)=\left[\begin{array}{ll}{R}_z\hfill & d_i\hfill \\ 0\hfill & 1\hfill \end{array}\right]$

35

Thus, if $\overline{p}={[x,y,z,1]}^T$ is the current location of the robotic manipulator, the transformation from the manipulator’s frame (m) to the robot’s base frame (b) is given by

$T_m^b\left(\theta \right)=T_y\left(\delta {\theta }_1^{*}\right)T_z\left(\delta {\theta }_2^{*}\right)T_z\left(\delta {\theta }_3^{*}\right)T_y\left(\delta {\theta }_4^{*}\right)T_z\left(\delta {\theta }_5^{*}\right)T_y\left(\delta {\theta }_6^{*}\right)\overline{p}$

36

Since the system is linearized, using small time increments would be the most advantageous for increased accuracy along the path.□

Theorem 4

An optimal solution always exists for the system $J\delta {\theta }_i^{*}=0$ .

Proof

Since any system $Ax=0$ has at least the solution $x=0$ , the system (35) always has a solution. Therefore, the inverse kinematic mapping based on the relative joint space method has no singularities.

Transforming error of the simple robot to the 6-DOF manipulator

What if, given the average error per unit time of the simplified robot, we would like to predict the average error of the 6-DOF manipulator? Based on equation (13), we define the error measured by the simple robot to be

${\epsilon }_{\text{simple}}\left(t\right)=\overline{q}\left(t\right)-T_{\text{simple}}\left(0\right)\overline{p}\left(0\right)-T_{\text{simple}}\left(t\right)\overline{p}\left(t\right)$

37

where $p\left(t\right)$ and $q\left(t\right)$ are the beginning and end points of the travel segment, $T_{\text{simple}}\left(t\right)$ is the homogeneous transformation from the end effector’s frame to the base frame of the simple robot. In particular, we have included the initial location of the simple robot’s end effector as $T_{\text{simple}}\left(0\right)\overline{p}\left(0\right)$ and subtracted it from $q\left(t\right)$ to “realign” the path traced by the physical robot with the commanded path given by the equation. Likewise, the error measured by the 6-DOF manipulator is given by

${\epsilon }_{6\text{DOF}}\left(t\right)=\overline{q}\left(t\right)-T_{6\text{DOF}}\left(0\right)\overline{p}\left(0\right)-T_{6\text{DOF}}\left(t\right)\overline{p}\left(t\right)$

38

Let ${\xi }_{\text{simple}}=T_{\text{simple}}\left(0\right)\overline{p}\left(0\right)$ and ${\xi }_{6\text{DOF}}=T_{6\text{DOF}}\left(0\right)\overline{p}\left(0\right)$ . With these substitutions, and dropping the $\left(t\right)$ from the defined variables for simplicity, we may solve equation (37) for the position $\overline{p}$ as

$\overline{p}=T_{\text{simple}}^{-1}\left(\overline{q}-{\xi }_{\text{simple}}-{\epsilon }_{\text{simple}}\right)$

39

Substituting equation (39) into equation (38) for $\overline{p}\left(t\right)$ gives

${\epsilon }_{6\text{DOF}}=\overline{q}-T_{6\text{DOF}}{T}_{\text{simple}}^{-1}\left(\overline{q}-{\xi }_{\text{simple}}-{\epsilon }_{\text{simple}}\right)-{\xi }_{6\text{DOF}}$

40

This transformation gives a mapping at time t from the error measured by the simple robot to the error measured by the 6-DOF manipulator. Since the transformations $T_{\text{simple}}$ and $T_{6\text{DOF}}$ are in the group $SE\left(3\right)$ , they are guaranteed to be invertible for all time t.

Cost function

Different configurations have dissimilar path-following attributes that may make certain configurations more desirable than others. For instance, the more a linear joint moves along its DOFs, the more actuation is needed to produce the motion. Also, the heavier a joint, the more stress on the actuators to overcome the joint’s inertia. Finally, the most important factor—how well the robot, using the inverse kinematic map, follows the path, or the accuracy of the robot along the path—must be considered in a cost function. In this section, we focus on formulating a cost function that encapsulates all of these concepts in order to aid us in identifying which configuration has the most efficient behavior along a given path.

The cost function is composed of the joint variables, the inertia penalty, and the path error. We choose an energy-like cost function so that it will be amenable to control design. To obtain the cost, we use

$C=\frac{1}{2}\left({\overrightarrow{x}}^T{M}_x\overrightarrow{x}+{\overrightarrow{y}}^T{M}_y\overrightarrow{y}+{\overrightarrow{z}}^T{M}_z\overrightarrow{z}+I_p+E_x^2+E_y^2+E_z^2\right)$

41

where $\overrightarrow{x}$ , $\overrightarrow{y}$ , and $\overrightarrow{z}$ are the $n\hspace{0.17em}\times \hspace{0.17em}1$ vectors containing the x–y–z variables for each joint of the robot, $M_x$ , $M_y$ , $M_z$ are positive semi-definite weighting matrices, where we have chosen to set each equal to

$M_x=M_y=M_z=\left[\begin{array}{llll}{m}_1\hfill & 0\hfill & 0\hfill & 0\hfill \\ 0\hfill & m_2\hfill & 0\hfill & 0\hfill \\ 0\hfill & 0\hfill & \ddots \hfill & 0\hfill \\ 0\hfill & 0\hfill & 0\hfill & m_n\hfill \end{array}\right]$

42

and $I_p$ is the inertia penalty. Let $\overrightarrow{\Theta }$ be the $n\hspace{0.17em}\times \hspace{0.17em}1$ vector containing the angular displacement for each of the n joints, let $I_1$ , $I_2$ , …, $I_n$ be the inertia of each joint about its own center of mass, and let

$\begin{array}{ccc}{I}_{r,1}\left(t\right)& =& m_1{\displaystyle \sum }_{k=1}^K{d}_{1,k}{(t)}^2\\ \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}\\ I_{r,2}\left(t\right)& =& m_2{\displaystyle \sum }_{k=1}^K{d}_{2,k}{(t)}^2\\ \hspace{0.17em}& \hspace{0.17em}& ⋮\\ I_{r,n}\left(t\right)& =& m_n{\displaystyle \sum }_{k=1}^K{d}_{n,k}{(t)}^2\end{array}$

43

where $d_{1,k}$ , $d_{2,k}$ ,…, $d_{n,k}$ are the time-varying distances of the joints from each of the K centers of rotation. Note that only joints between the rotary joint and the path will be affected by the rotary inertia. Placing all of these components together, we obtain the inertia penalty as

$I_p=\bar{\Theta }{(t)}^T\left[\begin{array}{llll}{I}_1\hfill & 0\hfill & 0\hfill & 0\hfill \\ 0\hfill & I_2\hfill & 0\hfill & 0\hfill \\ 0\hfill & 0\hfill & \ddots \hfill & 0\hfill \\ 0\hfill & 0\hfill & 0\hfill & I_n\hfill \end{array}\right]\bar{\Theta }\left(t\right)+\bar{\Theta }{(t)}^T\left[\begin{array}{llll}{I}_{r,1}\left(t\right)\hfill & 0\hfill & 0\hfill & 0\hfill \\ 0\hfill & I_{r,2}\left(t\right)\hfill & 0\hfill & 0\hfill \\ 0\hfill & 0\hfill & \ddots \hfill & 0\hfill \\ 0\hfill & 0\hfill & 0\hfill & I_{r,n}\left(t\right)\hfill \end{array}\right]\bar{\Theta }\left(t\right)$

44

$E_x$ , $E_y$ , and $E_z$ are the x-, y-, and z-projections of the path error, given by

$E_x=〈E,\widehat{x}〉$

45

$E_y=〈E,\widehat{y}〉$

46

$E_z=〈E,\widehat{z}〉$

47

where E is the difference between the actual path $p_{\text{actual}}$ and the path generated from the inverse kinematics joint variables $p_{\text{IK}}$ , given by

$E=p_{\text{actual}}-p_{\text{IK}}$

48

Results and discussion

In this section, we present the results of simulations using the proposed methodology for determining two-dimensional (2-D) path inverse kinematic maps, and for the case study of multiple configurations of a simplified robot whose joint angles have been mapped to a 6-DOF robotic arm. We implemented the algorithms in MATLAB (T. M. MATLAB 8.5, Inc., Natick, Massachusetts, United States). ⁶

2-D path generation along a planar sinusoid—error cost function

The planar path that we consider is given by

$\begin{array}{ccc}x\left(t\right)& =& \text{cos}\left(t\right)\\ y\left(t\right)& =& \text{sin}\left(t\right)\\ z\left(t\right)& =& 0\end{array}$

49

where $0\le t\le 2\pi$ in steps of size $\delta t$ . The simplified robot is in the (1234) configuration in Figure (3). Following the analysis in section “Materials and methods,” we assumed that the lengths of the three linear links are $L_x=L_y=L_z=0.1$ m.

The results are shown in Figure 4, and confirm that smaller time increments—which correspond to small δθ increments in Theorem 3—produce higher accuracy between the generated curve and the curve produced using the output of the inverse kinematic map.

Figure 4.

Tracing of a curved planar sinusoidal trajectory using the proposed inverse kinematic mapping. The time increments of $\delta t=0.1$ and $\delta t=0.001$ give average path errors of $e=0.1609$ and $e=0.0016$ , respectively. The error cost function is linearly proportional to the time increment, and both have the same order of magnitude.

2-D path generation along a planar sinusoid—Generalized cost function

To test the inverse kinematic mapping, we use a three stages “joints” robot. Each slot is stacked vertically (along the y-direction). We have two linear joints, the first with one redundant DOF (x) and the other with two DOFs (x–z), and a rotary joint. We make the simplifying assumption that the rotary joint is a disk, and the two linear joints are rectangular plates as shown in Figure 5. We will use this system to test the inverse kinematic mapping along a sinusoid. The joint dimensions, movement capabilities, and slot information are shown in Table 1.

Figure 5.

Physical representation of each stage “joint”.

Table 1.

Part and slot information.

Part no.	Range	Symbol	End effector	Dimensions	Mass
1	x = 0–1.0 m	a1	0.125 m	1:1.395 m; w:0 .225 n1; h:0.100 m	44.1 kg
2	x = 0–0.15 m; z = 0–0.15 m	a2; d2	0.1375 m	1:0.275 m; w:0.130 m; h:0.04 m	4.6 kg
3	θ = 0–2π rad	Θ3	0.0475 m	h:0.055 m; r:0.0475 m	2.0 kg
Slot no.	Rel. Loc. (m)	Symbol	Rel. Rot.	Relative to:
1	[0, 0, 1.950]	[S11, S12, S13]	I3 × 3	Origin
2	[0, 0.15, 0]	[S21, S22, S23]	I3 × 3	1
3	[0, 0.15, 0]	[S31, S32, S33]	I3 × 3	2

For this exercise, we study the robot in the 1-2-3 configuration, whose forward kinematic mapping is given by

$T_e^O(\theta )=\left[\begin{array}{llll}\text{cos}\left({\theta }_3\right)\hfill & 0\hfill & -\text{sin}\left({\theta }_3\right)\hfill & S_{11}+S_{21}+S_{31}+a_1+a_2+\frac{19}{400}\text{cos}\left({\theta }_3\right)+\frac{21}{80}\hfill \\ 0\hfill & 1\hfill & 0\hfill & S_{12}+S_{22}+S_{32}\hfill \\ \text{sin}\left({\theta }_3\right)\hfill & 0\hfill & \text{cos}\left({\theta }_3\right)\hfill & S_{13}+S_{23}+S_{33}+d_2+\frac{19}{400}\text{sin}\left({\theta }_3\right)\hfill \\ 0\hfill & 0\hfill & 0\hfill & 1\hfill \end{array}\right]$

50

We can already see that for the selection of joint types and their slot orientations, the path of the end effector will always be planar, regardless of configuration.

The general equation for the sinusoid that we use is

$z_t=p_{\text{start}}+A\text{sin}\left(\frac{\Omega }{B}{x}_t\right)$

51

where $p_{\text{start}}$ is the insertion point in the x–z plane of the sinusoid function, A is the amplitude of the sinusoidal wave, Ω is the frequency, and B is the wavelength. Using the position in $T_e^O\left(0\right)$ , we were able to identify values of $x_t$ and $z_t$ at $t=0$ (given by $p_{1,x}$ and $p_{1,z}$ , respectively) and solve for $p_{\text{start}}$ . With the parameters $A=0.01$ , $B=0.5$ , and $\Omega =2\pi$ , we determined $p_{\text{start}}$ to be given by the equation

$p_{\text{start}}=p_{1,z}-A\hspace{0.17em}\text{sin}\left(\frac{\Omega }{B}{p}_{1,x}\right)$

52

where $p_{1,z}$ is the z-component of the insertion point, and $p_{1,x}$ is the x-component.

Table 2 and Figure 6 compare the accuracy of the inverse kinematic map along the sinusoidal curve for increasing numbers of data points.

Table 2.

Comparison of accuracy with 10, 20, 40 and 100 data points.

Number of data points	1-Norm of position error
10	0.8609
20	0.3687
40	0.1268
100	0.0927

Figure 6.

(a) 40 data points and (b) 100 data points.

We use the 1-norm for comparison because using a single number allows us to quickly see the net effect of error over time. Taking the 1-norm of the difference between the desired and actual paths confirms that increasing the number of data points increases accuracy of the estimated path. Using the automatic-inverse kinematic map generation algorithm, we were able to compare n! configurations for the three-joint robot. The results are shown in Table 3.

Table 3.

Comparison of 1-norm of costs for all configurations when n = 3.

Permutation	1-Norm of cost
(3-1-2),(3-2-1),(1-3-2),(1-2-3),(2-1-3),(2-3-1)	546.3976
(1-3), (3-1)	550.9591
(2-3), (3-2)	12.3874
(2-1), (1-2)	556.1739
(3)	16.4216
(1)	537.3598
(2)	12.4337

Under the chosen cost function, all of the six configurations that use all three joints have similar costs, while marked differences in cost can be seen for the six configurations that use only two of the joints and the three configurations that use only one of the joints. The reason for this pattern in the six configurations that use all three joints may be that regardless of the joint sequence, the same net movement is undertaken by the robot. This seems to indicate that the order in which parts are assembled does not matter, and could be due to the commutative nature of many of the joints involved and the fact that the approximate inverses always have a “trivial” unique solution. Increases in the inertia penalty due to configurations that place the rotary joint in either the first slot or the second slot are leveraged by smaller translations made by the joints. The “sweeping motion” of the rotary joint, particularly when the rotary joint is further in position from the path, eliminates much of the needed translation.

The costs for two-joint configurations is largely dependent on the joints being used. For instance, configuration 1-3 has a much lower cost than configuration 2-3 because 2-3 cannot traverse the entire sinusoid, and therefore accrues much more error than configuration 1-3. Configurations 1-2 and 2-1 actually have fairly low costs, even without the help of the rotary joint, because they are independently able to traverse the entire sinusoidal curve.

3-D path generation along a curved sinusoid

To test the efficiency of the various configurations of the simple robot, we consider using a curved sinusoid, given by the parametric representation

$\begin{array}{l}x\left(t\right)=\text{cos}\left(t\right)\\ y\left(t\right)=\text{sin}\left(t\right)\\ z\left(t\right)=\text{sin}\left(15t\right)\end{array}$

53

where $\pi /3\le t\le 2\pi /3$ , with increments of $\delta t=0.001$ . We tested both the simple robot and the 6-DOF manipulator along the path. Starting with the simple robot of four groups of configurations shown in Figure 7, we chose to make it both massless and dimensionless, that is, all motion is concentrated at a single moving point. All 64 configurations were tested. Configuration sets that have the same cardinality (4, 3, 2, or 1 link(s)) also have the same average error, because of the small increments in movement taken. When infinitesimally small changes in the joint angles are used in majority of simplified manipulator, the end effector’s motion is independent of the sequence of links. The findings in Table 4 indicate that differences in average error between configurations depends on the number of links used. The same idea holds true for computational efficiency. The number of computational operations depends primarily on the implementation of the algorithm, and this number is independent of permutation when configurations of the same cardinality are used. Between cardinalities, the difference lies in the length of the chain of exponentials used to compute the forward kinematic map. For the multiplication of two matrices A and B of size $b\hspace{0.17em}\times \hspace{0.17em}n$ and $n\hspace{0.17em}\times \hspace{0.17em}p$ , matrix multiplication alone has $bp\hspace{0.17em}\left(2n-1\right)$ operations (often referred to as “flops”).

Figure 7.

Schematic of simple robot with four groups of configurations.(a) Group of one rotational joint, (b) group of two rotational joints, (c) group of three rotational joints, and (d) group of four rotational joints.

Table 4.

Four groups of configurations for the simple robot based on the number of links. The time increment is $\delta t=\mathrm{0.001.}$

Assigned number	Configurations (in order)	Average error (mm), computational operations
1,2,3,4,5,6,7,8,9,10,11,12,13; 14,15,16,17,18,19,20,21, 22,23 ;24	(1432),(1423),(1342),(1324), (1234),(1243), (4132),(4123), (4312),(4321),(4231),(4213), (3412),(3421),(3142),(3124), (3214),(3241), (2431),(2413), (2341),(2314),(2134),(2143)	0.000014344, 471
25,26,27,28,29,30,31,32, 33,34,35,36,35,37,38, 39,40,41,42,43,44,45,46,47,48	(341),(314),(431),(413),(143), (134),(241),(214),(421),(412), (142),(124),(231),(213),(321), (312),(132),(123),(234),(243), (324),(342),(432),(423)	0.002154348, 407
49,50,51,52,53,54,55,56,57,58,59,60	(41),(14),(31),(13),(21),(12), (34),(43),(24),(42),(23),(32)	0.009943061, 343
61,62,63,64	(1),(4),(3),(2)	0.009929444, 279

Depending on how the forward kinematics algorithms are implemented, matrix computations can be the most costly. Thus, the programmer must weigh the cost of production, judging between accuracy and computational efficiency.

The patterns in accuracy and computational efficiency of the simple robot can be extended to the 6-DOF manipulator when considered as separate cases. However, generating a correspondence relationship between the simple robot’s four links and the links of the 6-DOF manipulator is elusive at best, since several of the 6-DOF robot’s joint angles contribute to multiple global DOFs. To see an example of this, suppose that we set $\delta x_{\text{global}}$ to zero to represent the removal of the simple robot’s x-directional link. The unevaluated Jacobian matrix equation (31) will still be the same as if we had not. The only way to alter the Jacobian would be to set $\delta x_{\text{global}}$ and $\delta y_{\text{global}}$ to zero (removing the x- and y-directional links of the simple robot). Doing this would leave the equation ${\theta }_B-\delta {\theta }_6-\delta {\theta }_4-\delta {\theta }_1=0$ , and there are an infinite number of solutions to this equation.

Mapping error from the simple robot to the 6-DOF manipulator is more reliable. Using equation (40) and only considering the case when all four links of the simple robot and all six links of the 6-DOF robot are used, the results of the error mapping are shown in Figure 7. The inputs were the forward kinematic maps of the simple and 6-DOF robots, the error of the simple robot at time t, the end point $\overline{q}$ of a path segment at time t, and the initial positions ${\xi }_{\text{simple}}$ and ${\xi }_{6\text{DOF}}$ of each robot’s end effector. To generate Figure 8, the result of equation (40) was subtracted from the original path’s coordinates to obtain the actual path followed by the 6-DOF robot. As expected, since we use all four links of the simple robot, the error is relatively small, where there is no difference observed between the plotted curves in the graph.

Figure 8.

Comparison of path following given the original path $p\left(t\right)$ , the path generated by the 6-DOF manipulator’s inverse kinematics map, the path generated by $p\left(t\right)-{\epsilon }_{6\text{DOF}}\left(t\right)$ , and the path generated by $p\left(t\right)-{\epsilon }_{\text{simple}}\left(t\right)$ . 6-DOF: 6 degrees of freedom.

Conclusion

In this work, we have described a computationally efficient joint space algorithm that determines the configuration of a modular robotic system for a given path. We presented a framework in which the joint angles may be retrieved for motion along planar paths in Cartesian coordinates, where we reconstruct robot configurations using the change in the joint angles rather than the absolute values. To accomplish this, we first determine the inverse kinematic mapping based solely on the path traveled; then, applying the inverse kinematic map to all configurations of the robot, we compare each configuration based on accuracy in following the path, and minimal computation required. The developed algorithm allows users to input information about the joints of a robotic platform and receive an automatically generated forward kinematic mapping for any configuration of the joints, as well as the time-dependent joint variables along a given path. The generalized forward kinematic mapping can handle a set of distinct robotic joints as well as duplicate joints and still produce the appropriate mapping by the utilization of numerical part markers that the algorithm uses to query information about the joints considered. Using this inverse kinematic mapping has proven to render accurate reproductions of a path with bounded curvature and without singularity when the measured joint variables were placed through the forward kinematic mapping.

We also successfully map joint angles obtained from an existing 3-D inverse kinematics method to a physical 6DOF manipulator and set up the framework for comparing multiple configurations of this robot. Although such simple manipulator may have closed form solution, but our method is suitable for solving inverse kinematic of robotics system involving large number of joint variable, example may include self-assembly of molecular scale robots. This new comparison metric will enable designers and micro-machinists to build manufacturing platforms with prior knowledge of the configurations that produce the optimal output for the task, hence speeding up the rapid prototyping process significantly.

In future work, we will extend the two dimensional inverse kinematics scheme to three dimensions in order to completely specify a mechanism for the fast determination of optimal configurations in modular robotic platforms. We will also develop a cost function to more fully characterize the trade-off between computational efficiency and accuracy.