Using Lie Derivatives with Dual Quaternions for Parallel Robots

Abstract

We introduce the notion of the Lie derivative in the context of dual quaternions that represent rigid motions and twists. First we define the wrench in terms of dual quaternions. Then we show how the Lie derivative helps understand how actuators affect an end effector in parallel robots, and make it explicit in the two cases case of Stewart Platforms, and cable-driven parallel robots. We also show how to use Lie derivatives with the Newton-Raphson Method to solve the forward kinematic problem for over constrained parallel actuators. Finally, we derive the equations of motion of the end effector in dual quaternion form, which include the effect of inertia from the actuators.

1. Introduction

This paper is broadly about poses and/or rigid motions, and how they can be represented by unit dual quaternions. A pose is a description of a frame of reference with respect to a fixed frame of reference, and consists of an orientation, and a position. A rigid motion consists of a rotation followed by a translation. From a mathematical point of view, these can be considered to be the same thing, and so we use these terms interchangeably (but see [1] for a different point of view).

Along with poses and/or rigid motions, we have the notion of screws. First we have rates of change of poses/rigid bodies, which are angular and translational velocities, and are described by twists. Second we have descriptions of how to change the inertia of the rigid body, the wrench, that is, the torque and the force applied to the body. For more reading on rigid body kinematics and dynamics, including screw theory, we refer the reader to [2,3,4,5].

The notion of the dual quaternion, and its use to represent poses and rigid motions, seems to go back to McAulay [6], inspired by the earlier work by Clifford [7]. The notion of using dual quaternions to represent twists may be found in [8,9]. A basic introduction to dual quaternions may be found in [10,11], the latter also covering twists. Many authors have used dual quaternions to represent hierarchies of poses, that is, chains of manipulators [9,10,12,13,14]. Papers on representing kinematics or dynamics via dual quaternions include [8,15,16,17,18,19,20]. Dual quaternions have also found great use in computer graphics [21,22].

(The reader should be aware that [8,17] have incorrect formulas for the logarithm and exponential of dual quaternions—the correct formulas may be found in [23], and [24] for the exponential.)

The purpose of this paper is to introduce the notion of using Lie derivatives for dual quaternions. We show that these can be used to essentially automate the creation of rather complex formulas, which are required for forward kinematics, and for dynamic equations of motion.

The authors have successfully used these formulas, combine with the algorithm described in [25], to create software for controlling a cable-driven parallel robot, which was built by the Dynamic Systems Test Branch of the Software, Robotics, and Simulation Division (ER5) at the NASA Johnson Space Center. Because of Export Administration Regulations we are unable to provide many more details.

The paper is quite heavy with mathematical formulas. For this reason, the proofs of most of the statements, and many of the comments of a mathematical nature, relegated to the Appendix A.

2. Notation

Rotations can be represented by unit quaternions [26,27], which we briefly describe here. A quaternion is a quadruple of real numbers, written as $A=w+x\mathit{i}+y\mathit{j}+z\mathit{k}$, with the algebraic operations ${\mathit{i}}^2={\mathit{j}}^2={\mathit{k}}^2=\mathit{i}\mathit{j}\mathit{k}=-1$. Its conjugate is $A^{*}=w-x\mathit{i}-y\mathit{j}-y\mathit{k}$, its norm is $\left|A\right|={(w^2+x^2+y^2+z^2)}^{1/2}=\sqrt{A{A}^{*}}=\sqrt{A^{*}A}$, its normalization is $\widehat{A}=A/\left|A\right|$, its real part is $Re\left(A\right)=w=\frac{1}{2}(A+A^{*})$, and its imaginary part is $Im\left(A\right)=\mathit{i}x+\mathit{j}y+\mathit{k}z=\frac{1}{2}(A-A^{*})$. It is called a unit quaternion if $\left|A\right|=1$, a real quaternion if $Im\left(A\right)=0$, and a vector quaternion if $Re\left(A\right)=0$. If $A\ne 0$, the multiplicative inverse is given by $A^{-1}=A^{*}/{\left|A\right|}^2$. (Note that many sources use the word “pure” instead of “vector” in this context).

We identify three dimensional vectors with vector quaternions, by identifying $\mathit{i}$, $\mathit{j}$, and $\mathit{k}$ with the three standard unit vectors. A unit quaternion Q represents the rotation which takes the direction $\mathit{r}$ to $Q\mathit{r}{Q}^{*}$. A rotation by angle a about an axis $\mathit{s}$, where $\left|\mathit{s}\right|=1$, has two unit quaternion representations: $\pm (cos\left(\frac{1}{2}a\right)+\mathit{s}sin\left(\frac{1}{2}a\right))=\pm exp\left(\frac{1}{2}a\mathit{s}\right)$. Composition of rotations corresponds to multiplication of unit quaternions.

We can represent quaternions as four dimensional vectors, and give it the inner product

$$A·B=Re\left(A{B}^{*}\right)=Re\left(A^{*}B\right).$$

(1)

A dual quaternion is a pair of quaternions, written as $\eta =A+ϵB$, with the extra algebraic operation ${ϵ}^2=0$. We call $A=\mathcal{P}\left(\eta \right)$ the primary part of $\eta$, and $B=\mathcal{D}\left(\eta \right)$ the dual part of $\eta$.

The conjugate dual quaternion of $\eta =A+ϵB$ is ${\eta }^{*}=A^{*}+ϵB^{*}$. Conjugation reverses the order of multiplication:

$${\left({\eta }_1{\eta }_2\right)}^{*}={\eta }_2^{*}{\eta }_1^{*}.$$

(2)

There is another conjugation for dual quaternions: $\overline{A+ϵB}=A-ϵB$, but we have no cause to use it in this paper, except in Equation (7) below.

A unit dual quaternion $\eta =Q+ϵB$ is a dual quaternion such that ${\eta }^{*}\eta =1$, equivalently, that Q is a unit quaternion and $B·Q=0$. A vector dual quaternion $A+ϵB$ is a dual quaternion such that both A and B are vector quaternions.

If $\eta =A+ϵB$ is a dual quaternion with $A\ne 0$, then its multiplicative inverse can be calculated using the formula

$${\eta }^{-1}=A^{-1}-ϵA^{-1}B{A}^{-1}.$$

(3)

If $\eta$ is a unit dual quaternion, then there is a computationally much faster formula (see [8]):

$${\eta }^{-1}={\eta }^{*}.$$

(4)

We set $\mathbb{D}\mathbb{H}$ for the set of invertible dual quaternions (that is, $A+ϵB$ where $A\ne 0$), $\mathfrak{d}\mathfrak{h}$ for the set of dual quaternions, $\mathbb{S}\mathbb{D}\mathbb{H}$ for the set of unit dual quaternions, and $\mathfrak{s}\mathfrak{d}\mathfrak{h}$ for the set of vector dual quaternions.

A rigid motion of the form

$$\mathit{r}\mapsto Q\mathit{r}{Q}^{*}+\mathit{t},$$

(5)

where here $\mathit{t}$ is a translation, can be represented by the unit dual quaternion

$$\eta =Q+\frac{1}{2}ϵ\mathit{t}Q.$$

(6)

Composition of rigid motions corresponds to multiplication of unit dual quaternions, where the notation ${\eta }_1{\eta }_2$ means to apply first the rigid motion represented by ${\eta }_2$, and then by ${\eta }_1$, that is, the dual quaternion acts by left multiplication. If $\mathit{r}$ is a 3-vector, and $\mathit{s}$ is the image of $\mathit{r}$ under the action of the rigid motion $\eta =Q+ϵB$, then

$$1+ϵ\mathit{s}=\eta (1+ϵ\mathit{r}){\overline{\eta }}^{*},$$

(7)

but generally it is easier to use the formula

$$\mathit{s}=Q\mathit{r}{Q}^{*}+2B{Q}^{*}=(Q\mathit{r}+2B)Q^{*}.$$

(8)

For a dual quaternion, it is not really possible to mix its primary and dual parts additively. For example, for a unit dual quaternion that represents a rigid motion, the primary part is unitless, whereas the dual has units of length. For this reason, when measuring how large such a dual quaternion is, everything must be with respect to a characteristic length scale l. (For example, for a parallel robot, the characteristic length might be the width of the workspace of the end effector). The size of a dual quaternion is defined to be

$${\mathrm{size}}_l\left(\eta \right)={\left({\left|\mathcal{P}\left(\eta \right)\right|}^2+l^{-2}{\left|\mathcal{D}\left(B\right)\right|}^2\right)}^{1/2}.$$

(9)

A twist is the pair of vectors $(\mathit{w},\mathit{v})$ that describes the rate of change of pose or rigid motion, where $\mathit{w}$ is the angular velocity, and $\mathit{v}$ is the translational velocity. One can think of the twist $(\mathit{w},\mathit{v})$ as a rigid motion function of time t:

$$\mathit{r}\mapsto \mathit{r}+t(\mathit{w}\times \mathit{r}+\mathit{v})+O\left(t^2\right)\mathrm{as}t\to 0.$$

(10)

It has two possible meanings, depending upon whether the twist is understood to be with respect to the fixed frame, or with respect to the moving frame. If it is understood to be with respect to the moving frame, we have the formula

$$\frac{d}{dt}(Q\mathit{r}{Q}^{*}+\mathit{t})=Q(\mathit{w}\times \mathit{r}+\mathit{v})Q^{*}+\mathit{t},$$

(11)

and if it is understood to be with respect to the fixed frame

$$\frac{d}{dt}(Q\mathit{r}{Q}^{*}+\mathit{t})=\mathit{w}\times (Q\mathit{r}{Q}^{*}+\mathit{t})+\mathit{v}.$$

(12)

Then this twist can be represented by a vector dual quaternion [8,15]

$$\phi =\frac{1}{2}\mathit{w}+\frac{1}{2}ϵ\mathit{v},$$

(13)

where if we understand it to be with respect to the moving frame, we have the formula

$$\phi ={\eta }^{-1}\dot{\eta },\mathrm{or}\dot{\eta }=\eta \phi ,$$

(14)

and if we understand it to be with respect to the fixed frame

$$\phi =\dot{\eta }{\eta }^{-1},\mathrm{or}\dot{\eta }=\phi \eta .$$

(15)

In this paper, unless otherwise stated, we always understand the twist to be with respect to the moving frame.

We make the identification

$$\mathfrak{d}\mathfrak{h}\cong {\mathbb{R}}^8,$$

(16)

using the basis

$$({\beta }_1,{\beta }_2,{\beta }_3,{\beta }_4,{\beta }_5,{\beta }_6,{\beta }_7,{\beta }_8)=(\mathit{i},\mathit{j},\mathit{k},ϵ\mathit{i},ϵ\mathit{j},ϵ\mathit{k},1,ϵ),$$

(17)

and similarly, we make the identification

$$\mathfrak{s}\mathfrak{d}\mathfrak{h}\cong {\mathbb{R}}^6,$$

(18)

using the basis $({\beta }_1$, ${\beta }_2$, ${\beta }_3$, ${\beta }_4$, ${\beta }_5$, ${\beta }_6)$. With these identifications, we can define the dot product between two dual quaternions by transferring the usual definition of dot product on ${\mathbb{R}}^8$, that is

$$(A+ϵB)·(C+ϵD)=A·C+B·D.$$

(19)

In this way, every dual quaternion $\eta$ can be written in component form as

$$\eta =\sum _{i=1}^8{\eta }_i{\beta }_i,$$

(20)

and every vector dual quaternion $\theta$ as

$$\theta =\sum _{i=1}^6{\theta }_i{\beta }_i.$$

(21)

Finally, we give a few extra formulas. Let ${\phi }_m$ denote the twist with respect to the moving frame, and ${\phi }_f$ denote the twist with respect to the fixed frame. From Equations (14) and (15) we obtain

$${\phi }_f=\eta {\phi }_m{\eta }^{-1}.$$

(22)

Since in any algebra we have

$$\frac{d}{dt}{\eta }^{-1}=-{\eta }^{-1}\dot{\eta }{\eta }^{-1},$$

(23)

we obtain the surprisingly simple formula for the change of frame for acceleration:

$${\dot{\phi }}_f=\eta {\dot{\phi }}_m{\eta }^{-1}.$$

(24)

Note that this does not generalize to higher derivatives, for example, the formula for change of frame for jerk is

$${\ddot{\phi }}_f=\eta {\ddot{\phi }}_m{\eta }^{-1}+\eta (\phi \dot{\phi }-\dot{\phi }\phi ){\eta }^{-1}.$$

(25)

3. Dual Quaternions to Represent Wrenches

Let the pose $\eta$ represent the reference frame that moves with the end effector. It is not necessary (although it can simplify things) that the center of mass of the end effector coincides with the origin of the moving frame.

The wrench dual quaternion is defined to be

$$\tau =2\mathit{q}+2ϵ\mathit{p},$$

(26)

where $\mathit{q}$ and $\mathit{p}$ are the torque and force, respectively, applied to the end effector at the origin of the moving frame, measured with respect to the moving frame.

If ${\mathit{r}}_0$ is the center of mass of the end effector in the moving frame, then the twist about the center of mass is given by

$${\phi }_0=\phi +\frac{1}{2}ϵ\mathit{w}\times {\mathit{r}}_0,$$

(27)

where $\phi ={\eta }^{-1}\dot{\eta }$, and the wrench applied about the center of mass is

$${\tau }_0=\tau +2\mathit{p}\times {\mathit{r}}_0.$$

(28)

The reason for introducing the factor 2 in definition (26) is so that the rate of change of work done to the end effector is given by

$$\frac{d}{dt}(\mathrm{work}\mathrm{done})=\tau ·\phi ={\tau }_0·{\phi }_0.$$

(29)

(The second equality follows from vector identities).

See [2] for the origins of the term twist and wrench as pairs of 3-vectors, which are examples of screws. The ’work done’ formulas are also known as reciprocal screw relationships.

4. The Normalization of a Dual Quaternion

A dual number is anything of the form $a+ϵb$, where a and b are real numbers. The norm of a dual quaternion $\eta =A+ϵB$ is the dual number defined by the two steps:

$$\begin{array}{c}{\left|\eta \right|}^2={\eta }^{*}\eta =\eta {\eta }^{*}={\left|A\right|}^2+2ϵ(A·B),\end{array}$$

(30)

$$\begin{array}{c}\left|\eta \right|=\sqrt{{\left|\eta \right|}^2}=\left|A\right|+ϵ(A·B)/\left|A\right|.\end{array}$$

(31)

The norm preserves multiplication, that is, if ${\eta }_1$ and ${\eta }_2$ are two dual quaternions, then

$$|{\eta }_1{\eta }_2|=|{\eta }_1\left|\right|{\eta }_2|.$$

(32)

If $\eta =Q+ϵB$ is an invertible dual quaternion, then we define its normalization to be the unit dual quaternion

$$\widehat{\eta }={\left|\eta \right|}^{-1}\eta ={\eta \left|\eta \right|}^{-1}={Q/\left|Q\right|+ϵ(B/|Q|-(B·Q)Q/|Q|}^3).$$

(33)

(We remark that the normalization of an invertible dual quaternion is used in the computer graphics industry [21,22]). While this normalization formula might seem initially quite complicated, after thinking about it one can see that it is the simplest projection that enforces $\left|Q\right|=1$ and $B·Q=0$.

The normalization also satisfies the following properties.

• If $\eta$ is a unit dual quaternion, then $\widehat{\eta }=\eta$.
• Normalization preserves multiplication, that is, if ${\eta }_1$ and ${\eta }_2$ are two dual quaternions, then

  $$\widehat{{\eta }_1{\eta }_2}={\widehat{\eta }}_1{\widehat{\eta }}_2.$$

  (34)

5. Notation for Three by Three Matrices

Let $\mathsf{I}$ denote the $(3\times 3)$ identity matrix, and $\mathsf{0}$ denote the $(3\times 3)$ zero matrix. If $\mathit{r}$ is a 3-vector, then the Hodge star operator of $\mathit{r}$ is

$$\mathsf{S}\left(\mathit{r}\right)=\left[\begin{array}{ccc}0& -r_3& r_2\\ r_3& 0& -r_1\\ -r_2& r_1& 0\end{array}\right].$$

(35)

Note that

$$\mathsf{S}\left(\mathit{r}\right)\mathit{s}=\mathit{r}\times \mathit{s}.$$

(36)

If $\mathit{u}$ is a unit vector, then the projection onto the complement of the unit vector $\mathit{u}$ is defined by

$${\mathsf{P}}_{\mathit{u}}\mathit{x}=\mathit{x}-(\mathit{u}·\mathit{x})\mathit{u}.$$

(37)

Consistent with the identification (18), if $\mathsf{A}$, $\mathsf{B}$, $\mathsf{C}$, and $\mathsf{D}$ are $(3\times 3)$ matrices, and $\theta =\frac{1}{2}\mathit{a}+\frac{1}{2}ϵ\mathit{b}$, $\psi =\frac{1}{2}\mathit{c}+\frac{1}{2}ϵ\mathit{d}$ are vector dual quaternions, we have

$$\begin{array}{c}\left[\begin{array}{cc}\mathsf{A}& \mathsf{B}\end{array}\right]\theta =\frac{1}{2}\mathsf{A}\mathit{a}+\frac{1}{2}\mathsf{B}\mathit{b},\end{array}$$

(38)

$$\begin{array}{c}\left[\begin{array}{cc}\mathsf{A}& \mathsf{B}\\ \mathsf{C}& \mathsf{D}\end{array}\right]\theta =\frac{1}{2}\mathsf{A}\mathit{a}+\frac{1}{2}\mathsf{B}\mathit{b}+\frac{1}{2}ϵ\mathsf{C}\mathit{a}+\frac{1}{2}ϵ\mathsf{D}\mathit{b},\end{array}$$

(39)

$$\begin{array}{c}\theta ·\left[\begin{array}{cc}\mathsf{A}& \mathsf{B}\\ \mathsf{C}& \mathsf{D}\end{array}\right]\psi =\frac{1}{4}\mathit{a}·\mathsf{A}\mathit{c}+\frac{1}{4}\mathit{a}·\mathsf{B}\mathit{d}+\frac{1}{4}\mathit{b}·\mathsf{C}\mathit{c}+\frac{1}{4}\mathit{b}·\mathsf{D}\mathit{d}.\end{array}$$

(40)

6. Lie Derivatives

The notion of the Lie derivative, sometimes in our context called the directional derivative, is a combination of two ideas that may be found in the literature. First is the concept of a Lie derivative with respect to a vector field [28]. Secondly, the definition of the Lie algebra is that it is the vector space of vector fields that are invariant under left multiplication by elements of the Lie group [29]. In this way, we can define the Lie derivative of a function with respect to an element of the Lie algebra. One place in the literature where they are combined is in [30] (Equation (5), Chapter II).

These standard abstract definitions can be made more concrete in our special case where the Lie group is the set of unit dual quaternions, and the Lie algebra is the set of vector dual quaternions.

Suppose one has a quantity that is a function of pose $g\left(\eta \right)$, and whose range is any vector space. (But for intuition, consider the case when the range is the real numbers, and think of the vector valued case as the formulas below simply being applied component wise).

Then we usually think of the derivative of $g\left(\eta \right)$ as the Jacobian with respect to the components of $\eta$. But it really makes more sense to compute the derivative with respect to the components of the perturbation of $\eta$. The latter is the Lie derivative.

The definition is this. Given a differentiable function g whose domain is the unit dual quaternions, $\mathbb{S}\mathbb{D}\mathbb{H}$, and whose codomain is any vector space over the real numbers, we can extend it arbitrarily to a differentiable function whose domain is an open neighborhood of $\mathbb{S}\mathbb{D}\mathbb{H}$ in $\mathfrak{d}\mathfrak{h}$. Given a unit dual quaternion $\eta$ and a vector dual quaternion $\theta$, we define the Lie derivative of $g\left(\eta \right)$ in the direction of $\theta$ to be

$${\mathcal{L}}_{\theta }g=\underset{r\to 0}{lim}\frac{g\left(\eta \right(1+r\theta \left)\right)-g\left(\eta \right)}{r}={\left.\frac{d}{dr}g\left(\eta (1+r\theta )\right)\right|}_{r=0}.$$

(41)

Since $\eta (1+r\theta )$ isn’t necessarily a unit dual quaternion, it is not obvious that the definition of the directional Lie derivative doesn’t depend upon how the domain of g was extended from $\mathbb{S}\mathbb{D}\mathbb{H}$, but it is, as is shown in Lemma A2 below.

Given a generic function g whose domain is the invertible dual quaternions, $\mathbb{D}\mathbb{H}$, and whose codomain is any vector space over the real numbers, we define its Jacobian to be the dual quaternion

$$\frac{\partial g}{\partial \eta }=\sum _{i=1}^8\left[\frac{\partial }{\partial {\eta }_i}g\left(\eta \right)\right]{\beta }_i,$$

(42)

where the partial derivative ${\displaystyle \frac{\partial }{\partial {\eta }_i}}$ is interpreted using the identification of ${\eta }_i$ with the ${\beta }_i$ components of $\eta$ as described in (17).

If the domain of g is the vector dual quaternions, $\mathfrak{s}\mathfrak{d}\mathfrak{h}$, we have Equation (42) except with 8 replaced by 6.

We have the following formula, which is useful for explicitly calculating the Lie derivative if an extension of g to an open neighborhood of $\mathbb{S}\mathbb{D}\mathbb{H}$ in $\mathfrak{d}\mathfrak{h}$ is known:

$${\mathcal{L}}_{\theta }g=\frac{\partial g}{\partial \eta }·\left(\eta \theta \right):=\sum _{i=1}^8\left(\frac{\partial g}{\partial {\eta }_i}\right){\left(\eta \theta \right)}_i,$$

(43)

where the notation ${\left(\eta \theta \right)}_i$ means the ${\beta }_i$ component of $\eta \theta$, as described in (17).

We define the partial Lie derivatives to be

$${\mathcal{L}}_{i}g={\mathcal{L}}_{{\beta }_i}g,(1\le i\le 6),$$

(44)

and its full Lie derivative to be the vector dual quaternion (or if the range of g is a vector space, the tensor product of a vector dual quaternion with a vector)

$$\mathcal{L}g=\sum _{i=1}^6{\mathcal{L}}_{i}g{\beta }_i,$$

(45)

so that for all vector dual quaternions $\theta$ it satisfies:

$$\theta ·\mathcal{L}g={\mathcal{L}}_{\theta }g.$$

(46)

To gain some intuition, write

$$\theta =\frac{1}{2}\mathit{a}+\frac{1}{2}ϵ\mathit{b}.$$

(47)

Since we have that

$$\frac{\partial }{\partial \theta }=2\frac{\partial }{\partial \mathit{a}}+2ϵ\frac{\partial }{\partial \mathit{b}},$$

(48)

we see that $\theta$ represents a change in pose by an infinitesimal translation $\mathit{b}$ and an infinitesimal rotation $\mathit{a}$, measured in the moving frame of reference. Thus $\mathcal{L}g$ is a vector dual quaternion giving twice the change in g with respect to an infinitesimal rotation, plus $ϵ$ times twice the change of g with respect to an infinitesimal translation.

One important property of the Lie derivative is that if $\eta$ represents a pose, with twist $\phi$, then by Equation (41), with r replaced by t, we see that

$$\frac{d}{dt}\left[g\left(\eta \right)\right]={\mathcal{L}}_{\phi }g.$$

(49)

The Lie derivative satisfies various rules, which easily follow from either Equation (41) or (43), which are also useful for explicitly calculating the Lie derivative when g is known.

• If $g\left(\eta \right)$ is linear in $\eta$, then

  $${\mathcal{L}}_{\theta }g\left(\eta \right)=g\left(\eta \theta \right).$$

  (50)
• The product rule: if * is any binary operator which is bilinear over the real numbers, such as the product of real numbers, the inner product, the cross product, or the dual quaternion product, then

  $${\mathcal{L}}_{\theta }(g_1*g_2)=g_1*\left({\mathcal{L}}_{\theta }{g}_2\right)+\left({\mathcal{L}}_{\theta }{g}_1\right)*g_2.$$

  (51)
• The chain rule:

  $${\mathcal{L}}_{\theta }\left(h(g_1,g_2,\cdots ,g_m)\right)=\sum _{i=1}^m\frac{\partial }{\partial g_i}h(g_1,g_2,\cdots ,g_m){\mathcal{L}}_{\theta }{g}_i.$$

  (52)
• Let $\tilde{\mathit{s}}$ be a constant position vector, and $\tilde{\mathit{n}}$ be a constant direction. Let $\mathit{s}$ and $\mathit{n}$ be their corresponding values with respect to the moving frame. Then

  $$\begin{array}{cc}\hfill {\mathcal{L}}_{\theta }\mathit{s}& =2\left[\begin{array}{cc}\mathsf{S}\left(\mathit{s}\right)& -\mathsf{I}\end{array}\right]\theta ,\hfill \end{array}$$

  (53)

  $$\begin{array}{cc}\hfill {\mathcal{L}}_{\theta }\mathit{n}& =2\left[\begin{array}{cc}\mathsf{S}\left(\mathit{n}\right)& \mathsf{0}\end{array}\right]\theta .\hfill \end{array}$$

  (54)

To simplify the writing of application software, Using these rules, one can create a software library in C++ that performs automatic Lie differentiation. Since the domain, and hence range, of the Lie derivative can be any vector space, the sensible way to do this is using templates to allow for a variety of different data types. Note also that the product rule (51) has to be implemented for every product that is used, and similarly the chain rule (52) for every function h that is used.

7. Applications to Parallel Robots

To aid our description, we show the example of a cable-driven parallel robot in Figure 1 and Figure 2. In Figure 1, we show the end effector in green, controlled by eight cables shown in red, which go around pulleys shown in green. Each cable, and how it is attached to the pulley, is shown in Figure 2.

Suppose that the position of the end effector of a parallel robot is given by n actuators, described by quantities

$$\ell ={\left({\ell }_j\right)}_{1\le j\le n}.$$

(55)

For example, for a cable-driven parallel robot, these represent the lengths of the cables, and typically $n=8$. The numbers ${\ell }_j$ only need to be determined up to a fixed number. Thus, for example, in Figure 2, the number ${\ell }_j$ could represent the length of cable from the end effector attachment point, to the small cross marked at the top of the pulley.

For the Gough or Stewart Platforms [4], we have $n=6$, and ${\ell }_j$ is simply the distance between the end effector attachment point and the ground attachment point.

Let us also denote the force exerted by the actuators by

$$\mathit{f}={\left(f_j\right)}_{1\le j\le n},$$

(56)

defined so that the rate of change of work performed through the actuators is given by

$$\frac{d}{dt}(\mathrm{work}\mathrm{done})=\mathit{f}·\dot{\ell }.$$

(57)

Suppose we have a function $\mathsf{L}:\mathbb{S}\mathbb{D}\mathbb{H}\to {\mathbb{R}}^n$, which calculates the required actuator values, ℓ, from the pose $\eta$ of the end effector frame. This is the inverse kinematics function.

We also define the $(n\times 6)$ matrix $\mathrm{\Lambda }$ by

$$\mathrm{\Lambda }\theta ={\mathcal{L}}_{\theta }\mathsf{L}={\mathcal{L}}_{\theta }\ell ,$$

(58)

or more explicitly, by

$${\mathrm{\Lambda }}_{i,j}={\mathcal{L}}_j{\ell }_i.$$

(59)

From Equation (49) we obtain

$$\dot{\ell }=\mathrm{\Lambda }\phi .$$

(60)

There is also a ($6\times n$) matrix $\mathsf{T}$ that maps the actuator forces to the wrench dual quaternion:

$$\tau =\mathsf{T}\mathit{f}.$$

(61)

This can be computed by balancing the force and torque exerted upon the end-effector. But it can also be computed with the following important identity:

$$\mathsf{T}={\mathrm{\Lambda }}^T.$$

(62)

8. Second Lie Derivatives

If g is a function of dual quaternions with codomain any vector space over the real numbers, we define its Hessian to be the $(8\times 8)$ matrix

$$\left[\frac{{\partial }^{2}g}{\partial {\eta }^2}\right]={\left[\frac{{\partial }^{2}g}{\partial {\eta }_i\partial {\eta }_j}\right]}_{1\le i,j\le 8}.$$

(63)

Thus the expression $\left[\frac{{\partial }^{2}g}{\partial {\eta }^2}\right]\gamma$ should be interpreted as a matrix product with $\gamma$ treated as an eight dimensional column vector.

Second Lie derivatives will be used in Newton’s Method, as well as in our statements of the equations of motion, both described below.

We have that

$$\begin{array}{cc}\hfill {\mathcal{L}}_{\theta }{\mathcal{L}}_{\psi }g& =\left(\eta \psi \right)·\left[\frac{{\partial }^{2}g}{\partial {\eta }^2}\right]\left(\eta \theta \right)+\frac{\partial g}{\partial \eta }·\left(\eta \theta \psi \right)\hfill \\ \hfill & :=\sum _{i=1}^8\sum _{j=1}^8{\left(\eta \psi \right)}_i\left(\frac{{\partial }^{2}g}{\partial {\eta }_i\partial {\eta }_j}\right){\left(\eta \theta \right)}_j+\sum _{i=1}^8\left(\frac{\partial g}{\partial {\eta }_i}\right){\left(\eta \theta \psi \right)}_i.\hfill \end{array}$$

(64)

Another way one might try to define the second derivative is to use the formula ${\left.\frac{{\partial }^2}{\partial {\theta }^2}g\left(\eta (1+\theta )\right)\right|}_{\theta =0}$. Unfortunately, this definition doesn’t work, as it depends upon the choice of how to extend the domain of g to all dual quaternions. The obvious choice of extension is to use the normalization $\tilde{g}\left(\eta \right)=g\left(\widehat{\eta }\right)$. We have the following formula for the Hessian of $\tilde{g}$:

$${\left.\frac{{\partial }^2}{\partial {\theta }^2}\tilde{g}\left(\eta (1+\theta )\right)\right|}_{\theta =0}={\left[\frac{1}{2}({\mathcal{L}}_i{\mathcal{L}}_{j}g\left(\eta \right)+{\mathcal{L}}_j{\mathcal{L}}_{i}g\left(\eta \right))\right]}_{1\le i,j\le 6}.$$

(65)

9. The Examples of a Stewart Platform, and a Cable-Driven Parallel Robot

As an example, let us consider parallel robots such as cable-driven robots, or Stewart Platforms. In this case, the end effector has certain ’attachment points’ on it, ${\mathit{r}}_1$, ${\mathit{r}}_2,\cdots$, ${\mathit{r}}_n$ where the cables or legs attach, the cables or legs are attached at the other end to ${\mathit{s}}_k$, and unit vectors ${\mathit{u}}_1$, ${\mathit{u}}_2,\cdots$, ${\mathit{u}}_n$ which are directions the cables or legs come into the end effector, all of these being measured in the end effector’s frame of reference. Note that for $1\le k\le n$

$${\mathit{u}}_k=\frac{{\mathit{r}}_k-{\mathit{s}}_k}{|{\mathit{r}}_k-{\mathit{s}}_k|}.$$

(66)

Then

$$\mathcal{L}{\ell }_k=2{\mathit{r}}_k\times {\mathit{u}}_k+2ϵ{\mathit{u}}_k=2\left[\begin{array}{c}\mathsf{S}\left({\mathit{r}}_k\right)\\ \mathsf{I}\end{array}\right]{\mathit{u}}_k,$$

(67)

or representing vectors as columns, we have

$$\mathrm{\Lambda }=2\left[\begin{array}{cc}{({\mathit{r}}_1\times {\mathit{u}}_1)}^T& {\mathit{u}}_1^T\\ {({\mathit{r}}_2\times {\mathit{u}}_2)}^T& {\mathit{u}}_2^T\\ ⋮& ⋮\\ {({\mathit{r}}_n\times {\mathit{u}}_n)}^T& {\mathit{u}}_n^T\end{array}\right].$$

(68)

For calculating the second Lie derivative of ${\ell }_k$, we need only know the first Lie derivative of ${\mathit{u}}_k$. Note that in the simple case that the attachment point of the cable to the frame is fixed in the fixed frame (for example, as in a Stewart Platform), we have

$$\mathcal{L}{\mathit{u}}_k=\frac{1}{|{\mathit{r}}_k-{\mathit{s}}_k|}{\mathsf{P}}_{{\mathit{u}}_k}\left[\begin{array}{cc}-\mathsf{S}\left({\mathit{s}}_k\right)& \mathsf{I}\end{array}\right],$$

(69)

and therefore

$${\mathcal{L}}^2{\ell }_k=2\left[\begin{array}{c}\mathsf{S}\left({\mathit{r}}_k\right)\\ \mathsf{I}\end{array}\right]\frac{1}{|{\mathit{r}}_k-{\mathit{s}}_k|}{\mathsf{P}}_{{\mathit{u}}_k}\left[\begin{array}{cc}-\mathsf{S}\left({\mathit{s}}_k\right)& \mathsf{I}\end{array}\right].$$

(70)

Let us also describe a more complicated situation, which matches the cable-driven parallel robot the NASA Johnson Space Center described in the introduction. For simplicity of notation, we drop the subscripts k. See Figure 2 for reference. Let us suppose that the cable attaches via an assembly, which is free to rotate about an axis parallel to the unit vector $\mathit{n}$, and passing through the point $\mathit{x}$. Attached to the assembly, at a fixed distance w from the fixed point $\mathit{x}$, by a rod perpendicular to $\mathit{n}$, is the center of a pulley of radius p, which rotates in the plane containing the axis of rotation and the point on the end effector $\mathit{r}$. A cable passes over the pulley in the $\mathit{n}$ direction from the center of the pulley. All of the vectors are expressed in the moving frame of reference.

We work in the $(\mathit{m},\mathit{n})$ coordinate system in which the origin is $\mathit{r}$. Let $(h_1,h_2)$ be the coordinates of the center of the pulley. Then

$$\mathit{u}=u_1\mathit{m}+u_2\mathit{n},$$

(71)

where

$$(u_1,u_2)=-\frac{1}{(h_1^2+h_2^2)}\left[\begin{array}{cc}{h}_1& -h_2\\ h_2& h_1\end{array}\right]\left[\begin{array}{c}\sqrt{h_1^2+h_2^2-p^2}\\ p\end{array}\right].$$

(72)

The cable length, up to an additive constant, is

$$\ell =\sqrt{h_1^2+h_2^2-p^2}+p{tan}^{-1}\left(\frac{u_2}{u_1}\right).$$

(73)

To avoid singularities, the optimal way to compute the inverse tangent, at least with the configuration shown in Figure 2, is to use $\mathrm{atan}2(-u_2,-u_1)$, where the commonly available $\mathrm{atan}2(y,x)$ function solves for $\theta$ where $x=rcos\theta$, $y=rsin\theta$, $r>0$, $-\pi <\theta \le \pi$.

The various required quantities are

$$\begin{array}{c}\mathit{d}={\mathsf{P}}_{\mathit{n}}(\mathit{x}-\mathit{r}),\end{array}$$

(74)

$$\begin{array}{c}\mathit{m}=\widehat{\mathit{d}}=\frac{\mathit{d}}{\left|\mathit{d}\right|},\end{array}$$

(75)

$$\begin{array}{c}{h}_2=\mathit{n}·(\mathit{x}-\mathit{r}),\end{array}$$

(76)

$$\begin{array}{c}{h}_1=\mathit{m}·(\mathit{x}-\mathit{r})-w=\sqrt{{|\mathit{x}-\mathit{r}|}^2-h_2^2}-w.\end{array}$$

(77)

The Lie derivatives of these quantities can be calculated using the automatic Lie differentiation described at the end of Section 6. The only second Lie derivative required is that of ℓ, and since we have Equation (67), no automatic second Lie differentiation is required.

Singularity Analysis for Stewart Platforms

When operating a Stewart Platform, a singularity occurs when there are no viable cable forces that can create an arbitrary wrench, or equivalently, when there exists infinitesimal perturbations of the end effector pose that don’t require a leg length change. These singularities are often called bifurcations, because after a Stewart Platform encounters a singularity, the end effector is free to move in more than one direction. Encountering a singularity can cause great damage to the Stewart Platform.

From these considerations, it becomes clear that a singularity happens for a Stewart Platform if and only if $det\left(\mathrm{\Lambda }\right)=det\left(\mathsf{T}\right)=0$. This is in agreement with the results obtained by Gosselin and Angeles [31].

Note that when considering more complex parallel robots, that a singularities can happen in other situations as well (see, for example, [32]).

10. Forward Kinematics

In robotics, there are several ways to find the pose of the end effector. One method is to use an optical system, but this is not very accurate. Another is to use a proprioceptive sensor, where the pose is found by integrating the acceleration and angular velocity of the end effector, but this is subject to drift. It would be extremely helpful if the end effector could be computed from the numbers ${\ell }_k$. These numbers lengths can be found easily and with high sampling frequency using, for example, encoders.

Let the set of admissible actuator values, $\mathbb{L}\subset {\mathbb{R}}^n$, be the range of the function $\mathsf{L}$. Then the forward kinematics function is

$$\mathsf{Y}:\mathbb{L}\to \mathbb{S}\mathbb{D}\mathbb{H},$$

(78)

which is a left inverse to $\mathsf{L}$. Because of possible measurement errors, $\mathsf{Y}$ should produce decent answers even if the actuator values are merely close to $\mathbb{L}$.

The methods described here are essentially the Newton-Raphson Method, and are all iterative methods. Given a guess ${\eta }_k$, we create a new guess ${\eta }_{k+1}$:

$${\eta }_{k+1}={\eta }_k\widehat{(1+{\theta }_{k+1})},$$

(79)

and iterate until some criterion is met. We measure ${\mathrm{size}}_l({\eta }_{k+1}-{\eta }_k)$, and see when it is smaller than some preset value, like ${10}^{-16}$.

Here we merely describe the method. In Appendix A, we explain why they work.

We would also like to mention a different approach to using dual quaternions to solve forward kinematics problems in [33], although we don’t think it will cover the more complicated situation described in Section 9.

10.1. Forward Kinematics for Stewart Platforms

First we describe how to solve the exact-constrained problem, that is, when $n=6$. This would be the case for a Stewart Platfom. Typically this is solved by writing the pose using Euler angles, which provides a way to represent the pose using a 6-vector. However, in the opinion of the authors, this becomes a rather complicated set of equations, resulting in quite lengthy code.

Our algorithm is

$${\theta }_{k+1}=-\mathrm{\Lambda }{\left({\eta }_k\right)}^{-1}(\mathsf{L}\left({\eta }_k\right)-\ell ).$$

(80)

This is easy to code, certainly simpler than methods which use Euler angles.

10.2. Forward Kinematics for Over-Constrained Parallel Robots

Next we focus on the over-constrained problem, that is, when the number of actuators n is greater than 6.

This problem has been solved by many others, for example, [33,34]. But we feel that this is more easily solved using dual quaternions. For example, using the programming language C++, one can quickly build classes representing dual quaternions, and then these formulas can be applied without any real thought.

We compute

$${\delta }_k={\left(\sum _{m=1}^n{\mathrm{\Lambda }}_{m,i}({\mathsf{L}}_m\left({\eta }_k\right)-{\ell }_m)\right)}_{1\le i\le 6},$$

(81)

and

$${\mathsf{H}}_k={\left[\sum _{m=1}^n{\mathrm{\Lambda }}_{m,i}{\mathrm{\Lambda }}_{m,j}+\frac{1}{2}{\mathcal{L}}_i{\mathrm{\Lambda }}_{m,j}({\mathsf{L}}_m\left({\eta }_k\right))+\frac{1}{2}{\mathcal{L}}_j{\mathrm{\Lambda }}_{m,i}({\mathsf{L}}_m\left({\eta }_k\right))\right]}_{1\le i,j\le 6}.$$

(82)

The algorithm is

$${\theta }_{k+1}={\mathsf{H}}_k^{-1}{\delta }_k.$$

(83)

10.3. Results of Simulations for Forward Kinematics

The software used to check these algorithms described here is currently proprietary, and so we cannot give too many details. We hope to get permission to release the software at a later date, and make it available to everyone.

To test the simple Stewart Platform algorithm, we created 10,000 random poses. The orientation of each pose had an angle no greater than ${30}^{\circ }$ from the identity pose. For each pose, the cable lengths were calculated. The Newton-Raphson Method was then applied to the cable lengths, with a random initial guess pose.

Only one run out of 10,000 failed. The average number of required iterations was about 5.1. The run time for each forward kinematics calculation was a little under 100 microseconds, using a fairly modern but low-end laptop. Increasing the allowed angle to ${45}^{\circ }$ also gave a failure rate of one in 10,000.

The Newton-Raphson Method for the over-constrained robot given by the more complicated situation described in Section 9 was more delicate. In particular, since it is minimizing a loss function rather than directly solving the problem, it is possible that it might find local minima of the loss function which didn’t correspond to the solution.

For the first test, we created 1000 random poses, and computed their cable lengths. The Newton-Raphson Method was applied, with a random initial guess pose, and was allowed up to 50 iterations. If the loss function of the final answer was greater than ${10}^{-16}$, the Newton-Raphson Method was applied again.

All poses were found, but the average number of times the Newton-Raphson had to be applied was about 80. The average time to find a pose was about 5 ms.

For the second test, we again created 1000 random poses and computed their cable lengths. Then the Newton-Raphson was applied, with an initial guess that was 1% different from the original pose. Again success was measured by computing the loss function. But there were no second chances.

All poses were found, the average number of iterations required was 4.2, and the average time taken was a little under 100 microseconds.

Increasing the allowable difference between the initial guess and the original pose to 5% resulted in only about 89% of the runs being successful.

Note that in all these trials, since the poses were randomly created, it is quite likely that some of the poses were non-feasible for the parallel robots under consideration. (For example, it might be impossible to maintain that pose while keeping the cables under tension, or moving to that pose might require crossing cables).

When using this algorithm, to find the initial pose, we suggest to use the first method of trying 1000 different initial guesses, and choosing the pose with the lowest loss function. But thereafter, sample the cable lengths often, and use the previously measured pose as the initial guess. In our numerical simulations, the pose changes by about 0.3% per time step (about 1 ms), and the algorithm has never failed to solve the forward kinematics problem.

11. Dynamics of the End Effector

The dynamics equations of motion of rigid bodies is well known. But in this section, we also consider the additional effect of overcoming inertia from the actuators. Furthermore, we feel that it is nice to see this stated and derived in the context of dual quaternions.

Let us suppose that the kinetic energy of the parallel robot is given by

$$e=\frac{1}{2}\phi ·\mathsf{M}\phi ,$$

(84)

where $\mathsf{M}$ is a $(6\times 6)$ positive definite matrix, which depends only upon $\eta$, and which we call the effective mass of the parallel robot.

We define the no-load forces to be the actuator forces if there is no end effector present:

$${\mathit{f}}_0={\mathsf{M}}_0\ddot{\ell },$$

(85)

where ${\mathsf{M}}_0$ is a positive definite $(n\times n)$ matrix denoting what we shall call the effective no-load mass of the actuators. This might be caused by, for example, the reflected moment of inertia of the motor that drives each actuator, in which case ${\mathsf{M}}_0$ is simply a constant multiple of the identity. Since $\frac{1}{2}\dot{\ell }·{\mathsf{M}}_0\dot{\ell }=\frac{1}{2}\phi ·{\mathrm{\Lambda }}^T{\mathsf{M}}_0\mathrm{\Lambda }\phi$ is part of the kinetic energy, it follows that $\mathsf{M}-{\mathrm{\Lambda }}^T{\mathsf{M}}_0\mathrm{\Lambda }$ is a positive semi-definite matrix.

If $m_e$ is the mass of the end effector, ${\mathsf{M}}_e$ is the moment of inertia tensor of the end effector about its center of mass, and ${\mathit{r}}_0$ is the center of mass of the end effectors, all measured with respect to the moving frame, then

$$e=\frac{1}{2}{m}_e{|\mathit{v}+\mathit{w}\times {\mathit{r}}_0|}^2+\frac{1}{2}\mathit{w}·{\mathsf{M}}_e\mathit{w}+\frac{1}{2}\dot{\ell }·{\mathsf{M}}_0\dot{\ell },$$

(86)

that is

$$\mathsf{M}=4\left[\begin{array}{cc}{\mathsf{M}}_e-m_e\mathsf{S}{\left({\mathit{r}}_0\right)}^2& m_e\mathsf{S}\left({\mathit{r}}_0\right)\\ -m_e\mathsf{S}\left({\mathit{r}}_0\right)& m_e\mathsf{I}\end{array}\right]+{\mathrm{\Lambda }}^T{\mathsf{M}}_0\mathrm{\Lambda }.$$

(87)

Theorem 1. 

If the kinetic energy satisfies Equation (84) with Equation (87) holding, and the potential energy v is calculated in the usual manner from the mass of the end effector in a constant gravitational field whose value is $\mathit{g}$ measured with respect to the moving frame, then the equation of motion is

$$\tau =\mu +\mathsf{M}\alpha ,$$

(88)

where

$$\begin{array}{c}\mu =2\mathit{w}\times \left({\mathsf{M}}_e\mathit{w}\right)+2ϵm_e\mathit{w}\times \mathit{v}\hfill \\ +2m_e((\mathit{w}·{\mathit{r}}_0)(\mathit{w}\times {\mathit{r}}_0)+{\mathit{r}}_0\times (\mathit{w}\times \mathit{v})+ϵ\mathit{w}\times (\mathit{w}\times {\mathit{r}}_0))\\ \hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{\mathrm{\Lambda }}^T{\mathsf{M}}_0\left({\mathcal{L}}_{\phi }\mathrm{\Lambda }\right)\phi -2m_e({\mathit{r}}_0\times \mathit{g}+ϵ\mathit{g}),\end{array}$$

(89)

and

$$\mathsf{M}\alpha =2{\mathsf{M}}_e\dot{\mathit{w}}+2m_e({\mathit{r}}_0\times \dot{\mathit{v}})+2ϵm_e\dot{\mathit{v}}+2ϵm_e(\dot{\mathit{w}}\times {\mathit{r}}_0)+{\mathrm{\Lambda }}^T{\mathsf{M}}_0\mathrm{\Lambda }\alpha .$$

(90)

The various terms in Equation (89) can be interpreted as follows.

• ${\mathsf{M}}_e\dot{\mathit{w}}$ and $m_e\dot{\mathit{v}}$ are inertial resistance to change of angular and translational velocities.
• $m_e\mathit{w}\times \mathit{v}$ is the centripetal force required to rotate and move at the same time.
• $\mathit{w}\times \left({\mathsf{M}}_e\mathit{w}\right)$ is the precession torque (so that if the moment of inertia is not isotropic, then the body spins in a counter-intuitive manner, see, for example, [35]).
• $${\mathrm{\Lambda }}^T{\mathsf{M}}_0\left({\mathcal{L}}_{\phi }\mathrm{\Lambda }\right)\phi +{\mathrm{\Lambda }}^T{\mathsf{M}}_0\mathrm{\Lambda }\alpha =\mathsf{T}{\mathit{f}}_0$$

  (91)

  is the wrench required to move the actuators, where the no-load forces may be computed using

  $${\mathit{f}}_0={\mathsf{M}}_0{\mathcal{L}}_{\phi }\left(\mathrm{\Lambda }\phi \right)+{\mathsf{M}}_0\mathrm{\Lambda }\alpha .$$

  (92)
• $m_e\mathit{g}$ is the force due to gravity.
• All terms containing ${\mathit{r}}_0$ are corrections required since the center of gravity isn’t necessarily the same as the origin of the moving frame of reference. They could be derived by first finding the equations of motion when ${\mathit{r}}_0=0$, and then applying Equations (27) and (28).

Numerical Verification of the Dynamics Equations

The way we numerically verified these equations was by running a simulated motion, and then calculating the work done in three ways. The first method was to integrate the inner product of the twist and the wrench on the end effector. The second method was to integrate the sum of the actuator force times the rate of change of the length of the actuator. The third method was to calculate the total kinetic energy using Equation (86), and the potential energy using the standard mass times gravity time height formula. The numerical simulations gave the same results for all three methods up to machine precision.

12. Conclusions

This paper has given a comprehensive and consistent description of how to use dual quaternions to represent poses, rigid motions, twists, and wrenches. We have introduced the notion of the Lie derivative for dual quaternions. We have shown how these formula are helpful for first producing Newton-Raphson methods for solving the forward kinematics problem for parallel robots, and secondly for a self contained derivation for the dynamic equations of motion of the end effector that includes the inertia of the actuators.

Finally, in Equation (A1), we give an approximation of the normalization of a vector dual quaternion perturbation of the identity, which shows that it is equal up to the second order to the exponential of the vector dual quaternion. This equation was essential for calculating the Hessian in the forwards kinematics algorithms. We feel that this formula will be of independent interest to other researchers in the field of dual quaternions.