# Dynamics and robust control of a space manipulator system servicing a satellite with flexible appendages

Sofiane Kraïem, Mathieu Rognant, Sérgio Waitman, Idriss Chelikh, Jean-Marc Biannic

► To cite this version:

Sofiane Kraïem, Mathieu Rognant, Sérgio Waitman, Idriss Chelikh, Jean-Marc Biannic. Dynamics and robust control of a space manipulator system servicing a satellite with flexible appendages. IFAC Journal of Systems and Control, 2025, 33, 10.1016/j.ifacsc.2025.100328. hal-05214757

HAL Id: hal-05214757

https://hal.science/hal-05214757v1

Submitted on 19 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Full length article

# Dynamics and robust control of a space manipulator system servicing a satellite with flexible appendages

Sofiane Kraïem *, Mathieu Rognant, Sérgio Waitman, Idriss Chelikh, Jean-Marc Biannic

ONERA - The French Aerospace Lab, Information Processing and Systems Department, 31000, Toulouse, France

## ARTICLE INFO

**Article history:**

Received 19 February 2024

Received in revised form 7 July 2025

Accepted 14 July 2025

Available online 22 July 2025

**Keywords:**

Space manipulator system

Modeling

Flexible structures

Robust control

Nonlinear Dynamic Inversion

Disturbance compensation

## ABSTRACT

Space Manipulator Systems (SMS) are becoming pivotal in space exploitation and exploration, offering a versatile range of solutions from space debris capture to structure assembly. However, recent missions involving manipulators aboard satellites and space structures must contend with lightweight and large elements that exhibit flexible behaviors. Despite the challenges posed by flexible elements in their surroundings, enhancing the autonomy of SMS remains crucial to ensure their viability as solutions. For the pre-design of the SMS, path-planning applications, or controller design, there is a necessity for methods to assess the couplings between the manipulator, the SMS platform, and any flexible elements moved by a manipulator or attached to the platform. Moreover, recently proposed control strategies have demonstrated a keen interest in developing model-based controllers, which advantageously provide an efficient utilization of actuators and mitigation of internal disturbances within the system. This paper initially presents, for control purposes, the derivation of the kinematics and dynamics of a free-floating Space Manipulator System (SMS) with a flexible body at the end of a kinematic chain, employing a Lagrangian formalism. Subsequently, a robust joint-space control law is designed for the On-Orbit Servicing (OOS) of a satellite with flexible appendages. The control structure consists of Nonlinear Dynamic Inversion (NDI) for system linearization and decoupling. A structured $H_\infty$ controller synthesis is developed to provide robustness against flexible disturbances, model uncertainties, and sensor noise. To evaluate the effectiveness of this control strategy, it is implemented on a real-time simulation platform that ensures tight and high-fidelity space robot dynamics and flexible structures. This platform incorporates visual environment models and virtual sensors to provide an accurate representation of real-world conditions. The considered use-case involves the servicing of a satellite equipped with a flexible solar array using a free-floating dual-arm SMS with flexible appendages.

© 2025 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY license

(http://creativecommons.org/licenses/by/4.0/).

## 1. Introduction

Space Manipulator Systems (SMS) are predicted a bright future to address a wide array of challenges in the space exploitation and exploration. In particular SMS provide advantageous and versatile solutions to reduce the number of humanized flight missions. Active Debris Removal (ADR) or On-Orbit Servicing (OOS) missions, encompassing tasks like refueling, maintenance and deployment of structural components, are prime examples that stand to benefit significantly from robotic systems (Flores-Abad et al., 2014; Papadopoulos et al., 2022). Over the years, the Canadian Space Agency (CSA) has made substantial contributions to the development of robotic technologies. For in-space assembly, the first version of the Canadarm on the space shuttle played a key role in the installation of the MIR space station's docking

module in 1995 (Hiltz et al., 2001). The subsequent version, Canadarm2, is still employed for assembly tasks, thereby extending the scope of robotic applications. Likewise, the National Space Development Agency of Japan (NASDA) were pioneers in the development of autonomous control technologies for SMS. The Engineering Test Satellite No. 7 (ETS VII) in 1997 conducted an initial demonstration mission without a crew, enabling the experimental validation of technologies for both tele-operated and autonomous manipulator tasks (Oda, 2000). However, SMS have predominantly been tele-operated, either from the ground, the space shuttle, or a space station. This poses a significant challenge in ensuring the safe and precise operation of a manipulator to meet future requirements (Hirzinger et al., 2004; Penin, 2002). Democratizing the use of SMS will hinge on enhancing their autonomy. However, the requirements for highly precise and autonomous systems create multiple challenges (Flores-Abad et al., 2014; Papadopoulos et al., 2022).

The principal challenge when controlling space robot is to deal with the different couplings between manipulator(s) and

<footnote>* Corresponding author.
E-mail address: sofiane.kraiem@onera.fr (S. Kraïem).</footnote>
the spacecraft's platform. As defined by (Wilde et al., 2018), free-floating and free-flying systems can be distinguished. While free-flying SMS may compensate for manipulator motions with actuators providing external forces/torques (Fauré et al., 2022), such as thrusters, free-floating robots can only compensate for spacecraft base rotations with kinetic moment exchange devices, such as reaction-wheels or control moment gyroscopes (Wu et al., 2018). Free-floating systems will be favored for on-orbit assembly (Rognant et al., 2019) while use of thrusters will be unavoidable for capture scenarios (Fauré et al., 2022). Therefore, according to the OOS operation, different control strategies have been investigated. To mitigate the coupling effects between the manipulator and the spacecraft's base, early studies have suggested the utilization of the Generalized Jacobian Matrix (GJM) (Yoshida & Umetani, 1993) or coordinated base/manipulator control (Papadopoulos & Dubowsky, 1991). The GJM, introduced by Umetani (Umetani & Yoshida, 1989), is the matrix expressing the relation between the manipulator's end-effector and the joints in the Cartesian space. (Rybus et al., 2015) analyzed a control scheme based on the fixed-base Jacobian inverse with the inclusion of spacecraft velocity to reduce complexity associated with the GJM. The use of an GJM allows to perform manipulator motions minimizing the direct coupling between the end-effector and the base motions. (Pisculli & Gasbarri, 2015) adapted this strategy introducing a generalized mass Jacobian matrix to reduce the impact of a captured object on the servicer's base. Recently, kinematic indices derived from the GJM have been introduced to evaluate actuator participation, analyzing the feasibility of a motion (Rognant et al., 2019b). Such indices illustrated better improvement of possible motions from a coordinated control of the spacecraft base and manipulator's joint. As for capture applications, a separated control of the servicer base and the manipulator is common (Colmenarejo et al., 2018; Cumer et al., 2021), from manipulator's workspace adjustment (Giordano et al., 2017) or considering the SMS global Center of Mass (CoM) (Giordano et al., 2019) a simultaneous control of both the base and the manipulator offers efficient propellant consumption strategies. Likewise, to deal with low reaction-wheels saturation thresholds the combination with control moment gyroscopes (Wu et al., 2018) or with thrusters (Giordano et al., 2020) have been investigated.

Furthermore, a new challenge arise with the increase of light and large elements composing space structures that will introduce flexible behaviors within the overall space system. Current and future missions are involving multiple solar arrays, sunshields and antennas that are examples of sources for flexible dynamics (Rognant et al., 2019). As experienced with the Hubble telescope, in presence of flexible appendages, a poor modeling can lead to a possible loss of control and thus mission without intervention (Foster et al., 1995). Literature on flexible multi-body systems is rich (Banerjee, 2022; Bauchau, 2011; Bremer, 2008; Géradin & Cardona, 2001; Shabana, 1997; Simeon, 2013). However, introducing flexible dynamics in SMS requires relative knowledge of the flexible modes which represents a first modeling difficulty. A second difficulty is the increase of system states to account for the different flexible modes of each elements. For this reasons, recent studies have considered simple dynamic models to tackle the coupling between the manipulator and the flexible appendages prior developing manipulator control strategies while suppressing vibrations (Hirano et al., 2014; Meng et al., 2017). To overcome the fastidious modeling of flexible bodies among multi-body systems, robust control to unknown external/internal disturbances have been considered to reduce the couplings of flexible elements to the manipulator (Qiao et al., 2019; Seddaoui & Saaj, 2019). However, by leveraging the knowledge of the flexible modes of each appendage, alongside mitigating internal disturbances caused by vibrations using a model-based

observer (Kraïem et al., 2022), it is possible to achieve an effective and minimized use of the base actuators for executing manipulator tasks (Kraïem et al., 2021).

Improving SMS autonomy will be achieved by developing new model-based control strategies and an inherent effort to model new more and more complex space structure. Modeling of a space robot refers to the derivation of both dynamic and kinematic models. A kinematic model gives the relations between the manipulator's joint states and the end-effector in the satellite frame. A dynamic model provides the system's equations of motions in terms of accelerations, forces and torques. A distinctive aspect of space robots is that the inverse/forward kinematic problem is inherently dynamic (Flores-Abad et al., 2014). Focusing on the manipulation of flexible elements in a detumbling or deployment scenario, two main modeling methods can be developed to obtain a kinematic/dynamic model: the Newton-Euler and the Lagrangian formalism. The Newton-Euler approach, which is a direct interpretation of Newton's second law, formulates that the dynamic model by considering for each body of the system the force and moment applied. The recursivity of the method is a first advantage for its implementation (Carpenter et al., 2019; Virgili-Llop et al., 2016). Moreover, the independent consideration of each body allows more versatility in the modeling process and leads to the possibility of developing plug-in approaches (Alazard et al., 2008). Thus, for instance closed-kinematic chains and additional elements connected to the end-effector can be studied. The Lagrangian formalism provides a simpler and systematic modeling method based on the energies and works present in the system. Regarding energies of a rigid space robot, only the kinetic ones are considered which allows by developing a Lagrangian approach to obtain an easier method to derive a dynamic model (Wilde et al., 2018). With rigid systems and assuming no external interactions, the momentum conservation is a means of state parameters estimation to retrieve the physical properties of the space manipulator systems (Christidi-Loumpasefski et al., 2017). In this way, it is a common strategy to model a rigid multi-body system evolving in a free fall environment by initiating a Lagrangian formalism. Newton-Euler approaches have been adapted to integrate flexible dynamics as for each body a force/moment is listed so that a connection between two bodies is made with summing forces and torques on the junction. Such methods have been put forward by different authors (Rodrigues et al., 2022; Rui et al., 2014; Sanfedino et al., 2018). Nevertheless, a limit of the methods developed in the literature is the consideration of linearized model. Such models are well-suited for the development of control laws but limit the range of SMS applications. To alleviate assumptions while maintaining a high level of modeling fidelity, a Lagrangian formalism allows the inclusion of flexible appendages with an associated linear model, while accounting for nonlinear couplings between rigid and flexible elements. Such approach finds its interest to develop model-based control strategies (Kraïem et al., 2022). All modeling methods developed within the multibody context can leverage the robust and widely validated Finite Element Method (FEM), which is the most commonly used numerical technique in structural analysis. Two methods that utilize FEM results within a substructure approach are commonly employed: the Finite Element-Transfer Matrix (FE-TM) method (Leckie & Pestel, 1960) and the Component Modes Synthesis (CMS) (Rui et al., 2022). The FE-TM combines the finite element approach with transfer matrix models. By characterizing each flexible body with a transfer function between two state vectors at each side of the body and defined by the generalized displacements of the point and the corresponding generalized forces, the dynamic model of the structure is obtained by multiplication of the transfer functions of each body connected in series. A limitation of this
method lies in the derivation of the transfer matrix (Rui et al., 2022) which makes the approach less common than the CMS. The CMS approach (Craig & Bampton, 1968; Hurty, 1965) which consists in separating boundary and internal Degrees of Freedom (DoFs) has been employed to describe substructure architectures. To model multiple flexible body in chain-like configuration without losing the dynamic information corresponding to the free tip appendage points, the Two-Input Two-Output (TITOP) Port model has been introduced in (Alazard et al., 2008) and recently extended for larger angular velocities with a Lagrangian formalism (Rodrigues et al., 2025). It constructs the dynamic model of each substructure as a transfer between the accelerations and the external forces at each connection point in a state-space form.

This paper addresses the problem of integrating flexible bodies at the end of rigid multi-kinematic chain systems. After integrating rigidly connected flexible elements in the open-source spacecraft robotics toolkit (Virgili-Llop et al., 2016) in our previous study (Kraïem et al., 2022), we extended the open-source toolbox to facilitate the design of a model-based control that aims at linearizing and decoupling the system expressed in its joint-space representation. A feedback, design with an $H_\infty$ synthesis, allows to provide robustness against modeling uncertainties and sensor noises. The proposed solution's efficiency is demonstrated through the servicing of a satellite and its flexible solar array by a free-floating dual-arm SMS equipped with flexible solar arrays. Extensive testing is carried out on a real-time simulation platform that faithfully replicates realistic and high-fidelity conditions (Kraïem et al., 2024). The paper is structured as follows: Section 2 derives from a Lagrangian formalism the kinematics/dynamics of a multi-body system with flexible appendages ending kinematic chains; Section 3 develops the proposed control framework and its design; Section 4 illustrates the efficiency of the method through simulations. Finally, Section 5 concludes and states future perspectives.

## 2. Flexible-rigid dynamics of a multi-kinematic chains and multi-body space robot

In this section, we propose to develop a kinematic and dynamic model of a multi-body space robot with flexible bodies positioned at the end of rigid kinematic chains. This model aims at providing the principal couplings between flexible appendages and the rest of the SMS for controller synthesis, workspace analysis and path-planning applications. The integration of flexible dynamics is performed with a Lagrangian approach after a Finite Element Method providing the main flexible modes. This allows to consider couplings due to Coriolis and Centrifugal effects induced by flexible bodies on the structure (Rodrigues et al., 2025). Modeling assumptions consider in this approach are: rigid manipulators, small velocities of the system, orbital motions are neglected as well as external disturbances to focus on the impact of internal couplings. However, such assumptions are aligned with the considerations of On-Orbit Deployment of large and light structure, detumbling of a satellite with flexible appendages.

### 2.1. Introduction of the notations with Denavit-Hartenberg formalism

A kinematic chain comprises a series of rigid bodies interconnected by joints with one-DoF, such as a revolute or prismatic joint, or a fixed connection. The Denavit-Hartenberg (DH) convention (Denavit & Hartenberg, 1955) enables the characterization of the relative position of two solids using only four parameters. These parameters are defined for each body with a specific frame associated with it, as visualized in Fig. 1. Hence, a recursive formulation of the kinematic model can be established.

<figure><img src="image_2.png" /><figcaption>Fig. 1. Illustration of Denavit-Hartenberg parameters.</figcaption></figure>

Before introducing the method for identifying DH parameters, we will define the construction of the multi-kinematic chain system using the following notations. The common base for all kinematic chains is defined as the first rigid solid of the system and referred as $S_0$. From the robot's base, $n_q$ additional solids (or links) are considered, with $n_q \in \mathbb{N}$. It is assumed that the last solid of a unique kinematic chain system is the $n_{EE}^{th}$. To facilitate readability, these notations are consistently retained for multiple kinematic chains with more than one link. The subscript $EE_{c_i}$ will refer to the last solid of the $i$th kinematic chain that has more than one-DoF. The notation $S_i$ is used to refer at the $i$th solid, with $i \in [1, n_q]$. Between the bodies $S_{i-1}$ and $S_i$ a connecting joint, denoted $A_i$, is either considered fixed, prismatic or revolute. Moreover, for each joint and solid, a Cartesian body frame, $R_{i-1} = (O_{i-1}, \vec{x}_{i-1}, \vec{y}_{i-1}, \vec{z}_{i-1})$, is associated in respect to the following steps:

* The axis $\vec{z}_{i-1}$ is defined along the axis of joint $A_i$

* The origin $O_{i-1}$ corresponds to the intersection between $z_{i-1}$ and the normal to $z_{i-2}$ and $z_{i-1}$ that carries $\vec{x}_{i-1}$

* The vector $\vec{y}_{i-1}$ is deduced from the following cross product $\vec{x}_{i-1} = \vec{z}_{i-1} \times \vec{x}_{i-1}$

* The first Cartesian frame $R_0 = (O_0, \vec{x}_0, \vec{y}_0, \vec{z}_0) = R_{sat}$ is defined such that $\vec{z}_0$ is along the axis of $A_1$, $\vec{x}_0$ and $\vec{y}_0$ are arbitrarily chosen

* The last frame $R_{n_{EE}} = (O_{n_{EE}}, \vec{x}_{n_{EE}}, \vec{y}_{n_{EE}}, \vec{z}_{n_{EE}})$ is defined such that $\vec{x}_{n_{EE}}$ is normal to the axis of $A_{n_{EE}}$, $\vec{z}_{n_{EE}}$ and $\vec{y}_{n_{EE}}$ are arbitrarily chosen

Then the four DH parameters are defined with these referential frames such as:

* $d_i$ is the distance between $x_{i-1}$ and $x_i$ along axis $\vec{z}_{i-1}$, it corresponds to a prismatic DoF

* $\theta_i$ is the rotation around $\vec{z}_i$ between $\vec{x}_{i-1}$ et $\vec{x}_i$, it corresponds to a rotational DoF

* $a_i$ is the distance between $z_{i-1}$ and $z_i$ along axis $\vec{x}_i$, it corresponds to the common normal of $\vec{z}_i$ and $\vec{z}_{i-1}$

* $\alpha_i$ is the rotation around $\vec{x}_i$ between $\vec{z}_{i-1}$ and $\vec{z}_i$

Moreover, the geometry of the space robot with respect to an inertial frame, $R_{ine}$, can be described with the position vector, $\mathbf{r}_{A_i} \in \mathbb{R}^{3 \times 1}$, of each joint $A_i$ and the position vector of each Center-of-Mass (CoM), $\mathbf{r}_{S_i} \in \mathbb{R}^{3 \times 1}$, of link $S_i$ ($i \in [0, n_q]$). With the DH parameters, homogeneous transformations can be defined
to recursively express the pose and orientation of the ith link or joint. For this purpose, the matrix $\mathbf{T}_{A_{i-1}, A_i}$ is introduced as a homogeneous transformation matrix between the frames $\mathcal{R}_{i-1}$ and $\mathcal{R}_i$, which applies a rotation and a translation from the first frame to the second. The matrix $\mathbf{T}_{\text{trans}}^{\chi}(d)$ denotes a translation matrix along the $\bar{x}$ axis by a distance $d$, and $\mathbf{T}_{\text{rot}}^{\chi}(\theta)$ represents a rotational matrix around the $\bar{x}$ axis by an angle $\theta$. The matrix $\mathbf{T}_{A_{i-1}, A_i}$ is then defined with the DH parameters as follows:

$$ \mathbf{T}_{A_{i-1}, A_i} = \mathbf{T}_{\text{trans}}^{z_{i-1}}(d_i) \mathbf{T}_{\text{rot}}^{z_{i-1}}(\theta_i) \mathbf{T}_{\text{trans}}^{\chi}(a_i) \mathbf{T}_{\text{rot}}^{\chi}(\alpha_i) \quad (1) $$

where:

$$ \mathbf{T}_{\text{trans}}^{z_{i-1}}(d_i) = \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & d_i \\ 0 & 0 & 0 & 1 \end{bmatrix} \quad (2a) $$

$$ \mathbf{T}_{\text{rot}}^{z_{i-1}}(\theta_i) = \begin{bmatrix} c(\theta_i) & -s(\theta_i) & 0 & 0 \\ s(\theta_i) & c(\theta_i) & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix} \quad (2b) $$

$$ \mathbf{T}_{\text{trans}}^{\chi}(a_i) = \begin{bmatrix} 1 & 0 & 0 & a_i \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix} \quad (2c) $$

$$ \mathbf{T}_{\text{rot}}^{\chi}(\alpha_i) = \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & c(\alpha_i) & -s(\alpha_i) & 0 \\ 0 & s(\alpha_i) & c(\alpha_i) & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix} \quad (2d) $$

with $c(\cdot)$ and $s(\cdot)$ representing the cosine and sine functions, respectively. The homogeneous transformation matrix between $\mathcal{A}_{i-1}$ and $\mathcal{A}_i$ is then given by:

$$ \mathbf{T}_{A_{i-1}, i} = \begin{bmatrix} c(\theta_i) & -c(\alpha_i)s(\theta_i) & s(\alpha_i)s(\theta_i) & a_i c(\theta_i) \\ s(\theta_i) & c(\theta_i)c(\alpha_i) & -c(\theta_i)s(\alpha_i) & a_i s(\theta_i) \\ 0 & s(\alpha_i) & c(\alpha_i) & d_i \\ 0 & 0 & 0 & 1 \end{bmatrix} \quad (3) \\ = \begin{bmatrix} \mathbf{R}_{A_i} & (\mathbf{r}_{A_i} - \mathbf{r}_{A_{i-1}}) \\ \mathbf{0}_{1 \times 3} & 1 \end{bmatrix} $$

It is composed of a translation defined by the position, $\mathbf{r}_{A_i}$, of the joint $\mathcal{A}_i$ in $\mathcal{R}_{ine}$, and a rotation given by the Direct Cosine Matrix (DCM), $\mathbf{R}_{\mathcal{A}_i} \in \mathbb{R}^{3 \times 3}$, between $\mathcal{A}_{i-1}$ and $\mathcal{A}_i$ that specifies the direction of the rotation of $\mathcal{A}_i$ relative to the previous joint. A similar operation can be defined to describe the homogeneous transformation between the links $S_{i-1}$ and $S_i$. The frame $\mathcal{R}_{S_{i-1}}$ attached to the solid $S_{i-1}$ is defined such as its origin $O_{S_{i-1}}$ corresponds to the CoM of $S_{i-1}$, $\mathbf{r}_{S_i}$, and its axis are collinear to the ones defining the frame $\mathcal{R}_{A_i}$. The origin of the Cartesian base frame $\mathcal{R}_{sat}$ is chosen such as $\mathbf{r}_{S_i} = \mathbf{r}_{A_0}$. Therefore, the homogeneous transformation matrix between the links $S_{i-1}$ and $S_i$ is given with the DH parameters and (3) as:

$$ \mathbf{T}_{S_{i-1}, S_i} = \begin{bmatrix} \mathbf{R}_{S_i} & (\mathbf{r}_{S_i} - \mathbf{r}_{S_{i-1}}) \\ \mathbf{0}_{1 \times 3} & 1 \end{bmatrix} \quad (4) $$

The rotations of the rigid link $S_i$ and the joint $\mathcal{A}_{i+1}$ are identical which leads to the equality $\mathbf{R}_{S_i} = \mathbf{R}_{A_{i+1}}$. Finally, the homogeneous transformation matrix from $S_i$ to $\mathcal{A}_{i+1}$, $\mathbf{T}_{S_i, \mathcal{A}_{i+1}}$, is given by:

$$ \mathbf{T}_{S_i, \mathcal{A}_{i+1}} = \begin{bmatrix} \mathbf{R}_{\mathcal{A}_{i+1}} & (\mathbf{r}_{\mathcal{A}_{i+1}} - \mathbf{r}_{S_i}) \\ \mathbf{0}_{1 \times 3} & 1 \end{bmatrix} \quad (5) $$

Moreover, to express the homogeneous transformation from $\mathcal{A}_i$ to $\mathcal{A}_{i-1}$ in a decreasing recursive order, the following relation is provided:

$$ \mathbf{T}_{\mathcal{A}_{i}, \mathcal{A}_{i-1}} = \mathbf{T}_{\mathcal{A}_{i-1}, \mathcal{A}_{i}}^{-1} = \begin{bmatrix} \mathbf{R}_{\mathcal{A}_{i}}^{T} & -\mathbf{R}_{\mathcal{A}_{i}}^{T} (\mathbf{r}_{\mathcal{A}_{i}} - \mathbf{r}_{\mathcal{A}_{i-1}}) \\ \mathbf{0}_{1 \times 3} & 1 \end{bmatrix} \quad (6) $$

For control purposes, the DCM of the robot's base can be expressed as a function of Euler angles, $[\phi \ \theta \ \psi]^T$ (roll, pitch, yaw $\psi$). With the following rotation sequence $z-y-x$, it is written as:

$$ \begin{aligned} \mathbf{R}_{S_0} &= \mathbf{R}_{A_0} \\ &= \begin{bmatrix} c(\psi) & -s(\psi) & 0 \\ s(\psi) & c(\psi) & 0 \\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} c(\theta) & 0 & s(\theta) \\ 0 & 1 & 0 \\ -s(\theta) & 0 & c(\theta) \end{bmatrix} \begin{bmatrix} 1 & 0 & 0 \\ 0 & c(\phi) & -s(\phi) \\ 0 & s(\phi) & c(\phi) \end{bmatrix} \\ &= \begin{bmatrix} c(\psi)c(\theta) & c(\psi)s(\theta)c(\phi) - s(\psi)c(\theta)c(\phi) \\ s(\psi)c(\theta) & c(\psi)c(\phi) + s(\psi) \\ -s(\theta) & c(\theta)s(\phi) \end{bmatrix} \\ & \qquad \begin{bmatrix} s(\psi)s(\phi) + c(\psi)s(\theta)c(\phi) \\ -c(\psi)s(\phi) + s(\psi)s(\theta)c(\phi) \\ c(\theta)c(\phi) \end{bmatrix} \end{aligned} $$

The interest of using this formalism relies in the recursivity to express a body state. Denoting $\mathbf{T}_{\mathcal{A}_i, \mathcal{R}_{ine}}$ and $\mathbf{T}_{S_i, \mathcal{R}_{ine}}$ respectively the transformation of joint $\mathcal{A}_i$ and body $S_i$ expressed in $\mathcal{R}_{ine}$, a recursive transformation is given by:

$$ \mathbf{T}_{A_0, \mathcal{R}_{ine}} = \mathbf{T}_{S_0, \mathcal{R}_{ine}} \quad (7a) $$

$$ \begin{aligned} \mathbf{T}_{A_1, \mathcal{R}_{ine}} &= \mathbf{T}_{S_0, \mathcal{R}_{ine}} \mathbf{T}_{A_1, S_0} \\ &= \mathbf{T}_{A_0, \mathcal{R}_{ine}} \mathbf{T}_{A_1, A_0} \end{aligned} \quad (7b) $$

$$ \begin{cases} \mathbf{T}_{\mathcal{A}_i, \mathcal{R}_{ine}} = \mathbf{T}_{\mathcal{A}_{i-1}, \mathcal{R}_{ine}} \mathbf{T}_{\mathcal{A}_{i-1}, \mathcal{A}_i} \\ = \mathbf{T}_{A_1, \mathcal{R}_{ine}} \prod_{j=2}^{i} \mathbf{T}_{\mathcal{A}_{j-1}, \mathcal{A}_j}, \forall i \in [2, n_q] \end{cases} \quad (7c) $$

$$ \mathbf{T}_{S_i, \mathcal{R}_{ine}} = \mathbf{T}_{\mathcal{A}_i, \mathcal{R}_{ine}} \mathbf{T}_{\mathcal{A}_i, S_i}, \forall i \in [2, n_q] \quad (7d) $$

One will note that the transformation $\mathbf{T}_{\mathcal{A}_i, \mathcal{R}_{ine}}$ is composed of the overall DCM of the joint $\mathcal{A}_i$ expressed in the $\mathcal{R}_{ine}$, $\mathbf{R}_{\mathcal{A}_i, \mathcal{R}_{ine}}$, and its position given in $\mathcal{R}_{ine}$. This DCM is used to derive the direction of the joint rotation, $\mathbf{k}_i$, in the inertial frame as:

$$ \mathbf{k}_i = \mathbf{R}_{\mathcal{A}_i, \mathcal{R}_{ine}} \begin{bmatrix} 0 \\ 0 \\ 1 \end{bmatrix} = \mathbf{R}_{A_0} \prod_{j=1}^{i} \mathbf{R}_{\mathcal{A}_j} \begin{bmatrix} 0 \\ 0 \\ 1 \end{bmatrix} \quad (8) $$

To construct the kinematic model, a mass $m_i$ ($i \in [0, n_q]$) and an associated inertia $\mathcal{I}_{S_i}$ ($i \in [0, n_q]$) expressed in the corresponding body frame are defined for each body. With (8), the inertia matrix expressed in $\mathcal{R}_{ine}$ and denoted $\mathcal{I}_i$ is written as:

$$ \mathcal{I}_i = \mathbf{R}_{S_i, \mathcal{R}_{ine}} \mathcal{I}_{S_i} \mathbf{R}_{S_i, \mathcal{R}_{ine}}^T \quad (9) $$

Denoting $m_{tot}$ as the total mass of the system, the position of the CoM in $\mathcal{R}_{ine}$ is expressed as a function of each body CoM positions, $\mathbf{r}_{S_i}$, as follows:

$$ \mathbf{r}_{\text{CoM}} = \frac{\sum_{i=0}^{n_q} m_i \mathbf{r}_{S_i}}{\sum_{i=0}^{n_q} m_i} = \frac{\sum_{i=0}^{n_q} m_i \mathbf{r}_{S_i}}{m_{\text{tot}}} \quad (10) $$

And a recursive expression of the links CoM can be developed as:

$$ \mathbf{r}_{S_i} = \mathbf{r}_{S_0} + \sum_{j=1}^{i} (\mathbf{r}_{S_i} - \mathbf{r}_{S_{i-1}}), \forall i \in [0, n_q] \quad (11) $$

With the DH formalism introduced in this section, expressions for the pose and orientation of each element composing a multi-body system with one or more kinematic chains can be obtained recursively. The homogeneous transformations between two solids/joints are determined based on the system's DoFs and more generally, on the nature of the junctions.
## 2.2. Velocities and accelerations

To describe system velocities, twist vectors are introduced. The twist vector, $\mathbf{t}_i$, of the ith body gathers the angular and linear velocities, both expressed in $\mathcal{R}_{ine}$ and denoted respectively as $\omega_i \in \mathbb{R}^{3 \times 1}$ and $\dot{\mathbf{r}}_{S_i} \in \mathbb{R}^{3 \times 1}$. Therefore, for all $i \in [0, n_q]$, the twist of $S_i$ is given as:

$$ \mathbf{t}_i = \begin{bmatrix} \omega_i \\ \dot{\mathbf{r}}_{S_i} \end{bmatrix}, \quad \forall i \in [0, n_q] \qquad (12) $$

The kinematic model is obtained by expressing the twist vectors in function of the joint-space velocities. We introduce a generalized joint rate-variables vector, $\dot{\mathbf{q}} = [\dot{\mathbf{q}}_1^T \dots \dot{\mathbf{q}}_{n_q}^T]^T \in \mathbb{R}^{n_q \times 1}$, such that $\forall i \in [1, n_q]$, $\dot{\mathbf{q}}_i \in \mathbb{R}$ represents the velocity state of the ith DoF in the kinematic chain(s), and $\dot{\mathbf{q}}_0 \in \mathbb{R}^{6 \times 1}$ denotes the base linear and angular velocity joint state. With the DH parameters, one can express $\forall i \in [1, n_q]$, the velocity $\dot{\mathbf{q}}_i$ as:

* $\dot{\mathbf{q}}_i = \dot{\theta}_i$ if $A_i$ is a rotational joint

* $\dot{\mathbf{q}}_i = \dot{d}_i$ if $A_i$ is a prismatic joint

The joint-space is defined in function of the nature of the junctions and thus the system's DoFs. For control applications, the spacecraft angular velocity is preferably expressed in its body frame, $\mathcal{R}_{sat}$. Defining $\dot{\mathbf{q}}_0$ such as $\dot{\mathbf{q}}_0 = [\omega_0^{sat} \ \dot{\mathbf{r}}_{S_0}^T]^T$, with the angular velocity expressed in $\mathcal{R}_{sat}$ and linear velocity in $\mathcal{R}_{ine}$. Its expression allows to obtain the base's twist, $\mathbf{t}_0$, in $\mathcal{R}_{ine}$ operating the following base transformation:

$$ \mathbf{t}_0 = \begin{bmatrix} \omega_0 \\ \dot{\mathbf{r}}_{S_0} \end{bmatrix} = \begin{bmatrix} \mathbf{R}_{S_0} & \mathbf{0}_3 \\ \mathbf{0}_3 & \mathbf{I}_3 \end{bmatrix} \begin{bmatrix} \omega_0^{sat} \\ \dot{\mathbf{r}}_{S_0} \end{bmatrix} = \mathbf{P}_0 \dot{\mathbf{q}}_0 \qquad (13) $$

If not specified, all velocities and pose are expressed in $\mathcal{R}_{ine}$. Defining $\mathbf{t} = [\mathbf{t}_0^T \dots \mathbf{t}_{n_q}^T]^T \in \mathbb{R}^{(6+n_q) \times 1}$ as the generalized twist vector, a velocity transformation between the task-space and the joint-space is defined as (Angeles & Lee, 1988):

$$ \mathbf{t} = \begin{bmatrix} \mathbf{I}_6 & \mathbf{0}_6 & \dots & \mathbf{0}_6 \\ \mathbf{B}_{10} & \mathbf{I}_6 & \dots & \mathbf{0}_6 \\ \vdots & \vdots & \ddots & \vdots \\ \mathbf{B}_{n_q 0} & \mathbf{B}_{n_q 1} & \dots & \mathbf{I}_6 \end{bmatrix} \begin{bmatrix} \mathbf{P}_0 & \mathbf{0}_6 & \dots & \mathbf{0}_6 \\ \mathbf{0}_6 & \mathbf{p}_{m_1} & \dots & \mathbf{0}_6 \\ \vdots & \vdots & \ddots & \vdots \\ \mathbf{0}_6 & \mathbf{0}_6 & \dots & \mathbf{p}_{m_{n_q}} \end{bmatrix} \begin{bmatrix} \dot{\mathbf{q}}_0 \\ \dot{\mathbf{q}} \end{bmatrix} \qquad (14) $$

with:

* $\mathbf{B}_{kj} = \begin{bmatrix} \mathbf{I}_3 \\ (\mathbf{r}_{S_k} - \mathbf{r}_{S_j})^\times \end{bmatrix} \in \mathbb{R}^{6 \times 6}$, if $S_k$ and $S_j$ are in the same kinematic chain

* $\mathbf{B}_{kj} = \mathbf{0}_6$, if $S_k$ and $S_j$ are in different kinematic chains

* $\mathbf{p}_{m_k} = \begin{bmatrix} \mathbf{k}_k \\ \mathbf{k}_k \times (\mathbf{r}_{S_k} - \mathbf{r}_{A_k}) \end{bmatrix} \in \mathbb{R}^{6 \times 1}$, if $A_k$ is a revolute joint

* $\mathbf{p}_{m_k} = \begin{bmatrix} \mathbf{0}_{3 \times 1} \\ \mathbf{k}_k \end{bmatrix} \in \mathbb{R}^{6 \times 1}$, if $A_k$ is a prismatic joint

* $\mathbf{p}_{m_k} = \mathbf{0}_{6 \times 1}$, if $A_k$ is a fixed joint

where $\mathbf{x}^\times$ refers to the skew-symmetric matrix¹ of vector $\mathbf{x} \in \mathbb{R}^{3 \times 1}$, $\mathbf{p}_{m_k}$ is the twist propagation vector, $\mathbf{B}_{ij}$ is the twist propagation matrix and $\mathbf{N}$ is the velocity transformation matrix. Its expression is function of the geometry and configuration of the system. The matrix $\mathbf{N}$ can be decoupled into a lower and an upper triangular block matrix (Virgili-Llop et al., 2016). Moreover,

$$ \mathbf{1} \mathbf{x}^\times = \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix}^\times = \begin{bmatrix} 0 & -x_3 & x_2 \\ x_3 & 0 & -x_1 \\ -x_2 & x_1 & 0 \end{bmatrix} $$

thanks to vectors $\mathbf{p}_{m_k}$ and matrices $\mathbf{B}_{kj}$, the relation (14) enables to consider multiple kinematic chains with a common base. Likewise, with the recursivity of the formalism, from (14) the twist expression can be expressed as:

$$ \mathbf{t}_0 = \mathbf{P}_0 \dot{\mathbf{q}}_0 \qquad (15a) $$

$$ \mathbf{t}_i = \mathbf{B}_{i0} \mathbf{P}_0 \dot{\mathbf{q}}_0 + \sum_{k=1}^{i-1} (\mathbf{B}_{ik} \mathbf{p}_{m_k} \dot{\mathbf{q}}_k) + \mathbf{p}_{m_i} \dot{\mathbf{q}}_i \qquad (15b) $$

In space robotics, and more largely for robots that involve multi-body kinematic chains, the term Jacobian is employed to refer the transformation matrix between a body twist, or the twist of a point in the kinematic chain considered, and the generalized joint rate-variables vector. With the rows of $\mathbf{N}$, the twist of the ith body in the kinematic chain is given as function of the Jacobian matrix, $\mathbf{J}_i \in \mathbb{R}^{6 \times (6+n_q)}$, by:

$$ \mathbf{t}_i = \begin{bmatrix} \mathbf{B}_{i0} \mathbf{P}_0 & \mathbf{B}_{i1} \mathbf{p}_{m_1} & \dots & \mathbf{B}_{i(i-1)} \mathbf{p}_{m_{i-1}} & \mathbf{p}_{m_i} \\ \mathbf{J}_i(\mathbf{q}) & \mathbf{J}_i(\dot{\mathbf{q}}_0) & \dots & \mathbf{J}_i(\dot{\mathbf{q}}_i) & \mathbf{J}_i(\dot{\mathbf{q}}) \end{bmatrix} \begin{bmatrix} \dot{\mathbf{q}}_0 \\ \dot{\mathbf{q}} \end{bmatrix} \qquad (16) $$

Likewise, $\mathbf{J}_i$ can be decomposed into a Jacobian, $\mathbf{J}_{0i} \in \mathbb{R}^{6 \times 6}$, for the base contribution on $\mathbf{t}_i$ and a second Jacobian, $\mathbf{J}_{qi} \in \mathbb{R}^{6 \times n_q}$, for the kinematic chain(s) DoF contributions. This allows to re-write (16), $\forall i \in [1, n_q - 1]$, as:

$$ \mathbf{t}_i = \mathbf{B}_{i0} \mathbf{P}_0 \dot{\mathbf{q}}_0 + \begin{bmatrix} \mathbf{B}_{i1} \mathbf{p}_{m_1} & \dots & \mathbf{B}_{i(i-1)} \mathbf{p}_{m_{i-1}} & \mathbf{p}_{m_i} & \mathbf{0}_{6 \times (n_q - i)} \end{bmatrix} \begin{bmatrix} \dot{\mathbf{q}}_1 \\ \vdots \\ \dot{\mathbf{q}}_{n_q} \end{bmatrix} \qquad (17) $$

For the $n_p$ flexible bodies ending a kinematic chain, the ith flexible body referred as the $P_{fi}^{th}$ one, one can adapt Eq. (17) to express the twist of each flexible bodies in $\mathcal{R}_{ine}$ such that $\forall \mathbf{P}_{fi} \in [0, n_p]$ its expression is given by:

$$ \mathbf{t}_{q_{P_{fi}}} = \begin{bmatrix} \omega_{P_{fi}} \\ \dot{\mathbf{r}}_{P_{fi}} \end{bmatrix} = \mathbf{J}_{0P_{fi}} \dot{\mathbf{q}}_0 + \mathbf{J}_{q_{P_{fi}}} \dot{\mathbf{q}} \qquad (18) $$

Capitalizing on the generic formulation, all accelerations are obtained with the time-derivative of the twist vector (17) such that $\forall i \in [1, n_q]$:

$$ \dot{\mathbf{t}}_i = \dot{\mathbf{j}}_0 \dot{\mathbf{q}}_0 + \mathbf{J}_{0i} \ddot{\mathbf{q}}_0 + \dot{\mathbf{J}}_{qi} \dot{\mathbf{q}} + \mathbf{J}_{qi} \ddot{\mathbf{q}} \qquad (19) $$

where:

$$ \begin{cases} \dot{\mathbf{j}}_{0i} = \begin{bmatrix} \mathbf{I}_3 \\ (\mathbf{r}_{S_0} - \mathbf{r}_{S_i})^\times \end{bmatrix} \mathbf{0}_3 \\ \Omega_0 = \begin{bmatrix} \omega_0^\times & \mathbf{0}_3 \\ \mathbf{0}_3 & \mathbf{0}_3 \end{bmatrix} \end{cases} \qquad (20a) $$

$$ \begin{cases} \dot{\mathbf{B}}_{kj} = \begin{bmatrix} \mathbf{0}_3 \\ (\dot{\mathbf{r}}_{S_k} - \dot{\mathbf{r}}_{S_j})^\times \end{bmatrix} \mathbf{0}_3 \\ \dot{\mathbf{J}}_{qi} = \begin{bmatrix} \dot{\mathbf{J}}_{qi,1} & \dots & \dot{\mathbf{J}}_{qi,k} & \dots & \dot{\mathbf{J}}_{qi,n_q} \end{bmatrix} \end{cases} \qquad (20b) $$

$$ \dot{\mathbf{J}}_{qi,k} = \begin{bmatrix} \mathbf{I}_3 \\ (\mathbf{r}_{S_k} - \mathbf{r}_{S_i})^\times \end{bmatrix} \mathbf{0}_3 \qquad (20c) $$

$$ \dot{\mathbf{J}}_{qi,k} = \begin{bmatrix} \mathbf{I}_3 \\ (\mathbf{r}_{S_k} - \mathbf{r}_{S_i})^\times \end{bmatrix} \mathbf{0}_3 \mathbf{\Omega}_k \mathbf{p}_{m_k} + \dot{\mathbf{B}}_{ki} \mathbf{p}_{m_k} \qquad (20d) $$

$$ \mathbf{\Omega}_k = \begin{bmatrix} \omega_k^\times & \mathbf{0}_3 \\ \mathbf{0}_3 & \omega_k^\times \end{bmatrix} \qquad (20e) $$
## 2.3. Multi-body systems dynamics with flexible bodies

In order to develop the dynamic with a Lagrangian approach, the energy and forces/torques of the overall system are first resumed. The kinetic energy of the complete system, $\mathcal{T}$, can be detailed as:

$$ \mathcal{T} = \mathcal{T}_{hub} + \mathcal{T}_{actuators} + \mathcal{T}_{appendages} \quad (21) $$

where $\mathcal{T}_{hub}$, $\mathcal{T}_{actuators}$ and $\mathcal{T}_{appendages}$ respectively correspond to the hub (spacecraft platform and rigid element connected not part of a kinematic chain), the actuators and the appendages kinetic energies. Additionally, due to the presence of flexible elements a potential energy, $\mathcal{V}_\eta$, and a dissipative force, $\mathbf{F}_\eta$, are considered. Energies of rigid kinematic chains composing the SMS are given as:

$$ \mathcal{T}_{hub} = \frac{1}{2} \mathbf{t}_0^T \begin{bmatrix} \mathcal{I}_0 & \mathbf{0}_3 \\ \mathbf{0}_3 & m_0 \mathbf{I}_3 \end{bmatrix} \mathbf{t}_0 \quad (22a) $$

$$ \mathcal{T}_{actuators} = \frac{1}{2} \sum_{i=1}^{n_q} \left( \mathbf{t}_i^T \begin{bmatrix} \mathcal{I}_i & \mathbf{0}_3 \\ \mathbf{0}_3 & m_i \mathbf{I}_3 \end{bmatrix} \mathbf{t}_i \right) \quad (22b) $$

and with Eqs. (15) and (17) the rigid terms' kinetic energies can be expressed in function of the actuator generalized joint rate-variable, $\dot{\mathbf{q}}$, and the base velocity joint state, $\dot{\mathbf{q}}_0$, as:

$$ \mathcal{T}_{hub} = \frac{1}{2} \dot{\mathbf{q}}_0^T \mathbf{P}_0^T \begin{bmatrix} \mathcal{I}_0 & \mathbf{0}_3 \\ \mathbf{0}_3 & m_0 \mathbf{I}_3 \end{bmatrix} \mathbf{P}_0 \dot{\mathbf{q}}_0 \quad (23a) $$

$$ \begin{aligned} \mathcal{T}_{actuators} = & \frac{1}{2} \dot{\mathbf{q}}_0^T \sum_{i=1}^{n_q} \begin{bmatrix} \mathbf{J}_{0i}^T & \begin{bmatrix} \mathcal{I}_i & \mathbf{0}_3 \\ \mathbf{0}_3 & m_i \mathbf{I}_3 \end{bmatrix} \mathbf{J}_{0i} \end{bmatrix} \dot{\mathbf{q}}_0 \\ & + \frac{1}{2} \dot{\mathbf{q}}_0^T \sum_{i=1}^{n_q} \begin{bmatrix} \mathbf{J}_{qi}^T & \begin{bmatrix} \mathcal{I}_i & \mathbf{0}_3 \\ \mathbf{0}_3 & m_i \mathbf{I}_3 \end{bmatrix} \mathbf{J}_{qi} \end{bmatrix} \dot{\mathbf{q}}_0 \\ & + \frac{1}{2} \dot{\mathbf{q}}_0^T \sum_{i=1}^{n_q} \begin{bmatrix} \mathbf{J}_{qi}^T & \begin{bmatrix} \mathcal{I}_i & \mathbf{0}_3 \\ \mathbf{0}_3 & m_i \mathbf{I}_3 \end{bmatrix} \mathbf{J}_{0i} \end{bmatrix} \dot{\mathbf{q}}_0 \\ & + \frac{1}{2} \dot{\mathbf{q}}_0^T \sum_{i=1}^{n_q} \begin{bmatrix} \mathbf{J}_{0i}^T & \begin{bmatrix} \mathcal{I}_i & \mathbf{0}_3 \\ \mathbf{0}_3 & m_i \mathbf{I}_3 \end{bmatrix} \mathbf{J}_{qi} \end{bmatrix} \dot{\mathbf{q}}_0 \end{aligned} \quad (23b) $$

$$ \begin{aligned} \mathcal{T}_{actuators} = & \frac{1}{2} \dot{\mathbf{q}}_0^T \sum_{i=1}^{n_q} \begin{bmatrix} \mathbf{J}_{0i}^T & \begin{bmatrix} \mathcal{I}_i & \mathbf{0}_3 \\ \mathbf{0}_3 & m_i \mathbf{I}_3 \end{bmatrix} \mathbf{J}_{0i} \end{bmatrix} \dot{\mathbf{q}}_0 \\ & + \frac{1}{2} \dot{\mathbf{q}}_0^T \sum_{i=1}^{n_q} \begin{bmatrix} \mathbf{J}_{qi}^T & \begin{bmatrix} \mathcal{I}_i & \mathbf{0}_3 \\ \mathbf{0}_3 & m_i \mathbf{I}_3 \end{bmatrix} \mathbf{J}_{qi} \end{bmatrix} \dot{\mathbf{q}}_0 \\ & + \frac{1}{2} \dot{\mathbf{q}}_0^T \sum_{i=1}^{n_q} \begin{bmatrix} \mathbf{J}_{qi}^T & \begin{bmatrix} \mathcal{I}_i & \mathbf{0}_3 \\ \mathbf{0}_3 & m_i \mathbf{I}_3 \end{bmatrix} \mathbf{J}_{0i} \end{bmatrix} \dot{\mathbf{q}}_0 \\ & + \frac{1}{2} \dot{\mathbf{q}}_0^T \sum_{i=1}^{n_q} \begin{bmatrix} \mathbf{J}_{0i}^T & \begin{bmatrix} \mathcal{I}_i & \mathbf{0}_3 \\ \mathbf{0}_3 & m_i \mathbf{I}_3 \end{bmatrix} \mathbf{J}_{qi} \end{bmatrix} \dot{\mathbf{q}}_0 \end{aligned} \quad (23b) $$

In order to integrate the flexible dynamics at the end of a rigid kinematic chain, a similar interaction model as developed in Kraïem et al. (2022) is adapted in this study. Each interactions between a flexible body and the rigid spacecraft will occur on the common junction such that for an external force/torque applying or exciting the appendage a linear or angular displacement occurs and reciprocally any excitation lead to a force/torque response on the junction DoFs. Associating the modal displacement vector $\eta_{\mathbf{P}_{fi}} \in \mathbf{R}^{n_{\mathbf{P}_{fi}}} \times 1$ to the $\mathbf{P}_{fi}^{th}$ flexible body composed of the $n_{\mathbf{P}_{fi}}$ flexible DoFs, the second order equation of motions in Sanfedino et al. (2018) is adapted with an appropriate flexible base change and expressed in $\mathcal{R}_{ine}$ as:

$$ \begin{bmatrix} \bar{\mathbf{M}}_r & \mathbf{L}_{\eta \mathbf{P}_{fi}} \mathbf{P}_{fi} \\ \mathbf{L}_{\eta \mathbf{P}_{fi}}^T \mathbf{P}_{fi} & \mathbf{I}_{n_{\eta \mathbf{P}_{fi}}} \end{bmatrix} \begin{bmatrix} \dot{\mathbf{t}}_{fi} \\ \ddot{\mathbf{t}}_{fi} \end{bmatrix} + \begin{bmatrix} \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{c}_{\eta \mathbf{P}_{fi}} \end{bmatrix} \begin{bmatrix} \mathbf{t}_{fi} \\ \dot{\eta}_{\mathbf{P}_{fi}} \end{bmatrix} + \begin{bmatrix} \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{k}_{\eta \mathbf{P}_{fi}} \end{bmatrix} \begin{bmatrix} \star \\ \eta_{\mathbf{P}_{fi}} \end{bmatrix} = \begin{bmatrix} \mathbf{F}_{\mathbf{P}_{fi}} \\ \mathbf{F}_{\eta \mathbf{P}_{fi}} \end{bmatrix} \quad (24) $$

where:

* $\mathbf{P}_{fi}$ is the index designating the considered flexible body

* $\eta_{\mathbf{P}_{fi}}$ is the modal displacement vector

* $\mathbf{t}_{fi}$ is the twist vector of the fixed joint between the appendages and the satellite base

* $\mathbf{L}_{\eta \mathbf{P}_{fi}} \mathbf{P}_{fi}$ is the matrix of the participation factors

* $\bar{\mathbf{M}}_r$ is the condensed mass matrix of the rigid bodies composed of the junction and the appendage

* $\mathbf{c}_{\eta \mathbf{P}_{fi}}$ is the damping matrix assumed diagonal

* $\mathbf{k}_{\eta \mathbf{P}_{fi}}$ is the diagonal generalized stiffness matrix

* $\mathbf{F}_{\eta \mathbf{P}_{fi}}$ is the forces/torques imposed on the internal flexible DoFs

* $\mathbf{F}_{\mathbf{P}_{fi}}$ is the external forces/torques imposed on the appendage junction

Assuming no external forces/torques applying to the system and from Eq. (19), the interaction dynamic (24) can be adapted and generalized for the $\mathbf{P}_{fi}^{th}$ body as:

$$ \bar{\mathbf{M}} = \begin{bmatrix} \ddot{\mathbf{q}}_0 \\ \ddot{\mathbf{q}} \\ \ddot{\eta}_{\mathbf{P}_{fi}} \end{bmatrix} + \begin{bmatrix} \mathbf{0}_{6 \times 1} \\ \mathbf{0}_{n_q \times 1} \\ \mathbf{c}_{\eta \mathbf{P}_{fi}} \dot{\eta}_{\mathbf{P}_{fi}} \end{bmatrix} + \begin{bmatrix} \mathbf{0}_{6 \times 1} \\ \mathbf{0}_{n_q \times 1} \\ \mathbf{k}_{\eta \mathbf{P}_{fi}} \eta_{\mathbf{P}_{fi}} \end{bmatrix} = \begin{bmatrix} \mathbf{0}_{6 \times 1} \\ \mathbf{0}_{n_q \times 1} \\ \mathbf{0}_{n_{\eta \mathbf{P}_{fi}}} \times 1 \end{bmatrix} \quad (25) $$

with:

$$ \left\{ \begin{aligned} \bar{\mathbf{M}}(\mathbf{q}_0, \mathbf{q}) &= \begin{bmatrix} \mathbf{J}_{0P_{fi}}^T \bar{\mathbf{M}}_{P_{fi}} \mathbf{J}_{0P_{fi}} \\ * \\ * \\ \mathbf{J}_{qP_{fi}}^T \bar{\mathbf{M}}_{P_{fi}} \mathbf{J}_{qP_{fi}} \\ * \\ \mathbf{J}_{qP_{fi}}^T \bar{\mathbf{M}}_{P_{fi}} \mathbf{J}_{qP_{fi}} \end{bmatrix} \begin{bmatrix} \mathbf{J}_{0P_{fi}}^T \mathbf{L}_{\eta \mathbf{P}_{fi}}^T \mathbf{P}_{fi} \\ * \\ \mathbf{J}_{qP_{fi}}^T \mathbf{L}_{\eta \mathbf{P}_{fi}}^T \mathbf{P}_{fi} \\ * \\ \mathbf{I}_{n_{\eta \mathbf{P}_{fi}}} \end{bmatrix} \\ \bar{\mathbf{M}}_{P_{fi}} &= \begin{bmatrix} \mathbf{I}_{P_{fi}} & \mathbf{0}_3 \\ \mathbf{0}_3 & m_{P_{fi}} \mathbf{I}_3 \end{bmatrix} \end{aligned} \right. $$

Thus, the kinetic energy of a flexible appendage, or an element in a one-DoF kinematic chain, is given by:

$$ \mathcal{T}_{appendage} = \frac{1}{2} \begin{bmatrix} \dot{\mathbf{q}}_0^T & \dot{\mathbf{q}}^T & \dot{\eta}_{\mathbf{P}_{fi}}^T \end{bmatrix} \bar{\mathbf{M}}(\mathbf{q}_0, \mathbf{q}) \begin{bmatrix} \dot{\mathbf{q}}_0 \\ \dot{\mathbf{q}} \\ \dot{\eta}_{\mathbf{P}_{fi}} \end{bmatrix} \quad (26) $$

the associated dissipative force, and potential energy are defined as:

$$ \mathbf{F}_{\eta \mathbf{P}_{fi}} = -\mathbf{c}_{\eta \mathbf{P}_{fi}} \dot{\eta}_{\mathbf{P}_{fi}} \quad (27a) $$

$$ \mathcal{V}_{\eta \mathbf{P}_{fi}} = \frac{1}{2} \eta_{\mathbf{P}_{fi}}^T \mathbf{k}_{\eta \mathbf{P}_{fi}} \eta_{\mathbf{P}_{fi}} \quad (27b) $$

which can be generalized to $n_p$ appendages by introducing: the generalized flexible modal displacement vector,

$$ \eta_p \in \mathbf{R}^{n_p} = \sum_{i=1}^{n_p} n_{\eta \mathbf{P}_{fi}} \times 1, \quad \text{defined as } \eta_p = \begin{bmatrix} \eta_{\mathbf{P}_{fi}} & \cdots & \eta_{\mathbf{P}_{fnp}} \end{bmatrix}^T; $$

the generalized forces/torques applying on the appendages $\tau_{ext \eta \mathbf{P}_{fi}} = \begin{bmatrix} \tau_{ext \eta \mathbf{P}_{fi}}^T & \cdots & \tau_{ext \eta \mathbf{P}_{fnp}}^T \end{bmatrix}^T$; the matrix of participation factor to the joint, $\mathbf{L}_{\eta \mathbf{P}_{fi}}$, and the matrix of participation factor to the base, $\mathbf{L}_{\eta 0}$, respectively defined as:

$$ \left\{ \begin{aligned} \mathbf{L}_{\eta \mathbf{P}_{fi}} &= \begin{bmatrix} (\mathbf{L}_{\eta \mathbf{P}_{fi}} \mathbf{P}_{fi})^T & \cdots & (\mathbf{L}_{\eta \mathbf{P}_{fnp}} \mathbf{P}_{fnp})^T & \mathbf{J}_{qP_{fi}} \end{bmatrix}^T \\ \mathbf{L}_{\eta 0} &= \begin{bmatrix} (\mathbf{L}_{\eta \mathbf{P}_{fi}} \mathbf{P}_{fi})^T & \cdots & (\mathbf{L}_{\eta \mathbf{P}_{fnp}} \mathbf{P}_{fnp})^T & \mathbf{J}_{0P_{fi}} \end{bmatrix}^T \end{aligned} \right. $$

Thus, for $n_p$ appendages, Eq. (26) is adapted as:

$$ \mathcal{T}_{appendages} = \frac{1}{2} \begin{bmatrix} \dot{\eta}_{\mathbf{P}}^T & \dot{\mathbf{q}}_0^T & \dot{\mathbf{q}}^T & \mathbf{M}(\mathbf{q}_0, \mathbf{q}) \begin{bmatrix} \dot{\eta}_{\mathbf{P}} \\ \dot{\mathbf{q}}_0 \\ \dot{\mathbf{q}} \end{bmatrix} \end{bmatrix} \quad (28) $$
with:

$$
\mathbf{M}(\mathbf{q}_0, \mathbf{q}) =
\begin{bmatrix}
\mathbf{L}_{n_{\eta p}} & \mathbf{L}_{\eta 0} & \mathbf{L}_{\eta p_f} \\
\mathbf{L}_{\eta 0}^T & \sum_{P_{f_i} \in C_f} (\mathbf{J}_{0P_{f_i}}^T \bar{\mathbf{M}}_{P_{f_i}} \mathbf{J}_{0P_{f_i}}) & \sum_{P_{f_i} \in C_f} (\mathbf{J}_{0P_{f_i}}^T \bar{\mathbf{M}}_{P_{f_i}} \mathbf{J}_{qP_{f_i}}) \\
\mathbf{L}_{\eta p_f}^T & \sum_{P_{f_i} \in C_f} (\mathbf{J}_{qP_{f_i}}^T \bar{\mathbf{M}}_{P_{f_i}} \mathbf{J}_{0P_{f_i}}) & \sum_{P_{f_i} \in C_f} (\mathbf{J}_{qP_{f_i}}^T \bar{\mathbf{M}}_{P_{f_i}} \mathbf{J}_{qP_{f_i}})
\end{bmatrix}
\quad (29)
$$

and

$$
\mathbf{F}_{\eta p_f} = -\operatorname{diag}(\mathbf{c}_{\eta p_f}) \dot{\eta}_{\mathbf{P}} = -\mathbf{c}_{\eta p} \dot{\eta}_{\mathbf{P}} \quad (30a)
$$

$$
\begin{equation}
\left\{
\begin{aligned}
\mathbf{V}_{\eta p_f} &= \frac{1}{2} \eta_{\mathbf{P}}^T \operatorname{diag}(\mathbf{k}_{\eta p_f}) \eta_{\mathbf{P}} = \frac{1}{2} \eta_{\mathbf{P}}^T \mathbf{K}_{\eta p} \eta_{\mathbf{P}} && (30b) \\
\mathbf{F}_{\eta p_f} &= -\operatorname{diag}(\mathbf{c}_{\eta p_f}) \dot{\eta}_{\mathbf{P}} = -\mathbf{c}_{\eta p} \dot{\eta}_{\mathbf{P}} && (30c)
\end{aligned}
\right.
\end{equation}
$$

With Eqs. (23) and (28), the total kinetic energy of the system (21) can be expressed as:

$$
\tau = \frac{1}{2} \begin{bmatrix} \dot{\mathbf{q}}_0^T & \dot{\mathbf{q}}^T & \dot{\eta}_P^T \end{bmatrix} \underbrace{\begin{bmatrix} \mathbf{H}_0 & \mathbf{H}_{0q} & \mathbf{H}_{0\eta p} \\ \star & \mathbf{H}_q & \mathbf{H}_{q\eta p} \\ \star & \star & \mathbf{H}_{\eta p} \end{bmatrix}}_{\mathbf{H}(\mathbf{q}_0, \mathbf{q})} \begin{bmatrix} \dot{\mathbf{q}}_0 \\ \dot{\mathbf{q}} \\ \dot{\eta}_P \end{bmatrix} \quad (31)
$$

where the inertia matrix, **H**(q₀, q), is a positive matrix that varies
in function of the robot configuration. Each terms are detailed
with Eqs. (23) and (28) as:

$$
\mathbf{H}_0 = \mathbf{P}_0^T \begin{bmatrix} \mathcal{I}_0 & \mathbf{0}_3 \\ \mathbf{0}_3 & m_0 I_3 \end{bmatrix} \mathbf{P}_0 + \sum_{i=1}^{n_q} (\mathbf{J}_{0i}^T \bar{\mathbf{M}}_i \mathbf{J}_{0i}) \quad (32a)
$$

$$
\mathbf{H}_{0q} = \sum_{i=1}^{n_q} (\mathbf{J}_{0i}^T \bar{\mathbf{M}}_i \mathbf{J}_{qi}) \quad (32b)
$$

$$
\mathbf{H}_{0\eta p} = \mathbf{L}_{\eta 0} \quad (32c)
$$

$$
\mathbf{H}_q = \sum_{i=1}^{n_q} (\mathbf{J}_{q_i}^T \bar{\mathbf{M}}_i \mathbf{J}_{q_i}) \quad (32d)
$$

$$
\mathbf{H}_{q\eta p} = \mathbf{L}_{\eta p_f} \quad (32e)
$$

$$
\mathbf{H}_{\eta p} = \mathbf{I}_{n_{\eta p}} \qquad (32f)
$$

aßhe system's total energy is expressed with Eqs. (30) and (31)

$$
\mathcal{L} = \tau - V_{\eta p} \tag{33}
$$

System dynamics are then obtained evaluating the Euler–Lagrange
equations of (33) for each DoFs such as:

$$
\left\{
\begin{array}{l}
\displaystyle \frac{d}{dt} \left( \frac{\partial \mathcal{L}}{\partial \dot{\mathbf{q}}_0} \right) - \frac{\partial \mathcal{L}}{\partial \mathbf{q}_0} = \boldsymbol{\tau}_0 + \boldsymbol{\tau}_{\text{ext}_0} \hfill (34\text{a}) \\[1.5ex]
\displaystyle \frac{d}{dt} \left( \frac{\partial \mathcal{L}}{\partial \dot{\mathbf{q}}} \right) - \frac{\partial \mathcal{L}}{\partial \mathbf{q}} = \boldsymbol{\tau}_q + \boldsymbol{\tau}_{\text{ext}_q} \hfill (34\text{b})
\end{array}
\right.
$$

$$
\left\{
\begin{array}{l}
\displaystyle \frac{d}{dt} \left( \frac{\partial \mathcal{L}}{\partial \dot{\mathbf{q}}} \right) - \frac{\partial \mathcal{L}}{\partial \mathbf{q}} = \boldsymbol{\tau}_q + \boldsymbol{\tau}_{\text{ext}_q} \hfill (34\text{b}) \\[1.5ex]
\displaystyle \frac{d}{dt} \left( \frac{\partial \mathcal{L}}{\partial \dot{\boldsymbol{\eta}}_{\mathbf{P}}} \right) - \frac{\partial \mathcal{L}}{\partial \boldsymbol{\eta}_{\mathbf{P}}} = \boldsymbol{\tau}_{\eta_{\mathbf{P}}} + \boldsymbol{\tau}_{\text{ext}_{\eta_{\mathbf{P}}}} = \boldsymbol{\tau}_{\text{ext}_{\eta_{\mathbf{P}}}} \hfill (34\text{c})
\end{array}
\right.
$$

$$
\left\{
\begin{array}{l}
\displaystyle \frac{d}{dt} \left( \frac{\partial \mathcal{L}}{\partial \dot{\eta}_{\mathbf{P}}} \right) - \frac{\partial \mathcal{L}}{\partial \eta_{\mathbf{P}}} = \boldsymbol{\tau}_{\eta_{\mathbf{P}}} + \boldsymbol{\tau}_{\text{ext}_{\eta_{\mathbf{P}}}} = \boldsymbol{\tau}_{\text{ext}_{\eta_{\mathbf{P}}}} \quad (34c) \\[1.5ex]
\displaystyle \frac{d}{dt} \left( \frac{\partial \mathcal{L}}{\partial \dot{\mathbf{q}}} \right) - \frac{\partial \mathcal{L}}{\partial \mathbf{q}} = \boldsymbol{\tau}_q + \boldsymbol{\tau}_{\text{ext}_q} \quad (34b)
\end{array}
\right.
$$

where $\boldsymbol{\tau}_0 \in \mathbb{R}^{6 \times 1}$, $\boldsymbol{\tau}_q \in \mathbb{R}^{n_q \times 1}$ and $\boldsymbol{\tau}_{\eta_{\mathbf{P}}} \in \mathbb{R}^{n_{\eta p} \times 1}$ are respectively the base, the joints and flexible mode control torques/forces. In this study, no means are considered to actively control the flexible mode displacement, thus $\boldsymbol{\tau}_{\eta_{\mathbf{P}}} = \mathbf{0}_{n_{\eta p}} \times 1$. $\boldsymbol{\tau}_{\text{ext}_0} \in \mathbb{R}^{6 \times 1}$, $\boldsymbol{\tau}_{\text{ext}_q} \in \mathbb{R}^{n_q \times 1}$ and $\boldsymbol{\tau}_{\text{ext}_{\eta_{\mathbf{P}}}} \in \mathbb{R}^{n_{\eta p} \times 1}$ are respectively external torques/forces applying on the system's base, joints and flexible mode DoFs. Developing expressions (34) allows to express the

system dynamics under a compact form as:

$$
\left[
\begin{array}{c}
\mathbf{H}_0 \quad \mathbf{H}_{0q} \quad \mathbf{H}_{0\eta p} \\
\star \quad \mathbf{H}_q \quad \mathbf{H}_{q\eta p} \\
\star \quad \star \quad \mathbf{H}_{\eta p}
\end{array}
\right]
\left[
\begin{array}{c}
\ddot{\mathbf{q}}_0 \\
\ddot{\mathbf{q}} \\
\ddot{\eta}_{\mathbf{P}}
\end{array}
\right]
+
\underbrace{
\left[
\begin{array}{c}
\mathbf{C}_0 \quad \mathbf{C}_{0q} \quad \mathbf{C}_{0\eta p} \\
\mathbf{C}_{q0} \quad \mathbf{C}_q \quad \mathbf{C}_{q\eta p} \\
\mathbf{C}_{\eta p0} \quad \mathbf{C}_{\eta pq} \quad \mathbf{C}_{\eta p}
\end{array}
\right]
}_{\mathbf{C}(\mathbf{q}_0, \mathbf{q}, \dot{\mathbf{q}}, \dot{\mathbf{q}}, \dot{\eta}_{\mathbf{P}})}
\left[
\begin{array}{c}
\dot{\mathbf{q}}_0 \\
\dot{\mathbf{q}} \\
\dot{\eta}_{\mathbf{P}}
\end{array}
\right]
+
\left[
\begin{array}{c}
\mathbf{0}_{6 \times 1} \\
\mathbf{0}_{n_q \times 1} \\
\mathbf{K}_{\eta p} \eta_{\mathbf{P}}
\end{array}
\right]
=
\left[
\begin{array}{c}
\boldsymbol{\tau}_0 \\
\boldsymbol{\tau}_q \\
\mathbf{0}_{n_{\eta p} \times 1}
\end{array}
\right]
+
\left[
\begin{array}{c}
\boldsymbol{\tau}_{\text{ext}_0} \\
\boldsymbol{\tau}_{\text{ext}_q} \\
\boldsymbol{\tau}_{\text{ext}_{\eta p}}
\end{array}
\right]
\tag{35}
$$

with the convective matrix, **C**(q₀, q, ˙q₀, ˙q), corresponding to the Coriolis and Centrifugal terms. Matrix **c** is detailed as:

$$
\mathbf{c}_0 = \frac{\partial}{\partial \mathbf{q}_0} (\dot{\mathbf{q}}_0^T \mathbf{H}_0 + \dot{\mathbf{q}}^T \mathbf{H}_{0q}^T + \dot{\eta}_P^T \mathbf{H}_{0\eta p}^T) \quad (36a)
$$

$$
\mathbf{c}_{0q} = \frac{\partial}{\partial \mathbf{q}_0} (\dot{\mathbf{q}}_0^T \mathbf{H}_{0q} + \dot{\mathbf{q}}^T \mathbf{H}_q + \dot{\eta}_P^T \mathbf{H}_{q\eta p}^T) \quad (36b)
$$

$$
\mathbf{c}_{0\eta p} = \frac{\partial}{\partial \mathbf{q}_0} (\dot{\mathbf{q}}_0^T \mathbf{H}_{0\eta p} + \dot{\mathbf{q}}^T \mathbf{H}_{q\eta p} + \dot{\eta}_P^T \mathbf{H}_{\eta p}^T) \quad (36c)
$$

$$
\mathbf{c}_{q0} = \frac{\partial}{\partial \mathbf{q}} (\dot{\mathbf{q}}_0^T \mathbf{H}_0 + \dot{\mathbf{q}}^T \mathbf{H}_{0q}^T + \dot{\eta}_P^T \mathbf{H}_{0\eta p}^T) \quad (36d)
$$

$$
\mathbf{c}_q = \frac{\partial}{\partial \mathbf{q}} (\dot{\mathbf{q}}_0^T \mathbf{H}_{0q} + \dot{\mathbf{q}}^T \mathbf{H}_q + \dot{\eta}_P^T \mathbf{H}_{q\eta p}^T) \quad (36e)
$$

$$
\mathbf{c}_{q\eta p} = \frac{\partial}{\partial \mathbf{q}} (\dot{\mathbf{q}}_0^T \mathbf{H}_{0\eta p} + \dot{\mathbf{q}}^T \mathbf{H}_{q\eta p} + \dot{\eta}_P^T \mathbf{H}_{\eta p}^T) \quad (36f)
$$

$$
\mathbf{c}_{\eta p0} = \frac{\partial}{\partial \eta_{\mathbf{P}}} (\dot{\mathbf{q}}_0^T \mathbf{H}_0 + \dot{\mathbf{q}}^T \mathbf{H}_{0q}^T + \dot{\eta}_P^T \mathbf{H}_{0\eta p}^T) \quad (36g)
$$

$$
\mathbf{c}_{\eta p q} = \frac{\partial}{\partial \eta_{\mathbf{P}}} (\dot{\mathbf{q}}_0^T \mathbf{H}_{0q} + \dot{\mathbf{q}}^T \mathbf{H}_q + \dot{\eta}_P^T \mathbf{H}_{q\eta p}^T) \quad (36h)
$$

$$
\mathbf{c}_{\eta p} = \frac{\partial}{\partial \eta_{\mathbf{P}}} (\dot{\mathbf{q}}_0^T \mathbf{H}_{0\eta p} + \dot{\mathbf{q}}^T \mathbf{H}_{q\eta p} + \dot{\eta}_P^T \mathbf{H}_{\eta p}^T) \quad (36i)
$$

with the partial derivative expressions given in Appendix. The
time-derivative expressions of the inertia sub-matrices are ob-
tained with Eqs. (20) as:

$$
\dot{\mathbf{H}}_0 = \begin{bmatrix} \dot{\mathcal{I}}_0 & \mathbf{0}_3 \\ \mathbf{0}_3 & \mathbf{0}_3 \end{bmatrix} + \sum_{i=1}^{n_q} (\mathbf{J}_{0i}^T \bar{\mathbf{M}}_i \mathbf{J}_{0i}) \quad (37a)
$$

$$
+ \sum_{i=1}^{n_q} (\mathbf{J}_{0i}^T \dot{\bar{\mathbf{M}}}_i \mathbf{J}_{0i}) + \sum_{i=1}^{n_q} (\mathbf{J}_{0i}^T \bar{\mathbf{M}}_i \mathbf{J}_{0i}) \quad (37b)
$$

$$
\dot{\mathbf{H}}_{0q} = \sum_{i=1}^{n_q} (\dot{\mathbf{J}}_{0i}^T \bar{\mathbf{M}}_i \mathbf{J}_{qi}) + \sum_{i=1}^{n_q} (\mathbf{J}_{0i}^T \dot{\bar{\mathbf{M}}}_i \mathbf{J}_{qi}) + \sum_{i=1}^{n_q} (\mathbf{J}_{0i}^T \bar{\mathbf{M}}_i \dot{\mathbf{J}}_{qi}) \quad (37\text{b})
$$

$$
\dot{\mathbf{H}}_{0\eta p} = \begin{bmatrix} (\mathbf{L}_{\eta p f_1} \mathbf{P}_{f_1} \dot{\mathbf{J}}_{0 p_{f_1}}) \\ \vdots \\ (\mathbf{L}_{\eta p f_{np}} \mathbf{P}_{f_{np}} \dot{\mathbf{J}}_{0 p_{f_{np}}}) \end{bmatrix} \quad (37c)
$$

$$
\dot{\mathbf{H}}_q = \sum_{i=1}^{n_q} (\mathbf{J}_{qi}^T \bar{\mathbf{M}}_i \mathbf{J}_{qi}) + \sum_{i=1}^{n_q} (\mathbf{J}_{qi}^T \dot{\bar{\mathbf{M}}}_i \mathbf{J}_{qi}) + \sum_{i=1}^{n_q} (\mathbf{J}_{qi}^T \bar{\mathbf{M}}_i \dot{\mathbf{J}}_{qi}) \quad (37\text{d})
$$

$$
\dot{\mathbf{H}}_{q\eta p} = \begin{bmatrix} (\mathbf{L}_{\eta p f_1} \mathbf{P}_{f_1} \dot{\mathbf{J}}_{q p_{f_1}}) \\ \vdots \\ (\mathbf{L}_{\eta p f_{np}} \mathbf{P}_{f_{np}} \dot{\mathbf{J}}_{q p_{f_{np}}}) \end{bmatrix} \quad (37e)
$$

$$
\dot{\mathbf{H}}_{\eta p} = \mathbf{0}_{n_{\eta p}} \quad (37f)
$$
<figure><img src="image_3.png" /><figcaption>Fig. 2. Block diagram of the proposed control.</figcaption></figure>

with the time-derivative of the inertia matrices $I_i$ (9) in $R_{ine}$ given, $\forall i \in [0, n_q]$, by:

$$ \dot{\mathbf{I}}_i = \omega_i^\times \mathcal{R}_{S_i, \mathcal{R}_{ine}} \mathcal{I}_{S_i} \mathcal{R}_{S_i, \mathcal{R}_{ine}}^T + \mathcal{R}_{S_i, \mathcal{R}_{ine}} \mathcal{I}_{S_i} \mathcal{R}_{S_i, \mathcal{R}_{ine}}^T \omega_i^\times \quad (38) $$

Flexible modes are obtained from a linearization process and expressed in a normalized base which leads to the independence of inertia matrices to the flexible mode displacement. Thus, the following simplifications in the convective terms (36) can be made:

$$ \mathbf{c}_{\eta_p 0} = \mathbf{0}_{n_{\eta_p} \times 6} \quad (39a) $$

$$ \mathbf{c}_{\eta_p q} = \mathbf{0}_{n_{\eta_p} \times n_q} \quad (39b) $$

$$ \mathbf{c}_{\eta_p} = \mathbf{0}_{n_{\eta_p}} \quad (39c) $$

$$ \frac{\partial}{\partial \mathbf{q}_0} (\dot{\eta}_P^T \mathbf{H}_{\eta_p}) = \mathbf{0}_{6 \times n_{\eta_p}} \quad (39d) $$

$$ \frac{\partial}{\partial \mathbf{q}} (\dot{\eta}_P^T \mathbf{H}_{\eta_p}) = \mathbf{0}_{n_q \times n_{\eta_p}} \quad (39e) $$

With the hypothesis of a free-floating SMS, the equation of motions (35) not subject to any external torques/forces are reduced to:

$$ \mathbf{H}(\mathbf{q}, \mathbf{q}_0) \begin{bmatrix} \ddot{\mathbf{q}}_0 \\ \ddot{\mathbf{q}} \\ \ddot{\eta}_P \end{bmatrix} + \mathbf{C}(\mathbf{q}_0, \mathbf{q}, \dot{\mathbf{q}}_0, \dot{\mathbf{q}}, \dot{\eta}_P) \begin{bmatrix} \dot{\mathbf{q}}_0 \\ \dot{\mathbf{q}} \\ \dot{\eta}_P \end{bmatrix} \\ + \begin{bmatrix} \mathbf{0}_{6 \times 1} \\ \mathbf{0}_{n_q \times 1} \\ \mathbf{K}_{\eta_p} \mathbf{\eta}_P \end{bmatrix} = \begin{bmatrix} \mathbf{0}_{6 \times 1} \\ \tau_q \\ \mathbf{0}_{n_{\eta_p}} \end{bmatrix} \quad (40) $$

### 3. Control strategy

After deriving the equation of motions of a free-floating SMS with flexible appendages ending kinematic chains, a robust spacecraft's attitude and manipulator(s) control framework is proposed to perform OOS of a target with flexible appendages. The control strategy consists in an NDI to linearize and decouple the system dynamics, as illustrated in Fig. 2. To provide robustness towards modeling uncertainties, sensor noises and disturbances due to the vibrations of flexible appendages, an $H_\infty$ synthesis is presented. When implementing an NDI, the derivation of the inverse dynamics/kinematics presents the drawback of being highly time-consuming. In order to validate the proposed control strategy on a real-time simulation platform (Kraïem et al., 2024), the evaluation of the inverse dynamics/kinematics is provided by the open-source spacecraft robotics toolkit (Virgili-Llop et al.,

2016) extended as developed in (Kraïem et al., 2022). This advantageously provides a fast enough processing time to evaluate control torques at a 100 Hz rate. However, this toolbox implemented with a Newton-Euler approach neglects the nonlinear couplings provided by the differential derivatives terms (36). Hence, the $H_\infty$ feedback is designed such that the modeling errors, necessary to perform real-time control, are not impacting the performances. The control strategy is applied to a rotation-free-floating SMS with reaction wheels as such actuators are commonly used for attitude control (Cumer et al., 2021; Wilde et al., 2018) which both allows to simplify equations of motion and focus on the flexible disturbance rejection. Moreover, it is assumed that the flexible displacements are not measured.

#### 3.1. Open-loop dynamics

When considering only reaction wheels as base actuators, the equation of motion (40) can be simplified and detailed as:

$$ \begin{bmatrix} \mathbf{H}_0^\omega & \mathbf{H}_0^{\omega L} & \mathbf{H}_{0qr}^\omega & \mathbf{H}_{0qm}^\omega & \mathbf{H}_{0\eta_p}^\omega \\ \mathbf{\star} & \mathbf{H}_0^L & \mathbf{0} & \mathbf{H}_{0qm}^L & \mathbf{H}_{0\eta_p}^L \\ \mathbf{\star} & \mathbf{\star} & \mathbf{H}_{qr} & \mathbf{0} & \mathbf{0} \\ \mathbf{\star} & \mathbf{\star} & \mathbf{\star} & \mathbf{H}_{qm} & \mathbf{H}_{qm\eta_p} \\ \mathbf{\star} & \mathbf{\star} & \mathbf{\star} & \mathbf{\star} & \mathbf{H}_{\eta_p} \end{bmatrix} \begin{bmatrix} \dot{\omega}_0 \\ \ddot{\mathbf{r}}_{S0} \\ \dot{\mathbf{q}}_r \\ \dot{\mathbf{q}}_m \\ \ddot{\eta}_P \end{bmatrix} + \begin{bmatrix} \mathbf{C}_0^\omega & \mathbf{C}_0^{\omega L} & \mathbf{C}_{0qr}^\omega & \mathbf{C}_{0qm}^\omega & \mathbf{C}_{0\eta_p}^\omega \\ \mathbf{C}_0^{\omega L} & \mathbf{C}_0^L & \mathbf{0} & \mathbf{C}_{0qm}^L & \mathbf{C}_{0\eta_p}^L \\ \mathbf{C}_{q_r 0}^\omega & \mathbf{0} & \mathbf{C}_{qr} & \mathbf{0} & \mathbf{0} \\ \mathbf{C}_{q_m 0}^\omega & \mathbf{C}_{q_m 0}^L & \mathbf{0} & \mathbf{C}_{qm} & \mathbf{C}_{qm\eta_p} \\ \mathbf{C}_{\eta_p 0}^\omega & \mathbf{C}_{\eta_p 0}^L & \mathbf{0} & \mathbf{C}_{\eta_p q_m} & \mathbf{C}_{\eta_p} \end{bmatrix} \begin{bmatrix} \dot{\omega}_0 \\ \ddot{\mathbf{r}}_{S0} \\ \dot{\mathbf{q}}_r \\ \dot{\mathbf{q}}_m \\ \ddot{\eta}_P \end{bmatrix} + \begin{bmatrix} \mathbf{0}_{3 \times 1} \\ \mathbf{0}_{3 \times 1} \\ \mathbf{0}_{n_{qr} \times 1} \\ \mathbf{0}_{n_{qm} \times 1} \\ \mathbf{K}_{\eta_p} \mathbf{\eta}_P \end{bmatrix} = \begin{bmatrix} \mathbf{0}_{3 \times 1} \\ \mathbf{0}_{3 \times 1} \\ \tau_{qr} \\ \tau_{qm} \\ \mathbf{0}_{n_{\eta_p} \times 1} \end{bmatrix} \quad (41) $$

in which the following notations were introduced: subscript $m$ designates manipulators quantities, subscript $r$ designates reaction wheels quantities, superscript $\omega$ designates the base angular contribution of a quantity and superscript $L$ the linear contribution. This expression highlights that no direct couplings between the base actuators and the manipulator's joints are considered. Then the open-loop equation of motions describing the angular motions of the spacecraft base and the manipulator joints are obtained by a re-writing of Eq. (41) in which the linear acceleration of the base, $\ddot{\mathbf{r}}_0$, and the reaction wheels accelerations are
not present. First step is to isolate from (41) the angular base and manipulator joints accelerations:

$$
\begin{aligned}
& \begin{bmatrix} \mathbf{H}_0^\omega & \mathbf{H}_{0q_m}^\omega \\ \star & \mathbf{H}_{q_m}^\omega \end{bmatrix} \begin{bmatrix} \dot{\mathbf{\omega}}_0 \\ \ddot{\mathbf{q}}_m \end{bmatrix} + \begin{bmatrix} \mathbf{C}_0^\omega & \mathbf{C}_{0q_m}^\omega \\ \mathbf{C}_{q_m0}^\omega & \mathbf{C}_{q_m}^\omega \end{bmatrix} \begin{bmatrix} \mathbf{\omega}_0 \\ \dot{\mathbf{q}}_m \end{bmatrix} \\
& + \begin{bmatrix} \mathbf{H}_0^{\omega L} & \mathbf{H}_{0qr}^\omega \\ \star & \mathbf{H}_{qr}^\omega \end{bmatrix} \begin{bmatrix} \ddot{\mathbf{r}}_{S0} \\ \ddot{\mathbf{q}}_r \end{bmatrix} + \begin{bmatrix} \mathbf{C}_0^{\omega L} & \mathbf{C}_{0qr}^{\omega} \\ \mathbf{C}_{qr}^{\omega} & \mathbf{C}_{q_r}^{\omega} \end{bmatrix} \begin{bmatrix} \dot{\mathbf{r}}_0 \\ \dot{\mathbf{q}}_r \end{bmatrix} \\
& + \begin{bmatrix} \mathbf{H}_{0\eta_P}^\omega \\ \mathbf{H}_{q_m\eta_P}^\omega \end{bmatrix} \ddot{\mathbf{\eta}}_P + \begin{bmatrix} \mathbf{0} & \mathbf{C}_{0\eta_P}^\omega \\ \mathbf{0} & \mathbf{C}_{q_m\eta_P}^\omega \end{bmatrix} \begin{bmatrix} \dot{\eta}_P \\ \dot{\eta}_P \end{bmatrix} = \begin{bmatrix} \mathbf{0} \\ \tau_{q_m} \end{bmatrix}
\end{aligned}
\qquad (42)
$$

the linear base and reaction wheels accelerations:

$$
\begin{aligned}
& \begin{bmatrix} \mathbf{H}_0^L & \mathbf{0} \\ \mathbf{0} & \mathbf{H}_{qr} \end{bmatrix} \begin{bmatrix} \ddot{\mathbf{r}}_0 \\ \ddot{\mathbf{q}}_r \end{bmatrix} + \begin{bmatrix} \mathbf{C}_0^L & \mathbf{0} \\ \mathbf{0} & \mathbf{C}_{qr} \end{bmatrix} \begin{bmatrix} \dot{\mathbf{r}}_{S0} \\ \dot{\mathbf{q}}_r \end{bmatrix} \\
& + \begin{bmatrix} \mathbf{H}_0^{\omega L^T} & \mathbf{H}_{0q_m}^L \\ \mathbf{H}_{0qr}^{\omega^T} & \mathbf{0} \end{bmatrix} \begin{bmatrix} \dot{\mathbf{\omega}}_0 \\ \ddot{\mathbf{q}}_m \end{bmatrix} + \begin{bmatrix} \mathbf{C}_0^{\omega L} & \mathbf{C}_{0q_m}^L \\ \mathbf{C}_{q_m0}^{\omega} & \mathbf{0} \end{bmatrix} \begin{bmatrix} \mathbf{\omega}_0 \\ \dot{\mathbf{q}}_m \end{bmatrix} \\
& + \begin{bmatrix} \mathbf{H}_{0\eta_P}^L \\ \mathbf{0} \end{bmatrix} \ddot{\mathbf{\eta}}_P + \begin{bmatrix} \mathbf{0} & \mathbf{C}_{0\eta_P}^L \\ \mathbf{0} & \mathbf{0} \end{bmatrix} \begin{bmatrix} \dot{\eta}_P \\ \dot{\eta}_P \end{bmatrix} = \begin{bmatrix} \mathbf{0} \\ \tau_{qr} \end{bmatrix}
\end{aligned}
\qquad (43)
$$

and the flexible accelerations:

$$
\begin{aligned}
& \begin{bmatrix} \mathbf{H}_{0\eta_P}^L & \mathbf{0} \end{bmatrix} \begin{bmatrix} \ddot{\mathbf{r}}_0 \\ \ddot{\mathbf{q}}_r \end{bmatrix} + \begin{bmatrix} \mathbf{C}_{\eta_P 0}^L & \mathbf{0} \end{bmatrix} \begin{bmatrix} \dot{\mathbf{r}}_{S0} \\ \dot{\mathbf{q}}_r \end{bmatrix} \\
& + \begin{bmatrix} \mathbf{H}_{0\eta_P}^{\omega L^T} & \mathbf{H}_{q_m\eta_P}^T \end{bmatrix} \begin{bmatrix} \dot{\mathbf{\omega}}_0 \\ \ddot{\mathbf{q}}_m \end{bmatrix} + \begin{bmatrix} \mathbf{C}_{\eta_P 0}^{\omega} & \mathbf{C}_{q_m\eta_P}^L \end{bmatrix} \begin{bmatrix} \mathbf{\omega}_0 \\ \dot{\mathbf{q}}_m \end{bmatrix} \\
& + \mathbf{H}_{\eta_P} \ddot{\mathbf{\eta}}_P + \begin{bmatrix} \mathbf{K}_{\eta_P} & \mathbf{C}_{\eta_P} \end{bmatrix} \begin{bmatrix} \dot{\eta}_P \\ \dot{\eta}_P \end{bmatrix} = \mathbf{0}
\end{aligned}
\qquad (44)
$$

Combining the previous Eqs. (42), (43) and (44), the open-loop system can be expressed as:

$$
\begin{aligned}
& \mathbf{M}^{\circ}(\mathbf{q}, \mathbf{q}_0) \begin{bmatrix} \dot{\mathbf{\omega}}_0 \\ \ddot{\mathbf{q}}_m \end{bmatrix} + \mathbf{D}^{\circ}(\mathbf{q}, \mathbf{q}_0, \dot{\mathbf{q}}, \dot{\mathbf{q}}_0, \dot{\mathbf{\eta}}_P) \begin{bmatrix} \mathbf{\omega}_0 \\ \dot{\mathbf{q}}_m \end{bmatrix} \\
& + \mathbf{D}^{\circ}_{\mathrm{Lr}}(\mathbf{q}, \mathbf{q}_0, \dot{\mathbf{q}}, \dot{\mathbf{q}}_0, \dot{\mathbf{\eta}}_P) \begin{bmatrix} \dot{\mathbf{r}}_{S0} \\ \dot{\mathbf{q}}_r \end{bmatrix} + \mathbf{M}^{\circ}_{\eta}(\mathbf{q}, \mathbf{q}_0) \ddot{\mathbf{\eta}}_P \\
& + \mathbf{D}^{\circ}_{\eta}(\mathbf{q}, \mathbf{q}_0, \dot{\mathbf{q}}, \dot{\mathbf{q}}_0, \dot{\mathbf{\eta}}_P) \dot{\mathbf{\eta}}_P = \mathbf{J}_{\tau}^{\circ} \tau_q
\end{aligned}
\qquad (45)
$$

with:

$$
\left\{
\begin{aligned}
\mathbf{H}^{\circ} &= \begin{bmatrix} \mathbf{H}_{0}^{\omega L} & \mathbf{H}_{0q_r}^{\omega} \\ \mathbf{H}_{0q_m}^{L^T} & \mathbf{0} \end{bmatrix} \begin{bmatrix} \mathbf{H}_{0}^{L} & \mathbf{0} \\ \mathbf{0} & \mathbf{H}_{qr} \end{bmatrix}^{-1} \\
&= \begin{bmatrix} \mathbf{H}_{0}^{\omega L} \mathbf{H}_{0}^{L^{-1}} & \mathbf{H}_{0qr}^{\omega} \mathbf{H}_{qr}^{-1} \\ \mathbf{H}_{0q_m}^{L^T} \mathbf{H}_{0}^{L^{-1}} & \mathbf{0} \end{bmatrix} \\
\mathbf{M}^{\circ} &= \begin{bmatrix} \mathbf{H}_{0}^{\omega} & \mathbf{H}_{0q_m}^{\omega} \\ \mathbf{H}_{0q_m}^{\omega^T} & \mathbf{H}_{q_m} \end{bmatrix} - \mathbf{H}^{\circ} \begin{bmatrix} \mathbf{H}_{0}^{\omega L^T} & \mathbf{H}_{0q_m}^{L} \\ \mathbf{H}_{0qr}^{\omega^T} & \mathbf{0} \end{bmatrix} \\
\mathbf{D}^{\circ} &= \begin{bmatrix} \mathbf{C}_{0}^{\omega} & \mathbf{C}_{q_m0}^{\omega} \\ \mathbf{C}_{q_m0}^{\omega} & \mathbf{C}_{q_m} \end{bmatrix} - \mathbf{H}^{\circ} \begin{bmatrix} \mathbf{C}_{0}^{\omega L} & \mathbf{C}_{0q_m}^{L} \\ \mathbf{C}_{q_r0}^{\omega} & \mathbf{0} \end{bmatrix} \\
\mathbf{D}^{\circ}_{\mathrm{Lr}} &= \begin{bmatrix} \mathbf{C}_{0}^{\omega L} & \mathbf{C}_{0q_r}^{\omega} \\ \mathbf{C}_{q_m0}^{L} & \mathbf{0} \end{bmatrix} - \mathbf{H}^{\circ} \begin{bmatrix} \mathbf{C}_{0}^{L} & \mathbf{0} \\ \mathbf{0} & \mathbf{C}_{qr} \end{bmatrix} \\
\mathbf{M}^{\circ}_{\eta} &= \begin{bmatrix} \mathbf{H}_{0\eta_P}^{\omega} & \mathbf{H}_{q_m\eta_P}^{\omega} \end{bmatrix} - \mathbf{H}^{\circ} \begin{bmatrix} \mathbf{C}_{0\eta_P}^{L} \\ \mathbf{0} \end{bmatrix} \\
\mathbf{D}^{\circ}_{\eta} &= \begin{bmatrix} \mathbf{C}_{q_m\eta_P}^{\omega} & -\mathbf{H}_{0q_m}^{L^T} \mathbf{H}_{0}^{L^{-1}} \mathbf{C}_{0\eta_P}^{L} \\ \mathbf{C}_{q_m\eta_P}^{\omega} & -\mathbf{H}_{0q_m}^{L} \mathbf{H}_{0}^{L^{-1}} \mathbf{C}_{0\eta_P}^{L} \end{bmatrix} \\
\mathbf{J}_{\tau}^{\circ} &= \begin{bmatrix} -\mathbf{H}^{\circ} & \begin{bmatrix} \mathbf{0}_{3 \times n_r} \\ \mathbf{I}_{n_r} \end{bmatrix} \\ \begin{bmatrix} \mathbf{0}_{3 \times n_m} \\ \mathbf{I}_{n_m} \end{bmatrix} \end{bmatrix} \\
&= \begin{bmatrix} -\mathbf{H}_{0qr}^{\omega} \mathbf{H}_{qr}^{-1} & \mathbf{0}_{3 \times n_m} \\ \mathbf{0}_{n_m \times n_r} & \mathbf{I}_{n_m} \end{bmatrix}
\end{aligned}
\right.
$$

where $\mathbf{M}^{\circ}$ and $\mathbf{M}_{\eta}^{\circ}$ are inertia matrices including the mass participation of non-actuated elements that only depend on the system configuration; $\mathbf{D}^{\circ}$, $\mathbf{D}_{\mathrm{Lr}}^{\circ}$ and $\mathbf{D}_{\eta_P}^{\circ}$ are convective matrices that depend on both the system configuration and velocities, i.e. $(\mathbf{q}_0, \mathbf{q}, \dot{\mathbf{q}}_0, \dot{\mathbf{q}}, \dot{\mathbf{\eta}}_P)$; $\mathbf{J}_{\tau}^{\circ}$ is an equivalent Jacobian matrix that does not depend on the SMS configuration but only on the geometry of the reaction wheels. This expression highlights the internal disturbance torques which impact the actuators control.

A similar re-writing of (41) can be done such that the flexible dynamics can be expressed as a function of the base angular motions and manipulator joints. This provides a distinct perturbation dynamics impacting the control performances. Thus, by injecting (43) in (44), the disturbances due to flexible vibrations are expressed as:

$$
\begin{aligned}
& \mathbf{M}^{\circ}(\mathbf{q}, \mathbf{q}_0) \begin{bmatrix} \dot{\mathbf{\omega}}_0 \\ \ddot{\mathbf{q}}_m \end{bmatrix} + \mathbf{D}^{\circ}(\mathbf{q}, \mathbf{q}_0, \dot{\mathbf{q}}, \dot{\mathbf{q}}_0, \dot{\mathbf{\eta}}_P) \begin{bmatrix} \mathbf{\omega}_0 \\ \dot{\mathbf{q}}_m \end{bmatrix} \\
& + \mathbf{D}_{\mathrm{Lr}}^{\circ}(\mathbf{q}, \mathbf{q}_0, \dot{\mathbf{q}}, \dot{\mathbf{q}}_0, \dot{\mathbf{\eta}}_P) \begin{bmatrix} \dot{\mathbf{r}}_{S0} \\ \dot{\mathbf{q}}_r \end{bmatrix} + \mathbf{M}_{\eta}^{\circ}(\mathbf{q}, \mathbf{q}_0) \ddot{\mathbf{\eta}}_P \\
& + \mathbf{D}_{\eta}^{\circ}(\mathbf{q}, \mathbf{q}_0, \dot{\mathbf{q}}, \dot{\mathbf{q}}_0, \dot{\mathbf{\eta}}_P) \dot{\mathbf{\eta}}_P + \mathbf{K}_{\eta}^{\circ} \mathbf{\eta}_P = \mathbf{0}_{n_{\eta_P} \times 1}
\end{aligned}
\qquad (46)
$$

with:

$$
\left\{
\begin{aligned}
\mathbf{H}^{\circ} &= \mathbf{H}_{0\eta_P}^{L^T} \mathbf{H}_{0}^{L^{-1}} \\
\mathbf{M}^{\circ} &= \begin{bmatrix} \mathbf{H}_{0\eta_P}^{\omega^T} - \mathbf{H}^{\circ} \mathbf{H}_{0}^{\omega^T} & \mathbf{H}_{q_m\eta_P}^{T} - \mathbf{H}^{\circ} \mathbf{H}_{0q_m}^{L} \\ \mathbf{H}_{q_m\eta_P}^{\omega^T} & \mathbf{H}_{0q_m}^{L^T} - \mathbf{H}^{\circ} \mathbf{C}_{0}^{L\omega} \end{bmatrix} \\
\mathbf{D}^{\circ} &= \begin{bmatrix} \mathbf{C}_{\eta_P 0}^{\omega} & -\mathbf{H}^{\circ} \mathbf{C}_{0}^{L\omega} \\ \mathbf{C}_{\eta_P 0}^{\omega} & \mathbf{C}_{q_m0}^{L\omega} \end{bmatrix} \\
\mathbf{D}_{\mathrm{Lr}}^{\circ} &= \begin{bmatrix} \mathbf{C}_{\eta_P 0}^{L} & -\mathbf{H}^{\circ} \mathbf{C}_{0}^{L} \\ \mathbf{C}_{q_r0}^{L} & \mathbf{0} \end{bmatrix} \\
\mathbf{M}_{\eta}^{\circ} &= \mathbf{H}_{\eta_P} - \mathbf{H}^{\circ} \mathbf{H}_{0\eta_P}^{L} \\
\mathbf{D}_{\eta}^{\circ} &= \mathbf{C}_{\eta_P} - \mathbf{H}^{\circ} \mathbf{C}_{0\eta_P}^{L} \\
\mathbf{K}_{\eta}^{\circ} &= \mathbf{K}_{\eta_P}
\end{aligned}
\right.
$$

## 3.2. Control architecture

To both decouple and linearize the task-space, defined by the spacecraft's base attitude and the manipulator's joint states, an NDI framework is designed. However, as it is assumed that no flexible measures are available, one will introduce the notation $\hat{\mathbf{X}}$ to designate a quantity $\mathbf{X}$ obtained in the inverse kinematics/dynamics evaluation from available measures, i.e. $(\mathbf{q}_0, \mathbf{q}, \dot{\mathbf{q}}_0, \dot{\mathbf{q}})$. A control torque that linearizes and decouples the system for a given desired dynamic $\mathbf{v}$ is chosen as:

$$
\tau_{qc} = \mathbf{J}_{\tau}^{\circ} + \left( \hat{\mathbf{M}}^{\circ} \begin{bmatrix} \dot{\mathbf{\omega}}_0 \\ \ddot{\mathbf{q}}_m \end{bmatrix}_d + \hat{\mathbf{M}}^{\circ} \tau_{\infty} + \hat{\mathbf{D}}^{\circ} \begin{bmatrix} \mathbf{\omega}_0 \\ \ddot{\mathbf{q}}_m \end{bmatrix} + \hat{\mathbf{D}}^{\circ}_{\mathrm{Lr}} \begin{bmatrix} \dot{\mathbf{r}}_{S0} \\ \dot{\mathbf{q}}_r \end{bmatrix} \right) \qquad (47)
$$

With the superscript + designating the Moore-Penrose pseudo-inverse. Using the subscript $d$ to indicate a desired quantity, the imposed dynamic is defined as $\tau_{\infty} = \mathbf{K}_{\infty} [\epsilon_c^T \dot{\epsilon}_c^T]^T$ where $\epsilon_c = [\theta_0^T \mathbf{q}_m^T]_d - [\theta_0^T \mathbf{q}_m^T]^T$ is the control error signal of the base attitude $\theta_0$ and manipulator joints rotations/displacements; the control gain $\mathbf{K}_{\infty}$ the $(3+n_m) \times 2(3+n_m)$. Injecting (47) in (45) the control error dynamics are given as:

$$
\begin{aligned}
\ddot{\epsilon}_c = \mathbf{M}^{\circ -1} \left( -\hat{\mathbf{M}}^{\circ} \mathbf{K}_{\infty} \begin{bmatrix} \epsilon_c \\ \dot{\epsilon}_c \end{bmatrix} + \mathbf{M}_{\eta}^{\circ} \ddot{\eta}_P \right. \\
\left. + \mathbf{D}_{\eta}^{\circ} \dot{\eta}_P + (\mathbf{D}^{\circ} - \hat{\mathbf{D}}^{\circ}) \begin{bmatrix} \mathbf{\omega}_0 \\ \ddot{\mathbf{q}}_m \end{bmatrix} + (\mathbf{D}_{\mathrm{Lr}}^{\circ} - \hat{\mathbf{D}}_{\mathrm{Lr}}^{\circ}) \begin{bmatrix} \dot{\mathbf{r}}_{S0} \\ \dot{\mathbf{q}}_r \end{bmatrix} \right) \\
+ (\mathbf{M}^{\circ} - \hat{\mathbf{M}}^{\circ}) \begin{bmatrix} \dot{\mathbf{\omega}}_0 \\ \ddot{\mathbf{q}}_m \end{bmatrix}_d
\end{aligned}
\qquad (48)
$$

The $\mathcal{H}_{\infty}$ feedback aims at providing robustness against modeling uncertainties and the effect of flexible vibrations on the joint control through the couplings with the spacecraft base. We introduce
the notation $\Delta X$ to express the error between the actual quantity $X$ and its measure $\hat{X}$, such that $X = \hat{X} + \Delta X$. Then injecting (46) in (48) allows to express the error dynamics as:

$$ \ddot{\epsilon}_c = \mathbf{M}^{\circ-1} \left( -\hat{\mathbf{M}}^{\circ} \mathbf{K}_{\infty} \begin{bmatrix} \epsilon_c \\ \dot{\epsilon}_c \end{bmatrix} + \mathbf{M}_{\eta}^{\circ} \ddot{\eta}_p \right. \\ \left. + \mathbf{D}_{\eta}^{\circ} \dot{\eta}_p + \Delta \mathbf{D}^{\circ} \begin{bmatrix} \dot{\omega}_0 \\ \dot{\hat{q}}_m \end{bmatrix} + \Delta \mathbf{D}_{Lr}^{\circ} \begin{bmatrix} \dot{\hat{r}}_{S0} \\ \dot{\hat{q}}_r \end{bmatrix} + \Delta \mathbf{M}^{\circ} \begin{bmatrix} \dot{\omega}_0 \\ \ddot{\hat{q}}_m \end{bmatrix}_d \right) \quad (49) $$

In order to perform the gain synthesis only as a function of the desired dynamics, i.e a planned path, expressing the system momentum and artificially adding/subtracting the desired velocities in (49) allows us to reduce the conservatism of the synthesis. However, contrary to the assumption of a fully rigid system (Wilde et al., 2018), flexible elements introduce a dissipative force that results in a loss of kinematic momentum. To account for this momentum dissipation, the kinetic angular momentum, $\delta_{\omega} \in \mathbb{R}^{3 \times 1}$ can be expressed as:

$$ \begin{bmatrix} \mathbf{H}_0^{\omega} & \mathbf{H}_{0qm}^{\omega} \end{bmatrix} \begin{bmatrix} \dot{\omega}_0 \\ \dot{\hat{q}}_m \end{bmatrix} + \begin{bmatrix} \mathbf{H}_0^{\omega L} & \mathbf{H}_{0qr}^{\omega} \end{bmatrix} \begin{bmatrix} \dot{\hat{r}}_{S0} \\ \dot{\hat{q}}_r \end{bmatrix} + \mathbf{H}_{0\eta_p}^{\omega} \dot{\eta}_p = \delta_{\omega} \quad (50) $$

with $\delta_{\omega}$ assumed to be initially null and bounded as for its first time-derivative. Thus, by injecting (50) in (49) and artificially adding the desired angular base/manipulators velocities, one can express the control error dynamics as:

$$ \ddot{\epsilon}_c = \mathbf{M}^{\circ-1} \left( -\hat{\mathbf{M}}^{\circ} \mathbf{K}_{\infty} \begin{bmatrix} \epsilon_c \\ \dot{\epsilon}_c \end{bmatrix} + \mathbf{M}_{\eta}^{\circ} \ddot{\eta}_p \right. \\ \left. + \Delta \mathbf{M}^{\circ} \begin{bmatrix} \dot{\omega}_0 \\ \ddot{\hat{q}}_m \end{bmatrix}_d + \tilde{\mathbf{D}}_{\eta}^{\circ} \dot{\eta}_p + \tilde{\mathbf{D}}^{\circ} \begin{bmatrix} \dot{\omega}_0 \\ \ddot{\hat{q}}_m \end{bmatrix}_d - \tilde{\mathbf{D}}^{\circ} \dot{\epsilon}_c + \tilde{\mathbf{D}}_{\delta}^{\circ} \delta_{\omega} \right) \quad (51) $$

with:

$$ \left\{ \begin{array}{l} \tilde{\mathbf{D}}_{\eta}^{\circ} = \mathbf{D}_{\eta}^{\circ} - \Delta \mathbf{D}_{Lr}^{\circ} \begin{bmatrix} \mathbf{H}_0^{\omega L} & \mathbf{H}_{0qr}^{\omega} \end{bmatrix}^{+} \mathbf{H}_{0\eta_p}^{\omega} \\ \tilde{\mathbf{D}}^{\circ} = \Delta \mathbf{D}^{\circ} - \Delta \mathbf{D}_{Lr}^{\circ} \begin{bmatrix} \mathbf{H}_0^{\omega L} & \mathbf{H}_{0qr}^{\omega} \end{bmatrix}^{+} \begin{bmatrix} \mathbf{H}_0^{\omega} & \mathbf{H}_{0qm}^{\omega} \end{bmatrix} \end{array} \right. \quad (52a) $$

$$ \left\{ \begin{array}{l} \tilde{\mathbf{D}}^{\circ} = \Delta \mathbf{D}^{\circ} - \Delta \mathbf{D}_{Lr}^{\circ} \begin{bmatrix} \mathbf{H}_0^{\omega L} & \mathbf{H}_{0qr}^{\omega} \end{bmatrix}^{+} \begin{bmatrix} \mathbf{H}_0^{\omega} & \mathbf{H}_{0qm}^{\omega} \end{bmatrix} \\ \tilde{\mathbf{D}}_{\delta}^{\circ} = \Delta \mathbf{D}_{Lr}^{\circ} \begin{bmatrix} \mathbf{H}_0^{\omega L} & \mathbf{H}_{0qr}^{\omega} \end{bmatrix}^{+} \end{array} \right. \quad (52b) $$

$$ \tilde{\mathbf{D}}_{\delta}^{\circ} = \Delta \mathbf{D}_{Lr}^{\circ} \begin{bmatrix} \mathbf{H}_0^{\omega L} & \mathbf{H}_{0qr}^{\omega} \end{bmatrix}^{+} \quad (52c) $$

Similarly, the flexible dynamics (46) can be re-written in function of the desired dynamics and the momentum dissipation (50) such as:

$$ \mathbf{M}_{\eta}^{\circ} \ddot{\eta}_p = - \mathbf{M}^{\circ} \begin{bmatrix} \dot{\omega}_0 \\ \ddot{\hat{q}}_m \end{bmatrix}_d - \tilde{\mathbf{D}}^{\circ} \begin{bmatrix} \dot{\omega}_0 \\ \ddot{\hat{q}}_m \end{bmatrix}_d - \tilde{\mathbf{D}}_{\eta}^{\circ} \dot{\eta}_p - \mathbf{K}_{\eta}^{\circ} \eta_p - \tilde{\mathbf{D}}_{\delta}^{\circ} \delta_{\omega} \quad (53) $$

with:

$$ \left\{ \begin{array}{l} \tilde{\mathbf{D}}_{\eta}^{\circ} = \mathbf{D}_{\eta}^{\circ} - \mathbf{D}_{Lr}^{\circ} \begin{bmatrix} \mathbf{H}_0^{\omega L} & \mathbf{H}_{0qr}^{\omega} \end{bmatrix}^{+} \mathbf{H}_{0\eta_p}^{\omega} \\ \tilde{\mathbf{D}}^{\circ} = \mathbf{D}^{\circ} - \mathbf{D}_{Lr}^{\circ} \begin{bmatrix} \mathbf{H}_0^{\omega L} & \mathbf{H}_{0qr}^{\omega} \end{bmatrix}^{+} \begin{bmatrix} \mathbf{H}_0^{\omega} & \mathbf{H}_{0qm}^{\omega} \end{bmatrix} \end{array} \right. \quad (54a) $$

$$ \left\{ \begin{array}{l} \tilde{\mathbf{D}}_{\delta}^{\circ} = \mathbf{D}^{\circ} - \mathbf{D}_{Lr}^{\circ} \begin{bmatrix} \mathbf{H}_0^{\omega L} & \mathbf{H}_{0qr}^{\omega} \end{bmatrix}^{+} \begin{bmatrix} \mathbf{H}_0^{\omega} & \mathbf{H}_{0qm}^{\omega} \end{bmatrix} \\ \tilde{\mathbf{D}}_{\delta}^{\circ} = \mathbf{D}_{Lr}^{\circ} \begin{bmatrix} \mathbf{H}_0^{\omega L} & \mathbf{H}_{0qr}^{\omega} \end{bmatrix}^{+} \end{array} \right. \quad (54b) $$

$$ \tilde{\mathbf{D}}_{\delta}^{\circ} = \mathbf{D}_{Lr}^{\circ} \begin{bmatrix} \mathbf{H}_0^{\omega L} & \mathbf{H}_{0qr}^{\omega} \end{bmatrix}^{+} \quad (54c) $$

One can then define the block diagram for the $\mathcal{H}_{\infty}$ synthesis as given in Fig. 3 and detailed with Fig. 4 and Fig. 5. This presents the advantage of clearly defining the flexible disturbances on the control performances for a set of pre-defined motions.

In the current synthesis, no external disturbances are considered to put the focus on the rejection of flexible vibrations of a flexible elements being manipulated. However, any external torques/forces could be considered in the $\mathcal{H}_{\infty}$ by virtually adding an equivalent torque/force in (45). This term would be included in the synthesis as an additional input pre-multiplied by $\mathbf{M}^{\circ-1}$ in the expression of the control dynamics (49).

## 4. Illustration of the proposed method

The significance of the proposed control (referred as $C_1$) strategies is demonstrated through its implementation on a real-time

<figure><img src="image_4.png" /><figcaption>Fig. 3. $\mathcal{H}_{\infty}$ block diagram synthesis with details provided in Figs. 4 and 5.</figcaption></figure>

simulation platform (Kraïem et al., 2024) and compared with a more classical state-of-the-art separated base/manipulators control (Cumer et al., 2021) (referred as $C_2$). The platform provides a high-fidelity environment with actuator/sensor dynamics and noises as well as a detailed visual rendering of the scene. As detailed in the reports of the H2020 PULSAR project$^2$, in particular in the “D32.1-Detailed Design Document of In-Space Assembly in Simulation Demonstrator (dISAS)” document, the proposed formalism to integrate flexible elements has been validated with model analysis comparing with the Satellite Dynamic Toolbox (SDT) (Alazard et al., 2008). The servicing of a satellite with a flexible appendage is employed as an illustrative use-case, depicted in its initial and final configuration in Fig. 6. It is assumed that a joint-space path has been provided such that the manipulator approaches the target to the servicer base. The SMS platform should maintain a fixed attitude throughout the course of the manipulators’ motions. To address modeling uncertainties, sensor noise, and effectively counteract the undesired effects of appendage vibrations not accounted for in the path-planning method, we present the design of the $\mathcal{H}_{\infty}$ gain, $\mathbf{K}_{\infty}$, followed by simulation results.

### 4.1. Use-case

The considered SMS is composed of a 1960 kg base with a payload of 1440 kg on top, corresponding to the PULSAR telescope base (Garzon et al., 2022). A 6-Degree of Freedoms (DoF) vispa manipulators (Shilton et al., 2024) is used. Base attitude is controlled by six identical reaction wheels, each capable of delivering a maximum control torque of 0.12 $N \cdot m$. Solar arrays attached to the servicer are considered rigid to highlight the impact of a flexible body at the end of a kinematic chain. The serviced target is a scaled-down version of the base and a flexible panel with a total weight of 641 kg and an inertia matrix of diag(1115, 734, 734) kg $\cdot m^2$, which corresponds to one-third of the base’s. The flexible solar array exhibit three modes with a damping coefficient of 0.005 and natural frequencies equals to from 1.3069, 5.7177 and 8.8216 rad/s.

<footnote>$^2$ https://cordis.europa.eu/project/id/821858/results.</footnote>
<figure><img src="image_5.png" /><figcaption>Fig. 4. Detailed block diagram of the tracking control error dynamics (51) in Fig. 3.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. Detailed block diagram of the disturbance dynamics (53) in Fig. 3.</figcaption></figure>

Prior developing the $H_\infty$ gain synthesis, a manipulator motion is defined such that the target is significantly moved. In order to guarantee the feasibility of the maneuver while maintaining a fixed-base orientation, the manipulator motion is slowed down to stay within reaction-wheels capabilities to compensate for base-manipulator couplings (Rognant et al., 2019b).

## 4.2. $H_\infty$ gain synthesis

To demonstrate robustness towards modeling uncertainties/ errors, for the considered motion, the nominal dynamics are

obtained with the proposed formalism (40) while the uncertain dynamics is evaluated with the extension of the SPART toolbox (Virgili-Llop et al., 2016) as developed in Kraïem et al. (2023).

The MATLAB function `systune` is used to synthesize the gain of the structured $H_\infty$ controller. Due to the decoupled nature of (51) provided by the NDI, the structure of the gain $K_\infty$ is chosen as a concatenation of two diagonal square matrices of size $(3+n_m)$, representing a PD controller.

We begin by considering the hard constraints used on the design. Robustness is imposed via constraints on the gain (6 dB).
<figure><img src="image_7.png" /><figcaption>Fig. 6. Snapshots of the SMS servicing scenario.</figcaption></figure>

and phase (45°) margins at the command $\tau_{\infty}$ level. The response time and the static error are fixed through constraints on the sensitivity function at $\tau_{\infty}$ with the weights $W_{s,i}(s)$ defined as $W_{s,i} = (2s + \lambda_i \cdot 1.732e-6)/(s + \lambda_i \cdot 1.732e-3)$, for $i \in \{b, a\}$ representing the base and arm, respectively, and with $\lambda_b = 1$ and $\lambda_a = 100$. This represents required bandwidths of 0.001 Hz and 0.1 Hz for the base and arms, respectively. Finally, the control effort is tuned through constant weights on the transfer $T_{r \to u}$ between the reference $[\dot{\omega}_0^T \ \ddot{q}_m^T]^T$ and $\tau_{\infty}$ for the frequency range [1, 10] rad/s. The weight values, 0.0032 for the base and 398.1072 for the arms, are found by tuning in view of the amplitude of the commands obtained in closed loop.

In addition to these hard constraints, two soft constraints are also considered. First, a bound of 0.0032 is imposed on the transfer $T_{\delta \to \epsilon}$ between the disturbance and the tracking error, to limit the effect of the variation of the kinetic angular momentum due to the flexible elements. Finally, a bound of 31.6228 is imposed on the transfer $T_{r \to \dot{\epsilon}}$ between the reference and the velocity error for the frequency range [1e-10, 1e-3], to ensure a good tracking error at lower frequencies.

The nonconvex design is performed from 4 random initial conditions. To further constrain the problem and prevent the gain $K_{\infty}$ from becoming excessively large, we impose hard limits on its values. The limits used (proportional: 0.0022 for the base and 0.1540 for the arms; derivative 1.1e-6 for the base and 0.0028 for the arms) are obtained by trial-and-error during the controller tuning. The obtained design ensures a value of 1.8423 for the hard constraints, indicating that they are not completely satisfied. The simulation results show that the controller provides good performance and robustness, as seen in the following section.

## 4.3. Results and discussions

Results of the closed-loop simulations are presented in this section. One drawback of developing a Lagrangian formalism is the number of recursive function calls to establish the dynamic model as opposed to a Newton-Euler approach that presents a more optimized approach to compute kinematic and dynamic models (Virgili-Llop et al., 2016). For this reason, in order to reach computational times of 100 Hz on the SpeedGoat target, the extension of the SPART toolbox as detailed in Kraïem et al. (2022) is used. This modeling formalism will only account for couplings between the spacecraft's base and the link connected to a flexible element, as opposed to the formalism developed in this paper where couplings of a flexible element at the end of a kinematic chain are considered for all elements of this rigid chain. This leads to the necessity of introducing robustness to modeling uncertainties in a model based control approach such as

<figure><img src="image_8.png" /><figcaption>Fig. 7. End-effectors accurate pose tracking: upper subplot depicts the end-effector motion around its initial pose for both control strategies; lower subplot depicts the error signals $\epsilon_c(\mathbf{r}_{EE}) = \mathbf{r}_{EEd} - \mathbf{r}_{EE}$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 8. Evolution of the base angular velocities for both control strategies.</figcaption></figure>

the proposed NDI controller. After evaluation of with of the Task Execution Time of the Matlab Simulink running on the Speedgoat, one can conclude on the necessity to implement on the simulation platform (Kraïem et al., 2024) the extended Newton-Euler formalism developed in Virgili-Llop et al. (2016) and extended as in Kraïem et al. (2022) to reach computation times under 100 Hz on the Speedgoat, against more than two times slower simulations with the formalism proposed in this study. Robustness towards modeling errors, i.e. lack of couplings between the manipulator's joints and the flexible body at its end, is thus evaluated. Additionally, all sensors are considered noisy and inertia/mass random uncertainties of ±5% and ±20% are respectively considered for the servicer and the target's bodies.

The efficacy of the proposed robust control strategy is demonstrated in Fig. 7. In the upper subplot, a comparison of the evolution of the end-effector motion around its initial pose is displayed, while the lower subplot illustrates the tracking errors for both control strategies. Both strategies exhibit satisfactory performance for a substantial task involving nearly half a meter of target movement. However, the proposed strategy, $C_1$, showcases slightly less precise tracking than the separated control approach, $C_2$. Similar observations can be made regarding the precise control of the base velocity, as illustrated in Fig. 8. Acceptable and similar performances are obtained for both control strategies. Furthermore, Fig. 9 provides insight into the robustness of the designed and implemented robust control law. This figure offers a comparison of the flexible forces and torques acting on
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Comparison, for both control strategies, of the maximal norm of flexible forces and torques on all axes.</figcaption></figure>

the spacecraft base by depicting the evolution of the maximal forces/torques norm between each axes at every moment. It is evident that both control approaches result in attenuation of the norms. However, a notable observation is the lower level of flexible excitation obtained with $C_1$. This can be explained by the trade off with a precise tracking of the end-effector desired path. Besides illustrating the feasibility of implementing a model-based control strategy on our simulation platform, the designed controller exhibits robustness towards modeling uncertainties, sensor noises and system's internal disturbances torques. Such control framework enables simplification of path-planning methods by not requiring consideration of flexible dynamics.

## 5. Conclusions

This paper introduces a formalism for integrating flexible bodies at the end of rigid multi-kinematic chain systems for free-floating systems. The key advantage of this formalism lies in its explicit definition of the nonlinear couplings between flexible bodies and the rest of the free-floating system. Extending an open-source Matlab toolbox with the proposed formalism enables system analysis, path-planning, and controller design. This tool proves beneficial in designing a robust model-based control law. The proposed framework is based on NDI and structured $H_\infty$ control. The synthesis of the $H_\infty$ controller is performed by designing appropriate weighting functions using the extended Matlab tools. Control performance is assessed under real-time conditions for an OOS scenario with a simulation platform providing highly realistic conditions. The proposed control strategies demonstrate efficiency and robustness against internal disturbances, modeling uncertainties, and sensor noises.

Future work will concentrate on designing control strategies for detumbling targets with flexible behaviors and performing on-orbit deployment of large and reconfigurable structures. Additionally, the Lagrangian formalism will be expanded to address robots walking on flexible elements.

## CRediT authorship contribution statement

**Sofiane Kraïem:** Writing – original draft, Methodology, Formal analysis, Writing – review & editing, Validation, Investigation.
**Mathieu Rognant:** Writing – original draft, Investigation, Writing – review & editing, Methodology.
**Sérgio Waitman:** Writing – review & editing, Methodology, Writing – original draft.
**Idriss Chelikh:** Validation, Software.
**Jean-Marc Biannic:** Supervision, Methodology.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Appendix. Expressions of the partial derivative terms to compute the convective matrix in Eq. (35) and 36

State-derivative expressions are given to detail the convective terms in (35) and 36. To conserve a general formalism the generalized joint-state variable, $\mathbf{q} \in \mathbb{R}^{n_q \times 1}$, defined as $\mathbf{q} = [q_1 \dots q_{n_q}]$, corresponds to the kinematic chain(s) DoFs while $\mathbf{x}_0 \in \mathbb{R}^{6 \times 1}$, defined as $\mathbf{q}_0 = [\theta_{0x} \ \theta_{0y} \ \theta_{0z} \ r_{0x} \ r_{0y} \ r_{0z}]^T = [q_0 \ q_1 \ q_2 \ q_3 \ q_4 \ q_5 \ q_6]^T$, corresponds to the 6-DoFs of the base.

We also introduce the set, $C_{0,l}$, that gathers all links in the branch and prior the $l$th link starting from the original child of the kinematic chain.

### A.1. Base-derivative of inertia matrices

The convective terms $\mathbf{c}_0$, $\mathbf{c}_{0q}$ and $\mathbf{c}_{0\eta_p}$ (36) are detailed as:

$$ \mathbf{c}_0 = \begin{bmatrix} \dot{\mathbf{q}}_0^T \frac{\partial}{\partial \theta_{0x}} (\mathbf{H}_0) + \dot{\mathbf{q}}^T \frac{\partial}{\partial \theta_{0x}} (\mathbf{H}_{0q}^T) + \dot{\eta}_p^T \frac{\partial}{\partial \theta_{0x}} (\mathbf{H}_{0\eta_p}^T) \\ \dot{\mathbf{q}}_0^T \frac{\partial}{\partial \theta_{0y}} (\mathbf{H}_0) + \dot{\mathbf{q}}^T \frac{\partial}{\partial \theta_{0y}} (\mathbf{H}_{0q}^T) + \dot{\eta}_p^T \frac{\partial}{\partial \theta_{0y}} (\mathbf{H}_{0\eta_p}^T) \\ \dot{\mathbf{q}}_0^T \frac{\partial}{\partial \theta_{0z}} (\mathbf{H}_0) + \dot{\mathbf{q}}^T \frac{\partial}{\partial \theta_{0z}} (\mathbf{H}_{0q}^T) + \dot{\eta}_p^T \frac{\partial}{\partial \theta_{0z}} (\mathbf{H}_{0\eta_p}^T) \\ \dot{\mathbf{q}}_0^T \frac{\partial}{\partial r_{0x}} (\mathbf{H}_0) + \dot{\mathbf{q}}^T \frac{\partial}{\partial r_{0x}} (\mathbf{H}_{0q}^T) + \dot{\eta}_p^T \frac{\partial}{\partial r_{0x}} (\mathbf{H}_{0\eta_p}^T) \\ \dot{\mathbf{q}}_0^T \frac{\partial}{\partial r_{0y}} (\mathbf{H}_0) + \dot{\mathbf{q}}^T \frac{\partial}{\partial r_{0y}} (\mathbf{H}_{0q}^T) + \dot{\eta}_p^T \frac{\partial}{\partial r_{0y}} (\mathbf{H}_{0\eta_p}^T) \\ \dot{\mathbf{q}}_0^T \frac{\partial}{\partial r_{0z}} (\mathbf{H}_0) + \dot{\mathbf{q}}^T \frac{\partial}{\partial r_{0z}} (\mathbf{H}_{0q}^T) + \dot{\eta}_p^T \frac{\partial}{\partial r_{0z}} (\mathbf{H}_{0\eta_p}^T) \end{bmatrix} $$

$\mathbf{c}_{0q} =$

$$ \begin{bmatrix} \dot{\mathbf{q}}_0^T \frac{\partial}{\partial \theta_{0x}} (\mathbf{H}_{0q}) + \dot{\mathbf{q}}^T \frac{\partial}{\partial \theta_{0x}} (\mathbf{H}_q) + \dot{\eta}_p^T \frac{\partial}{\partial \theta_{0x}} (\mathbf{H}_{q\eta_p}) \\ \dot{\mathbf{q}}_0^T \frac{\partial}{\partial \theta_{0y}} (\mathbf{H}_{0q}) + \dot{\mathbf{q}}^T \frac{\partial}{\partial \theta_{0y}} (\mathbf{H}_q) + \dot{\eta}_p^T \frac{\partial}{\partial \theta_{0y}} (\mathbf{H}_{q\eta_p}) \\ \dot{\mathbf{q}}_0^T \frac{\partial}{\partial \theta_{0z}} (\mathbf{H}_{0q}) + \dot{\mathbf{q}}^T \frac{\partial}{\partial \theta_{0z}} (\mathbf{H}_q) + \dot{\eta}_p^T \frac{\partial}{\partial \theta_{0z}} (\mathbf{H}_{q\eta_p}) \\ \dot{\mathbf{q}}_0^T \frac{\partial}{\partial r_{0x}} (\mathbf{H}_{0q}) + \dot{\mathbf{q}}^T \frac{\partial}{\partial r_{0x}} (\mathbf{H}_q) + \dot{\eta}_p^T \frac{\partial}{\partial r_{0x}} (\mathbf{H}_{q\eta_p}) \\ \dot{\mathbf{q}}_0^T \frac{\partial}{\partial r_{0y}} (\mathbf{H}_{0q}) + \dot{\mathbf{q}}^T \frac{\partial}{\partial r_{0y}} (\mathbf{H}_q) + \dot{\eta}_p^T \frac{\partial}{\partial r_{0y}} (\mathbf{H}_{q\eta_p}) \\ \dot{\mathbf{q}}_0^T \frac{\partial}{\partial r_{0z}} (\mathbf{H}_{0q}) + \dot{\mathbf{q}}^T \frac{\partial}{\partial r_{0z}} (\mathbf{H}_q) + \dot{\eta}_p^T \frac{\partial}{\partial r_{0z}} (\mathbf{H}_{q\eta_p}) \end{bmatrix} $$

$$ \mathbf{c}_{0\eta_p} = \begin{bmatrix} \dot{\mathbf{q}}_0^T \frac{\partial}{\partial \theta_{0x}} (\mathbf{H}_{0\eta_p}) + \dot{\mathbf{q}}^T \frac{\partial}{\partial \theta_{0x}} (\mathbf{H}_{q\eta_p}) \\ \dot{\mathbf{q}}_0^T \frac{\partial}{\partial \theta_{0y}} (\mathbf{H}_{0\eta_p}) + \dot{\mathbf{q}}^T \frac{\partial}{\partial \theta_{0y}} (\mathbf{H}_{q\eta_p}) \\ \dot{\mathbf{q}}_0^T \frac{\partial}{\partial \theta_{0z}} (\mathbf{H}_{0\eta_p}) + \dot{\mathbf{q}}^T \frac{\partial}{\partial \theta_{0z}} (\mathbf{H}_{q\eta_p}) \\ \dot{\mathbf{q}}_0^T \frac{\partial}{\partial r_{0x}} (\mathbf{H}_{0\eta_p}) + \dot{\mathbf{q}}^T \frac{\partial}{\partial r_{0x}} (\mathbf{H}_{q\eta_p}) \\ \dot{\mathbf{q}}_0^T \frac{\partial}{\partial r_{0y}} (\mathbf{H}_{0\eta_p}) + \dot{\mathbf{q}}^T \frac{\partial}{\partial r_{0y}} (\mathbf{H}_{q\eta_p}) \\ \dot{\mathbf{q}}_0^T \frac{\partial}{\partial r_{0z}} (\mathbf{H}_{0\eta_p}) + \dot{\mathbf{q}}^T \frac{\partial}{\partial r_{0z}} (\mathbf{H}_{q\eta_p}) \end{bmatrix} $$
With the expression of the inertia matrices (32), the derivative in function of $q_{0l}$ ($l \in [1, 6]$) are given by:

$$ 
\left\{
\begin{aligned}
\frac{\partial}{\partial q_{0l}} (\mathbf{H}_0) &= \mathbf{P}_0^T \left[ \begin{array}{cc} \frac{\partial}{\partial q_{0l}} (\mathcal{I}_0) & \mathbf{0}_3 \\ \mathbf{0}_3 & \mathbf{0}_3 \end{array} \right] \mathbf{P}_0 \\
&\quad + \left[ \frac{\partial}{\partial q_{0l}} ((\mathcal{R}_{S_0})^T) \begin{array}{cc} \mathbf{0}_3 & \mathbf{0}_3 \end{array} \right] \bar{\mathbf{M}}_0 \mathbf{P}_0 \\
&\quad + \mathbf{P}_0^T \bar{\mathbf{M}}_0 \left[ \begin{array}{cc} \frac{\partial}{\partial q_{0l}} (\mathcal{R}_{S_0}) & \mathbf{0}_3 \\ \mathbf{0}_3 & \mathbf{0}_3 \end{array} \right] + \sum_{i=1}^{n_q} \left( \frac{\partial}{\partial q_{0l}} (\mathbf{J}_{0i})^T \bar{\mathbf{M}}_i \mathbf{J}_{0i} \right. \\
&\qquad \left. + \mathbf{J}_{0i}^T \frac{\partial}{\partial q_{0l}} (\bar{\mathbf{M}}_i) \mathbf{J}_{0i} + \mathbf{J}_{0i}^T \bar{\mathbf{M}}_i \frac{\partial}{\partial q_{0l}} (\mathbf{J}_{0i}) \right) \\
\frac{\partial}{\partial q_{0l}} (\mathbf{H}_{0q}) &= \sum_{i=1}^{n_q} \left( \frac{\partial}{\partial q_{0l}} (\mathbf{J}_{0i})^T \bar{\mathbf{M}}_i \mathbf{J}_{qi} \right. \\
&\qquad \left. + \mathbf{J}_{0i}^T \bar{\mathbf{M}}_i \frac{\partial}{\partial q_{0l}} (\mathbf{J}_{qi}) + \mathbf{J}_{0i}^T \frac{\partial}{\partial q_{0l}} (\bar{\mathbf{M}}_i) \mathbf{J}_{qi} \right) \\
\frac{\partial}{\partial q_{0l}} (\mathbf{H}_{0\eta_p}) &= \left[ \begin{array}{c} \mathbf{L}_{\eta p_f_1} \mathbf{P}_{f_1} \frac{\partial}{\partial q_{0l}} (\mathbf{J}_{0p_{f_1}}) \\ \vdots \\ \mathbf{L}_{\eta p_{f_{np}}} \mathbf{P}_{f_{np}} \frac{\partial}{\partial q_{0l}} (\mathbf{J}_{0p_{f_{np}}}) \end{array} \right] \\
\frac{\partial}{\partial q_{0l}} (\mathbf{H}_q) &= \sum_{i=1}^{n_q} \left( \frac{\partial}{\partial q_{0l}} (\mathbf{J}_{qi})^T \bar{\mathbf{M}}_i \mathbf{J}_{qi} \right. \\
&\qquad \left. + \mathbf{J}_{qi}^T \bar{\mathbf{M}}_i \frac{\partial}{\partial q_{0l}} (\mathbf{J}_{qi}) + \mathbf{J}_{qi}^T \frac{\partial}{\partial q_{0l}} (\bar{\mathbf{M}}_i) \mathbf{J}_{qi} \right) \\
\frac{\partial}{\partial q_{0l}} (\mathbf{H}_{q\eta_p}) &= \left[ \begin{array}{c} \mathbf{L}_{\eta p_f_1} \mathbf{P}_{f_1} \frac{\partial}{\partial q_{0l}} (\mathbf{J}_{qp_{f_1}}) \\ \vdots \\ \mathbf{L}_{\eta p_{f_{np}}} \mathbf{P}_{f_{np}} \frac{\partial}{\partial q_{0l}} (\mathbf{J}_{qp_{f_{np}}}) \end{array} \right] \\
\frac{\partial}{\partial q_{0l}} (\mathbf{H}_{\eta p}) &= \mathbf{0}_{n_{np}}
\end{aligned}
\right.
$$

### A.1.1. Inertia derivatives expressions

The inertia matrices $\mathcal{I}_0$ and $\mathcal{I}_i$ are expressed in the inertial frame $\mathcal{R}_{ine}$ as given by (9). With (7) and (8), one can detail the recursive expressions as:

$$ 
\left\{
\begin{aligned}
\mathcal{I}_0 &= \mathbf{R}_{A_0} \mathcal{I}_{S_0} \mathbf{R}_{A_0}^T \\
\mathcal{I}_i &= \mathbf{R}_{A_0} \prod_{j \in C_{1,i}} \mathbf{R}_{A_j} \mathcal{I}_{S_i} \left( \mathbf{R}_{A_0} \prod_{j \in C_{1,i}} \mathbf{R}_{A_j} \right)^T \\
\frac{\partial}{\partial q_{0l}} (\mathcal{I}_0) &= \frac{\partial}{\partial q_{0l}} (\mathbf{R}_{A_0}) \mathcal{I}_{S_0} \mathbf{R}_{A_0}^T + \mathbf{R}_{A_0} \mathcal{I}_{S_0} \frac{\partial}{\partial q_{0l}} (\mathbf{R}_{A_0}^T) \\
\frac{\partial}{\partial q_{0l}} (\mathcal{I}_i) &= \frac{\partial}{\partial q_{0l}} (\mathbf{R}_{A_0}) \prod_{j \in C_{1,i}} \mathbf{R}_{A_j} \mathcal{I}_{S_i} \left( \prod_{j \in C_{1,i}} \mathbf{R}_{A_j} \right)^T \mathbf{R}_{A_0}^T \\
&\quad + \mathbf{R}_{A_0} \prod_{j=1}^i \mathbf{R}_{A_j} \mathcal{I}_{S_i} \left( \prod_{j \in C_{1,i}} \mathbf{R}_{A_j} \right)^T \frac{\partial}{\partial q_{0l}} (\mathbf{R}_{A_0}^T)
\end{aligned}
\right.
$$

Therefore, the inertia derivative in function of a base state is reduced to the derivation of the rotation matrix $\mathbf{R}_{A_0}$. This leads to $\frac{\partial}{\partial r_{0x}} (\mathbf{R}_{A_0}) = \frac{\partial}{\partial r_{0y}} (\mathbf{R}_{A_0}) = \frac{\partial}{\partial r_{0z}} (\mathbf{R}_{A_0}) = \mathbf{0}_3$ and with the Euler Angles rotation order given in the derivative in function of the base attitude are:

$$ 
\left\{
\begin{aligned}
\frac{\partial}{\partial \theta_{0x}} (\mathbf{R}_{A_0}) &=
\begin{bmatrix}
0 & c(\theta_{0z})s(\theta_{0y})c(\theta_{0x}) + s(\theta_{0z})s(\theta_{0x}) & 0 \\
0 & -c(\theta_{0z})s(\theta_{0x}) + s(\theta_{0z})s(\theta_{0y})c(\theta_{0x}) & 0 \\
0 & c(\theta_{0y})c(\theta_{0x}) & -c(\theta_{0z})s(\theta_{0y})s(\theta_{0x}) + s(\theta_{0z})c(\theta_{0x}) \\
& & -c(\theta_{0z})c(\theta_{0x}) - s(\theta_{0z})s(\theta_{0y})s(\theta_{0x}) \\
& & -c(\theta_{0y})s(\theta_{0x})
\end{bmatrix} \\
\frac{\partial}{\partial \theta_{0y}} (\mathbf{R}_{A_0}) &=
\begin{bmatrix}
-c(\theta_{0z})s(\theta_{0y}) & c(\theta_{0z})c(\theta_{0y})s(\theta_{0x}) & 0 \\
-s(\theta_{0z})s(\theta_{0y}) & s(\theta_{0z})c(\theta_{0y})s(\theta_{0x}) & 0 \\
-c(\theta_{0y}) & -s(\theta_{0y})s(\theta_{0x}) & 0 \\
& & c(\theta_{0z})c(\theta_{0y})c(\theta_{0x}) \\
& & s(\theta_{0z})c(\theta_{0y})c(\theta_{0x}) \\
& & -s(\theta_{0y})c(\theta_{0x})
\end{bmatrix} \\
\frac{\partial}{\partial \theta_{0z}} (\mathbf{R}_{A_0}) &=
\begin{bmatrix}
-s(\theta_{0z})c(\theta_{0y}) & -s(\theta_{0z})s(\theta_{0y})s(\theta_{0x}) - c(\theta_{0z})c(\theta_{0x}) & 0 \\
c(\theta_{0z})c(\theta_{0y}) & -s(\theta_{0z})c(\theta_{0x}) + c(\theta_{0z})s(\theta_{0y})s(\theta_{0x}) & 0 \\
0 & 0 & c(\theta_{0z})s(\theta_{0x}) - s(\theta_{0z})s(\theta_{0y})c(\theta_{0x}) \\
& & s(\theta_{0z})s(\theta_{0x}) + c(\theta_{0z})s(\theta_{0y})c(\theta_{0x}) \\
& & 0
\end{bmatrix}
\end{aligned}
$$

The derivative of matrix $\mathbf{P}_0$ is then obtained by:

$$ \frac{\partial}{\partial q_{0l}} (\mathbf{P}_0) = \left[ \begin{array}{cc} \frac{\partial}{\partial q_{0l}} (\mathbf{R}_{A_0}) & \mathbf{0}_3 \\ \mathbf{0}_3 & \mathbf{0}_3 \end{array} \right] $$

J$_{0i}$ derivative expression: The derivative of J$_{0i}$ in function of $q_{0l}$ is given $\forall l \in [1, 6]$ by:

$$ \frac{\partial}{\partial q_{0l}} (\mathbf{J}_{0i}) = \frac{\partial}{\partial q_{0l}} (\mathbf{B}_{i0}) \mathbf{P}_0 + \mathbf{B}_{i0} \frac{\partial}{\partial q_{0l}} (\mathbf{P}_0) $$

Considering $\mathbf{r}_{S_i}$ with $i \in [1, n_q]$, the matrix $\mathbf{B}_{i0}$ (14) only depends on the position of the solid $S_i$ in the inertial frame. Thus, $\frac{\partial}{\partial \theta_{0x}} (\mathbf{B}_{i0}) = \frac{\partial}{\partial \theta_{0y}} (\mathbf{B}_{i0}) = \frac{\partial}{\partial \theta_{0z}} (\mathbf{B}_{i0}) = \mathbf{0}_3$ and with (11) $\frac{\partial}{\partial r_{0x}} (\mathbf{B}_{i0}) = \frac{\partial}{\partial r_{0y}} (\mathbf{B}_{i0}) = \frac{\partial}{\partial r_{0z}} (\mathbf{B}_{i0}) = \mathbf{0}_3$

J$_{q_i}$ derivative expression: With (17), the expression of the Jacobian matrix J$_{q_i}$ is given $\forall i \in [1, n_q]$ by:

$$ \mathbf{J}_{q_i} = [\mathbf{B}_{i1} \mathbf{p}_{m_1} \dots \mathbf{B}_{i(i-1)} \mathbf{p}_{m_{i-1}} \mathbf{p}_{m_i} \mathbf{0}_{6 \times (n_q-1)}] $$

As only $\mathbf{p}_{m_i}$ ($i \in [1, i-1]$) depends on the base state, $\frac{\partial}{\partial q_{0i}} (\mathbf{B}_{il}) = \mathbf{0}_6$ and thus:

$$ \frac{\partial}{\partial q_{0i}} (\mathbf{J}_{q_i}) = \left[ \mathbf{B}_{i1} \frac{\partial}{\partial q_{0i}} (\mathbf{p}_{m_1}) \dots \frac{\partial}{\partial q_{0i}} (\mathbf{p}_{m_i}) \mathbf{0}_{6 \times (n_q - i)} \right] $$
Then two cases are considered, either $S_l$ and $S_i$ are in the same kinematic chain or not. Thus, for $i \in [1, n_q]$, if $i \in C_l$ then:

$$ \frac{\partial}{\partial q_{0i}} (\mathbf{p}_{m_l}) = \begin{bmatrix} \frac{\partial}{\partial q_{0i}} (\mathbf{k}_l) \\ \frac{\partial}{\partial q_{0i}} (\mathbf{k}_l^\times) (\mathbf{r}_{S_l} - \mathbf{r}_{A_l}) \end{bmatrix} $$

if not:

$$ \frac{\partial}{\partial q_{0i}} (\mathbf{p}_{m_l}) = \begin{bmatrix} \mathbf{0}_{3 \times 1} \\ \frac{\partial}{\partial q_{0i}} (\mathbf{k}_l) \end{bmatrix} $$

With the recursive expression of $\mathbf{k}_l(8)$, its derivative expression is given $\forall l \in [1, i-1]$ by:

$$ \frac{\partial}{\partial q_{0i}} (\mathbf{k}_l) = \frac{\partial}{\partial q_{0i}} (\mathbf{R}_{A_0}) \prod_{j \in C_{0,i}} \mathbf{R}_{A_j} \begin{bmatrix} \mathbf{0} \\ \mathbf{0} \\ \mathbf{1} \end{bmatrix} $$

with the detail of $\frac{\partial}{\partial q_{0i}} (\mathbf{R}_{A_0})$ given in the previous section and the derivative for of $\mathbf{p}_{m_i}$ is given by:

$$ \frac{\partial}{\partial q_{0i}} (\mathbf{p}_{m_i}) = \begin{bmatrix} \frac{\partial}{\partial q_{0i}} (\mathbf{R}_{A_0}) \prod_{j \in C_{0,i}} \mathbf{R}_{A_j} \begin{bmatrix} \mathbf{0} \\ \mathbf{0} \\ \mathbf{1} \end{bmatrix} \\ + \mathbf{R}_{A_0} \prod_{j \in C_{0,i-1}} (\mathbf{R}_{A_j}) \frac{\partial}{\partial q_{0i}} (\mathbf{R}_{A_i}) \begin{bmatrix} \mathbf{0} \\ \mathbf{0} \\ \mathbf{1} \end{bmatrix} \\ \left( \frac{\partial}{\partial q_{0i}} (\mathbf{R}_{A_0}) \prod_{j \in C_{0,i}} \mathbf{R}_{A_j} \begin{bmatrix} \mathbf{0} \\ \mathbf{0} \\ \mathbf{1} \end{bmatrix} \right)^\times \left( + \mathbf{R}_{A_0} \prod_{j \in C_{0,i-1}} (\mathbf{R}_{A_j}) \frac{\partial}{\partial q_{0i}} (\mathbf{R}_{A_i}) \begin{bmatrix} \mathbf{0} \\ \mathbf{0} \\ \mathbf{1} \end{bmatrix} \right)^\times (\mathbf{r}_{S_l} - \mathbf{r}_{A_l}) \end{bmatrix} $$

### A.1.2. Joint-derivative of inertia matrices

The convective terms $\mathbf{c}_{q0}$, $\mathbf{c}_q$ and $\mathbf{c}_{q\eta_p}$ (36) are detailed as: This convective terms can be developed as:

$$ \begin{cases} \mathbf{c}_{q0} = \begin{bmatrix} \dot{\mathbf{q}}_0^T \frac{\partial}{\partial q_1} (\mathbf{H}_0) + \dot{\mathbf{q}}^T \frac{\partial}{\partial q_1} (\mathbf{H}_{0q}^T) + \dot{\eta}_P^T \frac{\partial}{\partial q_1} (\mathbf{H}_{0\eta_P}^T) \\ \vdots \\ \dot{\mathbf{q}}_0^T \frac{\partial}{\partial q_{n_q}} (\mathbf{H}_0) + \dot{\mathbf{q}}^T \frac{\partial}{\partial q_{n_q}} (\mathbf{H}_{0q}^T) + \dot{\eta}_P^T \frac{\partial}{\partial q_{n_q}} (\mathbf{H}_{0\eta_P}^T) \end{bmatrix} \\ \\ \mathbf{c}_q = \begin{bmatrix} \dot{\mathbf{q}}_0^T \frac{\partial}{\partial q_1} (\mathbf{H}_{0q}) + \dot{\mathbf{q}}^T \frac{\partial}{\partial q_1} (\mathbf{H}_q) + \dot{\eta}_P^T \frac{\partial}{\partial q_1} (\mathbf{H}_{q\eta_P}^T) \\ \vdots \\ \dot{\mathbf{q}}_0^T \frac{\partial}{\partial q_{n_q}} (\mathbf{H}_{0q}) + \dot{\mathbf{q}}^T \frac{\partial}{\partial q_{n_q}} (\mathbf{H}_q) + \dot{\eta}_P^T \frac{\partial}{\partial q_{n_q}} (\mathbf{H}_{q\eta_P}^T) \end{bmatrix} \\ \\ \mathbf{c}_{q\eta_P} = \begin{bmatrix} \dot{\mathbf{q}}_0^T \frac{\partial}{\partial q_1} (\mathbf{H}_{0\eta_P}) + \dot{\mathbf{q}}^T \frac{\partial}{\partial q_1} (\mathbf{H}_{q\eta_P}) \\ \vdots \\ \dot{\mathbf{q}}_0^T \frac{\partial}{\partial q_{n_q}} (\mathbf{H}_{0\eta_P}) + \dot{\mathbf{q}}^T \frac{\partial}{\partial q_{n_q}} (\mathbf{H}_{q\eta_P}) \end{bmatrix} \end{cases} $$

where, $\forall l \in [1, n_q]$:

$$ \left\{ \begin{aligned} & \frac{\partial}{\partial q_l} (\mathbf{H}_0) \\ &= \sum_{i=1}^{n_q} \left( \frac{\partial}{\partial q_l} (\mathbf{J}_{0i})^T \bar{\mathbf{M}}_i \mathbf{J}_{0i} + \mathbf{J}_{0i} \frac{\partial}{\partial q_l} (\bar{\mathbf{M}}_i) \mathbf{J}_{0i} + \mathbf{J}_{0i} \bar{\mathbf{M}}_i \frac{\partial}{\partial q_l} (\mathbf{J}_{0i}) \right) \\ & \frac{\partial}{\partial q_l} (\mathbf{H}_{0q}) \\ &= \sum_{i=1}^{n_q} \left( \frac{\partial}{\partial q_l} (\mathbf{J}_{0i})^T \bar{\mathbf{M}}_i \mathbf{J}_{qi} + \mathbf{J}_{0i} \frac{\partial}{\partial q_l} (\bar{\mathbf{M}}_i) \mathbf{J}_{qi} + \mathbf{J}_{0i} \bar{\mathbf{M}}_i \frac{\partial}{\partial q_l} (\mathbf{J}_{qi}) \right) \\ & \frac{\partial}{\partial q_l} (\mathbf{H}_{0\eta_P}) = \begin{bmatrix} \mathbf{L}_{\eta_P f_1} \mathbf{P}_{f_1} \frac{\partial}{\partial q_l} (\mathbf{J}_{0f_1}) \\ \vdots \\ \mathbf{L}_{\eta_P f_{n_P}} \mathbf{P}_{f_{n_P}} \frac{\partial}{\partial q_l} (\mathbf{J}_{0f_{n_P}}) \end{bmatrix} \\ & \frac{\partial}{\partial q_l} (\mathbf{H}_q) \\ &= \sum_{i=1}^{n_q} \left( \frac{\partial}{\partial q_l} (\mathbf{J}_{qi})^T \bar{\mathbf{M}}_i \mathbf{J}_{qi} + \mathbf{J}_{qi} \frac{\partial}{\partial q_l} (\bar{\mathbf{M}}_i) \mathbf{J}_{qi} + \mathbf{J}_{qi} \bar{\mathbf{M}}_i \frac{\partial}{\partial q_l} (\mathbf{J}_{qi}) \right) \\ & \frac{\partial}{\partial q_l} (\mathbf{H}_{q\eta_P}) = \begin{bmatrix} \mathbf{L}_{\eta_P f_1} \mathbf{P}_{f_1} \frac{\partial}{\partial q_l} (\mathbf{J}_{qf_1}) \\ \vdots \\ \mathbf{L}_{\eta_P f_{n_P}} \mathbf{P}_{f_{n_P}} \frac{\partial}{\partial q_l} (\mathbf{J}_{qf_{n_P}}) \end{bmatrix} \\ & \frac{\partial}{\partial q_l} (\mathbf{H}_{\eta_P}) = \mathbf{0}_{n_{\eta_P}} \end{aligned} \right. $$

In order to keep the generality, the transformation matrix $\mathbf{T}_{A_i, \mathcal{R}_{ine}}$ is differentiated with the recursive relation (7) re-written with (3) $\forall l, i \in [1, n_q]$ as:

$$ \mathbf{T}_{A_i, \mathcal{R}_{ine}} = \mathbf{T}_{A_{i-1}, \mathcal{R}_{ine}} \mathbf{T}_{A_i, i-1} \prod_{j \in C_{i+1,i}} \mathbf{T}_{A_{j-1}, A_j} $$

according to $q_l \forall l \in C_{1,i}$ is given by:

In function of the joint natures, the differentiation of $\mathbf{T}_{A_i, \mathcal{R}_{ine}}$

• If the lth joint is prismatic, then:

$$ \frac{\partial}{\partial \mathbf{q}_l} (\mathbf{T}_{A_i, \mathcal{R}_{ine}}) = \mathbf{T}_{A_{i-1}, \mathcal{R}_{ine}} \prod_{j \in C_{i+1,i}} \mathbf{T}_{A_{j-1}, A_j} $$

• If the lth joint is revolute then $\forall l \in [1, i]$:

$$ \frac{\partial}{\partial \mathbf{q}_l} (\mathbf{T}_{A_i, \mathcal{R}_{ine}}) = \mathbf{T}_{A_{i-1}, \mathcal{R}_{ine}} \begin{bmatrix} -s(q_l) & -c(\alpha_l)c(q_l) & s(\alpha_l)c(q_l) & -a_l s(q_l) \\ c(q_l) & -s(q_l)c(\alpha_l) & -s(q_l)s(\alpha_l) & a_l c(q_l) \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{bmatrix} \prod_{j \in C_{i+1,i}} \mathbf{T}_{A_{j-1}, A_j} $$

Else:

$$ \frac{\partial}{\partial \mathbf{q}_l} (\mathbf{T}_{A_i, \mathcal{R}_{ine}}) = \mathbf{0}_4 $$
Thus, $\forall l \in [1, n_q]$:

$$ \frac{\partial}{\partial \mathbf{q}_l} (\mathbf{T}_{A_i, \mathcal{R}_{ine}}) = \begin{bmatrix} \frac{\partial \mathbf{R}_{A_i, \mathcal{R}_{ine}}}{\partial \mathbf{q}_l} & \frac{\partial \mathbf{r}_{A_i}}{\partial \mathbf{q}_l} \\ \frac{\partial \mathbf{q}_l}{\partial \mathbf{0}_{1 \times 3}} & \mathbf{0} \end{bmatrix} $$

With (7), the differentiation of $\mathbf{T}_{S_i, \mathcal{R}_{ine}}$ is given by:

$$ \frac{\partial}{\partial \mathbf{q}_l} (\mathbf{T}_{S_i, \mathcal{R}_{ine}}) = \frac{\partial}{\partial \mathbf{q}_l} (\mathbf{T}_{A_i, S_i}) + \mathbf{T}_{A_i, \mathcal{R}_{ine}} \frac{\partial}{\partial \mathbf{q}_l} (\mathbf{T}_{A_i, S_i}) $$

where:

$$ \frac{\partial}{\partial \mathbf{q}_l} (\mathbf{T}_{A_i, S_i}) = \frac{\partial}{\partial \mathbf{q}_l} \left( \begin{bmatrix} \mathbf{R}_{A_i} & (\mathbf{r}_{S_i} - \mathbf{r}_{A_i}) \\ \mathbf{0}_{1 \times 3} & 1 \end{bmatrix} \right) \\ = \begin{bmatrix} \frac{\partial}{\partial \mathbf{q}_l} (\mathbf{R}_{A_i}) & \frac{\partial}{\partial \mathbf{q}_l} ((\mathbf{r}_{S_i} - \mathbf{r}_{A_i})) \\ \mathbf{0}_{1 \times 3} & 0 \end{bmatrix} $$

without detailing the differentiation of the matrix $\mathbf{T}_{A_i, S_i}$ depends on the nature of joints and can be expressed as:

$$ \frac{\partial}{\partial \mathbf{q}_l} (\mathbf{T}_{S_i, \mathcal{R}_{ine}}) = \begin{bmatrix} \frac{\partial}{\partial \mathbf{q}_l} (\mathbf{R}_{S_i, \mathcal{R}_{ine}}) & \frac{\partial}{\partial \mathbf{q}_l} (\mathbf{r}_{S_i}) \\ \mathbf{0}_{3 \times 1} & \mathbf{0} \end{bmatrix} $$

Inertia derivatives expressions: The expression of the derivative mass matrix, $\tilde{\mathbf{M}}_i$ in function of the joint $q_l$ is given by:

$$ \frac{\partial}{\partial q_l} (\tilde{\mathbf{M}}_i) = \begin{bmatrix} \frac{\partial}{\partial q_l} (\mathcal{I}_i) & \mathbf{0}_3 \\ \mathbf{0}_3 & \mathbf{0}_3 \end{bmatrix} $$

With the recursive expression of inertias expressed in $\mathcal{R}_{ine}$ (A.3), the derivatives in function of $q_l \forall l \in [1, n_q]$ are given by:

$$ \frac{\partial}{\partial q_l} (\mathcal{I}_i) = \mathbf{R}_{A_0} \frac{\partial}{\partial q_l} \left( \prod_{j \in C_{1,i}} \mathbf{R}_{A_j} \right) \mathcal{I}_{S_i} \left( \prod_{j \in C_{1,i}} \mathbf{R}_{A_j} \right)^T \mathbf{R}_{A_0}^T \\ + \mathbf{R}_{A_0} \prod_{j \in C_{1,i}} \mathbf{R}_{A_j} \mathcal{I}_{S_i} \frac{\partial}{\partial q_l} \left( \left( \prod_{j \in C_{1,i}} \mathbf{R}_{A_j} \right)^T \right) \mathbf{R}_{A_0}^T $$

with the differentiation of $\prod_{j \in C_{1,i}} \mathbf{R}_{A_j}$ previously detailed to obtain

$$ \frac{\partial}{\partial q_l} (\mathbf{J}_{0i}) = \frac{\partial}{\partial q_l} (\mathbf{B}_{i0}) \mathbf{P}_0 \qquad (A.8) $$

$$ \frac{\partial}{\partial q_l} (\mathbf{B}_{i0}) = \begin{bmatrix} \mathbf{0}_3 \\ \left( \frac{\partial (\mathbf{r}_{S_i} - \mathbf{r}_{S_0})}{\partial \mathbf{q}_l} \right)^\times \\ \mathbf{0}_3 \end{bmatrix} = \begin{bmatrix} \mathbf{0}_3 \\ \left( \frac{\partial \mathbf{r}_{S_i}}{\partial \mathbf{q}_l} \right)^\times \\ \mathbf{0}_3 \end{bmatrix} \qquad (A.9) $$

with the differentiation of $\mathbf{r}_{S_i}$ previously detailed to obtain.

J_qi derivatives expression: The Jacobian matrix J_m is differentiated $\forall l, i \in [1, n_q]$ as:

$$ \frac{\partial}{\partial q_l} (\mathbf{J}_{mi}) = \begin{bmatrix} \frac{\partial}{\partial q_l} (\mathbf{B}_{i1} \mathbf{p}_{m_1}) & \dots & \frac{\partial}{\partial q_l} (\mathbf{B}_{i(i-1)} \mathbf{p}_{m_{i-1}}) & \frac{\partial}{\partial q_l} (\mathbf{p}_{m_i}) & \mathbf{0}_{6 \times (n_q - i)} \end{bmatrix} $$

First the joint orientation vector, $\mathbf{k}_k$, (8) is differentiate $\forall k \in [1, i]$ with as:

$$ \frac{\partial}{\partial q_l} (\mathbf{k}_k) = \frac{\partial}{\partial q_l} (\mathbf{R}_{A_k, \mathcal{R}_{ine}}) \begin{bmatrix} \mathbf{0}_3 \\ \mathbf{0} \\ \mathbf{1} \end{bmatrix} $$

According to the nature of the joint $\mathcal{A}_k$ with $\forall k \in [1, i]$:

* If $\mathcal{A}_k$ is prismatic:

$$ \frac{\partial}{\partial q_l} (\mathbf{B}_{ik} \mathbf{p}_{m_k}) = \begin{bmatrix} \frac{\partial}{\partial q_l} (\mathbf{k}_k) \end{bmatrix} $$

* If $\mathcal{A}_k$ is revolute:

$$ \frac{\partial}{\partial q_l} (\mathbf{B}_{ik} \mathbf{p}_{m_k}) = \begin{bmatrix} \frac{\partial}{\partial q_l} (\mathbf{k}_k) \\ \left( \frac{\partial}{\partial q_l} ((\mathbf{r}_{S_i} - \mathbf{r}_{S_k})^\times) \mathbf{k}_k + (\mathbf{r}_{S_i} - \mathbf{r}_{S_k})^\times \frac{\partial}{\partial q_l} (\mathbf{k}_k) \right) \\ + \frac{\partial}{\partial q_l} (\mathbf{k}_k^\times) (\mathbf{r}_{S_i} - \mathbf{r}_{A_k}) + \mathbf{k}_k^\times \frac{\partial}{\partial q_l} (\mathbf{r}_{S_i} - \mathbf{r}_{A_k}) \end{bmatrix} $$

## Data availability

No data was used for the research described in the article.

## References

Alazard, D., Cumer, C., & Tantawi, K. (2008). Linear dynamic modeling of spacecraft with various flexible appendages and on-board angular momentum. In 7th international ESA conference on guidance, navigation & control systems.

Angeles, J., & Lee, S. K. (1988). The formulation of dynamical equations of holonomic mechanical systems using a natural orthogonal complement.

Banerjee, A. (2022). Flexible multibody dynamics: efficient formulations with applications. CRC Press.

Bauchau, O. A. (2011). Flexible multibody dynamics, vol. 176. Springer.

Bremer, H. (2008). Elastic multibody dynamics. Intelligent Systems, Control and Automation: Science and Engineering (ISCA), 13, 85–90.

Carpentier, J., Saurel, G., Buondonno, G., Mirabel, J., Lamiraux, F., Stasse, O., & Mansard, N. (2019). The Pinocchio C++ library: A fast and flexible implementation of rigid body dynamics algorithms and their analytical derivatives. In 2019 IEEE/SICE international symposium on system integration (pp. 614–619). IEEE.

Christidi-Loumpasefski, O. O., Nanos, K., & Papadopoulos, E. (2017). On parameter estimation of space manipulator systems using the angular momentum conservation. In 2017 IEEE international conference on robotics and automation (pp. 5453–5458). IEEE.

Colmenarejo, P., Branco, J., Santos, N., Serra, P., Telaar, J., Strauch, H., Fruhnert, M., Giordano, A., Stefano, M., Ott, C., et al. (2018). Methods and outcomes of the COMRADE project-Design of robust combined control for robotic spacecraft and manipulator in servicing missions: comparison between between Hinf and nonlinear Lyapunov-based approaches. In 69th international astronautical congress.

Craig, R. R., & Bampton, M. C. (1968). Coupling of substructures for dynamic analyses. AIAA Journal, 6(7), 1313–1319.

Cumer, C., Rognant, M., Biannic, J. M., & Roos, C. (2021). Modelling and attitude control design for autonomous in-orbit assembly. In ESA GNC 2021.

Denavit, J., & Hartenberg, R. S. (1955). A kinematic notation for lower-pair mechanisms based on matrices. American Society of Mechanical Engineers.

Fauré, M., Henry, D., Cieslak, J., Colmenarejo, P., & Ankersen, F. (2022). A H∞ control solution for space debris removal missions using robotic arms: the ESA e. Deorbit case. In 2022 UKACC 13th international conference on control (pp. 122–129). IEEE.

Flores-Abad, A., Ma, O., Pham, K., & Ulrich, S. (2014). A review of space robotics technologies for on-orbit servicing. Progress in Aerospace Sciences, 68, 1–26.

Foster, C. L., Tinker, M. L., Nurre, G. S., & Till, W. A. (1995). Solar-array-induced disturbance of the Hubble space telescope pointing system. Journal of Spacecraft and Rockets, 32(4), 634–644.

Géradin, M., & Cardona, A. (2001). Flexible multibody dynamics: a finite element approach.
Giordano, A. M., Dietrich, A., Ott, C., & Albu-Schäffer, A. (2020). Coordination of thrusters, reaction wheels, and arm in orbital robots. *Robotics and Autonomous Systems, 131*, Article 103564.

Giordano, A. M., Garofalo, G., & Albu-Schäffer, A. (2017). Momentum dumping for space robots. In *2017 IEEE 56th annual conference on decision and control* (pp. 5243–5248). IEEE.

Giordano, A. M., Ott, C., & Albu-Schäffer, A. (2019). Coordinated control of spacecraft's attitude and end-effector for space robots. *IEEE Robotics and Automation Letters, 4*(2), 2108–2115.

Hiltz, M., Rice, C., Boyle, K., & Allison, R. (2001). Canadarm: 20 years of mission success through adaptation. In *International symposium on artificial intelligence, robotics and automation, JSC-CN-6877*.

Hirano, D., Fujii, Y., Abiko, S., Lampariello, R., Nagaoka, K., & Yoshida, K. (2014). Simultaneous control for end-point motion and vibration suppression of a space robot based on simple dynamic model. In *2014 IEEE international conference on robotics and automation* (pp. 6631–6637). IEEE.

Hirzinger, G., Landzettel, K., Brunner, B., Fischer, M., Preusche, C., Reintsema, D., Albu-Schäffer, A., Schreiber, G., & Steinmetz, B. M. (2004). DLR's robotics technologies for on-orbit servicing. *Advanced Robotics, 18*(2), 139–174.

Hurty, W. C. (1965). Dynamic analysis of structural systems using component modes. *AIAA Journal, 3*(4), 678–685.

Kraïem, S., Rognant, M., Biannic, J. M., & Brière, Y. (2021). Robust control of rotation-floating space robots with flexible appendages for on-orbit servicing. *IFAC-PapersOnLine, 54*(20), 134–140.

Kraïem, S., Rognant, M., Biannic, J. M., & Brière, Y. (2022). Dynamics and robust control of a space manipulator with flexible appendages for on-orbit servicing. *CEAS Space Journal, 1*–20.

Kraïem, S., Rognant, M., Biannic, J. M., & Brière, Y. (2023). Dynamics and robust control of a space manipulator with flexible appendages for on-orbit servicing. *CEAS Space Journal, 15*(5), 681–700.

Kraïem, S., Rognant, M., & Waitman, S. (2024). Simulation platform to design and validate control laws for a space manipulator system performing on-orbit servicing. In *2024 international conference on space robotics* (pp. 98–103). http://dx.doi.org/10.1109/iSpaRo60631.2024.10687658.

Leckie, F., & Pestel, E. (1960). Transfer-matrix fundamentals. *International Journal of Mechanical Sciences, 2*(3), 137–167.

Meng, D., Wang, X., Xu, W., & Liang, B. (2017). Space robots with flexible appendages: dynamic modeling, coupling measurement, and vibration suppression. *Journal of Sound and Vibration, 396*, 30–50.

Oda, M. (2000). Experiences and lessons learned from the ETS-VII robot satellite. In *Proceedings 2000 ICRA. Millennium conference. IEEE international conference on robotics and automation. Symposia proceedings* (Cat. No. 00CH37065), vol. 1 (pp. 914–919). IEEE.

Papadopoulos, E., Aghili, F., Ma, O., & Lampariello, R. (2022). Robotic manipulation and capture in space. *Frontiers in Robotics and AI, 9*.

Papadopoulos, E., & Dubowsky, S. (1991). Coordinated manipulator/spacecraft motion control for space robotic systems. In *ICRA* (pp. 1696–1701).

Penin, L. F. (2002). Teleoperation with time delay-a survey and its use in space robotics. *Technical Report of National Aerospace Laboratory*.

Pisculli, A., & Gasbarri, P. (2015). A minimum state multibody/FEM approach for modeling flexible orbiting space systems. *Acta Astronautica, 110*, 324–340.

Qiao, J., Wu, H., & Yu, X. (2019). High-precision attitude tracking control of space manipulator system under multiple disturbances. *IEEE Transactions on Systems, Man, and Cybernetics: Systems*.

Roa Garzon, M. A., Koch, C., Rognant, M., Ummel, A., Letier, P., Turetta, A., Lopez, P., Trinh, S., Rodriguez Brena, I. V., Nottensteiner, K., et al. (2022). PULSAR: Testing the technologies for on-orbit assembly of a large telescope. In *16th symposium on advanced space technologies in robotics and automation*. ESA.

Rodrigues, R., Alazard, D., Sanfedino, F., Mauriello, T., & Iannelli, P. (2025). Modeling and analysis of a flexible spinning Euler-Bernoulli beam with centrifugal stiffening and softening: A linear fractional representation approach with application to spinning spacecraft. *Applied Mathematical Modelling, 137*, Article 115699.

Rodrigues, R., Preda, V., Sanfedino, F., & Alazard, D. (2022). Modeling, robust control synthesis and worst-case analysis for an on-orbit servicing mission with large flexible spacecraft. *Aerospace Science and Technology*, Article 107865.

Rognant, M., Cumer, C., Biannic, J. M., Roa, M., Verhaeghe, A., & Bissonnette, V. (2019). Autonomous assembly of large structures in space: a technology review. In *EUCASS 2019*.

Rognant, M., Kraiem, S., & Sasiadek, J. (2019). Kinematic indices of rotation-floating space robots for on-orbit servicing. In *IFToMM world congress on mechanism and machine science* (pp. 3107–3116). Springer.

Rui, X., Zhang, J., Wang, X., Rong, B., He, B., & Jin, Z. (2022). Multibody system transfer matrix method: the past, the present, and the future. *International Journal of Mechanical System Dynamics, 2*(1), 3–26.

Rui, X., Zhang, J., & Zhou, Q. (2014). Automatic deduction theorem of overall transfer equation of multibody system. *Advances in Mechanical Engineering, 6*, Article 378047.

Rybus, T., Barciński, T., Lisowski, J., & Seweryn, K. (2015). Analyses of a free-floating manipulator control scheme based on the fixed-base Jacobian with spacecraft velocity feedback. *Aerospace Robotics II, 59*–69.

Sanfedino, F., Alazard, D., Pommier-Budinger, V., Falcoz, A., & Boquet, F. (2018). Finite element based N-port model for preliminary design of multibody systems. *Journal of Sound and Vibration, 415*, 128–146.

Seddaoui, A., & Saaj, C. M. (2019). Combined nonlinear $H_\infty$ controller for a controlled-floating space robot. *Journal of Guidance, Control, and Dynamics, 42*(8), 1878–1885.

Shabana, A. A. (1997). Flexible multibody dynamics: review of past and recent developments. *Multibody System Dynamics, 1*, 189–222.

Shilton, M., Garland, M., Hackett, C., Allouis, E., Lisle, M., Meringolo, C., Paganini, D., & Hall, A. (2024). VISPA URDF. URL https://github.com/AirbusDefenceAndSpace/visa.

Simeon, B. (2013). Computational flexible multibody dynamics. In *A differential-algebraic approach. Differential-algebraic equations forum*. Heidelberg: Springer.

Umetani, Y., & Yoshida, K. (1989). Resolved motion rate control of space manipulators with generalized Jacobian matrix. *IEEE Transactions on Robotics and Automation, 5*(3), 303–314.

Virgili-Llop, J., Drew, D., Romano, M., Hobson, G., Wakefield, B., & Roberts, W. (2016). Spacecraft robotics toolkit: an open-source simulator for spacecraft robotic arm dynamic modeling and control. In *6th international conference on astrodynamics tools and techniques*.

Wilde, M., Kwok Choon, S., Grompone, A., & Romano, M. (2018). Equations of motion of free-floating spacecraft-manipulator systems: an engineer's tutorial. *Frontiers in Robotics and AI, 5*, 41.

Wu, Y., Han, F., Zheng, M., He, M., Chen, Z., Hua, B., & Wang, F. (2018). Attitude control for on-orbit servicing spacecraft using hybrid actuator. *Advances in Space Research, 61*(6), 1600–1616.

Yoshida, K., & Umetani, Y. (1993). Control of space manipulators with generalized Jacobian matrix. In *Space robotics: Dynamics and control* (pp. 165–204). Springer.