Finite-time state-dependent Riccati equation regulation of anthropomorphic dual-arm space manipulator system in free-flying conditions

Acta Astronautica 216 (2024) 504–517 Available online 23 January 2024 0094-5765/© 2024 The Author(s). Published by Elsevier Ltd on behalf of IAA. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). Contents lists available at ScienceDirect Acta Astronautica journal homepage: www.elsevier.com/locate/actaastro Research paper Finite-time state-dependent Riccati equation regulation of anthropomorphic dual-arm space manipulator system in free-flying conditions Alessandro Scalvini∗, Alejandro Suarez, Saeed Rafee Nekoo, Anibal Ollero The GRVC Robotics Laboratory, Escuela Técnica Superior de Ingeniería, Universidad de Sevilla, Camino de los Descubrimientos s/n, Seville, 41092, Spain ARTICLE INFO Keywords: Dual-arm space manipulator SDRE State-dependent coefficient online computation Optimal control Space robotics Nonlinear control ABSTRACT This paper introduces a novel approach for regulating the pose of a free-flying dual-arm anthropomorphic space manipulator system (SMS) using a finite-time state-dependent Riccati equation (SDRE) controller. The proposed system finds applications in on-orbit satellite inspection, servicing, space structure assembly, and debris manipulation. The dual-arm SMS presented in this work consists of two 7 degrees of freedom (DoF) robotic arms mounted on a free-flying spacecraft, resulting in a complex 20-DoF system. Due to the high number of DoFs, advanced controller design and efficient computations are necessary. The finite-time SDRE controller relies on the state-dependent coefficient (SDC) parameterization matrices, which are nonlinear apparent linearizations of the dynamics. Conventionally, the computation of SDC matrices is offline and relies on the a priori derivation of the analytical equations governing the dynamics of the system. However, this strategy becomes computationally impractical for high DoF plants. To overcome this issue and deliver a more viable solution, a numerical method to construct and update the SDC matrices at each time step is presented. This approach relies on a screw-theory-based recursive Newton–Euler algorithm designed to reconstruct the manipulator inertia and Coriolis matrices. These quantities are the building blocks of the SDC parameters used in the synthesis of the SDRE controller. Simulation results demonstrate the performances of the finite-time SDRE controller augmented with the online update of the state-dependent coefficients. 1. Introduction On-orbit operations present multiple risks and difficulties due to the inherent challenges of working in a harsh space environment. However, a sustainable program of space exploration and exploitation relies on the capability of performing complex tasks such as servicing [1,2], assembling [3], debris removal [4], and repairing directly in orbit [5]. Nowadays, a large number of the most complex operations are conducted by humans during spacewalks. These are hours-long, dangerous, and physically demanding missions in which the astronaut and the environment are separated only by a (relatively) thin space suite. In fact, the absence of gravity, extreme temperature fluctuations, vacuum, radiation, and micro-meteoroids pose significant risks for humans. Because of the mentioned reasons, resorting to space robots to perform Extra Vehicular Activities (EVA) is particularly appealing since they are potentially capable of effectively operating in this hostile environment without endangering human lives. In order to replicate the dexterous manipulation capabilities of human astronauts, it is considered the application of anthropomorphic dual-arm space manipulator systems (SMS) for the realization of maintenance operations on tumbling or cooperative satellites, as illustrated in Fig. 1. ∗Corresponding author. E-mail addresses: [email protected] (A. Scalvini), [email protected] (A. Suarez), [email protected] (S.R. Nekoo), [email protected] (A. Ollero). Motivated by the necessity of helping and gradually replacing humans in the most dangerous operations, the field of space robotics has been evolving since the Shuttle Remote Manipulator System, known as Canadarm I, was first introduced in 1981. Remarkable examples of the advancements in this field include the first on-orbit operation conducted by the ETS-VII [6], the Space Station Remote Manipulator System, also known as Canadarm II, and the Special Purpose Dextrous Manipulator. These last are extensively employed on the International Space Station. An additional step forward consists of employing a Space Manipulator System (SMS) capable of conducting On-Orbit Servicing (OOS) missions with limited (teleoperation) or without (fully autonomous) human supervision [7]. As a matter of fact, SMS can be remotely operated to perform challenging tasks such as replacing faulty components/batteries and assisting the astronauts during spacewalks [8]. Direct telemanipulation, achieved utilizing virtual reality and telepresence techniques (visual and haptic feedback) augments the capability of teleoperation systems by providing a more immersive feeling for the operator [9]. Additionally, advanced rendezvous [10] and docking systems [11] allow the servicing spacecraft to approach https://doi.org/10.1016/j.actaastro.2024.01.027 Received 26 September 2023; Received in revised form 26 December 2023; Accepted 18 January 2024 Acta Astronautica 216 (2024) 504–517 505 A. Scalvini et al. Fig. 1. Render of the Dual arm Space Manipulator System performing a maintenance task. Credits for background image: [16]. and birth to the target satellite, ensuring a secure connection for maintenance operations. With these current solutions, on-orbit servicing is poised to revolutionize space operations by extending the lifespan of satellites [12] and reducing space debris [13], thus improving the sustainability of space missions cost-effectively. The most recent technological advancements in this direction have been extensively reviewed in [1,14,15]. A common trend that can be identified is the necessity to develop a highly agile and dexterous space manipulator system that is capable of performing a wide range of tasks in the different phases of the OOS summarized in [1], namely: (i) rendezvous and closing maneuver, (ii) target identification, (iii) attitude synchronization, (iv) manipulator deployment, (v) capture and (iv) post-capture maneuvers. In this paper, specific attention will be devoted to the coordinated base-arms motion taking place after the target identification. The base and arms adopt a convenient configuration to grasp the target while ensuring it remains within reach of the end effectors. Thus, in this work, the system is required to change both the attitude and position of the base while deploying the dual-arm manipulator in order to reach an overall pose more suitable for capture. Common strategies used for similar scenarios include the Generalized Jacobian [17] and the Virtual Manipulator [18] to solve the inverse kinematics of the SMS (out of the scope of this paper) in order to obtain a reference for the control. Relatively recent techniques include Optimal control, Robust control [19], and Sliding Mode control [20]. The former tends to be excessively computationally demanding and the latter is likely to produce chattering [1]. What is proposed in this paper is a finite-time SDRE augmented with an online update of the state and input matrices. This differs from the prevalent method reported in the literature. Conventional strategies involve using symbolic computation tools for the offline derivation of the dynamics parameters to a priori assemble the analytical statedependent coefficients. The SDC parameterization is then numerically evaluated every time the gain matrix of the SDRE is computed. While the standard approach is well-suited for systems with a limited number of Degrees of Freedom (DoF) [21], it becomes exceedingly demanding as the system’s DoFs increase (𝑛≥10). The limitation arises from the length and complexity of the mass and Coriolis matrices entries, which results from the inherently strong nonlinear dynamics exhibited by space robots. This is summarized in Fig. 2. Another distinctive feature of this work, setting it apart from similar studies, lies in the modeling of the base as ‘‘free-flying’’ rather than ‘‘free-floating’’. Previous works considered the model reference adaptive SDRE and output feedback SDRE to address the control of space manipulators in point-to-point and trajectory tracking scenarios [22–24]. However, the authors restricted Fig. 2. When the system is characterized by a large number of DoF, the analytical form of the dynamics equations results in considerably long and complex symbolic expressions which can make it impractical to convert them in a useful form (as a function handles). The solution adopted is to implement a Newton–Euler (NE) based algorithm to reconstruct the mass and Coriolis matrices at each time instant. The quantities have been transformed into numerical values, making them considerably more manageable. the spacecraft’s motion to free-floating conditions. Furthermore, the model employed in the existing literature on SDRE and space robotics often makes use of strong simplifications to reduce the DoFs. Typically, the case study is limited to planar scenarios as in [23,24]. In this paper, a more complex 20-degree-of-freedom system with a fully actuated freeflying base is considered instead. This additional complexity is reflected in the resulting equations of motion, which become considerably large and difficult to manipulate. As a consequence, the derivation of statedependent coefficient matrices becomes quite challenging, highlighting once more the contribution of this paper. Consequently, this work demonstrates the potential of this nonlinear method for synthesizing the control of complex nonlinear plants while maintaining asymptotic stability. Even though the SDRE belongs to the Optimal control category, it has the great advantage of leading to a closed-loop solution to the optimal control [25–27]. The solution solves the differential Riccati equation (DRE) which can be reliably evaluated (even online) and thus the control policy can be readily retrieved without actually solving an optimization problem. There are several advantages to using this controller. Firstly, it has local asymptotic stability (global for some cases) based on the Lyapunov theorem and it is optimal in nature [25,28]. Moreover, the numerical solution to the DRE can be reliably computed while the solver might not be able to find a solution to Optimal and Robust control problems, especially when the optimization is nonconvex. Another advantage of using this controller is the relative ease of enforcing high-level constraints. This can be done by operating on the weighting matrices of the controller to tune its performances [29]. In conclusion, it is demonstrated the SDRE can be a valuable and versatile tool for controlling complex and highly non-linear systems such as SMSs. The main contributions of this paper can be summarized as follows: 1. Application of the algorithm for reconstructing the dynamics of free-flying systems, as presented in the Ref. [30], to a robot with multiple connected open chains. Implementation of a fully numerical method based on the Newton–Euler algorithm for the online update of the statedependent coefficient matrices. The usual approach based on the offline analytical derivation of the SDC would fail due to the complexity of the system equations of motion. The considerable amount of highly non-linear terms would greatly impair their manipulation, thus making the assembly of the SDC matrices inconvenient. By instead constructing them online in their numerical form, the problem is avoided entirely. Acta Astronautica 216 (2024) 504–517 506 A. Scalvini et al. Table 1 Notation table. GENERIC SYMBOLS Symbol Domain Definition  (⋅)[–] Estimated quantity  (⋅)[–] First time derivative  (⋅)[–] Second time derivative  (⋅)[–] Lie wedge operator (⋅)†[–] Right pseudoinverse {⋅}[–] Reference frame 𝑖NIndex of 𝑖-DoF/link 𝑗NIndex 𝑗th sub-chain 𝑛NNumber of DoF 𝑡R+Time ROBOT MECHANICS Symbol Domain Definition 𝑞𝑗 𝑖R𝑖th DoF ∈𝑗th sub-chain 𝐜(𝐪(𝑡), 𝐪(𝑡)) R𝑛Coriolis and centrifugal forces vector 𝐶(𝐪(𝑡), 𝐪(𝑡)) R𝑛×𝑛Coriolis and centrifugal forces matrix 𝑀(𝐪(𝑡)) R𝑛×𝑛Manipulator inertia matrix (also ‘mass matrix’) 𝐪(𝑡)R𝑛Vector of the generalized DoF (joint variables) {𝐴}𝐫{𝐵}R3Position of {𝐵}with respect to {𝐴} {𝐴}𝐑{𝐵}R3×3 Orientation of {𝐵}with respect to {𝐴} {𝐴}𝐓{𝐵}R4×4 Pose of {𝐵}with respect to {𝐴} 𝜉𝑖(𝑡)R6Twist coordinates of 𝑖th joint/link 𝜏(𝑡)R𝑛Vector of the generalized forces CONTROL THEORY Symbol Domain Definition 𝐴(𝐱(𝑡)) R2𝑛×2𝑛State dependent coefficient (state-related) 𝐵(𝐱(𝑡)) R2𝑛×𝑛State dependent coefficient (input-related) 𝐹R2𝑛×2𝑛Final state weight 𝑋(𝐱(𝑡)) R𝑚×2𝑛SDRE gain matrix 𝑄(𝐱(𝑡)) R2𝑛×2𝑛State trajectory weight 𝑅(𝐱(𝑡)) R𝑛×𝑛Input trajectory weight 𝐮(𝑡)R𝑛Control input vector 𝐱(𝑡)R2𝑛State vector: 𝐱(𝑡) ∶= [𝐪⊤(𝑡), 𝐪⊤(𝑡)]⊤ 2. Implementation of the finite-time SDRE to a 20 Degree of Freedom (DoF), anthropomorphic dual-arm SMS in free-flying mode for coordinated base-arms motion control. The rest of the paper is organized as follows. Section 2presents the system modeling and the algorithm to retrieve the dynamic model and Section 3describes the design of the SDRE controller used in this paper and its synthesis. Section 4provides implementation details of the simulated system and the controller, whereas Section 5presents the simulation results. Finally, the conclusions of this work are discussed in Section 6. 2. System modeling In this section, the adopted notation and main assumptions of the dual arm space manipulator model are presented. After an overview of the system, screw theory is used to formulate the kinematics and the dynamics of the robot since it scales well with the complexity of the model, especially when compared to the Denavit–Hartemberg standard parameterization [31]. Then, the algorithm used to retrieve the manipulator inertia and Coriolis matrices is reported and explained for completeness. 2.1. Notation Table 1 provides a concise summary of the most commonly used symbols. Due to the formal notation, many of these symbols’ definitions will be reiterated for clarity. 2.2. System overview The first critical assumption pertains to the spacecraft operating in a free-flying mode in the microgravity environment of space. Consequently, in accordance to the classification proposed [32], the SMS is capable of employing its actuators for the simultaneous coordinated control of both attitude and position. Consider now the dual arm SMS represented in Fig. 3. The space robot consists of two identical 7-DoF arms with anthropomorphic kinematics and a free-flying 6-DoF base, resulting in a total of 20-DoF system. In particular, the anthropomorphic kinematics of the lightweight dual arm developed in our previous work for aerial manipulation [33] is implemented in this paper as well. It consists of three joints at the shoulder, one at the elbow and three at the wrist, obtaining the following overall configuration: 1. Shoulder flexion/extension. 2. Shoulder adduction/abduction. 3. Medial/lateral rotation. 4. Elbow flexion/extension. 5. Wrist adduction/abduction. 6. Wrist pronation/supination. 7. Wrist flexion/extension. All the parts are assumed to be rigid bodies and all the DoFs can be actuated independently. The mass of the joints is assigned to the links and the latter are assumed to have uniform mass distribution (uniform density). The SMS is displayed in two configurations: Fig. 3(a) shows the home configuration in which all the joint variables are set to zero. This configuration is chosen for convenience in the derivation, as outlined in [31,34]. When the SMS adopts this pose, the center of mass (CoM) of the base is coincident with the origin of the inertial frame {𝑁} since the variables associated with the base DoFs (𝑞0 𝑖for 𝑖= 1,…,6) are considered part of the overall kinematic chain. Thus, they are null at home configuration. Fig. 3(b) shows instead a generic pose of the system. Here, the most important frames can be observed: the inertial frame {𝑁}, the SMS-base frame {0} and the arms-base frames {11} and {12}. The right superscript ‘𝑗’ indicates the sub-chain index. The quantities related to the base are indicated with 𝑗= 0, while those belonging to the right and left arm are indicated with 𝑗= 1 and 𝑗= 2, respectively. Consequently, with 𝑞𝑗 𝑖, is indicated the 𝑖th DoF of the 𝑗th sub-chain, and with 𝑃𝑗 𝑖is indicated the point on the 𝑖th joint axis of the 𝑗th sub-chain. Recall that the axes and the points lying on them are the building blocks used to derive the kinematics using screw theory [31,34]. Finally, Fig. 3(b) illustrates the position vectors of the origins of specific frames (each attached to a distinct component of the manipulator). These frames are configured with their origins coinciding with the CoM of the respective link and are aligned with the principal axes of the corresponding link. The vector pointing form the origin of frame {𝐴}to the origin of frame {𝐵}expressed in {𝐴}is represented with {𝐴}𝐫{𝐵}∈R3. Similarly, {𝐴}𝐑{𝐵}∈R3×3 indicates the rotation matrix expressing the orientation of frame {𝐵}with respect to {𝐴}. Recalling that base DoFs are considered part of the chain, as clarified by Fig. 4, the common chain connecting {𝑁}to {0} can be interpreted as composed of six links and six joints. The first five links of the chain are virtual1while the last is simply the SMS-base itself. From now on, a system with ‘hybrid dynamics’ and ‘tree-like’ kinematics chain will be referred to as an Equivalent Manipulator (EM) model. Note that this interpretation is convenient in the development of ground testing systems that employ robotic manipulators for simulating the 6- DoF motion of a space platform [14,35]. In the EM, the 6-DoF joint connecting the base with the inertial frame is imagined as a series of three virtual prismatic (P) joints followed by three virtual revolute (R) joints. Notice that using the 3R joints to describe the orientation is 1By ‘virtual’ it is meant that the component is point-like and mass-less. Acta Astronautica 216 (2024) 504–517 507 A. Scalvini et al. Fig. 3. (a) Home Configuration. Notice in this configuration the base frame {0} and inertial frame {𝑁}coincide since all the joint variables are null. The 𝑖th DoF belonging to sub-chain 𝑗is indicated with 𝑞𝑗 𝑖. (b) Generic Configuration. With {𝐴}r𝑗 {𝐵}is intended the distance of the origin of the link-frame {𝐵}from {𝐴}(expressed in {𝐴}). Furthermore, both link-frames belong to the 𝑗th sub-chain. The arms are detached from the base only for clarity of visualization. Fig. 4. Kinematic configuration of the anthropomorhpic dual arm space manipulator. The Equivalent Manipulator model employ 6virtual joints to represent the position (𝑞0 1, 𝑞0 2, 𝑞0 3) and attitude (𝑞0 4, 𝑞0 5, 𝑞0 6) of the spacecraft with respect to the Inertial Frame {𝑁}. completely equivalent to the {𝑋, 𝑌 , 𝑍}Euler angles parameterization of the rigid body attitude. Connecting the 3P joint in series with the 3R, as depicted in Fig. 4, the free-flying joint (F) is obtained. The main benefit of this interpretation is that it allows treating the base as six links connected by six joints of a fixed-base open-chain manipulator whose (virtual) base is attached to the inertial frame {𝑁}. Consequently, the screw theory version of the Newton–Euler algorithm presented in [34] (for fixed-based manipulator) and [30] (for a hybrid dynamics system with single open chain) can be readily extended to our dual-arm freeflying space manipulator (see Section 2.4). Finally, it is important to highlight that all effects resulting from Orbital Mechanics are neglected in this study. Consequently, the reference frame denoted as {𝑁}is assumed to be inertial. However, it is important to acknowledge that, given the intended operation in Low Earth Orbit (LEO) or Geostationary Orbit (GEO), non-inertial effects are indeed of significance. In fact, the origin of reference frame {𝑁}would actually coincide with the one of Local Vertical Local Horizontal ({𝐿𝑉 𝐿𝐻}) frame of the target satellite, which is intrinsically non-inertial [36]. The mentioned assumption is however grounded in the particular phase of rendezvous under consideration in this work. In fact, during the terminal phase of rendezvous leading to inspection or docking, the spacecraft trajectory is fully controlled, to the point of being completely arbitrary. In some cases, the trajectory of the chaser (the SMS) can even appear linear in the {𝐿𝑉 𝐿𝐻}of the target, e.g. 𝑉-bar approach [36]. This motion is achieved by momentarily counteracting the influences of orbital mechanics through the active propulsion of thrusters2. In conclusion, 2It is important to emphasize that this operational mode can only be sustained for a very limited duration, owing to its substantial fuel consumption a thorough treatment of the mentioned non-inertial effects is out of the scope of this paper and they are thus neglected by assuming to be able to effectively counteract them as proposed in the literature, see for instance [36–38]. 2.3. Kinematics and dynamics To derive the pose of each relevant mass a screw theory-based formulation of the kinematics is used. Despite the initial effort needed to understand the mathematical formalism, it leads to a quite simple and mechanical procedure [31]. Firstly, the system is placed in a convenient home configuration (see Fig. 3(a)) where each variable associated with a DoF is set to zero. This is a pose that facilitates the identification of the position of strategic points, axes, and orientations of frames with respect to the inertial one {𝑁}. Namely, the position of the CoM of each link, the direction of joint axes, the coordinates of a point lying in them, and the orientation of each link principal axes (PA) are computed with respect to {𝑁}. Since the mentioned quantities are all expressed with respect to the inertial frame {𝑁}only, a proper choice of home configuration greatly simplifies their identification. Once these elements are available, the homogeneous transformation matrix (HTM) {𝑁}𝐓{𝐾}(𝑡) ∈ R4×4, describing the pose of any frame of interest {𝐾}with respect to {𝑁}can be systematically derived by applying the product of exponential formula: {𝑁}𝐓{𝐾}(𝑡) = 𝑘 ∏ 𝑖=1 [𝑒 𝝃𝑖(𝑡)𝑞𝑖(𝑡)]{𝑁}𝐓{𝐾}(0).(1) A remarkable result presented in [31] reveals that there is no need to compute the Lagrange equations to derive the second-order Equation of Motion (EoM) for any system displaying a Lagrangian with the following structure: where 𝑒(⋅)is the exponential matrix,  𝝃𝑖(𝑡) ∈ R4×4 is the twist associated to the 𝑖th joint and 𝑞𝑖(𝑡)is the corresponding joint variable. The operator  (⋅)stands for the Lie-algebra ‘wedge’ [31] that maps the twist coordinates, which are packed into a six-entries vector 𝝃𝑖(𝑡) ∈ R6, to the twist element of the set 𝑠𝑒(3) itself [31]. The pose of {𝐾}relative to {𝑁}when the system is at home configuration is represented by the HTM {𝑁}𝐓{𝐾}(0) ∈ R4×4. The index 𝑘represents the last joint influencing the forward kinematics of frame {𝐾}, which is a generic frame whose pose is of interest (usually the frame attached to the end effector). Thus, index 𝑘indicates the joint immediately preceding frame {𝐾}in the kinematic tree. In order to increase readability, the superscript 𝑗is omitted in the above equation. However, and it is reserved for final stages of the rendezvous process, a scenario that precisely aligns with the focus of this paper. Acta Astronautica 216 (2024) 504–517 508 A. Scalvini et al. it is important to notice that the product of the exponential formula must be applied for each sub-chain 𝑗= 0,1 and 2 independently. A possible optimization could be storing the pose of base sub-chain 𝑗= 0 ({𝑁}𝐓{6}) since it is common to the following DoFs. Indeed the base DoFs influence both arms (thus kinematic chains 𝑗= 1 and 𝑗= 2), but the forward kinematics of one arm is independent of the joint variables of the other. The kinematics of the EM is the starting point for deriving the dynamics since, by applying Eq. (1) for each link and end-effector, the time history of the pose of each relevant mass is obtained. This is needed to compute the kinetic energy of the bodies and hence the Lagrangian. 𝐿=1 2 𝐪(𝑡)⊤𝐌(𝐪(𝑡))  𝐪(𝑡) − 𝑉(𝐪(𝑡)),(2) where 𝐪(𝑡) ∈ R𝑛is the generalized coordinate vector which collects the 𝑛joint variables,  𝐪(𝑡)its first time derivative, 𝑉(𝐪(𝑡)) ∈ Ris the potential energy. Since the operation is conducted in a micro-gravity environment, it is possible to immediately set 𝑉(𝐪(𝑡)) = 0. The quantity 𝐌(𝐪(𝑡)) ∈ R𝑛×𝑛is the symmetric positive definite manipulator inertia matrix defined as: 𝐌(𝐪(𝑡)) = 𝑛 ∑ 𝑖=1 𝐉𝑏 𝑖(𝐪(𝑡))⊤𝑖𝐉𝑏 𝑖(𝐪(𝑡)).(3) The terms of Eq. (3) are the generalized inertia matrix 𝑖∈ R6×6 of the 𝑖th rigid body of the robot and 𝐉𝑏 𝑖(𝐪(𝑡)) ∈ R𝑛×6 is the corresponding body manipulator Jacobian [31]. The role played by this quantity is to ‘rotate’ the matrix 𝑖, which is conveniently expressed in the PA frame of the 𝑖th rigid body, to the current configuration, thus taking into account the information provided by the forward kinematics. This operation can therefore be considered a generalization of Steiner’s theorem [34]. For any mechanical system producing Eq. (2) the second-order dynamics EoMs are readily retrieved: 𝐌(𝐪(𝑡)) 𝐪(𝑡) + 𝐂(𝐪(𝑡), 𝐪(𝑡))  𝐪(𝑡) = 𝝉(𝑡),(4) where 𝝉(𝑡)is the vector of the generalized forces, 𝐌(𝐪(𝑡)) is known from Eq. (3) while 𝐂(𝐪(𝑡), 𝐪(𝑡)) ∈ R𝑛×𝑛is the Coriolis and centrifugal forces matrix. Interested readers can find a comprehensive treatment of this approach in [31] and in [34]. Although the described Lagrangian formulation is very elegant, it leads to computational problems when used for systems with many DoFs, such as the one presented in this work. In particular, the matrices 𝐌(𝐪(𝑡)) and 𝐂(𝐪(𝑡), 𝐪(𝑡)) obtained with the Lagrangian method have a large number of highly complex entries. The complexity is due to the highly non-linear and strongly coupled nature of the system. This leads to products of a very large number of trigonometric functions in almost every term of the mentioned matrices. This hinders both the computational efficiency and the precision of the evaluation of these quantities. Moreover, when the number of DoF becomes sufficiently high (𝑛≥10), the complexity can make it impractical to perform the symbolic derivation, manipulation, and conversion of these matrix functions into a practical and usable form (such as a function handle). Nonetheless, retrieval of the dynamics is imperative for the implementation of a model-based controller. The solution proposed in this paper is to use an implementation of the Newton–Euler (NE) algorithm efficient enough to be run online. 2.4. Algorithm To obtain computational efficiency, the screw theory-based approaches presented in [34] and [30] are adapted to the dual-arm free-flying space manipulator examined in this work. The adopted procedure is composed of inner and outer parts. The inner solves the inverse dynamics with a NE-based approach for a given set of 𝐪(𝑡), 𝐪(𝑡) and  𝐪(𝑡). The outer part calls the inverse dynamics and sets the value of the mentioned quantities to compute the mass matrix 𝐌(𝐪(𝑡)) and the Coriolis vector 𝐜(𝐪(𝑡), 𝐪(𝑡)) ∈ R𝑛. This last is defined as: 𝐜(𝐪(𝑡), 𝐪(𝑡)) = 𝐂(𝐪(𝑡), 𝐪(𝑡))  𝐪(𝑡). Consequently, the inverse dynamics algorithm is run several times with: •𝐪(𝑡) = 𝐪meas(𝑡), 𝐪(𝑡) =  𝐪meas(𝑡)and  𝐪(𝑡) = 𝟎𝑛×1 to compute the Coriolis vector. •𝐪(𝑡) = 𝐪meas(𝑡), 𝐪(𝑡) = 𝟎𝑛×1 and 𝑞𝑖(𝑡) = 1 for 𝑖= 1,…, 𝑛 to iteratively compute the columns of the mass matrix. The subscript ‘meas’ indicates the output of the physical model. In the case of a hardware-in-the-loop simulation, these quantities are replaced by the actual sensor measurements. Additionally, in both cases, the wrench applied at the end effector is set to zero. In Algorithms 1 and 2a simplified pseudo-code is reported to show the strategy without the excessive complexity introduced by a formal notation. Notice that the measurement of the acceleration, which is normally quite noisy, is not needed in this algorithm. In fact, the vector  𝐪(𝑡)has either all components equal to zero or only one unitary (the remaining are still zeros). Algorithm 1 OUTER LOOP- Set quantities 1: Given qmeas(𝑡)and  qmeas(𝑡) 2: procedure SetValues(qmeas(𝑡), qmeas(𝑡)) 3: if M(q(𝑡)) then 4: q(𝑡)←qmeas(𝑡), 5:  q(𝑡)←𝟎𝑛×1 6:  q(𝑡)←𝟎𝑛×1 7: for 𝑖= 1,…, 𝑛 do 8: 𝑞𝑖(𝑡)←1 9: 𝑀(∶, 𝑖)←𝐼𝑛𝑣𝐷𝑦𝑛(q(𝑡), q(𝑡), q(𝑡)) ⊳ 𝑖-th column 10: return M(q(𝑡)) 11: end for 12: else if c(q(𝑡), q(𝑡)) then 13: q(𝑡)←qmeas(𝑡), 14:  q(𝑡)← qmeas(𝑡) 15:  q(𝑡)←𝟎𝑛×1 16: c(q(𝑡), q(𝑡)) ←𝐼𝑛𝑣𝐷𝑦𝑛(q(𝑡), q(𝑡), q(𝑡)) 17: return c(q(𝑡), q(𝑡)) 18: end if 19: end procedure Algorithm 2 INNER LOOP - Newton Euler Inverse Dynamics 1: Given q(𝑡), q(𝑡)and  q(𝑡) 2: procedure ForwardIterations(q(𝑡), q(𝑡), q(𝑡)) 3: for 𝑖= 1,…, 𝑛 do 4: Compute {𝑖}T{𝑖−1}(𝑡)⊳Pose link 𝑖wrt 𝑖− 1 5: Compute 𝝃𝑖(𝑡)⊳Twist link 𝑖 6: Compute  𝝃𝑖(𝑡)⊳Twist Rate link 𝑖 7: end for 8: end procedure 9: return all {𝑖}T{𝑖−1}(𝑡),𝝃𝑖(𝑡), 𝝃𝑖(𝑡) 10: procedure BackwardIterations({𝑖}T{𝑖−1}(𝑡),𝝃𝑖(𝑡), 𝝃𝑖(𝑡)) 11: for 𝑖= 1,…, 𝑛 do 12: Compute 𝑖(𝑡)⊳Wrench link 𝑖 13: Compute 𝜏𝑖(𝑡)⊳Force/Torque Joint 𝑖 14: end for 15: end procedure 16: return 𝝉(𝑡) It is important to notice that the computation of 𝐌(𝐪(𝑡)) is independent of 𝐜(𝐪(𝑡), 𝐪(𝑡)) since they are the two separate case of 1. To prove this consider Eq. (4). When the mass matrix is being computed with the underlined procedure, the joint rates must be set to zero. Thus, the NE Inverse Dynamics is called with  𝐪(𝑡) = 𝟎𝑛×1. This leads to the Coriolis term becoming null: 𝐜(𝐪(𝑡), 𝐪(𝑡)) = 𝐂(𝐪(𝑡), 𝐪(𝑡))  𝐪(𝑡) = 𝟎𝑛×1. Consequently, the left-hand-side of Eq. (4) reduces to the inertial term only: 𝐌(𝐪(𝑡)) 𝐪(𝑡) = 𝝉(𝑡).(5) Acta Astronautica 216 (2024) 504–517 509 A. Scalvini et al. Similarly, when the Coriolis vector is derived the joint accelerations vector that is nullified. Thus the dynamics equation simplifies to: 𝐜(𝐪(𝑡), 𝐪(𝑡)) = 𝝉(𝑡).(6) Furthermore, even the iterations to obtain the columns of the mass matrix are all independent from each other. Recall again that for a given measured configuration 𝐪meas(𝑡), the iterations to compute M(q(t)) require to call the inverse dynamics (Algorithm 2) setting 𝐪(𝑡) = 𝐪meas(𝑡)and  𝐪(𝑡) = 𝟎𝑛×1 except for one single unitary component 𝑞𝑖(𝑡)=1. Substituting into Eq. (5): 𝝉(𝑡) = 𝐌(𝐪meas(𝑡)) [0,…,1,0,…,0]⊤=𝐌(∶, 𝑖), where the dependence from 𝐪meas(𝑡)is omitted for clarity. Thus, at each iteration, the generalized torque output of the Inverse Dynamics (Algorithm 2) is actually the 𝑖th column of the mass matrix. An almost identical reasoning is applied to compute 𝐜(𝐪(𝑡), 𝐪(𝑡)). In conclusion, all these operations can be run in parallel, thus producing a very efficient procedure. In particular, the Inverse Dynamics can be called in parallel 𝑛+ 1 times. The first 𝑛are due to the mass matrix, while the last is for the Coriolis vector. The final step is obtaining the Coriolis matrix from vector 𝐜(𝐪(𝑡), 𝐪(𝑡)) utilizing the right pseudoinverse of the joint rates vector: 𝐂(𝐪(𝑡), 𝐪(𝑡)) = 𝐜(𝐪(𝑡), 𝐪(𝑡))  𝐪†(𝑡).(7) Once these matrices are available, it is possible to build the SDC matrices 𝐀(𝐱(𝑡)) and 𝐁(𝐱(𝑡)) needed for the synthesis of the SDRE, as described in Sections 2.5 and 3, respectively. 2.5. SDC parameterization At this point, there is still a missing link between the dynamics and control. Namely, the dynamics equation must be expressed in a suitable form for the SDRE. This is done by constructing a state-space- like formulation known as SDC parameterization. With this aim, Eq. (4) is first solved for the second derivative of the joint position:  𝐪(𝑡)=−𝐌−1(𝐪(𝑡))𝐂(𝐪(𝑡), 𝐪(𝑡))  𝐪(𝑡) + 𝐌−1(𝐪(𝑡))𝝉(𝑡).(8) By now defining the state vector 𝐱(𝑡) ∶= [𝐪⊤(𝑡), 𝐪⊤(𝑡)]⊤∈R2𝑛, the previous can be written in a quasi-state-space form:  𝐱(𝑡) = [𝟎𝑛×𝑛𝐈𝑛×𝑛 𝟎𝑛×𝑛−𝐌−1(𝐪(𝑡))𝐂(𝐪(𝑡), 𝐪(𝑡))]𝐱(𝑡) + [𝟎𝑛×𝑛 −𝐌−1(𝐪(𝑡))]𝝉(𝑡).(9) Recalling now the definition of a time-invariant non-linear controlaffine system:  𝐱(𝑡) = 𝐟(𝐱(𝑡)) + 𝐠(𝐱(𝑡),𝐮(𝑡)),(10) where state vector is denoted by 𝐱(𝑡) ∈ R2𝑛and the input one by 𝐮(𝑡) ∈ R𝑛. The equilibrium point of the system is zero 𝐟(𝟎) = 𝟎.𝐟(𝐱(𝑡)) ∶ R2𝑛→R2𝑛and 𝐠(𝐱(𝑡),𝐮(𝑡)) ∶ R2𝑛×R𝑛→R2𝑛are vector-valued smooth piecewise-continuous functions satisfying the local Lipschitz condition. The SDC parameterization of system (10) is: 𝐟(𝐱(𝑡)) =𝐀(𝐱(𝑡))𝐱(𝑡), 𝐠(𝐱(𝑡),𝐮(𝑡)) =𝐁(𝐱(𝑡))𝐮(𝑡),(11) where 𝐀(𝐱(𝑡)) ∶ R2𝑛→R2𝑛×2𝑛and 𝐁(𝐱(𝑡)) ∶ R2𝑛→R2𝑛×𝑛. By inspecting Eqs. (9) and (11) it is possible to recognize: 𝐀(𝐱(𝑡)) = [𝟎𝑛×𝑛𝐈𝑛×𝑛 𝟎𝑛×𝑛−𝐌−1(𝐪(𝑡))𝐂(𝐪(𝑡), 𝐪(𝑡))], 𝐁(𝐱(𝑡)) = [𝟎𝑛×𝑛 𝐌−1(𝐪(𝑡))]. (12) The parameters 𝐀(𝐱(𝑡)) and 𝐁(𝐱(𝑡)) are recomputed at each time step of the simulation by numerically updating the matrices 𝐌(𝐪(𝑡)) and 𝐂(𝐪(𝑡), 𝐪(𝑡)) with the NE-based algorithm described in Section 2.4. The SDC is then used to compute the gain 𝐗(𝐱(𝑡)) as described in the next section. 3. Controller design Once the SDC parameterization of the dynamics (Eq. (4)) is available, it is possible to formulate the control problem. The objective of the SDRE control design is to minimize the cost function 𝐽=1 2𝐱⊤(𝑇f)𝐅 𝐱(𝑇f)+ 1 2∫𝑇f 0 {𝐱⊤(𝑡)𝐐(𝐱(𝑡))𝐱(𝑡) + 𝐮⊤(𝑡)𝐑(𝐱(𝑡))𝐮(𝑡)} d𝑡, (13) to obtain a trade-off between the error of the system and the input effort of the control law. In the cost function integral (13),𝐐(𝐱(𝑡)) ∶ R2𝑛→ R2𝑛×2𝑛represents the weighting matrix of states and 𝐑(𝐱(𝑡)) ∶ R2𝑛→ R𝑛×𝑛denotes the weighting matrix for the inputs, symmetric positive semi-definite and definite, respectively. 𝐅∈R2𝑛×2𝑛sets the penalty of states at final time 𝑇f[s]. Two conditions must be met to find the solution to the SDRE: controllability and observability conditions. Condition 1. The SDC set of {𝐀(𝐱(𝑡)),𝐁(𝐱(𝑡))} is a controllable pair of matrices of system (10) for all 𝐱(𝑡) ∈ R2𝑛in 𝑡∈R+[23]. Condition 2. The SDC set of {𝐀(𝐱(𝑡)),𝐐1∕2(𝐱(𝑡))} is an observable pair of matrices of system (10) for all 𝐱(𝑡) ∈ R2𝑛in 𝑡∈R+in which 𝐐1∕2(𝐱(𝑡)) is the Cholesky decomposition of 𝐐(𝐱(𝑡)) [23]. Constructing the Hamiltonian function and applying the stationary condition on that results in the control law: 𝐮(𝑡)=−𝐗(𝐱(𝑡))𝐱(𝑡),(14) The control gain is 𝐗(𝐱(𝑡)) = 𝐑−1(𝐱(𝑡))𝐁⊤(𝐱(𝑡))𝐊(𝐱(𝑡)) in which 𝐊(𝐱(𝑡)) ∶ R2𝑛→R2𝑛×2𝑛is symmetric positive-definite solution to the state-dependent differential Riccati equation (SDDRE): 𝐀⊤(𝐱(𝑡))𝐊(𝐱(𝑡)) + 𝐊(𝐱(𝑡))𝐀(𝐱(𝑡)) −𝐊(𝐱(𝑡))𝐁(𝐱(𝑡))𝐑−1(𝐱(𝑡))𝐁⊤(𝐱(𝑡))𝐊(𝐱(𝑡)) + 𝐐(𝐱(𝑡)) = −  𝐊(𝐱(𝑡)),(15) with final boundary condition 𝐊(𝐱(𝑇f)) = 𝐅. The derivation of the SDRE controller is based on the regulation of the system to zero equilibrium point. However, it is possible to regulate any workspace point by defining the deviation from a constant reference state: 𝛿𝐱(𝑡) = 𝐱ref −𝐱(𝑡),(16) where 𝐱ref is the constant reference state. Consequently, it is possible to redefine the control input vector for the regulation of a state different from zero: 𝐮(𝑡) = 𝐗(𝑡)𝛿𝐱(𝑡).(17) Notice that the minus sign originally present in Eq. (14) is already included in the definition of 𝛿𝐱(𝑡). A block-diagram representation of the computation of the 𝐮(𝑡)is included in Fig. 5. The solution to (15) in this work is chosen the ‘‘Lyapunov-based method’’ with a corresponding negative definite solution to Riccati equation [39]. By considering 𝑡→ ∞[39], the solution to an infinite-horizon version of (15) is: 𝐀⊤(𝐱(𝑡))𝐊− ss(𝐱(𝑡)) + 𝐊− ss(𝐱(𝑡))𝐀(𝐱(𝑡)) −𝐊− ss(𝐱(𝑡))𝐁(𝐱(𝑡))𝐑−1(𝐱(𝑡))𝐁⊤(𝐱(𝑡))𝐊− ss(𝐱(𝑡) + 𝐐(𝐱(𝑡)) = 𝟎,(18) in which 𝐊− ss(𝐱(𝑡)) is the steady-state symmetric negative-definite solution to (18). Subtracting (15) from (18) and using new variable: 𝜅−1(𝐱(𝑡)) = 𝐊(𝐱(𝑡)) − 𝐊− ss(𝐱(𝑡)), result in state-dependent differential Lyapunov equation: 𝜅(𝐱(𝑡)) = 𝐀cl(𝐱(𝑡))𝜅(𝐱(𝑡)) + 𝜅(𝐱(𝑡))𝐀⊤ cl(𝐱(𝑡)) − 𝐁(𝐱(𝑡))𝐑−1(𝐱(𝑡))𝐁⊤(𝐱(𝑡)),(19) where 𝐀cl(𝐱(𝑡)) = 𝐀(𝐱(𝑡)) − 𝐁(𝐱(𝑡))𝐑−1(𝐱(𝑡))𝐁⊤(𝐱(𝑡))𝐊− ss(𝐱(𝑡)), Acta Astronautica 216 (2024) 504–517 510 A. Scalvini et al. Fig. 5. Simplified block diagram of simulation model. with final boundary condition 𝜅(𝐱(𝑇f)) = {𝐅−𝐊− ss(𝐱(𝑡))}−1. The closed-form solution to (19) is found by 𝜅(𝐱(𝑡)) = 𝐄(𝐱(𝑡)) + exp{𝐀cl(𝐱(𝑡))(𝑡−𝑇f)}[𝐏(𝜅(𝑇f)) − 𝐄(𝐱(𝑡))]exp{𝐀⊤ cl(𝐱(𝑡))(𝑡−𝑇f)}, where 𝐄(𝐱(𝑡)) is the solution to algebraic Lyapunov equation: 𝐀cl(𝐱(𝑡))𝐄(𝐱(𝑡)) + 𝐄(𝐱(𝑡))𝐀⊤ cl(𝐱(𝑡)) = 𝐁(𝐱(𝑡))𝐑−1(𝐱(𝑡))𝐁⊤(𝐱(𝑡)). Finally, the solution to (15) is defined by 𝐊(𝐱(𝑡)) = 𝐊− ss(𝐱(𝑡)) + 𝜅−1(𝐱(𝑡)). The application of this controller is described in Section 4. 4. Implementation In this section, the SDC parameterization of the dynamics is derived and analyzed in terms of observability and controllability. Moreover, the SIMULINK-Simscape implementation used in this work is described using the simplified block diagram scheme represented in Fig. 5. 4.1. Observability and controllability To guarantee the existence of a solution to the state-dependent differential Riccati equation (Eq. (15)), it is necessary and sufficient to demonstrate that the system is observable and controllable, as stated by Conditions 1 and 2presented in Section 3. These conditions are satisfied if and only if the Kalman observability and controllability matrices are full rank. For the sake of completeness, recall that the observability (𝐱(𝑡)) ∶ R2𝑛→R2𝑛2×2𝑛and controllability (𝐱(𝑡)) ∶ R2𝑛→R2𝑛×2𝑛2matrices are defined as: (𝐱(𝑡)) = ⎡⎢⎢⎢⎢⎢⎢⎢⎣ 𝐐1∕2 𝐐1∕2𝐀 𝐐1∕2𝐀2 ⋮ 𝐐1∕2𝐀𝑛−1 ⎤⎥⎥⎥⎥⎥⎥⎥⎦ , (𝐱(𝑡)) = [𝐁 𝐀𝐁 𝐀2𝐁…𝐀𝑛−1𝐁]. (20) In the observability matrix, 𝐐=𝐐(𝐱(𝑡)) is used in place of the output matrix since a full-state feedback regulation case is considered [23]. If these matrices are always full rank it implies the linear independence of their columns, thus guaranteeing that they will never be singular. These conditions can be checked throughout the real-world operation by making sure that they are not singular at time 𝑡= 0 and that their ranks never change. Even though these requirements could be computed online during the mission, they would not be of actual help if a singular (unobservable/uncontrollable) configuration is reached while the SMS is operating. This would prevent the computation of the gain matrix 𝐗(𝐱(𝑡)), thus leaving the system uncontrolled and potentially causing a failure of the mission. The mere knowledge of having lost one of these conditions is likely not enough to avoid instability. Consequently, it is of paramount importance to have some guarantees of not being affected by this scenario. In this work, it is decided instead to exploit the SIMULINK-Simscape model for validation through simulation. Indeed, the rank of the matrices (𝐱(𝑡)) and (𝐱(𝑡)) is checked at each time step of the simulation when the SDDRE is solved. A hint on the reason why observability and controllability should be constant properties of the system can be found in the block-structure of the SDC matrices 𝐀(𝐱(𝑡)) and 𝐁(𝐱(𝑡)), defined in Eq. (12). By recalling that the generalized inertia matrix 𝐌(𝐪(𝑡)) is symmetric and positive definite, thus its diagonal component will always be greater than zero. By looking at the lower block 𝐁(𝐱(𝑡)) can be immediately implied that control authority over the system ( 𝐪(𝑡)more specifically) will never be lost. The only remaining issue to be checked is the possibility of −𝐌−1(𝐪(𝑡))𝐂(𝐪(𝑡), 𝐪(𝑡)) to be singular. Given the observations made about 𝐌(𝐪(𝑡)), if a problem is to exist it must lie in the matrix 𝐂(𝐪(𝑡), 𝐪(𝑡)). During the time frame of the simulation, this is never encountered. Even though the discussion above is limited to controllability, the same observation could be extended to observability as well. In conclusion, even though the intuition here described is not as strong as an ‘if and only if’ condition, it can be considered acceptable when united with successful simulation results. 4.2. SIMULINK implementation The SMS dynamics and control are implemented and simulated within SIMULINK environment. In particular, the Simscape built-in application is leveraged to model and simulate multi-body and multidomain physical systems. In this work, Simscape is used to easily assemble and simulate the physical model of the space robot, depicted in Fig. 3. A simplified block diagram describing the overall signal flow as implemented is shown in Fig. 5. The simulation starts at the initial state 𝐱0, which coincides with the Home Configuration of Fig. 3(a). Using the algorithm explained in Section 2.4, the generalized inertia and Coriolis matrices of the SMS are retrieved for the initial conditions and the SDC are assembled as described in Section 2.5. These matrices are then used for the synthesis of the SDRE gain matrix outlined in Section 3. The output of the Finite Time SDRE block is thus the matrix 𝐗(𝐱)described in Eq. (14). More specifically, the MATLAB built-in commands care and lyap are used to numerically solve the Riccati and Lyapunov equations, respectively. Notice that this is the only computationally inefficient step of the proposed controller, not the online computation of 𝐌(𝐪(𝑡)) and 𝐂(𝐪(𝑡), 𝐪(𝑡)). The reason for this bottleneck is that the mentioned functions are not ‘naturally’ implemented in SIMULINK but run on MATLAB engine instead, thus reducing substantially the performances. However, this can be readily fixed by writing the same algorithm in a more efficient language when implemented in real hardware. After the gain matrix 𝐗(𝐱(𝑡)) is available, it is right-multiplied by the current state error 𝛿(𝐱(𝑡)) and the resulting control input 𝐮(𝐱(𝑡)) is given to the physical model. SIMULINK/Simscape assembles and integrates the equations of motions of the actuated system and thus the new current state 𝐱(𝑡)is extracted and fed back to the NE-based algorithm for parameter retrieval. The cycle is then repeated until the final time instant. 5. Simulation results In this section, the results of the simulations are reported and analyzed. Firstly the numerical stability of the algorithm presented in Section 2.4 is evaluated by introducing an open loop case as a benchmark. Then, the system is simulated in a representative case while controlled using the finite time SDRE presented in Section 3to evaluate Acta Astronautica 216 (2024) 504–517 511 A. Scalvini et al. Fig. 6. Open loop joints accelerations. The solid lines represent the time evolution of 𝑞𝑗 𝑖, while the dotted lines stand for the desired accelerations. The reference acceleration is given for the first 5(𝑠)of the simulation which is kept running until 10(𝑠)to show that no significant errors occur afterwards. its performance. Finally, the proposed method is compared against the LQR controller to demonstrate its potentiality and it is tested (without re-tuning) in a scenario with increased dynamic coupling to prove its adaptability. 5.1. Testing the algorithm: Open loop feedback linearization To evaluate the numerical stability of the algorithm employed to retrieve 𝐌(𝐪(𝑡)) and 𝐂(𝐪(𝑡), 𝐪(𝑡)), it is decided to implement an open loop version of the feedback linearization (FL) controller. Consider the following control law: 𝝉FL( 𝐪des(𝑡), 𝑡) =  𝐌(𝐪(𝑡)) 𝐪des(𝑡) +  𝐂(𝐪(𝑡), 𝐪(𝑡))  𝐪(𝑡),(21) where  𝐪des(𝑡)collects the desired time histories of the joints’ accelerations and the symbol  (⋅)indicates that the quantity is an estimation of the real one. Substituting (21) in Eq. (4) results in: 𝐌(𝐪(𝑡)) 𝐪(𝑡) + 𝐂(𝐪(𝑡), 𝐪(𝑡))  𝐪(𝑡) =  𝐌(𝐪(𝑡)) 𝐪des(𝑡) +  𝐂(𝐪(𝑡), 𝐪(𝑡))  𝐪(𝑡).(22) If the dynamics is correctly retrieved by the algorithm a full cancellation of the Coriolis and inertial terms should be achieved, thus obtaining the simple double integrator equation:  𝐪(𝑡) =  𝐪des(𝑡).(23) This is a simple and reliable performance metric: the system can be simulated applying at each time instant the control input 𝝉FL( 𝐪des(𝑡), 𝑡) for a desired acceleration vector and check afterward the extent to which the one obtained from the physical model matches the reference. An implicit assumption of this strategy is assuming the Simscape physical model as the ground truth that the algorithm under analysis aims to reconstruct. This allows us to test if this method can successfully estimate the dynamic properties from a partially known model and the measurements of 𝐪(𝑡)and  𝐪(𝑡). In this work, it is decided to construct  𝐪des(𝑡)such that the corresponding joint rates profiles are parabolic in time, which implies that the joints’ accelerations will be linear in 𝑡. Fig. 7. Open loop joints accelerations error. The solid lines represent the time evolution of the error on the joint acceleration 𝛿 𝑞𝑗 𝑖. The reference acceleration is given for the first 5(𝑠)of the simulation which is kept running until 10(𝑠)to show that no unexpected errors occur afterwards. Thus, the 𝑖th desired joint rate of the 𝑗th sub-chain takes the following simple polynomial form: 𝑞𝑗 𝑖, des(𝑡) = 𝑎0+𝑎1𝑡+𝑎2𝑡2.(24) Among all the possible desired trajectories, the one beginning and ending at rest is selected. Thus, by setting null rate at initial (𝑡= 0) and final time (𝑡=𝑇f), the following conditions apply: 𝑞𝑗 𝑖, des(0) = 0 ⟹𝑎0= 0, 𝑞𝑗 𝑖, des(𝑇f) = 0 ⟹𝑎2= −𝑎1𝑇f.(25) Integrating the Eq. (24) it possible to impose the final amount of displacement of the 𝑖th joint of 𝑗th sub-chain 𝜃𝑗 𝑖, thus finding the last free parameter 𝑎1: 𝑞𝑗 𝑖, des(𝑇f) = 𝜃𝑗 𝑖⟹𝑎1= −6 𝜃𝑗 𝑖 𝑇3 f .(26) Consequently, substituting the coefficients in Eq. (24) and taking the time derivative, the (linear) acceleration profile of the 𝑖th joint variable is obtained: 𝑞𝑗 des,𝑖(𝑡) = −6 𝜃𝑗 𝑖 𝑇3 f (1−2𝑡).(27) In Fig. 6 the results of this test are displayed. The reference is tracked quite accurately since the curves almost overlap and only a small deviation is produced. The results are confirmed by analyzing the evolution of the error of the acceleration displayed in Fig. 7. The first two plots show that the base linear and angular errors remain negligible throughout the simulation. The acceleration error experienced by the arms joints is greater but remains low and almost constant until about 3(𝑠). Afterwards, some components start growing linearly. Under the normal assumption of never operating in an open loop for several seconds, this outcome ensures that the controller’s feedback will prevent the propagation of errors over time. Acta Astronautica 216 (2024) 504–517 512 A. Scalvini et al. Table 2 Final pose joint positions. The parameter 𝜃𝑗 𝑖indicates final amount of displacement of the 𝑖th joint of 𝑗th sub-chain. Final pose joint values Base Arm - R Arm - L Joint Value Joint Value Joint Value 𝜃0 11 [m] 𝜃1 1−60 [deg] 𝜃2 1−60 [deg] 𝜃0 20 [m] 𝜃1 2−25 [deg] 𝜃2 225 [deg] 𝜃0 31 [m] 𝜃1 310 [deg] 𝜃2 3−10 [deg] 𝜃0 440 [deg] 𝜃1 4−45 [deg] 𝜃2 4−45 [deg] 𝜃0 50 [deg] 𝜃1 5−15 [deg] 𝜃2 515 [deg] 𝜃0 6−25 [deg] 𝜃1 625 [deg] 𝜃2 6−25 [deg] NA NA 𝜃1 70 [deg] 𝜃2 70 [deg] Table 3 Weight matrices of the finite time SDRE cost functional. Weighting Symbol Value State 𝐐diag{𝟏𝟎6×1,𝟖14×1,𝟏20×1} Final State 𝐅20𝐐 Control Effort 𝐑diag{𝟏6×1 × 10−3,𝟏𝟎14×1} 5.2. Regulation case Once the reliability of the algorithm presented in Section 2.4 is proven, the loop is closed to validate the capability of the finite time SDRE described in Section 3to regulate the SMS pose. The space robot is required to move from the home configuration discussed in Section 2 and depicted in Fig. 3(a) to a final pose whose values are reported in Table 2. Since a regulation scenario is under analysis, the desired final rates are simply equal to zero: 𝑞𝑗 𝑖(𝑇f) =  𝜃𝑗 𝑖= 0 ∀ 𝑖, 𝑗. (28) In this paper, the simulation time interval is set to 𝑡= [0,10](s) and is carried out using the weighting matrices (Eq. (13)) reported in Table 3. The tuning process of the Finite Time SDRE controller involves adjusting these three key coefficient matrices: 𝐐,𝐅, and 𝐑. Iterative refinement is often necessary to enhance the controller’s performance [40]. Usually, these matrices are selected as diagonal, in order to be able to gain some a priori intuition of their influence on the SDRE. In this case, the 𝐐matrix is linked to the state and demands an increase in its elements to minimize errors in regulation or tracking. Conversely, reducing the elements of the 𝐑matrix enhances the precision of the controller. The first half of the diagonal elements of 𝐐pertain to the joints’ positions, while the second half relates to their rates. Thus, an analogy can be drawn with the influence of proportional and derivative gains in a PD controller. Finally, the 𝐅matrix, weighting on the final state, is typically chosen to be proportional to 𝐐, with the coefficient expressing the relative importance of 𝐅compared to 𝐐in shaping the controller’s behavior. Careful adjustment of these matrices is crucial for achieving adequate system performance. In this work, the standard approach described before is used to realize the controller, but newer strategies are under development. For instance, in [41] the implementation of (SDRE) control algorithms through the use of artificial neural networks is presented with the aim of reducing the computational cost of this controller. Also in [42] it is proposed a combination of fuzzy logic and neural networks is leveraged to approximate the SDRE control. The evolution of the DoFs belonging to the base (𝑗= 0) during regulation is shown in Fig. 8. In the left plots, it is possible to notice that each joint variable approaches and meets the reference at 𝑡=𝑇f. The lack of overshooting and the negligible steady state error indicate that Fig. 8. Closed loop base joints positions and rates. The solid lines represent the time evolution of 𝑞0 𝑖and 𝑞0 𝑖, while the dotted lines stand for the reference values. Fig. 9. Closed loop arms joints positions and rates. The solid lines represent the time evolution of 𝑞𝑗 𝑖and 𝑞𝑗 𝑖, while the dotted lines stand for the reference values. the controller is capable of producing precise and well-damped closedloop dynamics. In the plots on the right, the trajectories of the joint rates can instead be appreciated. As it is possible to see, they start and end at zero following a bell-shaped trajectory. This shows that the SMS base has reached a stable final pose. Similar observations can be made when the arm joints are studied. In Fig. 9 the closed loop evolution in time of the arms’ joint angles and