Jéromine Dumon
email: jeromine.dumon@isae.fr

Yannick Bury
email: yannick.bury@isae.fr

Nicolas Gourdain
email: nicolas.gourdain@isae.fr

Laurent Michel
email: laurent.michel@isae.fr

NUMERICAL AND EXPERIMENTAL INVESTIGATIONS OF BUFFET ON A DIAMOND AIRFOIL DESIGNED FOR SPACE LAUNCHERS APPLICATIONS

Keywords: buffet, URANS, LES, wind tunnel test, fluid structure interaction, composite

The capability to reuse space launchers for new missions requires to better understand flow phenomena in the transonic regime, such as buffet, and its interaction with the structure. Indeed, the mechanical integrity of the launcher can be compromised by shock/boundary layer interactions, that induce lateral forces responsible for plunging and pitching moments. This work reports some numerical and experimental investigations about the aerodynamic and aero-elastic behavior of a diamond airfoil, designed for microsatellite-dedicated launchers, with a particular interest for the fluid/structure interaction during buffeting. Experiments have been conducted, based on Schlieren visualizations, and compared with numerical predictions obtained with unsteady RANS and Large-Eddy Simulation. Finally, the effect of buffeting on the composite aileron is studied by solving the equation of the dynamics, showing that the aerodynamic response of the airfoil tends to damp the structural displacement, and thus limit the effect of buffeting.

Introduction

Microsatellite-dedicated launchers, for which payload is lower than 50kg, are of paramount importance for future space missions. Among the difficulties encountered for the design of such launchers, the aero-elastic behavior of the ailerons in transonic regimes still remains nowadays partly unknown. Such ailerons are originally designed to provide stability to the launcher, especially when lateral winds are experienced. The PERSEUS' project (French acronym for Academic and Scientific European Student project for Space research) is an initiative of CNES, the French Space Agency, to promote the emergence of innovative technical solutions for space launchers. The present work takes place in the frame of the development of SERA (Supersonic Experimental Rocket ARES) series launchers (Figure 1). To improve the stability of the rocket, SERA is equipped with three ailerons, composed of diamond airfoils made in composite material. With the objective to increase their reliability in turbulent transonic flows, it is necessary to better understand the interaction between the unsteady flow, including shock induced flows, and the composite walls. Moreover, such information is relevant in the context of reusable launchers, where the number of cycles that can be accomplished by the aileron must be accurately predicted. It is thus necessary to predict the loads that occurs during buffeting, to ensure that structure components and subsystems possess adequate strength, stress and fatigue margins in regard to the structural dynamic response. Buffeting is a well-known instability, that occurs in the transonic regime. Buffet is characterized by fluctuating pressures resulting from flow-induced turbulence, flow separation, wake effects, and shock oscillations. The interaction between the shock wave and the separated boundary layer causes the inception of instabilities responsible for a selfsustained periodic motion of the shock over the surface of the airfoil. In a typical flight of a SERA rocket, the transonic regime occurs during less than two seconds during its atmospheric phase as shown in Figure 2. Transonic buffet is observed in many aeronautical applications, including internal flows (e.g. compressor passages) and external flows (e.g. aircraft wings). This phenomenon has been extensively studied, see for example the works of Percey [START_REF] Pearcey | The interaction between local effect at the shock and rear separation -a source of a significant scale effects in wind-tunnel tests on airfoils and wings[END_REF], Tijdeman [START_REF] Tijdeman | Investigations of the transonic flow around oscillating airfoils[END_REF] and Lee [START_REF] Lee | Self-sustained shock oscillations on airfoils at transonic speeds[END_REF], thanks to experimental campaigns on reference geometries [START_REF] Mc Devitt | Static and dynamic pressure measurements on a NACA0012 airfoil at AMES high Reynolds number facility[END_REF][START_REF] Jacquin | Experimental Study of Shock Oscillation over a Transonic Supercritical Profile[END_REF] or numerical simulations [START_REF] Thiery | Numerical Prediction of Shock Induced Oscillations over a 2D Airfoil: Influence of Turbulence Modelling and Test Section Walls[END_REF] [7] [START_REF] Brunet | Zonal-Detached Eddy Simulation of Transonic Buffet on a Civil Aircraft Type Configuration[END_REF]. Based on this extensive knowledge, it has been possible to delay or alleviate buffeting in such geometries [START_REF] Corre | Transonic Flow Control Using a Navier-Stokes Solver and a Multi-Objective Genetic Algorithm[END_REF] [10] [START_REF] Gao | Active control of transonic buffet flow[END_REF]. Unfortunately, the detailed mechanisms that are responsible for the buffet inception and its dynamics are still debated. Moreover, contrary to classical supercritical profiles for civil aircraft, there is a lack of studies for diamond airfoils, adapted to supersonic flows, which are the target of this work.

The first part of this paper deals with the experimental and numerical methods that have been used to study buffet. In a second part, aerodynamic data are compared and analyzed, in order to highlight some of the mechanisms related to buffet for a non-moving airfoil. In the last part of the paper, numerical simulations are conducted considering a moving airfoil, that dynamically responds to aerodynamic forces. Finally, some conclusions and perspectives are drawn.

Methods

2.1.

EXPERIMENTAL SETUP

The ISAE-SUPAERO transonic wind tunnel has a 130mm-by-80-mm rectangular slotted test section, Fig. 3(a). It is powered by four vacuum pumps and provides flow Mach numbers ranging from 0.7±0.05 to 1.3±0.1. Here the Mach number is determined from the stagnation pressure P i and the test section static pressure P according to equation ( 1)), with Ma the Mach number and γ the heat capacity ratio,

" " # = 1 + '() * 𝑀𝑎 * ( - -./ . ( 1
)
Time-resolved Schlieren visualizations are recorded using a high speed Photron camera. Two sets of data can be recorded: 1) 704×512 pixels' image with an acquisition frequency of 20,000 frames per second or 2) a 512×272 pixels' image with an acquisition frequency of 50,000 frames per second. However, Schlieren technique intrinsically integrates 3D information into a 2D image.

This complicates the analysis of the images when the flow naturally exhibits 3D structures, and is responsible for a hard-to-quantify inaccuracy.

(a) (b) In order to determine the potential occurrence of a coupling between the characteristic frequencies of the flow (in particular with the oscillating shock waves) and a specific vibratory frequency of the aileron, a modal analysis is conducted on a vibrating pot (Figure 5). The setup of the aileron on the vibrating pot is chosen similar to its setup in the wind tunnel model, taking into account both the fixing beam of the aileron to the wind tunnel structure and the window in close contact with the aileron. The modal analysis shows three main natural frequencies 𝑓 ) = 196 𝐻𝑧, 𝑓 * = 226 𝐻𝑧 and 𝑓 6 = 850 𝐻𝑧. The frequency f 3 is associated with the fixation of the window to the vibrating pot, f 1 and f 2 correspond to the two first flexion modes of the aileron.

The predicted aerodynamics frequencies are away of the above-mentioned structural frequencies, which ensures that i) the potential occurrence of pressure fluctuations on the surface of the aileron due to the buffeting phenomenon and to other flow unsteadinesses during the wind tunnel tests will not be induced or enhanced by the structural deformation of the aileron, ii) the aileron will not experience severe deformation promoted by the aerodynamic excitation and its coupling with the structural deformation of the model under resonant effects.

NUMERICAL SETUP

This section presents the numerical analysis, led with unsteady Reynolds Averaged Navier Stokes (RANS) and Large Eddy Simulation (LES) approaches, considering the operating conditions of the wind tunnel (including wind tunnel walls). As the zone located downstream of the shock wave is subsonic, it is expected that perturbations generated in the boundary layer and in the close wake of the aileron travel upstream and impact the shock development in the zone of the lambda-shaped shock pedestal, especially as the expected buffeting results from an interaction between separated boundary layer and the shock wave.

The numerical model used in the URANS approach is purely two dimensional, in order to reduce computational time effort, corresponding to the section of the aileron. The dimensions of the computational domain are similar to those of the wind tunnel test section (130mm high and 30 chords long). The center of the model is positioned 10 chords downstream of the inlet of the domain. For the LES computations, the 2D section of the actual wind tunnel model is extruded in the spanwise direction, with a span corresponding to 25% of the chord, in order to ensure uncorrelated turbulence.

The URANS simulations are performed using STAR-CCM+ v11.02. The k-ω SST-Menter turbulence model [START_REF] Menter | Two-Equation Eddy-Viscosity Turbulence Models for Engineering Applications[END_REF] is used for modeling the turbulence. A compressible solver is used with a second order Runge-Kutta scheme for the time discretization. The spatial discretization of the convective fluxes is performed by a third-order MUSCL scheme. Regarding the grid, an unstructured polyhedral 2D mesh is generated with StarCCM+, based on prism layers close to the airfoil walls and polyhedral cells in the rest of the computational domain. The size of the mesh is highly refined close to the aileron surface, and in the zone where the shock waves are expected to develop. The prism layers were set on the aileron and wall surfaces to better capture the flow gradients in the boundary layer. The size of the first layer was chosen to impose y + <0.5 on the aileron surfaces. A grid convergence study was performed, based on the evolution of both drag and lift coefficients, showing that a 1 Million cell mesh is sufficient to ensure convergence on the efforts while capturing the buffeting phenomenon.

As shown in Fig. 7, for URANS, the inlet and outlet flow conditions are modelled as freestream and the walls are considered as adiabatic with a no-slip condition. For LES, total pressure and temperature are imposed at the inlet, static pressure at the outlet, and walls are considered as adiabatic with a no-slip condition.

Figure 7. Schematic of the computational domain, mesh grid refinement as a function of grey levels and boundary conditions

LES computations are performed using the CharLES X solver [START_REF] Bermejo-Moreno | Solving the compressible Navier-Stokes equations on up to 1.97 million cores and 4.1 trillion grid points[END_REF], which solves the spatially filtered compressible Navier-Stokes equations using a finite volume formulation on unstructured hexahedral meshes. A fourth-order central scheme is used for the computation (2 nd order on stretched volumes as in the present study). An explicit third-order Runge-Kutta (RK3) scheme is used for time integration with the Vreman subgrid-scale (SGS) model [START_REF] Vreman | An eddy-viscosity subgrid-scale model for turbulent shear flow: Algebraic theory and applications[END_REF]. The approach relies on the combination of a non-dissipative centered numerical scheme and an essentially nonoscillatory (ENO) second-order shock-capturing scheme, with a shock sensor [START_REF] Bermejo-Moreno | Solving the compressible Navier-Stokes equations on up to 1.97 million cores and 4.1 trillion grid points[END_REF]. Two grid strategies have been considered. The first one relies on a wall-modelling approach, with y + ≈15, x + ≈ 30 and z + ≈ 50, leading to a 30 million cells grid. The second method relies on a wall resolved approach [START_REF] Kawai | Wall modeling in large eddy simulation: length scales, grid resolution and accuracy[END_REF], with y + ≈ 1, x + ≈ 30 and z + ≈ 20, that leads to a 120 million cells grid. Beyond the mesh size reduction, the main interest with the wall-modelling approach is the possibility to increase the time step by a factor of 10 in contrast to the wall resolved approach. Indeed, the cost ratio here between wall resolved and wall modelling approaches is around 40.

Analysis of results at wind tunnel conditions

The spectral analysis presented in section 2 is applied for the numerical data on the density gradient captures. Moreover, aerodynamic loads (drag and lift) are analyzed in order to highlight the link between loads and the position of shock waves, using Fast Fourier Transform (FFT).

Shock wave and separated boundary layer oscillations are observed in both experimental and numerical simulations. The analysis of the Schlieren-based data shows an aperiodic three-dimensional, time-varying deformation of the shock wave in the spanwise direction, Figure 9.

A secondary oscillation of the shock wave is observed in both numerical and experimental data: the local boundary layer separation induces an oblique weak shock wave in front of the strong shock wave, resulting in a classical delta shaped pedestal. While the dominant frequency results from the interaction between the strong shock wave and the boundary layer separation, the secondary oblique shock wave also oscillates due to the local periodic flow separation in front of the strong shock wave.

An important observation is that experimental data exhibit a chaotic behavior, Fig. 10. Two states are alternatively distinguished: the "flag state" and the "delta state". The flag state corresponds to a 3D shock, oscillating in the spanwise direction (as a wind-flapping flag). Due to the 2D spatial integration of the Schlieren technique, it cannot be confirmed if it is or not periodic. The delta state is a state where a sinusoidal oscillation of the shock is observed. In this state, both the oblique weak shock and the strong shock oscillate at the same frequency. For this state, the flow can reasonably be considered as two dimensional.

With an angle of attack (2 degrees), both flag and delta states are present. But contrary to a zero angle of attack, the scenario doesn't switch alternatively from one state to the other state. The flag state is present on the suction side of the aileron and the delta state on the pressure side.

The frequency of the shock oscillation is evaluated to f=620 Hz (URANS) and f=310 Hz (LES) with the spectral analysis. However, the results should be considered with caution for LES due to the limited amount of time available. As mentioned, the analysis of the Schlieren-based data remains difficult due to the presence of a three-dimensional, time-varying deformation of the shock wave in the spanwise direction. The magnitude of the shock wave oscillations varies depending on the method used (LES, URANS or wind tunnel tests). It is worth to mention that the angle of attack in the experimental setup is not exactly zero, as revealed by a non-strictly symmetric distribution of the shock wave footprint. Relatively to the wind tunnel results, the amplitude of the shock wave oscillations is better predicted with LES than with URANS. However, the mean position of the shock is better predicted with URANS than with LES. The shock position is greatly dependent on the boundary layer thickness on the lateral walls for a given Mach number, which is a difficult parameter to predict.

The aerodynamic loads are unsteady due to the shock oscillations, as shown in Figure 13. Indeed, the lift and drag oscillations are driven by the buffet phenomenon. The FFT in Figure 13 highlights a fundamental peak at 516 Hz and its harmonics. The Strouhal number, based on the upstream velocity and the chord, is equal to 0,089.

The analysis of load signals is not straightforward. Interestingly, correlations between instantaneous flow fields and lift signals reveal that the instant at which the lift is null does not correspond to the instant where the shock waves are symmetric on both sides of the airfoil. This behavior is explained by the dynamics of the shock wave, which is different depending on its direction of displacement (upstream or downstream). This shock dynamic is thus associated with both lift and drag signals which are not sinusoidal, which in turn explains why the zero lift is not achieved when the positions of the shocks are symmetric. This analysis also shows that the drag frequency is twice the lift frequency.

The lift and drag predicted with the wall-resolved LES are shown in Figure 14. Only one period of the shock oscillation has been simulated at the moment, which is insufficient to fully analyzed the spectral content of the signal. However, these results highlight two main frequency: one low frequency associated to buffeting and one higher frequency, related to the oblique shock oscillation. The magnitude order in terms of lift and drag are similar to URANS predictions.

Analysis of results at real operating conditions

To complete the numerical simulations, a study has been performed at Reynolds condition close to the one encountered by the SERA rocket at real dimensions. The URANS approach is used with the same methodology, as previously described. The 2D geometry considered here is based on a slice of the real aileron, at mid-span of the real swept aileron. At this location, the chord of the airfoil is 200mm and the relative thickness is 12% of the chord, Fig. 15.

The numerical simulations are performed for several Mach numbers so that the transonic range is entirely covered, from Ma = 0.77 to Ma = 0.83. Figure 15 depicts three instants of the simulation which allows to identify buffeting thanks to the shock wave and the boundary layer oscillations. Buffeting is observed for Mach numbers ranging from 0.78 to 0.81. As illustrated in Figure 16, the normalized frequency of buffeting is constant, corresponding to a Strouhal number (based on upstream velocity and axial chord) of 0.083, which is still very close to the Strouhal number reported for the 1:4 th scaled model-based URANS computations. The amplitude of oscillations increases as a function of the Mach number.

Fluid-structure interaction: impact of buffeting on aileron

The objective is now to study the effect of aerodynamic forces on the dynamic response of the profile. As a first approximation, only rigid movements of the aileron are considered, such as bending and torsion, as shown in Fig. 17(a). The 3D aileron is reduced to its 2D extruded shape. The bending is modeled by a pure vertical translation in the plane while the torsion is modeled by a rotation in the plane, as in Fig. 17 The aileron is made of a sandwich composite material. The core is an epoxy foam and the skin is a carbon laminate. The aileron is considered as a beam with a thin web cross section. Moments of inertia are then calculated geometrically, for a diamond cross section of diagonals equal to 220mm and 12.5 mm, with a 1 mm thick web corresponding to the laminate skin. The Young's The dynamic response of the aileron is driven by the fundamental equation of dynamics. The ideal solution is to simulate the solid displacement in a fully coupled way with the flow. However, the large difference between the characteristic times of the fluid and the solid makes this approach impracticable in the present case. The method relies thus on a separation between the numerical simulation of the flow and the resolution of the aileron dynamics. As a first step, only flapping is considered, that is driven by:

(

) 2
with M is the matrix of mass, D is the matrix of damping, K is the matrix of stiffness and F A are the aerodynamic forces. In order to simplify the problem, an equivalent homogeneous material is considered, so matrices are reduced to a scalar.

This equation is time marched by means of a classical four steps Runge-Kutta scheme. At each time step, the aerodynamic forces are composed of two components: one related to the phenomenon of buffet and one due to the reaction of the aerodynamic force, induced by the profile displacement. The buffet force is extracted from the numerical simulations (URANS and LES database) and the aerodynamic response is modelled. Due to the low thickness of the profile, and the fact that displacement velocity is small compared to the fluid velocity, the lift force is modelled using the thin profile theory, with a compressibility correction [START_REF] Anderson | Fundamentals of Aerodynamics[END_REF]:

(3)

The angle of the flow seen by the profile at each instant of time is estimated as

(4)
Since the flow is responding in phase with the solid, the aerodynamic response is a positive damping term in the equation of the dynamic, Eq. ( 2). The main limitation of this approach is that the aerodynamic response to the profile displacement is instantaneous, while a lag is observed in practice. For this reason, this method is valid only when the ratio between the flow velocity and the displacement velocity of the solid is large (so the time lag becomes negligible). The results are expressed in terms of non-dimensional parameters U * , representing the ratio between the fluid velocity and the displacement velocity of the solid and m * , representing the mass ratio between the solid and the fluid:

U ⇤ = U 1 1 2⇡ q k m .C F lif t = 2⇡(↵ ↵ 0 ) p 1 M 2 ⇥ 1 2 ⇢ 1 U 2 1 .S M ẍ + D ẋ + Kx = X F A tan(↵) = ẋ U 1 m ⇤ = m ⇢SL
with k the stiffness, C the chord of the aileron and m the mass of the solid. In the case of the aileron at real dimensions, the value of the normalized velocity U * is estimated to be close to 100. Two cases are considered: without aerodynamic coupling (only the force coming from buffet is applied) and with coupling (aerodynamic response of the profile is added). The evolution of the normalized displacement, y * =y/C, with respect to the normalized velocity U * is represented in Fig. 18(a) (using URANS data) and Fig. In Fig. 18(a), the peak close to U * =10 corresponds to the resonance between the natural frequency of the solid and the aerodynamic excitation (buffet). In the case of LES, the buffet is associated to a more complex signal in term of spectral content, so two frequencies of resonance are shown in Fig. 18(b). In the uncoupled case, the maximum displacement increased with U * (corresponding to a reduction of the stiffness). In the coupled case, as expected, the aerodynamic forces act as a damping term, which reduced the amplitude of the displacement, especially at the resonance frequency. A plateau is also reached for value of U * higher than 500. At such velocity ratios, the periodic excitation due to buffet is no longer seen by the solid (since the time-average force of buffet is null, the solid does not react to this phenomenon). The effect of buffet predicted by LES has a weaker effect on the solid displacement compared to URANS.

Two normalized masses are compared: m * =80 (light hollow aileron) and m * =1600 (heavy plain aileron). The previous conclusions are globally unchanged with the lower mass ratio parameter (m * =80), as shown in Fig. 19. However, as expected, the maximum displacement is increased (especially at low to moderate values of U * ) and a smoothening of the resonance peak, due to an increase of the aerodynamic damping (velocity displacement of the solid is increased compared to m * =1600, and so the value of the flow angle is also increased). In all cases, the value of the angle stays below 7-8 degrees, which is still acceptable for the thin profile theory.

Conclusions

Both experimental and numerical investigations have been performed to understand the buffeting phenomenon and its impact on a composite aileron. This study can be summarized by the following points:

-First, buffeting appears at transonic speed on a diamond aileron at zero angle of attack, as verified with both numerical simulation and measurements, -Experimental approach reveals a 3D oscillation of the shock, but at this step there is no clear evidence of the physical phenomena promoting this 3D mode, -LES identifies an oscillation of the width of the shock foot (λ-width) which has a significant impact on the aerodynamics load, but this phenomenon is not seen with the URANS simulation, -at last, the experiments highlighted a chaotic behavior of the shock wave that shifts from an oscillatory periodic state (delta state) to an erratic 3D flapping state (flag state), transiently separated by transition states, where the shock moves from a close-to-trailing edge position (delta state) to a close-to-dihedron position (flag state).

The results reported in this paper still highlight a lack of deciphering of the buffeting origin. A 3D numerical simulation (both URANS and LES) should be conducted in order to better understand the spanwise oscillation of the shock observed experimentally. More information from the wind tunnel is also needed. The use of unsteady pressure sensors on the aileron will provide a more accurate validation of the numerical simulations in the future. Fluid/Structure interaction has been investigated, considering only the bending mode. The effect of torsion should now be investigated since it should have more impact on the flow due to the increase of incidence that amplifies buffeting effects. A provision for future work will also consist in using a time-dependent stiffness, in order to represent the influence of progressive damages of unsteady loads on the aileron (fatigue mode).

Figure 1 .

 1 Figure 1. Global view of the SERA launcher and details of the aileron shape and dimensions [mm]

Figure 2 .

 2 Figure 2. Evolution of the Mach number, as recorded during a typical flight of the SERA rocket

Figure 3 .

 3 Figure 3. (a) View of the wind tunnel equipped with Schlieren test bench and (b) detailed view of the aileron in the test section

Figure 4 .

 4 Figure 4. Schematic view and dimensions of the wind tunnel aileron

Figure 5 .Figure 6 .

 56 Figure 5. (a) Close view of the vibrating pot and (b) details of the laser scanner used for the modal analysisOn the basis of the time-resolved Schlieren visualizations (Figure6), a spectral analysis of the shock oscillation is proposed, based on a three-step process:-A one pixel-height sensor line is selected in the shock oscillation area, -A time series of the grey level signal is then extracted, -The Power Spectral Density (PSD) of this signal is computed (if the signal is periodic).

Figure 8 .

 8 Figure 8. Schlieren visualizations at 2 distinct instants revealing the oscillations of the shock associated with the buffeting: (a-b) measurements, (c-d) URANS and (e) LES.

Figure 9 .Figure 10 .

 910 Figure 9. 3D supposed shock wave visualization on the aileron confronted to the 2D Schlieren vizualisation

Figure 11 .

 11 Figure 11. Temporal history of the shock state changes (4 corresponds to the delta state and 1 to the flag state)

Figure 12 .

 12 Figure 12. Shock magnitude oscillation for (a) LES, (b) URANS and (c) wind tunnel test

Figure 13 .

 13 Figure 13. aerodynamics efforts (left: lift, right: drag) as a function of time [s] (URANS results) (left) and FFT of lift (up right) and drag (down right) from the URANS simulation

Figure 14 .Figure 15 .

 1415 Figure 14. Aerodynamics loads (left: lift, right: drag) as a function of time [s] (LES results)

Figure 17 .

 17 Figure 17. Example of the (a) aileron deformation (bending and torsion) and (b) 2D modelling

Figure 16 .

 16 Figure 16. Lift coefficient and Strouhal number evolution with respect to the Mach number

Figure 18 .

 18 Figure 18. Evolution of the maximal displacement and maximal velocity displacement, with respect to the parameter U * : (a) URANS and (b) LES (m * =1600)

Figure 19 .

 19 Figure 19. Influence of the mass ratio parameter on the displacement and displacement velocity at U * =100 (URANS data)

Table 1 :

 1 modulus E and the shear modulus G are material properties, determined by the Classical Laminate Theory[START_REF] Berthelot | Matériaux composites, « comportement mécanique et analyse des structures[END_REF]. The stiffness results are presented in Table1. Estimation of the stiffness of the real aileron

		𝐾	𝑓 [𝐻𝑧]
	Flexion	1373 N/m	13,2
	Torsion	346 Nm	55,1

Acknowledgements

This work is partially funded by the French Space Agency (CNES), this support is greatly acknowledged. The authors particularly thank Jean Oswald and Florent Puel from DLA (Direction des LAnceurs) at CNES for their help on this study. The experimental campaign was conducted with the technical team of ISAE-Supaero (special thanks to Emmanuel Rivet, Patrick Cheze and Marc Grellet). The simulations performed were achieved using CALMIP computing means, under project p1425 and GENCI A0042A07178.


Julio C Sanchez

Francisco Gavilan

Rafael Vazquez
email: rvazquez1@us.es

Christophe Louembet

A Flatness-Based Predictive Controller for Six-Degrees of Freedom Spacecraft Rendezvous

Keywords: Impulsive rendezvous, Attitude control, Model predictive control, Flatness theory Nomenclature B Input matrix e Rotation axis H Angular momentum Id Identity matrix

This work presents a closed-loop guidance algorithm for six-degrees of freedom spacecraft rendezvous with a passive target flying in an eccentric orbit. The main assumption is that the chaser vehicle has an attitude control system, based on reaction wheels, providing the necessary torque to change its orientation whereas the number of thrusters is arbitrary. The goal is to design fuel optimal manoeuvres while satisfying operational constraints and rejecting disturbances. The proposed method is as follows; first, the coupled translational and angular dynamics are transformed to equivalent algebraic relations using the relative translational states transition matrix and the attitude flatness property. Then, a direct transcription method, based on B-splines parameterization and discretization of time continuous constraints, is developed to obtain a tractable static program. Finally, a Model Predictive Controller, based on linearization around the previously computed solution, is considered to handle disturbances. Numerical results are shown and discussed.

Autonomous spacecraft rendezvous and docking is becoming a more important topic in the space industry as access to space continues increasing. From the first rendezvous attempts (Gemini missions) to the Rosetta mission in 2014, the rendezvous manoeuvre has played a key role in different kinds of space missions such as Apollo, ISS, Hubble, etc. After decades of development, many approaches to achieve rendezvous for different mission profiles have been used, see [START_REF] Woffinden | Navigating the Road to Autonomous Orbital Rendezvous[END_REF] for an historical review or [START_REF] Fehse | Automated Rendezvous and Docking of Spacecraft[END_REF] for the basics. Nowadays, an increasing interest to demonstrate autonomous rendezvous and flight formation operations for lightweight and low-power spacecraft is arising with CPOD, PRISMA and PROBA-3 missions as examples, see [START_REF] Roscoe | Overview and GNC design of the Cubesat Proximity Operations Demonstration (CPOD) Mission[END_REF][START_REF] Persson | PRISMA -An Autonomous Formation Flying Mission[END_REF][START_REF] Castellani | PROBA-3 mission[END_REF].

Typically, the rendezvous problem has been widely studied just considering orbit control making the assumption that translational and rotational motions are decoupled. This problem has been usually tackled by means of direct transcription methods which transform the optimal control problem into a discrete optimization problem as in [START_REF] Lu | Autonomous Trajectory Planning for Rendezvous Proximity Operations by Conic Optimization[END_REF][START_REF] Breger | Safe Trajectories for Autonomous Rendezvous of Spacecraft[END_REF][START_REF] Richards | Spacecraft Trajectory Planning with Avoidance Constraints Using Mixed-Integer Linear Programming[END_REF][START_REF] Gavilan | Chance-constrained model predictive control for spacecraft rendezvous with disturbance estimation[END_REF][START_REF] Vazquez | Pulse-width predictive control for LTV systems with application to spacecraft rendezvous[END_REF] among others. The main advantage of these methods, against indirect ones, is that several kinds of constraints can be easily added to the problem such as approach corridors through the docking axis (V-bar or R-bar guidance), way-points, thrust direction inhibition, obstacle avoidance or fault-tolerant trajectories.

However, orbit and attitude control subsystems are mutually coupled, which is mainly due to the dependence of the thrusters orientation on the relative attitude between target and pursuer (at least in the short-term). Spacecraft attitude planning for direction reorientation manoeuvres, which are the ones needed to point the thrusters in an adequate way, is a topic with a vast literature. Reference [START_REF] Leve | Spacecraft Momentum Control Systems[END_REF] proposed two dimensional attitude profiles with time derivatives saturation up to the jerk. Model Predictive Control (MPC) techniques based on linearization around a set point have been used in the works of [START_REF] Øyvind Hegrenaes | Spacecraft Attitude Control using Explicit Model Predictive Control[END_REF][START_REF] Guiggiani | Fixed-Point Constrained Model Predictive Control of Spacecraft Attitude[END_REF]. A remarkable approach is the one followed by [START_REF] Louembet | Path planning for satellite slew manoeuvres: a combined flatness and collocationbased approach[END_REF][START_REF] Caubet | A Motion Planning Method for Spacecraft Attitude Maneuvers Using Single Polynomials[END_REF] which is based on the attitude dynamics flatness property (see [START_REF] Fliess | Flatness and defect of non linear systems: introductory theory and examples[END_REF] for more details about flatness theory) that allowed them to transform the attitude dynamics into algebraic relations avoiding the need of numerical integration.

Regarding previous works on six-degrees of freedom relative motion, adaptive tracking controllers based on feedback has been considered by [START_REF] Filipe | Adaptive Position and Attitude-Tracking Controller for Satellite Proximity Operations Using Dual Quaternions[END_REF] for rendezvous and by [START_REF] Wong | Output Feedback Control for Spacecraft Formation Flying with Coupled Translation and Attitude Dynamics[END_REF][START_REF] Wang | Sun, 6-DOF robust adaptive terminal sliding mode control for spacecraft formation flying[END_REF] for flight-formation while strategies based on backstepping control have been employed by [START_REF] Kristiansen | Spacecraft coordination control in 6DOF: Integrator backstepping vs passivity-based control[END_REF] for flight-formation and by [START_REF] Zhang | Integrated relative position and attitude control of spacecraft in proximity operation missions with control saturation[END_REF][START_REF] Yan | Integrated Translational and Rotational Control for Rendezvous and Docking on Ellipse Orbits[END_REF] for rendezvous operations. Sliding mode control has also been explored by [START_REF] Terui | Position and Attitude Control of a Spacecraft by Sliding Mode Control[END_REF]. References [START_REF] Naasz | Integrated Orbit and Attitude Control for a Nanosatellite with Power Constraints[END_REF][START_REF] Siva | Coordinated 6-dof control of dual spacecraft formation[END_REF][START_REF] Moon | Quaternion Based Attitude Control and Suboptimal Rendezvous Guidance on Satellite Proximity Operation[END_REF][START_REF] Wu | Relative Motion Decoupled Control for Spacecraft Formation Translational and Rotational Dynamics[END_REF] proposed a two stages approach, first they used an optimal control method for the translational motion, LQR in [START_REF] Naasz | Integrated Orbit and Attitude Control for a Nanosatellite with Power Constraints[END_REF][START_REF] Siva | Coordinated 6-dof control of dual spacecraft formation[END_REF][START_REF] Moon | Quaternion Based Attitude Control and Suboptimal Rendezvous Guidance on Satellite Proximity Operation[END_REF] or convex optimization in [START_REF] Wu | Relative Motion Decoupled Control for Spacecraft Formation Translational and Rotational Dynamics[END_REF], and then they designed an attitude controller to obtain the orientations demanded by the translational plan. Reference [START_REF] Biggs | Motion planning on a class of 6-D Lie groups via a covering map[END_REF] proposed a covering map to convert the original differential equations into two coupled equations evolving on a 3-D Lie group, this map is applicable to 6-DOF rendezvous. These authors proposed in [START_REF] Sanchez | A Flatness-Based Trajectory Planning Algorithm for Rendezvous of Single-Thruster Spacecraft[END_REF] a method based on the translational state transition matrix and the attitude flatness property to solve the integrated optimal control rendezvous problem. Amongst the previous works, dual quaternions, which contain information of both translational and rotational states, were used in [START_REF] Filipe | Adaptive Position and Attitude-Tracking Controller for Satellite Proximity Operations Using Dual Quaternions[END_REF][START_REF] Wang | Sun, 6-DOF robust adaptive terminal sliding mode control for spacecraft formation flying[END_REF]. Concerning the number of thrusters, [START_REF] Filipe | Adaptive Position and Attitude-Tracking Controller for Satellite Proximity Operations Using Dual Quaternions[END_REF][START_REF] Kristiansen | Spacecraft coordination control in 6DOF: Integrator backstepping vs passivity-based control[END_REF][START_REF] Wang | Sun, 6-DOF robust adaptive terminal sliding mode control for spacecraft formation flying[END_REF] assumed a pair of them available on each direction, six in total, whereas [START_REF] Naasz | Integrated Orbit and Attitude Control for a Nanosatellite with Power Constraints[END_REF] considered four thrusters and [START_REF] Siva | Coordinated 6-dof control of dual spacecraft formation[END_REF][START_REF] Wu | Relative Motion Decoupled Control for Spacecraft Formation Translational and Rotational Dynamics[END_REF][START_REF] Moon | Quaternion Based Attitude Control and Suboptimal Rendezvous Guidance on Satellite Proximity Operation[END_REF][START_REF] Yan | Integrated Translational and Rotational Control for Rendezvous and Docking on Ellipse Orbits[END_REF][START_REF] Sanchez | A Flatness-Based Trajectory Planning Algorithm for Rendezvous of Single-Thruster Spacecraft[END_REF] studied the case of single-thruster operations. The results of [START_REF] Biggs | Motion planning on a class of 6-D Lie groups via a covering map[END_REF] are applied to both a classical six thrusters configuration and a single thruster one. The previous works assumed that torque is provided by an independent ACS system whereas [START_REF] Zhang | Integrated relative position and attitude control of spacecraft in proximity operation missions with control saturation[END_REF] considered six thrusters in a cuboid layout configuration providing both force and torque. Apart from rendezvous and flight-formation operations, coupled motion has also been studied for geostationary satellites station-keeping [START_REF] Weiss | Model Predictive Control for Spacecraft Rendezvous and Docking Strategies for Handling Constraints and Case Studies[END_REF] and solar sails control [START_REF] Gong | Coupled attitude-orbit dynamics and control for displaced solar orbits[END_REF].

In this paper, we consider a spacecraft equipped with reaction wheels and an arbitrary number of thrusters which seeks to rendezvous with a target flying in an eccentric orbit. The employed formulation allows to consider the coupled problem in an optimal way without any assumptions on the number of available thrusters which increase the applicability of the algorithm to different types of missions. In a similar way as [START_REF] Urbina | Guidance and robust control methods for the approach phase between two orbital vehicles with coupling between translational and rotational motions[END_REF], a hybrid system is considered where the propulsive action is modelled as impulses but the attitude control is time continuous.

The proposed solution method transforms the time-continuous dynamics into algebraic relations by means of the translational state transition matrix and the attitude flatness property. Then, this equivalent optimal control problem is parameterized and discretized to obtain a finite tractable static program. Once an open-loop solution is obtained, a closed-loop MPC scheme, see [START_REF] Camacho | Model Predictive Control[END_REF], based on linearization around the previously computed solution, is developed to reject disturbances and cope with unmodelled dynamics.

The structure of this paper is as follows. Section 2 describes the coupled translational and angular motion for spacecraft rendezvous. Next, Section 3 presents the time-continuous rendezvous problem and its conversion to an equivalent problem. Section 4 describes the employed methodology to solve this equivalent problem by means of parameterization and discretization. Section 5 presents the linearized close-loop MPC scheme. Section 6 shows results for cases of interest. Finally, Section 7 closes this paper with some additional considerations.

Model of Spacecraft Rendezvous

In this section, a six-degrees of freedom model for spacecraft rendezvous is presented. Firstly, the translational relative motion between the two vehicles is derived; secondly, the chaser angular motion model is described; and finally, both translational and angular motions are coupled.

Translational motion

There is a considerable number of translational dynamic models for spacecraft rendezvous; the one to be chosen depends on the objectives and constraints on the mission. For instance, if the target vehicle is orbiting in a closed Keplerian orbit, the linearised equations of the relative position between an active chaser spacecraft and a passive target vehicle can be expressed in a cartesian reference frame as in [START_REF] Tschauner | Rendezvous zu einemin elliptischer Bahn umlaufenden Ziel[END_REF], leading to the well known Tschauner-Hempel equations, or by means of its relative orbital elements as in [START_REF] D'amico | Proximity Operations of Formation-Flying Spacecraft Using an Eccentricity/Inclination Vector Separation[END_REF]. In this work, a cartesian reference frame is used

ẍ = νz + 2 ν ż + ν2 x - µx r 3 t + n T p=1 F x,p m , (1) ÿ
= - µy r 3 t + n T p=1 F y,p m , (2)
z = -νx -2 ν ẋ + ν2 z + 2µz r 3 t + n T p=1 F z,p m , (3)
where x, y and z denote the position of the chaser in a local-vertical/localhorizontal (LVLH) frame of reference fixed on the center of gravity of the target vehicle (see Fig. 1), in which z refers to the radial position (positive pointing towards the centre of the Earth), y to the cross-tack position (opposite to the orbit angular momentum) and x closes the right-handed system (note that x is not necessarily aligned with the target velocity due to eccentricity). The velocity of the chaser in the LVLH frame is given by ẋ, ẏ and ż; the variables F x,p , F y,p and F z,p are the projections on the LVLH frame of the thrust force exerted by each one of the n T thrusters; and m is the spacecraft mass which, for close enough rendezvous operations, is considered constant. The variables r t and ν are the target radius and true anomaly along its orbit, which are a function of time and its orbital elements (semi-major axis and eccentricity).

The gravitation parameter of the Earth is µ=398600.4 km 3 /s 2 . The state space model, x(t)=[x, y, z, ẋ, ẏ, ż] T , is governed by the linear time varying system (LTV) given by Eq.( 1)-(3). The independent variable of this LTV system can be changed from time to true anomaly leading to the simplified Tschauner-Hempel equations, see [START_REF] Tschauner | Rendezvous zu einemin elliptischer Bahn umlaufenden Ziel[END_REF]. A formal solution of the Tschauner-Hempel model by means of its state transition matrix, known as the Yamanaka-Ankersen matrix, was proposed by [START_REF] Yamanaka | New State Transition Matrix for Relative Motion on an Arbitrary Elliptical Orbit[END_REF]. This transition matrix is computed by means of its fundamental matrix and inverse without need of numerical integration. In this work, following [START_REF] Vazquez | Pulse-width predictive control for LTV systems with application to spacecraft rendezvous[END_REF] (note that the axes are not the same as in this work), the Yamanaka-Ankersen state transition matrix is expressed by means of the eccentric anomaly E,

Φ Φ Φ(t, t 0 ) = Y E(t) Y -1 E(t0) . (4)
Note that a one-to-one relation exist between time and eccentric anomaly through the Kepler equation

n(t -t p ) = E -e sin E, (5)
where t p is the time at periapsis and is used as a reference point to measure E. The time t p is chosen such that it is equal or less than the starting manoeuvre time denoted by t 0 (substracting, if necessary, any number of orbital periods). Kepler's equation ( 5) is not analytically invertible, but its inverse can be found numerically with any desired degree of precision (see any Orbital Mechanics reference, such as [START_REF] Wie | Space Vehicle Dynamics and Control[END_REF]).

Using the aforementioned state transition matrix, Eq.( 4), and considering, as a simplification, an impulsive model, Eq.( 7), the translational states transition equation is given by

x(t) = Φ Φ Φ(t, t 0 )x(t 0 ) + k i=0 Φ Φ Φ(t, t i )Bu(t i ), t k ≤ t < t k+1 , (6)
where the input matrix is B=[Θ 3×3 , Id] T , and the propulsive control signal u(t) is modelled as impulses (i.e. instantaneous changes of velocity) which describe with adequate accuracy a typical chemical thruster

u(t) = n T p=1 Np k=0 lim ∆t k →0 t k +∆t k t k F p (t) m dt = n T p=1 Np k=0 ∆V p (t k )δ(t -t k ), (7)
being ∆V p ∈R 3 the velocity increment given by the thruster p and N p +1∈N the number of thruster firings during the manoeuvre.

Angular motion

In this section, the attitude representation parameter is chosen and some of their properties are presented. Then, the angular dynamics of a spacecraft considering only internal torques, which are the ones produced by reaction wheels, is derived. Finally, the attitude flatness property of the resulting angular dynamics is introduced (this property will be then exploited in Section 3).

Attitude representation and angular dynamics

In this work, the modified Rodrigues parameters (MRP) representation (see [START_REF] Marandi | A preferred coordinate system and the associated orientation representation in attitude dynamics[END_REF][START_REF] Schaub | Stereographic Orientation Parameters for Attitude Dynamics: A Generalization of the Rodrigues Parameters[END_REF] for more details about MRP) is chosen rather than the widely used attitude quaternion. The modified Rodrigues parameters have the advantage of being a minimal attitude representation and are easier to linearize than attitude quaternions (incremental addition does not work for quaternions). Moreover, the unit-norm constraint of attitude quaternions is avoided in the problem formulation. The counterpart is that MRP suffer singularities when representing 3D rotations. The MRP are denoted as σ σ σ=[σ 1 , σ 2 , σ 3 ] T and its relation with the rotation angle, θ rot , and axis, e, is

σ σ σ = etan(θ rot /4), (8)
where singularities arise when θ rot =±2π(2j-1)π with j∈N. However, they can be avoided by constraining θ rot ∈(-2π, 2π). The rotation matrix to change a vector from one reference frame to another is given by

R(σ σ σ) = Id + 8σ σ σ × σ σ σ × -4(1 -σ σ σ 2 2 )σ σ σ × 1 + σ σ σ 2 2 2 , (9)
where σ σ σ × ∈R 3×3 is the cross product matrix, see [START_REF] Wie | Space Vehicle Dynamics and Control[END_REF]. The attitude evolution of the chaser is defined by the kinematic and dynamic equations. The translational equations are expressed on a local frame so it is of interest to work with the attitude of the body frame with respect to the LVLH frame as in [START_REF] Wie | Space Vehicle Dynamics and Control[END_REF]. The kinematics are given by

σ σ σ(t) = C(σ σ σ(t))[ω ω ω(t) -R(σ σ σ(t))ω ω ω L/I (t)], (10)
being ω ω ω=[ω 1 , ω 2 , ω 3 ] T the angular velocity of the chaser body frame with respect to the inertial frame and ω ω ω L/I =[0, -ν, 0] T the angular velocity of the LVLH frame with respect to the inertial frame expressed on the local frame. The matrix C has the following expression

C(σ σ σ) = 1 4   1 + σ 2 1 -σ 2 2 -σ 2 3 2(σ 1 σ 2 -σ 3 ) 2(σ 1 σ 3 + σ 2 ) 2(σ 1 σ 2 + σ 3 ) 1 -σ 2 1 + σ 2 2 -σ 2 3 2(σ 2 σ 3 -σ 1 ) 2(σ 1 σ 3 -σ 2 ) 2(σ 2 σ 3 + σ 1 ) 1 -σ 2 1 -σ 2 2 + σ 2 3   . ( 11
)
Additionally, the following equation describes the angular momentum variation, expressed on the chaser body frame, when the only considered torques are internal to the system (the ACS consists of reaction wheels)

I ω ω ω(t) + Ḣrw (t) + ω ω ω(t) × H tot = 0, (12)
where I∈R 3×3 and H tot ∈R 3 are, respectively, the moment inertia matrix and the angular momentum of the spacecraft whereas H rw ∈R 3 is the angular momentum of the reaction wheels. Note that Eq.( 10)- [START_REF] Øyvind Hegrenaes | Spacecraft Attitude Control using Explicit Model Predictive Control[END_REF] give the attitude evolution of the body frame, B, with respect to the inertial frame, I. From the fact that no external torques are applied, the spacecraft angular momentum is constant

H tot = H b (t) + H rw (t) ≡ constant, (13)
where H b (t)=Iω ω ω(t) is the angular momentum of the platform. The attitude control signal is the exerted torque by the reaction wheels through its angular momentum variation, Ḣrw (t).

Attitude flatness property

The angular motion given by Eq.( 10) and Eq.( 12) is non-linear, hence accounting for them in the resolution of an optimal control problem usually require numerical integration, see [START_REF] Desai | Six-Degree-of-Freedom Trajectory Optimization Using a Two-Timescale Collocation Architecture[END_REF]. However, the considered angular dynamics has the flatness property and it is called a flat system, see [START_REF] Louembet | Path planning for satellite slew manoeuvres: a combined flatness and collocationbased approach[END_REF].

Remark 1: a flat system has a flat output which can be used to explicitly express all states and inputs in terms of the flat output and a finite number of its derivatives, see [START_REF] Fliess | Flatness and defect of non linear systems: introductory theory and examples[END_REF].

Following [START_REF] Louembet | Path planning for satellite slew manoeuvres: a combined flatness and collocationbased approach[END_REF], the attitude representation parameter σ σ σ(t) is chosen as the flat output. The differential equations of the angular motion, Eq.( 10) and Eq.( 12), can be transformed into algebraic relations, as a function of the flat output and its derivatives. Solving the angular velocity in Eq.( 10) and deriving the obtained expression with respect to time

ω ω ω(t) = C -1 (σ σ σ) σ σ σ + R(σ σ σ)ω ω ω L/I , (14)
ω ω ω(t) = C -1 (σ σ σ) σ σ σ + Ċ(σ σ σ)R(σ σ σ)ω ω ω L/I + Ṙ(σ σ σ)ω ω ω L/I + R(σ σ σ) ω ω ω L/I -Ċ(σ σ σ)ω ω ω, (15)
and introducing Eq.( 14)-( 15) into Eq.( 12), the angular momentum variation of the reaction wheels is explicitly obtained as

Ḣrw (t) = -I C -1 (σ σ σ) σ σ σ + Ċ(σ σ σ)R(σ σ σ)ω ω ω L/I + Ṙ(σ σ σ)ω ω ω L/I +R(σ σ σ) ω ω ω L/I -Ċ(σ σ σ)ω ω ω -C -1 (σ σ σ) σ σ σ + R(σ σ σ)ω ω ω L/I × H tot . (16)
Using the angular momentum conservation, Eq.( 13), the reaction wheels angular momentum can also be expressed as a function of the flat output and its derivatives

H rw (t) = H tot -IC -1 (σ σ σ) σ σ σ, (17)
Note that time dependencies have been omitted at the right-hand side of Eq.( 14)-( 17) for clarity.

Coupling between translational and angular motion

Now, the translational and angular motion coupling between the previous models is presented. The velocity increment given by each thruster p on the LVLH frame, denoted by L, is

∆V p (t k ) = R T (σ σ σ(t k ))w p u p (t k ), u p (t k ) ≥ 0, ( 18
)
where w p ∈R 3 is a unit-norm vector representing the p thruster orientation on the pursuer body frame and u p (t k )∈R is the impulse amplitude of the thruster p at time t k and R(σ σ σ)∈R 3×3 is the rotation matrix between the chaser body frame and the LVLH frame. Introducing Eq.( 18) into Eq.( 7)

u(t) = n T p=1 Np k=0 R T (σ σ σ(t k ))w p u p (t k )δ(t -t k ). ( 19
)
The coupling between translational and angular motion arises when the translational control input given by Eq.( 19) is introduced into the translational states transition equation given by Eq.( 6) leading to

x(t) = Φ Φ Φ(t, t 0 )x(t 0 ) + k i=0 n T p=1 Φ Φ Φ(t, t i )BR T (σ σ σ(t i ))w p u p (t i ), t k ≤ t < t k+1 , (20)
Note that the propulsive action projected on the LVLH frame, u(t), depends on the vehicle attitude in a non-linear way by means of the rotation matrix between the pursuer body frame and the LVLH frame, see Eq.( 9). The angular motion is not affected by the translational motion (gravity-gradient effects are neglected), hence Eq.( 10) and Eq.( 12) still hold for the coupled model.

Rendezvous planning problem

In this section, the objective function and constraints are presented. In a generic form, the rendezvous optimal control problem states as follows minimize

up(t k ), Ḣrw(t) J(u p (t k ), Ḣrw (t)), subject to v(t) = -2ω ω ω × v -ω ω ω × r -ω ω ω × (ω ω ω × r) -µ/ r t 3 2 r -3 r T t r/ r t 2 2 r t , ṙ(t) = v(t) + Np k=0 n T p=1 R T (σ σ σ(t))w p u p (t k )δ(t -t k ), ω ω ω(t) = -I -1 Ḣrw (t) + ω ω ω(t) × H tot , σ σ σ(t) = C(σ σ σ(t))[ω ω ω(t) -R(σ σ σ(t))ω ω ω L/I (t)], g 1 (r(t)) ≤ 0, LOS region g 2 (u p (t k )) ≤ 0, ∆V bounds g 3 (H rw (t), Ḣrw (t)) ≤ 0, ACS bounds g 4 (r(t 0 ), v(t 0 ), σ σ σ(t 0 ), ω ω ω(t 0 )) = 0, t 0 conditions g 5 (r(t f ), v(t f ), σ σ σ(t f ), ω ω ω(t f )) = 0, t f conditions ( 21
)
where

r(t)=[x(t), y(t), z(t)] T , v(t)=[ ẋ(t), ẏ(t), ż(t)] T and r t (t)=[0, 0, -r t (t)] T .
Note that time dependencies have been omitted at the right hand side of the translational dynamics equation for clarity. The control inputs are the thrusters impulses amplitudes at the firing times (which are known beforehand) and the reaction wheels angular momentum variation. Next the objective function and constraints appearing in [START_REF] Zhang | Integrated relative position and attitude control of spacecraft in proximity operation missions with control saturation[END_REF] are detailed in (3.1) and (3.2) respectively.

Objective function

The chosen objective function seeks to minimize fuel consumption, which is equivalent to minimize the L 1 -norm of the applied impulses

J = Np k=0 n T p=1 u p (t k ). ( 22
)
Note that the absolute value symbol is not needed since u p (t k ) are always positive because of Eq.( 18). Moreover, reaction wheels use electrical power and therefore their associated cost do not appear in Eq.( 22).

Constraints of the problem

Three sets of constraints are considered in this paper. Firstly, path constraints on the relative translational states (g 1 ); secondly, the control variables (impulses amplitudes and reaction wheels angular momentum) are bounded (g 2 and g 3 ); and finally, initial and terminal states values are prescribed (g 4 and g 5 ).

Path constraints

For sensing purposes (see [START_REF] Breger | Safe Trajectories for Autonomous Rendezvous of Spacecraft[END_REF]), it is required that the chaser vehicle remains inside a line of sight (LOS) area from the docking port, thus guaranteeing that the chaser spacecraft is at all time visible from the docking port. The LOS region can be defined by the equations x≥c y (y -y 0 ), x≥-c y (y + y 0 ), x≥c z (z -z 0 ), x≥-c z (z + z 0 ) and x≥0; these equations limit the relative translational state space by five planes as shown in Fig. 2. One can define the LOS constraint algebraically, at any instant t, as A L x(t) ≤ b L , where

A L =       -1 c y 0 0 0 0 -1 -c y 0 0 0 0 -1 0 c z 0 0 0 -1 0 -c z 0 0 0 -1 0 0 0 0 0       , b L =       c y y 0 c y y 0 c z z 0 c z z 0 0       . ( 23
)

Control bounds

Regarding the thrusters performance, it is assumed that the impulse amplitude provided by each thruster is bounded above (and below by zero)

0 ≤ u p (t k ) ≤ u p,max , p = 1 . . . n T . (24)
Note that u p can take any value in the allowed interval (it is assumed that thrusters opening times can be adjusted to produce the exact impulse amount). On the other hand, each one of the reaction wheels saturates when it stores a certain amount of angular momentum which is equivalent to consider that reaction wheels velocities are limited. Moreover limits on the angular acceleration exist for each wheel

-H i,max ≤ H i,rw (t) ≤ H i,max , i = 1, 2, 3, (25)
-Ḣi,max ≤ Ḣi,rw (t) ≤ Ḣi,max , i = 1, 2, 3.

(26)

Boundary constraints

The chaser is assumed to depart from a given point and velocity with a given orientation and angular velocity

x(t 0 ) = x 0 , σ σ σ(t 0 ) = σ σ σ 0 , ω ω ω(t 0 ) = ω ω ω 0 , (27)
and it has to met prescribed states at the end of the manoeuvre

x(t f ) = x f , σ σ σ(t f ) = σ σ σ f , ω ω ω(t f ) = ω ω ω f , (28)
where the last three components of x f are null and ω ω ω f =R(σ σ σ f )ω ω ω L/I (t f ) to have no relative angular velocity between the body and LVLH frame.

Equivalent rendezvous planning problem

The aim of the optimal control problem ( 21) is to guarantee rendezvous with the target along a prescribed approach region (LOS) while respecting control bounds and minimizing fuel consumption. Using the coupled transition equation for the translational states, see Eq. [START_REF] Kristiansen | Spacecraft coordination control in 6DOF: Integrator backstepping vs passivity-based control[END_REF], and the algebraic relations derived from the attitude flatness property, developed through Eq.( 14)-( 17), we formulate an equivalent planning problem, where differential equations are replaced by algebraic ones (as a function of the flat output and its derivatives), without losing any information, minimize

up(t k ), σ σ σ(t) Np k=0 n T p=1 u p (t k ), subject to x(t) = Φ Φ Φ(t, t 0 )x(t 0 ) + k i=0 n T p=1 Φ Φ Φ(t, t i )BR T (σ σ σ(t i ))w p u p (t i ), t k ≤ t < t k+1 , A L x(t) ≤ b L , 0 ≤ u p (t) ≤ u p,max , p = 1 . . . n T , -H i,max ≤ H i,rw (σ σ σ(t), σ σ σ(t)) ≤ H i,max , i = 1, 2, 3, -Ḣi,max ≤ Ḣi,rw (σ σ σ(t), σ σ σ(t), σ σ σ(t)) ≤ Ḣi,max , i = 1, 2, 3, x(t 0 ) = x 0 , σ σ σ(t 0 ) = σ σ σ 0 , ω ω ω(t 0 ) = ω ω ω 0 , x(t f ) = x f , σ σ σ(t f ) = σ σ σ f , ω ω ω(t f ) = ω ω ω f . (29)
Beside the fact that the equivalent planning problem ( 29) is integration free, it is still infinite dimensional. In the next section, it is shown how to make this problem tractable by means of parameterization and discretization.

Optimal control computation

In this section, the resolution method to the equivalent optimal control problem (29) is presented. The proposed methodology is based on a B-spline parameterization of the flat output (MRP) and the discretization of the infinite dimensional constraints. The result is a tractable NLP problem which needs an initial guess to be solved.

Non-linear programming description

First of all, it is considered that the N p +1 impulses application times are equally spaced through the manoeuvre time, t∈[t 0 , t f ], with timespan T =(t ft 0 )/N p , hence t k =t 0 +kT for k=0 . . . N p . These firing times will be denoted as nodes.

B-splines parameterization of the flat output

The attitude flatness property allows any kind of MRP time evolution parameterization. In this work, following [START_REF] Louembet | Path planning for satellite slew manoeuvres: a combined flatness and collocationbased approach[END_REF], B-splines, see [START_REF] Kress | Numerical Analysis[END_REF] for more details about them, are chosen to parameterize the flat output since they define flexible trajectories with a high degree of differentiability using a low number of parameters

σ σ σ(t) = nc j=1 a j B j,q (t), (30)
where the B j,q (t) are qth order B-splines built on the knots sequence, t knots ∈ R n knots , and the a j ∈R 3 coefficients are called control points.

Remark 2: The B-splines intrinsically assure continuity up to C q . Given the order q and the number of coefficients n c , the number of knots must satisfy n knots =n c +q+1.

The attitude profile has to be continuous up to its second derivative, hence, q≥2. On the other hand, it is chosen to have at least one control point to represent the attitude at each node plus four additional control points to impose σ σ σ and σ σ σ at the beginning and end of the manoeuvre. The previous consideration leads to n c =N p +5, therefore n knots =N p +q+6. The knots are selected as the nodes t k , augmented at left and right by repeating t 0 and t f t knots = [t 0 , . . . , t 0 , t 1 , . . . , t Np-1 , t f , . . . , t f ] T .

(31)

Discretization of time continuous constraints

The time continuous constraints are the path constraint related to the LOS region, see Eq.( 23), and the bounds on the reaction wheels angular momentum and its variation, Eq.( 25)- [START_REF] Moon | Quaternion Based Attitude Control and Suboptimal Rendezvous Guidance on Satellite Proximity Operation[END_REF]. Each one of these constraints is discretized with a time grid within each interval, k. The LOS constraint is gridded with n L equally spaced subintervals of duration T L =T /n L at which the constraint is imposed

A L x(t k,l ) ≤ b L , t k,l = t 0 + (k -1)T + lT L , k = 1 . . . N p , l = 1 . . . n L , (32)
whereas the reaction wheels constraints are gridded with n M equally spaced subintervals of duration

T M =T /n M -H i,max ≤ H i (t k,m , a j ) ≤ H i,max , i = 1, 2, 3, -Ḣi,max ≤ Ḣi (t k,m , a j ) ≤ Ḣi,max , i = 1, 2, 3, t k,m = t 0 + (k -1)T + mT M , m = 0 . . . n M . (33)

Discrete optimization problem

To ease the notation, following [START_REF] Gavilan | Chance-constrained model predictive control for spacecraft rendezvous with disturbance estimation[END_REF], a compact formulation of the discrete problem is developed. Defining the following stack vectors x S ∈R 6n L Np , u Sp ∈ R Np+1 and a S ∈R 3nc as

x S = [x T 1,1 , . . . , x T 1,n L , x T 2,1 , . . . , x T 2,n L , x T 3,1 , . . . . . . , x T Np,n L ] T , (34)
u Sp = [u p,0 , u p,1 , . . . , u p,Np ] T , (35)
a S = [a T 1 , a T 2 , . . . , a T nc ] T , (36)
and the stack matrices F∈R 6n L Np×6 and G p ∈R 6n L Np×(Np+1)

F = [Φ Φ Φ T (t 1,1 , t 0 ), . . . , Φ Φ Φ T (t 1,n L , t 0 ), Φ Φ Φ T (t 2,1 , t 0 ), . . . . . . , Φ Φ Φ T (t Np,n L , t 0 )] T , (37)
G p =                     Φ Φ Φ(t 1,1 , t 0 )BR T (σ σ σ 0 )w B p Θ 6×1 . . . Θ 6×1 . . . . . . . . . . . . Φ Φ Φ(t 1,n L , t 0 )BR T (σ σ σ 0 )w B p Θ 6×1 . . . Θ 6×1 Φ Φ Φ(t 2,1 , t 0 )BR T (σ σ σ 0 )w p Φ Φ Φ(t 2,1 , t 1 )BR T (σ σ σ 1 )w p . . . Θ 6×1 . . . . . . . . . . . . Φ Φ Φ(t 2,n L , t 0 )BR T (σ σ σ 0 )w p Φ Φ Φ(t 2,n L , t 1 )BR T (σ σ σ 1 )w p . . . Θ 6×1 Φ Φ Φ(t 3,1 , t 0 )BR T (σ σ σ 0 )w p Φ Φ Φ(t
                    , (38)
where σ σ σ k =σ σ σ(t k , a S ) represents the attitude at each node. The relation between the stack vectors and matrices defined in Eq.( 34)-( 38) is given by

x S = Fx 0 + n T p=1 G p (a S )u Sp . (39)
Now, the infinite dimensional problem (29) boils down to NLP, expressed with the compact formulation, by means of the above parameterization and continuous constraints discretization minimize

u Sp , a S n T p=1 u Sp 1 , subject to A LS n T p=1 G p (a S )u Sp ≤ b LS -A LS Fx 0 , 0 ≤ u Sp ≤ u Sp,max , p = 1 . . . n T , -H i,max ≤ H i (t k,m , a S ) ≤ H i,max , i = 1, 2, 3, -Ḣi,max ≤ Ḣi (t k,m , a S ) ≤ Ḣi,max , i = 1, 2, 3, A rend n T p=1 G p (a S )u Sp = x f -A rend Fx 0 , σ σ σ(t 0 , a S ) = σ σ σ 0 , σ σ σ(t 0 , a S ) = 0, σ σ σ(t 0 , a S ) = 0, σ σ σ(t f , a S ) = σ σ σ f , σ σ σ(t f , a S ) = 0, σ σ σ(t f , a S ) = 0, (40)
where A LS ∈R 5n L Np×6n L Np and b LS ∈R 5n L Np stack the LOS matrix (diagonally) and vector, see Eq.( 23), respectively. The parameters u Sp,max ∈R Np+1 are stack vectors whose components are all equal to u p,max . The matrix A rend = [Θ 6×6(n L Np-1) , Id 6×6 ] is employed to impose the rendezvous condition. It has been considered that reaction wheels kinetic momentum variation at initial and final time shall be zero which constrains σ σ σ. A NLP solver is required to obtain a solution of the static program (40).

Initial guess computation (hotstart)

Any NLP solver needs an initial guess to compute the optimal solution of problem [START_REF] Desai | Six-Degree-of-Freedom Trajectory Optimization Using a Two-Timescale Collocation Architecture[END_REF]. In this case, the process is composed of two steps; first, a traditional six-thrusters spacecraft model with three-degrees of freedom is employed to formulate and solve a linear programming (LP) problem; and then, this obtained LP solution is converted to NLP decision variables, u Sp and a S .

Six-thrusters problem formulation

Considering a pair of thrusters available for each LVLH axis, the control can be expressed at each node as u k =[∆V x,k , ∆V y,k , ∆V z,k ] T , hence, the translational states transition equation is linear

x(t) = Φ Φ Φ(t, t 0 )x 0 + k i=0 Φ Φ Φ(t, t i )Bu i , t k ≤ t < t k+1 , (41)
and the LP problem is posed as minimize

u k Np k=0 u k 1 , subject to A L x ( t k,l ) ≤ b L , -max(u p,max )/ √ 3 ≤ u k ≤ max(u p,max )/ √ 3, x Np = x f , A u0 u 0 = 0, A u Np u Np = 0, (42)
where the bounds on the impulse amplitude for each direction have been conservatively chosen to not overpass the upper bound of the thruster with more available impulse amplitude when all thrusters saturate (thus the use of √ 3). The purpose of the last linear constraints, expressed by means of the matrices A u0 and A u Np , is to make the initial and final orientations compatibles with the initial, σ σ σ 0 , and desired final attitude, σ σ σ f , respectively. It should be noted that the L 1 -norm term in the objective function is non-linear because u k can take both positive or negative values. However, this issue is avoided by adding optimization slack variables only allowed to take positive values.

Six-thrusters solution transformation to a NLP solution

Once the LP problem ( 42) is solved, the impulses amplitudes on each thruster are chosen as u 1,k = u k 2 and u p =1,k =0. The thruster labelled with p=1 is the one with higher u p,max .

The B-spline control points, a S , are obtained matching the N p +1 demanded orientations at the nodes by the LP solution. The MRP at the nodes can be obtained with the aid of the rotation angle and axis. First, denote by k * i , where the subscript i refers to the number of required thruster firings ( u k * i 2 >0), the nodes at which a non-null impulse amplitude is demanded or an attitude has to be reached (instant t Np ) and then compute the unitary vector z k * i representing the velocity increment orientation, expressed on the inertial frame since attitude is defined between the chaser body frame and the inertial frame, at these nodes

z k * i = [u x,k , u y,k , u z,k ] T / u k 2 , k * i = k, if u k 2 > 0. ( 43
)
Using z k * i , it is possible to obtain the rotation MRP, σ σ σ rot , between consecutive orientations. For the nodes without thruster firings ( u k 2 =0), the attitude at this node k is chosen as the value of the interpolated MRP, between the nodes k * i-1 and k * i , evaluated at the instant t k . The rotation MRP between t k-1 and t k is

σ σ σ rot k/k-1 = e k * i tan(s k θ k * i /4), s k = k -k * i-1 k * i -k * i-1 , t k-1 , t k ∈ [t k * i-1 , t k * i ]. ( 44
)
where the rotation angle and axis of Eq.( 44) are obtained using the previously computed orientations, see Eq.( 43)

θ k * i = acos(z k * i • z k * i-1 ), ( 45
)
e k * i = z k * i × z k * i-1 z k * i × z k * i-1 2 . ( 46
)
This way, smooth attitude transitions are obtained when some nodes do not have burnings. Since θ mi ∈[-π, π], no singularities arise when computing σ σ σ rot . Once the rotation MRP is obtained, it is possible to compute the MRP at each node t k . The MRP at the nodes are determined applying the attitude composition rule given by

σ σ σ k = (1 -σ σ σ rot k/k-1 2 2 )σ σ σ k-1 + (1 -σ σ σ k-1 2
2 )σ σ σ rot k/k-1 + 2σ σ σ k-1 × σ σ σ rot k/k-1 1 + ( σ σ σ rot k/k-1 2 σ σ σ k-1 2 ) 2 -2σ σ σ rot k/k-1 • σ σ σ k-1 . ( 47
)
The last step is to compute the control points for this nodes sequence. Imposing null σ σ σ and σ σ σ at t 0 and t f a linear system of 3n c equations with 3n c unknowns (remember that n c =N p +5) can be easily solved to obtain the initial guess Bsplines control points a S as it is proposed in [START_REF] Kress | Numerical Analysis[END_REF].

MPC scheme

Once the NLP problem ( 40) is solved, an open-loop solution for the rendezvous manoeuvre is available. However, disturbances, unmodelled dynamics, etc., will perturb the planned path while the spacecraft is manoeuvring, hence a MPC scheme, based on linearization around this previously computed solution, is developed in this section. The trajectory is recomputed on-line, in a sliding horizon framework, by solving a quadratic programming problem after each sampling interval which eases the computational burden (compared to the NLP) and does not need an initial guess. The terminal constraints are relaxed, considering them as terminal costs instead of constraints, to prevent feasibility issues and augment stability.

Linearized model

Allowing small increments of the decision variables, ∆u p,k and ∆a j , the translational transition Eq.( 20) states as follows

x(t) =Φ Φ Φ(t, t r )x r + r+k i=r n T p=1 Φ Φ Φ(t, t i )B R T (σ σ σ(t i ))w p (u p,i + ∆u p,i ) + r+nc j=r+1 ∆R p,aj (σ σ σ(t i ))u p,i ∆a j , t ∈ [t r+k , t r+k+1 ), (48)
where r=1 . . . N p is the current MPC step and the matrix ∆R p,aj ∈R 3×3 has the following expression

∆R (p,aj ),αβ = 3 β=1 ∂ R T αβ w p,β ∂σ β ∂σ β ∂a j,β σ σ σi,aj , α = 1, 2, 3, β = 1, 2, 3. ( 49
)
Following with the compact formulation developed through Eq.( 34)-( 39), one can define the following stack vectors ∆u Sp (r)∈R Np+1 , ∆a S (r)∈R 3nc

∆u Sp (r) = [∆u p,r , ∆u p,r+1 , . . . , ∆u p,r+Np ] T ,

S (r) = [∆a T r+1 , ∆a T r+2 , . . . , ∆a T r+nc ] T , (50) ∆a
and the stack matrix

H p (r)∈R 6n L Np×3nc H p (r) =              Φ Φ Φ(t r+1
Φ Φ Φ(t r+Np,n L , t i )B∆R p,ar+1 u p,i . . . r+Np i=r Φ Φ Φ(t r+Np,n L , t i )B∆R p,ar+n c u p,i             
, where one should note that many of the matrices ∆R p,a =Θ 3×3 because by definition the interval between two consecutive B-spline knots has at most q+1 non-null coefficients, see [START_REF] Kress | Numerical Analysis[END_REF]. Using the stack vectors and matrices, the following linearized translational states transition equation is obtained in compact form

x S (r) = Fx r + n T p=1 G p (a S )(u Sp + ∆u Sp ) + H p (u Sp , a S )∆a S , (52)
where the dependence with r has been omitted at the right-hand side of Eq.(52) for clarity.

Linearized planning problem

The linearized planning problem seeks the same objectives as the NLP problem [START_REF] Desai | Six-Degree-of-Freedom Trajectory Optimization Using a Two-Timescale Collocation Architecture[END_REF] but now the terminal constraints are included in the cost function as in [START_REF] Gavilan | Chance-constrained model predictive control for spacecraft rendezvous with disturbance estimation[END_REF] to encode a prescribed arrival time. This improves feasibility when considering disturbances, see [START_REF] Chen | A Quasi-Infinite Horizon Nonlinear Predictive Control Scheme with Guaranteed Stability[END_REF], and can improve asymptotic stability properties without needing terminal constraints, see [START_REF] Limon | On the Stability of Constrained MPC Without Terminal Constraint[END_REF]. At each MPC step, r, the linearized optimization problem to solve is minimize

∆u Sp ,∆a S n T p=1 r+Np i=r ∆u p,i + γ x r k=k0 (x Np+k -x f ) T Q x (x Np+k -x f ) + γ v r k=k0 x T Np+k Q v x Np+k + γ σ r k=k0 (σ σ σ Np+k -σ σ σ f ) T Id(σ σ σ Np+k -σ σ σ f ) + γ ω T 2 r k=k0 σ σ σ T Np+k Id σ σ σ Np+k , k 0 = max(0, r -N p ), subject to A LS n T p=1 G p (r)(u Sp (r) + ∆u Sp (r)) + A LS H(r)∆a S (r) ≤ b LS -A LS F(r)x r , 0 ≤ u Sp (r) + ∆u Sp (r) ≤ u Sp,max , |H i,rw (t r+k,m , a S (r)) + ∆H i,rw (t r+k,m , a S (r), ∆a S (r))| ≤ H i,max , | Ḣi,rw (t r+k,m , a S (r)) + ∆ Ḣi,rw (t r+k,m , a S (r), ∆a S (r))| ≤ Ḣi,max , ∆σ σ σ(t r , ∆a S (r)) = σ σ σ r -σ σ σ(t r , a S (r)), ∆ σ σ σ(t r , ∆a S (r)) = σ σ σ r -σ σ σ(t r , a S (r)), ∆σ σ σ(t r , ∆a S (r)) = σ σ σ r -σ σ σ(t r , a S (r)), -∆u Sp,max ≤ ∆u Sp ≤ ∆u Sp,max , -∆a S,max ≤ ∆a S ≤ ∆a S,max , (53)
where the variables σ σ σ r , σ σ σ r and σ σ σ r are the measured attitude and its derivatives at the end of each sampling interval. These values are necessary because the desired attitude path could suffer some deviations caused by reaction wheels saturation since the employed local reduction technique only guarantees time continuous constraint satisfaction at some discrete times. The increment on the reaction wheels kinetic momentum and its variation are

∆H i,rw (t r+k,m , a S (r), ∆a S (r)) = r+nc j=r+1 3 α=1 ∂H i,rw ∂σ α ∂σ α ∂a j,α σ σ σ r+k ,aj ∆a j,α , (54)
∆ Ḣi,rw (t r+k,m , a S (r), ∆a S (r)) = r+nc j=r+1

3 α=1 ∂ Ḣi,rw ∂σ α ∂σ α ∂a j,α σ σ σ r+k .aj ∆a j,α , (55)
The matrices associated with the terminal translational states costs terms are

Q x = Id 3×3 Θ 3×3 Θ 3×3 Θ 3×3 , Q v = Θ 3×3 Θ 3×3 Θ 3×3 Id 3×3 , (56)
and γ x , γ v , γ σ and γ ω are positive scalars that weight the relative cost of each one of the terminal conditions with respect to fuel consumption. Since both the translational states propagation, see Eq.( 52), and the flat output and its derivatives relation with the B-spline control points is linear, see Eq.( 30), the proposed objective function is quadratic, hence the optimization problem ( 53) is a QP problem.

MPC scheme

Summarizing the development of previous sections, the MPC scheme expressed as pseudocode is as follows The steps 2-4 are computed off-line while the vehicle is performing stationkeeping around the departure point waiting the command to start the manoeuvre so no hard real-time requirements appear when computing this solution. However, the steps 8-14 within the while loop are performed on-line during the manoeuvre which require a fast computation. That is the main reason why a QP problem based on linearization around a previously computed solution has been developed instead of solving the NLP problem at each step.

p w p u p,max [m/s] p w p u p,max [m/s] 1 [1, 0, 0] T 1 6 [0, 0, -1] T 1 2 [-1, 0, 0] T 1 7 [ √ 2, √ 2, 0] T /2 1 3 [0, 1, 0] T 1 8 [ √ 2, - √ 2, 0] T /2 1 4 [0, -1, 0] T 1 9 [- √ 2, √ 2, 0] T /2 1 5 [0, 0, 1] T 1 10 [- √ 2, - √ 2, 0] T /2 1
T m r f [2, 0, 0] T m v 0 [1, 1, -1] T m/s v f [0, 0, 0] T m/s θ θ θ 0 [0 • , 0 • , 0 • ] T θ θ θ f [90 • , 90 • , 90 • ] T ω ω ω B/L,0 [0 • , 0 • , 0 • ] T s -1 ω ω ω B/L,f [0 • , 0 • , 0 • ] T s -1
Table 2: Scenario 1 boundary conditions

Controller parameters

Regarding controller parameters, for both cases, the B-splines order is chosen to be quintic which is equivalent to take q=5 in Eq. [START_REF] Weiss | Model Predictive Control for Spacecraft Rendezvous and Docking Strategies for Handling Constraints and Case Studies[END_REF]. The discrete grids sizes, to evaluate the time continuous constraints, are chosen as n L =2 and n M =12, while The objective function weights are taken as γ x =10, γ v =5, γ σ =2 and γ ω =1. The LOS parameters for both cases are c y =c z =1/tan(π/4) and y 0 =z 0 =2.5 m.

Satellite with 10 thrusters

In this scenario, a conventional cargo satellite with 10 thrusters has to rendezvous with a target flying in an eccentric low Earth orbit with e=0.1, h p =600 km and ν(t 0 )=π/4. Table 1 shows the characteristics of the considered propulsive layout. On the other hand, the chaser inertia matrix is chosen to be similar to the russian Progress cargo spacecraft, see [START_REF] Fehse | Automated Rendezvous and Docking of Spacecraft[END_REF],

I =   31 0 0 0 31 0 0 0 5   • 10 3 kg • m 2 , (61)
whereas the bounds of the reaction wheels angular momentum and its variation are taken as H i,max =500 N•m•s and Ḣi,max =20 N•m, respectively. At the beginning, the angular momentum of the system is considered to be null H tot =0. The manoeuvre boundary conditions are given by Table 2. The considered disturbance parameters for this simulation are δ δ δ=0.0175, Σ δ,ij =0.0175δ ij , δ ūp =0.02 and σ δup =0.05.

Satellite with 2 thrusters

In this scenario, a lightweight satellite with only 2 available thrusters has to rendezvous with a target flying in an eccentric low Earth orbit with e=0.5, 3. For this case, the chaser inertia matrix is chosen to be the one corresponding to the CNES small satellite MYRIADE, see [START_REF] Louembet | Path planning for satellite slew manoeuvres: a combined flatness and collocationbased approach[END_REF]

T m r f [2, 0, 0] T m v 0 [1, 1, -1] T m/s v f [0, 0, 0] T m/s θ θ θ 0 [0 • ,0 • , 0 • ] T θ θ θ f [90 • , 90 • , 90 • ] T • ω ω ω B/L,0 [0 • , 0 • , 0 • ] T s -1 ω ω ω B/L,f [0 • , 0 • , 0 • ] T s -1
I =   40 -3 -0.5 -3 28 -1 -0.5 -1 45   kg • m 2 , (62)
whereas the bounds of the reaction wheels angular momentum and its variation are taken as H i,max =1 N•m•s and Ḣi,max =0.05 N•m. At the beginning the angular momentum of the system is considered to be null, H tot =0. In this case, the manoeuvre boundary conditions are shown in Table 4 The considered disturbance parameters for this simulation are δ δ δ=0, Σ δ,ij =0.0175δ ij , δ ūp =0 and σ δup =0.01.

Simulation results

For each scenario, 100 realizations for the chosen disturbance parameters are simulated. Then, the obtained results are shown and discussed.

Scenario with 10 thrusters

First, analyse the scenario with 10 thrusters. For all the realizations, the linear QP program is feasible and the chaser reaches the proximity of the target without trespassing the LOS region, see Fig. 4. A typical 3D path of a realisation is shown in Fig. 3 while the attitude profile is shown in Fig. 5. For the shown realization the desired orientation is met at the end while the angular velocity is driven to a quasi-null value due to the considered uncertainties. More details on the terminal accuracy for this scenario are given in Table 5 where δ measures the mismatch between the obtained and the desired terminal value. Regarding the planned impulses, for the plotted realisation, see Fig. 10 significantly along the manoeuvre. Regarding the cost, the NLP program reduces fuel consumption in a 21.054% compared to the converted solution from the LP problem, see Table 5. Finally, in Fig. 7 it is shown that the reaction wheels have saturations (both on angular velocity and acceleration) at the initial and final instants of the manoeuvre, but then desaturate immediately and keep providing torque.

Scenario with 2 thrusters

Analysing the second scenario with 2 thrusters, similar conclusions with the first scenario still holds, see Fig. 8. Note that the spacecraft in this case is underactuated in translational control. Moreover, the desired final orientation is not favourable at all to brake the spacecraft since the thruster 1 nozzle will end pointing to the +x axis and the thruster 2 nozzle to the -z axis. The terminal accuracy is shown in Table 6 (it shows higher accuracy than the 10 thrusters scenario due to the lighter perturbations). For the plotted realization of Fig. 10, it is shown that the final braking impulse has to be advanced one interval due to the non favourable last orientation. In this case, there is not an improvement in fuel consumption when compared to the obtained LP solution but the reaction wheels saturation peak has been lowered from 5.5608 N•m•s (LP solution converted to NLP solution) to 1 N•m•s, see Table 6.

Concluding remarks

This paper has presented a predictive guidance and control algorithm for sixdegrees of freedom spacecraft rendezvous based on the translational state transition matrix, the attitude flatness property, discretization and a MPC scheme based on linearization. One of the main contributions of the proposed algorithm is its ability to consider several chaser spacecraft configurations which not only reduces fuel consumption but also allows to consider propulsive and ACS constraints in an integrated framework.

The numerical experiments shown in Section 6 have validated the method for two different spacecraft configurations. Additionally, the simulations have demonstrated convergence of the proposed MPC to the desired final state even in the presence of disturbances. However, a formal proof of stability has not been addressed and is left as future work.

Possible future research lines include the following. First, it will be of great interest to consider on/off thrusters as it is done in [START_REF] Vazquez | Pulse-width predictive control for LTV systems with application to spacecraft rendezvous[END_REF]. This will cause continuous coupling between translational and rotational motion since the vehicle will be spinning when the thrusters are fired. Second, to improve the robustness of the underactuated case, robust MPC techniques in the spirit of [START_REF] Gavilan | Chance-constrained model predictive control for spacecraft rendezvous with disturbance estimation[END_REF] could be considered. Finally, another possible line is to consider more advanced techniques, that does not rely on discretization, to handle the time continuous constraints of the problem.

Figure 1 :

 1 Figure 1: LVLH frame

Figure 2 :

 2 Figure 2: LOS region

Algorithm 1 :begin 2 3 4 5 6 Initialize the current MPC step r=1; 7 while r ≤ N p do 8 9 u p,Np+r =0, 10 σ

 12345678910 MPC scheme 1 Obtain a solution of the LP problem (42); Transform the LP solution to NLP decision variables using Eq. (43)-(47); Obtain a solution of the NLP problem (40), u Sp and a S ; Apply u p,0 and Ḣrw (t) for t∈[t 0 , t 1 ); Prescribe reference controls at N p + r: σ σ(t Np+r , a S (r))=σ σ σ f , σ σ σ(t Np+r , a S (r))=0, σ σ σ(t Np+r , a S (r))=0 -→ a S (r); 11 Obtain the solution of the QP linearized rendezvous problem 53, ∆u Sp (r) and ∆a S (r) ; 12 Update the decision variables u Sp (r)=u Sp (r)+∆u Sp (r), a S (r)=a S (r)+∆a S (r); 13 Apply ∆V p,r and Ḣrw (t) for t∈[t r , t r+1 ); 14 Update the current MPC step, r=r+1;

x

  , the thrusters {1,2,6,9} have relevant firings while thrusters {3,4,5,7,8,10} are not operated 0

Figure 3 :

 3 Figure 3: Chaser 3D path of scenario 1 for the first random realization

Figure 4 : 3 Figure 5 : 10 Figure 6 :

 435106 Figure 4: Chaser trajectory on the target orbital plane of scenario 1 for all random realizations

Figure 7 :

 7 Figure 7: Reaction wheels angular momentum and its variation of scenario 1 for the first random realization

Figure 8 :

 8 Figure 8: Chaser trajectory on the target orbital plane of scenario 2 for all random realizationsResults (intrinsic Euler angles sequence 3→1→3)J LP =3.49 m/s, J N LP =3.40 m/s E[J M P C ]=3.43 m/s, s[J M P C ]=0.029 m/s E[ δr(t f ) 2 ] 0.82 m s[ δr(t f ) 2 ] 0.36 m E[ δv(t f ) 2 ] 1.34 cm/s s[δv(t f )] 0.64 cm/s E[θ θ θ(t f )][97.39 • , 93.73 • , 89.10 • ] s[θ θ θ(t f )] [4.57 • , 2.86 • , 2.51 • ] E[ ω ω ω(t f ) 2 ] 0.31 • /s s[ ω ω ω(t f ) 2 ] 0.13 • /s

3 Figure 9 : 2 Figure 10 :Figure 11 :

 3921011 Figure 9: Chaser attitude of scenario 2 for the first random realization

  [START_REF] Roscoe | Overview and GNC design of the Cubesat Proximity Operations Demonstration (CPOD) Mission[END_REF],n L , t 1 )BR T (σ σ σ 1 )w p . . .

				Θ 6×1
	. . .	. . .	. . .	. . .
	. . .	. . .	. . .	. . .

Φ Φ Φ(t Np,n L , t 0 )BR T (σ σ σ 0 )w p Φ Φ Φ(t Np,n L , t 1 )BR T (σ σ σ 1 )w p . . . BR T (σ σ σ Np )w p

Table 1 :

 1 Thrusters configuration for scenario 1

	Boundary conditions (intrinsic Euler angles sequence 3→1→3)
	t 0	0 s	t f	900 s
	r 0	[400, -250, -200]

Table 3 :

 3 Thrusters configuration for scenario 2

	Boundary conditions (intrinsic Euler angles sequence 3→1→3)
	t 0	0 s	t f	900 s
	r 0	[350, 200, 200]

Table 4 :

 4 Scenario 2 boundary conditions h p =400 km and ν(t 0 )=π. The thrusters are mounted in an orthogonal configuration as shown by Table

Table 5 :

 5 Scenario 1 terminal results

Table 6 :

 6 m/s E[ δr(t f ) 2 ] 0.82 m s[ δr(t f ) 2 ] 0.36 m E[ δv(t f ) 2 ] 1.34 cm/s s[δv(t f )] 0.64 cm/s E[θ θ θ(t f )][97.39 • , 93.73 • , 89.10 • ] s[θ θ θ(t f )] [4.57 • , 2.86 • , 2.51 • ] E[ ω ω ω(t f ) 2 ] 0.31 • /s s[ ω ω ω(t f ) 2 ] 0.13 • /s Scenario 2 terminal results

Acknowledgements

The authors gratefully acknowledge Universidad de Sevilla for funding part of this work under its V-PPI US. Rafael Vazquez acknowledges financial support of the Spanish Ministerio de Economía y Competitividad under grant MTM2015-65608-P.

Simulation results

Since the employed formulation does not make any assumptions on the chaser number of thrusters, two different scenarios will be considered. The first one will correspond to a heavy rendezvous satellite equipped with 10 thrusters while the other one corresponds to a low-power spacecraft equipped with 2 thrusters. The simulations of this sections have been obtained using MATLAB routines with Gurobi optimization package, see [START_REF]Gurobi optimizer reference manual[END_REF], as LP and QP solver whereas the IPOPT optimization package, see [START_REF]Introduction to IPOPT: A tutorial for downloading, installing and using IPOPT, COIN-OR[END_REF], is used as NLP solver.

Rendezvous model

It is important to remark that although a linear model, see Eq.( 1)-( 3), is used to compute the control sequence, the plant is considered to be dominated by the following non-linear relative motion dynamics, see [START_REF] Wie | Space Vehicle Dynamics and Control[END_REF],

(59)

Disturbances model

In a similar way as [START_REF] Gavilan | Chance-constrained model predictive control for spacecraft rendezvous with disturbance estimation[END_REF] (note that there the disturbance is considered in an additive way), a disturbance on each of the thrusters performance is added to test the capabilities of the MPC scheme developed in Section 5. This disturbance is modelled in the chaser body axes as

where u p is the commanded output computed by the control laws, δθ δθ δθ ∼ N 3 (E[δθ δθ δθ], Σ Σ Σ δ ) is a vector of random small angles and δu p ∼ N (E[δu p ], s[δu p ]) is a random scalar. These disturbances model several physical aspects. First, the attitude control of the chaser will not be perfect, so one can expect some alignment errors, modelled by Ω(δθ δθ δθ) in a simplified way. On the other hand, with δu p one can model thrust level disturbances.

Simulation scenarios

To test the capabilities of the proposed algorithm, two scenarios for different pursuer architectures are considered. For the first scenario, a heavy spacecraft equipped with 10 thrusters is considered while for the second scenario a lightweight satellite with a limited propulsion plant, with only 2 thrusters is simulated.