Equations of Motion of Free-Floating Spacecraft-Manipulator Systems: An Engineer's Tutorial

The paper provides a step-by-step tutorial on the Generalized Jacobian Matrix (GJM) approach for modeling and simulation of spacecraft-manipulator systems. The General Jacobian Matrix approach describes the motion of the end-effector of an underactuated manipulator system solely by the manipulator joint rotations, with the attitude and position of the base-spacecraft resulting from the manipulator motion. The coupling of the manipulator motion with the base-spacecraft are thus expressed in a generalized inertia matrix and a GJM. The focus of the paper lies on the complete analytic derivation of the generalized equations of motion of a free-floating spacecraft-manipulator system. This includes symbolic analytic expressions for all inertia property matrices of the system, including their time derivatives and joint-angle derivatives, as well as an expression for the generalized Jacobian of a generic point on any link of the spacecraft-manipulator system. The kinematics structure of the spacecraft-manipulator system is described both in terms of direction-cosine matrices and unit quaternions. An additional important contribution of this paper is to propose a new and more detailed definition for the modes of maneuvering of a spacecraft-manipulator. In particular, the two commonly used categories free-flying and free-floating are expanded by the introduction of five categories, namely floating, rotation-floating, rotation-flying, translation-flying, and flying. A fully-symbolic and a partially-symbolic option for the implementation of a numerical simulation model based on the proposed analytic approach are introduced and exemplary simulation results for a planar four-link spacecraft-manipulator system and a spatial six-link spacecraft manipulator system are presented.


INTRODUCTION
Small spacecraft equipped with robotic manipulators will be the backbone of any robotic on-orbit servicing, asteroid mining, or active space debris removal missions. Examples for such missions are discussed in more detail by Yoshida (2003), Shoemaker and Wright (2004), Reintsema et al. (2010), Barnhart et al. (2013). The use of spacecraft-mounted robotic manipulator systems for the assembly and maintenance of spacecraft has a long and successful history with the Space Shuttle and International Space Station programs. The Space Shuttle orbiters were equipped with the Shuttle Remote Manipulator System, colloquially known as Canadarm (Sallaberger, 1997). The SRMS was successfully used to capture the Hubble space telescope and other satellites during servicing missions, to position astronauts during extra-vehicular activities, and to assemble and resupply the ISS (Goodman, 2006;Hale et al., 2011). The ISS is currently equipped with two manipulator systems, the Space Station Remote Manipulator System (Stieber et al., 1997), also known as Canadarm 2, and the Japanese Experiment Module Remote Manipulator System (Sato and Wakabayashi, 2001). The Space Station Remote Manipulator System is used to capture and berth H-II Transfer Vehicle (HTV), Dragon, and Cygnus vehicles, and to position supplies and astronauts (Dreyer, 2009;Bain, 2010;Ueda et al., 2010). With the Special Purpose Dexterous Manipulator (also called Dextre) as end-effector, the Space Station Remote Manipulator System is capable of fine manipulation (Coleshill et al., 2009). This capability is, for example, used in the NASA Robotic Refueling Mission, which demonstrated the feasibility of accessing and refueling a typical satellite fuel port (Cepollina, 2013). The Japanese Experiment Module Remote Manipulator System is mostly used to service the Japanese Experiment Module Kibo and is also equipped with a dexterous end-effector. The manipulator suite of the ISS is planned to be completed with the installation of the European Robotic Arm on the Russian module (Patten et al., 2002). The use of robotic manipulator arms mounted on small servicing spacecraft to capture and manipulate a client satellite has successfully been demonstrated with the Engineering Test Satellite VII (ETS-VII) and Orbital Express missions (Yoshida, 2003;Kennedy, 2008). The upcoming onorbit servicing demonstration missions Restore-L and Robotic Servicing of Geostationary Satellites will also feature robotic manipulators as core component of the capture and servicing system (Reed et al., 2016;Roesler et al., 2017).
Operating a robotic manipulator on a spacecraft results in a complex system overlapping the disciplines of robotics and aerospace engineering. Since the dynamics between the manipulator and the base-spacecraft are coupled, the system requires an integrated control system to meet the capture or manipulation goals and ensure mission success. Moreover, the complex dynamics of the spacecraft-manipulator system must be accounted for in the maneuver planning and in its overall maneuver timeline. The effects of manipulator operations on the orientation and position of the Shuttle orbiter and the ISS have been small and could be managed through operational procedures (Sargent, 1984). For fast-moving manipulators mounted on small base-spacecraft, the position and orientation disturbances due to manipulator motion become critical, as demonstrated on ETS-VII and Orbital Express (Oda, 2000;Kennedy, 2008). Hence, the engineers developing the spacecraft control system, specifying the sensor systems to be used, developing the communications system, and developing the operations plan must have an understanding of the complex dynamics arising from the multi-body system. However, multibody dynamics and modeling of robotic systems are not part of a typical undergraduate aerospace engineering curriculum and are rarely covered in aerospace engineering graduate programs. Therefore, aerospace engineers must resort to the academic literature to fill this knowledge gap. Available literature focusses on the application and performance analysis of various approaches to describing the coupled dynamics of a spacecraftmanipulator system (Moosavian and Papadopulos, 2007;Flores-Abad et al., 2014), but does not provide a complete description of the modeling approach that would allow an aerospace engineer to access the topic without consulting a combination of research publications and textbooks.
There are two common methods for modeling the dynamics and deriving the equations of motion of multibody systems: the recursive Newton-Euler method and the Lagrangian method. A general description of the Newton-Euler method is given in Siciliano et al. (2010). Stoneking (2007) provides a detailed presentation of the use of Newton-Euler dynamics to develop the equations of motion for a multibody spacecraft. Applications in the field of spacecraft-manipulator dynamics are shown in Longman et al. (1987) and Mukherjee and Nakamura (1992). In the Newton-Euler method, the equations of motion of the multibody system are computed from the equilibria of forces and torques acting on each link of the system. From this, a recursive algorithm can be developed. In the forward recursion through the structure of the multibody system, the linear and angular link accelerations and velocities are computed. The forces and moments acting on the links are then computed in the backward recursion. To develop the equations of motion of a system of flexible links, the Direct Path Method was developed (Ho, 1977;Hughes, 1979). In the Direct Path Method, the point of reference of the equations of motion is moved from the system centerof-mass to a fixed point in one of the bodies, which is typically selected to be the center-of-mass of the base spacecraft. The structure of the system is then described following the most direct path through the links. The torques acting on the links are taken about the joints instead of the link centers-of-mass, thus eliminating constraint forces and torques between the links.
The Lagrangian method develops the equations of motion of a multibody system from its kinetic and potential energies, using a set of generalized coordinates describing the positions of the links (Siciliano et al., 2010). Following Siciliano et al. (2010), the Lagrangian method is advantageous in it being systematic and easily comprehensible and in providing the equations of motion in a compact analytical form facilitating control systems design. The fundamental advantage of the Newton-Euler approach is its computational efficiency as a recursive algorithm. The Lagrangian method is thus used to explain the development of the equations of motion in the present paper, due to its systematic nature.
When controlling the motion of the spacecraft-manipulator system, the dynamic coupling of the base spacecraft and the manipulator becomes a concern. Since the base spacecraft in one of the five maneuvering modes described in section Detailed Classification of Spacecraft-Manipulator System Maneuvering, and thus not fixed in space, any motion of the manipulator will cause a rotation and translation of the base spacecraft. A comprehensive overview of methods to account for the dynamic coupling in controlling the position and orientation of both the end-effector of the manipulator and the base-spacecraft is provided in Flores-Abad et al. (2014). Three of these methods shall be highlighted here: The Virtual Manipulator (VM) approach, the Dynamically Equivalent Manipulator (DEM) approach, and the Generalized Jacobian Matrix (GJM) approach. The VM approach replaces the physical spacecraft-manipulator system with a dynamically consistent VM system (Vafa and Dubowsky, 1987). The base of the VM is a spherical joint located at the center-of-mass of the physical spacecraft-manipulator system. The orientation of this joint equals the orientation of the base-spacecraft in inertial space. Since, in the absence of external forces, the system center-of-mass remains stationary, the complex free-floating system is replaced by a dynamically consistent fixed-base system. The VM is a purely kinematic computational model in the sense that it is a massless kinematic chain with constant dimensions calculated from the geometry and mass properties of the spacecraft-manipulator system. It is typically used for workspace analysis for floating spacecraftmanipulator systems (Vafa and Dubowsky, 1987;Dubowsky and Papadopoulos, 1993). In particular, the VM cannot be represented by a physical manipulator. However, this can be achieved by the DEM approach (Liang et al., 1998). The base joint of the DEM is also a spherical joint at the system centerof-mass, and the DEM is geometrically identical to the VM for the same spacecraft-manipulator system. Differing from the VM, the DEM links have masses and moments-of-inertia matrices that are calculated from the mass distributions of the real system. Therefore, the DEM computation model can be reproduced in a real mechanical system, and thus be used in experimental systems.
The Generalized Jacobian approach was originally proposed by Yoshida and Umetani (1993), expanding the dynamic analysis previously introduced by Dubowsky and Papadopoulos (1993) and in the dissertation of Papadopoulos (1990). The Generalized Jacobian approach was used successfully in running simulations and developing control algorithms for the ETS-VII demonstrator mission (Yoshida, 2003). The approach was also further generalized to serve as a description for any under-actuated manipulator system (Yoshida and Nenchev, 1998). The Generalized Jacobian method formed the basis for the development of non-holonomic path planning algorithms (Nakamura and Mukherjee, 1989;Xu et al., 2008), target capture algorithms (Yoshida and Nakanishi, 2003), spacecraftmanipulator control strategies (Marchesi, 1997;Caccavale and Siciliano, 2001;Aghili, 2009), hardware-in-the-loop simulation of space robotic systems (Wei et al., 2006), Reaction Null-Space Control algorithms  and contact dynamics models . The method is advantageous for a symbolic analytic description of the dynamics of a complex spacecraft-manipulator system, which is important for the synthesis of guidance and control laws, as well as to guide the selection of proximity operations sensor systems and the design of pointing mechanisms for antennas, solar arrays, etc.
Notwithstanding the widespread use of the GJM approach, a complete description for the computation of all inertial parameters of the spacecraft-manipulator system is missing in the literature for the general case of a spatial N-link manipulator mounted on a six degrees-of-freedom (DOF) base-spacecraft. However, this is critical for its use in the design of the base spacecraft and the overall rendezvous and capture system and mission. This paper aims at filling this significant gap by presenting a complete symbolic, step-bystep derivation of the coupled equations of motion of a spacecraft-manipulator system following the GJM approach, including the time and joint angle derivatives of the generalized inertia matrix. The tutorial presented here uses the Lagrangian method for a single manipulator mounted on a base-spacecraft, under the assumption of zero linear and angular momenta, which makes it applicable to the description of the majority of current spacecraft-manipulator systems. The kinematics of the spacecraft-manipulator system are described both using direction cosine matrices and quaternions. For the extension of this approach to multi-manipulator systems, the reader is referred to Moosavian and Papadopoulos (2004), while the equations of motion of a spacecraft-manipulator system with non-zero angular momentum is covered in Nanos and Papadopoulos (2011). This complete discussion of the GJM approach enables the computation of symbolic expressions of the spacecraft-manipulator system equations of motion, which can then be used for the formulation of guidance and control laws in addition to numerical simulations.
A second important contribution of the paper is the introduction of a new, more accurate and detailed categorization for spacecraft-manipulator control modes. In particular, the two commonly used categories (free-flying and free-floating) are expanded by the introduction of five categories (namely floating, rotation-floating, rotation-flying, translation-flying, and flying).
A third contribution of the paper is the definition of the generalized Jacobian for an arbitrary point of the spacecraftmanipulator system, which is useful for guidance and control, obstacle avoidance, and collision detection and modeling.
To the authors' knowledge, this is the first time that the Generalized Jacobian approach is derived in completeness, drawing from multiple references in literature and the authors' own research. The paper thus forms a valuable resource for aerospace engineers to understand and model the complex dynamics of a robotic manipulator operated on a spacecraft.
The paper is organized as follows. Section Detailed Classification of Spacecraft-Manipulator System Maneuvering provides a detailed classification of the maneuvering modes of a spacecraft-manipulator system. Section Kinematics of a Spacecraft-Manipulator System develops expressions for the kinematics of a spacecraft-manipulator system based on a customized version of Denavit-Hartenberg parameters. Using these expressions, the dynamics of a floating spacecraftmanipulator system are described in section Dynamics of a Floating Spacecraft-Manipulator System. The paper proceeds to develop the generalized form of the equations of motion of a floating spacecraft-manipulator system from section Generalized Form of the Equations of Motion for a Floating Spacecraft-Manipulator System, followed by the generalized form of the generic Jacobians of a spacecraft-manipulator system in section Generalized Jacobian. Section Implementation of a Simulation Model for a Floating Spacecraft-Manipulator System then discusses the implementation options for a computer simulation model based on the generalized approach. The results of sample simulations are reported in section Numerical Simulations. Finally, section Conclusion provides concluding remarks. a Apart from the effect of external contact forces (e.g., with the target object). b The rotation of the base-spacecraft is controlled only by momentum-exchange devices (reaction wheels or control-moment gyroscopes) (Hughes, 2004). c The rotation of the base-spacecraft is controlled only by external torques (typically provided by reaction-jet thrusters).

DETAILED CLASSIFICATION OF SPACECRAFT-MANIPULATOR SYSTEM MANEUVERING
The spacecraft-manipulator system can maneuver in different modes, typically designated by the terms free-floating and free-flying (Umetani and Yoshida, 1989;Dubowsky and Papadopoulos, 1993). To arrive at a more detailed and complete classification of spacecraft maneuvering modes, the authors propose to add three modes, thus fully covering all possible spacecraft maneuvers (see Table 1). The new modes are defined for an isolated spacecraft-manipulator system operating in pure weightlessness and in the absence of friction.

Conserved Linear and Angular Momentum
The maneuvering cases here defined as floating and rotationfloating have in previous literature been referred to as freefloating. In both cases the total linear momentum and the total angular momentum of the spacecraft-manipulator system are not subject to any external control forces or torques and are thus conserved.

Floating
A spacecraft-manipulator system is here defined to be floating when maneuvering in a six DOF under-actuated mode in which only the manipulator joints are actively controlled. The system moves only under the effect of the internal reactions due to the actuation of the manipulator's joint motors. As robotic manipulators mounted on spacecraft typically only use revolute joints, these internal reactions are typically torques. Therefore, the 3 DOF of orientation of the base-spacecraft and the 3 DOF of translation of the system's center-of-mass are not actively controlled.

Rotation-Floating
A spacecraft-manipulator system is here defined to be rotationfloating when maneuvering in a 3 DOF under-actuated mode, in which both the DOF at the manipulator's joints and the three DOF of orientation of the base-spacecraft are controlled by internal torques only. The rotation-floating maneuver case thus differs from floating in that the attitude of the basespacecraft is actively controlled by momentum-exchange devices [i.e., reaction wheels or control-moment gyroscopes (Hughes, 2004)]. The three DOF of translation of the system's center-ofmass are not actively controlled.

Controlled Linear and Angular Momentum
All of the three maneuvering cases here defined as rotationflying, translation-flying, and flying have been defined in previous literature as free-flying. In all three cases, the total angular momentum, total linear momentum, or both are actively controlled by external forces and/or torques, e.g., by a reaction control system.

Rotation-Flying
A spacecraft-manipulator system is here defined to be rotationflying when the DOF at the manipulator joints are actively controlled by joint motor torques and the three DOF of orientation of the base-spacecraft are actively controlled by external torques only. This is typically achieved by reaction-jet thrusters firing in couples, thus generating a pure torque with total null force. The three DOF of translation of the system center-of-mass are not actively controlled. Therefore, the system's total linear momentum is in this case constant while the angular momentum is time-varying.

Translation-Flying
A spacecraft-manipulator system is here defined to be translation-flying when both the DOF at the manipulator joints and the three DOF of translation of the base-spacecraft are actively controlled. The manipulator DOF are controlled by joint motor torques. The base-spacecraft translation is controlled by external forces, provided typically by reaction-jet thrusters. Furthermore, the three DOF of orientation of the base-spacecraft are either not actively controlled or controlled only by momentum-exchange devices. Therefore, the system's total angular momentum is in this case constant while the linear momentum is time-varying.

Flying
A spacecraft-manipulator system is here defined to be flying when maneuvering in a mode in which all of the DOF at the manipulator joints are actively controlled by joint motor torques and the six DOF of motion of the base-spacecraft are actively controlled by external forces, provided typically by reaction jet thrusters. In this case, both the system's total angular momentum and the linear momentum are time-varying. The classification above is rigorously valid only for an isolated spacecraft-manipulator system. In reality, a spacecraftmanipulator system is never isolated but orbiting an extended body (e.g., the Earth) under its gravitational attraction. However, FIGURE 1 | Geometry of a spacecraft-manipulator system. the classification above can still be used in an approximate sense, due to the weightless (e.g., free-falling) condition of the system center-of-mass (due to the balancing of gravitational attraction and centrifugal force) and neglecting the effect of environmental torques (typically dominated by gravity-gradient torque, atmospheric torque, and solar radiation-pressure torque ) and non-gravitational environmental forces (typically dominated by atmospheric drag and solar radiation-pressure). The analysis and simulation of floating, rotation-floating and rotation-flying maneuvering modes can be typically conducted with good accuracy as if the system was isolated. In those three cases, a coordinate system centered at the center-of-mass of the orbiting spacecraft-manipulator system and having axes oriented in a fixed way with respect to an inertial frame (i.e., having zero absolute angular velocity) can be considered as equivalent to an inertial coordinate system for the description of the spacecraftmanipulator system motion.

KINEMATICS OF A SPACECRAFT-MANIPULATOR SYSTEM
The geometry of a base-spacecraft equipped with a single Nlink manipulator is illustrated in Figure 1. The base-spacecraft is defined as link 0, with joint 0 residing at the base-spacecraft center-of-mass. The end-effector (EE) of the manipulator is part of link N and located at the end of that link. Therefore, its location can be treated as a virtual joint N+1. The geometry of a spacecraft-manipulator system with respect to the origin of an inertial reference frame J can be described unequivocally by the position vectors of the joints, p i (0 ≤ i ≤ N +

DH parameter
Geometric meaning Distance along the common normal between z J i and z J i+1 1), and by the position vectors of the centers-of-mass of the links, r j (0 ≤ j ≤ N). In this paper, the word "vector" is rigorously reserved for Gibbsian vectors (see, for instance, Hughes, 2004). Three types of Cartesian coordinate systems are used in the analysis: an inertial coordinate system J , a set of N+2 jointfixed coordinate systems The base-spacecraft coordinate system J 0 ≡ L 0 has its origin at the center-of-mass of the base-spacecraft, and has an arbitrary orientation (for instance, corresponding to its principal directions of inertia). Each joint coordinate system J i (1 ≤ i ≤ N + 1) is built using of the Denavit-Hartenberg (DH) convention as specified in Siciliano et al. (2010). The axis z J i of J i is parallel to the joint rotation axisk i for 1 ≤ i ≤ N. The origin ̺ 1 of J 1 as well as the directions of x J 1 and y J 1 can be selected as needed for a particular geometry. The same is true for the end-effector coordinate system J N+1 . Furthermore, for 2 ≤ i ≤ N, the origin ̺ i of J i is at the intersection of z J i and the common normal between z J i−1 and z J i . Whenever z J i−1 and z J i are parallel, the common normal is uniquely defined and ̺ i is selected for convenience (e.g., passing through the centerof-mass of the link i-1). The axis x J i points from ̺ i along the direction of the common normal defined above. The axis y J i completes the right-handed Cartesian coordinate system. The link coordinate system L i is parallel to J i+1 (for 1 ≤ i ≤ N), but with the origin at the center-of-mass (CM i ) of link i. This way of defining the joint coordinate systems follows Siciliano with the following difference: Siciliano's link i coordinate system is identical to our joint i+1 coordinate system J i+1 . Furthermore, the origins of the link coordinate systems L i are placed into the link centers-of-mass to comply with Yoshida's approach (Yoshida and Umetani, 1993).
The geometrical relationship between subsequent joint coordinate systems J i and J i+1 (1 ≤ i ≤ N) is described by the four DH parameters d i , θ i , α i , and c i (see Table 2 and Figure 2). For a link with a rotary joint, θ i is the only variable parameter and is identical to the joint angle q i . The parameters d i , α i , and c i represent fixed geometric properties of the manipulator link. For a prismatic joint, the parameter d i is variable and identical to the joint extension. Since prismatic joints are not commonly used in orbital robotics applications, only rotary joints are considered in this paper.
The homogenous transformation matrix J i T J i+1 from mathcalJ i to J i+1 coordinates can be expressed as a function of the DH parameters as follows: where the DH transformation matrix function A (β, u, γ , w) is defined as: Therefore, the position and orientation of each joint coordinate system J i (2 ≤ i ≤ N + 1) can be expressed recursively by a product of DH transformation matrices. The resulting homogenous transformation matrix J T J i contains the direction-cosine matrix (DCM) J R J i and the position vector J p i : Analogously, the homogenous transformation matrix J T L i contains the DCM J R L i and position J r i . Notably, because of the fact that L i is oriented as J i+1 , it yields: The joint DCM are used to derive the directionsk i of the joint rotation axes in the inertial frame. In the DH convention, the rotation axis of a rotary joint is defined as the z axis in the corresponding link frame. So,k i is defined by: Furthermore, it is useful to define the following additional vectors. The position vector a i connects the origin ̺ i of the joint i coordinate system J i to the center-of-mass of link i (CM i , origin of link i coordinate system L i ). The position vector b i connects CM i and ̺ i+1 . Notably, if CM i lies on a straight line connecting ̺ i and ̺ i+1 , which is a valid assumption for most common spacecraft manipulators, a i and b i can be substituted by their scalar components a i and b i on the x axis of L i and J i+1 . Notably, since the base-spacecraft is designated as link 0 and its center-of-mass is joint 0, it is r 0 = p 0 , L 0 = J 0 , J R L 0 = J R J 0 , and J T L 0 = J T J 0 . In particular, we have chosen to express the DCM for the base-spacecraft in Euler angles (roll angle φ, pitch angle θ , yaw angle ψ) in the following 1-2-3 rotation sequence: Therefore: cos ψ cos θ cos ψ sin θ sin ψ − sin ψ cos ψ cos ψ sin θ cos ψ + sin ψ sin ψ sin ψ cos θ sin ψ sin θ sin ψ + cos ψ cos ψ sin ψ sin θ cos ψ − cos ψ sin ψ I r 0 −sin θ cos θ sin ψ cos θ cos ψ The joint coordinate system J 1 has a constant position 0 b 0 and orientation L 0 R J 1 , with respect to the base-spacecraft coordinate system. By expressing L 0 R J 1 in Euler angles with 1-2-3 sequence (ǫ, ζ , and η), it yields: cos η cos ζ cos η sin ζ sin ǫ − sin η cos ǫ cos η sin ζ cos ǫ + sin η sin ǫ sin η cos ζ sin η sin ζ sin ǫ + cos η cos ǫ sin η sin ζ cos ǫ − cos η sin ǫ 0 b 0 −sin ζ cos ζ sin ǫ cos ζ cos ǫ Frontiers in Robotics and AI | www.frontiersin.org For all of the subsequent joint coordinate systems J i , by taking into account Equation (1), it yields: where the transformation matrices for the link coordinate systems, I T L i , are derived as: The relative position between the base-spacecraft center-of-mass and the link centers-of-mass, J r 0i , and the position of the system center-of-mass, J r C , can be expressed as: where m tot = N i = 0 m i . Furthermore, the position of the system center-of-mass with respect to the base-spacecraft center-of-mass is given by: Equations (12) and (13) yield: Given the geometry of the spacecraft-manipulator system, the linear velocity J v X i of any point X i on the i-th link of the manipulator and the angular velocity J ω i of the i-th link can be expressed as: The notation a × is used to denote the skew-symmetric matrix (Hughes, 2004) associated with any vector a, so that the vector product c = a × b can be expressed in the matrix-vector operation c = a × b: Alternative Expression of Transformations Using the Unit Quaternion Instead of the DCM The use of DCM in developing the transformation matrices for the spacecraft-manipulator system is familiar and intuitive. DCM, or equivalently rotation matrices, offer a unique and singularity-free parameterization of the orientation. As such, the DCM method is used above to describe the DH formalism, and DCM-based notation is used throughout this paper. However, the parameterization of the rotation using three parameters, e.g., Euler angle combinations, creates singularities, thus making the resulting DCM non-invertible. In addition, the multiplication of transformation matrices results in a large number of operations, 27 for two [3 × 3] rotation matrices and 64 for two [4 × 4] homogenous transformation matrices, thus making the approach computationally inefficient. The alternative is to use Euler symmetric parameters, also called unit quaternions (Wertz, 1978;Hughes, 2004). Unit quaternions are used to describe rotations following the Euler axis-angle theorem, with the vector part ε of the quaternion describing the rotation axis, and the scalar part η describing the rotation angle. A unit quaternion q expressing a rotation about any axisû by any angle α can be written as a four-vector: The magnitude of q is q = 1, hence q is called a unit quaternion. The unit quaternion thus corresponds to a rotation matrix R (Hughes, 2004): A sequence of rotations, R 3 = R 1 R 2 , is then expressed by the quaternion product (Campa and Camarillo, 2008): Therefore, any vector a, expressed as four-vector a = a 1 a 2 a 3 0 T , can be rotated using a unit quaternion (Kuipers, 2000): with the conjugate quaternion q * = −ε T η T . Therefore, the application of the DH convention results in an unit quaternion J i q rJ i+1 for the rotation between coordinate systems J i and J i+1 , and a four-vector J i t J i+1 for the translation between origins of J i and J i+1 (Campa and Camarillo, 2008): The rotation from the base joint J 1 to the end-effector J N+1 can thus be computed from a sequence of rotations using Equation (20): The sequence of translations results from Equation (21): 1 t E = 1 t N+1 = 1 t 2 + + 1 q r 2 ⊗ 2 t 3 ⊗ 1 q * r2 + +( 1 q r 2 ⊗ 2 q r 3 ) ⊗ 3 t 4 ⊗ 1 q r 2 ⊗ 2 q r 3 * + . . . +

DYNAMICS OF A FLOATING SPACECRAFT-MANIPULATOR SYSTEM
The equations of motion of the spacecraft-manipulator system illustrated in Figure 1 are here derived using the Lagrangian approach for the case of floating maneuvering mode (see section Detailed Classification of Spacecraft-Manipulator System Maneuvering and Table 1). In this case, the potential energy of the system is zero, and the Lagrangian equals the kinetic energy T: By collecting the mass and inertia properties of the spacecraftmanipulator system into the [6 × 6] base-spacecraft inertia matrix H 0 , the [N × N]manipulator inertia matrix H m , and the [6 × N] dynamic-coupling inertia matrix H 0m , after some algebraic steps the kinetic energy can be expressed as (Yoshida and Umetani, 1993): where Jẋ 0 = J v 0 T J ω T 0 T is the combined and angular velocity matrix of the base-spacecraft. The base-spacecraft inertia matrix, manipulator inertia matrix, and the dynamiccoupling inertia matrix are thoroughly explained in the following subsections.

Base-Spacecraft Inertia Matrix
The [6 × 6] base-spacecraft inertia matrix H 0 is resulting to be: where I 3,3 is the [3 × 3] identity matrix and the [3 × 3] submatrix H S collects the moments-of-inertia of the spacecraft-manipulator system about the base-spacecraft center-of-mass, expressed in the inertial frame J : Frontiers in Robotics and AI | www.frontiersin.org The i-th link moments-of-inertia matrix in the inertial coordinate system is derived from the moments-of-inertia matrix in the i-th link coordinate system as:

Dynamic-Coupling Inertia Matrix
The [6 × N] dynamic-coupling inertia matrix H 0m expresses the contribution of the dynamic coupling between the manipulator and the base-spacecraft to the kinetic energy of the spacecraftmanipulator system. In particular, it is: where the [3 × N] submatrix J TS collects the contribution to the system kinetic energy of the combination of the effect of the manipulator joint rateq and the base-spacecraft linear velocity J v 0 . In detail, it is: where the [3 × N] matrix J Ti represents the linear velocity Jacobian for the center-of-mass of the i-th link (Siciliano et al., 2010), and is given by: Finally, the [3 × N] submatrix H Sq contains the contribution to the system kinetic energy of the combination of the effect of the manipulator joint rateq and the base-spacecraft angular velocity J ω 0 . In detail, it is: where the [3 × N] matrix J Ri represents the angular velocity Jacobian for the i-th link (Siciliano et al., 2010):

Manipulator Inertia Matrix
The [N × N] inertia matrix for the manipulator, H m , is identical to that of any fixed-base manipulator (Siciliano et al., 2010) and expressed as:

Equations of Motion of a Floating Spacecraft-Manipulator System
The generalized coordinates for the description of the spacecraftmanipulator system were chosen to be the manipulator joint angles q and the base-spacecraft combined linear and angular position J x 0 , resulting in the following two matricial Lagrangian equations. In the following, the coordinate system superscript will be omitted for better readability: where the [N × 1] matrix τ contains the manipulator joint torques. It would be possible to include the effects of the internal torques due to the presence of momentum-exchange devices on the right-hand side of Equation (37), following the procedure outlined in Wie (2008). That would extend the analysis to the case of rotation-floating (see Table 1). By substituting the kinetic energy expressed in Equation (27) into Equations (37) and (38), by computing the derivatives, and by properly rearranging the terms, the matricial equations of motion (corresponding to N+6 scalar equations) for the floating spacecraft-manipulator system result as: where the [6 × 1] matrix c 0 and the [N × 1] matrix c m are defined as:

GENERALIZED FORM OF THE EQUATIONS OF MOTION FOR A FLOATING SPACECRAFT-MANIPULATOR SYSTEM
The equations of motion in Equation (39) govern the dynamics of a floating spacecraft-manipulator system in terms of coupled base-spacecraft motion (first six scalar equations) and manipulator motion (last N scalar equations). In this section, the N+6 scalar equations of motion are rewritten into a set of N scalar generalized equations, as typically done for under-actuated systems (Yoshida and Nenchev, 1998). The linear momentum, P, and the angular momentum, L, are given by (Umetani and Yoshida, 1989): where M 0 = P L T is the [61] combined matrix of initial linear and angular momentum about the base-spacecraft center-ofmass, expressed in the inertial frame J . Because of the hypothesis of floating maneuvering in the absence of any external forces or moments, the location of the system center-of-mass remains constant, the momenta are conserved, and therefore it yields: d dt Equation (42) is solved forẋ 0 . The result is inserted into Equation (43), which is then solved forẍ 0 . The resulting expressions are finally inserted into the last N scalar equations of Equation (39). By considering the symmetry of H −1 0 (which will be discussed below), the resulting matricial equation (corresponding to N scalar equations) is: where the [N × N] generalized inertia matrix H ⋆ is defined to be: The further treatment of this generalized equation of motion in the presence on non-zero angular momentum is described in Nanos and Papadopoulos (2011). For the purposes of this tutorial, we assume from now on that the initial momentum of the floating spacecraft-manipulator system is zero, i.e.: M 0 = 0 6,1 .
Therefore, Equation (44) simplifies to: where the [N × 1] matrix c ⋆ is defined as: The velocity and position dependent terms in the generalized equations of motion can be combined into a single [N × 1] matrix C ⋆ : This permits the familiar representation of the generalized equations of motion which closely resembles the equations of motion for a fixed-base manipulator (Siciliano et al., 2010): The thorough symbolic expression of the terms of Equation (47) is provided in the following subsections.

Inverse Base-Spacecraft Inertia Matrix
The calculation of the symbolic expression of H ⋆ from Equation (45) By following the approach taken by Mukherjee and Nakamura (1992), the inverse of H 0 is determined using the Banachiewicz inversion formula (Baksalary and Styan, 2002) and the Schur complement S U of the non-singular submatrix U: where the Schur complement of U for a matrix U V W X is defined as Baksalary and Styan (2002): In detail, S U is a symmetric [3 × 3] matrix with the elements:

Time Derivative of the Generalized Inertia Matrix
The time derivative of the generalized inertia matrix,Ḣ ⋆ , is calculated by taking the time derivative of Equation (45), resulting in:Ḣ The derivative of a matrix inverse with respect to a scalar can be expressed as: as it can be immediately obtained from where Q is any square-invertible matrix and x is a scalar. This allows computing the time derivative of the generalized inertia matrix without the need for calculating the time derivative of the inverse of the base-spacecraft inertia matrix, H −1 0 : The following subsections discuss the calculation of the individual matrix derivatives in Equation (65).

Time Derivative of the Manipulator Inertia Matrix
The only new values to be calculated are the time derivatives of the angular and linear motion Jacobians,J Ri andJ Ti , as well as the time derivative of the link moments-of-inertia matrices, . For the calculation of these terms, the time derivatives of the joint and link direction-cosine matrices, JṘ J i and JṘ L i are required. The time derivative of any joint directioncosine matrix J R J i is given by the Darboux equation (Hughes, 2004): The angular velocity of any joint coordinate system with respect to the inertial frame can be computed recursively from: Due to the definition of the joint and link coordinate systems in section Kinematics of a Spacecraft-Manipulator System: Furthermore, the time derivative of Equation (30) yields: The time derivative of the angular velocity Jacobian defined in Equation (35) is: The time derivative of the linear velocity Jacobian defined in Equation (33) is computed from: where:

Time Derivative of the Base-Spacecraft Inertia Matrix
The time derivative of Equation (28) is given as: By using Equations (11) and (13), the upper right and lower left sub-matrices of Equation (74) can be obtained from the following: The lower right sub-matrix of Equation (74) is computed by using Equation (29):

Time Derivative of the Dynamic-Coupling Inertia Matrix
The time derivative of the dynamic-coupling matrix H 0m in Equation (31) results in:Ḣ where from Equations (32) and (34): Expression of the Matrix c ⋆ From Equation (48), the elements of c ⋆ can be expressed as: For the planar case, Papadopoulos (1990) provides an expression for c ⋆ . For the general case, the symbolic computation of c ⋆ requires the derivatives of H ⋆ with respect to each joint angle q k . By following the analogous procedure as used for the time derivative of H ⋆ (see Equation 65), it yields: As for the time derivative above, the computation of the derivative of the generalized inertia matrix requires the separate computations shown in the following paragraphs.

Joint-Angle Derivative of the Manipulator Inertia Matrix
The joint-angle derivatives of the manipulator inertia matrix from Equation (36) are given by: The calculation of this expression requires the joint-angle derivatives of the joint and link transformation matrices J T J i and J T L i , as well as the joint angle derivatives of the linear and angular motion Jacobians. From Equation (9), the joint transformation matrices are: where the homogeneous transformation matrix J T J k and the DH matrices A(q l , d l , α l , l a l + l b l ) (k + 1 ≤ l ≤ i − 1) are independent of joint angle q k . Therefore, only a single DH transformation matrix in the multiplication chain depends on q k , with the joint angle derivative of a DH matrix given by: −sin q k − cos q k cos α k cos q k sin α k −c k sin q k cos q k − sin q k cos α k sin q k sin α k c k cos q k 0 Therefore, three cases must be considered: (1) If k > i−1, the joint transformation matrix is independent of q k : (2) If k = i − 1, the final factor in the product of transformation matrices depends on q k , and the derivative changes to: (3) For k < i − 1, one of the contributing matrices depends on q k : A q l , d l , α l , l a l + l b l .
With (10), the joint angle derivatives of the link transformation matrices become: As for the link transformation matrices J T Li , there are three cases to consider: (1) For k > i, the link transformation matrix is independent of q k .
(2) For k = i, the final DH transformation matrix must be differentiated for q k : (3) If k < i, we can reuse the derivative of the joint transformation matrix developed above: The joint-angle derivatives of the transformation matrices contain the derivatives of both the direction-cosine matrix and the position. Therefore, the derivatives of the angular and linear motion Jacobians can now be determined from Equations (35) and (33) as follows: where:

Joint-Angle Derivative of the Dynamic-Coupling Matrix
The derivatives of the dynamic-coupling matrix are computed from Equation (31) as: where Since r 0 does not depend directly on any joint angle, r 0i (see Equation 11) is differentiated as:

Joint-Angle Derivative of the Base-Spacecraft Inertia Matrix
The joint-angle derivatives of the base-spacecraft inertia matrix H 0 from Equation (28) result in: The joint-angle derivatives of the relative center-of-mass position vector r 0C from Equation (13) are calculated in: The system moments-of-inertia matrix from Equation (29) is differentiated by: where from Equation (30):

GENERALIZED JACOBIAN
For control purposes, a transformation between joint-space and task-space requires the use of a Jacobian (Siciliano et al., 2010). Typically, the Jacobian is used to map generalized joint velocities to a spatial velocity of the end-effector. However, for reasons of collision avoidance and contact dynamics, Jacobians must also be expressed for the combined velocity of any arbitrary point X i on the i-th link of the manipulator. From Equations (15) and (16), it follows: where the contribution of the manipulator joint rates to the combined velocity of X i is expressed by the [6 × N] manipulator Jacobian J m X i : with x i being the position of point X i . The contribution of the base-spacecraft combined velocity is expressed by the [6 × 6] base-spacecraft Jacobian J 0 X i : where When using the generalized form of the equations of motion, a [6 × N] generalized Jacobian J ⋆ X i can be formulated, such that: Therefore: Using Equation (42) and assuming M 0 = 0 6,1 ,ẋ 0 can be expressed as:ẋ The generalized Jacobian is thus defined as: If x i coincides with the position of the end-effector (p E = p N+1 ), the resulting generalized Jacobian Jacobian defined by Equation (110), for an arbitrary point, becomes the same as the Generalized Jacobian J ⋆ of the floating manipulator as originally defined by Yoshida and Umetani (1993). Therefore, the joint rates required to have the end-effector move at a linear velocity ν E and angular velocity ω E can be calculated from:

IMPLEMENTATION OF A SIMULATION MODEL FOR A FLOATING SPACECRAFT-MANIPULATOR SYSTEM
By using the generalized equations in symbolic form introduced above, a computer model can be implemented for the simulation of the floating spacecraft-manipulator system. This computer model can be based on two implementation approaches, which we call here full symbolic implementation and partial symbolic implementation (see Figure 3). The numerical simulation can thus be run by evaluating at each integration step static functions containing symbolic expressions, without executing iterative procedures such as the recursive Newton-Euler algorithm as in the available literature on the subject (refer to, for instance Mukherjee and Nakamura, 1992;Carignan and Akin, 2000). In the partial symbolic implementation, the numerical values for matrices J T J i , J T L i , H 0 , H 0m , H m , and H -1 0 are calculated from their symbolic expressions. These numerical values are then used to find the generalized matrices H ⋆ , C ⋆ , and J ⋆ . In the full symbolic implementation, H ⋆ , C ⋆ , J ⋆ are calculated directly from their symbolic expressions.
The symbolic implementations are easily adaptable to any number of joints and any structure of the kinematic chain of the spacecraft-manipulator system.

NUMERICAL SIMULATIONS
In order to verify the proposed analytical approach, a simulation model was implemented for two case studies: (1) a planar spacecraft-manipulator system with a 4 DOF manipulator,   as commonly used in ground experimentation of spacecraftmanipulator systems, and (2) a spatial spacecraft-manipulator system with a spatial 6 DOF manipulator, similar to the one described in Yoshida's original discussion of the generalized Jacobian approach (Yoshida and Umetani, 1993). The simulation model was implemented in Matlab/Simulink by following the partial symbolic implementation approach outlined in Figure 3. The block diagram of the overall simulated control architecture is illustrated in Figure 4.
With the complete knowledge of the generalized equations of motion of the spacecraft-manipulator system, the system can be controlled using a standard control scheme called Computed Torque Control (Siciliano et al., 2010). In particular, the commanded torque for the joint motors is computed from the desired joint angular acceleration, angular rate, and angular position by means of a direct dynamics operation. The joint acceleration control input u is computed by means of a proportional-derivative (PD) control law: In general, the controller gains can be designed based on the solutions of the harmonic oscillator: with T max,i being the maximum motor torque for each joint. The control input is finally fed into the dynamic model as acceleration to derive the resulting reaction torque. This torque then serves as the input to the real manipulator system. With a manipulator employed on an orbiting spacecraft, there is always a substantial level of uncertainty when it comes to the inertial properties of the spacecraft-manipulator system, mostly due to the settling state of the fuel tanks and the associated fuel sloshing. Therefore, the computed torque control is based on estimates of the inertia properties matrices,H ⋆ andC ⋆ : In the numerical simulation presented here, perfect knowledge of the inertia properties of the system was assumed, thusH

Simulations With the Planar Four-Link Spacecraft-Manipulator System
The simulated planar spacecraft-manipulator system consists of a base-spacecraft and four identical manipulator links. The mass and inertia properties of the system are given in Table 3, whereas the properties of the kinematic chain are summarized in Table 4. For illustration purposes, the gains of the controller were set to K D i = 1 and K P i = 1 for all joints i. The manipulator is initially fully extended along the x axis, which corresponds to zero angle for each joint. In the sample maneuver sequence one joint at a time accelerates to 0.1 rad/s for 10 s, then rests for 2 s, rotates at −0.1 rad/s for 20 s, rests again for 2 s, then rotates back to the initial condition.
The base-spacecraft is initially located at the origin of the inertial coordinate system, with the three Euler angles being all zero. As shown in Figure 5A, the dynamic coupling between the base-spacecraft and the manipulator generates compensatory motion of the base-spacecraft. As expected, the base-spacecraft has an angular rate component which has opposite sign with respect to the joint rates, and is smaller in magnitude, since the base-spacecraft has higher inertia than the manipulator. Of note is also that the base-spacecraft reaction becomes smaller as the joint maneuver sequence proceeds "outward", since the outer joints move less mass than the inner joints.
Similar behavior is evident in Figure 5B, which illustrates the base-spacecraft position deviation depending on the joint angles. Since the total center-of-mass of the system remains stationary, the base-spacecraft center-of-mass is pulled forward and pushed back to its initial x position, and is pushed to the left and right in the y direction as the manipulator moves to the left and right. The magnitude of the position deviation depends on the combination of distance between base-spacecraft center-of-mass and the center-of-mass of the actuated links, and the total mass or the actuated links. The location of the total center of mass of the multi-body system remains constant.
Since the manipulator joints are in this example purely actuated by internal torques, both the linear momentum and the angular momentum of the spacecraft-manipulator system remain constant, and in this case zero, as expected. Therefore, the linear and angular momenta of the base-spacecraft must be opposite in sign but equal in magnitude to the manipulator momenta. This expected behavior is confirmed in Figures 6A,B. Since the spacecraft-manipulator system is only moving in the x-y plane, there is no linear momentum along the z axis, and no angular momentum about the x and y axes. 5 | Spatial spacecraft-manipulator system with six-link manipulator: mass and inertia properties (Yoshida and Umetani, 1993

Simulations With the Spatial Six-Link Spacecraft-Manipulator System
The simulated spatial spacecraft-manipulator system consists of a base-spacecraft and six manipulator links. The mass and inertia properties of the system are given in Table 5, the DH parameters are summarized in Table 6. For illustration purposes, the gains of the controller were set to K D i = 1 and K P i = 1 for all joints i. The manipulator is initially fully extended along the x axis, which corresponds to zero angle for each joint. In the sample maneuver sequence, one joint at a time accelerates to 0.1 rad/s for 10 s, then rests for 2 s, rotates at −0.1 rad/s for 20 s, rests again for 2 s, then rotates back to zero angle. The base-spacecraft is initially located at the origin of the inertial coordinate system, with the three Euler angles being zero. The rotation axis of joint 1 is along the x axis of the spacecraft coordinate system. Therefore, the rotation of joint 1 causes a 3 DOF rotation of the base-spacecraft, see Figure 7A. Furthermore, it causes the center-of-mass of the base-spacecraft to deviate from its initial position in x, y and z direction during the joint 1 cycle, and only in x and z direction when joints 2-6 are active, as shown in Figure 7B. The location of the total center of mass of the multi-body system remains constant. The linear and angular momenta of base-spacecraft and manipulator compensate each other (see Figure 8) as the total momentum is conserved.

CONCLUSION
This paper presents a complete, step-by-step analytic derivation of the equations of motion of a floating spacecraft-manipulator system, using the Generalized Jacobian approach for modeling the dynamics of a spacecraft-manipulator system. This includes symbolic analytic expressions for all inertia matrices of the spacecraft-manipulator system, as well as their time derivatives and joint-angle derivatives. It also includes a general expression for the Jacobian of a generic point on any of the links of the spacecraft-manipulator system, which is required for the analysis and simulation of collision-avoidance systems and contact dynamics. Furthermore, new and more detailed definitions are proposed in this paper (see Table 1) for the possible modes of maneuvering a spacecraft-manipulator system. In particular, the two commonly used categories (free-flying and free-floating) are  a planar four-link spacecraft-manipulator system, and a spatial six-link spacecraft manipulator system, show that the modeling approach is effective and consistent with the physical principles.
The description of the geometry of the spacecraft-manipulator system uses DH parameters, which allows complete generality. The approach is thus, in principle, extendable to multiple manipulators with open-chain configuration. A fundamental assumption that was used in the description of the kinematics and dynamics of the spacecraft-manipulator system is that of straight links with the center-of-mass of the link being located on the straight line through the origins of the adjacent Cartesian joint coordinate systems. While this is a restriction of generality, it also reflects the type of robotic manipulators commonly used in orbital robotics systems. Furthermore, only rotary joints were considered, since prismatic joints are not commonly used in orbital robotics applications. However, the description of the kinematics and dynamics of the system can be, in principle, straightforwardly extended to include both irregularly shaped links and prismatic joints. Gears in the joints were not considered but they could be added in the analysis by using, for instance, the approach detailed in Siciliano et al. (2010). Euler angles were used in the description of the base-spacecraft orientation: in principle it would be straightforward to use quaternions to avoid any possible orientation singularity problem.
The goal is for this detailed presentation of the Generalized Jacobian approach to serve as a tutorial to build a complete analytic model of the complex dynamics of a spacecraftmanipulator system. While this tutorial is mostly aimed at aerospace engineers faced with the challenge of modeling a robotic system while designing a spacecraft for a rendezvous and capture mission, it is also a good reference for robotic engineers, since it condenses material from distributed sources into one complete presentation.

AUTHOR CONTRIBUTIONS
MW is the principal author. He performed the literature research into the generalized Jacobian approach and derived the mathematical equations to fill the gaps in published literature. He also implemented the Matlab/Simulink model used to generate the examples shown in the paper. He also had the lead in writing the paper. SK worked with MW on checking the mathematical derivations for correctness and on clarifying the presentation of the Denavit-Hartenberg convention. AG worked with MW on the implementation in Matlab/Simulink and on producing the data plots used in the paper. MR provided expertise on the Lagrangian method used to derive the equations of motion and on matrix/vector mathematics.
coordinate system origin between system A and system B T = kinetic energy τ max,i = maximum motor torque for joint i B T A = [4 × 4] homogenous transformation matrix from the coordinate system A to the coordinate system B u = [N × 1] matrix of control inputs for Computed Torque Control v 0 = linear velocity of the base-spacecraft with respect to the inertial frame J v E = linear velocity of the end-effector with respect to the inertial frame J v X i = linear velocity of arbitrary point X i located on link i of the manipulator with respect to the inertial frame J