Dual Quaternion Framework for Modeling of Spacecraft-Mounted Multibody Robotic Systems

This paper lays out a framework to model the kinematics and dynamics of a rigid spacecraft-mounted multibody robotic system. The framework is based on dual quaternion algebra, which combines rotational and translational information in a compact representation. Based on a Newton-Euler formulation, the proposed framework sets up a system of equations in which the dual accelerations of each of the bodies and the reaction wrenches at the joints are the unknowns. Five different joint types are considered in this framework via simple changes in certain mapping matrices that correspond to the joint variables. This differs from previous approaches that require the addition of extra terms that are joint-type dependent, and which decouple the rotational and translational dynamics.


INTRODUCTION
The interest to operate servicing spacecraft in space has led to wide-ranging research in academia, governmental agencies, and private companies. The space servicing market is growing, and with it, also the need for effective and easy-to-use tools to model the different phases of the mission. One tool of particular interest that has garnered attention for proximity operations, during which not only the attitude, but also the position of a spacecraft has to be precisely controlled, are dual quaternions, see for example Filipe and Tsiotras (2013a), Seo (2015), and Filipe et al. (2015). We add to this body of literature, having as a goal to provide an intuitive development of the multibody dynamics of a spacecraft-mounted manipulator system in dual quaternion algebra using a Newton-Euler approach. The aim is to provide a unified mathematical framework in which to model the different phases of a servicing mission.

Multibody Dynamics for Space Applications
When it comes to mounting a robotic manipulator on a spacecraft, the development of the equations of motion is not as straightforward as ground-based robotic applications, due to the complex interaction between reaction forces that arise at the joints and the conservation of angular momentum. This is especially important for relatively small spacecraft with large manipulators, as the stationarity assumption of the base is not longer valid. In such scenarios (that become increasingly popular in practice) the combined base-manipulator motion has to be accounted for. In terms of prior work in spacecraft equipped with manipulators, we can mention Hooker (1970), who first derived the equations of motion for an n-body satellite. In his derivation, the reaction forces and torques at the joints are not explicit in the formulation and aims to expose the body axes so that it is convenient to incorporate control laws, internal forces and other disturbance forces into the model that would not be straightforward to introduce using a Lagrangian formulation. Hooker's approach is based on the addition of the independent equations of motion for each of the bodies to cancel the reaction forces, and the cancellation of reaction torques through a clever manipulation of the equations of motion. This leads to a system of equations where the unknowns are the angular acceleration of the base, and the generalized accelerations at the joints. Longman et al. (1987) proposed a model for the operation of a robotic arm mounted on the Space Shuttle when attitude control is enabled. The authors in Longman et al. (1987) developed a forward and inverse kinematic model based on an initial determination of where the center of mass of the system is. This allows for identification of where the satellite-body is in inertial space as a function of joint angles, enabling a custom-derived solution of the forward and inverse kinematics problem. The authors then provided an approach to extract the reaction forces and torques applied on the satellite base due to the robotic arm through the extension of results derived using a fixed-base approach. Umetani and Yoshida (1989) introduced the equations of motion for systems with revolute joints and an uncontrolled (not actively controlled) base. The authors introduce an innovative representation of the task-space motion of the end-effector using two Jacobians-one capturing the effect of the spacecraft's motion, and another one to capture the effect ofjoint rates. This framework, however, did not account for external forces or torques. Dubowsky et al. (1989) dealt with the problem of thruster, or joint actuator saturation as an integral part of path-planning for the manipulators. Their model consists of a nine-generalized-coordinate system, and they derive the equations of motion using the Lagrangian approach, involving a 9 × 9 × 9 tensor to compute the Coriolis-like term. Papadopoulos and Dubowsky (1991a) rewrite the equations of motion of the satellite-mounted robot arm, but this time include actuation of the satellite base, and embed them in a quasi-Lagrangian approach. Dubowsky (1990, 1991b) succinctly describe the equations of motion for a robotic arm on a satellite under the assumption of zero initial angular momentum using the Routhian and a compact representation of the kinetic energy of the system. The authors proceed to argue that fixed-base and space-based manipulators can almost always be controlled using the same control algorithms, given the structural similarities between the corresponding model matrices. Xu and Shum (1991) developed a dynamical model for a robotic arm mounted on a satellite base in the absence of thruster jets. This implies that the motion of the system obeys the conservation of linear and angular momentum, a fundamental fact in their derivation. Walker and Weel (1991) provided the equations of motion for a six degree of freedom robotic arm on a satellite base. The method incorporates three reaction wheels and the equations are derived using a Lagrangian formulation. They eliminate the velocity of the satellite base from this formulation, given the constraint of no external forces on the system, without necessarily assuming that the initial momenta of the system are zero. Their formulation leads to a complicated system of equations that relies on the pre-computation of a significant amount of partial derivatives. Carignan and Akin (2000) proposed a recursive Newton-Euler algorithm that is easy to implement, intuitive, and has been well adopted by the engineering community. As an example, Dubanchet et al. (2015) hinges on this dynamics framework to implement H ∞ control on a linearized version of the plant with the objective of designing a debris collection robotic manipulator in space. Stoneking (2007) uses a Newton-Euler approach which exposes the reaction forces of the system, solved for by a matrix inversion that also yields linear and angular accelerations. The same author proposes a decoupling of the equations for users not interested in the reaction forces at the joints. Furthermore, he provides a formulation for the case in which the joints are not given by a simple primitive (revolute or prismatic), such as is the case of a gimbaled joint. Bishop et al. (2014) used this method for path planning and control during rapid maneuvering of a robotic arm mounted on a spacecraft. Stoneking (2013) also proposed an approach based on Kane's equations of motion, in which the generalized coordinates appear as part of a minimal representation. In this case, extracting knowledge about the reaction forces and torques, which are particularly important during design phases, becomes a significantly more complicated task. Jain (1991) and Rodriguez et al. (1991Rodriguez et al. ( , 1992) provided a numerically efficient multibody dynamics framework based on spatial operator algebra. Featherstone and Orin (2000) and Featherstone (2008), provided generalizable algorithms to model multibody dynamics. In particular, in section 9.3 of Featherstone (2008) the author specializes his algorithm to free-floating bases. Saha (1999), Mohan and Saha (2007), and Saha et al. (2013) provided another numerical algorithm for recursive dynamics, which claims to be even more efficient than the one by Featherstone, and it relies on using projection matrices to eliminate reaction forces and torques. Software has also been developed to model general dynamical systems. For example, Moosavian and Papadopoulos (2004) describes SPACEMAPLE, a tool that uses an analytical formulation of the Lagrangian equations of motion. At Tohoku University in Japan, Yoshida (1999) and his research team developed the SpaceDyn toolbox. The toolbox uses a recursive Newton-Euler approach, method further explored by Carignan and Akin (2000). Other open source toolboxes available online include SPART by Virgili-Llop et al. (2017), developed specifically for spacecraftmounted manipulators, and DART by Lee et al. (2018), which is aimed for general multibody systems, among many others. Commercial software packages also exist. Among these, SD/FAST by Sherman and Rosenthal (2001) is a commonly used software package for spacecraft modeling.
In this wide literature for dynamic modeling of spacecraftmounted robotic manipulators, dual quaternions are mentioned and used only a handful of times. In particular, Dooley and McCarthy (1993) proposed using dual quaternions as generalized coordinates, and Brodsky and Shoham (1999) proposed a rigorous dual-number based methodology that resulted in a Lagrangian-like framework. Brodsky and Shoham did draw parallelisms with a Newton-Euler-type equation, but these were always projected onto the dual axes of motion for the cases concerning serial manipulators, obscuring any potential insight into the reaction forces and torques at the joints.

Operation Definition
Addition Multiplication by a scalar λa = (λa 0 , λa) The lack of previous work using dual quaternions in a classical Newton-Euler framework to model serial manipulator systems on a spacecraft motivated the work of this paper. In particular, the highly generalizable dynamics framework presented herein aims exploits the versatility of dual quaternions to capture coupled rotational and translational dynamic quantities, and to capture joint kinematic constraints at both, the velocity and the acceleration levels. The framework is developed using dual quaternions, an extension of the well-known quaternions, a mathematical language that is familiar to the practitioner in the field of spacecraft dynamics and control. Additionally, the proposed framework consists of a non-recursive approach that solves a well-defined system of equations for a satellite with a tree-like architecture. By providing a simple-to-follow algorithm, the proposed work aims at increasing the accessibility of the uninitiated into the realm of multibody dynamics.

Quaternions
The group of quaternions as defined by Hamilton in 1843 extends the well-known imaginary unit j, which satisfies the equation j 2 = −1. This non-abelian group is defined by the presentation Q 8 {−1, i, j, k : i 2 = j 2 = k 2 = ijk = −1}. The algebra constructed from Q 8 over the field of real numbers is the quaternion algebra, and it gives rise to the set H. We define quaternions as H {q = q 0 + q 1 i + q 2 j + q 3 k : i 2 = j 2 = k 2 = ijk = −1, q 0 , q 1 , q 2 , q 3 ∈ R}. This defines an associative, non-commutative, division algebra.
In practice, quaternions are often referred to by their scalar and vectors parts as q = (q 0 , q), where q 0 ∈ R and q = [q 1 , q 2 , q 3 ] T ∈ R 3 . The properties of quaternion algebra are summarized in Table 1. Previous literature has defined quaternion multiplication as the multiplication between a 4 × 4 matrix and a vector in R 4 , see for example Filipe and Tsiotras (2013b).
Since any rotation can be described by three parameters, the unit norm constraint is imposed on quaternions for attitude representation. Unit quaternions are closed under multiplication, but not under addition. A quaternion describing the orientation of frame X with respect to frame Y, q X/Y , satisfies q * X/Y q X/Y = q X/Y q * X/Y = 1, where 1 = (1, 0 3×1 ). This quaternion can be constructed as q X/Y = (cos(φ/2),n sin(θ/2)), wheren and θ are the unit Euler axis, and Euler angle of the rotation, respectively. It is worth emphasizing that q * Y/X = q X/Y , and that q X/Y and −q X/Y represent the same rotation. Furthermore, given quaternions q Y/X and q Z/Y , the quaternion describing the rotation from X to Z is given by q Z/X = q Y/X q Z/Y . This equation for composition of rotations corresponds to a Hamilton product convention as opposed to a Shuster convention, both of which are described in Sommer et al. (2018). Three-dimensional vectors can also be interpreted as special cases of quaternions. Specifically, givens X ∈ R 3 , the coordinates of a vector expressed in frame X, its quaternion representation is given by s Filipe, 2014 for further information). The change of the reference frame on a vector quaternion is achieved by the adjoint operation, and is given by where The three-dimensional attitude kinematics evolve aṡ where ω Z X/Y (0, ω Z X/Y ) ∈ H v and ω Z X/Y ∈ R 3 is the angular velocity of frame X with respect to frame Y expressed in Z-frame coordinates.

Dual Quaternions
The group of dual quaternion elements can be defined as Dual quaternion algebra arises as the algebra of the dual quaternion group Q d over the field of real numbers, and is denoted as H d . When dealing with the modeling of mechanical systems, it is convenient to define this algebra as

Operation Definition
Addition where ǫ is the dual unit. We call q r the real part, and q d the dual part of the dual quaternion q. Filipe and Tsiotras (2014) and Filipe (2014) have laid out much of the groundwork in terms of the notation and the basic properties of dual quaternions. The main properties of dual quaternion algebra are listed in Table 2. Filipe and Tsiotras (2013b) also conveniently define a multiplication between matrices and dual quaternions that resembles the well-known matrix-vector multiplication by simply representing the dual quaternion coefficients as a vector in R 8 .
Analogous to the set of vector quaternions H v , we can define the set of vector dual quaternions as H v d {q = q r + ǫq d : q r , q d ∈ H v }. For vector dual quaternions we will define the skew-symmetric operator For dual quaternions a = a r + ǫa d and b = b r + ǫb d ∈ H d , the left and right dual quaternion multiplication where The following lemma deals with the transformation invariance of the dual quaternion cross product operation.
LEMMA 1. The dual quaternion cross product is invariant to frame transformations. Specifically, PROOF. From the definition of the dual quaternion cross product given in Table 2, we have that The following lemma recasts the identity operation on a dual quaternion in terms of the left and right dual quaternion multiplication operations.
LEMMA 2. Given unit q ∈ H d , the left and right dual quaternion multiplication matrix operators satisfy the following identities: PROOF. To prove the first equality, let us apply the left-hand-side on the generic dual quaternion a ∈ H d and apply the definition of the multiplication matrix operators given in Equation (8) as and since qq * = qq * = 1, the result follows. The second equality is proven analogously.
Since rigid body motion has six degrees of freedom, a dual quaternion needs two constraints to parameterize it. The dual quaternion describing the relative pose of frame B relative to I is given by q B/I = q B/I,r + ǫq B/I,d = q B/I + ǫ 1 2 q B/I r B B/I , where r B B/I is the position quaternion describing the location of the origin of frame B relative to that of frame I, expressed in B-frame coordinates. It can be easily observed that q B/I,r · q B/I,r = 1 and q B/I,r · q B/I,d = 0, where 0 = (0,0), providing the two necessary constraints. Thus, a dual quaternion representing a pose transformation is a unit dual quaternion, since it satisfies q · q = q * q = 1, where 1 1 + ǫ0. Additionally, we also define 0 0 + ǫ0.
Similar to the standard quaternion relationships, two dual quaternion transformations can be composed to yield a third one via dual quaternion multiplication as q Z/X = q Y/X q Z/Y . Finally, the dual quaternion inverse is obtained using the conjugation operation, denoted as q * Y/X = q X/Y . A useful equation is the generalization of the velocity of a rigid body in dual form, which contains both the linear and angular velocity components. The dual velocity of the Y-frame with respect to the Z-frame, expressed in X-frame coordinates, is defined as where ω X Y/Z = (0,ω X Y/Z ) and v X Y/Z = (0,v X Y/Z ),ω X Y/Z andv X Y/Z ∈ R 3 are, respectively, the angular and linear velocity of the Y-frame with respect to the Z-frame expressed in X-frame coordinates, and r X X/Y = (0,r X X/Y ), wherer X X/Y ∈ R 3 is the position vector from the origin of the Y-frame to the origin of the X-frame expressed in X-frame coordinates. In particular, the dual velocity of a rigid body assigned to frame i with respect to the inertial frame, expressed in i -frame coordinates will be denoted ω i i /I = ω i i /I + ǫv i i /I . In general, the dual quaternion kinematics can be expressed as Filipe and Tsiotras (2013b)

Wrench and Twists and Their Transformations Using Dual Quaternions
In order to take full advantage of the potential of dual quaternions in the context of dynamic modeling of multibody systems, we have to specify how forces and torques are shifted from one frame to another. This will allow us, for example, to easily shift the application of a reaction force at a joint onto the center of mass of a given body, among other applications. A wrench W Z (O p ) ∈ H v d expressed in Z-frame coordinates can be expressed in terms of its components as where the extra torque term is due to the moment arm from point O q to point O p , captured by the position vector r Z p/q . Applying a frame transformation operation on a wrench about point O X FIGURE 1 | Wrench interpretation. expressed in X-frame coordinates, given by W X (O X ) = f X + ǫτ X , yields the following expression and by the definition of the cross product of two pure quaternion quantities given in Table 2, we get that The transformation described by Equation (18) implies that, given the dual force (e.g., force and torque) applied on a body at location O X , the equivalent wrench about a different location O Y can be computed by using a simple frame transformation operation, commonly known as the shifting law. As expected, the transformation changes the reference frame in which the original (X-frame) force and torque are being expressed, but it also adds a torque term that arises due to the lever of the force f X with respect to the new reference point O Y . Equivalently, the following transformation of W Y (O Y ) = f Y + ǫτ Y can be easily derived: Finally, when using wrenches, subscripts will be used to denote the source of, or a descriptor for, the wrench. For example, W X ext (O p ) denotes that the source of the wrench is "ext, " which for our case denotes an external force and torque applied at the end effector of the robotic arm. It is worth emphasizing that the wrench transformation can be used to merely change the orientation of the frame on which the wrench is expressed, or to simply translate the origin, without re-orientating the axes.
The frame transformation relationships we have just derived not only apply to wrenches, but also to twists. Therefore, given the twist s X = s X r + ǫs X d the adjoint transformation can be described by Frontiers in Robotics and AI | www.frontiersin.org Equivalently, given s Y = s Y r + ǫs Y d , the inverse adjoint transformation is described by where m i ∈ R is the mass of the i-th body,Ī i ∈ R 3×3 is the rigid body mass inertia matrix of the i-th body, and I 3×3 is the 3-by-3 identity matrix.
The dual momentum of body i computed about its center of mass and expressed in local frame i can be defined as where the ⋆ operator can be interpreted as conventional matrixvector multiplication when the dual quaternion is represented as a vector in R 8 . The kinetic energy for a rigid body can be computed as We can also define the matrix operator H (·) : R 8×8 → R 8×8 to eliminate the swap operation in a matrix-dual quaternion multiplication. Applied on a matrix multiplying a dual quaternion a ∈ H d , we have that In block form, this operator acts on M ∈ R 8×8 , composed of blocks M 1 , M 2 ∈ R 8×4 , as follows and, in particular, it acts on the dual inertia matrix M i as Therefore, we can also write the dual momentum as The following lemma deals with the invertibility of H M i .
For a multibody system S, with B rigid bodies whose centers of mass are located at i , Equation (23) can be generalized to yielding the dual momentum of the system computed about the origin of the inertial frame and expressed in I-frame coordinates. The kinetic energy of Equation (24) can be generalized as From Equation (23), we can compute the 6-DOF dynamic equations of motion of body i aṡ or equivalently, where W i i (O i ) = f i + ǫτ i is the net wrench applied on body i about its center of mass.

MULTIBODY SYSTEM MODELING USING DUAL QUATERNIONS
This section aims to provide a generalized dual quaternion framework to model the kinematics and the dynamics of a multibody system that contains joints of the following types: The approach is aimed toward characterizing spacecraft with one or more serial robotic arms having varying lengths. The framework, in fact, will hold for robotic arms that branch out themselves, while preserving a rooted tree structure, with the satellite base being the root.
As in previous sections, we will use roman variables for frames, subscripts and superscripts of physical quantities. We will use standard math font for the labeling of physical components, like bodies and joints. For example, body i will have its center of mass at i .

Variable Definition and Conventions
We will model the spacecraft as a graph G(v, e), where v is the number of vertices, and e represents the number of edges. This graph, in particular, will correspond to that of a directed and rooted tree with arborescent branching, that is, a graph with tree structure where direction of the edges matters, and these in general point away from the root.
For our specific application, the nodes of the graph will be the different rigid bodies composing the serial manipulator(s), and the edges will be the different joints of the manipulator(s). Figure 2A shows an example of the labeling for the different rigid bodies composing a two-arm configuration on a satellite. The same configuration is shown in Figure 2B with the labeling of the vertices (nodes) and edges accordingly.
As is common in graph theory, matrices will aid in the description of the system's topology. Two matrices will be particularly useful in this generalization: the incidence matrix, denoted by C, and the branch termination vector, denoted as T. The incidence matrix contains information about the connectivity between the joints and the bodies. The columns of the incidence matrix represent rigid bodies, while the rows represent joints. Thus, entry C ij indicates the relationship between joint i and rigid body j as follows if joint i is not connected to body j, −1, if joint i is distal, body j is proximal, where the relative positions are with respect to the satellite base.
The branch termination vector, T denotes whether the given body is the end of a branch. The body will most likely be an end-effector and external wrenches due to interaction with the environment may be applied on it. We define the vector T as (T) i = t i 1, if body i ends a branch, 0, otherwise.
We will define the functions N(·), P(·) as follows. Given a row or column of matrix C, or vector T, they output the indices of the "−1" entries, and the indices of the "+1" entries, respectively. Additionally, we will use the notation C(:, j) to identify the j-th column of C, C(i, :) to identify the i-th row of matrix C. It is worth emphasizing that each row will contain exactly one "−1" entry FIGURE 3 | Robotic arm configuration on a satellite base. and exactly one "+1" entry, although, in general, columns can have several "−1" or "1" entries 1 . Figure 3 are given by

EXAMPLE 1. The incidence and branch termination matrices for the architecture shown in
As example of the usage of the functions N(·) and P(·), we have P(C(1, :)) = P(row 1 of matrix C) = {2}, 1 The column corresponding to the satellite base will only have "−1" values, since no joint is proximal. The columns corresponding to end-effector bodies will only possess "+1" values since end-effectors are all distal with respect to their corresponding joint.
Let N i be the length of branch i, d i be the degrees of freedom of joint J i , J the total number of joints, and B the total number of rigid bodies. Therefore, B = 1 + J, and J = i ∈ Branches N i .
Using this notation, matrix C ∈ R J×B and vector T ∈ R B . We will define D as the total number of degrees of freedom added by the joints, which can be computed as D = i ∈ Joints d i .
Exploiting the duality between degrees of freedom at a joint, d i , and the dimensionality of the reaction wrench, r i , we will define R = i ∈ Joints The vector y ∈ R 8B is defined as the collection of dual velocities, given by The vector of generalized coordinates Ŵ ∈ R D represents the generalized coordinates of the joints and it is defined as where the form of Ŵ J i is dependent upon the type of joint J i . Table 3 lists the parametrization used for each type of joint. Here it is worth noting that the generalized coordinates parametrize the motion of the i frame (fixed to the distal body with respect to the joint) with respect to the proximal body, which is captured by the index k, where k = N(C(i, :)). In particular, S joints are modeled with an Eulerian 3-2-1 (yaw ψ, pitch θ , roll φ) rotation sequence for uniformity with other types of joints, even though these can be better parametrized by quaternions to avoid issues with singularities during the evaluation of the kinematics. Thus, the state vector for any given spacecraft-robotic arm configuration will be given by where q 1 /I ∈ H d is the pose of the base. Figure 4 shows joint J i with its associated frame i; the frame i+1 , which has the same orientation as frame i but its origin is at the center of mass of body i + 1; and the frame at the center of mass of the proximal body denoted by k , where k = N(C(i, :)). The origin of the i frame is positioned at the physical interface between the two adjoining bodies. Figure 4 also shows three types of wrenches. The reaction and actuation wrenches appear at the joint, with their point of application being the origin of the joint frame O i , and their coordinates expressed in the i frame. We additionally show the body wrench W i+1 i+1 (O i+1 ). Joint actuation wrenches W i act,i (O i ) induce motion about the degrees of freedom of the joint. Reaction wrenches W i i+1/k (O i ) arise due to physical constraints at the joints, and they are dual in nature to the joint actuation wrenches. Body wrenches, which are assumed to act at the center of mass of the body, come from control sources or other natural phenomena such as gravitational effects, or atmospheric drag, all appropriately transformed to the center of mass through the shifting law. For an example on the use of the naming convention for frames and wrenches, the reader is referred to the Appendix.
It will be assumed that the degrees of freedom of the joints are along the Z i -axis, which is a common assumption in the field of robotics, while the X i and Y i axes can be selected according to any predetermined set of rules, such as those laid out in Chapter 5 of Jazar (2010). The exceptions are the Cartesian and spherical joints, both of which have three degrees of freedom, and for which an orientation of the axes must be assumed a priori. For the cartesian joint, the local coordinate system is defined such that it is parallel to the physical axes of motion. For the spherical joint, one suggestion is to define the X i pointing toward the i+1th rigid body, while the Y i and Z i complete the orthogonal axis system.
We will define T ∈ R R , the collection of reduced reaction wrenches, as . . .
whereW i i+1/N(C(i,:)) ∈ R r i is obtained from W i i+1/N(C(i,:)) ∈ H d by eliminating the entries that correspond to the generalized coordinate of the joint, since there are no reaction forces or torques applied on the bodies about that generalized coordinate. In general, we can obtainW i i+1/N(C(i,:)) from W i i+1/N(C(i,:)) using the relationship W i i+1/N(C(i,:)) = V iW i i+1/N(C(i,:)) , the form of the matrix V i ∈ R 8×r i depending on the type of joint. Table 4 lists the general wrench W i i+1/N(C(i,:)) , the reduced wrench W i i+1/N(C(i,:)) , and the mapping matrix V i for each of the joints considered. The matrix E π (1,2,3,4,5,6,7,8;i) is formed by removing rows π(1, 2, 3, 4, 5, 6, 7, 8; i) from the eight-by-eight identity matrix, I 8×8 . The function π(·; i) selects an ordered subset of {1, 2, 3, 4, 5, 6, 7, 8} based on the type of joint i. The matrices i are provided for compactness, as they will be used in a future section as a way of eliminating a degree of freedom from a constraint equation for a given type of joint. Also, for completion purposes, we provide the form of the actuation wrenches in Table 5 and its corresponding mapping matrix from reduced actuation wrenches, identified by V act,i .

Kinematics
The kinematics of the system are fully characterized by the kinematics of the satellite base, and the kinematics of the joint generalized coordinates. The pose of the satellite base evolves aṡ q 1 /I = 1 2 q 1 /I ω 1 1 /I .
The joint dual velocity expressed in joint coordinates can be determined from while the generalized coordinates of the joints can be determined to evolve aṡ The matrix L J i depends on the joint type, and these are listed in Table 6.    Furthermore, from Equation (46), we can derive an acceleration-level relationship at each joint given bẏ

Joint type
resulting in where we have used the fact that iω i i/ k = 0, by construction of i , defined in Table 4.

Dynamics
We will now generalize the rigid body Newton-Euler dynamics to that of a spacecraft with multiple serial robotic manipulators. We will show that the equations of motion can be cast in the form We will define each of the blocks S 11 ∈ R 8B×8B , S 12 ∈ R 8B×R , S 21 ∈ R R×8B , S 22 ∈ R R×R , B 1 ∈ R 8B , and B 2 ∈ R R independently.
The block S 11 is composed of the dual inertia matrix for each of the bodies. It is given by Notice that since this matrix is block diagonal, its inverse can be easily computed as the inverse of its sub-blocks, which exist as proven in Lemma 3. Thus, in cases when there are no moving mechanical components, fluid slosh, or fuel consumption, the inverse of S 11 can be pre-computed and stored in memory to speed up computations. The block S 22 ∈ R R×R represents the effect of the reaction wrenches on the constraint equations. Since wrenches do not appear in the constraint equations, this block is composed of zeros. Explicitly, this block is given by The block S 12 ∈ R 8B×R couples the reaction wrenches with the dynamics of each body. These wrenches initially appear on the right-hand side of the Newton-Euler equation and are moved to the left-hand side as an unknown. The point of application of the wrench and the frame of reference are shifted to the center of mass of the body for which the equation is being derived. The matrix is composed of blocks of size (S 12 ) ij ∈ R 8×r j , corresponding to the attachment of body i to joint j, where each of these blocks is specified as The form of matrix V j depends on the type of joint as was detailed in Table 4. The block S 21 ∈ R R×8B introduces the dual accelerations of each body into the constraint equations. The matrix is composed of blocks (S 21 ) ij ∈ R r i ×8 , corresponding to the constraint at joint i and its relationship with body j as described by Equation (49). These sub-blocks are specified as [ [[q j /i The form of matrix i depends on the type of joint and it is provided in Table 4. The vector B 1 ∈ R 8B corresponds to the right hand side of the Newton-Euler equation. In particular, it contains the nonlinear term ω × (M ⋆ ω s ), the known wrenches applied at the center of mass, and the wrenches due to joint actuation. If the body ends a branch, it is assumed that it can interact with the environment at a specific point in the body. This is included in B 1 as well through "external" wrenches. External wrenches for branch i will be assumed to act at frame G i , the frame assigned to the end-effector of branch-terminating body i, and they will be denoted by W G i ext,i (O G i ). The vector is composed of sub-vectors (B 1 ) i ∈ R 8 given by (B 1 The vector B 2 ∈ R R corresponds to the right-hand-side of the constraint equations for each of the joints. In particular, it contains a cross term of dual velocities that arises when taking the derivative of the dual velocity constraint to yield a dual acceleration constraint, detailed in Equation (49). The vector is composed of sub-vectors (B 2 ) i ∈ R r i , given by where in the last equality we used the invariance of the dual quaternion cross product, proven in Lemma 1. Finally, since S 11 is always invertible and S 22 = 0 R×R , we can avoid inverting the large matrix on the left-hand-side of Equation (50) we define the Schur complement of block S 11 as S/S 11 −S 21 S −1 11 S 12 . Therefore, the inverse of S is given by 11 + S −1 11 S 12 (S/S 11 ) −1 S 21 S −1 11 −S −1 11 S 12 (S/S 11 ) −1 −(S/S 11 ) −1 S 21 S −1 11 (S/S 11 ) −1 .
Hence, we can solve for the unknowns as which upon expansion, yields

Locking or Prescribing Joint Motion
In some instances, it is desirable to lock a certain degree a freedom or prescribe its generalized coordinate, while still being able to determine the reaction wrenches produced by this motion. Additionally, knowledge of the required actuation wrench can provide insight into the holding torque that a given motor must provide, or exert during specific smaneuvers. A straight-forward modification of the equations provided herein can yield this information.
Let the admissible dual velocity and acceleration of the prescribed-motion for joint J i be given by The generalized speed is still mapped as followṡ Assuming knowledge of the proximal body's dual acceleratioṅ ω k k /I , which must be solved for in tandem with all other dual accelerations and reaction wrenches, and since all velocitylevel quantities are known, the distal body's dual velocity and acceleration are fully described by the kinematic relationships both of which can be easily derived from Equations (46) and (48). Since the dual accelerationω i+1 i+1 /I is no longer an unknown, we must remove the corresponding equations from the system of equations presented in Equation (50). To do this, we removė ω i+1 i+1 /I from the vector of unknownsẏ, and block-matrices (S 11 ) {:,i+1} , (S 21 ) {:,i+1} , which are the corresponding coefficients ofω i+1 i+1 /I that appear in both Newton-Euler, and constraint equations. For the sake of exposition, let us rename these modified variables asŷ,Ŝ 11 , andŜ 21 .
Next, we need to manipulate the modified Newton-Euler equation for bodies i + 1 and k = N(C(i, :)), since both are connected to joint J i , to include the actuation wrench as part of the vector of unknowns. In general terms, this equation is given by and where we have defined By manipulating Equations (65) and (66) Here we have that while S act,i,1 ∈ R (8(B−1)+d i )×d i is described by The vectorsB 1 and B 1 are identical, except the (i + 1)-th and k-th entries, which are computed as Additionally, the block diagonal matrix ϒ is described as It is worth emphasizing that the resulting matrix belongs to R (8(B−1)+R+d i )×(8(B−1)+R+d i ) and thus, it is square and invertible.

Framework Summary
Algorithm 1 provides a detailed description of how to implement the kinematics and dynamics framework introduced in the previous sections.
The deviation of the center of mass of the system with respect to its initial position is shown in Figure 5A. The total kinetic energy of the system is shown in Figure 5B, and the condition number for matrix S, described in Equation (57), is plotted in Figure 5C at every time step. The equations of motion were derived for the same architecture using classical Newton-Euler techniques as in the framework proposed by Stoneking Stoneking (2007), where the rotational dynamics are decoupled from the translational dynamics. The implementation required deriving the constraint equations for revolute joints, since these are not explicitly addressed by Stoneking. The numerical performance differences between the dual quaternion approach (DQ), and the decoupled formulation (Decoupled) of the dynamics, were evaluated for the same set of inputs. Figure 6A shows the comparison of the norm of the change of the center of mass of the system with respect to its initial position as a function of time. Next, the conservation of the linear and angular momenta of both systems is compared as shown in Figures 6B,C. As expected, the dual quaternion formulation possesses a numerical advantage since it more naturally accounts for the coupling between the rigid bodies' translational and rotational motion.

CONCLUSION
In this paper we have provided an intuitive approach to derive the dynamics of a satellite with a rooted-tree configuration with different joint types, including revolute, prismatic, spherical, cylindrical, and cartesian joints using dual quaternions. The approach exploits the structure of the Newton-Euler form of the dynamical equations of motion for a rigid body in dual quaternion form, allowing for the determination of the reaction wrenches at the joints. The different nature of the joints is captured by simple changes in the mapping matrices associated with each joint, and not through a fundamental change in the form of the equations -an advantage provided by the coupled nature of the kinematic and dynamic relationships expressed in terms of dual quaternions. The proposed framework can be particularly beneficial during proximity operations of a robotic servicing mission. Combining existing dual quaternion-based pose-tracking controllers with the proposed dual quaternion framework for the modeling of the multibody robotic servicer allows for the use of a unified algebra to model the different phases, ranging from navigation, to grappling and servicing.

AUTHOR CONTRIBUTIONS
The content of this article was developed as part of AV's Ph.D. thesis, in close collaboration with the PT.

FUNDING
This research was performed, in part, at the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space