Suboptimal LQR-based spacecraft full motion control: Theory and experimentation

.


Introduction
As spacecraft technology has evolved, robust and efficient automated control has become an essential mission capability.Over the last fifty years, in fact, worldwide aerospace research environments have addressed their work to the optimization of guidance, navigation and control performances, also paying attention to the propellant consumption and time-to-launch costs.
It is clear, then, how high efficiency controls, ensuring both position accuracy and propellant optimization, have become a critical issue in the control systems scope and specifically in the aerospace sector.
The spacecraft six degrees of freedom optimal control problem has been widely analyzed in the literature, lately focusing on spacecraft relative motion, usually requiring numerical methods [1].Spacecraft formation and the optimization of relative maneuvers are becoming increasingly important topics of investigation.This is due to the benefits Contents lists available at ScienceDirect journal homepage: www.elsevier.com/locate/actaastro in cost, responsiveness, and flexibility of a multi-spacecraft system versus the classical monolithic satellite.
Particularly new is the use of continuous on-off engines, appearing on small spacecraft.This control constraint adds new complexity to finding the optimal solution.In fact, most of the literature on spacecraft optimal control assumes that the thrust can be finely modulated, partially to mitigate the aforementioned problem.Unfortunately, this is not the case with real engines, which are usually limited to some sustained value for thrust.As a partial response to this problem, a new methodology has been presented in [2] with the aim to control spacecraft rendezvous maneuvers assuming multilevel continuous thrusters and impulsive thrusters on the same vehicle.Furthermore, very recently, the interests of the Department of Defense have been focusing on time/propellant optimal rendezvous and capture maneuvers of a noncooperative target satellite [3][4][5].This research has been pushing the envelope with regards to fast computation of practical optimal/sub-optimal trajectories.
The recent numerical approaches to the optimal control problem could be in short classified in two main categories: 1. the indirect methods which employ the calculus of variations to obtain the first-order optimality conditions [6], where the resulting boundary-value problem is sometimes impossible or time-consuming to solve; 2. direct methods which approximate the trajectory via parameterization, and transform the cost functional into a cost function [7][8][9].
These second methods appear to be a viable tool for real-time spacecraft optimal control.However, the major issues with direct methods relate to the difficulties of defining parameters to represent feasible trajectories.chosen parameterization may lead, for example, to unsatisfactory final conditions, or, it may satisfy the final conditions but violate some of the constraints on the controls.Adding penalties to the cost function may help, but it provides no guarantee of improvement as in [5,9].Also, the trajectory computation must be as fast as possible, since it is iterated by an outer optimizer looping on the parameters.The rendezvous problem to a tumbling object analyzed in [5] is a very good example, where a direct method may not converge to feasible solutions or may converge in an unacceptable amount of CPU time for on board implementation.Addressing the linear quadratic control problem, LQRbased algorithms were widely used to control space systems, specifically satellites and spacecraft assemblies.Yang proved the effectiveness of a quaternion-based LQR method for the design of nonlinear spacecraft control systems in [10], demonstrating how the designed controller globally stabilized the nonlinear spacecraft system, whereas it locally optimized the spacecraft performance.
Beatty, in [11], successfully demonstrated the superiority of a linear quadratic regulator with respect to a proportional-derivative (PD) controller.His work consisted of a comparison of spacecraft attitude control methodologies that use reaction wheels for torque actuation and star trackers to infer spacecraft orientation and angular rate.
Another proof of LQR reliability and effectiveness was given by Walker and Spencer [12]: this work focused on a system for relative navigation and automated proximity operations for a microsatellite using continuous thrust propulsion and low-cost visible and infrared imagers.In this case the state error was employed, together with the thrust vector, to define the cost function to be minimized by the LQR optimal gains; moreover the parametric weighting matrices defined as functions of the actual distance to the goal and of the maximum available thrust, while Pulse-Width-Modulation was used to determine thrust control.
Bevilacqua et al. worked on multiple spacecraft control by developing an autonomous distributed control algorithm for close proximity operations of multiple spacecraft systems, including rendezvous and docking scenarios.In order to validate the proposed control approach, both theoretical simulations [13] and experimental tests [14][15][16] were carried out, choosing the LQR as the system controller.
The LQR control effort served as the attractive force toward goal positions, while APF repulsive functions provided collision avoidance for both fixed and moving obstacles.Previous experiments, by the above authors, assumed that each spacecraft is equipped with an attitude control system, thus focusing only on the translational control.Regarding the control strategy, which represents the main subject of the present paper, as in [12], the weighting matrices were defined introducing a parametric structure depending on the maximum allowable values of the states and control effort.Simulations proved the LQR/ APF to be both effective and efficient in conducting simultaneous spacecraft missions and in disturbance rejection.
The previous results obtained in [13] were then implemented and verified in [14], where both theoretical developments and experimental validation of the hybrid LQR/ APF were presented.In this case propellant consumption was sub-optimized in real-time through re-computation of the LQR at each sample time, while the APF performed collision avoidance and a high level decisional logic.
There is a gap that has not been addressed directly in the past literature: real-time techniques for optimal control of both attitude and position.This is partly due to the fact that attitude and position are commonly dealt with independently; satellites primarily control their attitude to communicate to ground stations, to point to targets, to collect measurements or absorb energy from the Sun, whereas translational control is activated infrequently as it is only required during specific maneuvers such as orbit insertion and station-keeping.Nevertheless, current spacecraft rendezvous and docking maneuvers require the cooperation of multiple spacecraft deliberately designed to work together; hence both attitude and position need to be controlled constantly and simultaneously to maintain the formation and avoid collisions.
The present paper, in keeping with this context, illustrates a control algorithm taking into account both position and attitude, thus allowing a global control of a full sixdegree-of-freedom spacecraft with a unique controller, in view of an implementation on multiple spacecraft assembly for the execution of rendezvous and docking maneuvers.
The adoption of a single controller would simplify the system's architecture, reducing the computational burden, without sacrificing required performances.Specifically, the control strategy is based on a LQR, whose gains are recomputed at each time sample, hence generating a suboptimized solution of the six-degree-of-freedom problem.A comparison with an optimal control problem solver, based on a Radau pseudospectral Gaussian quadrature method, constitutes the theoretical validation of the proposed algorithm, evaluating the real optimization level.Finally, the presented LQR-based control system was also experimentally verified through a series of tests carried out on the test bed at the Advanced Autonomous Multiple Spacecraft (ADAMUS) laboratory [17,18].This work particularly showcased the feasibility of a real-time approach.
The main contribution of this research to the state-of-theart on spacecraft optimal relative motion control consists in the application of a SDARE approach to obtain a real-time control algorithm for six-degree-of-freedom spacecraft.SDARE approaches have been successfully applied in the past to several engineering problems [19][20][21][22] but, to the authors' knowledge, the complete spacecraft control problem has not received the attention it deserves.The paper completely characterizes the methodology via a comparison of its performances versus those of a general purpose software for solving nonlinear optimal control problems: GPOPS-II.GPOPS-II is based on an adaptive Radau pseudospectral Gaussian quadrature method.Furthermore, this work presents experimental validation of the proposed optimizer, through hardware-in-the-loop experimentation on a full sixdegree-of-freedom test bed, showing real-time feasibility.
This paper is organized as follows.Section 2 illustrates the equations of the six-degree-of-freedom spacecraft motion, the dynamics linearization and deals with the control strategy by means of a linear quadratic regulator.Section 3 introduces the comparison with the optimal control software GPOPS-II, paying specific attention to the transformation of the results in order to have a meaningful matching with those from the LQR controller, thus ensuring an effective and truthful comparison.Section 4 illustrates the robotic spacecraft simulator employed for the experimental verification.Section 5 presents the results of the experimental tests.Section 6 concludes the paper.

Theory
The complete dynamics of spacecraft relative motion, as it will be expanded in (2.1) and (2.2), consists of the Euler equations and the Clohessy-Wiltshire equations, describing the spacecraft attitude and translation respectively.When dealing with the spacecraft simulator, used for the experiments, the same equations have been used, simplifying the Clohessy-Wiltshire equations to a double integrator, thus obtaining a direct correspondence between an orbiting spacecraft's dynamics and the robot dynamics itself.The resulting dynamics is represented by linear equations for the translation and nonlinear Euler's equations for the rotation.In the next subsections we will simplify the dynamics to represent the spacecraft simulator and we will use the terms spacecraft (S/C) simulator or robot interchangeably.
Turning to detail, the position and attitude motion is represented by the following twelve component vector: three relative position coordinates: x, y, z; the relative three linear velocity components: a set of three Euler angles: θ x , θ y , θ z ; three angular velocity components associated with the rotation rate about the three axes constituting the body reference frame: ω x , ω y , ω z .
Fig. 1 summarizes the conventional reference frames (body B, Local Vertical/Local Horizontal LVLH, inertial E) commonly employed to describe the relative motion of a spacecraft with respect to a generic target and expresses the sets of translational and rotational coordinates mentioned above.
The most common solution employed to control spacecraft entails a set of body-mounted thrusters providing unidirectional forces.This requires thruster mapping strategies, converting the tri-axial forces and torques, obtained from the equations of motion, into a series of unidirectional forces correspondent to each employed thruster.

Translational dynamics
The translational equations of motion of the spacecraft simulator have been derived from the more general Clohessy-Wiltshire equations describing the relative motion of a chase vehicle with respect to a target vehicle lying on a circular orbit of radius R 0 and orbital rate n where m is the mass of the vehicle, that is, in the present work, the spacecraft simulator, g is the gravitational acceleration and F thrust is the control force, [23] illustrates how to derive the linear equations of relative motion, obtaining the well known expressions: The spacecraft simulator is not orbiting, thus implying the coincidence of the local reference frame with the inertial one in the laboratory: hence the terms related to the angular rate n, clearly null, vanish, finally obtaining: In the sequel, the forces will be assumed to be expressed in the coordinate system fixed to the spacecraft, requiring a rotation matrix to obtain then in the inertial reference frame.This choice is convenient to simplify the procedure of thruster mapping.

Attitude dynamics
When the reference frame is the laboratory's floor, the spacecraft simulator attitude dynamics can be modeled employing the Euler's rotational equation of absolute motion in vector/dyadic form, Since the principal axes are also a central triad for the robot, that is the reference frame is baricentric and the inertial tensor is diagonal, expanding along the three axial directions all the products of inertia cancel out and Eq. ( 6) becomes The kinematics equation, using the Euler angles parametrization with chosen rotation sequence yxz, is [23] where e j i are the basis vectors of different reference frames (b body frame) the subscript of which defines the rotation axis, whereas the superscript refers to the partial configuration obtained by each rotation.Expanding along the three axial components x, y, z, and then introducing the matrix notation, the previous equation assumes the form where ci and si represent the shorten form of cos ðθ i Þ and sin ðθ i Þ respectively.

The linearized dynamics
A linearization process has been applied to the nonlinear rotational dynamics of the robot to obtain the matrices representing the required inputs of the LQR to determine the optimal Kalman gain.The following derivations address decoupled translational and rotational motions, and the LQR solver generates independently the sub-optimal force and torque 3-component vectors to be mapped to the body mounted thrusters.Despite the ability of the LQR solver to generate directly the 12-component vector of the forces to be generated by the thrusters, thus addressing the coupled rototranslational dynamics, the authors chose to work with decoupled dynamics to retain the ability of a comparison with GPOPS-II.In fact, GPOPS-II has proven to have several convergence difficulties with the coupled dynamics, and to be very sensitive to the desired initial and final conditions.Given the twelve component state vector X expressed by where x ¼ fx; y; zg T , _ x ¼ f_ x; _ y; _ zg T , θ ¼ fθ x ; θ y ; θ z g T and ω ¼ fω x ; ω y ; ω z g T , the decoupled rotational and translational equation of motion can be summarized with the following vectorial expression: where B is a constant matrix.This represents a multi-input multi-output (MIMO) nonlinear system, a function both of the 12-dimensional state vector X and the 6-dimensional input F , including forces F and torques M. The system is then linearized at the desired state, at every time step.Proceeding in this way for every time step, the whole trajectory will finally consist of a sequence of linearized sections, and will approach the non-linear one as the time sample tends to zero.In light of these considerations the linearized equations are expressed as where the subscript des specifies the desired state terms and the matrix A is the state Jacobian matrix computed at the linearization point: Introducing the actual error vector as Thus, given the rotation sequence yxz, a symbolic calculation tool was used to derive the parametric expression of the matrix A.
In particular the state matrix A is composed of a translational and a rotational term, then gathered generating the global 12 Â 12 state matrix: where The expression of the state matrix rotational term has been further partitioned due to complexity of its elements, allowing us, to focus on those terms depending on the yxz convention, that is on the Euler angles: A rot 1 ¼ The 12 Â 6 input matrix B (where 12 represents the number of the states, while 6 represents the number of the force and torque inputs) remains unvaried, and is given by where Moreover, it was assumed that all the states are observed output.This model provides the inputs of the LQR that given the problem dimensions (12 states and 6 inputs) provides an optimal Kalman gain as a 6 Â 12 matrix to be multiplied by the 12 Â 1 state error X err , finally obtaining the optimal solution F LQR .This represents the optimal input generalized force, including both forces and torques, necessary to reach the desired state X des , while minimizing the cost functional.However the spacecraft simulator is moved by twelve cold gas thrusters, placed on the upper part of the simulator; hence the generalized optimal force has been converted into a set of twelve unidirectional thruster commands by means of a mapping matrix H, as proposed in [24].According to Fig. 2, thrusters are placed on the edge of the four arms symmetrically connected to the basis of the attitude stage and extended upwards and downwards respectively: specifically eight thrusters are placed with an inclination of 45°w ith respect to the sides of the imaginary square (side 2d ¼ 0:64 m) circumscribed to the basis, while the remaining four (thrusters 3, 6, 9, and 12) provide pure vertical motion when the basis is placed horizontally.The thrust distribution matrix is then given by where H f and H m represent the force and torque terms respectively.
with p ¼ cos ð451Þ ¼ sin ð451Þ and d ¼ 0:32 m.The twelve resulting normalized continuous inputs, U cont , have been then computed according to the following expression: pinvðÞ is the pseudoinverse operator, whereas u a is a dimensionless parameter whose magnitude corresponds to the nominal thrust force F t ¼ 0:3 N, introduced to normalize the values.After the pseuoinverse is computed a pulsewidth-modulation strategy is used to finally derive a set of twelve dimensionless binary inputs, u 10 , to be sent to the onboard on/off thrusters.It is worth underlining that the pseudoinverse solution is applied to a constant matrix, posing no computational or solution existence issues.

Linear quadratic control
The control strategy is sub-optimized since the LQR problem is solved at each time step, using dynamically sized weighting matrices Q and R, adapting gains with respect to the actual position reached by the robot.Thus the negative effects deriving from the linearization will be minimized since the whole trajectory will consist of a sequence of linearized sections and will approach the nonlinear one as the sample time tends to zero.
The definition of the control law requires the introduction of a state feedback gain whose design is a trade-off between the transient response and the control effort.The optimal control approach [25] to this design trade-off is to define and minimize a performance index; furthermore, since generally the robot has to reach non-zero target position and attitude, a non-zero set point optimal control [26] has been considered, shifting the actual state X by the desired quantity X des , obtaining the error, defined in Eq. ( 34): Consequently, the shifted optimal regulation problem becomes with A and B from Eqs. ( 18) and ( 27) respectively, and where the state weighting matrix Q is assumed symmetric and positive semi-definite and the input weighting matrix R symmetric and positive definite.In this case, the LQR generates the optimal solution U Ã ðtÞ (all the following functions are listed as functions of time, as they depend on X des which depends on time) where KðtÞ is often referred to as Kalman gain and PðtÞ is the solution of the Riccati differential equation, obtaining the following linear state feedback control: The presented control strategy has been simulated employing a Simulink integrated MATLAB function, which is based on the original code employed in the 'lqr' MATLAB command [27].It requires the following input matrices: A (dynamic matrix), B (control matrix), C (state-output mapping matrix), D (control-output mapping matrix or feedforward matrix), and Q and R weighting matrices.Lastly, it provides the optimal gain matrix K, the solution of the Riccati equation S and the closed-loop eigenvalues of the matrix ðA À BKÞ, allowing us to verify whether or not the system has been stabilized, obtaining negative real part eigenvalues.The presented function executes the routine, in both linux and win32/64 environment where the correspondent C code is generated for compilation under RTAI linux.
Concerning the weighting matrices, the tuning technique is based on the use of constant elements q i and r i for both the matrices Q and R, as part of a trial and error iterative process where each weight and the factor ρ are gradually adjusted with the aim of obtaining the desired performances, that is the reaching of the target position.
Moreover, two additional gains K pos and K rot are introduced in the Q , and the weighting factor ρ becomes variable throughout the simulations as well as the experimental test, depending on these two new gains, thus increasing the control sensitivity to the system evolution.The deriving matrices have the following structure: Specifically, two controls, over the spacecraft position and attitude respectively, allow these new gains to switch between two calibrated values, on the basis of the actual position and attitude: once the i-th component of the position or attitude gets into an expressly chosen band, the boundaries of which have been defined as percentages of the relative target component (i.e.720% of the target position), a counter variable starts to run; when the counter reaches a predefined value the rela- tive K gain (K transl or K rot ) switches to a second value appropriately calibrated for the steady state phase (robot near the target position); these new K gains are applied to the position and atti- tude terms of the Q (Eqs.( 42) and ( 43)), thus forcing a change of the weights as the state approaches the target; when both the position and the attitude gains K pos and K rot have changed, the factor ρ is also allowed to change, hence switching to a new value which is more suitable for the steady state phase.

Comparison with GPOPS-II
A validation of the proposed sub-optimal controller, implemented in Matlab and Simulink, is carried out by comparing it with an optimizer to highlight how well the proposed controller approximates the best solution.At the time of writing, several open-source, freeware and commercial optimal control and nonlinear programming interfaces are available.However, results accuracy, advancements in mesh refinement and generality of problem formulation represent some of the reasons behind the choice of GPOPS-II as the basis for comparison with the Simulink model.GPOPS-II available at [28] has been developed at the University of Florida in cooperation with the U.S. Office of Naval Research (ONR) and the U.S. Defense Advanced Research Projects Agency (DARPA).Turning to detail, it employs an adaptive Radau pseudospectral Gaussian quadrature method and a sparse finite-differencing is used to estimate all first and second derivatives required by the nonlinear programming (NLP) solver.Moreover, it has been designed to work with the NLP solvers Sparse Nonlinear OPTimizer (SNOPT) and Interior Point OPTimizer (IPOPT) and to be extremely flexible, allowing a user to formulate a wide variety of applications including engineering, economics, and medicine.
The optimal control problem is then stated as a NLP problem (refer to Eqs. ( 12) and (44), whose differential algebraic equations are collocated using nodes obtained from a Gaussian quadrature, whereas the state and the control are parameterized using Legendre polynomials and their linear combinations.As stated in [29], the combination of Legendre polynomials with Gaussian quadrature ensures exponential convergence for smooth solution problems.More specifically, in order to find the optimal solution, that is a combination of state and input vector solutions minimizing the objective function in the Bolza form [25], the general purpose optimal control software GPOPS-II uses a set of Legendre Gauss Radau (LGR) collocation points: this set is defined on the domain ½À1; 1, but contains only one of the endpoints.According to [29,30], the whole iterative procedure, employed by GPOPS-II to determine the optimal solution, could be efficaciously divided into five steps representing partial results of the whole method: 1. identification of the first-order optimality conditions of the continuous Bolza problem, on the basis of the Pontryagin Minimum Principle [25]; 2. Radau pseudospectral discretization of the continuoustime first-order optimality conditions of the continuous Bolza problem; 3. Radau pseudospectral discretization of the continuous time optimal control problem, resulting in a discrete NLP; 4. statement of the Karush-Kuhn-Tucker (KKT) conditions related to the NLP; 5. costate estimation obtained from the results of steps 3 and 4.
A new model of the spacecraft simulator has then been developed, introducing the same rotational and translational equation of motion (12), but following, this time, according to the GPOPS-II logic design [31], a different problem formulation.Two phases have been considered, in order to analyze both the transient response and the steady state one, introducing the appropriate linkage constraints in the endpoint function to ensure continuity of the solution and set time and event criteria, defining the transit from the previous phase to the next one.In this way, during the mesh iterations, when an optimal solution is found, the mesh is analyzed in each of these phases, verifying if the mesh error tolerance is satisfied; if not, the solver automatically increases the number of collocation points, updating the mesh refinement and generating a non uniform adaptive grid designed to reduce the error.
Finally, the following aspects have been considered to obtain a correct comparison between the approach proposed herein and GPOPS one.
1.The comparison of the LQR method and GPOPS has been focused on the global cost required by the controller.Specifically the global dimensionless cost J, that is the cost or objective functional to be minimized both by the LQR and the optimal solver, consists of two terms P and F, each related to the generalized position and velocity accuracy, and the control effort respectively.Since the Simulink model employs a fixed time step of 0.02 s and a discrete solver, the cost function has been computed using the following discrete expression, to match the same type of cost function generated by GPOPS.This summation approximates the integral cost function.
In this analysis more importance has been attributed to the propellant cost, setting the tuning parameters in order to save as much propellant as possible.2. Since both the optimal solver GPOPS and the LQR solver cannot handle binary variables, the propellant term has been computed using the continuous control vector U cont ; hence, during these simulations, the block used in the Simulink model to transform the continuous control into a binary series of on/off commands has been bypassed.3.For the sake of simplicity the constant tuning technique (refer to Section 2.4) has been applied in both models, using also the same numeric values, thus to observe the results on equal weighting terms.
One last expedient has been devised to make the comparison even more truthful: the phase of postprocessing transforms the results obtained with GPOPS, interpolating data and then extracting values on the basis of the equally spaced time vector used in Simulink.This technique allows a point by point comparison, hence a more detailed analysis of the results and their discrepancy.
The most important data affecting the outcome of this comparison, hence the Simulink model performances, were the simulation costs.A set of eight simulations has then been executed both using the Simulink and the GPOPS software, generating Table 1, which summarizes the cost values as well as the initial and desired condition for each of these eight tests.It is worth reminding that the translational dynamics is linear, while the rotational is nonlinear.The initial and final values were randomly chosen to cover a wide range of maneuvers.
From a statistical analysis it can be inferred that on average the sub-optimal controller requires 15% additional cost.However, although the GPOPS controller seems to be the best alternative, some observations have to be taken into account before expressing the final evaluation: Both Simulink and GPOPS simulations have been run using a laptop with standard performances, at the time of writing (Processor Intel Core i3, RAM 2 GB).Given a fixed time simulated duration of 200 s in either case, the optimal solver definitely incurred longer computation time; specifically, setting a mesh tolerance between 1e-5 and 1e-3 (suggested values 1e-6, default value 1e-3), each simulation required on average 600-700 s, three times the simulated time, while, using Simulink, simulations required 200-500 s.Even though the LQR appears slightly faster in simulation, this comparison does not show such a substantial gap between the two techniques yet.There is, in fact, an additional fundamental observation that needs to be made: while GPOPS needs to compute an entire trajectory before returning a solution, the proposed approach can solve the LQR problem in real time, i.e., it decides how to actuate for the immediate next time step.
The LQR solution constitutes a linear state feedback control, generating a closed loop system, while the multi-purpose optimal control software GPOPS-II provides a guidance, that is the determination of an entire desired trajectory, which implies incapability to execute in real time.In fact, a 600 s convergence time for a 200 s simulation clearly shows inability for real-time implementation.
To summarize, the Simulink theoretical controller represents the best alternative: in particular the feedback nature and the higher level of flexibility of the controller allow the use of the sub-optimal solution for real-time implementation; furthermore the previous cost analysis ensures that the same sub-optimal controller maintains a high optimization level, since the costs discrepancy is 15%, compensated by implementability with a real spacecraft.4. The six-degree-of-freedom spacecraft simulator testbed Once the reliability and the effective optimization level of the Simulink model were proven, a series of experimental tests were conducted on the ADAMUS test bed [32,33].
This test bed consists of a six-degree-of-freedom spacecraft simulator based on an air bearing technology that allows it to move as it would in space, i.e., torque and force free.The so-called moving platform (MP) is controlled by twelve cold gas thrusters, appropriately placed on the upper part of the simulator.The MP is powered by two Lithium-ions batteries which are connected to an on board power management system.Moreover, a 13 Â 15 ft (3.96 Â 4.57 m) flat epoxy floor (Fig. 3) is being used as the main base upon which to move the robot.
The overall system represented in Fig. 4 consists of two main stages: The translational stage (TS) provides near frictionless planar translational motion by means of three linear air bearings, to transform operative conditions (1-g) into nearmicrogravity ones; the near gravity free vertical translation is instead obtained introducing a system of near-frictionless air bearing pulleys and variable-mass counterbalances.
The attitude stage (AS) represents the real spacecraft simulator and is connected to the TS through a spherical segment air bearing; it gives rise to the roll, pitch, yaw degree-of-freedom and hosts the onboard computer and position tracking system housing together with the thrusters supports.This test bed belongs to the most complex category of simulators, the category of combination systems, which integrates the capabilities of translational systems with those of rotational systems.
The ADAMUS robot, unlike all the other existing six-degreeof-freedom systems [34][35][36][37][38][39][40], ensures more realistic simulations, hence results, thanks to the dynamical reproduction of motion along the full 6 degrees of freedom and in particular along the vertical translational degree of freedom.Another interesting feature of the ADAMUS platform is the flexibility: by simply substituting the attitude stage, different categories of spacecraft could be tested, considerably extending its scope of applications.Additionally, special attention has been given to the mass balancing system: a Balancing Platform (BP) allows the robot to modify, by means of three linear motors, the AS center of mass position.

Translational stage
The translational state (TS), built by Guidance Dynamics Corporation s (GDC), is made up of an horizontal basis, connected to the linear air bearings providing the near-frictionless contact with the epoxy floor, and of a column with a sleeve for the vertical translational motion, at the top of which a spherical air bearing cup and its correspondent spherical segment ball allow the TS to be connected with the AS.Moreover, the TS contains two groups of tanks: the first group, placed on the lower base of the TS, supplies the compressed air to the air bearings, while the second one, placed on the intermediate base, feeds the air pulleys for the near gravity free vertical motion.

Attitude stage
The attitude stage (AS), designed and built by the ADAMUS lab, provides 740°about the pitch and roll axes and 360°of yaw motion.The main body of the attitude stage consists of a discoid basis of composite material (fiber glass and high density foam) containing the slots of the onboard computer, the power management system, the motor drives and the controller card.Four arms are symmetrically connected to the basis and extend two upwards and two downwards respectively.Three thrusters Fig. 3.The reference frame in the ADAMUS laboratory [32].Fig. 4. The two stages and the main components of the moving platform [32].
per arm are located on their edges, along three mutually orthogonal directions.
Furthermore, both tanks, providing the compressed air to the thrusters, and Lithium-ion batteries are connected to the lower arms for the sake of the AS stability.Lastly, 6 LEDs are distributed over the 4 arms and the AS base; these LEDs, together with the puck, are essential components of the PhaseSpace Impulse System, an optical motion tracker, designated to the position tracking.

Real-time experimentation
The guidance, navigation and control algorithms were created in Simulink, generating an executable file for Real-Time Application Interface (RTAI) Linux.The Simulink file used to program the experiments is based on the ideal Simulink file employed for the comparison with GPOPS, but has several differences listed in the following.
A set of s-functions was used to interface algorithms and hardware on the robot, providing a way to store, save   or modify the desired variables and to send the suboptimal control vector to the actuators' drivers.
A linear Kalman filter (LKF) and an extended Kalman filter (EKF) [41,42] provided a way to estimate the translational state (position and linear velocities) and the rotational state (Euler angles and the angular velocities) respectively, in such a manner that the error is minimized statistically by means of two-step algorithms (predict/    Lastly, previous model employed to control the ADA-MUS spacecraft simulator used angles error threshold between 8 and 5°, while in the present work the threshold was set to only 2°for the angles and 2 cm for the positions; these assumptions require higher control performances, but, at the same time, allow the controller to increase position accuracy reducing the steady state oscillations.Table 2 summarizes the initial and desired conditions of the first full motion test, Table 3 specifies the weighting matrices values, while Fig. 5 describes its main results.
Since the vertical motion range is small on this first test, an additional representation of the vertical degree of freedom y is provided separately in Fig. 6: A visualization of the experiment is presented in Figs.7 and 8.
Figs. 9-15 describe the complete response, including the error which correctly reaches the zero value for each  component of the state vector or remains within the relative threshold.This set of full motion tests allowed the authors to study how the real controller could handle the translational and rotational dynamics deriving differences: the translational dynamics is much slower than the rotational one, reaching the target position in about eighty seconds.On the other hand, this behavior influences also the attitude control: much more thrust is in fact required to let the robot translate, hence stressing the attitude control that in this way easily tends to oscillate or possibly overshoot.Consequently the attitude control comes out to be the most critical component; nevertheless, thanks to the last tuning adjustments, overshoot and oscillations have been greatly contained obtaining satisfactory results.
In this case, in fact, the three controlled angles reach the desired values in about twenty seconds, exhibiting only minimal steady state oscillation.Figs.11 and 12 illustrate the translational velocity and angular velocity.Since the desired values of both velocity terms have been set to zero, the velocity error comes out to be equal to the actual value except for the effect of the error threshold and thus have been omitted for the sake of conciseness.
Lastly since the present test requires an increment in the vertical position of the robot's attitude stage, thrusters 3 and 12, providing positive vertical motion, remain activated for most of the time window, after a transient phase; these results represent the proof of a successful thruster mapping allowing to convert the triaxial forces and torques into a binary configuration of twelve on/off thrusters.The second example tests the response to a negative vertical translation, starting from different initial conditions, as shown in Table 4, while Table 5 summarizes the weighting matrices used values.
Following the same order of the previous test, Figs.16 and 17 summarize the main results of this second test.
The minimal representation and the comparison with the real robot motion are described in Figs.18 and 19.
Figs. 20-26 provide a more detailed illustration of the complete response, including the actual state components as well as the error components.Similar conclusions can be drawn here as those provided in the previous test.
In this second case the vertical thrusters present an opposite configuration with respect to the previous test: the vertical position, in fact, must decrease, activating, after a transient phase, thrusters 6 and 9 providing negative vertical motion.
Finally, an additional analysis concerning the propellant cost has been conducted, with the aim to demonstrate the discrepancy between the Simulink model and the experiment: the propellant cost term has been computed as a net propellant expenditure as defined in the following equation: On average the propellant cost gap between the simulated solution and the real one corresponds to 3-5% of the real cost.However with equal propellant cost the simulated solution provides a better control, while with equal control the simulated solution requires less propellant.
The consistency between the simulations and the experimental tests could be further increased by incorporating the dual effects of the friction, slowing down the translational motion but at the same time reducing the oscillations, and improving the mass balancing system in order to ensure higher stability in the attitude stage.
The coupling between translation and rotation, by means of the same control inputs, is such that the attitude control is highly affected by the translational one, the latter  requiring much greater thrust values, which tend to generate oscillations or overshoot in the rotational components.It should be emphasized, however, that these irregularities never produced instability and that, as a result of tuning techniques and gain adjustments, these effects have been strongly reduced obtaining satisfactory results.

Conclusion
This work presented a real time, LQR-based sub-optimal approach for full six-degree-of-freedom spacecraft control.Previous works proved the effectiveness and reliability of linear quadratic optimal control, while this research has  illustrated a way to quantify the trade-off between realtime feasibility and optimality.This consisted of a direct comparison with the general purpose optimal control software GPOPS-II.A deepened cost analysis has shown a discrepancy of 15%, although the sub-optimal solution is characterized by easier real-time implementability, reduced computational burden and hence lower simulation time.The proposed control approach has been tested on the ADAMUS test bed with a series of hardware-in-the-loop experiments that constituted the final validation of the  Both the error and the actual state response have been analyzed, together with the thrusters profile, for a complete study of the controller performances.The results showed that both position and attitude reach the desired values, although some differences occur: the Euler angles plots, in fact, lack of the same regularity as the translational components, presenting occasional oscillations or slight overshoots.This imbalance directly derives from the existing differences between the translational dynamics and the rotational one, the former being much slower than the latter, and represents the only limit of the proposed controller, which is the main subject of the present research.

Fig. 1 .
Fig. 1.Vectorial representation of the chaser position with respect to the target, including the inertial frame E, the local vertical local horizontal frame LVLH and the body frame B; attitude representation through the Euler angles parametrization, according to the sequence yxz (box).

Fig. 2 .
Fig. 2. Axonometric view of the attitude stage (AS) and geometric representation of the thrusters configuration: dotted lines represent the inertial principal axes, black thin lines refer to the AS with the relative upper/lower arms, thick arrows illustrate the thrust direction, light thin lines symbolize the imaginary square circumscribed to the AS base.

Fig. 18 .
Fig. 18.Experimental results.Second test: representation characterized by the data illustrated in Table4.

Table 1
Costs comparison.

Table 2
First experimental test: simulation data.

Table 3
First experimental test: weighting matrices data.

Table 4
Second experimental test: simulation data.

Table 5
Second experimental test: weighting matrices data.