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

- Dynamics and Control Systems Laboratory, School of Aerospace Engineering, Georgia Institute of Technology, Atlanta, GA, United States

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.

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

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

Papadopoulos and 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. (1991, 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 spacecraft-mounted 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 spacecraft-mounted 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.

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.

## 2. Mathematical Preliminaries

### 2.1. 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}\triangleq \left\{-1,i,j,k:{i}^{2}={j}^{2}={k}^{2}=ijk=-1\right\}$. The algebra constructed from *Q*_{8} over the field of real numbers is the quaternion algebra, and it gives rise to the set ℍ. We define quaternions as $\mathbb{H}\triangleq \left\{q={q}_{0}+{q}_{1}i+{q}_{2}j+{q}_{3}k:{i}^{2}={j}^{2}={k}^{2}=ijk=-1,\text{}{q}_{0},{q}_{1},{q}_{2},{q}_{3}\in \mathbb{R}\right\}$. 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},\overline{q})$, where *q*_{0}∈ℝ and $\overline{q}={\left[{q}_{1},{q}_{2},{q}_{3}\right]}^{\text{T}}\in {\mathbb{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 ℝ^{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}_{\text{X}/\text{Y}}^{*}{q}_{\text{X}/\text{Y}}={q}_{\text{X}/\text{Y}}{q}_{\text{X}/\text{Y}}^{*}=1$, where 1 = (1, 0_{3 × 1}). This quaternion can be constructed as ${q}_{\text{X}/\text{Y}}=(cos(\varphi /2),\stackrel{\u0304}{n}sin(\theta /2))$, where $\stackrel{\u0304}{n}$ and θ are the *unit* Euler axis, and Euler angle of the rotation, respectively. It is worth emphasizing that ${q}_{\text{Y}/\text{X}}^{*}={q}_{\text{X}/\text{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, given ${\stackrel{\u0304}{s}}^{\text{X}}\in {\mathbb{R}}^{3}$, the coordinates of a vector expressed in frame X, its quaternion representation is given by ${s}^{\text{X}}=(0,{\stackrel{\u0304}{s}}^{\text{X}})\in {\mathbb{H}}^{v}$, where ℍ^{v} is the set of *vector* quaternions defined as ${\mathbb{H}}^{v}\triangleq \left\{({q}_{0},\overline{q})\in \mathbb{H}:\text{}{q}_{0}=0\right\}$ (see 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 ${s}^{\text{Y}}={q}_{\text{Y}/\text{X}}^{*}{s}^{\text{X}}{q}_{\text{Y}/\text{X}}$. Additionally, given *s*∈ℍ^{v}, we can define the operation [·]^{×}:ℍ^{v} → ℝ^{4 × 4} as

For quaternions $a=({a}_{0},\overline{a})$ and $b=({b}_{0},\overline{b})\in \mathbb{H}$, the left and right quaternion multiplication operators ${\left[\left[\xb7\right]\right]}_{\text{L}},{\left[\left[\xb7\right]\right]}_{\text{R}}:\mathbb{H}\to {\mathbb{R}}^{4\times 4}$ will be defined as

where

The three-dimensional attitude kinematics evolve as

where ${\omega}_{\text{X}/\text{Y}}^{\text{Z}}\triangleq (0,{\overline{\omega}}_{\text{X}/\text{Y}}^{\text{Z}})\in {\mathbb{H}}^{v}$ and ${\overline{\omega}}_{\text{X}/\text{Y}}^{\text{Z}}\in {\mathbb{R}}^{3}$ is the angular velocity of frame X with respect to frame Y expressed in Z-frame coordinates.

### 2.2. Dual Quaternions

The group of dual quaternion elements can be defined as

Dual quaternion algebra arises as the algebra of the dual quaternion group ℚ_{d} over the field of real numbers, and is denoted as ℍ_{d}. When dealing with the modeling of mechanical systems, it is convenient to define this algebra as ℍ_{d} = {** q** =

*q*

_{r}+ ϵ

*q*

_{d}:

*q*

_{r},

*q*

_{d}∈ ℍ}, 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 ℝ^{8}.

Analogous to the set of vector quaternions ℍ^{v}, we can define the set of vector dual quaternions as ${\mathbb{H}}_{d}^{v}\triangleq \left\{\text{q}={q}_{r}+\u03f5{q}_{d}:\text{}{q}_{r},{q}_{d}\in {\mathbb{H}}^{v}\right\}$. For vector dual quaternions we will define the skew-symmetric operator ${\left[\xb7\right]}^{\times}:{\mathbb{H}}_{d}^{v}\to {\mathbb{R}}^{8\times 8}$ as

For dual quaternions ** a** =

*a*

_{r}+ ϵ

*a*

_{d}and

**=**

*b**b*

_{r}+ ϵ

*b*

_{d}∈ ℍ

_{d}, the left and right dual quaternion multiplication operators ${\u301a\xb7\u301b}_{\text{L}},{\u301a\xb7\u301b}_{\text{R}}:{\mathbb{H}}_{d}\to {\mathbb{R}}^{8\times 8}$ are defined as

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 ∈ ℍ_{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** ∈ ℍ_{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 ${\mathit{q}}_{\text{B}/\text{I}}={q}_{\text{B}/\text{I},r}+\u03f5{q}_{\text{B}/\text{I},d}={q}_{\text{B}/\text{I}}+\u03f5{\scriptscriptstyle \frac{1}{2}}{q}_{\text{B}/\text{I}}{r}_{\text{B}/\text{I}}^{\text{B}}$, where ${r}_{\text{B}/\text{I}}^{\text{B}}$ 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,\stackrel{\u0304}{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 ${\text{q}}_{\text{Y}/\text{X}}^{*}={\text{q}}_{\text{X}/\text{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 ${\omega}_{\text{Y}/\text{Z}}^{\text{X}}=(0,{\stackrel{\u0304}{\omega}}_{\text{Y}/\text{Z}}^{\text{X}})$ and ${v}_{\text{Y}/\text{Z}}^{\text{X}}=(0,{\stackrel{\u0304}{v}}_{\text{Y}/\text{Z}}^{\text{X}})$, ${\stackrel{\u0304}{\omega}}_{\text{Y}/\text{Z}}^{\text{X}}$ and ${\stackrel{\u0304}{v}}_{\text{Y}/\text{Z}}^{\text{X}}\in {\mathbb{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}_{\text{X}/\text{Y}}^{\text{X}}=(0,{\stackrel{\u0304}{r}}_{\text{X}/\text{Y}}^{\text{X}})$, where ${\stackrel{\u0304}{r}}_{\text{X}/\text{Y}}^{\text{X}}\in {\mathbb{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 .

In general, the dual quaternion kinematics can be expressed as Filipe and Tsiotras (2013b)

#### 2.2.1. 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 ${\text{W}}^{\text{Z}}({O}_{p})\in {\mathbb{H}}_{d}^{v}$ expressed in Z-frame coordinates can be expressed in terms of its components as

where ${f}^{\text{Z}}=(0,{\stackrel{\u0304}{f}}^{\text{Z}}),\text{}{\tau}^{\text{Z}}=(0,{\stackrel{\u0304}{\tau}}^{\text{Z}})\in {\mathbb{H}}^{v}$ represent force and torque quaternions applied at point *O*_{p} as shown in Figure 1. Equivalently, we can describe the effect of *f*^{Z} and τ^{Z} about another point *O*_{q} 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}_{\text{p}/\text{q}}^{\text{Z}}$. Applying a frame transformation operation on a wrench about point *O*_{X} expressed in X-frame coordinates, given by ${\text{W}}^{\text{X}}({O}_{\text{X}})={f}^{\text{X}}+\u03f5{\tau}^{\text{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 ${\text{W}}^{\text{Y}}({O}_{\text{Y}})={f}^{\text{Y}}+\u03f5{\tau}^{\text{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, ${\text{W}}_{\text{ext}}^{\text{X}}({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 ${\text{s}}^{\text{X}}={s}_{r}^{\text{X}}+\u03f5{s}_{d}^{\text{X}}$ the adjoint transformation can be described by

Equivalently, given ${\text{s}}^{\text{Y}}={s}_{r}^{\text{Y}}+\u03f5{s}_{d}^{\text{Y}}$, the inverse adjoint transformation is described by

#### 2.2.2. Dual Inertia Matrix, Dual Momentum and 6-DOF Rigid Body Dynamics

The dual inertia matrix for a rigid body computed about its center of mass can be as follows Filipe and Tsiotras (2013b)

where *m*_{i} ∈ ℝ is the mass of the *i*-th body, 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 matrix-vector multiplication when the dual quaternion is represented as a vector in ℝ^{8}. The kinetic energy for a rigid body can be computed as

We can also define the matrix operator H(·):ℝ^{8 × 8} → ℝ^{8 × 8} to eliminate the swap operation in a matrix-dual quaternion multiplication. Applied on a matrix multiplying a dual quaternion **a** ∈ ℍ_{d}, we have that

In block form, this operator acts on *M* ∈ ℝ^{8 × 8}, composed of blocks ${M}_{1},{M}_{2}\in {\mathbb{R}}^{8\times 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}).

**LEMMA 3**. *The inverse of* H(*M*_{i}), H(*M*_{i})^{−1}, *exists and is given by*

**PROOF**. Through evaluation, .

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* as

or equivalently,

where is the net wrench applied on body *i* about its center of mass.

## 3. 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:

1. Revolute (R)

2. Prismatic (P)

3. Spherical (S)

4. Cylindrical (C)

5. Cartesian (U).

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

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

**Figure 2. (A)** Conceptual spacecraft architecture. **(B)** Architecture as rooted tree with labeled joints.

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

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

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 and exactly one “+1” entry, although, in general, columns can have several “−1” or “1” entries^{1}.

**EXAMPLE 1**. The incidence and branch termination matrices for the architecture shown in Figure 3 are given by

As example of the usage of the functions *N*(·) and *P*(·), we have

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 $\text{J}={\sum}_{i\text{}\in \text{}\mathrm{\text{Branches}}}{N}_{i}$. Using this notation, matrix *C* ∈ ℝ^{J × B} and vector *T* ∈ ℝ^{B}. We will define D as the total number of degrees of freedom added by the joints, which can be computed as $\text{D}={\sum}_{i\text{}\in \text{}\mathrm{\text{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 $\text{R}={\sum}_{i\text{}\in \text{}\mathrm{\text{Joints}}}{r}_{i}={\sum}_{i\text{}\in \text{}\mathrm{\text{Joints}}}6-{d}_{i}$.

The vector **y** ∈ ℝ^{8B} is defined as the collection of dual velocities, given by

The vector of generalized coordinates Γ ∈ ℝ^{D} represents the generalized coordinates of the joints and it is defined as

where the form of Γ_{Ji} 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} ∈ ℍ_{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 . Joint actuation wrenches ${\text{W}}_{\text{act},\text{i}}^{\text{i}}({O}_{\text{i}})$ induce motion about the degrees of freedom of the joint. Reaction wrenches ${\text{W}}_{\text{i}+1/\text{k}}^{\text{i}}({O}_{\text{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.

**Figure 4**. Body frame labeling and wrench definition at joint *J*_{i} between bodies *i*+1 and *k* = *N*(*C*(*i*, :)).

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}\in {\mathbb{R}}^{\text{R}}$, the collection of reduced reaction wrenches, as

where ${\stackrel{~}{\text{W}}}_{\text{i}+1/\text{N}(\text{C}(\text{i},:))}^{\text{i}}\in {\mathbb{R}}^{{r}_{i}}$ is obtained from ${\text{W}}_{\text{i}+1/\text{N}(\text{C}(\text{i},:))}^{\text{i}}\in {\mathbb{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 obtain ${\stackrel{~}{\text{W}}}_{\text{i}+1/\text{N}(\text{C}(\text{i},:))}^{\text{i}}$ from ${\text{W}}_{\text{i}+1/\text{N}(\text{C}(\text{i},:))}^{\text{i}}$ using the relationship ${\text{W}}_{\text{i}+1/\text{N}(\text{C}(\text{i},:))}^{\text{i}}={V}_{i}{\stackrel{~}{\text{W}}}_{\text{i}+1/\text{N}(\text{C}(\text{i},:))}^{\text{i}}$, the form of the matrix ${V}_{i}\in {\mathbb{R}}^{8\times {r}_{i}}$ depending on the type of joint. Table 4 lists the general wrench ${\text{W}}_{\text{i}+1/\text{N}(\text{C}(\text{i},:))}^{\text{i}}$, the reduced wrench ${\stackrel{~}{\text{W}}}_{\text{i}+1/\text{N}(\text{C}(\text{i},:))}^{\text{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. 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}.

### 3.2. 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 as

The joint dual velocity expressed in joint coordinates can be determined from

while the generalized coordinates of the joints can be determined to evolve as

The matrix *L*_{Ji} 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 by

resulting in

where we have used the fact that , by construction of Λ_{i}, defined in Table 4.

### 3.3. Dynamics

We will now generalize the rigid body Newton-Euler. We will show that the equations of motion can be cast in the form

We will define each of the blocks ${{S}}_{11}\in {\mathbb{R}}^{8\text{B}\times 8\text{B}}$, ${{S}}_{12}\in {\mathbb{R}}^{8\text{B}\times \text{R}}$, ${{S}}_{21}\in {\mathbb{R}}^{\text{R}\times 8\text{B}}$, ${{S}}_{22}\in {\mathbb{R}}^{\text{R}\times \text{R}}$, ${{B}}_{1}\in {\mathbb{R}}^{8\text{B}}$, and ${{B}}_{2}\in {\mathbb{R}}^{\text{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}\in {\mathbb{R}}^{\text{R}\times \text{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}\in {\mathbb{R}}^{8\text{B}\times \text{R}}$ couples the reaction wrenches with the dynamics of each body. These wrenches initially appear on the right-hand side of the Newton-Euler. 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}\in {\mathbb{R}}^{8\times {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}\in {\mathbb{R}}^{\text{R}\times 8\text{B}}$ introduces the dual accelerations of each body into the constraint equations. The matrix is composed of blocks ${({{S}}_{21})}_{ij}\in {\mathbb{R}}^{{r}_{i}\times 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

The form of matrix Λ_{i} depends on the type of joint and it is provided in Table 4.

The vector ${{B}}_{1}\in {\mathbb{R}}^{8\text{B}}$ corresponds to the right hand side of the Newton-Euler. In particular, it contains the non-linear 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}_{ext,i}^{{G}_{i}}({O}_{{G}_{i}})$. The vector is composed of sub-vectors ${({{B}}_{1})}_{i}\in {\mathbb{R}}^{8}$ given by

The vector ${{B}}_{2}\in {\mathbb{R}}^{\text{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}\in {\mathbb{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}_{\text{R}\times \text{R}}$, we can avoid inverting the large matrix on the left-hand-side of Equation (50) by using the Schur complement. Thus, if

we define the Schur complement of block ${{S}}_{11}$ as ${S}/{{S}}_{11}\triangleq -{{S}}_{21}{{S}}_{11}^{-1}{{S}}_{12}$. Therefore, the inverse of ${S}$ is given by

Hence, we can solve for the unknowns as

which upon expansion, yields

### 3.4. 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 follows

Assuming knowledge of the proximal body's dual acceleration , which must be solved for in tandem with all other dual accelerations and reaction wrenches, and since all velocity-level 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 is no longer an unknown, we must remove the corresponding equations from the system of equations presented in Equation (50). To do this, we remove from the vector of unknowns **ẏ**, and block-matrices ${({{S}}_{11})}_{\left\{:,i+1\right\}}$, ${({{S}}_{21})}_{\left\{:,i+1\right\}}$, which are the corresponding coefficients of that appear in both Newton-Euler. For the sake of exposition, let us rename these modified variables as $\widehat{\u1e8f}$, ${\widehat{{S}}}_{11}$, and ${\widehat{{S}}}_{21}$.

Next, we need to manipulate the modified Newton-Euler. In general terms, this equation is given by

and

where we have defined

By manipulating Equations (65) and (66), we obtain

and

Further manipulation of Equation (68) allows clearing ${\stackrel{~}{\text{W}}}_{\text{act},\text{i}}^{\text{i}}({O}_{\text{i}})$ of transformations as

where we have used Lemma 2 and the fact that ${V}_{\text{act},\text{i}}^{\text{T}}{V}_{\text{act},\text{i}}={I}_{{d}_{i}\times {d}_{i}}$ for ${\stackrel{~}{\text{W}}}_{\text{act},\text{i}}^{\text{i}}({O}_{\text{i}})\in {\mathbb{R}}^{{d}_{i}}$.

The resulting system of equations will be of the form

Here we have that

while ${{S}}_{\text{act},\text{i},1}\in {\mathbb{R}}^{(8(\text{B}-1)+{d}_{i})\times {d}_{i}}$ is described by

The vectors ${\widehat{{B}}}_{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 ${\mathbb{R}}^{(8(\text{B}-1)+R+{d}_{i})\times (8(\text{B}-1)+R+{d}_{i})}$ and thus, it is square and invertible.

### 3.5. Framework Summary

Algorithm 1 provides a detailed description of how to implement the kinematics and dynamics framework introduced in the previous sections.

## 4. Evaluation of Numerical Performance

We studied the performance of the algorithm on the satellite shown in Figure 3 without the end-effector. The inertias for the four different bodies were chosen as *M*_{1} = diag (1, 10, 10, 10, 1, 50, 50, 50) [kg, kg · m^{2}], *M*_{2} = diag (1, 5, 5, 5, 1, 2, 2, 1) [kg, kg · m^{2}], *M*_{3} = diag (1, 5, 5, 5, 1, 1, 2, 2) [kg, kg · m^{2}], and *M*_{4} = diag (1, 5, 5, 5, 1, 1, 2, 2) [kg, kg · m^{2}]. The geometry of the system was chosen as , , , , , , and , where the orientation of the frames for each of the bodies can be found in Chapter 5 of Valverde (2018).

The simulation was run using MATLAB R2017a's ODE45. The integrator's option AbsTol (absolute tolerance) was set to 1 × 10^{−14} and RelTol (relative tolerance) was set to 2.220 × 10^{−14}; the final time was set to *t*_{f} = 70 *s*. To evaluate center of mass, linear momentum, and angular momentum conservation, only internal (joint) wrenches were applied. The generalized forces on the wrenches

were set to

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.

**Figure 5. (A)** Center of mass deviation from initial position, **(B)** kinetic energy of the system, and **(C)** condition number of ${S}$.

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.

**Figure 6. (A)** Change in center of mass location, **(B)** linear, and **(C)** angular momentum comparison between decoupled and dual quaternion formulations.

## 5. 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 Administration, and funded through the internal Research and Technology Development program under award number SP.17.0004.008 RSA. This research was also performed under National Science Foundation award NRI1426945.

## Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

## Acknowledgments

The authors would like to thank David S. Bayard, from the G&C section at the Jet Propulsion Laboratory for valuable conversations that improved the content of this article.

## Footnotes

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.

## References

Bishop, B., Gargano, R., Sears, A., and Karpenko, M. (2015). Rapid maneuvering of multi-body dynamic systems with optimal motion compensation. *Acta Astronaut.* 117, 209–221. doi: 10.1016/j.actaastro.2015.07.035

Brodsky, V., and Shoham, M. (1999). Dual numbers representation of rigid body dynamics. *Mech. Mach. Theory* 34, 693–718. doi: 10.1016/S0094-114X(98)00049-4

Carignan, C. R., and Akin, D. L. (2000). The reaction stabilization of on-orbit robots. *IEEE Control Syst. Mag.* 20, 19–33. doi: 10.1109/37.887446

Dooley, J. R., and McCarthy, J. M. (1993). “On the geometric analysis of optimum trajectories for cooperating robots using dual quaternion coordinates,” in *Proceedings 1993 IEEE International Conference on Robotics and Automation*, Vol. 1 (Atlanta, GA), 1031–1036.

Dubanchet, V., Saussié, D., Alazard, D., Bérard, C., and Peuvédic, C. L. (2015). Modeling and control of a space robot for active debris removal. *CEAS Space J.* 7, 203–218. doi: 10.1007/s12567-015-0082-4

Dubowsky, S., Vance, E. E., and Torres, M. A. (1989). “The control of space manipulators subject to spacecraft attitude control saturation limits,” in *Proceedings of the NASA Conference on Space Telerobotics, Vol. 4* (Pasadena, CA), 409–418.

Featherstone, R., and Orin, D. (2000). “Robot dynamics: equations and algorithms,” in *Proceedings 2000 IEEE International Conference on Robotics and Automation*, Vol. 1 (San Francisco, CA), 826–834.

Filipe, N. (2014). *Nonlinear Pose Control and Estimation for Space Proximity Operations: An Approach Based on Dual Quaternions*. PhD thesis, Georgia Institute of Technology.

Filipe, N., Kontitsis, M., and Tsiotras, P. (2015). Extended Kalman filter for spacecraft pose estimation using dual quaternions. *J. Guidance Control Dyn.* 38, 1625–1641. doi: 10.2514/1.G000977

Filipe, N., and Tsiotras, P. (2013a). “Rigid body motion tracking without linear and angular velocity feedback using dual quaternions,” in *European Control Conference* (Zürich), 329–334.

Filipe, N., and Tsiotras, P. (2013b). “Simultaneous position and attitude control without linear and angular velocity feedback using dual quaternions,” in *Proceedings of the 2013 American Control Conference* (Washington, DC), 4815–4820.

Filipe, N., and Tsiotras, P. (2014). Adaptive position and attitude-tracking controller for satellite proximity operations using dual quaternions. *J. Guidance Control Dyn.* 38, 566–577. doi: 10.2514/1.G000054

Hooker, W. W. (1970). A set of r dynamical attitude equations for an arbitrary n-Body satellite having r rotational degrees of freedom. *AIAA J.* 8, 1205–1207. doi: 10.2514/3.5873

Jain, A. (1991). Unified formulation of dynamics for serial rigid multibody systems. *J. Guidance Control Dyn.* 14, 531–542. doi: 10.2514/3.20672

Jazar, R. N. (2010). *Theory of Applied Robotics: Kinematics, Dynamics, and Control*. New York, NY: Springer Science & Business Media.

Lee, J., Grey, M. X., Ha, S., Kunz, T., Jain, S., Ye, Y., et al. (2018). DART: dynamic animation and robotics toolkit. *J. Open Source Softw*. 3:500. doi: 10.21105/joss.00500

Longman, R. W., Lindbergt, R. E., and Zedd, M. F. (1987). Satellite-mounted robot manipulators - new kinematics and reaction moment compensation. *Int. J. Robot. Res.* 6, 87–103. doi: 10.1177/027836498700600306

Mohan, A., and Saha, S. K. (2007). A recursive, numerically stable, and efficient simulation algorithm for serial robots. *Multibody Syst. Dyn.* 17, 291–319. doi: 10.1007/s11044-007-9044-8

Moosavian, S. A. A., and Papadopoulos, E. (2004). Explicit dynamics of space free-flyers with multiple manipulators via spacemaple. *Adv. Robot.* 18, 223–244. doi: 10.1163/156855304322758033

Papadopoulos, E., and Dubowsky, S. (1990). “On the nature of control algorithms for space manipulators,” in *Proceedings 1990 IEEE International Conference on Robotics and Automation*, Vol. 2 (Cincinnati, OH), 1101–1108.

Papadopoulos, E., and Dubowsky, S. (1991a). “Coordinated manipulator/spacecraft motion control for space robotic systems,” in *Proceedings 1991 IEEE International Conference on Robotics and Automation* (Sacramento, CA), 1696–1701.

Papadopoulos, E., and Dubowsky, S. (1991b). On the nature of control algorithms for free-floating space manipulators. *IEEE Trans. Robot. Autom.* 7, 750–758. doi: 10.1109/70.105384

Rodriguez, G., Jain, A., and Kreutz-Delgado, K. (1991). A spatial operator algebra for manipulator modeling and control. *Int. J. Robot. Res.* 10, 371–381. doi: 10.1177/027836499101000406

Rodriguez, G., Jain, A., and Kreutz-Delgado, K. (1992). Spatial operator algebra for multibody system dynamics. *J. Astronaut. Sci.* 40, 27–50.

Saha, S. K. (1999). Dynamics of serial multibody systems using the decoupled natural orthogonal complement matrices. *J. Appl. Mech.* 66, 986–996. doi: 10.1115/1.2791809

Saha, S. K., Shah, S. V., and Nandihal, P. V. (2013). Evolution of the DeNOC-based dynamic modelling for multibody systems. *Mech. Sci.* 4, 1–20. doi: 10.5194/ms-4-1-2013

Seo, D. (2015). Fast adaptive pose tracking control for satellites via dual quaternion upon non-certainty equivalence principle. *Acta Astronaut.* 115, 32–39. doi: 10.1016/j.actaastro.2015.05.013

Sherman, M., and Rosenthal, D. (2001). *SD/FAST*. Available online at: https://support.ptc.com/support/sdfast/index.html

Sommer, H., Gilitschenski, I., Bloesch, M., Weiss, S., Siegwart, R., and Nieto, J. (2018). Why and how to avoid the flipped quaternion multiplication. *Aerospace* 5, 1–15. doi: 10.3390/aerospace5030072

Stoneking, E. (2007). “Newton-Euler dynamic equations of motion for a multi-body spacecraft,” in *AIAA Guidance, Navigation and Control Conference and Exhibit* (Hilton Head, SC).

Stoneking, E. (2013). “Implementation of Kane's method for a spacecraft composed of multiple rigid bodies,” in *AIAA Guidance, Navigation, and Control (GNC) Conference* (Boston, MA).

Umetani, Y., and Yoshida, K. (1989). Resolved motion rate control of space manipulators with generalized jacobian matrix. *IEEE Trans. Robot. Autom.* 5, 303–314. doi: 10.1109/70.34766

Valverde, A. (2018). *Dynamic Modeling and Control of Spacecraft Robotic Systems using Dual Quaternions*. PhD thesis, Georgia Institute of Technology.

Virgili-Llop, J. (2017). *SPART: Spacecraft Robotics Toolkit*. Available online at: https://github.com/NPS-SRL/SPART

Walker, M., and Weel, L.-B. (1991). “An adaptive control strategy for space based robot manipulators,” in *Proceedings of the IEEE International Conference on Robotics and Automation* (Sacramento, CA).

Xu, Y., and Shum, H.-Y. (1991). Dynamic control of a space robot system with no thrust jets controlled base. Technical Report CMU-RI-TR-91-33, Robotics Institute, Pittsburgh, PA.

Yoshida, K. (1999). “The SpaceDyn: a MATLAB toolbox for space and mobile robots,” in *Proceedings of the 1999 IEEE/RSJ International Conference on Intelligent Robots and Systems* (Osaka), 1633–1638.

## Appendix:

### Labeling Example On Two Arm Manipulator

We present an example of a spacecraft with two robotic arms, each with five links and three different types of joints. The example aims to further familiarize the reader with the proposed notation. The architecture of the satellite is shown in Figure 7. Figure 8 shows a schematic of the coordinate frames and the wrenches defined during implementation of the proposed framework. Reaction wrenches and actuation wrenches are assumed positive as applied on the body on which they are shown, and negative on the proximal body relative to the joint.

Keywords: multibody dynamics, dual quaternions, space operations, robotic servicing, Newton-Euler

Citation: Valverde A and Tsiotras P (2018) Dual Quaternion Framework for Modeling of Spacecraft-Mounted Multibody Robotic Systems. *Front. Robot. AI* 5:128. doi: 10.3389/frobt.2018.00128

Received: 22 March 2018; Accepted: 30 October 2018;

Published: 21 November 2018.

Edited by:

Evangelos G. Papadopoulos, National Technical University of Athens, GreeceReviewed by:

S. Ali A. Moosavian, K.N.Toosi University of Technology, IranVincent Dubanchet, Thales Alenia Space, France

Copyright © 2018 Valverde and Tsiotras. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Alfredo Valverde, avalverde3@gatech.edu