The Natural-CCD Algorithm,a Novel Method to Solve the Inverse Kinematics of Hyper-redundant and Soft Robots

Abstract

This article presents a new method to solve the inverse kinematics problem of hyper-redundant and soft manipulators. From an engineering perspective, this kind of robots are underdetermined systems. Therefore, they exhibit an infinite number of solutions for the inverse kinematics problem, and to choose the best one can be a great challenge. A new algorithm based on the cyclic coordinate descent (CCD) and named as natural-CCD is proposed to solve this issue. It takes its name as a result of generating very harmonious robot movements and trajectories that also appear in nature, such as the golden spiral. In addition, it has been applied to perform continuous trajectories, to develop whole-body movements, to analyze motion planning in complex environments, and to study fault tolerance, even for both prismatic and rotational joints. The proposed algorithm is very simple, precise, and computationally efficient. It works for robots either in two or three spatial dimensions and handles a large amount of degrees-of-freedom. Because of this, it is aimed to break down barriers between discrete hyper-redundant and continuum soft robots.

Introduction

The concept of degree-of-freedom (DOF) is very important in robotics. It refers to the number of movements or independent parameters that define the configuration of a robot.¹ At least 6 DOFs are required to properly position and orient a robot in a three-dimensional (3D) space: 3 to determine the position and 3 to determine the orientation. Hyper-redundant manipulators or snake-arm robots rely on a number of DOFs higher than the minimum required to perform a particular task. They have been widely studied due to their potential benefits: they are able to perform better obstacle avoidance than the conventional ones, they are fault-tolerant,² and also offer mechanic and kinematic advantages, such as maximization of the transmissibility of force, dexterity, rigidity, energy efficiency, or range availability of the joints. In other words, they are characterized for having a high versatility and can perform adequately even after a mechanical failure of some joints. They can be used for many applications such as surgical environments,^3,4 servicing underground tanks,⁵ nuclear core reactors,⁶ or aeronautical and space exploration.⁷ Although the study of hyper-redundant robots has more than three decades of history, their use is mostly limited to the research field. This is because of three main reasons. In first place, its mechanical design and implementation are often perceived as unnecessarily complex. In second place, the classical kinematic modeling techniques have not been particularly efficient or well adapted to the needs of hyper-redundant robots.⁸ And in third place, a large number of DOFs result in a high cost of design and production.

The main purpose of this research is to solve the inverse kinematics problem of hyper-redundant robots in an efficient manner. This problem, from a robotic point of view, aims at determining the joint angles, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\bar q \in { \mathbb{R}^n}$$ \end{document} , which yield a desired end-effector position and orientation, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${ { \bar x} _{ee}} \in { \mathbb{R}^m}$$ \end{document} . If the manipulator is redundant, it will satisfy the Equation (1).⁹ \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} n = \dim \left( { \mathop { \bar q} \limits^{} } \right) > \dim \left( {{{ \mathop { \bar x} \limits^{} }_{ee}}} \right) = m \tag{1} \end{align*} \end{document}

Because of their redundancy, they can be studied as underdetermined systems. Therefore, they exhibit an infinite number of joint configurations that satisfy the inverse kinematics problem. It is possible to find several approaches in the specific literature that try to choose the most appropriate one. In recent years, exhaustive methods have been replaced by optimization algorithms because the computational cost of the first ones grows faster versus the number of joints. Thus, the most recent optimization methods studied are Pseudo-Inverse Jacobian,¹⁰ Pattern Search,¹¹ Global Search,¹¹ Genetic Algorithms,¹¹ Simulated Annealing,¹¹ Artificial Neural Networks,¹¹ and Particle-Swarm Optimization.¹² The common limitation of these methods is that none of them provides a commitment to the desired criteria for solving the inverse kinematics problem: low computational times for a large number of DOF, high precision, and a good-quality solution.

According to their morphological structure and actuation, hyper-redundant robots can be either discrete or continuum.^13–16 There is a subfield of hyper-redundant continuum robots, developed with soft and deformable materials or compliant mechanical parts, denominated soft robots. The advantages of using those materials include a considerable reduction in the harm that could be caused by traditional robotic systems, increasing their potential to interact with humans. Compliant materials also adapt easier than others to various objects, simplifying tasks such as grasping.¹⁷ However, the movement and redundancy of soft robots often rely on their inherent structure. Therefore, they are difficult to model and hard to control.

The development of computationally efficient and accurate models for soft robots needs further study owing to both material and geometric nonlinearities. Some of the studies on this matter focus on biomechanical models such as squid tentacles,¹⁸ octopus arms,¹⁹ or elephant trunks.²⁰ To solve the inverse kinematics problem for soft robots, existing literature propose machine-learning algorithms,²¹ real-time finite element methods,²² visual servo control,²³ new optimization approaches,²⁴ or closed-form methods.²⁵ Also, the kinematics for multisection continuum robots in real time, for various types of actuation such as tendons or pneumatic pressures, has been proposed.²⁶ The main limitation of existing approaches for soft robots is that currently neither the whole body nor the pose of the end effector (which may be important for manipulation, sensing, etc.) is contemplated in the solution.¹⁷

The main contribution of this work is to break down barriers between discrete and continuum hyper-redundant robots and to overcome the current limitations to solve the inverse kinematics problem. According to some authors,¹³ models that approximate continuum manipulators by finite-DOF hyper-redundant ones can be appropriate. Thus, the natural-cyclic coordinate descent (CCD) algorithm is proposed below to solve the inverse kinematics problem of discrete hyper-redundant and soft robots. In addition, this article originally presents its application to continuous trajectories and robots with both rotational and prismatic joints. It also studies the fault tolerance in hyper-redundant robots and their performance in complex environments with obstacles.

In contrast to the previously mentioned state-of-the-art methods, the natural-CCD algorithm combines whole-body movements, higher precisions, lower computational times for real-time applications, and a management of a much larger number of DOFs. However, more importantly, it proposes one solution that can be extrapolated to soft robotic manipulators. After all, as discussed later, they are both inspired by nature.

Materials and Methods

In this section, the cyclic coordinate descent (CCD) will be explained because it is the basis of the following method. Then, the natural-CCD algorithm will be presented as a modification of its predecessor, the CCD. Also, it will be demonstrated why this algorithm is, to some extent, biomimetic or based on nature. Afterward, an explanation of how this work is applied to soft robotics is provided.

Cyclic coordinate descent

This article presents an original algorithm based on the CCD method.²⁷ The CCD is a derivative-free optimization algorithm that cyclically finds a local minimum of a function, f, using line search along one coordinate direction, q, at the current point in each iteration, k (Eq. 2). Some authors have shown that, under the assumption of minimizing the coordinates one by one, the method converges to minimize the function.²⁸ \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} q_i^{ \left( {k + 1} \right) } \in \arg \mathop { \min } \limits_{y \in \mathbb{R}} f \left( {q_1^{k + 1} , \ldots , q_{i - 1}^{k + 1} , y , q_{i + 1}^k , \ldots , q_N^k} \right) \tag{2} \end{align*} \end{document}

For the inverse kinematics problem, the purpose of the CCD is to minimize the function that relates the current position of the robot's end effector and the desired position to reach under a certain precision, by moving the robot joints in succession. To do this, the algorithm has to find in each iteration which movement of the current joint minimizes this function. It is applicable to robots either in two or three spatial dimensions. For example, only for rotational joints: being p_c the position of the current joint to be moved, p_e the current end-effector position, and p_f the desired end-effector position, there exists an angle \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\theta$$ \end{document} (Eq. 3) and a direction \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\vec d$$ \end{document} (Eq. 4) for which it is possible to rotate the joint. Vectors \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\overrightarrow {{p_c} \;{p_e}}$$ \end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\overrightarrow {{p_c} \;{p_f}}$$ \end{document} will be aligned, and the \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\left\vert { \overrightarrow {{p_e} \;{p_f}} } \right\vert$$ \end{document} distance will be minimized,²⁹ as shown in Figure 1. \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} & \theta = { \cos ^ { - 1 } } \left( { { \frac { { p_e } - { p_c } } { \| { { p_e } - { p_c } } \| } } \, \cdot \, { \frac { { p_f } - { p_c } } { \| { { p_f } - { p_c } } \| } } } \right) \\ & \\ \tag { 3 } \end{align*} \end{document} \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} \vec d = \; { \frac { { p_e } - { p_c } } { \| { { p_e } - { p_c } } \| } } { \rm { x } } { \frac { { p_f } - { p_c } } { \| { { p_f } - { p_c } } \| } } \tag { 4 } \end{align*} \end{document}

FIG. 1.

The CCD algorithm applied for solving the inverse kinematics problem of a planar robot with rotational joints. (a) Movement of a given joint. (b) Movement of the consecutive joint. CCD, cyclic coordinate descent; p_c, position of the current joint; p_e, current end-effector position; p_f, desired end-effector position; θ, rotation angle; d, direction vector.

During the optimization process, the position and orientation can be weighted in different manners, so as to reach the goal requirement that includes both parameters. Furthermore, the algorithm works for either rotational or prismatic joints. This has been studied in the mentioned work²⁷ and will be analyzed afterward.

The natural-CCD and selective coordinate descent algorithms

The CCD algorithm presents three main limitations. First, kinematic singularities are not managed. Second, it does not propose any kinematic constraint to avoid self-collisions. Third, since the \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\theta$$ \end{document} angle can take high values (Eq. 3), the movements of each joint can be abrupt, originating convoluted robot configurations. The solutions to these problems lead to the proposed algorithm: the natural-CCD. These modifications are explained below.

In first place, when vectors \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\overrightarrow {{p_c} \;{p_f}}$$ \end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\overrightarrow {{p_c} \;{p_e}}$$ \end{document} are parallel, the mentioned \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\vec d$$ \end{document} direction results in a null vector (Eq. 4). In this case, it is possible to decide not to move the corresponding joint, because the movement of the previous ones may help avoiding that singularity. However, there are some situations in which all the joints are blocked by singularities; for example, when all of them are collinear with the desired position. In those special cases, the use of a random direction \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\vec d$$ \end{document} (Eq. 4) and angle \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\theta$$ \end{document} (Eq. 3) is recommended to avoid the initial singularity.

To solve the second limitation, the existence of self-collisions, a maximum-angle upper bound is proposed (Eq. 5). This angle, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${ \theta _{max}}$$ \end{document} , is selected as the supplementary of a regular polygon's interior angle, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\gamma$$ \end{document} , with as many sides as links having the robot, N. Although the algorithm has been generalized in a 3D environment, that angle is the extreme case to ensure that any hyper-redundant robot will not self-collide. As a result, in an extreme case, the robot could enclose itself in an N-sided polygon. This angle can be slightly undervalued with a security parameter, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\varepsilon$$ \end{document} , for security reasons when the physical dimensions of the robot are taken into account. \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} { \theta _ { max } } \left( { rad } \right) = \varepsilon \; \left( { \pi - \gamma } \right) = \varepsilon \; \left( { \pi - \frac { { \pi \; \left( { N - 2 } \right) } } { N } } \right) = \; \frac { { 2 \pi \varepsilon } } { N } \tag { 5 } \end{align*} \end{document}

The third problem, the abruptness in the joint movements, has been solved using a variable called k-factor to reduce the value of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\theta$$ \end{document} (Eq. 6). Some authors have proposed similar approaches to coordinate descent algorithms.³⁰ In this study, this variable has been taken as the inverse of the Euclidean distance between the end-effector position and the desired position. The reason for this is that the convergence of the initial CCD algorithm was highly dependent on this distance. The k-factor has an upper limit equal to 1 to ensure the convergence of the method. In contrast with the CCD algorithm, the inclusion of this parameter allows obtaining not only the final configuration of the robot but also a set of valid intermediate configurations that lead to the appropriate solution. \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} k = { \frac { 1 \; } { \left\vert { \overrightarrow{ { { p_e } \; { p_f } } \; }} \right\vert } } \ { \rm where } \ 0 > k \ge 1\ \, { \rm and } \ { \theta ^* } = k \; \theta \tag { 6 } \end{align*} \end{document}

Moreover, it is possible to modify the k-factor to obtain a specific number of cycles for a given movement (Eq. 7). In other words, to manage the step of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\theta$$ \end{document} regardless of the number of joints, N, with the following strategy. \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} k = & \;{C \over {N \cdot \left\vert { \overrightarrow {{p_e} \;{p_f}} { \rm{ \;}}} \right\vert }} \ \;{ \rm{where \ C \ is \ a \ positive \ constant , }}\\ \quad &\ 0 > k \ge 1 \ { \rm{and}} \ { \theta ^*} = k \; \theta \tag{7}\ \end{align*} \end{document}

The main structure of the natural-CCD algorithm can be found in Table 1. Some authors have proposed other coordinate descent approaches, such as the adaptive coordinate descent³¹ that uses adaptive encoding to gradually build a transformation of the coordinate system. With this method, the new coordinates are as decorrelated as possible with respect to the objective function. Another technique worth mentioning is the randomized coordinate descent,³⁰ which selects randomly with uniform probability the coordinate to minimize.

Table 1.

The Natural-CCD Algorithm

while (distance between (p_e, p_f) > precision)
for (q_i = 1, q_i ≤ N, q_i = q_i +1)
Calculate (

\documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${ \rm{ \theta }}$$ \end{document}

)
Solve the abruptness in movements,

if vectors (

) are parallel then manage the singularity end if
if (

≥

) then

end if
Rotate q_i using (θ,

)
end for
Update distance (p_e, p_f)
end while

Natural-CCD, natural-cyclic coordinate descent; p_e, current end-effector position; p_f, desired end-effector position; q_i, current joint; N, number of joints; \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\theta$$ \end{document} , rotation angle; \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${ \rm{ \vec d}}$$ \end{document} , direction vector; k, variable to reduce the value of the rotation angle; p_c, position of the current joint; \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${ \theta _c}$$ \end{document} , angle between the current link to be moved and the previous one; \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${ \theta _{max}}$$ \end{document} , maximum angle upper bound.

In this article, a novel alternative called selective coordinate descent (SCD) is also presented. The aim is to use this new approach in underdetermined systems such as hyper-redundant robots (Eq. 1). The main difference between CCD and SCD is that the first one minimizes all the coordinates or variables uniformly, whereas the second one does not. The SCD selects a candidate to be moved among all the elements so as to meet a specific requirement. For example, in the inverse kinematics problem of hyper-redundant robots, the joint to be moved, in terms of minimizing energy, could be chosen in each cycle. This can be also combined with the modifications presented in the natural-CCD algorithm as shown in Table 2. Also, a comparison between both techniques is shown in Figure 2.

FIG. 2.

Movement of a hyper-redundant robot with 10 rotational joints (a) from an upright configuration to a desired position (b) using the natural-CCD algorithm (c) and the SCD algorithm. SCD, selective coordinate descent.

Table 2.

The Selective Coordinate Descent Algorithm

while (distance between (p_e, p_f) > precision)
Select q_i between [1,N] to meet a specific requirement such as energy efficiency
Manage local minima
Calculate (

)
Solve the abruptness in movements,

if vectors (

) are parallel then manage the singularity end if
if (

≥

) then

end if
Rotate q_i using (θ,

)
Update distance (p_e, p_f)
end while

p_e, current end-effector position; p_f, desired end-effector position; q_i, current joint; N, number of joints; \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\theta$$ \end{document} , rotation angle; \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${ \rm{ \vec d}}$$ \end{document} , direction vector; k, variable to reduce the value of the rotation angle; p_c, position of the current joint; \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${ \theta _c}$$ \end{document} , angle between the current link to be moved and the previous one; \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${ \theta _{max}}$$ \end{document} , maximum angle upper bound.

Biomimetics in the natural-CCD algorithm

It can be said that some trajectories resulting from the natural-CCD algorithm are biomimetic or similar to those observed in nature. Biomimicry derives from the Ancient Greek: bios (life) and mimesis (imitation). It is the science that studies nature as a source of inspiration and new innovative technologies to solve human problems, which the nature itself has already solved through system models, processes, and elements that are inspired therein.³²

Figure 3 shows the movement of a two-dimensional (2D) robot with ten rotational joints and ten equal-length links. It starts from an upright configuration, with a desired end-effector position close to the anchored base and with no orientation restrictions. It can be seen that the final configuration is close to a decagon. In this particular case, it is not a result of the Equation (5) but of the Equation (6), as it will be explained.

FIG. 3.

Movement of a planar robot with 10 rotational joints using the natural-CCD algorithm. (a) Initial upright configuration. (b) Intermediate configurations. (c) Final configuration and the end-effector path representation.

The final end-effector path has been interpolated one time per cycle, only when all joints have had the chance to be moved, as other authors have done in previous works.²⁷ In this case, the resulting paths within a 2D space are often quite particular. They tend to shapes commonly present in nature: logarithmic spirals.

A logarithmic spiral, equiangular spiral, or growth spiral is a self-similar spiral curve for which the radius grows exponentially with the angle. It has the following mathematical definition in Cartesian coordinates (Eq. 8) and in Polar coordinates (Eq. 9), where r is the distance from the origin, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\alpha$$ \end{document} is the angle from the x-axis, and a and b are arbitrary constants. \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} \begin{split}x = r \cos { \rm{ \alpha }} = { \rm{ \;}}a{e^{{ \rm{b \alpha }}}} \cos { \rm{ \alpha }} \\y = r \sin { \rm{ \alpha }} = a{e^{{ \rm{b \alpha }}}} \sin \rm{ \alpha }\end{split} \tag{8} \end{align*} \end{document} \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} r = a{e^{{ \rm{b \alpha }}}} \tag{9} \end{align*} \end{document}

In differential geometric terms, the logarithmic spiral can also be defined with the constant angle \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\beta$$ \end{document} (Eq. 10), measured from the tangent and radial vectors at every single point of the curve ( \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$r , { \rm{ \; \alpha }}$$ \end{document} ). The derivative of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$r \left( { \rm{ \alpha }} \right)$$ \end{document} is proportional to the b parameter. In other words, it determines how closely and in which direction the spiral wraps. In the extreme case, where b is equal to zero, the spiral would be a circle of radius a, and if it approached to infinity, the spiral would tend to be a ray. \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} \beta = { \cos ^ { - 1 } } { \frac { \left\langle { r \left( { \rm { \alpha } } \right) , \;r \prime \left( { \rm { \alpha } } \right) } \right\rangle } { \| { r \left( { \rm { \alpha } } \right) } \| \, \, \| { \;r \prime \left( { \rm { \alpha } } \right) } \| } } = { \rm { \; } } { \tan ^ { - 1 } } \frac { 1 } { b } \tag { 10 } \end{align*} \end{document}

When the robot has an initial upright configuration, all of its structure temporally belongs to a unidimensional space. On that basis, when the robot is moved using the natural-CCD algorithm and with no orientation restrictions, the joints will belong to a plane. In this scenario, some of the resulting end-effector paths tend to be logarithmic spirals. This can be tested using the Cesàro equation. In geometry of plane curves, this equation provides an indicator of the curvature, ҡ, which is the reciprocal of the radius. It can be said that two curves are congruent if they have the same Cesàro equation. In geometry, two figures or objects are congruent if they have the same shape and size, or if one has the same shape and size as the mirror image of the other. The Cesàro equations of some types of curves are the ones as follows (Eqs. 11–14), where R is the radius, C is a constant, and s is the arc length³³: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} { \rm{Straight \ line :}} \;{ { {\upkappa} }} =0 \tag{11} \end{align*} \end{document} \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} { \rm { Circle : } } \; { \rm { { \upkappa } } } = \frac { 1 } { R } \tag { 12 } \end{align*} \end{document} \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} { \rm { Circular \ involute : } } \; { { { \upkappa } } } = \frac { C } { { \sqrt s } } \tag { 13 } \end{align*} \end{document} \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} { \rm { Logarithmic \ spiral : } } \; { { { \upkappa } } } = \frac { C } { s } \tag { 14 } \end{align*} \end{document}

If the curve produced by the end effector is parameterized, it is possible to represent the relationship between the curvature, ҡ, with the arc length, s, as shown in Figure 4. For example, the particular movement of the previous Figure 3 tends to a logarithmic spiral since the arc length and the curvature are inversely proportional (Eq. 14).

FIG. 4.

Curvature-arc length test results for the movement of a robot with 10 links using the natural-CCD algorithm from an initial upright configuration and with a desired end-effector position close to the anchored base.

Moreover, it has been found that this particular path is, in fact, a specific type of logarithmic spiral: a golden spiral. This spiral is characterized by having a growth factor equal to \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\varphi$$ \end{document} , the golden ratio or divine proportion. In mathematics, two quantities are in the golden ratio if their proportion is the same as the ratio of their sum to the larger of the two quantities (Eqs. 15–16). \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} \frac { { u + v } } { u } = \frac { u } { v } \; \; { \rm { with } } \ u > v \tag { 15 } \end{align*} \end{document} \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} & \varphi = {u \over v} \\ & 1 + { \varphi ^{ - 1}} = \varphi \\ & { \varphi ^2} - \varphi - 1 = 0 \\ & \varphi \; = \;1.61803398 \ldots\ \tag{16} \end{align*} \end{document}

A widely known approximation of the golden spiral is the Fibonacci spiral. This curve has numerous applications in number theory, applied mathematics, computer science, physics, and biology.³⁴ Its popularity comes from its appearance in numerous figures of nature such as branches of trees, arrangement of leaves on the stem, flowers of artichokes and sunflowers,³⁵ inflorescences of Romanesco broccoli, horns and claws of some mammals, phalanges of human hands, sea waves, tails of comets, or even arms of galaxies. The Fibonacci spiral is generated drawing circular arcs that connect the opposite corners of several squares beside each other, having as sides the values of the Fibonacci sequence (Eq. 17). \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} x \left( {k + 2} \right) = x \left( {k + 1} \right) + x \left( k \right) \tag{17} \end{align*} \end{document}

Owing to this concept, the golden spiral and consequently its approximation, the Fibonacci spiral, can be defined according to the Equations (10, 18–20). \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} {e^{{b_f}{{ \rm{ \alpha }}_{right}}}} = \varphi \tag{18} \end{align*} \end{document} \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} \left\vert { { b_f } } \right\vert = \; { \frac { \ln \varphi } { \pi / 2 } } \; = 0.306349 \;rad \tag { 19 } \end{align*} \end{document} \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} { \beta _f } = { \tan ^ { - 1 } } \frac { 1 } { { { b_f } } } = { 72.9676^o } \tag { 20 } \end{align*} \end{document}

It can be demonstrated that the resulting path from Figure 3 is a Fibonacci spiral by polynomial curve fitting. It has been found that this appears for a sufficiently redundant number of equal-length links according to the Equation (20).

It can be seen in Figure 5 that for a number higher than 10 links, the resulting path tends to be a golden spiral (Fig. 6). There exists a slight difference between the stabilization of the theoretical constant angle \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${ \beta _f}$$ \end{document} and the obtained from the resulting path \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\beta$$ \end{document} , because the robot is discretized in a finite number of serial links, instead of being infinitely bendable on all of its points. Due to the application of the k-factor presented in Equation (6), all the joints are rotating nearly the same small amount in each cycle, converging to the desired position. Because of this, from a kinematics perspective, it is possible to set forth the following statement.

FIG. 5.

Comparison of the Fibonacci spiral constant angle β_f and the obtained angle β from the resulting path (Fig. 3c) versus the number of links of the robot.

FIG. 6.

(a) The golden spiral. (b) The Fibonacci spiral. (c) The resulting end-effector path using the natural-CCD algorithm.

A serial manipulator robot with a sufficiently redundant number of equal-length links can determine a biomimetic end-effector path, or Fibonacci spiral, by using the Natural-CCD algorithm, when their rotational joints are cyclically and uniformly moved from an upright configuration bit by bit to a position close to its base.

The demonstration of the biomimetic movement of the robot in Figure 3 gives reason to think that the harmonious movements originated by the natural-CCD algorithm are, to some extent, biomimetic or related to nature. To generalize that theory, the proposal is to imagine a robot with an infinite number of links. In fact, this could be a soft robot since it is continuous. In that context, the same statement can be expressed in geometric terms as follows.

If a straight segment with a fixed length and anchored on one end is cyclically, uniformly and infinitesimally bended on all of its points along the same direction, the other end will move along a Golden spiral until enclosing the initial straight segment to a circle.

It is expected that this property can help to clarify the common appearance of the Fibonacci spiral in numerous figures of nature. After all, it is related to the behavior of a very simple and familiar movement. Another way to visualize this phenomenon is by approximating the robot configurations to arcs of circumferences. Initially, since the configuration is a straight segment, it belongs to a circumference with an infinite radius. The length of the robot remains invariable with the movement, but it does vary the radius and the central angle of the circumferences that contain the robot on each cycle (Fig. 7).

FIG. 7.

The Fibonacci spiral path and representation of the circumferences around the robot configurations.

Application to soft robotics

Previously, a new method called natural-CCD has been presented to solve the inverse kinematics problem of hyper-redundant robots. It has been shown that this particular algorithm presents solutions that are, to some extent, biomimetic or related to nature. Soft robots are also inspired by nature.³⁶ Usually, they should be taken as continuous because of their intrinsic morphology, materials, and actuation. Some authors affirm that soft robotic manipulators are kinematically similar to hyper-redundant manipulators with a large number of DOFs, so models that approximate continuum manipulators by finite DOF hyper-redundant ones can be appropriate.¹³ The natural-CCD provides solutions for discrete hyper-redundant robots so far. Now, a transformation of those discrete solutions to continuous is introduced.

There are some studies that propose techniques to discretize a continuous curve to fit discrete robots.^37,38 However, in this section, we discuss the opposite methodology, which is similar to an interpolation problem: to find a continuous curve from a set of discrete points. The positions and orientations of each joint of the discrete and hyper-redundant robot are obtained in real time from the pose provided by the natural-CCD algorithm. With that information, a backbone curve can be defined.³⁹ It is a continuous curve that captures the important macroscopic geometry features of a hyper-redundant robot. Related to this, a backbone curve reference set consists of a backbone curve parameterization and an associated set of orthonormal frames that provide a complete representation of this kind of robots. Using this strategy, it is possible to convert the discrete solution to a continuous one for soft robots.

When working with hyper-redundant or soft robots, the pose of the whole body is usually as important as the pose of the end effector. A limitation of current studies is that neither the whole body nor the pose of the end effector is considered in the solution.¹⁷ In fact, whole-body movements for soft robots^40,41 are usually difficult to implement owing to their inherent compliance.

The natural-CCD algorithm offers a complete whole-body and end-effector solution for the inverse kinematics of soft robots. To achieve this, a new strategy is proposed. For a given desired and reachable whole-body configuration, q_d, if it is continuous it should be discretized to a sufficient number of links. Then, following previous analyses (Eqs. 3 and 4), it is proposed to change the values of p_e and p_f, while p_c remains the same (Eqs. 21–23). \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} {p_c} = \;p ( {q_i} ) \tag{21} \end{align*} \end{document} \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} {p_e} = \;p ( {q_{i + 1}} ) \tag{22} \end{align*} \end{document} \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} {p_f} = p ( {q_i} ) + \left( {p ( {q_{d , i + 1}}} \right) - p ( {q_{d , i}} ) ) \tag{23} \end{align*} \end{document}

where p_c is still the current position of the joint to be moved; p_e now takes the value of the current position of the following joint and p_f is now variable taken as the virtual projection of the desired joint position. This strategy allows vector \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\overrightarrow {{p_c} \;{p_f}}$$ \end{document} to be independent of the current joint and only dependent on the desired pose. Thus, it can be said that the natural-CCD algorithm provides a solution for the inverse kinematics problem and a set of whole-body movements that can be transformed from discrete to continuous for soft robotic manipulators under the assumption of having a representative number of joints.

Results

Description of experiments

The CCD algorithm provided a desired position and orientation error minimization for either rotational or prismatic joints.²⁷ This has been tested using the natural-CCD as well. Figure 8 shows the movement of a robot for a desired position using the natural-CCD algorithm and the variation of some important parameters.

FIG. 8.

Movement of a robot with three rotational joints using the natural-CCD from an upright configuration to a given desired point, where {q₁,q₂,…,q_N} are the robot joints from the base to the end effector. (a) Representation of the intermediate configurations. (b) k-Factor variation. (c) Convergence of the end effector to the desired position measured as a Euclidean distance. (d) Angle variation for each joint.

For a given position and orientation, the natural-CCD algorithm stands out for providing very harmonious and coherent solutions (Fig. 9).

FIG. 9.

Movement of a hyper-redundant robot with 10 rotational joints and 3 DOFs per joint using the natural-CCD. (a) Initial configuration and desired position and orientation for the end effector. (b) Final configuration in a 3D perspective. (c) Top view of the final configuration. 3D, three dimensional; DOF, degree-of-freedom.

In addition, the combined use of rotational and prismatic joints has also been studied, as shown in Figure 10.

FIG. 10.

Movement of a robot with five rotational joints and two rotational and prismatic joints (second and fifth from the base) from an upright configuration to a desired position using the natural-CCD. (a) Initial upright configuration. (b) Intermediate configurations. (c) Final configuration.

In addition to this, other important features have originally been studied in this work. For example, since it does not matter in what order the coordinates are minimized using the natural-CCD, it is possible to disable the movement of a given joint in real time to simulate a mechanical failure and exploit the error tolerance capabilities of the hyper-redundant robot (Fig. 11).

FIG. 11.

Simulation of the failure of 4 joints (fifth to the eighth from the base) over a total of 12 rotational joints to study the error tolerance of the natural-CCD algorithm for hyper-redundant robots. (a) Initial upright configuration. (b) Intermediate configurations. (c) Final configuration.

Furthermore, the inverse kinematics algorithm can be applied for a continuous trajectory. This trajectory can be discretized as a set of desired positions that have to be reached sequentially under a certain precision, instead of having a single desired position (Fig. 12). In this context, it is recommended to use the SCD algorithm instead of using the natural-CCD. With that strategy, the joint with more manipulability to move along a particular trajectory can be selected in each iteration, allowing to reduce significantly the number of iterations.

FIG. 12.

Trajectory planning of a robot with 20 rotational joints and 3 DOFs per joint (60 DOFs) using the SCD. (a) Desired straight-line trajectory. (b) Final configuration and end-effector path for the line trajectory. (c) Desired helical trajectory. (d) Final configuration and end-effector path for the helical trajectory.

Moreover, a study of the obstacle-avoidance capabilities of hyper-redundant robots using the natural-CCD has been made. The natural-CCD algorithm lets the end effector of the robot be attracted by the desired location. This can be combined with a robotic motion planning method, such as potential fields,⁴² to repel the obstacles in each iteration (Fig. 13). All the obstacles have been virtually swollen a slightly higher amount than half of the girth of the manipulator to prevent collisions. The main drawback of this strategy is that the algorithm can fall into a local minimum while minimizing the coordinates. To solve this issue, a tangential direction is determined to move the robot contouring the obstacle. This strategy can remind a Bug Algorithm,⁴³ a complete path planning method for mobile robots.

FIG. 13.

Combined use of natural-CCD and potential fields to avoid five obstacles modeled as consecutive spheres. (a) Initial robot configuration, obstacles, and desired position. (b) Final configuration and path.

A collision-free trajectory from the initial configuration to the desired position has been calculated when the obstacles and the environment are more complex, such as introducing the robot through a hole in a wall (Fig. 14). Thus, this trajectory will guide the end effector, while the rest of the robot is repelled from the obstacles if necessary. This trajectory cannot be any collision-free trajectory since the end effector is not a point mobile robot and it is linked to the joints toward its base. This approach avoids the use of the off-line calculation of the Configuration Space or C-Space,⁴⁴ which results in a high computational cost for hyper-redundant manipulators.

FIG. 14.

Motion planning for hyper-redundant robots in complex environments. No collisions occurred during the execution. (a) Initial robot configuration in an environment with two complex obstacles and representation of the desired position and a valid collision-free trajectory for the end effector. (b) Final robot configuration and path for the previous scenario. (c) Initial configuration, desired position, and collision-free trajectory for the end effector in other scenarios with walls and holes. (d) Final robot configuration and path for the previous scenario.

To conclude the simulation results, it has been tested all the previous capabilities of hyper-redundant robots at once to study their performance (Fig. 15). For example, a hyper-redundant robot with nine joints is studied. A mechanical failure in joints 2 and 3 has been introduced. The end-effector joint is both rotational and prismatic. In addition, nine consecutive spheres have been introduced in the working environment as obstacles that the robot has to avoid. The goal is to follow a linear trajectory with a desired orientation in all of its points with a precision of 0.1%. Despite all these constraints, the algorithm works successfully for the given task.

FIG. 15.

Use of the natural-CCD algorithm for resolution of a complex task with a lot of constraints. No collisions occurred during the execution. (a) Initial configuration of a robot with 9 rotational joints and 3 DOFs per joint. A mechanical failure is simulated in 2 joints (second and third from the base) and the last joint is also prismatic. The workspace includes 9 consecutive obstacles that have to be avoided. Then, the robot must conduct a linear trajectory with a given orientation. (b) Final robot configuration in a 3D perspective. (c) Final robot configuration in a front view.

Performance analysis

An application has been developed in MATLAB 2016a in Windows 10 to test the performance of the natural-CCD. A computer featuring an Intel^® Core i7-4720HQ (2.60 GHz) processor and 16GB of RAM memory were used to execute the algorithm. As previously mentioned in the Introduction, the most important criteria to analyze the effectiveness of an inverse kinematics algorithm for hyper-redundant robots are precision, computational time, and quality of the solution.

Precision

Using the natural-CCD algorithm, the end effector of the robot will always converge to the desired location, if reachable, under a certain precision. This is why this parameter has been taken as an input variable requirement, whereas the computational time is the critical one to measure. It is worth noting that the precision should be relative to the robot length. For example, it is not the same for a robot of 10 cm and a robot of 1 m length to reach a desired position with a precision of 1 mm. The precision in the first case should be 1% and in the second case 0.1%. Precisions up to 10⁻⁷% have been defined without significantly increasing the computational time when the natural-CCD algorithm is applied.

Computational time

This is usually the most critical variable when a hyper-redundant robot relies on a large number of DOFs. One of the main reasons to have chosen the CCD as starting point is because it is a very simple and computationally efficient method. The Monte Carlo method has been used to study the computational time against the number of joints. A set of 100 random initial robot configurations and desired locations have been taken inside the reachable 3D workspace—a solid sphere of radius the length of the robot—to measure the computational time. This has been done for robots with different number of joints and DOFs per joint. The precision was set 0.1% for all cases. Figure 16 shows that the computational times for 1, 2, or 3 DOFs per joint are quite similar. The differences among them are a result of the randomized input data. However, it is remarkable that the computational times are <1 s for robots of up to 300 DOFs for both desired positions and desired positions and orientations. When robots of up to 1000 DOFs are managed, it can be deduced that the computational time grows lineally with an increase of the number of joints, which is a huge advantage.

FIG. 16.

Mean computational time of the inverse kinematics problem resolution versus the number of joints for 1, 2, or 3 DOFs per joint (a) for given random positions (b) and for given random positions and orientations.

Quality of the solution

A good quality solution is one without self-collisions that maximizes the transmissibility of force, dexterity, rigidity, energy efficiency, or availability of joint range. This is, in other words, the one that best conveys the strength to the end effector, keeps each joint as close as possible to its central position, minimizes the speed, momentum, and inertia of the joints, avoids singularities, and maximizes the skill of the end effector. First of all, the natural-CCD ensures that the robot does not have any singularity or collide with itself (Eq. 5). Also, it provides an incremental and uniform movement for all the joints of the robot (Eq. 6). As a result, the robot shows an ability to easily move in a way that is biomimetic and fluid, giving place to coherent final configurations.

Discussion

The natural-CCD performance can be compared to other similar studies from the state-of-the-art both qualitatively and quantitatively. The other methods are as follows: Exhaustive methods,¹¹ Pattern Search,¹¹ Global Search,¹¹ Genetic Algorithms,¹¹ Simulated Annealing,¹¹ and Particle-Swarm Optimization.¹² The CCD algorithm has not been included in the following comparisons because it does not offer acceptable solutions in terms of singularities, self-collisions, and abrupt movements. First of all, due to the simplicity of the natural-CCD, it manages a much higher number of DOFs with acceptable computational times. It is a method with a very simple theoretical basis but with a high rate of convergence against other previous works. Despite the state-of-the-art methods, the precision of the presented algorithm is an input value that can be defined as a system requirement. It studies and manages the robot singularities, constraints, and the joint movements that lead to appropriate and coherent solutions. About the computational time, the critical variable, the natural-CCD, provides solutions for the inverse kinematics problem in <1 s for robots of up to 300 DOFs. Compared to the state-of-the-art methods that present computational times in order of seconds that increase exponentially against the number of joints, the natural-CCD presents far lower ones that grow lineally against the number of joints (Fig. 17).

FIG. 17.

Computational time comparison of the state-of-the-art methods and the natural-CCD against the number of DOFs in a logarithmic scale. PASO, PAth planning with Swarm Optimization.

In summary, the strengths of the algorithm compared to the state-of-the-art are simplicity, low computational cost, ability to establish highly desired precisions, and the quality of the solutions obtained. All of this could eventually lead to numerous potential benefits. First of all, the natural-CCD can be a viable solution to solve the inverse kinematics problem of hyper-redundant robots with real-time requirements. Second, it is based on biomimetic movements that help to understand some shapes of nature such as the Fibonacci spiral. And third, it proposes a method to convert the obtained discrete solutions to continuous for soft robotic applications.

Conclusion

In this document, soft and hyper-redundant robots have been studied because of their potential benefits against conventional manipulators. They can perform better obstacle avoidance, are error-tolerant, and offer mechanic and kinematic advantages. The main problem of this particular kind of robots is to effectively solve the inverse kinematics, since there exist an infinite number of possible solutions that satisfy that problem. A new algorithm based on the CCD and entitled as natural-CCD is proposed to solve this issue. This algorithm shows coherent and harmonious solutions, result of the graceful movements of the joints. Using this algorithm, a robot with a sufficiently redundant number of equal-length links can lead to an end-effector path that is very present in nature: a Fibonacci spiral. This theory has been extrapolated in geometrical terms to clarify the origin of many shapes in nature that tend to that specific spiral. Also, while the state-of-the-art methods show high computational times that grow exponentially against the number of DOFs, the natural-CCD presents far lower ones that could allow real-time applications, even for a large number of joints. Furthermore, the precisions using this method are a user input value and not part of a conditioned result. Finally, the discrete solutions obtained by this method have been extrapolated to continuous ones for soft robotics. Unlike the state-of-the-art methods, it allows whole-body and end-effector movements for soft robots. Because all of this, it can be stated that the natural-CCD algorithm is eventually a simple, efficient, and viable solution to solve the inverse kinematics problem of hyper-redundant and soft robots.

Footnotes

Acknowledgments

This work is the result of the research activities carried out at the Centre of Automation and Robotics (CAR-UPM-CSIC) of the Technical University of Madrid (Spain), within the Robotics and Cybernetics research group (RobCib). This work is supported by the project RoboCity2030-III-CM (Robótica aplicada a la mejora de la calidad de vida de los ciudadanos. Fase III; S2013/MIT-2748) and also by the Project DPI2014-56985-R (Protección Robotizada de Infraestructuras Críticas), funded by the Ministry of Economy and Competitiveness of the Government of Spain.

Author Disclosure Statement

No competing financial interests exist.

References

Barrientos

, et al. Fundamentos de Robótica. Madrid, Spain: McGraw-Hill, 2007.

Kimura

, et al. A fault-tolerant control algorithm having a decentralized autonomous architecture for space hyper-redundant manipulators. IEEE Trans Syst Man Cybern A Syst Hum, 1998; 28:521–527.

, Simaan

. Actuation compensation for flexible surgical snake-like robots with redundant remote actuation. In 2006 IEEE International Conference on Robotics and Automation (ICRA), Vols 1–10, Orlando, FL, 2006, pp. 4148–4154.

Freire

MAL

, et al. Development and implementation of a singularity management strategy for a cooperative redundant robotic system destined to assistance during surgical interventions. RIAI, 2015; 12:80–91.

Sofge

, Chiang

. Design, implementation, and cooperative coevolution of an autonomous/teleoperated control system for a serpentine robotic manipulator. In Amer. Nucl. Soc. 9th Top. Meeting on Robotics and Remote Systems, 2001, arXiv: 0706.1061, 2007.

, Hirose

, Yoshinada

. Development of a hyper-redundant multijoint manipulator for maintenance of nuclear-reactors. Adv Robot, 1995; 9:281–300.

Sujan

, Lichter

, Dubowsky

. Lightweight hyper-redundant binary elements for planetary exploration robots. In Proceedings of the 2001 IEEE/ASME International Conference on Advanced Intelligent Mechatronics, 2001, Como, Italy, 2001.

Chirikjian

, Burdick

. A hyper-redundant manipulator. IEEE Robot Autom Mag, 1994; 1:22–29.

Chirikjian

. Theory and applications of hyper-redundant robotic manipulators (PhD dissertation). California Institute of Technology, 1992.

10.

Marques

, et al. 3D hyper-redundant robot. In 11th Spanish Portuguese Conference on Electrical Engineering, Zaragoza, Spain, 2009.

11.

Espinoza

, et al. Inverse kinematics of a 10 DOF modular hyper-redundant robot resorting to exhaustive and error-optimization methods: a comparative study. In Robotics Symposium and Latin American Robotics Symposium (SBR-LARS), 2012 Brazilian, 2012.

12.

Collins

, Shen

W-M

. PASO: an integrated, scalable PSO-based optimization framework for hyper-redundant manipulator path planning and inverse kinematics. Information Sciences Institute, University of Southern California, Tech. Rep. ISI-TR-697, 2015.

13.

Trivedi

, et al. Soft robotics: Biological inspiration, state of the art, and future research. Appl Bionics Biomech, 2008; 5:99–117.

14.

Chirikjian

. A continuum approach to hyper-redundant manipulator dynamics. In Proceedings of the 1993 IEEE/RSJ International Conference on Intelligent Robots and Systems' 93, IROS'93, 1993.

15.

Chirikjian

. Hyper-redundant manipulator dynamics: a continuum approximation. Adv Robot, 1994; 9:217–243.

16.

Chirikjian

. Inverse kinematics of binary manipulators using a continuum model. J Intell Robot Syst, 1997; 19:5–22.

17.

Rus

, Tolley

. Design, fabrication and control of soft robots. Nature, 2015; 521:467.

18.

Van Leeuwen

, Kier

. Functional design of tentacles in squid: linking sarcomere ultrastructure to gross morphological dynamics. Philos Trans R Soc Lond B Biol Sci, 1997; 352:551–571.

19.

Yekutieli

, et al. Dynamic model of the octopus arm. I. Biomechanics of the octopus reaching movement. J Neurophysiol, 2005; 94:1443–1458.

20.

Wilson

, et al. Flexible robot manipulators and grippers: relatives of elephant trunks and squid tentacles. In: Dario

, Sandini

, Aebischer

(eds). Robots and Biological Systems: Towards a New Bionics?. Berlin, Heidelberg: Springer, 1993, pp. 475–494.

21.

Reinhart

, Steil

. Hybrid mechanical and data-driven modeling improves inverse kinematic control of a soft robot. Proc Technol, 2016; 26:12–19.

22.

Duriez

. Control of elastic soft robots based on real-time finite element method. In 2013 IEEE International Conference on Robotics and Automation (ICRA), Karlsruhe, Germany, 2013.

23.

Wang

, et al. Visual servo control of cable-driven soft robotic manipulator. In 2013 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), Tokyo, Japan. 2013.

24.

Mutlu

, Alici

, Li

. An effective methodology to solve inverse kinematics of electroactive polymer actuators modelled as active and soft robotic structures. Mech Mach Theory, 2013; 67: 94–110.

25.

Neppalli

, et al. Closed-form inverse kinematics for continuum manipulators. Adv Robot, 2009; 23:2077–2091.

26.

Jones

, Walker

. Kinematics for multisection continuum robots. IEEE Trans Robot, 2006; 22:43–55.

27.

Wang

L-C

, Chen

C-C

. A combined optimization method for solving the inverse kinematics problems of mechanical manipulators. IEEE Trans Robot Autom, 1991; 7:489–499.

28.

Saha

, Tewari

. On the finite time convergence of cyclic coordinate descent methods. arXiv preprint arXiv: 1005.2146, 2010.

29.

Muller-Cajar

, Mukundan

. Triangualation: a new algorithm for inverse kinematics. In Proceedings of Image and Vision Computing New Zealand 2007, New Zealand, 2007, pp. 181–186.

30.

Wright

. Coordinate descent algorithms. Math Program, 2015; 151:3–34.

31.

Loshchilov

, Schoenauer

, Sebag

. Adaptive coordinate descent. In Proceedings of the 13th annual conference on Genetic and Evolutionary Computation, Dublin, Ireland, 2011, pp. 885–892.

32.

Hargroves

, Smith

. Innovation inspired by nature: Biomimicry. ECOS, 2006; 129:27–29.

33.

Lawrence

. A Catalog of Special Plane Curves. New York: Courier Corporation, 2013.

34.

Conway

, Guy

. The Book of Numbers. New York: Copernicus, 1996.

35.

Jean

. Phyllotaxis: A Systemic Study in Plant Morphogenesis. New York: Cambridge University Press, 2009.

36.

Kim

, Laschi

, Trimmer

. Soft robotics: a bioinspired evolution in robotics. Trends Biotechnol, 2013; 31:287–294.

37.

Hatton

, Choset

. Generating gaits for snake robots: annealed chain fitting and keyframe wave extraction. Auton Robots, 2010; 28:271–281.

38.

Andersson

. Discretization of a continuous curve. IEEE Trans Robot, 2008; 24:456–461.

39.

Burdick

, Chirikjian

. The kinematics of hyper-redundant robots. In Baillieul J, Sastry SS, Sussmann HJ (eds). Essays on Mathematical Robotics. New York: Springer, 1998, pp. 61–89.

40.

Marchese

, et al. Design and control of a soft and continuously deformable 2d robotic manipulation system. In 2014 IEEE International Conference on Robotics and Automation (ICRA), Hong Kong, China, 2014.

41.

Marchese

, Katzschmann

, Rus

. Whole arm planning for a soft and highly compliant 2d robotic manipulator. In 2014 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS 2014), Chicago, IL, 2014.

42.

Khatib

. Real-time obstacle avoidance for manipulators and mobile robots. Int J Robot Res, 1986; 5:90–98.

43.

Lumelsky

, Stepanov

. Path-planning strategies for a point mobile automaton moving amidst unknown obstacles of arbitrary shape. Algorithmica, 1987; 2:403–430.

44.

Lozano-Perez

. Spatial planning: a configuration space approach. IEEE Trans Comput, 1983; 100:108–120.