USING LIE DERIVATIVES WITH DUAL QUATERNIONS FOR PARALLEL ROBOTS

. 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.


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 [7] 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 [5,6,13,36].
The notion of the dual quaternion, and its use to represent poses and rigid motions, seems to go back to McAulay [25], inspired by the earlier work by Clifford [8].The notion of using dual quaternions to represent twists may be found in [1,31].A basic introduction to dual quaternions may be found in [22,29], the latter also covering twists.Many authors have used dual quaternions to represent hierarchies of poses, that is, chains of manipulators [22,31,34,35,38].Papers on representing kinematics or dynamics via dual quaternions include [1,2,10,16,23,39,40]. Dual quaternions have also found great use in computer graphics [19,20].
(The reader should be aware that [1,16] have incorrect formulas for the logarithm and exponential of dual quaternions -the correct formulas may be found in [28], and [37] 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 [27], to create software for controlling a cabledriven 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, Section 12.

Notation
Rotations can be represented by unit quaternions [3,33], which we briefly describe here.A quaternion is a quadruple of real numbers, written as A = w + xi + yj + zk, with the algebraic operations its real part is Re(A) = w = 1 2 (A + A * ), and its imaginary part is Im(A) = ix+jy +kz = 1 2 (A−A * ).It is called a unit quaternion if |A| = 1, a real quaternion if Im(A) = 0, and a vector quaternion if Re(A) = 0.If A ̸ = 0, the multiplicative inverse is given by A −1 = A * /|A| 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 i, j, and k with the three standard unit vectors.A unit quaternion Q represents the rotation which takes the direction r to QrQ * .A rotation by angle a about an axis s, where |s| = 1, has two unit quaternion representations: ±(cos( 12 a) + s sin( 1 2 a)) = ± exp( 1 2 as).Composition of rotations corresponds to multiplication of unit quaternions.
We can represent quaternions as four dimensional vectors, and give it the inner product A dual quaternion is a pair of quaternions, written as η = A + ϵB, with the extra algebraic operation ϵ 2 = 0. We call A = P(η) the primary part of η, and B = D(η) the dual part of η.
The conjugate dual quaternion of η = A + ϵB is η * = A * + ϵB * .Conjugation reverses the order of multiplication: There is another conjugation for dual quaternions: A + ϵB = A − ϵB, but we have no cause to use it in this paper, except in equation ( 7) below.
A unit dual quaternion η = Q + ϵB is a dual quaternion such that η * η = 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 η = A+ϵB is a dual quaternion with A ̸ = 0, then its multiplicative inverse can be calculated using the formula If η is a unit dual quaternion, then there is a computationally much faster formula (see [1]): ) We set DH for the set of invertible dual quaternions (that is, A + ϵB where A ̸ = 0), dh for the set of dual quaternions, SDH for the set of unit dual quaternions, and sdh for the set of vector dual quaternions.
A rigid motion of the form where here t is a translation, can be represented by the unit dual quaternion Composition of rigid motions corresponds to multiplication of unit dual quaternions, where the notation η 1 η 2 means to apply first the rigid motion represented by η 2 , and then by η 1 , that is, the dual quaternion acts by left multiplication.If r is a 3-vector, and s is the image of r under the action of the rigid motion η = Q + ϵB, then but generally it is easier to use the formula 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 A twist is the pair of vectors (w, v) that describes the rate of change of pose or rigid motion, where w is the angular velocity, and v is the translational velocity.One can think of the twist (w, v) as a rigid motion function of time t: 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 and if it is understood to be with respect to the fixed frame Then this twist can be represented by a vector dual quaternion [1,2] where if we understand it to be with respect to the moving frame, we have the formula φ = η −1 η, or η = ηφ, (14) and if we understand it to be with respect to the fixed frame In this paper, unless otherwise stated, we always understand the twist to be with respect to the moving frame.We make the identification using the basis and similarly, we make the identification using the basis (β 1 , β 2 , β 3 , β 4 , β 5 , β 6 ).With these identifications, we can define the dot product between two dual quaternions by transferring the usual definition of dot product on R 8 , that is In this way, every dual quaternion η can be written in component form as and every vector dual quaternion θ as Finally, we give a few extra formulas.Let φ m denote the twist with respect to the moving frame, and φ f denote the twist with respect to the fixed frame.From equations ( 14) and ( 15) we obtain Since in any algebra we have we obtain the surprisingly simple formula for the change of frame for acceleration: φf = η φm η −1 .

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 where q and 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 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 where φ = η −1 η, and the wrench applied about the center of mass is 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 d dt (The second equality follows from vector identities.)See [5] for the origins of the term twist and wrench as pairs of 3vectors, which are examples of screws.The 'work done' formulas are also known as reciprocal screw relationships.

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 η = A + ϵB is the dual number defined by the two steps: The norm preserves multiplication, that is, if η 1 and η 2 are two dual quaternions, then If η = Q + ϵB is an invertible dual quaternion, then we define its normalization to be the unit dual quaternion (We remark that the normalization of an invertible dual quaternion is used in the computer graphics industry [19,20].)While this normalization formula might seem initially quite complicated, after thinking about it one can see that it is the simplest projection that enforces The normalization also satisfies the following properties.
• Normalization preserves multiplication, that is, if η 1 and η 2 are two dual quaternions, then

Notation for three by three matrices
Let I denote the (3 × 3) identity matrix, and 0 denote the (3 × 3) zero matrix.If r is a 3-vector, then the Hodge star operator of r is Note that If u is a unit vector, then the projection onto the complement of the unit vector u is defined by Consistent with the identification (18), if A, B, C, and D are (3 × 3) matrices, and θ = 1 2 a + 1 2 ϵb, ψ = 1 2 c + 1 2 ϵd are vector dual quaternions, † Note this formula is incorrect in the published version [30].The authors are grateful to Brent Koogler for bringing this to our attention.This correction makes only a small difference to the results reported in Subsection 10.3.
we have

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 [41].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 [24].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 [17,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(η), 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(η) 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 g whose domain is the unit dual quaternions, SDH, 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 SDH in dh.Given a unit dual quaternion η and a vector dual quaternion θ, we define the Lie derivative of g(η) in the direction of θ to be Since η(1+rθ) 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 SDH, but it is, as is shown in Lemma 2 below.Given a generic function g whose domain is the invertible dual quaternions, DH, and whose codomain is any vector space over the real numbers, we define its Jacobian to be the dual quaternion where the partial derivative ∂ ∂η i is interpreted using the identification of η i with the β i components of η as described in (17).
If the domain of g is the vector dual quaternions, sdh, 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 SDH in dh is known: where the notation (ηθ) i means the β i component of ηθ, as described in (17).
We define the partial Lie derivatives to be 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) so that for all vector dual quaternions θ it satisfies: To gain some intuition, write Since we have that we see that θ represents a change in pose by an infinitesimal translation b and an infinitesimal rotation a, measured in the moving frame of reference.Thus Lg 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 η represents a pose, with twist φ, then by equation ( 41), with r replaced by t, we see that The Lie derivative satisfies various rules, which easily follow from either equations ( 41) or (43), which are also useful for explicitly calculating the Lie derivative when g is known.
• If g(η) is linear in η, then • 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 • The chain rule: • Let s be a constant position vector, and ñ be a constant direction.Let s and n be their corresponding values with respect to the moving frame.Then 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.

Applications to parallel robots
To aid our description, we show the example of a cable-driven parallel robot in Figures 1 and 2. In Figure 1, we show the end effector in green, controlled by eight cables shown in red, which go around pulleys shown Suppose that the position of the end effector of a parallel robot is given by n actuators, described by quantities For example, for a cable-driven parallel robot, these represent the lengths of the cables, and typically n = 8.The numbers ℓ j only need to be determined up to a fixed number.Thus, for example, in Figure 2, the number ℓ 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 [13], we have n = 6, and ℓ 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 defined so that the rate of change of work performed through the actuators is given by d dt Suppose we have a function L : SDH → R n , which calculates the required actuator values, ℓ, from the pose η of the end effector frame.This is the inverse kinematics function.
We also define the (n × 6) matrix Λ by or more explicitly, by From equation ( 49) we obtain There is also a (6 × n) matrix T that maps the actuator forces to the wrench dual quaternion: 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:

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 × 8) matrix Thus the expression ∂ 2 g ∂η 2 γ should be interpreted as a matrix product with γ 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 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 g to all dual quaternions.The obvious choice of extension is to use the normalization g(η) = g(η).We have the following formula for the Hessian of g: . 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, r 1 , r 2 , . . ., r n where the cables or legs attach, the cables or legs are attached at the other end to s k , and unit vectors u 1 , u 2 , . . ., 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 Then or representing vectors as columns, we have For calculating the second Lie derivative of ℓ k , we need only know the first Lie derivative of 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 and therefore 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 n, and passing through the point x.Attached to the assembly, at a fixed distance w from the fixed point x, by a rod perpendicular to 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 r.A cable passes over the pulley in the n direction from the center of the pulley.All of the vectors are expressed in the moving frame of reference.
We work in the (m, n) coordinate system in which the origin is r.Let (h 1 , h 2 ) be the coordinates of the center of the pulley.Then where The cable length, up to an additive constant, is To avoid singularities, the optimal way to compute the inverse tangent, at least with the configuration shown in Figure 2, is to use atan2(−u 2 , −u 1 ), where the commonly available atan2(y, x) function solves for θ where x = r cos θ, y = r sin θ, r > 0, −π < θ ≤ π.
The various required quantities are 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.9.1.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(Λ) = det(T) = 0.This is in agreement with the results obtained by Gosselin and Angeles [15].
Note that when considering more complex parallel robots, that a singularities can happen in other situations as well (see, for example, [26]).

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 ℓ k .These numbers lengths can be found easily and with high sampling frequency using, for example, encoders.
Let the set of admissible actuator values, L ⊂ R n , be the range of the function L. Then the forward kinematics function is which is a left inverse to L. Because of possible measurement errors, Y should produce decent answers even if the actuator values are merely close to L.
The methods described here are essentially the Newton-Raphson Method, and are all iterative methods.Given a guess η k , we create a new guess η k+1 : and iterate until some criterion is met.We measure size l (η k+1 − η k ), and see when it is smaller than some preset value, like 10 −16 .
Here we merely describe the method.In Section 12, we explain why they work.
We would also like to mention a different approach to using dual quaternions to solve forward kinematics problems in [42], 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 This is easy to code, certainly simpler than methods which use Euler angles.
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, [32,42].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 and .
(82) The algorithm is 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, 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 • 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.
All runs were successful.The average number of required iterations was about 4.8.The run time for each forwards kinematics calculation was a little under 100 microseconds, using a fairly modern but low-end laptop.Increasing the allowed angle to 45 • gave a failure rate of 2 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 1,000 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 milliseconds.
For the second test, we again created 1,000 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 88% 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 1,000 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 millisecond), and the algorithm has never failed to solve the forward kinematics problem.

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 where M is a (6 × 6) positive definite matrix, which depends only upon η, 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: where M 0 is a positive definite (n × 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 M 0 is simply a constant multiple of the identity.Since If m e is the mass of the end effector, M e is the moment of inertia tensor of the end effector about its center of mass, and r 0 is the center of mass of the end effectors, all measured with respect to the moving frame, then that is 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 g measured with respect to the moving frame, then the equation of motion is where and The various terms in equation ( 89) can be interpreted as follows.
• M e ẇ and m e v are inertial resistance to change of angular and translational velocities.
• m e w × v is the centripetal force required to rotate and move at the same time.• w × (M e w) is the precession torque (so that if the moment of inertia is not isotropic, then the body spins in a counterintuitive manner, see, for example, [21]).
is the wrench required to move the actuators, where the no-load forces may be computed using • m e g is the force due to gravity.
• All terms containing 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 r 0 = 0, and then applying equations ( 27) and (28).
11.1.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. Appendix: Proofs Lemma 1.If θ is a vector dual quaternion, then We remark that since the exponential of a unit dual quaternion θ [40] satisfies then by equation (93) we have Let us clarify the big-Oh notation.We say that if there exists a constant c > 0 such that if size l (θ) is sufficiently small, then size l (η(θ)) ≤ c size l (γ(θ)).
(97) Let us show that this definition does not depend upon the characteristic length l.In [28], we show how to define the 'functional calculus' of dual quaternions.That is, given a continuously differentiable function f : C → C, satisfying f (z) = f (z), we can make sense of f (η) for any dual quaternion η = A + ϵB.First we can define f on quaternions by realizing that any unit vector n satisfies n 2 = −1, and hence one merely treats n as an imaginary unit.Next, if B is decomposed as where B 1 commutes with A, and b 2 is a vector quaternion that anti-commutes with Im(A), and setting then we have This means that for any functions f and g, that size l (f (A + ϵB)) ≤ size l (g(A + ϵB)) for all quaternions B if and In particular, the comparison of size does depend upon which characteristic length l is used.
Proof of Lemma 1.Since θ * = −θ, we have Hence using Taylor's series from which it follows that □ Lemma 2. The definition of L θ g in equation (41) does not depend upon the extension of g from SDH to a neighborhood of SDH in DH.
Proof.First note that from equation (93) we have Let g 1 and g 2 be two extensions of g from SDH to a neighborhood of Then where the second equality follows from the chain rule and equation (105).Similarly for g 2 .□ Proof that Definition (41) implies Equation (43).Using the chain rule for partial derivatives, we obtain Now set r = 0. □ Proof that Equation(47) implies Equation (48).Equation (47) can be written as where θ i is the β i component of θ as described in (17).If f is any function of the vector dual quaternions, we have

□
Proof of Equations (53) and (54).We have or that is Remembering a * = −a, and that as − sa = 2 a × s, we obtain

□
Proof of Equation (62).The rate of change of work done on the parallel robot can be computed in two different ways, either using equation ( 29), or (57).Substituting in equations ( 60) and (61), we obtain the last equality being a standard formula for the transpose of a matrix.Since this is true for arbitrary actuator forces f and end effector twists φ, the result follows.□ Proof of Equation (64).Applying the rules given in Section 6, we obtain

□
As a corollary to equation (64), we obtain the well known identity: which implies that Lie derivatives do not necessary commute.
Proof of Equation (65).We wish to find the Jacobian δ and the Hessian H of b at the origin, where We can find this by considering its Taylor series expansion b(θ) = b(0) + where Using the Taylor series, and using equation ( 93), one obtains Now by comparing coefficients, and considering equations ( 43) and (64), we obtain

□
Proof of Equation (67).Suppose that the end effector is moving with pose η = Q + ϵB, and twist ϕ = 1 2 w + 1 2 ϵv, then with respect to the fixed frame of reference, the velocity of the attachment point r k is Q(v + w × r k )Q * .And if a force f k is applied along the direction u k , then the force applied to the end-effector is f k Qu k Q * with respect to the fixed frame of reference.
Hence computing the rate of change of virtual work, we obtain (133) The result follows by equation (60).□ For complex situations, where the cables might pass through pulleys, the simplicity of equation ( 67) can be a bit surprising.This might best be intuitively understood by seeing that the involute of the curve describing the shape of the pulley is the curve traced by the end effector attachment point when the cable length of that actuator is kept fixed, and that the evolute is the opposite process.
Proof of Equation (69).Note that Now apply equation (53).□ Next, we justify the Newton-Raphson Methods for forward kinematics.The method as usually states only applies to linear vector spaces, whereas we are working on the non-linear manifold of unit dual quaternions.Thus given η k , we need to define a map from the vector space of vector dual quaternions to unit dual quaternions close to η k .See, for example, [18].
Most papers on the Newton-Raphson Method on manifolds construct this map using the so called exp function [9,11,12].So the map is 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 11.Another exp map is to follow a one-parameter subgroup, or equivalently, equation (94).We do not take these approaches.Our approach is to normalize: (However, this is numerically close to the second approach, as is shown by equation (95).Using normalization instead of the exp map slightly reduces the complexity of the calculations as transcendental trigonometric functions are not required.But if one wants to use the exp map, simply replace (1 + θ) by exp(θ) throughout.) The main substantive difference between these methods, and the standard method that is used on linear vector spaces, is that the map from the linear space of vector dual quaternions to the manifold of unit dual quaternions changes with each iteration, since the map depends upon η k .But the theory that the Newton-Raphson Method converges with quadratic order is based upon examination of each iteration separately, so this shouldn't pose a great issue.
Justification of Equation (80).If F : R 6 → R 6 is a function for which we wish to solve for F(x) = 0, the method is to iterate Our approach, then, is to consider the map The prior guess is then θ k = 0. We have that The result follows.□ Justification of Equation (83).We seek to find the pose η so that L(η) is close as possible to ℓ.This is performed by minimizing the loss function b (140) The standard Newton-Raphson Method for optimizing the real valued quantity F (x), where x is an element of a vector space, is to iterate In our case x is θ, ) and our previous iterate is θ k = 0.The Jacobian is and the Hessian is H k , which by equation ( 65) is noting that

□
Now we work on proving Theorem 1.We use the Euler-Lagrange equations.(However, one could also use standard formulas for rotating bodies, and Newtonian physics, to obtain the same result.) Define the cross product of two vector dual quaternions α = a + ϵb and β = c + ϵd by Define the adjoint products of vector dual quaternions by which can also be defined by the property that for all vector dual quaternions α, β, and γ we have (Note that in the Lie algebra literature, the map β → α × β is often denoted ad α .Thus the map γ → α ⋊ γ is the formal dual of ad α .) Theorem 2. If the kinetic energy e satisfies equation (84), and v = v(η) denotes the potential energy, then the equation of motion is where ) and for any constant vector dual quaternion ψ we have that is, (156) The Lagrangian of the parallel robot is We use the local coordinate system θ given by equation ( 154).The Euler-Lagrange Equation [4,14] tells us We suppose that η 0 = η(t 0 ), and θ(t 0 ) = 0, and from now on in this proof, all equations are stated assuming the condition t = t 0 .Thus we only prove our results when t = t 0 .But since t 0 is arbitrary, this is not a limitation.However, it is important that derivatives are calculated before setting t = t 0 .In particular, this means that for any function f of η that ∂f ∂θ = Lf.For the parts coming from the kinetic energy, using linearity, it is sufficient to prove it for the additive parts of M. The part not involving m 0 is proved using Theorem 2, various vector identities, and remembering equation (36).So we only need to prove the kinetic energy portion in the case M = Λ T M 0 Λ.The easiest way to show this is to simply differentiate M 0 l = M 0 Λφ with respect to time.To do it directly from the formulas is more complicated, as we now show.For any constant vector dual quaternion ψ, we have where we used equation (125).Then it is simply a matter of collecting terms.□

Conclusion
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 (93), 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.