Inverse Kinematics of a Humanoid Robot with Non-Spherical Hip: A Hybrid Algorithm Approach

Abstract

This paper describes an approach to solve the inverse kinematics problem of humanoid robots whose construction shows a small but non negligible offset at the hip which prevents any purely analytical solution to be developed. Knowing that a purely numerical solution is not feasible due to variable efficiency problems, the proposed one first neglects the offset presence in order to obtain an approximate “solution” by means of an analytical algorithm based on screw theory, and then uses it as the initial condition of a numerical refining procedure based on the Levenberg-Marquardt algorithm. In this way, few iterations are needed for any specified attitude, making it possible to implement the algorithm for real-time applications. As a way to show the algorithm's implementation, one case of study is considered throughout the paper, represented by the SILO2 humanoid robot.

1. Introduction

The kinematic modeling of humanoid robots is in fact a more difficult problem when compared to the classical industry manipulators, due to the locomotion that these robots are expected to perform while maintaining, at the same time, their postural stability. It is a formidable computational challenge, especially for real-time applications, due to the high number of degrees of freedom, the complex kinematics configurations that may lead to non-analytical solutions and the constantly changing point of support that is not rigidly attached to the floor [1]. In order to control the motion of these humanoid robots, specific trajectories in task space are required; that is, a desired sequence of attitudes for the feet and the body are required for accomplishing the desired stable locomotion task.

In the context of rigid multi-body systems, as is the case of robots, task space is different from joint space where motor commands are issued. Hence, movement planning requires either a preset sequence generated with play-back techniques or a suitable coordinate transformation from task to joint space. The former method is very simple as it does not require a complex mathematical analysis; however, it is not systematic as it is based on empirical tests for each movement, which make it tedious and non practical. The latter one, based on coordinates transformation, is supposed to be found by solving the so called Inverse Kinematics (IK) problem, which arises from the fact that the inverse transformation is often complex and ill-posed [2].

It is known that the IK problem cannot be solved by means of systematic processes as it is highly dependent on the robot's kinematic structure, and that it also either yields to infinite, multiple, one or none solutions. Only in simple cases, as in the case when the kinematic chain has no more than 6 degrees-of-freedom (dof) and three joint axes intersect, the IK problem for position and orientation can be decoupled, as stated by the so popular method proposed by Pieper in his PhD thesis [3, 4]. This is the case of many industrial manipulators and some of the most popular humanoid robots, i.e. ASIMO of Honda, HRP-2 of AIST and Kawada Industries or HOAP-2 of Fujitsu [5]. However, even for singularity-free mechanisms, the presence of offsets prevents this fact. Therefore, closed-form solutions are difficult, if not impossible, to find [6]. Or if they exist, they are very specific for a given kinematic chain [7, 8, 9].

For this reason, numerical methods have been commonly used to solve the IK problem, by using the pseudoinverse of the Jacobian (or a modification of it) [2, 10], the Cyclic Coordinate Descent (CCD) method [11, 12] or a procedure based on iterative triangulation [13, 14]. Some authors also use neural approaches [15, 16] or computational approaches based on conformal geometric algebra [17] or symbolic preprocessing [18]. However, these ones cannot be as efficient as the analytical ones, as they may require several iterations (not known a priori) in order to reach a particular precision, leading in many cases to computationally expensive solutions, which can be numerically unstable [19, 20]. Not to mention that a numerical algorithm may converge to a local minima instead of the desired global one, represented by the desired solution [21].

This makes it even more difficult to develop real-time applications, as these ones require fast and efficient algorithms for computing the gait patterns required to achieve dynamical stability. Therefore, finding an analytical solution to the IK problem, or a solution that takes a minimum of iterations, makes it possible to develop controllers able to compensate dynamical effects by changing motion patterns in real-time.

2. SILO2 Humanoid Robot

SILO2 stands for Locomotion System based on two limbs. This humanoid robot belongs to a family of walking prototypes that have been developed in the Department of Automatic Control at the Institute of Industrial Automation in Madrid, Spain [22].

This one in particular is currently on development and it is planned to be upgraded until getting a complete humanoid kinematic structure of 23 dof. However, in the current stage only the locomotor unit has been entirely developed (Figure 1).

Figure 1.

SILO2 Humanoid Robot [23]

As can be seen on this picture, this robot possesses 6 dof per leg: 3 at the hip, 1 at the knee, and 2 at the ankle. That is, each leg's mechanism is determined, so that each foot (or the body, depending on the basement link's election) can reach any desired attitude inside its workspace, defined by geometrical and mechanical constraints.

Joints are driven with two different types of actuators, as can be also noticed: Standard DC servomotors placed at the joint axis and non-linear transmission ratio actuators, used to improve the mechanical design by decreasing the energy consumption. The latter ones, patented as SMART (Special Mechatronic Actuator for Robot JoinTs) by the Department of Automatic Control at the Institute of Industrial Automation, are implemented as standard DC servomotors and four-bar linkages which provide the so called non-linear transmission ratio between the input angle driven by the servo and the output angle related to the joint variable.

Actually, three joints on each leg are driven by SMART actuators: One at the hip (abduction / adduction), one at the knee (flexion / extension), and one at the toe (flexion / extension), as can be seen on Figure 2a, 2b, and 2c.

Figure 2.

SMART actuators at (a) the hip, (b) the knee, (c) the ankle [23]

3. Humanoid Forward Kinematics

3.1. Kinematic diagram

In order to build the kinematic diagram of this robot it is worth to notice that the presence of SMART actuators does not add complexity to the kinematic model, as could be thought in a first sight. The four-bar linkages are just transmission mechanisms placed between the servomotors and the actual joints, such that the joint variable is indeed a function of the corresponding servomotor's angle.

Keeping this in mind, the kinematic structure of this robot can be represented with the diagram shown in Figure 3. In this diagram, joints are represented by cylinder-shaped icons, whose arrow represents its positive sense of motion. Each joint is labeled according to the following nomenclature: ${\mathrm{J}}_{si},$ where

Figure 3.

SILO2's kinematic diagram

$\mathrm{s}:\mathrm{l}\left(left\right)\text{\hspace{0.17em}or\hspace{0.17em}\hspace{0.17em}}\mathrm{r}\left(right\right).$

$\mathrm{i}:$ A number starting at 1 from the joint next to the body frame (base frame) and increasing until reaching the last joint in each kinematic chain.

Joint variables are labeled as ${\theta }_{si}$ according to the corresponding joint ${\mathrm{J}}_{si},$ such that each leg has a configuration determined by the vector ${\theta }_s=\left\{{\theta }_{\mathrm{s}1},{\theta }_{\mathrm{s}2},\dots ,{\theta }_{\mathrm{s}6}\right\}.$

Links follow a similar nomenclature, ${\mathrm{L}}_{si},$ where ${\mathrm{L}}_{\mathrm{s}0}$ refers to the base link (the body). In this way, joint ${\mathrm{J}}_{si}$ directly drives link ${\mathrm{L}}_{si}.$

Link dimensions are displayed along with the corresponding constant denoted by ${\mathrm{l}}_{\mathrm{i}},$ where side information is omitted due to symmetry reasons.

The body frame $\left\{\mathrm{D}\right\}$ of the robot is oriented in such a way that the y axis points up and the x axis points in its movement direction. The z axis is then established to complete a right-hand frame.

At the reference attitude $\left({\theta }_s=0\right),$ the feet frames (referred as $\left\{{\mathrm{S}}_{\mathrm{F}}\right\}$ ) have the same orientation of the body frame, so that their configuration with respect to the frame $\left\{\mathrm{D}\right\},{\bar{g}}_{\mathrm{d},sf}\left(0\right)$ [24], can be described as

${\bar{g}}_{\mathrm{d},sf}\left(0\right)=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& -\left({\mathrm{l}}_1+{\mathrm{l}}_2+{\mathrm{l}}_3+{\mathrm{l}}_4+{\mathrm{l}}_{\mathrm{f}}\right)\\ 0& 0& 1& \zeta {\mathrm{l}}_0\\ 0& 0& 0& 1\end{array}\right]$ (1)

Where

$\left(\mathrm{s}=\mathrm{r}\right)\Rightarrow \left(\zeta =1\right)\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\left(\mathrm{s}=\mathrm{l}\right)\Rightarrow \left(\zeta =-1\right)$ (2)

3.2. Screw motion

Forward kinematics can be described in a more geometric way if a screw motion is associated to every joint [24]. A screw motion represents a rotation by an amount of θ about an axis l followed by a translation by an amount d parallel to that axis, as can be seen on Figure 4. This axis is represented as a directed line ω through a point $q;$ such that $l=\left\{q+\lambda \omega |\lambda \in ℝ\right\}.$ The ratio $\mathrm{h}=\mathrm{d}/\theta$ is known as the pitch of the screw.

Figure 4.

Generalized Screw Motion [24]

In this way, if a point p is operated with the equivalent rigid motion transformation $\mathrm{g},$ its final location (Figure 4) is given by

$\mathrm{g}p=q+{\mathrm{e}}^{\stackrel{⌢}{\omega }\theta }\left(p-q\right)+\mathrm{h}\theta \omega$ (3)

Where ${\mathrm{e}}^{\stackrel{⌢}{\omega }\theta }$ is the exponential matrix of the skew-symmetric matrix $\stackrel{⌢}{\omega }$ corresponding to the unit vector $\omega ;$ that is, ${\mathrm{e}}^{\stackrel{⌢}{\omega }\theta }$ is just another way for representing the rotation matrix $\mathrm{R}\left(\omega ,\theta \right).$

By representing (3) in homogeneous coordinates, we have

$\bar{g}\hspace{0.17em}\bar{p}=\left[\begin{array}{cc}{\mathrm{e}}^{\stackrel{⌢}{\omega }\theta }& \left(\mathrm{I}-{\mathrm{e}}^{\stackrel{⌢}{\omega }\theta }\right)q+\mathrm{h}\theta \omega \\ 0^T& 1\end{array}\right]\left[\begin{array}{c}p\\ 1\end{array}\right]$ (4)

The transformation $\bar{g}$ in (4) can also be represented as an exponential matrix, $\bar{g}={\mathrm{e}}^{\stackrel{⌢}{\xi }\theta },$ of a matrix defined as

$\stackrel{⌢}{\xi }=\left[\begin{array}{cc}\stackrel{⌢}{\omega }& -\omega \times q+\mathrm{h}\omega \\ 0^T& 1\end{array}\right]$ (5)

In this way, a screw motion associated to a rotational joint $\left(\mathrm{h}=0\right)$ with joint variable θ is given by

${\mathrm{e}}^{\stackrel{⌢}{\xi }\theta }=\left[\begin{array}{cc}{e}^{\stackrel{⌢}{\omega }\theta }& \left(\mathrm{I}-{\mathrm{e}}^{\stackrel{⌢}{\omega }\theta }\right)q\\ 0^T& 1\end{array}\right]$ (6)

3.3. Product of exponentials formula

The transformation ${\bar{g}}_{\mathrm{d},sf}\left(0\right)$ shown in (1) describes the attitude of the corresponding foot with respect to the body frame.

Suppose that we fix every joint except the last one $\left({\mathrm{J}}_{\mathrm{s}6}\right),$ which directly affects the foot's attitude, then the new attitude can be calculated as

${\bar{g}}_{\mathrm{d},sf}\left({\theta }_{\mathrm{s}6}\right)={\mathrm{e}}^{{\stackrel{⌢}{\xi }}_{\mathrm{s}6}{\theta }_{\mathrm{s}6}}{\bar{g}}_{\mathrm{d},sf}\left(0\right)$ (7)

In the same way, we can fix every joint except ${\mathrm{J}}_{\mathrm{s}5},$ knowing that the foot's attitude has been modified by the ${\mathrm{J}}_{\mathrm{s}6}\text{'}\mathrm{s}$ action, then

$\begin{array}{l}{\bar{g}}_{\mathrm{d},sf}\left({\theta }_{\mathrm{s}5},{\theta }_{\mathrm{s}6}\right)={\mathrm{e}}^{{\stackrel{⌢}{\xi }}_{\mathrm{s}5}{\theta }_{\mathrm{s}5}}{\bar{g}}_{\mathrm{d},sf}\left({\theta }_{\mathrm{s}6}\right)\\ {\bar{g}}_{\mathrm{d},sf}\left({\theta }_{\mathrm{s}5},{\theta }_{\mathrm{s}6}\right)={\mathrm{e}}^{{\stackrel{⌢}{\xi }}_{\mathrm{s}5}{\theta }_{\mathrm{s}5}}{\mathrm{e}}^{{\stackrel{⌢}{\xi }}_{\mathrm{s}6}{\theta }_{\mathrm{s}6}}{\bar{g}}_{\mathrm{d},sf}\left(0\right)\end{array}$ (8)

By composition, after taking into account every joint in the corresponding kinematic chain the foot configuration becomes

${\bar{g}}_{\mathrm{d},sf}\left({\theta }_s\right)={\mathrm{e}}^{{\stackrel{⌢}{\xi }}_{\mathrm{s}1}{\theta }_{\mathrm{s}1}}{\mathrm{e}}^{{\stackrel{⌢}{\xi }}_{\mathrm{s}2}{\theta }_{\mathrm{s}2}}\dots {\mathrm{e}}^{{\stackrel{⌢}{\xi }}_{\mathrm{s}6}{\theta }_{\mathrm{s}6}}{\bar{g}}_{\mathrm{d},sf}\left(0\right)$ (9)

The expression in (9) is called the product of exponentials formula for the humanoid forward kinematics; that is, it describes the foot's attitude as a function of every joint in the kinematic chain [24].

3.4. Forward Kinematics parameters

As can be seen from expressions (9) and (6), in order to calculate the forward kinematics it is necessary to specify each exponential matrix ${\mathrm{e}}^{{\stackrel{⌢}{\xi }}_{si}{\theta }_{si}};$ that is, apart from the joint variable we need to define each axis direction ${\omega }_{si}$ and any point on each axis $q_{si},$ as described in the body's reference frame $\left\{\mathrm{D}\right\}.$

From the kinematic diagram shown in Figure 3, each axis direction is defined as

${\omega }_{s1}={\left[\begin{array}{ccc}0& 1& 0\end{array}\right]}^{\mathrm{T}}$ (10) ${\omega }_{s2}={\left[\begin{array}{ccc}0& 0& -\zeta \end{array}\right]}^{\mathrm{T}}$ (11) ${\omega }_{s3}={\left[\begin{array}{ccc}-1& 0& 0\end{array}\right]}^{\mathrm{T}}$ (12) ${\omega }_{s4}={\left[\begin{array}{ccc}0& 0& \zeta \end{array}\right]}^{\mathrm{T}}$ (13) ${\omega }_{s5}={\left[\begin{array}{ccc}0& 0& \zeta \end{array}\right]}^{\mathrm{T}}$ (14) ${\omega }_{s6}={\left[\begin{array}{ccc}1& 0& 0\end{array}\right]}^{\mathrm{T}}$ (15)

And a point on each axis is proposed as

$q_{s1}={\left[\begin{array}{ccc}0& 0& \zeta {\mathrm{l}}_0\end{array}\right]}^{\mathrm{T}}$ (16) $q_{s2}={\left[\begin{array}{ccc}0& -{\mathrm{l}}_1& 0\end{array}\right]}^{\mathrm{T}}$ (17) $q_{s3}={\left[\begin{array}{ccc}0& -\left({\mathrm{l}}_1+{\mathrm{l}}_2\right)& \zeta {\mathrm{l}}_0\end{array}\right]}^{\mathrm{T}}$ (18) $q_{s4}={\left[\begin{array}{ccc}0& -\left({\mathrm{l}}_1+{\mathrm{l}}_2+{\mathrm{l}}_3\right)& 0\end{array}\right]}^{\mathrm{T}}$ (19) $q_{s5}={\left[\begin{array}{ccc}0& -\left({\mathrm{l}}_1+{\mathrm{l}}_2+{\mathrm{l}}_3+{\mathrm{l}}_4\right)& 0\end{array}\right]}^{\mathrm{T}}$ (20) $q_{s6}={\left[\begin{array}{ccc}0& -\left({\mathrm{l}}_1+{\mathrm{l}}_2+{\mathrm{l}}_3+{\mathrm{l}}_4\right)& \zeta {\mathrm{l}}_0\end{array}\right]}^{\mathrm{T}}$ (21)

3.5. Orientation coordinates

Consider the representation of the humanoid's foot (and its frame) as shown in Figure 5. The chosen orientation coordinates are three angles, taken as follows: First, consider a rotation of ${\theta }_{\mathrm{s}}$ about the Z axis (pitch), followed by a rotation of ${\theta }_{\mathrm{f}}$ about the X axis of the current frame (roll). Finally, consider a rotation of ${\theta }_{\mathrm{t}}$ about the Y axis (yaw) of the current frame; that is, about a normal vector to the foot's sole. In this way, the resulting rotation matrix ${\mathrm{R}}_{foot}$ is given by

Figure 5.

Foot's frame definition

$\begin{array}{l}{\mathrm{R}}_{\text{foot}}=\mathrm{R}\left(\mathrm{Z},{\theta }_{\mathrm{s}}\right)\mathrm{R}\left(\mathrm{X},{\theta }_{\mathrm{f}}\right)\mathrm{R}\left(\mathrm{Y},{\theta }_{\mathrm{t}}\right)\\ {\mathrm{R}}_{\text{foot}}=\left[\begin{array}{ccc}{\mathrm{c}}_{\mathrm{s}}{\mathrm{c}}_{\mathrm{t}}-{\mathrm{s}}_{\mathrm{s}}{\mathrm{s}}_{\mathrm{f}}{\mathrm{s}}_{\mathrm{t}}& -{\mathrm{s}}_{\mathrm{s}}{\mathrm{c}}_{\mathrm{f}}& {\mathrm{c}}_{\mathrm{s}}{\mathrm{c}}_{\mathrm{t}}+{\mathrm{s}}_{\mathrm{s}}{\mathrm{s}}_{\mathrm{f}}{\mathrm{c}}_{\mathrm{t}}\\ {\mathrm{s}}_{\mathrm{s}}{\mathrm{c}}_{\mathrm{t}}-{\mathrm{c}}_{\mathrm{s}}{\mathrm{s}}_{\mathrm{f}}{\mathrm{s}}_{\mathrm{t}}& {\mathrm{c}}_{\mathrm{s}}{\mathrm{c}}_{\mathrm{f}}& {\mathrm{s}}_{\mathrm{s}}{\mathrm{s}}_{\mathrm{t}}-{\mathrm{c}}_{\mathrm{s}}{\mathrm{s}}_{\mathrm{f}}{\mathrm{c}}_{\mathrm{t}}\\ -{\mathrm{c}}_{\mathrm{f}}{\mathrm{s}}_{\mathrm{t}}& {\mathrm{s}}_{\mathrm{f}}& {\mathrm{c}}_{\mathrm{f}}{\mathrm{c}}_{\mathrm{t}}\end{array}\right]\end{array}$ (22)

Such that, from an arbitrary rotation matrix these rotation coordinates can be calculated as

${\theta }_{\mathrm{f}}=\arctan 2\left({\mathrm{r}}_{32},\pm \sqrt{1-{\mathrm{r}}_{32}^2}\right)$ (23) ${\theta }_{\mathrm{s}}=\arctan 2\left(-{\mathrm{r}}_{12},{\mathrm{r}}_{22}\right)$ (24) ${\theta }_{\mathrm{t}}=\arctan 2\left(-{\mathrm{r}}_{31},{\mathrm{r}}_{33}\right)$ (25)

4. Humanoid Inverse Kinematics

By looking at the kinematic diagram shown in Figure 3, it is possible to notice that the implementation of the humanoid's hip is not spherical; that is, the first three axes do not intersect in one point. There is an offset represented by the link ${\mathrm{L}}_{\mathrm{s}2}$ which prevents the usage of the kinematic decoupling method proposed by Pieper [3] to solve the IK problem. In this way, alternative geometric analytical solutions were also sought without success. However, a purely numerical solution is not feasible either due to the existence of a variable efficiency on its implementation.

For this reason, a hybrid approach is proposed. It consists on first finding an approximate analytical “solution” for the inverse kinematics problem by neglecting the presence of the offset at the hip, as it is small when compared to the actual size of the robot, and then to use this “solution” as the initial condition of a numerical method, so that the number of iterations needed be drastically reduced.

This idea is similar to the one proposed in [6]. However, the last one only deals with the case of a double universal joint placed at the wrist. On the other hand, the following procedure is not restricted to any particular kinematic configuration, as long as the offset is small compared to the size of the robot.

4.1. Approximate analytical “solution”

First, suppose that the three axes at the hip intersect in just one point, represented by the actual intersection of the first two axes. Likewise, suppose that the length of the third link is not ${\mathrm{l}}_3,$ but ${\mathrm{l}}_2+{\mathrm{l}}_3,$ in order to diminish the error.

By doing this, a purely analytical “solution” exists, so that it can be found by using the so-called Paden-Kahan methodology [24], which is based on the product of exponentials formula shown in (9) [1, 25].

The equation to be solved is

${\bar{g}}_{\mathrm{d},sf}\left(\theta \right)={\mathrm{e}}^{{\stackrel{⌢}{\xi }}_{\mathrm{s}1}{\theta }_{\mathrm{s}1}}{\mathrm{e}}^{{\stackrel{⌢}{\xi }}_{\mathrm{s}2}{\theta }_{\mathrm{s}2}}\dots {\mathrm{e}}^{{\stackrel{⌢}{\xi }}_{\mathrm{s}6}{\theta }_{\mathrm{s}6}}{\bar{g}}_{\mathrm{d},sf}\left(0\right):={\bar{g}}_{\mathrm{d}}$ (26) ${\mathrm{e}}^{{\stackrel{⌢}{\xi }}_{\mathrm{s}1}{\theta }_{\mathrm{s}1}}{\mathrm{e}}^{{\stackrel{⌢}{\xi }}_{\mathrm{s}2}{\theta }_{\mathrm{s}2}}\dots {\mathrm{e}}^{{\stackrel{⌢}{\xi }}_{\mathrm{s}6}{\theta }_{\mathrm{s}6}}={\bar{g}}_{\mathrm{d}}{\bar{g}}_{\mathrm{d},sf}^{-1}\left(0\right):={\bar{g}}_1$ (27)

To do that, let us propose the following four points based on the reference attitude of the robot shown in Figure 3:

• A point $p_H$ located at the position of the supposed spherical hip; that is, at the intersection of the axes of the first two joints:

$p_H:={\left[\begin{array}{ccc}0& -{\mathrm{l}}_1& {\mathrm{l}}_0\end{array}\right]}^{\mathrm{T}}$ (28)

• A point $p_A$ located at the position of the ankle; that is, at the intersection of the axes of the two last joints:

$p_A:={\left[\begin{array}{ccc}0& -\left({\mathrm{l}}_1+{\mathrm{l}}_2+{\mathrm{l}}_3+{\mathrm{l}}_4\right)& {\mathrm{l}}_0\end{array}\right]}^{\mathrm{T}}$ (29)

• A point $p_L$ located on the axis of the first joint:

$p_L:={\left[\begin{array}{ccc}0& 0& {\mathrm{l}}_0\end{array}\right]}^{\mathrm{T}}$ (30)

• A point $p_D$ located at the origin of the body frame:

$p_D:={\left[\begin{array}{ccc}0& 0& 0\end{array}\right]}^{\mathrm{T}}$ (31)

Having done this, the joint variables can be calculated as it is explained in the following [24]:

Calculation of ${\theta }_{s4}$ Let us apply both sides of (27) to the point $p_A$ and then, find the Euclidian distance between the transformation's result and the point $p_H.$ Notice that the point $p_A$ is located on the axes of the joints ${\mathrm{J}}_{\mathrm{s}5}$ and ${\mathrm{J}}_{\mathrm{s}6},$ so that the corresponding operators have no effect on this point and can be omitted. That is,

$\begin{array}{l}\Vert {\mathrm{e}}^{{\stackrel{⌢}{\xi }}_{\mathrm{s}1}{\theta }_{\mathrm{s}1}}\dots {\mathrm{e}}^{{\stackrel{⌢}{\xi }}_{\mathrm{s}5}{\theta }_{\mathrm{s}5}}{\mathrm{e}}^{{\stackrel{⌢}{\xi }}_{\mathrm{s}6}{\theta }_{\mathrm{s}6}}{p}_A-p_H\Vert =\Vert {\bar{g}}_1{p}_A-p_H\Vert \\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\Vert {\mathrm{e}}^{{\stackrel{⌢}{\xi }}_{\mathrm{s}1}{\theta }_{\mathrm{s}1}}\dots {\mathrm{e}}^{{\stackrel{⌢}{\xi }}_{\mathrm{s}4}{\theta }_{\mathrm{s}4}}{p}_A-p_H\Vert =\delta \end{array}$ (32)

In the same way, given that we have a similar situation at the hip

$p_H={\mathrm{e}}^{{\stackrel{⌢}{\xi }}_{\mathrm{s}1}{\theta }_{\mathrm{s}1}}{\mathrm{e}}^{{\stackrel{⌢}{\xi }}_{\mathrm{s}2}{\theta }_{\mathrm{s}2}}{\mathrm{e}}^{{\stackrel{⌢}{\xi }}_{\mathrm{s}3}{\theta }_{\mathrm{s}3}}{p}_H$ (33)

Equation (32) can be rewritten as

$\Vert {\mathrm{e}}^{{\stackrel{⌢}{\xi }}_{\mathrm{s}1}{\theta }_{\mathrm{s}1}}{\mathrm{e}}^{{\stackrel{⌢}{\xi }}_{\mathrm{s}2}{\theta }_{\mathrm{s}2}}{\mathrm{e}}^{{\stackrel{⌢}{\xi }}_{\mathrm{s}3}{\theta }_{\mathrm{s}3}}\left({\mathrm{e}}^{{\stackrel{⌢}{\xi }}_{\mathrm{s}4}{\theta }_{\mathrm{s}4}}{p}_A-p_H\right)\Vert =\delta$ (34)

Now then, recalling that a rigid body transformation preserves distances between points, the joint variable ${\theta }_{\mathrm{s}4}$ can be found by recognizing the equivalent expression as the third Paden-Kahan sub problem (P-K-3) [24], whose solution is already known:

$\Vert {\mathrm{e}}^{{\stackrel{⌢}{\xi }}_{\mathrm{s}4}{\theta }_{\mathrm{s}4}}{p}_A-p_H\Vert =\delta \hspace{0.17em}\underrightarrow{\mathrm{P}-\mathrm{K}-3\hspace{0.17em}}\hspace{0.17em}\hspace{0.17em}{\theta }_{\mathrm{s}4}$ (35)

Calculation of ${\theta }_{s5}$ and ${\theta }_{s6}$ Let us reorder (27) and apply both sides to the point $p_H,$ considering the fact shown in (33). In this way,

$\begin{array}{l}{\mathrm{e}}^{{\stackrel{⌢}{\xi }}_{\mathrm{s}5}{\theta }_{\mathrm{s}5}}{\mathrm{e}}^{{\stackrel{⌢}{\xi }}_{\mathrm{s}6}{\theta }_{\mathrm{s}6}}{g}_1^{-1}{p}_H={\mathrm{e}}^{-{\stackrel{⌢}{\xi }}_{\mathrm{s}4}{\theta }_{\mathrm{s}4}}{\mathrm{e}}^{-{\stackrel{⌢}{\xi }}_{\mathrm{s}3}{\theta }_{\mathrm{s}3}}\dots {\mathrm{e}}^{-{\stackrel{⌢}{\xi }}_{\mathrm{s}1}{\theta }_{\mathrm{s}1}}{p}_H\\ {\mathrm{e}}^{{\stackrel{⌢}{\xi }}_{\mathrm{s}5}{\theta }_{\mathrm{s}5}}{\mathrm{e}}^{{\stackrel{⌢}{\xi }}_{\mathrm{s}6}{\theta }_{\mathrm{s}6}}{g}_1^{-1}{p}_H={\mathrm{e}}^{-{\stackrel{⌢}{\xi }}_{\mathrm{s}4}{\theta }_{\mathrm{s}4}}{p}_H\end{array}$ (36)

Now then, as ${\theta }_{\mathrm{s}4}$ is already known, it is possible to define the following points:

$p:={\mathrm{g}}_1^{-1}{p}_H$ (37) $q:={\mathrm{e}}^{-{\stackrel{⌢}{\xi }}_{\mathrm{s}4}{\theta }_{\mathrm{s}4}}{p}_H$ (38)

So that, both ${\theta }_{\mathrm{s}5}$ and ${\theta }_{\mathrm{s}6}$ can be found by recognizing the expression in (36) as the second Paden-Kahan sub problem (P-K-2), after the proper substitution of $p\text{\hspace{0.17em}and\hspace{0.17em}}q$ [24]; that is,

${\mathrm{e}}^{{\stackrel{⌢}{\xi }}_{\mathrm{s}5}{\theta }_{\mathrm{s}5}}{\mathrm{e}}^{{\stackrel{⌢}{\xi }}_{\mathrm{s}6}{\theta }_{\mathrm{s}6}}p=q\hspace{0.17em}\underrightarrow{\mathrm{P}-\mathrm{K}-2\hspace{0.17em}}\hspace{0.17em}\hspace{0.17em}\left\{{\theta }_{\mathrm{s}5},{\theta }_{\mathrm{s}6}\right\}$ (39)

Calculation of ${\theta }_{s2}$ and ${\theta }_{s3}$ Let us reorder again (27) and apply both sides to the point $p_L.$

Remember that this point is located on the axis of ${\mathrm{J}}_{\mathrm{s}1},$ so that the joint has no effect on this point. In this way,

$\begin{array}{l}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\mathrm{e}}^{{\stackrel{⌢}{\xi }}_{\mathrm{s}2}{\theta }_{\mathrm{s}2}}{\mathrm{e}}^{{\stackrel{⌢}{\xi }}_{\mathrm{s}3}{\theta }_{\mathrm{s}3}}{\mathrm{e}}^{{\stackrel{⌢}{\xi }}_{\mathrm{s}4}{\theta }_{\mathrm{s}4}}\dots {\mathrm{e}}^{{\stackrel{⌢}{\xi }}_{\mathrm{s}6}{\theta }_{\mathrm{s}6}}{g}_1^{-1}{p}_L={\mathrm{e}}^{-{\stackrel{⌢}{\xi }}_{\mathrm{s}1}{\theta }_{\mathrm{s}1}}{p}_L\\ {\mathrm{e}}^{{\stackrel{⌢}{\xi }}_{\mathrm{s}2}{\theta }_{\mathrm{s}2}}{\mathrm{e}}^{{\stackrel{⌢}{\xi }}_{\mathrm{s}3}{\theta }_{\mathrm{s}3}}\left({\mathrm{e}}^{{\stackrel{⌢}{\xi }}_{\mathrm{s}4}{\theta }_{\mathrm{s}4}}\dots {\mathrm{e}}^{{\stackrel{⌢}{\xi }}_{\mathrm{s}6}{\theta }_{\mathrm{s}6}}{g}_1^{-1}{p}_L\right)=p_L\end{array}$ (40)

Now then, given that ${\theta }_{\mathrm{s}4},{\theta }_{\mathrm{s}5},$ and ${\theta }_{\mathrm{s}6}$ are known, it is possible to define

$p={\mathrm{e}}^{{\stackrel{⌢}{\xi }}_{\mathrm{s}4}{\theta }_{\mathrm{s}4}}{\mathrm{e}}^{{\stackrel{⌢}{\xi }}_{\mathrm{s}5}{\theta }_{\mathrm{s}5}}{\mathrm{e}}^{{\stackrel{⌢}{\xi }}_{\mathrm{s}6}{\theta }_{\mathrm{s}6}}{g}_1^{-1}{p}_L$ (41)

In such a way that ${\theta }_{\mathrm{s}2}\text{\hspace{0.17em}y\hspace{0.17em}\hspace{0.17em}}{\theta }_{\mathrm{s}3}$ can be found by recognizing the expression in (40) as the second Paden-Kahan sub problem (P-K-2) again, after the proper substitution of p [24]; that is,

${\mathrm{e}}^{{\stackrel{⌢}{\xi }}_{\mathrm{s}2}{\theta }_{\mathrm{s}2}}{\mathrm{e}}^{{\stackrel{⌢}{\xi }}_{\mathrm{s}3}{\theta }_{\mathrm{s}3}}p=p_L\underrightarrow{\mathrm{P}-\mathrm{K}-2\hspace{0.17em}\hspace{0.17em}}\hspace{0.17em}\left\{{\theta }_{\mathrm{s}2},{\theta }_{\mathrm{s}3}\right\}$ (42)

Calculation of ${\theta }_{s1}$ Reorder again (27) and apply both sides to the point $p_D.$ In this way,

${\mathrm{e}}^{{\stackrel{⌢}{\xi }}_{\mathrm{s}1}{\theta }_{\mathrm{s}1}}{p}_D={\mathrm{g}}_1{\mathrm{e}}^{-{\stackrel{⌢}{\xi }}_{\mathrm{s}6}{\theta }_{\mathrm{s}6}}\dots {\mathrm{e}}^{-{\stackrel{⌢}{\xi }}_{\mathrm{s}2}{\theta }_{\mathrm{s}2}}{p}_D$ (43)

Now then, as all the other angles are already known, it is possible to define

$q={\mathrm{g}}_1{\mathrm{e}}^{-{\stackrel{⌢}{\xi }}_{\mathrm{s}6}{\theta }_{\mathrm{s}6}}\dots {\mathrm{e}}^{-{\stackrel{⌢}{\xi }}_{\mathrm{s}2}{\theta }_{\mathrm{s}2}}{p}_D$ (44)

In such a way that ${\theta }_{\mathrm{s}1}$ can be finally found by recognizing the expression in (43) as the first Paden-Kahan sub problem (P-K-1), after the substitution of q [24]; that is,

${\mathrm{e}}^{{\stackrel{⌢}{\xi }}_{\mathrm{s}1}{\theta }_{\mathrm{s}1}}{p}_D=q\underrightarrow{\mathrm{P}-\mathrm{K}-1\hspace{0.17em}}\hspace{0.17em}\hspace{0.17em}{\theta }_{\mathrm{s}1}$ (45)

It is worth to remark that this procedure may lead to a maximum of eight different joint configurations, due to the multiple solutions of P-K-2 and P-K-3 for equations (35), (39), and (42) (each one has two different solutions). However, it is possible to get rid of four among these ones, as it is natural to pre-select the knee-front configuration, and then to notice that the remaining four are not all valid ones, as not every constraint is considered in every calculation (the procedure required to get rid of information on constraints in order to apply the Paden-Kahan sub-problems).

By choosing the knee-front configuration only two solutions are valid, which can be selected by calculating the four candidate joint configurations and measuring the attitude's error generated by each one. Ideally, the two valid ones must lead to a zero attitude error; but, as this method is based on an approximation, only a minimum error may be achieved. In this way, the joint configuration which lead to lowest error is then selected as the “solution” of the problem.

4.2. Numerical refining

The approximate “solution” obtained by supposing a spherical hip generates significant errors that can not be neglected when dealing with gait generation. However, these ones are small compared to the actual size of the robot. This fact can be depicted in Figure 6, where both the real kinematic structure and the supposed one are shown having the same joint values.

Figure 6.

Real kinematic structure (lighter one) vs. supposed one (darker one)

As can be seen, for ${\theta }_{\mathrm{s}3}=16.35°$ (which is actually the joint limit), there is an error of 2.787 cm with respect to the end-effector's position, whereas the robot's size is approximately 1 m. In this way, if such approximate “solution” is used as the initial condition of a numerical method, the algorithm will likely converge to the desired solution instead of getting trapped in a local minima. Furthermore, the number of iterations needed for accomplishing a negligible error can be drastically reduced in most of the cases. These statements rely in the fact that the convergence of any numerical method almost entirely depends on the selection of the initial values [20].

Numerical algorithms

A numerical algorithm for this problem basically increments (or decrements) the joint configuration θ by means of an increment vector $\mathrm{\Delta }\theta ,$ in a certain number of iterations; that is, $\theta :=\theta +\mathrm{\Delta }\theta$ in each one. To compute this $\mathrm{\Delta }\theta$ there are many available methods which primarily use the Jacobian $\mathrm{J},$ by taking advantage of the fact that [19]

$\mathrm{\Delta }x\approx \mathrm{J}\mathrm{\Delta }\theta$ (46)

Where $\mathrm{\Delta }x$ is an increment of the end-effector's attitude.

One way of doing this is using the Moore-Penrose's pseudoinverse of the Jacobian ${\mathrm{J}}^{\text{†}},$ such that

$\mathrm{\Delta }\theta ={\mathrm{J}}^{\text{†}}e$ (47)

Where e is the error vector between the desired attitude and the actual one.

The pseudoinverse provides the “best fit” solution to a system of linear equations that lacks of unique solution. However, it is prone to be unstable in the neighborhood of singularities, as the singular values of ${\mathrm{J}}^{\text{†}}$ are calculated from the inverse of the singular values of $\mathrm{J},$ which are small in the neighborhood of singularities [19].

The Levenberg-Marquardt's method (also known as the damped least squares' method) improves the pseudoinvere's behavior near the singularities, resulting in a stable way for computing $\mathrm{\Delta }\theta .$ This is done by including a damping factor $\lambda ,$ such that [19]

$\mathrm{\Delta }\theta ={\left({\mathrm{J}}^{\mathrm{T}}\mathrm{J}+{\lambda }^2\mathrm{I}\right)}^{-1}{\mathrm{J}}^{\mathrm{T}}\mathrm{e}$ (48)

The reason for such improvement is because the singular values ${\mathrm{f}}_{\mathrm{i}}$ of the matrix ${\left({\mathrm{J}}^{\mathrm{T}}\mathrm{J}+{\lambda }^2\mathrm{I}\right)}^{-1}{\mathrm{J}}^{\mathrm{T}}$ can also be expressed as

${\mathrm{f}}_{\mathrm{i}}=\frac{{\sigma }_{\mathrm{i}}}{{\sigma }_{\mathrm{i}}^2+{\lambda }^2}$ (49)

Where ${\sigma }_{\mathrm{i}},\mathrm{i}=1,\dots ,\rho \left\{\mathrm{J}\right\}$ are singular values of $\mathrm{J}.$

So, given that λ is small, if ${\sigma }_{\mathrm{i}}\gg \lambda ,$ then $\frac{{\sigma }_{\mathrm{i}}}{{\sigma }_{\mathrm{i}}^2+{\lambda }^2}\approx \frac{1}{{\sigma }_{\mathrm{i}}},$ such that both methods behave in a similar way. But, if ${\sigma }_{\mathrm{i}}\approx \lambda ,$ then $\forall \lambda \hspace{0.17em}\text{>}\hspace{0.17em}0,\frac{{\sigma }_{\mathrm{i}}}{{\sigma }_{\mathrm{i}}^2+{\lambda }^2}\to 0$ as ${\sigma }_{\mathrm{i}}\to 0,$ behaving in a stable manner [19]. This fact gives us the clue to calculate λ as

$\lambda =\min \left\{{\sigma }_{\mathrm{i}},{\lambda }_{\max }\right\}$ (50)

Attitude's error computation

In order to apply the Levenberg-Marquardt's method it is first necessary to define the error vector $e\in {ℝ}^6,$ which is obtained by a concatenation of the positioning error $e_t\in {ℝ}^3$ and the orientation error $e_R\in {ℝ}^3;$ that is

$e={\left[\begin{array}{cc}{e}_t{}^{\mathrm{T}}& e_R{}^{\mathrm{T}}\end{array}\right]}^{\mathrm{T}}$ (51)

Recall that the Jacobian relates differential increments in position and orientation of the end-effector to the differential increments of joint variables. Position's differential increments are effectively vectors with differential magnitude pointing in the movement's direction. In the other hand, orientation's ones are also vectors with differential magnitude, but pointing in the direction of the axis of rotation (unitary angular velocity vector).

Numerical algorithms deal with finite increments (errors), instead of differential ones (as an approximation). In this way, position and orientation errors are vectors with a finite magnitude, that point in the direction of the increment in the case of the former ones or in the direction of the axis of rotation in the case of latter ones.

So, whereas position error's calculation is straightfor-ward $e_t={\left[\mathrm{\Delta }\mathrm{X}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{\Delta }\mathrm{Y}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{\Delta }\mathrm{Z}\right]}^{\mathrm{T}},$ orientation error's one is not so obvious. This is because the chosen orientation coordinates are not the Cartesian components of a rotational increment vector.

In order to obtain the desired orientation error from the chosen orientation coordinates error, it is necessary first to obtain an expression that relates the rate of change of those coordinates $\left({\stackrel{.}{\theta }}_{\mathrm{s}},{\stackrel{.}{\theta }}_{\mathrm{f}},{\stackrel{.}{\theta }}_{\mathrm{t}}\right)$ with respect to the angular velocity vector Ω of the corresponding foot.

Consider the matrix in (22), the correspondent skew-symmetric matrix of $\mathrm{\Omega },$ written as $\stackrel{⌢}{\mathrm{\Omega }},$ can be calculated from [24]

$\stackrel{⌢}{\mathrm{\Omega }}={\stackrel{.}{R}}_{\text{foot}}{\mathrm{R}}_{\text{foot}}^{\mathrm{T}}$ (52)

From which the vector Ω can be extracted, so that

$\mathrm{\Omega }=\left[\begin{array}{ccc}0& {\mathrm{c}}_{\mathrm{s}}& -{\mathrm{s}}_{\mathrm{s}}{\mathrm{c}}_{\mathrm{f}}\\ 0& {\mathrm{s}}_{\mathrm{s}}& {\mathrm{c}}_{\mathrm{s}}{\mathrm{c}}_{\mathrm{f}}\\ 1& 0& {\mathrm{s}}_{\mathrm{f}}\end{array}\right]\left[\begin{array}{c}{\stackrel{.}{\theta }}_{\mathrm{s}}\\ {\stackrel{.}{\theta }}_{\mathrm{f}}\\ {\stackrel{.}{\theta }}_{\mathrm{t}}\end{array}\right]$ (53)

Having done this, the orientation error can be obtained by arbitrarily replacing the orientation coordinates rate of change with the corresponding error to build $e_R;$ that is,

$e_R=\left[\begin{array}{ccc}0& {\mathrm{c}}_{\mathrm{s}}& -{\mathrm{s}}_{\mathrm{s}}{\mathrm{c}}_{\mathrm{f}}\\ 0& {\mathrm{s}}_{\mathrm{s}}& {\mathrm{c}}_{\mathrm{s}}{\mathrm{c}}_{\mathrm{f}}\\ 1& 0& {\mathrm{s}}_{\mathrm{f}}\end{array}\right]\left[\begin{array}{c}\mathrm{\Delta }{\theta }_{\mathrm{s}}\\ \mathrm{\Delta }{\theta }_{\mathrm{f}}\\ \mathrm{\Delta }{\theta }_{\mathrm{t}}\end{array}\right]$ (54)

4.3. Workspace

Let us consider two types of workspace: the ideal one, ${\mathrm{W}}_{\mathrm{i}},$ defined by geometrical constraints; that is, by considering non-limited joints, and the real one, ${\mathrm{W}}_{\mathrm{r}},$ which takes into account the joint limits, so that ${\mathrm{W}}_{\mathrm{r}}\subset {\mathrm{W}}_{\mathrm{i}}.$

By having an analytical algorithm for the IK problem, it is possible to check a priori if a given attitude is inside the ideal workspace or not, in order to assure that a “solution” exists and avoid failures at the time of calculation.

In this case, an approximate workspace can be found by noticing that only the solution for P-K-3 considers a restricted-domain function, represented by the arccosine, arccos(.): $\left[-1,1\right]\to \left[0,\pi \right].$ So that if an unreachable attitude is specified, the argument of the arccosine function $\left(\mathrm{D}\right)$ will be a number outside its domain. In this way, if $-1\text{\hspace{0.17em}<\hspace{0.17em}}\mathrm{D}\text{\hspace{0.17em}<\hspace{0.17em}}1$ an approximate analytical “solution” exists.

5. Simulation results

In order to visually validate the algorithm proposed in this paper, a kinematic simulator was developed under the name of ARMS (Advanced Robot Motion Simulator) (Figure 7) [7].

Figure 7.

ARMS Kinematic Simulator

This simulator uses a VRML (Virtual Reality Modeling Language) model of the humanoid robot, which is driven by specifying its joint coordinates. These ones are directly specified by the user, or calculated from a desired attitude which is specified in task space.

For testing purposes a series of attitudes were specified for each foot, by considering the body as the base link. Some attitudes were proposed as lying inside the real workspace and others, outside, or as lying only inside the ideal workspace but not inside the real one, in order to test the algorithm's behavior. One example of each case is reported, along with some screenshots of the corresponding results.

5.1. Attitude inside the real workspace

Figure 8 shows the attitude reached by the model of the robot (frontal and sagittal views) after specifying that the right foot's sole should be located 75 cm below with respect to the pelvis and 20 cm just in front of it, while maintaining the reference attitude's distance between this foot and the sagittal plane of symmetry (17.3 cm). The right foot's orientation is also specified with a pitch angle of −20°. In the other hand, the location of the left foot's sole is specified to be 80 cm below with respect to the pelvis, 10 cm in front of it and 17.3 cm away from the sagittal plane of symmetry, while its orientation is specified with a roll angle of −10°, followed by a yaw angle of 30°.

Figure 8.

Resulting configuration after specifying an attitude inside the real workspace

As can be seen with the aid of superimposed labels and the screenshot of the specified attitude and the obtained one (Figure 9) the virtual robot not only reaches the desired goal, but it can be achieved with only 2 iterations.

Figure 9.

Control and feedback panels associated With Figure 8

5.2. Attitude inside the ideal workspace, but outside the real one defined by mechanical constraints

Figure 10 shows again an attitude reached by the model of the robot. In this case, the right foot's sole was specified to be located 75 cm below with respect to the pelvis and 20 cm just in front of it as before, but now it is also desired a distance of 30 cm between this foot and the sagittal plane of symmetry. In the other hand, the left foot's location is specified to be a the same height but 20 cm behind the pelvis, at a distance of 5 cm from the sagittal plane of symmetry. Both feet are expected to maintain the reference orientation; that is, 0° for all the orientation coordinates.

Figure 10.

Resulting configuration after specifying an attitude inside the ideal workspace, but outside the real one

However, the desired attitude cannot be reached as it can be confirmed by noticing that the right foot's sole is a little below of the desired one (marked with the aid of superimposed labels), and because the orientation of the left foot is not the reference one. This can also be confirmed in Figure 11, where the difference between the desired attitude and the obtained one is noticeable, due to the saturation of joint's values. The process needed 3 iterations in this case.

Figure 11.

Control and feedback panels associated with Figure 10

5.3. Attitude outside the ideal workspace

Consider now the case when a desired attitude is outside the ideal workspace, represented by the case when the feet are specified to reach a distance impossible to achieve, even at full extension of the legs. In this case the kinematic simulator performs no action at all, as indicated by the label “Ideally impossible attitude” (Figure 12).

Figure 12.

The attitude is not inside the ideal workspace

6. Evaluation

The performance of the proposed hybrid algorithm was evaluated with respect to the number of iterations needed to accomplish a desired precision by means of statistical comparison with a purely numerical solution, as well as with respect to the mean time required for performing the requested computations in both cases, as a way to measure the achieved improvement of performance.

This evaluation was carried out over different sets of data consisting on 1000 random 6-tuples of joint variables $\left\{{\theta }_{\mathrm{s}1},{\theta }_{\mathrm{s}2},\dots ,{\theta }_{\mathrm{s}6}\right\}$ each one. The leg was also chosen in a random manner.

Each 6-tuple, along with the leg's choice, was used to obtain the corresponding foot's attitude by means of applying the forward kinematics algorithm. The resulting attitude was then taken as the input of the IK algorithm. In this way, it was assured that this input would lead to a solution, which was already known for testing and comparison purposes.

When testing the proposed algorithm, the foot's attitude was used as the input of the approximate analytical algorithm in such a way that the approximate “solution” was taken as the initial condition of the Levenberg-Marquardt method.

When testing the purely numerical solution, the initial condition of the Levenberg-Marquardt method was set arbitrarily as ${\theta }_s=0.$

In both cases, an error vector magnitude less than 0.1 or a maximum of 1500 iterations were used as stop conditions. The maximum value for the parameter ${\lambda }_{\max }$ was also set to 0.01 (as this value led to the best performance through a set of tests). Then, the error between the desired attitude and the obtained one was computed and recorded, as well as the number of iterations required and the mean computation time.

This evaluation consisted on two different statistical populations described as follows:

1. The first one considered randomly generated attitudes inside the ideal workspace, but not necessarily inside the real one. These ones were obtained by constraining the random joint variables' values to −90^° and 90^°.

2. The second one considered randomly generated attitudes only inside the real workspace. In this case, the joint variables' values were constrained to the actual joint limits of the robot (Table 1).

Table 1.

Range of motion of each joint of SILO2 [23]

Joint		Range of motion
Hip rotation	LeftRight	−180° ≤ θl1 ≤ 15°−15° ≤ θr1 ≤ 180°
Hip flexion/extension	LeftRight	−70° ≤θl2 ≤ 70°−70° ≤ θr2 ≤ 70°
Hip abduction/adduction	LeftRight	73° ≤ θl3 ≤ 106°−106° ≤ θr3 ≤ 73°
Knee	LeftRight	−64° ≤ θl4 ≤−3°−64° ≤ θr4 ≤ −3°
Ankle flexion/extension	LeftRight	−15° ≤ θl5 ≤ 20°−15° ≤ θr5 ≤ 20°
Ankle abduction/adduction	LeftRight	−12° ≤ θl6 ≤ 12°−12° ≤ θr6 ≤ 12°

This evaluation tests were implemented in Matlab R2012a, and performed on a Laptop Dell XPS L702X, which has an Intel(R) Core(TM) i7-2670QM CPU @ 2.20 GHz, 8 GB in RAM, and a 64-bit Operating System (Windows 7 Utimate).

6.1. Random attitudes generated inside the ideal workspace, but not necessarily inside the real one

Figure 13a shows the histogram generated by the purely numerical approach, whereas Figure 13b shows the one generated by the proposed hybrid algorithm, obtained by constraining the random joint variables' values to −90° and 90°. In both cases, the histogram's range is set to 1 – 50 or more iterations (samples over 50 iterations are accumulated over 50 in the histogram).

Figure 13.

Computation efficiency comparison between (a) a purely numerical approach and (b) the proposed hybrid method, considering randomly-generated attitudes inside the ideal workspace

Table 2 also shows some numerical data associated with the evaluation; that is, the percentage of samples that succeeded to converge (within the specified precision) in 1500, 50, and 10 iterations, as well as the average error in position and orientation between the desired attitude and the obtained one (taking into account only the samples that succeeded to converge) and the mean time required per sample in both cases.

Table 2.

Performance comparison of both algorithms considering randomly-generated attitudes inside the ideal workspace

Indicator	Purely numerical	Hybrid solution
Convergence in 1500 it.	97.5%	99.5 %
Convergence in 50 it.	39.8%	92.6 %
Convergence in 10 it.	13.5%	91.4 %
Mean position error	0.0176 cm	0.0011 cm
Mean orientation error	0.0284°	0.0023°
Mean CPU time	71.663 ms	17.057 ms

As can be seen, the proposed algorithm shows an outstanding performance over the purely numerical one used in [19] or [21], as more than 91.4% of the samples are able to converge in 10 or less iterations, in contrast with the latter one which may need 50, 100, 1000 or more iterations to converge. In addition, the proposed algorithm is able to achieve a lower average error (with respect to position or orientation) while being almost 4 times faster.

6.2. Random attitudes generated inside the real workspace only

By taking into account the actual joint limits of the robot, both algorithms were compared once again. Figure 14a shows the histogram generated by the purely numerical approach, whereas Figure 14b shows the one generated by the proposed hybrid algorithm. The histogram's range is set again to 1 – 50 or more iterations. Table 2 shows the corresponding numerical data associated with this evaluation.

Figure 14.

Computation efficiency comparison between (a) a purely numerical approach and (b) the proposed hybrid method, considering randomly-generated attitudes inside the real workspace

Surprisingly, the proposed algorithm not only shows an outstanding performance over the purely numerical one, but a very good one feasible for implementation on real-time, as 99.5% of the samples converged in less than 10 iterations (in fact, 97.2% did it within 2 iterations). Notice that the average error is also lower than the previous case, and that the hybrid algorithm performs almost 5 times faster.

Table 3.

Performance comparison of both algorithms considering randomly-generated attitudes inside the real workspace

Indicator	Purely numerical	Hybrid solution
Convergence in 1500 it.	99.3%	99.9%
Convergence in 50 it.	61.8%	99.8 %
Convergence in 10 it.	45.0%	99.5 %
Mean position error	0.0025 cm	0.0006 cm
Mean orientation error	0.0042°	0.0002°
Mean CPU time	41.046 ms	7.762 ms

7. Conclusions

This paper proposed a semi-analytical algorithm for computing the IK of humanoid robots whose construction shows a small but non negligible offset at the hip (compared to its actual height) which prevents any purely analytical solution to be developed.

As a case of study, a robot with a non-typical kinematic structure, as it is the SILO2 prototype, is analyzed. However, the algorithm can be applied to any structure that possesses negligible offsets, so that an approximate analytical “solution” can be found after neglecting them, and then be used as the initial condition of a stable numerical method, ensuring in this way that the number of iterations needed will be drastically reduced (compared with the required ones if an arbitrary initial condition is proposed).

Besides that, as the approximate analytical “solution” is in the neighborhood of the desired target in the joint configuration space, it will be more likely that the numerical method reaches the desired solution, as the right minima is very near to the initial condition. So, this method will be as efficient as desired so long as the concerned offsets are small enough. This fact was also evaluated, as this algorithm needed only very few iterations:

Less than 10 in most of the cases, or a maximum of 2 if the desired attitudes belonged to a typical gait; that is, attitudes that do not require big joint ranges to be performed. However, even non-typical motions can be achieved in an affordable computation time.

Concerning to the SILO2 robot, this method effectively solved its IK problem, so that it was possible to calculate precise gaits in real-time, as they were required by the control laws that had been previously developed by the team work at the Institute of Industrial Automation. Furthermore, the kinematic simulator developed during this work represented a useful tool for testing desired attitudes on a virtual model before doing it on the real prototype, avoiding in this way failures that may cause severe damages to the robot if it falls down. The simulator was very useful to test attitudes in this way, as many of them cannot be realized due to mechanical constraints, even when they were feasible during a certain gait. By recognizing these attitudes, the gait could be modified until getting a realizable and useful one, or they may also be used to provide guidelines for modifying the SMART actuators, by redesigning the size of the crank and the connecting road in every four-link mechanism, in order to improve the design of the humanoid.