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 on-orbit 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, multi-body 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 spacecraft-manipulator 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 center-of-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 spacecraft-manipulator 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 center-of-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), spacecraft-manipulator 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 (Nenchev et al., 1999) and contact dynamics models (Nenchev and Yoshida, 1999). 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-by-step 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 spacecraft-manipulator 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 spacecraft-manipulator 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.

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.

TABLE 1

Table 1. Classification of modes of maneuvering for an isolated spacecraft-manipulator system.

Conserved Linear and Angular Momentum

The maneuvering cases here defined as floating and rotation-floating have in previous literature been referred to as free-floating. 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 rotation-floating 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 base-spacecraft 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-of-mass are not actively controlled.

Controlled Linear and Angular Momentum

All of the three maneuvering cases here defined as rotation-flying, 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 rotation-flying 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 spacecraft-manipulator system is never isolated but orbiting an extended body (e.g., the Earth) under its gravitational attraction. However, 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 spacecraft-manipulator system motion.

Kinematics of a Spacecraft-Manipulator System

The geometry of a base-spacecraft equipped with a single N-link 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 + 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).

FIGURE 1

Figure 1. Geometry of a spacecraft-manipulator system.

Three types of Cartesian coordinate systems are used in the analysis: an inertial coordinate system ${J}$, a set of N+2 joint-fixed coordinate systems ${{J}}_i$ (0 ≤ i ≤ N + 1), and a set of N + 1 link-fixed coordinate systems ${{L}}_i$ (0 ≤ i ≤ N), with ${{J}}_0\equiv {{L}}_0$ and ${{J}}_{N+1}={{J}}_{EE}$. The base-spacecraft coordinate system ${{J}}_0\equiv {{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 axis ${\widehat{k}}_i$ for 1 ≤ i ≤ N. The origin ϱ₁ 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 center-of-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.

TABLE 2

Table 2. Denavit-Hartenberg parameters and their geometric meaning.

FIGURE 2

Figure 2. Customized Denavit-Hartenberg kinematic parameters.

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:

$\begin{array}{ll}{ }^{{{J}}_i}{T}_{{{J}}_{i+1}}=A\left({\theta }_i,d_i,{\alpha }_i,c_i\right),& (1)\end{array}$

where the DH transformation matrix function A (β, u, γ, w) is defined as:

$\begin{array}{ll}A\left(\beta ,u,\gamma ,w\right)\triangleq \left[\begin{array}{cccc}\mathit{\text{ cos }}\beta & -\mathit{\text{ sin }}\beta \mathit{\text{ cos }}\gamma & \mathit{\text{ sin }}\beta \mathit{\text{ sin }}\gamma & w\mathit{\text{ cos }}\beta \\ \mathit{\text{ sin }}\beta & \mathit{\text{ cos }}\beta \mathit{\text{ cos }}\gamma & -\mathit{\text{ cos }}\beta \mathit{\text{ sin }}\gamma & w\mathit{\text{ sin }}\beta \\ 0& \mathit{\text{ sin }}\gamma & \mathit{\text{ cos }}\gamma & u\\ 0& 0& 0& 1\end{array}\right].& (2)\end{array}$

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$:

$\begin{array}{ll}{ }^{{J}}{T}_{{J}}{ }_i=\left[\begin{array}{cc}{ }^{{J}}{R}_{{{J}}_i}& { }^{{J}}{p}_i\\ 0_{1,3}& 1\end{array}\right],& (3)\end{array}$

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:

$\begin{array}{ll}{ }^{{J}}{R}_{{{L}}_i}{=}^{{J}}{R}_{{{J}}_{i+1}}.& (4)\end{array}$

The joint DCM are used to derive the directions ${\widehat{k}}_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, ${\widehat{k}}_i$ is defined by:

$\begin{array}{ll}{ }^{{J}}{\widehat{k}}_i{=}^{{J}}{R}_{{{J}}_i}\left(\begin{array}{c}0\\ 0\\ 1\end{array}\right).& (5)\end{array}$

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₀ = p₀, ${{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:

$\begin{array}{ll}{}^{{J}}{R}_{{{L}}_0}{=}^{{J}}{R}_{{{J}}_0}=\left[\begin{array}{ccc}\mathit{\text{ cos }}\psi & -\mathit{\text{ sin }}\psi & 0\\ \mathit{\text{ sin }}\psi & \mathit{\text{ cos }}\psi & 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{ccc}\mathit{\text{ cos }}\theta & 0& \mathit{\text{ sin }}\theta \\ 0& 1& 0\\ -\mathit{\text{ sin }}\theta & 0& \mathit{\text{ cos }}\theta \end{array}\right]\left[\begin{array}{ccc}1& 0& 0\\ 0& \mathit{\text{ cos }}\varphi & -\mathit{\text{ sin }}\varphi \\ 0& \mathit{\text{ sin }}\varphi & \mathit{\text{ cos }}\varphi \end{array}\right].& (6)\end{array}$

Therefore:

$\begin{array}{ll}{}^{{J}}{T}_{{J}}{ }_0{=}^{{J}}{T}_{{L}}{ }_0=\left[\begin{array}{cccc}\mathit{\text{ cos }}\psi \mathit{\text{ cos }}\theta & \mathit{\text{ cos }}\psi \mathit{\text{ sin }}\theta \mathit{\text{ sin }}\psi -\mathit{\text{ sin }}\psi \mathit{\text{ cos }}\psi & \mathit{\text{ cos }}\psi \mathit{\text{ sin }}\theta \mathit{\text{ cos }}\psi +\mathit{\text{ sin }}\psi \mathit{\text{ sin }}\psi & \\ \mathit{\text{ sin }}\psi \mathit{\text{ cos }}\theta & \mathit{\text{ sin }}\psi \mathit{\text{ sin }}\theta \mathit{\text{ sin }}\psi +\mathit{\text{ cos }}\psi \mathit{\text{ cos }}\psi & \mathit{\text{ sin }}\psi \mathit{\text{ sin }}\theta \mathit{\text{ cos }}\psi -\mathit{\text{ cos }}\psi \mathit{\text{ sin }}\psi & { }^I{r}_0\\ -\mathit{\text{ sin }}\theta & \mathit{\text{ cos }}\theta \mathit{\text{ sin }}\psi & \mathit{\text{ cos }}\theta \mathit{\text{ cos }}\psi & \\ 0& 0& 0& 1\end{array}\right].& (7)\end{array}$

The joint coordinate system ${{J}}_1$ has a constant position ${ }^{\mathbf{\text{0}}}{b}_{\mathbf{\text{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:

$\begin{array}{ll}{}^{{J}}{T}_{{J}}{ }_1{=}^{{J}}{T}_{{J}}{0}^{{J}}{ }_0{T}_{{J}}1={=}^{{J}}{T}_{{J}}{ }_0\left[\begin{array}{llll}\hfill \mathit{\text{ cos }}\eta \mathit{\text{ cos }}\zeta \hfill & \hfill \mathit{\text{ cos }}\eta \mathit{\text{ sin }}\zeta \mathit{\text{ sin }}ϵ-\mathit{\text{ sin }}\eta \mathit{\text{ cos }}ϵ\hfill & \hfill \mathit{\text{ cos }}\eta \mathit{\text{ sin }}\zeta \mathit{\text{ cos }}ϵ+\mathit{\text{ sin }}\eta \mathit{\text{ sin }}ϵ\hfill & \hfill \hfill \\ \hfill \mathit{\text{ sin }}\eta \mathit{\text{ cos }}\zeta \hfill & \hfill \mathit{\text{ sin }}\eta \mathit{\text{ sin }}\zeta \mathit{\text{ sin }}ϵ+\mathit{\text{ cos }}\eta \mathit{\text{ cos }}ϵ\hfill & \hfill \mathit{\text{ sin }}\eta \mathit{\text{ sin }}\zeta \mathit{\text{ cos }}ϵ-\mathit{\text{ cos }}\eta \mathit{\text{ sin }}ϵ\hfill & \hfill { }^0{b}_0\hfill \\ \hfill -\mathit{\text{ sin }}\zeta \hfill & \hfill \mathit{\text{ cos }}\zeta \mathit{\text{ sin }}ϵ\hfill & \hfill \mathit{\text{ cos }}\zeta \mathit{\text{ cos }}ϵ\hfill & \hfill \hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \end{array}\right].\hfill & (8)\end{array}$

For all of the subsequent joint coordinate systems ${{J}}_i$, by taking into account Equation (1), it yields:

$\begin{array}{l}{ }^{{J}}{T}_{{J}}{ }_i{=}^{{J}}{T}_{{J}}{{ }_{i-1}}^{{{J}}_{i-1}}{T}_{{{J}}_i}{=}^{{J}}{T}_{{J}}{ }_{i-1}A\left(q_{i-1},d_{i-1},{\alpha }_{i-1},a_{i-1}+b_{i-1}\right),\\ \forall \left\{2\le i\le N\right\},& (9)\end{array}$

where the transformation matrices for the link coordinate systems, ${ }^I{T}_{{{L}}_i}$, are derived as:

$\begin{array}{ll}{ }^{{J}}{T}_{{L}}{ }_i{=}^{{J}}{T}_{{J}}{{ }_{{ }_i}}^{{{ }_{{J}}}_{{ }_i}}{T}_{L_i}{=}^{{J}}{T}_{{J}}{ }_{i}A\left(q_i,d_i,{\alpha }_i,a_i\right),\forall \left\{1\le i\le N\right\}.& (10)\end{array}$

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:

$\begin{array}{ll}{ }^{{J}}{r}_{0i}{=}^{{J}}{r}_i{-}^{{J}}{r}_0,\left\{1\le i\le N\right\},& (11)\end{array}$

$\begin{array}{l}{ }^{{J}}{r}_C=\frac{1}{m_{tot}}{\displaystyle \sum _{i=0}^N{m}_i{}^{{J}}{r}_i}=\frac{1}{m_{tot}}{\displaystyle \sum _{i=0}^N{m}_i\left({ }^{{J}}{r}_0+{}^{{J}}{r}_{0i}\right)}\\ ={}^{{J}}{r}_0+\frac{1}{m_{tot}}{\displaystyle \sum _{i=1}^N{m}_i{}^{{J}}{r}_{0i}},& (12)\end{array}$

where $m_{tot}=\sum _{i=0}^N{m}_i$.

Furthermore, the position of the system center-of-mass with respect to the base-spacecraft center-of-mass is given by:

$\begin{array}{ll}{ }^{{J}}{r}_{0C}={}^{{J}}{r}_C-{}^{{J}}{r}_0.& (13)\end{array}$

Equations (12) and (13) yield:

$\begin{array}{ll}{m}_{tot}{ }^{{J}}{r}_{0C}={\displaystyle \sum _{i=1}^N{m}_i{ }^{{J}}{r}_{0i}}.& (14)\end{array}$

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}}{\omega }_i$ of the i-th link can be expressed as:

$\begin{array}{ll}{ }^{{J}}{v}_{x_i}{=}^{{J}}{v}_0{+}^{{J}}{\omega }_0^{\times }\left({ }^{{J}}{x}_i{-}^{{J}}{r}_0\right)+{\displaystyle \sum _{k=1}^i\left\{\left[{ }^{{J}}{\widehat{k}}_k^{\times }\left({ }^{{J}}{x}_i{-}^{{J}}{p}_k\right)\right]{\stackrel{.}{q}}_k\right\}},& (15)\end{array}$

$\begin{array}{ll}{ }^{{J}}{\omega }_{{L}}{ }_i{=}^{{J}}{\omega }_0+{\displaystyle \sum _{k=1}^i\left({ }^{{J}}{\widehat{k}}_k{\dot{q}}_k\right)}.& (16)\end{array}$

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:

$\begin{array}{ll}{a}^{\times }=\left[\begin{array}{ccc}0& -a_z& a_y\\ a_z& 0& -a_x\\ -a_y& a_x& 0\end{array}\right].& (17)\end{array}$

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 $\overline{q}$ expressing a rotation about any axis $\widehat{u}$ by any angle α can be written as a four-vector:

$\begin{array}{ll}\overline{q}=\left[\begin{array}{c}\epsilon \\ \eta \end{array}\right]=\left[\begin{array}{c}{\epsilon }_1\\ {\epsilon }_2\\ {\epsilon }_3\\ \eta \end{array}\right]=\left[\begin{array}{c}{u}_1\mathit{\text{ sin }}\left(\alpha /2\right)\\ u_2\mathit{\text{ sin }}\left(\alpha /2\right)\\ u_3\mathit{\text{ sin }}\left(\alpha /2\right)\\ \mathit{\text{ cos }}\left(\alpha /2\right)\end{array}\right].& (18)\end{array}$

The magnitude of $\overline{q}$ is $\left|\overline{q}\right|=1$, hence $\overline{q}$ is called a unit quaternion. The unit quaternion thus corresponds to a rotation matrix R (Hughes, 2004):

$\begin{array}{ll}R=\left[\begin{array}{ccc}1-2\left({\epsilon }_2^2+{\epsilon }_3^2\right)& 2({\epsilon }_1{\epsilon }_2-{\epsilon }_3\eta )& 2({\epsilon }_1{\epsilon }_3+{\epsilon }_2\eta )\\ 2({\epsilon }_1{\epsilon }_2+{\epsilon }_3\eta )& 1-2({\epsilon }_1^2+{\epsilon }_3^2)& 2({\epsilon }_2{\epsilon }_3-{\epsilon }_1\eta )\\ 2({\epsilon }_1{\epsilon }_3-{\epsilon }_2\eta )& 2({\epsilon }_2{\epsilon }_3+{\epsilon }_1\eta )& 1-2({\epsilon }_1^2+{\epsilon }_2^2)\end{array}\right].& (19)\end{array}$

A sequence of rotations, R₃ = R₁R₂, is then expressed by the quaternion product (Campa and Camarillo, 2008):

$\begin{array}{ll}{\overline{q}}_3={\overline{q}}_1\otimes {\overline{q}}_2=\left[\begin{array}{c}{\eta }_1{\epsilon }_2+{\eta }_2{\epsilon }_1+{\epsilon }_1^{\times }{\epsilon }_2\\ {\eta }_1{\eta }_2-{\epsilon }_1^T{\epsilon }_2\end{array}\right].& (20)\end{array}$

Therefore, any vector a, expressed as four-vector $\text{a}={\begin{array}{cccc}{[a}_1& a_2& a_3& 0]\end{array}}^T$, can be rotated using a unit quaternion (Kuipers, 2000):

$\begin{array}{ll}a\prime =\overline{q}\otimes a\otimes {\overline{q}}^{*},& (21)\end{array}$

with the conjugate quaternion ${\overline{q}}^{*}={\begin{array}{cc}-[{\epsilon }^T& \eta ]\end{array}}^T$. Therefore, the application of the DH convention results in an unit quaternion ${ }^{{{J}}_i}{{\overline{q}}_r}_{{{J}}_{i+\mathit{\text{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):

$\begin{array}{ll}{ }^{{{J}}_i}{\overline{q}}_r{ }_{{{J}}_{i+1}}=\left[\begin{array}{c}\text{ cos }\frac{{\theta }_i}{2}\mathit{\text{ sin }}\frac{{\alpha }_i}{2}\\ \\ \mathit{\text{ sin }}\frac{{\theta }_i}{2}\text{ cos }\frac{{\alpha }_i}{2}\\ \text{ cos }\frac{{\theta }_i}{2}\text{ cos }\frac{{\alpha }_i}{2}\end{array}\right],& (22)\end{array}$

$\begin{array}{ll}{ }^{{{J}}_i}{t}_{{{J}}_{i+1}}=\left[\begin{array}{c}{c}_i\text{ cos }{\theta }_i\\ c_i\mathit{\text{ sin }}{\theta }_i\\ d_i\\ 0\end{array}\right].& (23)\end{array}$

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):

$\begin{array}{ll}{ }^1{\overline{q}}_r{ }_E{=}^1{\overline{q}}_r{ }_{N+1}{=}^1{\overline{q}}_r{ }_2{\otimes }^2{\overline{q}}_r{ }_3\otimes \dots {\otimes }^{N-1}{\overline{q}}_r{ }_N\otimes {}^N{\overline{q}}_r{ }_{N+1}.& (24)\end{array}$

The sequence of translations results from Equation (21):

$\begin{array}{l}{ }^1{t}_E{=}^1{t}_{N+1}{=}^1{t}_2+\\ {+}^1{\overline{q}}_r{ }_2{\otimes }^2{t}_3{\otimes }^1{\overline{q}}_{r_2}^{*}+\\ +{(}^1{\overline{q}}_r{ }_2{\otimes }^2{\overline{q}}_r{ }_3){\otimes }^3{t}_4\otimes {\left({ }^1{\overline{q}}_r{ }_2{\otimes }^2{\overline{q}}_r{ }_3\right)}^{*}+\dots +\\ +{(}^1{\overline{q}}_r{ }_2\otimes \dots {\otimes }^{N-1}{\overline{q}}_r{ }_N){\otimes }^N{t}_{N+1}\\ \otimes {\left({ }^1{\overline{q}}_r{ }_2\otimes \dots {\otimes }^{N-1}{\overline{q}}_r{ }_N\right)}^{*}.& (25)\end{array}$

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:

$\begin{array}{ll}L=T=\frac{1}{2}{\displaystyle \sum _{i=0}^N\left({ }^{{J}}{\omega }_i{{ }^T}^{{J}}{I}_i{ }^{{J}}{\omega }_i+m_i{ }^{{J}}{\stackrel{.}{r}}_i{{ }^T}^{{J}}{\stackrel{.}{r}}_i\right)}.& (26)\end{array}$

By collecting the mass and inertia properties of the spacecraft-manipulator system into the [6 × 6] base-spacecraft inertia matrix H₀, 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):

$\begin{array}{l}T=\frac{1}{2}\left[{ }^{{J}}{\dot{x}}_0^T{\dot{q}}^T\right]\left[\begin{array}{cc}{H}_0& H_{0m}\\ H_{0m}^T& H_m\end{array}\right]\left[\begin{array}{c}{ }^{{J}}{\stackrel{.}{x}}_0\\ \stackrel{.}{q}\end{array}\right]\\ ={\frac{1}{2}}^{{J}}{\stackrel{.}{x}}_0^T{H}_0{ }^{{J}}{\stackrel{.}{x}}_0+\frac{1}{2}{\dot{q}}^T{H}_m\dot{q}\\ +{\frac{1}{2}}^{{J}}{\dot{x}}_{0}T{H}_{0m}\dot{q}+\frac{1}{2}{\stackrel{.}{q}}^T{H}_{0m}^T{ }^{{J}}{\stackrel{.}{x}}_0,& (27)\end{array}$

where ${ }^{{J}}{\stackrel{.}{x}}_0={\left[\begin{array}{cc}{ }^{{J}}{v}_0{ }^T& { }^{{J}}{\omega }_0^T\end{array}\right]}^T$ is the combined and angular velocity matrix of the base-spacecraft. The base-spacecraft inertia matrix, manipulator inertia matrix, and the dynamic-coupling inertia matrix are thoroughly explained in the following subsections.

Base-Spacecraft Inertia Matrix

The [6 × 6] base-spacecraft inertia matrix H₀ is resulting to be:

$\begin{array}{ll}{H}_0=\left[\begin{array}{cc}{m}_{tot}{𝕀}_{3,3}& -m_{tot}{ }^{{J}}{r}_{0C}^{\times }\\ m_{tot}{ }^{{J}}{r}_{0C}^{\times }& H_S\end{array}\right],& (28)\end{array}$

where 𝕀_{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}$:

$\begin{array}{ll}{H}_S={\displaystyle \sum _{i=1}^N\left({ }^{{J}}{I}_i-m_i{ }^{{J}}{r}_{0i}^{\times }{ }^{{J}}{r}_{0i}^{\times }\right){+}^{{J}}{I}_0}.& (29)\end{array}$

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:

$\begin{array}{ll}{ }^{{J}}{I}_i{=}^{{J}}{R}_{{{L}}_i}{ }^{{{L}}_i}{I}_i{ }^{{J}}{R}_{{{L}}_i}^T.& (30)\end{array}$

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 spacecraft-manipulator system. In particular, it is:

$\begin{array}{ll}{H}_{0m}=\left[\begin{array}{c}{J}_{TS}\\ H_{Sq}\end{array}\right],& (31)\end{array}$

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 rate $\dot{q}$ and the base-spacecraft linear velocity ${ }^{{J}}{v}_0$. In detail, it is:

$\begin{array}{ll}{J}_{TS}={\displaystyle \sum _{i=1}^N\left(m_i{J}_{Ti}\right)},& (32)\end{array}$

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:

$\begin{array}{l}{J}_{Ti}=\left[\begin{array}{cccc}{ }^I{\widehat{k}}_1^{\times }\left({ }^{{J}}{r}_i{-}^{{J}}{p}_1\right)& \cdots & { }^I{\widehat{k}}_i^{\times }\left({ }^{{J}}{r}_i{-}^{{J}}{p}_i\right)& 0_{3,N-i}\end{array}\right],\\ \forall \left(1\le i\le N\right).& (33)\end{array}$

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 rate $\dot{q}$ and the base-spacecraft angular velocity ${ }^{{J}}{\omega }_0$. In detail, it is:

$\begin{array}{ll}{H}_{Sq}={\displaystyle {\displaystyle \sum }_{i=1}^N}\left({ }^{{J}}{I}_i{J}_{Ri}+m_i{ }^{{J}}{r}_{0i}^{\times }{J}_{Ti}\right),& (34)\end{array}$

where the [3 × N] matrix ${J}$_Ri represents the angular velocity Jacobian for the i-th link (Siciliano et al., 2010):

$\begin{array}{ll}{J}_{Ri}=\left[\begin{array}{cccc}{ }^{{J}}{\widehat{k}}_1& \cdots & { }^{{J}}{\widehat{k}}_i& 0_{3,N-i}\end{array}\right],\forall \left(1\le i\le N\right).& (35)\end{array}$

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:

$\begin{array}{ll}{H}_m={\displaystyle \sum _{i=1}^N\left(J_{Ri}^T{ }^{{J}}{I}_i{J}_{Ri}+m_i{J}_{Ti}^T{J}_{Ti}\right)}.& (36)\end{array}$

Equations of Motion of a Floating Spacecraft-Manipulator System

The generalized coordinates for the description of the spacecraft-manipulator 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:

$\begin{array}{ll}\frac{d}{dt}\left(\frac{\partial T}{\partial {\dot{x}}_0}\right)-\frac{\partial T}{\partial x_0}=0_{6,1},& (37)\end{array}$

$\begin{array}{ll}\frac{d}{dt}\left(\frac{\partial T}{\partial \dot{q}}\right)-\frac{\partial T}{\partial q}=\tau ,& (38)\end{array}$

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:

$\begin{array}{ll}\left[\begin{array}{cc}{H}_0& H_{0m}\\ H_{0m}^T& H_m\end{array}\right]\left[\begin{array}{c}{\stackrel{\mathrm{..}}{x}}_0\\ \stackrel{\mathrm{..}}{q}\end{array}\right]+\left[\begin{array}{cc}{\stackrel{.}{H}}_0& {\stackrel{.}{H}}_{0m}\\ {\stackrel{.}{H}}_{0m}^T& {\stackrel{.}{H}}_m\end{array}\right]\left[\begin{array}{c}{\stackrel{.}{x}}_0\\ \stackrel{.}{q}\end{array}\right]+\left[\begin{array}{c}{c}_0\\ c_m\end{array}\right]=\left[\begin{array}{c}{0}_{6,1}\\ \tau \end{array}\right],& (39)\end{array}$

where the [6 × 1] matrix c₀ and the [N × 1] matrix c_m are defined as:

$\begin{array}{ll}{c}_0=-\frac{\partial T}{\partial x_0}=-\frac{1}{2}\frac{\partial }{\partial x_0}\left({\stackrel{.}{x}}_0^T{H}_0{\stackrel{.}{x}}_0+{\stackrel{.}{q}}^T{H}_m\stackrel{.}{q}+{\stackrel{.}{x}}_0^T{H}_{0m}\stackrel{.}{q}+{\stackrel{.}{q}}^T{H}_{0m}^T{\stackrel{.}{x}}_0\right),& (40)\end{array}$

$\begin{array}{ll}{c}_m=-\frac{\partial T}{\partial q}=-\frac{1}{2}\frac{\partial }{\partial q}\left({\stackrel{.}{x}}_0^T{H}_0{\stackrel{.}{x}}_0+{\stackrel{.}{q}}^T{H}_m\stackrel{.}{q}+{\stackrel{.}{x}}_0^T{H}_{0m}\stackrel{.}{q}+{\stackrel{.}{q}}^T{H}_{0m}^T{\stackrel{.}{x}}_0\right).& (41)\end{array}$

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):

$\begin{array}{ll}\left[\begin{array}{c}P\\ L\end{array}\right]=H_0{\stackrel{.}{x}}_0+H_{0m}\stackrel{.}{q}=M_0,& (42)\end{array}$

where $M_0={\left[\begin{array}{cc}P& L\end{array}\right]}^T$ is the [61] combined matrix of initial linear and angular momentum about the base-spacecraft center-of-mass, 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:

$\begin{array}{ll}\frac{d}{dt}\left(\begin{array}{c}P\\ L\end{array}\right)=H_0{\stackrel{\mathrm{..}}{x}}_0+H_{0m}\stackrel{\mathrm{..}}{q}+{\stackrel{.}{H}}_0{\stackrel{.}{x}}_0+{\stackrel{.}{H}}_{0m}\stackrel{.}{q}=0.& (43)\end{array}$

Equation (42) is solved for ${\dot{x}}_0$. The result is inserted into Equation (43), which is then solved for ${\ddot{x}}_0$. The resulting expressions are finally inserted into the last N scalar equations of Equation (39). By considering the symmetry of $H_0^{-\mathbf{\text{1}}}$ (which will be discussed below), the resulting matricial equation (corresponding to N scalar equations) is:

$\begin{array}{l}{H}^{\star }\stackrel{\mathrm{..}}{q}+H^{·\star }\dot{q}-\frac{1}{2}\frac{\partial }{\partial q}\left({\stackrel{.}{q}}^T{H}^{\star \stackrel{.}{q}}\right)=\tau -\frac{d}{dt}\left(H_{0m}^T{H}_0^{-1}\right)M_0\\ +\frac{1}{2}\frac{\partial }{\partial q}\left(M_0^T{H}_0^{-T}{M}_0\right),& (44)\end{array}$

where the [N × N] generalized inertia matrix H^⋆ is defined to be:

$\begin{array}{ll}{H}^{\star }=H_m-H_{0m}^T{H}_0^{-1}{H}_{0m}.& (45)\end{array}$

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

$\begin{array}{ll}{M}_0=0_{6,1}.& (46)\end{array}$

Therefore, Equation (44) simplifies to:

$\begin{array}{ll}{H}^{\star }\stackrel{\mathrm{..}}{q}+H^{·\star }\stackrel{.}{q}+c^{\star }=\tau ,& (47)\end{array}$

where the [N × 1] matrix c^⋆ is defined as:

$\begin{array}{ll}{c}^{\star }=-\frac{1}{2}\frac{\partial }{\partial q}\left({\stackrel{.}{q}}^T{H}^{\star }\stackrel{.}{q}\right)=-\frac{1}{2}{\stackrel{.}{q}}^T\frac{\partial H^{\star }}{\partial q}\stackrel{.}{q}.& (48)\end{array}$

The velocity and position dependent terms in the generalized equations of motion can be combined into a single [N × 1] matrix C^⋆:

$\begin{array}{ll}{C}^{\star }\left(q,\dot{q}\right)=H{\stackrel{.}{ }}^{\star }\dot{q}+c^{\star }.& (49)\end{array}$

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):

$\begin{array}{ll}{H}^{\star }\stackrel{\mathrm{..}}{q}+C^{\star }(q,\dot{q})=\tau .& (50)\end{array}$

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) requires the inverse of the [6 × 6] base-spacecraft inertia matrix H₀. The matrix H₀ is a partitioned block matrix of the form:

$\begin{array}{ll}{H}_0=\left[\begin{array}{cc}{m}_{tot}{𝕀}_{3,3}& -m_{tot}{r}_{0C}^{\times }\\ m_{tot}{r}_{0C}^{\times }& H_S\end{array}\right]=\left[\begin{array}{cc}U& V\\ W& X\end{array}\right].& (51)\end{array}$

By following the approach taken by Mukherjee and Nakamura (1992), the inverse of H₀ is determined using the Banachiewicz inversion formula (Baksalary and Styan, 2002) and the Schur complement S_U of the non-singular submatrix U:

$\begin{array}{l}{H}_0^{-1}=\left[\begin{array}{cc}{𝕀}_{3,3}& -U^{-1}V\\ 0_{3,3}& {𝕀}_{3,3}\end{array}\right]\left[\begin{array}{cc}{U}^{-1}& 0_{3,3}\\ 0_{3,3}& S_U^{-1}\end{array}\right]\left[\begin{array}{cc}{𝕀}_{3,3}& 0_{3,3}\\ -W{U}^{-1}& {𝕀}_{3,3}\end{array}\right]\\ =\left[\begin{array}{cc}\frac{1}{m_{tot}}{𝕀}_{3,3}-r_{0C}^{\times }{S}_U^{-1}{r}_{0C}^{\times }& r_{0C}^{\times }{S}_U^{-1}\\ S_U^{-1}{r}_{0C}^{\times }& S_U^{-1}\end{array}\right],& (52)\end{array}$

where the Schur complement of U for a matrix $\left[\begin{array}{cc}U& V\\ W& X\end{array}\right]$ is defined as Baksalary and Styan (2002):

$\begin{array}{ll}{S}_U=X-W{U}^{-1}V=H_S+m_{tot}{r}_{0C}^{\times }{r}_{0C}^{\times }.& (53)\end{array}$

In detail, S_U is a symmetric [3 × 3] matrix with the elements:

$\begin{array}{ll}{s}_{1,1}={\displaystyle \sum _{i=1}^N\left(I_{x{x}_i}+m_i\left(r_{0i_y}^2+r_{0i_z}^2\right)\right)}+I_{x{x}_0}-m_{tot}\left(r_{0C_y}^2+r_{0C_z}^2\right),& (54)\end{array}$

$\begin{array}{ll}{s}_{1,2}=s_{2,1}={\displaystyle \sum _{i=1}^N\left(I_{x{y}_i}-m_i{r}_{0i_x}{r}_{0i_y}\right)}+I_{x{y}_0}+m_{tot}{r}_{0C_x}{r}_{0C_y},& (55)\end{array}$

$\begin{array}{ll}{s}_{1,3}=s_{3,1}={\displaystyle \sum _{i=1}^N\left(I_{x{z}_i}-m_i{r}_{0i_x}{r}_{0i_z}\right)}+I_{x{z}_0}+m_{tot}{r}_{0C_x}{r}_{0C_z},& (56)\end{array}$

$\begin{array}{ll}{s}_{2,2}={\displaystyle \sum _{i=1}^N\left(I_{y{y}_i}+m_i\left(r_{0i_y}^2+r_{0i_z}^2\right)\right)+I_{y{y}_0}-m_{tot}\left(r_{0C_x}^2+r_{0C_z}^2\right)},& (57)\end{array}$

$\begin{array}{ll}{s}_{2,3}=s_{3,2}={\displaystyle \sum _{i=1}^N\left(I_{y{z}_i}-m_i{r}_{0i_y}{r}_{0i_z}\right)}+I_{y{z}_0}+m_{tot}{r}_{0C_y}{r}_{0C_z},& (58)\end{array}$

$\begin{array}{ll}{s}_{3,3}={\displaystyle \sum _{i=1}^N\left(I_{z{z}_i}+m_i\left(r_{0i_y}^2+r_{0i_z}^2\right)\right)}+I_{z{z}_0}-m_{tot}\left(r_{0C_x}^2+r_{0C_y}^2\right).& (59)\end{array}$

The determinant of S_U is given by:

$\begin{array}{l}\det \left(S_U\right)=-s_{13}^2{s}_{22}+2s_{12}{s}_{13}{s}_{23}-s_{11}{s}_{23}^2-s_{12}^2{s}_{33}\\ \text{+ }{s}_{\text{11}}{s}_{\text{22}}{s}_{\text{33}}\text{.}& (60)\end{array}$

Therefore, S_U can be inverted symbolically, resulting in the symmetric [3 × 3] matrix:

$\begin{array}{l}{S}_U^{-1}=\\ \frac{1}{det\left(S_U\right)}\left[\begin{array}{ccc}{s}_{22}{s}_{33}-s_{23}^2& s_{13}{s}_{23}-s_{12}{s}_{33}& s_{12}{s}_{23}-s_{13}{s}_{22}\\ s_{13}{s}_{23}-s_{12}{s}_{33}& s_{11}{s}_{33}-s_{13}^2& s_{12}{s}_{13}-s_{11}{s}_{23}\\ s_{12}{s}_{23}-s_{13}{s}_{22}& s_{12}{s}_{13}-s_{11}{s}_{23}& s_{11}{s}_{22}-s_{12}^2\end{array}\right]& (61)\end{array}$

This symbolic expression for S_U allows the calculation of the symbolic expression of $H_0^{-\mathbf{\text{1}}}$ by using Equation (52).

Time Derivative of the Generalized Inertia Matrix

The time derivative of the generalized inertia matrix, ${\dot{H}}^{\star }$, is calculated by taking the time derivative of Equation (45), resulting in:

$\begin{array}{ll}{\stackrel{.}{H}}^{\star }=\frac{d}{dt}\left(H_m-H_{0m}^T{H}_0^{-1}{H}_{0m}\right).& (62)\end{array}$

The derivative of a matrix inverse with respect to a scalar can be expressed as:

$\begin{array}{ll}\frac{d{Q}^{-1}}{dx}=-Q^{-1}\frac{dQ}{dx}{Q}^{-1},& (63)\end{array}$

as it can be immediately obtained from

$\begin{array}{ll}\frac{d\left(Q{Q}^{-1}\right)}{dx}=\frac{dE}{dx}=0=\frac{dQ}{dx}{Q}^{-1}+Q\frac{d{Q}^{-1}}{dx},& (64)\end{array}$

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_0^{-\mathbf{\text{1}}}$:

$\begin{array}{l}{\stackrel{.}{H}}^{\star }={\stackrel{.}{H}}_m-({\stackrel{.}{H}}_{0m}^T{H}_0^{-1}{H}_{0m}+H_{0m}^T{H}_0^{-1}{\stackrel{.}{H}}_{0m}\\ -H_{0m}^T{H}_0^{-1}{\stackrel{.}{H}}_0{H}_0^{-1}{H}_{0m}).& (65)\end{array}$

The following subsections discuss the calculation of the individual matrix derivatives in Equation (65).

Time Derivative of the Manipulator Inertia Matrix

By taking the time derivative of Equation (36), it yields:

$\begin{array}{l}{\stackrel{.}{H}}_m={\displaystyle \sum _{i=1}^N({\stackrel{.}{J}}_{Ri}^T{ }^{{J}}{I}_i{J}_{Ri}+J_{Ri}^T{}^{{J}}{\stackrel{.}{I}}_i{J}_{Ri}+J_{Ri}^T{ }^{{J}}{I}_i{\stackrel{.}{J}}_{Ri}}\\ +m_i\left({\stackrel{.}{J}}_{Ti}^T{J}_{Ti}+J_{Ti}^T{\stackrel{.}{J}}_{Ti}\right))& (66)\end{array}$

The only new values to be calculated are the time derivatives of the angular and linear motion Jacobians, ${\dot{J}}_{Ri}$ and ${\dot{J}}_{Ti}$, as well as the time derivative of the link moments-of-inertia matrices, ${ }^{{J}}{\dot{I}}_i$ (1 ≤ i ≤ N). For the calculation of these terms, the time derivatives of the joint and link direction-cosine matrices, ${ }^{{J}}{\dot{R}}_{{{J}}_i}$ and ${ }^{{J}}{\dot{R}}_{{{L}}_i}$ are required. The time derivative of any joint direction-cosine matrix ${ }^{{J}}{R}_{{{J}}_{\mathbf{\text{i}}}}$ is given by the Darboux equation (Hughes, 2004):

$\begin{array}{ll}{ }^{{J}}{\stackrel{.}{R}}_{{{J}}_i}{=}^{{J}}{\omega }_{{{J}}_i}^{\times }{ }^{{J}}{R}_{{{J}}_i}.& (67)\end{array}$

The angular velocity of any joint coordinate system with respect to the inertial frame can be computed recursively from:

$\begin{array}{l}{ }^{{J}}{\omega }_{{{J}}_i}{=}^{{J}}{\omega }_{{{J}}_{i-1}}{+}^{i-1}{\omega }_{{{J}}_i}{=}^{{J}}{\omega }_{{{J}}_{i-1}}{+}^{{J}}\widehat{k}{}_{{{J}}_{i-1}}{\stackrel{.}{q}}_{i-1},\\ \forall \left(1\le i\le N+1\right).& (68)\end{array}$

Due to the definition of the joint and link coordinate systems in section Kinematics of a Spacecraft-Manipulator System:

$\begin{array}{ll}{ }^{{J}}{\stackrel{.}{R}}_{{{L}}_i}{=}^{{J}}{\stackrel{.}{R}}_{{{J}}_{i+1}},\forall \left(1\le i\le N\right).& (69)\end{array}$

Furthermore, the time derivative of Equation (30) yields:

$\begin{array}{l}{}^{{J}}{\stackrel{.}{I}}_i{=}^{{J}}{\stackrel{.}{R}}_{{{L}}_i}{}^{{{L}}_i}{I}_i{}^{{J}}{R}_{{{L}}_i}^T{+}^{{J}}{R}_{{{L}}_i}{}^{{{L}}_i}{I}_i{}^{{J}}{\stackrel{.}{R}}_{{{L}}_i}^T,\\ \forall \left(1\le i\le N\right).& (70)\end{array}$

The time derivative of the angular velocity Jacobian defined in Equation (35) is:

$\begin{array}{l}{\stackrel{.}{J}}_{Ri}=\left[\begin{array}{ccccc}{ }^{{J}}{\stackrel{.}{R}}_{{{J}}_1}\left(\begin{array}{c}0\\ 0\\ 1\end{array}\right)& { }^{{J}}{\stackrel{.}{R}}_{{{J}}_2}\left(\begin{array}{c}0\\ 0\\ 1\end{array}\right)& \dots & { }^{{J}}{\stackrel{.}{R}}_{{{J}}_i}\left(\begin{array}{c}0\\ 0\\ 1\end{array}\right)& 0_{3,N-1}\end{array}\right],\\ \forall \left(1\le i\le N\right).& (71)\end{array}$

The time derivative of the linear velocity Jacobian defined in Equation (33) is computed from:

$\begin{array}{ll}{\stackrel{.}{J}}_{Ti}=[\begin{array}{ccccc}{\stackrel{.}{j}}_{T{i}_1}& {\stackrel{.}{j}}_{T{i}_2}& \cdots & {\stackrel{.}{j}}_{T{i}_i}& 0_{3,N-1}\end{array}],\forall \left(1\le i\le N\right),& (72)\end{array}$

where:

$\begin{array}{ll}{\stackrel{.}{j}}_{T{i}_k}=\begin{array}{c}{\left[{ }^{{J}}{\stackrel{.}{R}}_{{{J}}_k}\left(\begin{array}{c}0\\ 0\\ 1\end{array}\right)\right]}^{\times }\left(r_i-p_k\right)+{\widehat{k}}_k^{\times }\left({\stackrel{.}{r}}_i-{\stackrel{.}{p}}_k\right),\forall \left(k\le i\right)\end{array}.& (73)\end{array}$

Time Derivative of the Base-Spacecraft Inertia Matrix

The time derivative of Equation (28) is given as:

$\begin{array}{ll}{\stackrel{.}{H}}_0=\left[\begin{array}{cc}{0}_{3,3}& -m_{tot}{\stackrel{.}{r}}_{0C}^{\times }\\ m_{tot}{\stackrel{.}{r}}_{0C}^{\times }& {\stackrel{.}{H}}_S\end{array}\right].& (74)\end{array}$

By using Equations (11) and (13), the upper right and lower left sub-matrices of Equation (74) can be obtained from the following:

$\begin{array}{l}{m}_{tot}{\stackrel{.}{r}}_{0C}^{\times }={\displaystyle \sum _{i=1}^N{m}_i{\stackrel{.}{r}}_{0i}^{\times }}={\displaystyle \sum _{i=1}^N{m}_i{\left({\stackrel{.}{r}}_i-{\stackrel{.}{r}}_0\right)}^{\times }}={\displaystyle \sum _{i=1}^N{m}_i}\\ {\left({\omega }_0{ }^{\times }\left(r_i-r_0\right)+{\displaystyle \sum _{j=1}^i\left\{\left[{\left[{ }^{{J}}{R}_{{{J}}_j}\left(\begin{array}{c}0\\ 0\\ 1\end{array}\right)\right]}^{\times }\left(r_i-p_j\right)\right]{\dot{q}}_j\right\}}\right)}^{\times }.& (75)\end{array}$

The lower right sub-matrix of Equation (74) is computed by using Equation (29):

$\begin{array}{ll}{\stackrel{.}{H}}_S={\displaystyle \sum _{i=1}^N\left({ }^{{J}}{\stackrel{.}{I}}_i-m_i\left({\stackrel{.}{r}}_{0i}^{\times }{r}_{0i}^{\times }+r_{0i}^{\times }{\stackrel{.}{r}}_{0i}^{\times }\right)\right)}{+}^{{J}}{\stackrel{.}{I}}_0.& (76)\end{array}$

Time Derivative of the Dynamic-Coupling Inertia Matrix

The time derivative of the dynamic-coupling matrix H_0m in Equation (31) results in:

$\begin{array}{ll}{\stackrel{.}{H}}_{0m}=\left[\begin{array}{c}{\stackrel{.}{J}}_{TS}\\ {\stackrel{.}{H}}_{Sq}\end{array}\right],& (77)\end{array}$

where from Equations (32) and (34):

$\begin{array}{ll}{\stackrel{.}{J}}_{TS}={\displaystyle \sum _{i=1}^N{m}_i{\stackrel{.}{J}}_{Ti}},& (78)\end{array}$

$\begin{array}{ll}{\stackrel{.}{H}}_{Sq}={\displaystyle \sum _{i=1}^N\left({ }^{{J}}{\stackrel{.}{I}}_i{J}_{Ri}{+}^{{J}}{I}_i{\stackrel{.}{J}}_{Ri}+m_i\left({\stackrel{.}{r}}_{0i}^{\times }{J}_{Ti}+r_{0i}^{\times }{\stackrel{.}{J}}_{Ti}\right)\right)}.& (79)\end{array}$

Expression of the Matrix c^⋆

From Equation (48), the elements of c^⋆ can be expressed as:

$\begin{array}{ll}{c}_k^{\star }=-\frac{1}{2}{\stackrel{.}{q}}^T\frac{\partial H^{\star }}{\partial q_k}\stackrel{.}{q},\forall \left(1\le k\le N\right),& (80)\end{array}$

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:

$\begin{array}{l}\frac{\partial H^{\star }}{\partial q_k}=\frac{\partial H_m}{\partial q_k}-{\frac{\partial H_{0m}}{\partial q_k}}^T{H}_0^{-1}{H}_{0m}-H_{0m}^T{H}_0^{-1}\frac{\partial H_{0m}}{\partial q_k}\\ +H_{0m}^T{H}_0^{-1}\frac{\partial H_0}{\partial q_k}{H}_0^{-1}{H}_{0m},\forall \left(1\le k\le N\right)& (81)\end{array}$

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:

$\begin{array}{l}\frac{\partial H_m}{\partial q_k}={\displaystyle \sum { }_{i=1}^N}({\frac{\partial J_{Ri}}{\partial q_k}}^T{ }^{{J}}{I}_i{J}_{Ri}+J_{Ri}^T\frac{{\partial }^{{J}}{I}_i}{\partial q_k}{J}_{Ri}+J_{Ri}^T{ }^{{J}}{I}_i\frac{\partial J_{Ri}}{\partial q_k}\\ +m_i{\frac{\partial J_{Ti}}{\partial q_k}}^T{J}_{Ti}+m_i{J}_{Ti}{ }^T\frac{\partial J_{Ti}}{\partial q_k}),\forall \left(1\le k\le N\right)& (82)\end{array}$

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:

$\begin{array}{l}{ }^{{J}}{T}_{{J}}{ }_i{=}^{{J}}{T}_{{J}}{ }_k•A\left(q_k,d_k,{\alpha }_k{,}^k{a}_k{+}^k{b}_k\right)\\ •{\displaystyle \prod _{l=k+1}^{i-1}A\left(q_l,d_l,{\alpha }_l{,}^l{a}_l{+}^l{b}_l\right)}.& (83)\end{array}$

where the homogeneous transformation matrix ${ }^{{J}}{T}_{{J}}{ }_k$ and the DH matrices A($q_l,d_l,{\alpha }_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:

$\begin{array}{l}\frac{\partial A\left(q_k,d_k,{\alpha }_k{,}^k{c}_k\right)}{\partial q_k}=\\ \left[\begin{array}{cccc}-\mathit{\text{ sin }}{q}_k& - \mathit{\text{ cos }}{q}_k\mathit{\text{ cos }}{\alpha }_k& \mathit{\text{ cos }}{q}_k\mathit{\text{ sin }}{\alpha }_k& -c_k\mathit{\text{ sin }}{q}_k\\ \mathit{\text{ cos }}{q}_k& - \mathit{\text{ sin }}{q}_k\mathit{\text{ cos }}{\alpha }_k& \mathit{\text{ sin }}{q}_k\mathit{\text{ sin }}{\alpha }_k& c_k\mathit{\text{ cos }}{q}_k\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right].& (84)\end{array}$

Therefore, three cases must be considered:

(1) If k > i − 1, the joint transformation matrix is independent of q_k:

$\begin{array}{ll}{\frac{{\partial }^{{J}}{T}_{{J}}{ }_i}{\partial q_k}|}_{k>i-1}=0_{4,4},& (85)\end{array}$

(2) If k = i − 1, the final factor in the product of transformation matrices depends on q_k, and the derivative changes to:

$\begin{array}{ll}{\frac{{\partial }^{{J}}{T}_{{J}}{ }_i}{\partial q_k}|}_{k=i-1}{=}^{{J}}{T}_{{J}}{ }_{i-1}•\frac{\partial A\left(q_{i-1},d_{i-1},{\alpha }_{i-1}{,}^{i-1}{a}_{i-1}{+}^{i-1}{b}_{i-1}\right)}{\partial q_{i-1}},& (86)\end{array}$

(3) For k < i − 1, one of the contributing matrices depends on q_k:

$\begin{array}{l}{\frac{{\partial }^{{J}}{T}_{{J}}{ }_i}{\partial q_k}|}_{k<i-1}{=}^{{J}}{T}_{{J}}{ }_k•\frac{\partial A\left(q_k,d_k,{\alpha }_k{,}^k{a}_k{+}^k{b}_k\right)}{\partial q_k}\\ •{\displaystyle \prod _{l=k+1}^{i-1}A\left(q_l,d_l,{\alpha }_l{,}^l{a}_l{+}^l{b}_l\right)}.& (87)\end{array}$

With (10), the joint angle derivatives of the link transformation matrices become:

$\begin{array}{ll}\frac{{\partial }^{{J}}{T}_{{L}}{ }_i}{\partial q_k}=\frac{{\partial }^{{J}}{T}_{{J}}{ }_i}{\partial q_k}A\left(q_i,d_i,{\alpha }_i{,}^i{a}_i\right)+{}^{{J}}{T}_{{J}}{ }_i\frac{\partial A\left(q_i,d_i,{\alpha }_i{,}^i{a}_i\right)}{\partial q_k}.& (88)\end{array}$

As for the link transformation matrices ${ }^{{J}}{T}_{{L}}{ }_i$, there are three cases to consider:

(1) For k > i, the link transformation matrix is independent of q_k.

$\begin{array}{ll}{\frac{{\partial }^{{J}}{T}_{{L}}{ }_i}{\partial q_k}|}_{k>i}=0_{4,4},& (89)\end{array}$

(2) For k = i, the final DH transformation matrix must be differentiated for q_k:

$\begin{array}{ll}{\frac{{\partial }^{{J}}{T}_{{L}}{ }_i}{\partial q_k}|}_{k=i}={}^{{J}}{T}_{{J}}{ }_i•\frac{\partial A\left(q_i,d_i,{\alpha }_i{,}^i{a}_i\right)}{\partial q_i},& (90)\end{array}$

(3) If k < i, we can reuse the derivative of the joint transformation matrix developed above:

$\begin{array}{ll}\frac{{\partial }^{{J}}{T}_{{L}}{ }_i}{\partial q_k}=\frac{{\partial }^{{J}}{T}_{{J}}{ }_i}{\partial q_k}A\left(q_i,d_i,{\alpha }_i{,}^i{a}_i\right).& (91)\end{array}$

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:

$\begin{array}{ll}\frac{\partial J_{Ri}}{\partial q_k}=\left[\frac{{\partial }^{{J}}{R}_{{{J}}_1}}{\partial q_k}\begin{array}{cccc}\left(\begin{array}{c}0\\ 0\\ 1\end{array}\right)& \cdots & \frac{{\partial }^{{J}}{R}_{{{J}}_i}}{\partial q_k}\left(\begin{array}{c}0\\ 0\\ 1\end{array}\right)& 0_{3,N-i}\end{array}\right],\begin{array}{c}\forall \left(1\le i\le N\right)\\ \forall \left(1\le k\le N\right),\end{array}& (92)\end{array}$

$\begin{array}{ll}\frac{\partial J_{Ti}}{\partial q_k}=\left[\begin{array}{ccccc}\frac{\partial J_{T_i{ }_1}}{\partial q_k}& \frac{\partial J_{T_i{ }_2}}{\partial q_k}& \cdots & \frac{\partial J_{T_i{ }_i}}{\partial q_k}& 0_{3,N-i}\end{array}\right],\begin{array}{c}\forall \left(1\le i\le N\right)\\ \forall \left(1\le k\le N\right).\end{array}& (93)\end{array}$

where:

$\begin{array}{ll}\frac{\partial J_{T_{i_1}}}{\partial q_k}={\left[\frac{{\partial }^{{J}}{R}_{{{J}}_1}}{\partial q_k}\left(\begin{array}{c}0\\ 0\\ 1\end{array}\right)\right]}^{\times }\left(r_i-p_l\right)+{\widehat{k}}_l^{\times }\left(\frac{\partial r_i}{\partial q_k}-\frac{\partial p_l}{\partial q_k}\right),\begin{array}{c}\forall \left(l\le i\right)\\ \forall \left(1\le k\le N\right).\end{array}& (94)\end{array}$

Joint-Angle Derivative of the Dynamic-Coupling Matrix

The derivatives of the dynamic-coupling matrix are computed from Equation (31) as:

$\begin{array}{ll}\frac{\partial H_{0m}}{\partial q_k}=\left[\begin{array}{c}\frac{\partial J_{TS}}{\partial q_k}\\ \frac{\partial H_{Sq}}{\partial q_k}\end{array}\right], \forall (1\le k\le N),& (95)\end{array}$

where

$\begin{array}{ll}\frac{\partial J_{TS}}{\partial q_k}={\displaystyle \sum _{i=1}^N\left(m_i\frac{\partial J_{Ti}}{\partial q_k}\right)},\forall \left(1\le k\le N\right),& (96)\end{array}$

$\begin{array}{l}\frac{\partial H_{Sq}}{\partial q_k}={\displaystyle \sum _{i=1}^N\left\{\frac{{\partial }^{{J}}{I}_i}{\partial q_k}{J}_{Ri}{+}^{{J}}{I}_i\frac{\partial J_{Ri}}{\partial q_k}+m_i\left[{\left(\frac{{\partial }^{{J}}{r}_{0i}}{\partial q_k}\right)}^{\times }{J}_{Ti}+r_{0i}{ }^{\times }\frac{\partial J_{Ti}}{\partial q_k}\right]\right\}},\\ \forall \left(1\le k\le N\right).& (97)\end{array}$

Since r₀ does not depend directly on any joint angle, r_0i (see Equation 11) is differentiated as:

$\begin{array}{ll}\frac{\partial r_{0i}}{\partial q_k}=\frac{\partial r_i}{\partial q_k},\begin{array}{c}\forall \left(1\le i\le N\right)\\ \forall \left(1\le k\le N\right).\end{array}& (98)\end{array}$

Joint-Angle Derivative of the Base-Spacecraft Inertia Matrix

The joint-angle derivatives of the base-spacecraft inertia matrix H₀ from Equation (28) result in:

$\begin{array}{ll}\frac{\partial H_0}{\partial q_k}=\left[\begin{array}{cc}{0}_{3,3}& -m_{tot}{\left(\frac{\partial r_{0C}}{\partial q_k}\right)}^{\times }\\ m_{tot}{\left(\frac{\partial r_{0C}}{\partial q_k}\right)}^{\times }& \frac{\partial H_S}{\partial q_k}\end{array}\right].& (99)\end{array}$

The joint-angle derivatives of the relative center-of-mass position vector r_0C from Equation (13) are calculated in:

$\begin{array}{l}\frac{\partial r_{0C}}{\partial q_k}=\frac{1}{m_{tot}}{\displaystyle \sum _{i=1}^N\left(m_i\frac{\partial r_{0i}}{\partial q_k}\right)}=\frac{1}{m_{tot}}{\displaystyle \sum _{i=1}^N\left(m_i\frac{\partial r_i}{\partial q_k}\right)},\\ \forall \left(1\le k\le N\right).& (100)\end{array}$

The system moments-of-inertia matrix from Equation (29) is differentiated by:

$\begin{array}{l}\frac{\partial H_S}{\partial q_k}={\displaystyle \sum _{i=1}^N\left\{\frac{{\partial }^{{J}}{I}_i}{\partial q_k}-m_i\left[{\left(\frac{\partial r_i}{\partial q_k}\right)}^{\times }{r}_{0i}{ }^{\times }+r_{0i}{ }^{\times }{\left(\frac{\partial r_i}{\partial q_k}\right)}^{\times }\right]\right\}},\\ \forall \left(1\le k\le N\right),& (101)\end{array}$

where from Equation (30):

$\begin{array}{l}\frac{{\partial }^{{J}}{I}_i}{\partial q_k}=\frac{{\partial }^{{J}}{R}_{{{L}}_i}}{\partial q_k}{ }^{{{L}}_i}{I}_i{ }^{{J}}{R}_{{{L}}_i}^T{+}^{{J}}{R}_{{{L}}_i}{ }^{{{L}}_i}{I}_i{\frac{{\partial }^{{J}}{R}_{{{L}}_i}}{\partial q_k}}^T\\ \forall \left(1\le k\le N\right).& (102)\end{array}$

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:

$\begin{array}{ll}\left[\begin{array}{c}{\nu }_{X_i}\\ {\omega }_i\end{array}\right]=J_{0_{X_i}}{\stackrel{.}{x}}_0+J_{m_{X_i}}\stackrel{.}{q},& (103)\end{array}$

where the contribution of the manipulator joint rates to the combined velocity of X_i is expressed by the [6 × N] manipulator Jacobian J_mXi:

$\begin{array}{ll}{J}_{m_{X_i}}=\left[\begin{array}{cccc}{\widehat{k}}_1^{\times }\left(x_i-p_1\right)& \dots & {\widehat{k}}_i^{\times }\left(x_i-p_i\right)& 0_{3,N-i}\\ {\widehat{k}}_1& \dots & {\widehat{k}}_i& 0_{3,N-i}\end{array}\right],& (104)\end{array}$

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_0Xi:

$\begin{array}{ll}{J}_{0_{X_i}}=\left[\begin{array}{cc}E& -x_{0i}^{\times }\\ 0_{3,3}& E\end{array}\right],& (105)\end{array}$

where

$\begin{array}{ll}{x}_{0i}=x_i-p_0.& (106)\end{array}$

When using the generalized form of the equations of motion, a [6 × N] generalized Jacobian $J_{X_i}^{\star }$ can be formulated, such that:

$\begin{array}{ll}\left[\begin{array}{c}{ }^{{J}}{\nu }_{X_i}\\ { }^{{J}}{\omega }_i\end{array}\right]=J_{X_i}^{\star }\stackrel{.}{q}.& (107)\end{array}$

Therefore:

$\begin{array}{ll}{J}_{X_i}^{\star }\stackrel{.}{q}=J_{0_{X_i}}{\stackrel{.}{x}}_0+J_{m_{X_i}}\stackrel{.}{q}.& (108)\end{array}$

Using Equation (42) and assuming M₀ = 0_{6, 1}, ${\dot{x}}_0$ can be expressed as:

$\begin{array}{ll}{\stackrel{.}{x}}_0=-H_0^{-1}{H}_{0m}\stackrel{.}{q}.& (109)\end{array}$

The generalized Jacobian is thus defined as:

$\begin{array}{ll}{J}_{X_i}*=J_{m_{X_i}}-J_{0_X{ }_i}{H}_0^{-1}{H}_{0m}& (110)\end{array}$

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:

$\begin{array}{ll}\stackrel{.}{q}=J^{{\star }^{-1}}\left[\begin{array}{c}{\nu }_E\\ {\omega }_E\end{array}\right]& (111)\end{array}$

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₀, H_0m, H_m, and $H_0^{\mathbf{\text{-1}}}$ 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.

FIGURE 3

Figure 3. Pseudo-code structure illustrating the two proposed implementation variants of a floating spacecraft-manipulator system simulation model based on the Generalized Jacobian approach.

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

FIGURE 4

Figure 4. Block diagram of the overall simulated control architecture used for the numerical simulations.

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 $\overline{u}$ is computed by means of a proportional-derivative (PD) control law:

$\begin{array}{ll}\overline{u}=K_D\left({\stackrel{.}{q}}_D-\stackrel{.}{q}\right)+K_P\left(q_D-q\right)& (112)\end{array}$

In general, the controller gains can be designed based on the solutions of the harmonic oscillator:

$\begin{array}{ll}{K}_{P_i}=\frac{{\tau }_{max,i}}{{\left(q_D-q\right)}_{max}},& (113)\end{array}$

$\begin{array}{ll}{K}_{D_i}=\sqrt{2}{K}_{P_i},& (114)\end{array}$

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, ${\widetilde{\mathit{\text{H}}}}^{\star }$ and ${\widetilde{\mathit{\text{C}}}}^{\star }$:

$\begin{array}{ll}\tau ={\tilde{H}}^{\star }\overline{u}+{\tilde{C}}^{\star }(q,q\stackrel{.}{ }).& (115)\end{array}$

In the numerical simulation presented here, perfect knowledge of the inertia properties of the system was assumed, thus ${\widetilde{\mathit{\text{H}}}}^{\star }={\mathit{\text{H}}}^{\star }$ and ${\widetilde{\mathit{\text{C}}}}^{\star }={\mathit{\text{C}}}^{\star }$.

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_Di = 1 and K_Pi = 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.

TABLE 3

Table 3. Planar spacecraft-manipulator system with four-link manipulator: mass and inertia properties.

TABLE 4

Table 4. Planar spacecraft-manipulator system with four-link manipulator: customized Denavit-Hartenberg parameters.

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.

FIGURE 5

Figure 5. Planar spacecraft-manipulator system with four-link manipulator: (A) joint rates and the resulting base-spacecraft angular velocity components; (B) joint angles and the resulting base-spacecraft position deviation.

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.

FIGURE 6

Figure 6. Planar spacecraft-manipulator system with four-link manipulator: (A) linear momentum; (B) angular momentum.

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_Di = 1 and K_Pi = 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.

TABLE 5

Table 5. Spatial spacecraft-manipulator system with six-link manipulator: mass and inertia properties (Yoshida and Umetani, 1993).

TABLE 6

Table 6. Spatial spacecraft-manipulator system with six-link manipulator: customized Denavit-Hartenberg parameters.

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.

FIGURE 7

Figure 7. Spatial spacecraft-manipulator system with six-link manipulator: (A) joint rates and the resulting base-spacecraft angular velocity components; (B) joint angles and the resulting base-spacecraft position deviation.

FIGURE 8

Figure 8. Spatial spacecraft-manipulator system with six-link manipulator: (A) linear momentum; (B) angular momentum.

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 expanded by the introduction of five categories (namely floating, rotation-floating, rotation-flying, translation-flying, and flying). The authors believe that the adoption of the newly proposed definitions would contribute to increased clarity, advantageous to both students and researchers in the field. Furthermore, the paper introduces a full-symbolic implementation and a partial-symbolic implementation of computer simulation models based on the complete generalized Jacobian approach. These implementation options allow the development of efficient numerical simulations. Example simulations, for 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 spacecraft-manipulator 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.

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

This research received no specific grant from any funding agency in the public, commercial, or not-for-profit sectors.

References

Aghili, F. (2009). “Coordination control of a free-flying manipulator and its base attitude to capture and detumble a cooperative satellite,” in Proceedings of 2009 IEEE/RSJ International Conference on Intelligent Robots and Systems (St. Louis, MO), 2365–2372.

Google Scholar

Bain, M. E. (2010). “Cygnus: back to the future - applying commercial program lessons learned,” in AIAA SPACE 2010 Conference & Exposition (Anaheim, CA: AIAA), 1–8.

Google Scholar

Baksalary, J. K., and Styan, G. P. H. (2002). Generalized inverses of partitioned matrices in Banachiewicz-Schur form. Linear Algeb. Appl. 354, 41–47. doi: 10.1016/S0024-3795(02)00334-8

CrossRef Full Text | Google Scholar

Barnhart, D., Sullivan, B., Hunter, R., Bruhn, J., Fowler, E., Hoag, L., et al. (2013). “Phoenix project status 2013,” in AIAA Space 2012 Conference and Exhibition (San Diego, CA: AIAA), 1–17.

Caccavale, F., and Siciliano, B. (2001). “Quaternion-based kinematic control of redundant spacecraft/manipulator systems,” in Proceedings of the 2001 IEEE International Conference on Robotics and Automation (Seoul: IEEE), 435–440.

Google Scholar

Campa, R., and Camarillo, K. (2008). “Unit quaternions: a mathematical tool for modeling, path planning and control of robot manipulators,” in Robot Manipulators, ed M. Ceccarelli (Rijeka: InTech), 23–48.

Google Scholar

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

CrossRef Full Text | Google Scholar

Cepollina, F. (2013). Robotic Refueling Mission. NASA (Accessed January 6, 2014). Available online at: http://www.nasa.gov/mission_pages/station/research/experiments/778.html.

Coleshill, E., Oshinowo, L., Rembala, R., Bina, B., Daniel, R., and Shelley, S. (2009). Dextre: improving maintenance operations on the International Space Station. Acta Astronaut. 64, 869–874. doi: 10.1016/j.actaastro.2008.11.011

CrossRef Full Text | Google Scholar

Dreyer, L. (2009). “Latest Developments on SpaceX's Falcon 1 and Falcon 9 Launch Vehicles and Dragon Spacecraft,” in IEEE Aerospace Conference (Big Sky, MT: IEEE).

Google Scholar

Dubowsky, S., and Papadopoulos, E. (1993). The kinematics, dynamics, and control of free-flying and free-floating space robotics systems. IEEE Trans. Robotics Autom. 9, 531–543.

Google Scholar

Flores-Abad, A., Ma, O., Pham, K., and Ulrich, S. (2014). A review of space robotics technologies for on-orbit servicing. Progr. Aerospace Sci. 68, 1–26. doi: 10.1016/j.paerosci.2014.03.002

CrossRef Full Text | Google Scholar

Goodman, J. L. (2006). History of space shuttle rendezvous and proximity operations. J. Spacecr. Rockets 43, 944–959. doi: 10.2514/1.19653

CrossRef Full Text | Google Scholar

Hale, W., Lane, H., Lulla, K., and Chapline, G. (2011). Wings in Orbit: Scientific and Engineering Legacies of the Space Shuttles 1971-2010. Washington, DC: NASA.

Google Scholar

Ho, J. Y. L. (1977). Direct Path method for flexible multibody spacecraft dynamics. J. Spacecr. Rockets 14, 102–110.

Google Scholar

Hughes, P. C. (1979). Dynamics of a chain of flexible bodies. J. Astronautical Sci. 27, 359–380.

Hughes, P. C. (2004). Spacecraft Attitude Dynamics. Mineola, NY: Dover Publications.

Google Scholar

Kennedy, F. G. (2008). “Orbital express: accomplishments and lessons learned,” in Proceedings of the AAS Guidance and Control Conference (Breckenridge, CO: Univelt), 575–586.

Google Scholar

Kuipers, J. B. (2000). “Quaternions and Rotation Sequences,” in Proceedings of the International Conference on Geometry, Integrability and Quantization (Sofia: Coral Press Scientific Publishing), 127–143.

Google Scholar

Liang, B., Xu, Y., and Bergerman, M. (1998). Mapping a space manipulator to a dynamically equivalent manipulator. J. Dyn. Syst. Meas. Control 120, 1–7.

Google Scholar

Longman, R. W., Lindberg, R. E., and Zedd, M. F. (1987). Satellite-mounted manipulators - new kinematics and reaction moment compensation. Int. J. Rob. Res. 6, 87–103.

Google Scholar

Marchesi, M. (1997). “Control strategy of a free-flying space manipulator”, in Proceedings of the 8th International Conference on Advanced Robotics (ICAR '97) (Monterey, CA), 665–670

Google Scholar

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

CrossRef Full Text | Google Scholar

Moosavian, S. A.A., and Papadopulos, E. (2007). Free-flying robots in space: an overview of dynamics modeling, planning and control. Robotica 25, 537–547. doi: 10.1017/S0263574707003438

CrossRef Full Text | Google Scholar

Mukherjee, R., and Nakamura, Y. (1992). Formulation of efficient computation of inverse dynamics of space robots. IEEE Trans. Robot. Autom. 8, 400–406.

Google Scholar

Nakamura, Y., and Mukherjee, R. (1989). “Nonholonomic path planning of space robots,” in Proceedings of 1989 International Conference on Robotics and Automation (Scottsdale, AZ: IEEE), 1050–1055.

Google Scholar

Nanos, K., and Papadopoulos, E. (2011). On the use of free-floating space robots in the presence of angular momentum. Intell. Serv. Robot. 4, 3–15. doi: 10.1007/s11370-010-0083-2

CrossRef Full Text | Google Scholar

Nenchev, D. N., and Yoshida, K. (1999). Impact analysis and post-impact motion control issues of a free-floating space robot subject to a force impulse. IEEE Trans. Robot. Autom. 15, 548–557.

Google Scholar

Nenchev, D. N., Yoshida, K., Vichitkulsawat, P., and Uchiyama, M. (1999). Reaction null-space control of flexible structure mounted manipulator systems. IEEE Trans. Robot. Autom. 15, 1011–1023.

Google Scholar

Oda, M. (2000). “Experiences and lessons learned from the ETS-VII robot satellite,” in Proceedings of the 2000 IEEE International Conference on Robotics & Automation (San Francisco, CA: IEEE), 914–919.

Google Scholar

Papadopoulos, E. G. (1990). On the Dynamics and Control of Space Manipulators. Ph.D. Dissertation, Cambridge, MA: Massachusetts Institute of Technology.

Google Scholar

Patten, L., Evansa, L., Oshinowo, L., Ochisor, M., Kazuharu, N., Lodewijk, A., et al. (2002). “International Space Station Robotics: A Comparative Study of ERA, JEMRMS and MSS,” in 7th ESA Workshop on Advanced Space Technologies for Robotics and Automation ‘ASTRA 2002' (Noordwijk: ESA), 1–8.

Google Scholar

Reed, B. B., Smith, R. C., Naasz, B., Pellegrino, J., and Bacon, C. (2016). “The restore-L servicing mission,” in AIAA Space Forum (Long Beach, CA: AIAA), 1–8.

Google Scholar

Reintsema, D., Thaeter, J., Rathke, A., Naumann, W., Rank, P., and Sommer, J. (2010). “DEOS - the german robotics approach to secure and de-orbit malfunctioned satellites from low-earth orbits,” in Proceedings of International Symposium on Artificial Intelligence, Robotics and Automation in Space (i-SAIRAS) (Sapporo,: JAXA), 244–251.

Google Scholar

Roesler, G., Jaffe, P., and Henshaw, G. (2017). Orbital mechanics. IEEE Spectrum 54, 44–50. doi: 10.1109/MSPEC.2017.7864756

CrossRef Full Text | Google Scholar

Sallaberger, C. (1997). Canadian space robotic activities. Acta Astronaut. 41, 239–246.

Google Scholar

Sargent, D. G. (1984). “The Impact of Remote Manipulator Structurel Dynamics on Shuttle On-Orbit Flight Control,” in 17th Fluid Dynamics, Plasma Dynamics, and Lasers Conference (Seattle, WA: AIAA).

Google Scholar

Sato, N., and Wakabayashi, Y. (2001). “JEMRMS design features and topics from testing,” in Proceeding of the 6th International Symposium on Artificial Intelligence and Robotics & Automation in Space: i-SAIRAS 2001 (St-Hubert, QC: CSA).

Google Scholar

Shoemaker, J., and Wright, M. (2004). “Orbital express space operations architecture program,” in Proceedings of SPIE 5088, Space Systems Technology and Operations (Orlando, FL), 56–65.

Google Scholar

Siciliano, B., Sciavicco, L., Villani, L., and Oriolo, G. (2010). Robotics. London: Springer-Verlag.

Google Scholar

Stieber, M. E., Trudel, C. P., and Hunter, D. G. (1997). “Robotic systems for the international space station,” in Proceedings of the 1997 IEEE International Conference on Robotics and Automation (Albuquerque, NM: IEEE), 3068–3073.

Google Scholar

Stoneking, E. (2007). “Newton-euler dynamic equations of motion for a multi-body spacecraft,” in AIAA Guidance, Navigation and Control Conference and Exhibit, Guidance, Navigation, and Control and Co-located Conferences (Hilton Head, SC: AIAA), 1–13.

Google Scholar

Ueda, S., Kasai, T., and Uematsu, H. (2010). “HTV rendezvous technique and GN&C design evaluation based on 1st Flight on-orbit operation result,” in AIAA/AAS Astrodynamics Specialist Conference, Guidance, Navigation, and Control and Co-located Conferences (Toronto, ON: AIAA), 1–12.

Umetani, Y., and Yoshida, K. (1989). Resolved motion rate control of space manipulators with generalized Jacobian Matrix. IEEE Trans. Robotics Autom. 5, 303–314.

Google Scholar

Vafa, Z., and Dubowsky, S. (1987). “On the dynamics of manipulators in space using the virtual manipulator approach,” in IEEE International Conference on Robotics and Automation (Raleigh, NC).

Google Scholar

Wei, R., Jin, M. H., Xia, J. J., Xie, Z. W., Shi, J. X., and Liu, H. (2006). “High fidelity distributed hardware-in-the-loop simulation for space robot,” in Proceedings of the 2006 IEEE International Conference on Mechatronics and Automation (Luoyang: IEEE), 2150–2155.

Google Scholar

Wertz, J. R. (1978). Spacecraft Attitude Determination and Control. Boston, MA: Kluwer Academic Publishers.

Google Scholar

Wie, B. (2008). Space Vehicles Dynamics and Control 2nd Edn. Reston, VA: AIAA.

Xu, W., Liu, Y., Liang, B., Xu, Y., Li, C., and Qiang, W. (2008). Non-holonomic path planning of a free-floating space robotic system using genetic algorithms. Adv. Robot. 22, 451–476. doi: 10.1163/156855308X294680

CrossRef Full Text | Google Scholar

Yoshida, K. (2003). Engineering test satellite VII flight experiments for space robot dynamics and control: theories on laboratory test beds ten years ago, now in orbit. Int. J. Rob. Res. 22, 321–335. doi: 10.1177/0278364903022005003

CrossRef Full Text | Google Scholar

Yoshida, K., and Nakanishi, H. (2003). “Impedance matching in capturing a satellite by a space robot,” in Proceedings of the 2003 IEEE/RSJ Intl. Conference on Intelligent Robots and Systems. (Las Vegas, NV: IEEE), 3059–3064.

Google Scholar

Yoshida, K., and Nenchev, D. N. (1998). “A general formulation for under-actuated manipulators,” Eighth International Symposium on Robotics Research (Shonan).

Google Scholar

Yoshida, K., and Umetani, Y. (1993). “Control of space manipulators with generalized Jacobian,” in Space Robotics: Dynamics and Control, eds Y. Xu and T. Kanade (Boston, MA: Springer), 165–204.

Google Scholar

Nomenclature

${J}$ = Inertial Cartesian coordinate system orbiting with the spacecraft-manipulator system

${{J}}_i$ = Cartesian coordinate system attached to joint i

${{L}}_i$ = Cartesian coordinate system attached to link i

α_i = Denavit-Hartenberg parameter for rotation from link coordinate system axis z_i−1 to axis z_i about axis x_i

ω_E = angular velocity of the end-effector with respect to inertial frame ${J}$

ω_i = angular velocity of link i with respect to inertial frame ${J}$

ω₀ = base-spacecraft angular velocity with respect to inertial frame ${J}$

Ω = origin of the inertial Cartesian coordinate system ${J}$

${\left[\begin{array}{c}ϵ\&\zeta \&\eta \end{array}\right]}^T$ = Euler angles (sequence 1-2-3) expressing the orientation of the coordinate system attached to joint 1 with respect to the coordinate system at the spacecraft center-of-mass CM₀

${\left[\begin{array}{c}\varphi \&\theta \&\Psi \end{array}\right]}^T$ = Euler angles (sequence 1-2-3) expressing the orientation of the base-spacecraft (link 0) with respect to the inertial coordinate system ${J}$

ϱ_i = origin of the Cartesian coordinate system attached to joint i

τ = [N × 1] matrix of joint torques (each element τ_i expresses the scalar value of the torque component along ${\widehat{k}}_i$ due to the motor driving joint i)

θ_i = Denavit-Hartenberg parameter for rotation between link coordinate system axis x_i−1 and x_i about axis z_i−1

${ }_a^{ℬ}$ = [3 × 1] matrix of components of the vector a along the Cartesian coordinate system ${B}$

a ^× = [3 × 3] skew-symmetric matrix associated with vector a

a_i = position vector from ϱ_i to center-of-mass of link i

A = [4 × 4] Denavit-Hartenberg transformation matrix

b₀ = position vector from CM₀ to ϱ₁

b_i = position vector from CM_i to ϱ_i+1

c_i = Denavit-Hartenberg parameter for displacement between link coordinate system axis z_i−1 and axis z_i along axis x_i

c^⋆ = [N × 1] matrix of generalized non-linear terms of the equations of motion

C^⋆ = [N × 1] matrix of generalized centripetal and Coriolis terms

${\tilde{C}}^{\star }$ = [N × 1] matrix of estimated generalized centripetal and Coriolis terms

c₀ = [6 × 1] matrix of derivatives for base-spacecraft position

c_m = [N × 1] matrix of derivatives for manipulator joint angles

CM_i = center-of-mass of link i

CM_tot = spacecraft-manipulator system center-of-mass

d_i = Denavit-Hartenberg parameter for displacement between link coordinate system axis x_i−1 and axis x_i along axis z_i−1

ε = [3 × 1] vector part of unit quaternion q

η = Scalar part of unit quaternion q

𝕀_{3, 3} = [3 × 3] identity matrix

H^⋆ = [N × N] generalized inertia properties matrix

${\tilde{H}}^{\star }$ = [N × N] estimated generalized inertia properties matrix

H₀ = [6 × 6] base inertia properties matrix

H_0m = [6 × N] base/manipulator dynamic-coupling matrix

H_m = [N × N] manipulator inertia properties matrix

H_S = [3 × 3] system moments-of-inertia matrix

H_Sq = [3 × N] angular motion energy contribution matrix

${}^{{{L}}_i}I{ }_i$ = [3 × 3] moments-of-inertia matrix of link i assembly with respect to the Cartesian coordinate system attached to link i

${ }^{{J}}{I}_i$ = [3 × 3] moments-of-inertia matrix of link i assembly with respect to the inertial coordinate system ${J}$

J^⋆ = [6 × N] generalized Jacobian for the motion of the end-effector

$J_{X_i}^{\star }$ = [6 × N] generalized Jacobian for the motion of X_i

J_0Xi = [6 × 6] Jacobian for the contribution of the base-spacecraft motion to the velocity of X_i

$J_k^{\star }$ = [6 × N] generalized Jacobian for point of action of external force or torque F_k

J_mXi = [6 × N] Jacobian for the contribution of the manipulator motion to the velocity of X_i

J_Ri = [3 × N] angular motion Jacobian for the center-of-mass of link i

J_Ti = [3 × N] linear motion Jacobian for the center-of-mass of link i

j_Tik = [3 × 1] k-th column of matrix J_Ti

J_TS = [3 × N] linear motion energy contribution matrix

K_D = [N × 1] matrix of elements K_Di

K_Di = derivative controller gain

${\widehat{k}}_i$ = unit vector along the z axis of the coordinate system attached to joint i

K_P = [N × 1] matrix of elements K_Pi

K_Pi = proportional controller gain

L = Lagrangian

L = angular momentum vector about CM₀ expressed in frame ${J}$

m_i = mass of link i assembly (including connected parts of joints i and i+1)

${ }^{{J}}{M}_0$ = [6 × 1] combined matrix of initial momentum (linear and angular) with respect to the base-spacecraft center-of-mass

m_tot = mass of the spacecraft-manipulator system

N = number of manipulator joints

0_{i, j} = [i × j] zero matrix

P = linear momentum vector about CM₀ expressed in frame ${J}$

p_i = position vector from point Ω to the point ϱ_i

q = [N × 1] matrix of elements q_i

q_D = [N × 1] matrix of elements q_D

q = [4 × 1] unit quaternion

q_Di = desired angular displacement of joint i

q_i = angular displacement of joint i

${}^{{B}}q{ }_{r{A}}$ = [4 × 1] unit quaternion expressing the rotation between coordinate system ${A}$ and coordinate system ${B}$

r_C = position vector from Ω to CM_tot expressed in frame ${J}$

r_i = position vector from point Ω to point CM_i expressed in frame ${J}$

r_0C = position vector from CM₀ to CM_tot expressed in frame ${J}$

r_0i = position vector from CM₀ to CM_i expressed in frame ${J}$

${}^{{B}}R{ }_{{A}}$ = [3 × 3] direction-cosine matrix from the coordinate system ${A}$ to the coordinate system ${B}$

S_U = [3 × 3] Schur complement matrix of upper left sub-matrix of H_b

s_{j, k} = row j, column k element of S_U

${}^{{B}}t{ }_{{A}}$ = [4 × 1] four-vector expressing the translation of the 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₀ = 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_Xi = linear velocity of arbitrary point X_i located on link i of the manipulator with respect to the inertial frame ${J}$

${ }^{{J}}{x}_0$ = [6 × 1] combined position matrix of linear position components and attitude parameters of the base-spacecraft, (${ }^{{J}}{x}_b$ = ${\left[\begin{array}{cccc}{ }^{{J}}{r}_0{ }^T& \varphi & \theta & \psi \end{array}\right]}^T$)

${ }^{{J}}{\dot{x}}_0$ = [6 × 1] combined velocity matrix of linear velocity components and angular velocity components of the base-spacecraft,

$\begin{array}{l}{(}^{{J}}{\stackrel{.}{x}}_0={\left[\begin{array}{cc}{ }^{{J}}{v}_0{ }^T& { }^{{J}}{\omega }_0{ }^T\end{array}\right]}^T={\left[\begin{array}{cc}{ }^{{J}}{v}_0{ }^T& {\left(\left[\begin{array}{ccc}1& 0& -\text{ sin }\theta \\ 0& \text{ cos }\varphi & \text{ cos }\theta \text{ sin }\varphi \\ 0& -\text{ sin }\varphi & \text{ cos }\varphi \text{ cos }\theta \end{array}\right]\left[\begin{array}{c}\dot{\varphi }\\ \dot{\theta }\\ \dot{\psi }\end{array}\right]\right)}^T\end{array}\right]}^T)\end{array}$

x_0i = position vector from CM₀ to an arbitrary point on link i

X_i = arbitrary point on link i

x_i = position vector from Ω to an arbitrary point on link i

DCM = Direction-Cosine Matrix

DEM = Dynamically Equivalent Manipulator

DH = Denavit-Hartenberg

DOF = Degree of Freedom

EE = End-Effector

ETS-VII = Engineering Test Satellite VII

GJM = Generalized Jacobian Matrix

HTV = H-II Transfer Vehicle

ISS = International Space Station

JEMRMS = Japanese Experiment Module Remote Manipulator System

PD = proportional-derivative controller

SPDM = Special Purpose Dexterous Manipulator

SRMS = Shuttle Remote Manipulator System

VM = Virtual Manipulato.