DeepAI

# Using dual quaternions in robotics

We advocate for the use of dual quaternions to represent poses and twists for robotics. We show how to represent torques and forces using dual quaternions. We introduce the notion of the Lie derivative, and explain how it can be used to calculate the behavior of actuators. We show how to combine dual quaternions with the Newton-Raphson method to compute forward kinematics for parallel robots. We derive the equations of motion in dual quaternion form. This paper contains results we have not seen before, which are listed in the conclusion.

• 2 publications
• 2 publications
03/23/2022

### An introduction to using dual quaternions to study kinematics

We advocate for the use of dual quaternions to represent poses, twists, ...
10/05/2022

### ProxNLP: a primal-dual augmented Lagrangian solver for nonlinear programming in Robotics and beyond

Mathematical optimization is the workhorse behind several aspects of mod...
03/07/2022

### Minimax principle for right eigenvalues of dual quaternion matrices and their generalized inverses

Dual quaternions can represent rigid body motion in 3D spaces, and have ...
05/31/2022

### Hadamard matrices related to a certain series of ternary self-dual codes

In 2013, Nebe and Villar gave a series of ternary self-dual codes of len...
05/03/2019

### Asymmetric Dual-Arm Task Execution using an Extended Relative Jacobian

Dual-armed coordination is a classical problem in robotics, which gained...
10/25/2019

### DQ Robotics: a Library for Robot Modeling and Control Using Dual Quaternion Algebra

Dual quaternion algebra and its application to robotics have gained cons...
01/11/2018

### Review of theory and implementation of hyper-dual numbers for first and second order automatic differentiation

In this review we present hyper-dual numbers as a tool for the automatic...

## 1. Introduction

We assume that the reader is familiar with the notion of pose and rigid motion, and the use of quaternions to represent rotations. We use the terms pose and rigid motion interchangeably, since a pose can be considered as a rigid motion relative to a fixed reference frame, and the mathematical operations to manipulate them are identical. Thus we also use the terms rotation and orientation interchangeably, and similarly with translation and position. We refer the reader to [5, 10].

The use of dual quaternions to represent poses is well established, especially in the graphic card industry: [1, 12, 6, 16, 17, 19, 26]

. However, in the robotics industry, it is more common to represent a pose using a quaternion (for the rotation) and a three dimensional vector (for the translation). In this paper, we advocate for the use of dual quaternions.

Let us start with a general discussion of why people might want to represent rotations by unit quaternions instead of three by three orthogonal matrices with determinant . One reason given is that quaternions can be represented with only four numbers, whereas the matrix representation requires nine numbers. Another representation commonly used are Euler angles, and these only require three numbers. However, with Euler angles, the formula for composition of rotations is difficult to calculate. Also, Euler angles suffer from gimbal lock, in which if the rotation takes certain values, the choice of Euler angles isn’t unique, and discontinuities and singularities arise. Quaternions suffer from no such problems.

But in modern times, computer memory is very cost effective, and communication is fast, thus the aformentioned advantages might seem trivial at first. This begs the question: what other advantages do quaternions offer over matrices? Quaternions represent a quantity with three degrees of freedom by a four dimensional vector, whereas matrices use a nine dimensional vector. Suppose one obtains an approximation to a matrix or quaternion? For example, this would be common when using Newton-Raphson to find a rotation satisfying some formula, or to interpolate between rotations. So then one has to find a projection which takes the vector to an admissible vector, that is, a vector that actually represents a rotation. And projecting from a nine dimensional quantity to a three dimensional quantity is going to be much harder than the similarly projection from the four dimensional vector. (The projection from quaternions to unit quaternions is called the

normalization, and is defined in Section2. A suggestion of what such a projection would be for matrices is given at the end of Section 16.)

For representing poses, the usual representation is by a three by three orthogonal matrix with determinant

for the rotation, and a three dimensional vector for the translation, twelve dimensional in all. This has a natural representation as a four by four matrix, as we explain Section 4. A great advantage of this representation is the relationship between pose and twist:

 ddt(pose)=(pose)×(twist). (1)

where the multiplication is matrix multiplication, and a twist is a combination of angular velocity and translational velocity measured with respect to the moving frame.

Another natural representation used a lot is to use a quaternion and a three dimensional vector. But this new representation has problems in that it is not clear how the above equation would work in this situation. However, if one introduces the notion of the dual quaternion, which is an eight dimensional vector, equation (1) is still valid if one uses dual quaternion multiplication.

A challenge with dual quaternions is that they are abstract, which can make them harder to understand. Mathematicians find value behind abstract computational ideas, but Engineers usually find the end result, ‘black box approach,’ more useful.

A useful approach to developing intuition for dual quaternions is to understand the rules with their structure, and the mathematical operations needed to manipulate them, (that is, normalization, basis vector multiplication, pose differentiation). Once one becomes comfortable with their form and function, the advantage of duel quaternions over other methods becomes more apparant. For example, one realizes that computing the projection of an eight dimensional vector to the six degrees of freedom of a pose, is also quite simple. (This projection is also called the normalization, and is defined in Section 3.)

This makes it very useful if one is using the Newton-Raphson Method to find a pose satisfying certain properties. We show how the dual quaternion representation efficiently solves the forward kinematics problem.

Finally we consider the dynamics of poses, that is, the equations of motion of a rigid spinning body. In any representation this can be quite difficult, because of both centripetal and precessional effects. But from a theoretical perspective, using the dual quaternion representation is quite slick.

In summary, while the usual representation using Euler angles and translations is conceptually easy to understand, performing calculations can be difficult. Whereas with dual quaternions, while the are initially conceptually hard to understand, when it comes to calculations they are much easier to work with.

## 2. Motivation

A pose, or rigid motion, is a rotation followed by a translation. (Semantically, a pose is the position of the end effector of the robot, and thus different from a rigid motion, but mathematically they are the same if we consider a pose to be a rigid motion with respect to the stationary reference frame.) The rotation is represented by a orthogonal matrix with determinant . A translation is a vector in . These represent the pose :

 r↦Rr+t. (2)

Composition of poses is written from right to left. Thus

 (R1,t1)∘(R2,t2)=(R1R2,t1+R1t2). (3)

One way to represent a pose is using a four by four matrix

 [Rt01]=⎡⎢ ⎢ ⎢⎣R11R12R13t1R21R22R23t2R31R32R33t30001⎤⎥ ⎥ ⎥⎦ (4)

Then composition of poses may be computed simply by multiplying matrices of the form (4). In the literature, the Lie group of rotations is often denoted . The Lie group of poses is often denoted .

The problem with this representation is that twelve numbers are required to represent a pose. Other problems, which we discuss below, is that of interpolating a sequence of poses, and of computing forward kinematics.

One common way to solve both of these problems is to represent the rotation by a unit quaternion [2, 23], which we briefly describe here. A quaternion is a quadruple of real numbers, written as , with the algebraic operations . Its conjugate is , its norm is , its normalization is , its real part is , and its imaginary part is . It is called a unit quaternion if , a real quaternion if , and a pure quaternion if . Note the multiplicative inverse is given by .

We identify three dimensional vectors with pure quaternions, by identifying , , and with the three standard unit vectors. A unit quaternion represents the rotation . A rotation by angle about an axis , where , has two unit quaternion representations: . Composition of rotations corresponds to multiplication of unit quaternions. With some practice, it becomes easier to read a rotation from a quaternion than it does from Euler angles, the standard method of representing rotations.

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

 A⋅B=Re(A¯¯¯¯B)=Re(¯¯¯¯AB). (5)

Then a way to represent poses is as a pair . This is the way poses are represented internally in the Robotic Operating System [25]. Then the composition rule is

 (Q1,t1)∘(Q2,t2)=(Q1Q2,t1+Q1t2¯¯¯¯Q1). (6)

A difficulty with this is that the composition is not a bilinear operation as it was with the representation (4). (This lack of bilinear representation becomes particularly problematic when representing twists, which we describe below.) However there is a way to resolve this as follows. Write and . Then equation (6) can be rewritten as

 (Q1,2B1Q−11)∘(Q2,2B2Q−12) (7) =(Q1Q2,2B1Q−11+2Q1B2Q−12Q−11) =(Q1Q2,2(B1Q2+Q1B2)(Q1Q2)−1).

This leads to representing the pose by a unit dual quaternion, which we now describe.

## 3. Dual quaternions to represent poses

A dual quaternion is a pair of quaternions, written as , with the extra algebraic operation .

The conjugate dual quaternion of is . Conjugation reverses the order of multiplication:

 ¯¯¯¯¯¯¯¯¯¯η1η2=¯¯¯η2¯¯¯η1. (8)

There is another conjugation for dual quaternions: , but we have no cause to use it in this paper, except in equation (22) below.

A unit dual quaternion is a dual quaternion such that , equivalently, that is a unit quaternion and . A pure dual quaternion is a dual quaternion such that both and are pure quaternions.

If is a dual quaternion with , then its multiplicative inverse can be calculated using the formula

 η−1=Q−1−ϵQ−1BQ−1. (9)

If is a unit dual quaternion, then there is a computationally much faster formula:

 η−1=¯¯¯η. (10)

For a dual quaternion , it is not really possible to mix and additively. The quaternion is unitless, whereas the quaternion has units of length. For this reason, when measuring how large a dual quaternion is, everything must be with respect to , the characteristic length scale. The size of a dual quaternion is defined to be

 ∥A+ϵB∥=(|A|2+l−2|B|2)1/2. (11)

A dual number is anything of the form , where and are real numbers. The norm of a dual quaternion is the dual number defined by the two steps:

 |η|2=¯¯¯ηη=η¯¯¯η=|Q|2+2ϵ(B⋅Q), (12) |η|=√|η|2=|Q|+ϵ(B⋅Q)/|Q|. (13)

The norm preserves multiplication, that is, if and are two dual quaternions, then

 |η1η2|=|η1||η2|. (14)

If is any dual quaternion with , then we define its normalization to be the unit dual quaternion

 ˆη=|η|−1η=η|η|−1=Q/|Q|+ϵ(B−(B⋅Q)Q/|Q|2). (15)

(We remark that the normalization of a dual quaternion is used in the computer graphics industry [16, 17].) While this normalization formula might seem initially quite complicated, after thinking about it one can see that it is the simplest projection that enforces and .

The normalization also satisfies the following properties.

• If is a unit dual quaternion, then .

• Normalization preserves multiplication, that is, if and are two dual quaternions, then

 ˆη1η2=ˆη1ˆη2. (16)
• If is pure, then

 ˆ(1+θ)=1+θ+12θ2+O(θ3). (17)

(See Section 15 for the proof of this.)

Note that

 η−1=|η|−2¯¯¯η, (18)

and this could have been the definition of the multiplicative inverse of a dual quaternion, except that this definition is circular.

Let us clarify the big-Oh notation. We say that

 η(θ)=O(γ(θ)) (19)

if there exists a constant such that if is sufficiently small, then

 ∥η(θ)∥≤c∥γ(θ)∥. (20)

We set for the set of dual quaternions, for the set of unit dual quaternions, and for the set of pure dual quaternions.

If we are given a pose represented by , then the pose is also represented by the unit dual quaternion

 η=Q+12ϵtQ. (21)

As we have shown above, composition of poses corresponds to multiplication of unit dual quaternions. If is a 3-vector, and is the image of under the action of the pose , then

 1+ϵs=η(1+ϵr)¯¯¯η∗, (22)

but generally it is easier to use the formula

 s=Qr¯¯¯¯Q+2B¯¯¯¯Q=(Qr+2B)¯¯¯¯Q. (23)

The introduction of the factors and is unnecessary, but it slightly improves formulas for the twist, as we see below.

It is difficult for a human to look at the second part of a unit dual quaternion, and from that read what the translational part of the pose should be. So dual quaternions work better for an internal representations of poses rather than human readable representations.

One more notion is the exponential of a dual quaternion [26]:

 exp(θ)=∞∑k=0θkk!, (24)

The exponential map takes pure dual quaternions to unit dual quaternions. In particular, if is a pure dual quaternion, then by equation (17) we have

 exp(θ)=ˆ(1+θ)+O(θ3). (25)

We show how to calculate in Section 14.

## 4. Kinematics: dual quaternions to represent twists

The translational velocity of a translating reference frame is given by . The angular velocity of a rotating frame is given by the differential equation

 ˙R=Rw⋆, (26)

or , where is the Hodge star operator of [27]:

 w⋆=⎡⎢⎣0−w3w2w30−w1−w2w10⎤⎥⎦. (27)

(Equation (26) gives the angular velocity with respect to the rotating reference frame. If one wants the angular velocity with respect to the stationary reference frame, use instead .)

Note that

 (w⋆)r=w×r, (28)

where denotes the cross product.

A twist is the pair of vectors that describe the change of pose in the moving reference frame, that is:

 ddt[Rt01]=[Rt01][w⋆v00]. (29)

(Note in particular that , so that a continual application of a twist can result in a motion following a straight line, an arc of a circle, or a spiral.) The set of angular velocities, and the set of twists are both Lie algebras, often denoted by and in the literature.

The quaternion-translation formulation doesn’t have a representation like equation (1):

 ddt(Q,t)=(12Qw,Qv¯¯¯¯Q). (30)

But with the unit dual representation, because the composition of poses is given by a bilinear multiplication, it follows that the twist is represented by the dual quaternion satisfying the differential equation.

 φ=η−1˙η,or˙η=ηφ. (31)

It will be shown in Section 15 (see also [12]) that is the pure dual quaternion

 φ=12w+12ϵv. (32)

Note that if is a unit dual quaternion that is a function of time , then is always a pure dual quaternion. (Note that is a pure dual quaternion equal to the twist relative to the fixed frame of reference. Note also that the factor attached to in equation (32) is a consequence of the introduction of the factor in equation (21), and thus is essentially arbitrary.)

## 5. Dual quaternions to represent wrenches

Let the pose 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

 τ=2q+2ϵp, (33)

where and 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 is the center of mass of the end effector in the moving frame, then the twist about the center of mass is given by

 φ0=φ+12ϵw×r0, (34)

where , and the wrench applied about the center of mass is

 τ0=τ+2p×r0. (35)

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

 ddt(work done)=τ⋅φ=τ0⋅φ0. (36)

See [4] for the origins of the term twist and wrench as pairs of 3-vectors, which are examples of screws.

## 6. Interpolation of poses

Suppose that we are given a sequence of times and rotations, and . We would like to find a function from to rotations such that the third derivative of is bounded, and . If the quantities were merely vectors, we could use cubic splines. But if we perform a similar calculation on the matrices, we cannot guarantee that the matrix is a rotation matrix.

So what we want is to find a map that projects general matrices onto rotation matrices. But it is not obvious how to find such a map that is computationally efficient. It is here that using quaternions really has great advantages over matrices, because an obvious map from general quaternions to unit quaternions is to normalize.

Let be the quaternion representation of respectively. Let be the cubic spline interpolation such that . Then

 F(t)=1|G(t)|G(t) (37)

provides a computationally fast, low jerk, method of interpolating the rotations.

This can be extended to interpolate poses. Given times and unit dual quaternions , first interpolate to get a function such that , and then normalize it to get a function . This gives a computationally fast, low jerk, method for interpolating poses.

It is perhaps better to interpolate the rotation quaternions and the translations separately. This is because the center of gravity typically travels in a straight line rather than the arc of a curve, but naive interpolation of dual quaternions follows arcs of curves if the angular rotation is non-zero. (Interpolation of dual quaternions is more appropriate for ‘skinning’ [16, 17].)

If one wants to find a function that passes through these points, and has low jerk even at the end points and , this can be done by adding a ‘ramp up’ and ‘ramp down’ at the ends of the trajectory. Create a free cubic spline on and . Then define a function

 s(t)=⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩12(t0+t1)if t

where

 f(t)=12+12(2−t)t3. (39)

Since , , and , it follows that has bounded third derivative, for , for , and for . Thus

 η(t)=ˆγ(s(t)) (40)

provides a low jerk function, even at the end points, that passes through the appropriate points. it is important that and be large enough to allow the trajectory to ‘ramp up’ and ‘ramp down’ smoothly.

## 7. Perturbations of poses

Suppose one has a pose that varies in time, and a fixed time . Then we would like to find a good representation of the perturbation of from a reference pose when is close to . One way to do this is to choose a new reference frame in which , the identity pose, or rather, to think about how behaves.

It can be seen that for close to that there is a pure dual quaternion such that

 η=ηrˆ(1+θ)≈ηr(1+θ). (41)

Thus

 θ≈η−1rη−1=¯¯¯ηrη−1 (42)

gives a great way to measure the perturbation.

We prefer instead to use what we shall call the Lie difference, which we denote using , which is close to (42), but is always a pure dual quaternion:

 θ≈η△ηr=12(¯¯¯ηrη−¯¯¯ηηr)=Im(¯¯¯ηrη). (43)

Note that if and are pure dual quaternions, and is an invertible dual quaternion, then

 (η(1+θ1))△(η(1+θ2)) (44) =(1+θ1)△(1+θ2) =θ1−θ2.

Thus these give a great way to map perturbations of non-linear poses to linear pure dual quaternions. These can be very effective in control theory.

## 8. Dual quaternions as vectors

We make the identifications

 D≅R8, (45)

using the basis

 {β1,β2,β3,β4,β5,β6,β7,β8}={i,j,k,ϵi,ϵj,ϵk,1,ϵ}, (46)

and similarly, we make the identification

 sd≅R6, (47)

using the basis , , , , , . With these identifications, we can define the dot product between two dual quaternions by transferring the usual definition of dot product on , that is

 (A+ϵB)⋅(C+ϵD)=A⋅C+B⋅D. (48)

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

 η=8∑i=1ηiβi, (49)

and every pure dual quaternion as

 θ=6∑i=1θiβi. (50)

## 9. 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, 30]. 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 [20, 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 [13, 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 pure dual quaternions.

If one has a quantity that is a function of pose , then we usually think of its derivative as the Jacobian with respect to the components of . But it really makes more sense to compute the derivative with respect to the components of the perturbation of . The latter is the Lie derivative.

The definition is this. Given a differentiable function whose domain is the unit dual quaternions, , we can extend it arbitrarily to a differentiable function whose domain is an open neighborhood of in . Given a unit dual quaternion and a pure dual quaternion , we define the Lie derivative of in the direction of to be

 Lθg=limr→0g(η(1+rθ))−g(η)r=ddrg(η(1+rθ))∣∣∣r=0. (51)

Since 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 was extended from , but it is, as is shown in Lemma 1 below.

Given a generic function whose domain is the dual quaternions, , we define its Jacobian to be the dual quaternion

 ∂g∂η=8∑i=1∂g∂ηiβi. (52)

If its domain is the pure dual quaternions, , we have the same formula except with replaced by .

Using the chain rule for partial derivatives, we obtain the following formula, which is useful for explicitly calculating the Lie derivative if the function

is known.

 Lθg=∂g∂η⋅(ηθ). (53)

We define the partial Lie derivatives to be

 Lig=Lβig,(1≤i≤6), (54)

and its full Lie derivative to be the pure dual quaternion

 Lg=6∑i=1Ligβi, (55)

so that for all pure dual quaternions it satisfies:

 θ⋅Lg=Lθg. (56)

To gain some intuition, write

 θ=12a+12ϵb, (57)

so that

 ∂∂θ=2∂∂a+2ϵ∂∂b. (58)

Then we see that represents a change in pose by an infinitesimal translation and an infinitesimal rotation , measured in the moving frame of reference. Thus is a pure dual quaternion giving twice the change in with respect to an infinitesimal rotation, plus times twice the change of with respect to an infinitesimal translation.

One important property of the Lie derivative is that if represents a pose, with twist , then

 ddt[g(η)]=Lφg. (59)

The Lie derivative satisfies various rules, which are also useful for explicitly calculating the Lie derivative when is known.

• If is linear in , then

 Lθg(η)=g(ηθ). (60)
• The product rule: if is any product, such as inner product or dual quaternion product, then

 Lθ(g1×g2)=g1×(Lθg2)+(Lθg1)×g2. (61)
• The chain rule:

 Lθ(h(g1,g2,…,gm))=m∑i=1∂∂gih(g1,g2,…,gm)Lθgi. (62)

## 10. Applications to parallel robots

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

 ℓ=(ℓj)1≤j≤n. (63)

For example, for a cable-driven parallel robot, these represent the lengths of the cables, and typically . For the Gough or Stewart Platforms [10], we have .

Let us also denote the force exerted by the actuators by

 f=(fj)1≤j≤n, (64)

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

 ddt(work done)=f⋅˙ℓ. (65)

Suppose we have a function , which calculates the required actuator values, , from the pose of the moving frame. This is the inverse kinematics function. If it is a simple cable-driven parallel robot, this can be computed using Pythagoras’s Theorem, from the positions of the actuators on the fixed frame, and the positions of where the cables are attached to end effector. But this formula could be more complicated in other situations.

We also define the matrix by

 Λθ=LθL. (66)

From equation (59) we obtain

 ˙ℓ=Λφ. (67)

There is also a () matrix that maps the actuator forces to the wrench dual quaternion:

 τ=Tf. (68)

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:

 T=ΛT. (69)

This is because the rate of change of work done on the parallel robot can be computed in two different ways, either using equation (36), or (65). Substituting in equations (67) and (68), we obtain

 Tf⋅φ.=f⋅Λφ. (70)

Since this is true for arbitrary actuator forces and end effector twists , the result follows.

## 11. Second lie derivatives

If is a function of dual quaternions, we define its Hessian to be the matrix

 [∂2g∂η2]=[∂2g∂ηi∂ηj]1≤i,j≤8. (71)

Thus the expression should be interpreted with treated as an eight dimensional vector.

We also need to compute second Lie derivatives. Calculations show that

 LθLψg=(ηψ)⋅[∂2g∂η2](ηθ)+∂g∂η⋅(ηθψ). (72)

As a corollary, we obtain the well known identity:

 LθLψg−LψLθg=L(θψ−ψθ)g, (73)

which implies that Lie derivatives do not necessary commute.

Another way one might try to define the second derivative is to use the formula . Unfortunately, this definition doesn’t work, as it depends upon the choice of how to extend the domain of to all dual quaternions. The obvious choice of extension is to use the normalization:

 ~g(η)=g(^η). (74)

We have the following formula for the Hessian of :

 ∂2∂θ2~g(η(1+θ))∣∣∣θ=0=[12(LiLjg(η)+LjLig(η))]1≤i,j≤6. (75)

This equation is proved in Section 15.

## 12. Forward kinematics for parallel robots

Let the set of admissible actuator values, , be the range of the function . Then the forward kinematics function is

 Y:L→SD, (76)

which is a left inverse to . Because of possible measurement errors, should produce decent answers even if the actuator values are merely close to .

We will focus on the over-constrained problem, that is, when the number of actuators is greater than . If , a much simpler approach is possible by simply applying the Newton-Raphson method directly to . In that case, only first derivatives of are required.

This problem has been solved by many others, for example, [22, 31]. But we feel that this is much more easily solved using dual quaternions. The function takes actuator values , and seeks to find the pose so that is close as possible to . We do this by seeking to minimize the loss function

 b(η)=12|L(η)−ℓ|2 (77)

We use the Newton-Raphson Method, which given close to a minimum of the loss function, finds which is much closer to the same minimum of the loss function. Because this method as usually stated works only on linear vector spaces, we first have to define a map , with , from the linear vector space of pure dual quaternions to the manifold of unit dual quaternions [15].

Most papers on the Newton-Raphson Method on manifolds construct this map using the so called exp function [7, 8, 9]. So the map is

 ~bk(θ)=b(ηkexp(θ)). (78)

The exp map in these papers is following the path of a geodesic on the manifold, and this is equivalent to using the equations of motion of the end effector as described in Section 13. Another exp map is to follow a one-parameter subgroup, or equivalently, equation (24). We do not take these approaches, as they can be computationally expensive. Our approach is to normalize:

 ~bk(θ)=b(ηkˆ(1+θ)). (79)

This is numerically close to the second approach, as is shown by equation (25).

In order to apply the Newton-Raphson, we need to compute the Jacobian

 δk=∂~bk∂θ∣∣∣θ=0=(Lib(ηk))1≤i≤6, (80)

and its Hessian , which by equation (75) is

 Hk=[12(LiLjb(ηk)+LjLib(ηk))]1≤i,j≤6. (81)

The Newton-Raphson method we use is to iterate:

 ηk+1=ηkˆ(1−H−1kδk). (82)

For the criterion of when to terminate the Newton-Raphson Method, we measure .

This method works well in practice, because typically we measure the pose frequently, and thus the current pose won’t be much different than the previously measured pose. Thus the previously measured pose is a good starting value for using the Newton-Raphson method for finding the current pose. We found that only about two or three iterations of the Newton-Raphson method were generally required, and that it easily ran at a thousand times per second on a modern computer.

## 13. Dynamics of the end effector

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

 f0=M0¨ℓ, (83)

where is a positive definite matrix denoting what we shall call the effective no-load mass of the actuators.

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

 e=12φ⋅Mφ, (84)

where is a positive definite matrix, which depends only upon , and which we call the effective mass of the parallel robot. If is the mass of the end effector,

is the moment of inertia tensor of the end effector about its center of mass, and

is the center of mass of the end effectors, all measured with respect to the moving frame, then

 e=12me|v+w×r0|2+12w⋅Mew+12˙ℓ⋅M0˙ℓ, (85)

that is

 M=4[Me−me(r⋆0)2mer⋆0−me(r⋆0)meI]+ΛTM0Λ, (86)

where is the identity matrix.

The next result is proved in Section 15.

###### Theorem 1.

If the kinetic energy satisfies equation (84) with equation (86) holding, and the potential energy is calculated in the usual manner from the mass of the end effector in a constant gravitational field , then the equation of motion is

 τ=μ+Mα, (87)

where

 μ=2w×(Mew)+2ϵmew×v+2me((w⋅r0)(w×r0)+r0×(w×v)+ϵw×(w×r0))+ΛTM0(LφΛ)φ−2me(r0×~g+ϵ~g), (88)

with

 ~g=¯¯¯¯QgQ, (89)

and

 Mα=2Me˙w+2me(r0×˙v)+2ϵme˙v+ΛTM0Λα. (90)

The various terms in equation (88) can be interpreted as follows.

• and are inertial resistance to change of angular and translational velocities.

• is the centripetal force required to rotate and move at the same time.

• 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, [18]).

•  ΛTM0(LφΛ)φ+ΛTM0Λα=Tf0 (91)

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

 f0=M0Lφ(Λφ)+M0Λα. (92)
• is the force due to gravity, translated into the moving frame of reference.

• All terms containing 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 , and then applying equations (34) and (35).

## 14. Spherical linear interpolation of dual quaternions

Suppose we are given two unit dual quaternions, and ? The problem is to find the dual quaternion valued function of time, denoted , such that at , it takes the value , that at , it takes the value , and such that its twist, , is a constant pure dual quaternion. The name of the function, ‘slerp’, comes from the equivalent problem for orientations [24].

This requires two functions. The first is the exponential function, defined in equation (24). The function as function of time is the pose traveled if it is the identity pose at , and maintains a constant twist . The second is a logarithm function, which given a unit dual quaternion , finds such that . There are infinitely many possibilities for , so we will just describe the principal value, that is, the one for which the norm of the angular velocity is minimized. Then

 slerp(η1,η2,t) =η1exp(tlog(¯¯¯η1η2)) (93) =exp(tlog(η2¯¯¯η1))η1.

If one is merely interested in interpolating between the poses represented by these unit dual quaternions, multiply one of them by so that non-dual part of has non-negative real part.

The results of this section may be found in [21]. We have seen similar results in [26], but we believe our formulas to be more explicit.

Given any dual quaternion

 η=c+d+ϵ(x+y), (94)

by decomposing into a vector parallel to and perpendicular to

, we can suppose without loss of generality that there exists orthonormal vectors

and such that is a linear combination of , , , , and .

###### Theorem 2.

If and are orthonormal vectors, then if

 exp(12wta+(12v1ta+12v2tb)ϵ) (95) =(cos(12wt)+sin(12wt)a)(1+12v1taϵ) =+v2wsin(12wt)bϵ,

and if

 exp(12vtaϵ)=1+12vtaϵ. (96)

If , this represents a pose that rotates counterclockwise with angular velocity around , and translates by

 v1ta+v2w(sin(wt)b+(1−cos(wt))(a×b)),

that is, the sum of a point on a line in the direction of moving with speed , and a point on a circle in the plane perpendicular to of radius which is traversed once every time units. If , this represents a pose with the identity rotation and translating in a straight line with velocity .

###### Theorem 3.

Suppose and are orthonormal vectors, and that

 η=c+sa+ϵ(x+y1a+y2b) (97)

is a unit dual quaternion, with . Let

 t=\rm atan2(s,c), (98)

that is, the angle part of the polar coordinates of , with . Then if then

 log(η)=ta+ϵ(cy1−sx)a+ϵty2sb, (99)

where if we assume , and we set .

It is not possible to assign a canonical principal value for , because it could be for any unit vector .

## 15. Proofs

###### Proof of Equation (17).

Since , we have

 |1+θ|2=1−θ2. (100)

Hence using Taylor’s series

 |1+θ|−1=1+12θ2+O(θ3) (101)

from which it follows that

 ˆ(1+θ)=(1+θ)(1+θ)−1=1+θ+12θ2+O(θ3). (102)

###### Proof of Equation (32).

First, if we differentiate , we obtain , that is, is pure.

Suppose the position is the image of the constant position under the pose represented by and also represented by . Differentiating equation (2) we obtain

 ˙s=Rw⋆r+Rv=R(w×r+v). (103)

Differentiating equation (22) we obtain

 ϵ˙s=ddt(η(1+ϵr)¯¯¯η∗)=η(φ(1+ϵr)+(1+ϵr)¯φ∗)¯¯¯η∗. (104)

After some manipulation, and setting , we obtain

 ˙s=Q(a×r+b)¯¯¯¯Q. (105)

Comparing equations (103) and (105), we obtain

 w×r+v=a×r+b, (106)

and since this holds for all , the result follows. ∎

###### Lemma 1.

The definition of in equation (51) does not depend upon the extension of from to a neighborhood of in .

###### Proof.

Let and be two extensions of from to a neighborhood of in