Trajectory Planning and Tracking for a Re-Entry Capsule with a Deployable Aero-Brake

: In the last decade, the increasing use of NanoSats and CubeSats has made the re-entry capsule an emerging research ﬁeld needing updates in conﬁguration and technology. In particular, the door to advancements in terms of efﬁciency and re-usability has been opened by the introduction of inﬂatable and/or deployable aerodynamic brakes and the use of on-board electronics for active control. Such technologies allow smaller sizes at launch, controlled re-entries, and safe recovery. This paper deals with the design of a guidance and control algorithm for the re-entry of a capsule with a deployable aero-brake. A trajectory optimization model is used both in the mission planning phase to design the reference re-entry path and during the mission to update the trajectory in case of major deviations from the prescribed orbit, thanks to simpliﬁcations aimed at reducing the computational burden. Successively, a trajectory tracking controller, based on Nonlinear Model Predictive Control (NMPC), is able to modulate the opening of the aero-brake in order to follow the planned trajectory towards the target. A robustness analysis was carried out, via numerical simulations, to verify the reliability of the proposed controller in the presence of model uncertainties


Introduction
The reorientation of NASA's priorities after shuttle retirement and the current increasing number of NanoSats and CubeSats have brought the re-entry vehicle technology back in the news.The decrease of size, mass, and power results in a reduction of the overall mission costs, increasing the accessibility to space [1].However, despite the remarkable research and technology developments over time, re-entry capsules have been characterized by the same basic design concept.As a consequence of the miniaturization, more solutions are needed to offer the possibility to recover a payload and potential information for post-flight analyses coming from space [2,3].
Furthermore, several research projects have worked on the design of on-board electronics for small re-entry vehicles, increasing the efficiency and the re-usability, by enabling actively controlled re-entry and safe recovery [4].
To this purpose, inflatable and deployable aerodynamic brakes have been proposed and developed in the last few years.Umbrella-like or inflatable re-entry capsules are usually lighter and less expensive than conventional ones.They allow small sizes at launch and a sufficiently large surface area in re-entry.Thanks to their low ballistic coefficient, reducing the peak heat flux and the mechanical load, these kinds of capsule are usually based on commercial off-the-shelf (COTS) components and materials, minimizing development, building, and launching costs.
In the literature, several papers deal with inflatable or deployable aero-brake devices for capsule re-entry missions.Reference [5] introduced the concept of the parashield with application to a variety of capsules, including advanced manned spacecraft.In [6], the authors focused on the heat shield design, the flight dynamics, and the thermal loads in the re-entry mission of BREM-SAT 2, where the parashield resembles a reinforced umbrella, increasing the front area by a maximum factor of twelve.Reference [7] presented a drag device module, capable of de-orbiting a 15 kg CubeSat from a 700 km circular orbit in under 25 years.The geometry of such a drag device is also able to provide a three-axis attitude stabilization of the host CubeSat.Reference [8] focused on a flare-type membrane aeroshell sustained by an inflatable torus for an experimental vehicle of 15.6 kg.In [9], the authors analysed the aerodynamic longitudinal stability of a suborbital re-entry demonstrator in several conditions, in order to implement a proper re-entry strategy.Reference [10] showed the structural analysis and testing of the Inflatable Re-entry Vehicle Experiment (IRVE) [11] inflatable structure.In [12], the authors showed the design of a semi-rigid, deployable hypersonic decelerator system for heavy vehicles for manned missions on Mars.In this case, the architecture is similar to an umbrella, where the aerosurface is a flexible thin skin, used as thermal protection, which is draped over a structure of high-strength ribs.
Several challenges need to be addressed in guidance and control, due to the nonlinearity of the flight dynamics model, the constraints on the state and input variables, and the parametric uncertainties of the flight.Constraints about the heating rate and structural loads on the capsule can be considered in the optimal control problem as in [13].Reference [14] considered flight time constraints and no-fly zones in the optimization problem, in order to propose an analytical solution for a formation of Hypersonic Glide Vehicles (HGVs).The guidance logic in [15] computes the bank reversals by evaluating information from the reference cross-range profile, current cross-range, and estimated actual lift-todrag ratio; an noniterative numerical predictor, in the last phase of the trajectory, chooses whether a final bank reversal is needed to make null the heading error.In [16], the authors divided the re-entry mission into two phases and focused their attention on the second one, being the atmospheric re-entry, adding heating constraints to the trajectory optimization problem, where the design variable is the bank angle.
Optimal control methods were presented in [17,18].More specifically, in [17], the author compared five different types of optimum guidance and control algorithms, while in [18], the Pontryagin principle was used to generate an autonomous slew trajectory, and its performance was compared with a sinusoidal approach, whereas in [19], the authors improved the previous trajectory planning algorithms, with autonomous collision avoidance schemes for the re-entry.
Indeed, the re-entry of a capsule is a critical phase of the mission, where model uncertainties and environment disturbances can strongly affect the landing point on the ground.Non-lifting vehicles go through a ballistic trajectory, where the only control ability is achieved by changing the drag coefficients, retracting or deployingthe aero-brake [20].
Usually, the re-entry phase is divided into two main phases: de-orbit at very high altitude and atmospheric entry.The former is a long stage in a low-density atmosphere, where the decaying trajectory planning is fundamental to find the right starting point, while the latter is the most critical stage for the structure, being stressed by very high thermal flows [8,21].Therefore, in the atmospheric entry, the deployable aero-brake acts also as a thermal shield, resulting in reduced manoeuvring capabilities.
For all these reasons, trajectory tracking control becomes a challenging problem, and the uncertainties in air density and aerodynamic drag suggest the development of adaptive, parameter varying, or optimization-based closed-loop techniques aimed at minimizing the tracking error [22].
The actuated aero-brake has the ability to control the trajectory and eventually reaching a predetermined landing position by changing its shape, and consequently its wetted surface, during the descent flight [9,23,24].Reference [25] presented a de-orbit and guidance scheme to control the trajectory of a vehicle with a passive drag-based system, defined as the Exo-brake.Reference [26] dealt with a method for drag-controlled re-entry, applying to a retractable drag de-orbit device that can be attached to CubeSats and larger satellite structures.In [27], the authors designed a simple Proportional-Integrative-Derivative (PID) controller for satellite formation keeping, modulating differential drag.The scope of [28] was to mitigate casualty risk due to an undesired re-entry, by modulating drag in the last orbit revolutions, with the objective of increasing the probability to impact in the Southern Hemisphere.Reference [29] proposed an algorithm that computes the drag profile needed to de-orbit in a desired location.
Unfortunately, although open-loop control techniques were used in successful missions, they cannot react to environmental unpredictable uncertainties, turning out to be often unreliable.The trajectory tracking performance can be strongly improved by means of a closed-loop control, in terms of both reliability and tracking error.
In the literature, several approaches can be found to address the problem.The most typical approach is based on PID trajectory tracking controllers [30] and the use of gain scheduling techniques, where the trajectory is computed in advance and controller parameters can change during the decay depending on the state of the capsule [31][32][33].While [34] proposed a time-varying feedback controller, in [35][36][37], the re-entry strategy was based on fuzzy logic to track the drag reference profile.Gain scheduling was proposed also with optimal control techniques in [38][39][40].
Another approach to track the reference drag profile during the re-entry phase was proposed in [41] with a nonlinear PID control law able to guarantee, in the absence of control saturation, the global asymptotic stability.A sliding mode observer to estimate and track the variation of the drag coefficient was proposed in [42], while [43,44] was based on feedback linearization.The nonlinear dynamic inversion technique was used in [45,46].
In this paper, a unified approach to plan and control the re-entry mission of the capsule with a deployable aero-brake is proposed.This is based on an optimization strategy that covers both the off-line calculation of the optimal trajectory and the on-line possible path re-planning, taking into account the system dynamic response.The objective of the optimization is to find a (sub)optimal trajectory to reach a predefined target point at the Terminal Area Management (TAEM) fixed at an altitude around 30 km.The flight path must be compliant with the system dynamics and constraints, as well as the maximum admissible value of the heat flux, in order to guarantee a safe re-entry trajectory for the capsule.
While the off-line flight path planner relies on a Genetic Algorithm (GA) able to fully exploit the search domain, the on-line planner is based on a Sequential Quadratic Programming (SQP) algorithm to update the planned trajectory in case of deviation from the prescribed path computed off-line.SQP reduces the on-board computational burden and provides a sub-optimal solution in a finite and predictable time.
Once the re-entry has begun, the planned path is frozen and a Nonlinear-Model-Predictive-Control (NMPC)-based trajectory tracking algorithm modulates the aero-brake, to minimize the tracking error and reach the desired target position at a given altitude, ensuring the feasibility of the trajectory under suitable constraints [47,48].
The use of MPC is not new in the control of re-entry vehicles, and it was introduced in [49,50].In [51], the authors presented preliminary results on an atmospheric re-entry MPC controller, which does not make use of any pre-computed nominal trajectory, focusing directly on the target point.The present paper extends the use of MPC to the whole descent phase, including de-orbiting and considering an adaptive objective function to be minimized, whose weights depend on the actual vehicle state.
To demonstrate the effectiveness of the proposed strategy, a campaign of numerical simulations was executed making a sensitivity analysis on the main parameters of the control system and, so, guiding the designer to their optimal selection for balancing performance and computational burden.
Once the parameters were chosen, an extensive robustness analysis was carried out in the presence of model uncertainties, orbital perturbations, and measurement noises, to test further the capabilities of the proposed control scheme.
The paper is structured as follows: Section 2 deals with the problem statement, with details about the mission and the model; in Section 3, the proposed guidance system is described, with details the path planning procedure, the NMPC-based control scheme, and the adaptive strategy in the different phases of the re-entry; Section 4 resumes results about the path planning, sensitivity, and robustness analysis to prove the effectiveness of the proposed controller.

Mathematical Model
The satellite configuration consists of a detachable capsule and a service module, containing every subsystem needed during the in-orbit operations.The main feature of the capsule is an umbrella-like front structure, which can be used both as an aero-braking device and as a thermal shield.
Such a design allows at the same time slowing down the satellite and protecting the capsule from the extremely high temperatures reached during the atmospheric entry.
The re-entry mission is made-up of four main phases: • De-orbiting: At an altitude H 1 = 300 km, the satellite leaves its nominal orbit and begins the slow decay towards the Earth's surface, modulating the aero-brake in order to control the trajectory of the satellite.• Atmospheric approach: When the altitude H 2 = 150 km is reached, the service module releases the capsule, which, after the separation, approaches the atmospheric entry through a ballistic trajectory, controlled by the modulation of the deployable front structure.

•
Atmospheric entry: At the height H 3 = 100 km, the capsule enters the atmosphere and follows a purely ballistic path until the Terminal Area Energy Management (TAEM) interface at H 4 = 30 km.

•
Terminal area phase: In this phase, the capsule follows an uncontrolled vertical descent path until H 5 = 10 km, when the parachute is deployed.
In the whole mission, mainly two kinds of forces act on the capsule: gravity and aerodynamics.
Usually, the gravity acceleration between Earth and the satellite is considered as the gradient (∇ operator) of the gravitational potential W [52]: Let G be the universal gravitational constant and V ⊕ and σ ⊕ the Earth's volume and mass density, respectively.In the Earth-Centred Inertial (ECI) reference frame, the sum of the contributions of each element of the Earth's mass dM ⊕ is the gravitational potential: where the symbol • is the Euclidean norm, R c is the position vector of the individual mass elements dM ⊕ = σ ⊕ dV from the centre of the reference system, and R ECI = [X ECI , Y ECI , Z ECI ] T is the position vector of the satellite in the ECI frame.Consider a spherical coordinate system (r, φ, λ) such that: where r = R ECI is the distance of the capsule from the Earth's centre, λ (assumed positive towards the east) represents the longitude, and φ is the geocentric latitude.
In empty space (r > R c ), the gravity potential satisfies the Laplace equation; consequently, it can be expressed as an expansion in series of spherical harmonics: where R ⊕ is the Earth's equatorial radius.P lm = P lm (sinφ) are the associated Legendre functions of degree l and order m, and C lm and S lm are the Stokes coefficients [52] describing the dependence of gravity potential on the Earth's shape.Geopotential coefficients with l = 0 are called zonal coefficients, since they describe the part of the potential that does not depend on the longitude.In this case, the coefficients S l0 are equal to zero and the coefficients C l0 are denoted as −J l .
Sectoral coefficients are obtained with m = l.They describe the part of the gravity potential that depends only on the longitude.
C lm and S lm with m = l are called tesseral coefficients and describe the part of the gravity potential that depends both on the longitude and latitude.
Equation ( 4) can be efficiently approximated by considering only the first four zonal coefficients: with µ = GM ⊕ .In ( 5), the dominating term related to the first zonal coefficient J 1 was extracted from the summation.
Considering the transformation from polar coordinates to ECI Cartesian coordinates and using ( 5), the gravity acceleration in the ECI frame can be approximated as follows: where a J 2 (R ECI ), a J 3 (R ECI ), and a J 4 (R ECI ) are the perturbation accelerations related to zonal harmonics J 2 , J 3 , and J 4 .
The terms a J 2 , a J 3 , and a J 4 , as well as the other ones related to neglected zonal, tesseral, and sectoral coefficients are due to Earth's shape: the oblateness causes an asymmetry of the gravitational field, which results in a variety of perturbations in the satellite orbit, affecting satellites at low altitudes more.
The atmospheric drag represents the largest non-gravitational perturbations acting on low-altitude satellites.The lift contribution of a blunt body is negligible, whereas the aerodynamic drag plays a significant role more so in the atmospheric entry, due to the presence of the Earth's atmosphere.The drag force can be written as: where: • ρ is the air density, depending by the altitude H; T is the satellite velocity vector in the ECI reference frame; S is the surface of the deployable aero-brake.
The drag coefficient C D describes the interaction of the atmospheric constituents with the satellite surface.More specifically, it depends on the satellite surface material, the chemical constituents of the atmosphere and the molecular weight and temperature of the impinging particles.In the free molecular flow regime, characterized by a Knudsen number greater than 1 (Kn > 1), the particles re-emitted from the satellite do not interfere with the incident molecules.At lower altitudes, with Kn < 0.1, the re-emitted molecules partially shield the satellite from the incident flow.This phenomena causes a drag decrease [52,53].
The model of the Earth's atmosphere is not simple, and an accurate description would require many parameters.
The largest source of uncertainties in the evaluation of the atmospheric density at the highest altitudes is related to the solar activity.The Sun always emits a different amount of energy, changing the portion of energy penetrating and heating the atmosphere.The effect is a continuous change in density, extremely hard to predict, Solar activity being an almost unpredictable phenomenon.
In this paper, the USA Standard Atmosphere model was assumed [54]: where ρ is the air density at any altitude above sea level H, ρ 0 is the air density at some specified reference altitude H 0 , and SH is the atmospheric scale height at H 0 .
In the exponential model, the solar activity is approximately taken into account through the atmospheric scale height.In fact, the atmosphere can been divided into 28 sections, from 0 up to 1000 km, each having its adequate scale height [23].
Assuming an Earth-Centred Inertial (ECI) frame, the nonlinear motion equations are defined as follows: where M is the satellite mass.Equation ( 9) creates a relationship between the ballistic coefficient C b = M C D S and the velocity vector of the satellite; consequently, the control of the trajectory can then be achieved through the modulation of the aero-brake surface.
The aero-brake has also the role of a thermal shield in the hottest phase of re-entry.
Chapman's and Sutton-Grave's simplified methods allow computing the heating rate experienced by a re-entry vehicle [55]: where K is a suitable constant, ||V ECI || is the vehicle speed, and l nose is the capsule nose radius, which measures the "bluntness" of capsule's shape.
The altitude where the maximum heating rate is reached can be found using the following relation [55]: where γ is the flight-path angle and C b = M C d S is the ballistic coefficient of the capsule.

Guidance System Design
The scheme of the proposed guidance and control algorithm is shown in Figure 1, with the off-board and on-board parts.Before the whole mission, a path planning block is able to compute the optimal trajectory for the overall descent path.During the mission, before the de-orbiting phase, such a trajectory can be updated in order to take into account possible deviations from the prescribed orbital parameters by sending a replanning trigger to the local path planning subsystem.Once the re-entry has begun, the trajectory is considered frozen and the re-entry mission control block becomes active.Two MPC-based trajectory tracking algorithms are considered, with different optimization parameters fitting the model to the specific de-orbiting and atmospheric entry phases.

Path Planning
The goal of the path planning algorithm is to minimize the final position error with respect to a predefined target, by finding the optimal values of the control signal, i.e., the values assumed by the aero-brake surface along the trajectory, and the de-orbit point, i.e., the satellite on-orbit initial position, identified by the true anomaly.
The control signal, S(t), is assumed piecewise constant to avoid damages of the aerobrake actuation system caused by overloads due to consecutive openings and closures.
Given N h time ranges during the re-entry, the design variable vector can be defined as follows: where: • ν is the true anomaly of satellite when the descent phase starts; • S i with i ∈ 1, 2, . .., N h − 3 is the aero-brake surface to be required by the control system in the time range [t i , t i+1 ]; • S 0 , S N h −2 and S N h −1 are the aero-brake surfaces at 300, 150, and 100 km.
The choice of relating the last two aero-brake surface values to the altitude in the atmospheric entry depends on the strong change in the drag coefficient and mass after the capsule leaves the hub.
The guidance algorithm is based on the following optimization problem: subject to: with the orbital period P and a maximum descent time t max , evaluated by Cowell's method [56].It is worth noting that, in order to avoid unrealistic values of the aero-brake surface, the constraint S i ∈ S takes into account the mechanical limits and the quantization of the actuation system.Practically, such a set S is defined by loading a table containing all the values of the aero-brake surface actually achievable.
The objective function F to be minimized is: where: • R * ECEF is the desired position vector in the ECEF reference frame at the altitude H * = 30 km; • R * ECEF is the effective position in the ECEF reference frame at the altitude H * = 30 km; • qmax is the value of the heat flux peak along the trajectory; • qmax is the maximum admissible value of the heat flux peak, assumed as 600 kW m 2 ; • Q r is a suitable definite positive weights matrix and Q q is a positive scalar.
The objective function essentially is the sum of two contributions.The former is the position error with respect to the target point at altitude H * = 30 km; the latter is a term assessing the difference between the heat flux peak reached along the trajectory and its maximum allowable value.
The defined optimization model is strongly nonlinear due to the capsule dynamics.However, the off-line path planning phase can trust on ground-based equipment, so, to achieve a global search of the optimum, a Genetic Algorithm (GA) procedure was implemented, with the integration of Equation ( 9) with a fixed step size solver to compute the objective function and solve the optimization problem (13).
The on-board guidance and control system is provided with a local path planning block to update the trajectory whenever the navigation system estimates a deviation from the prescribed mission orbital parameters, before the de-orbiting phase begins.In this case, a SQP-based nonlinear local search [57] is carried out, starting from the reference trajectory, leading to a new sub-optimal solution.

MPC-Based Trajectory Tracking
The goal of the control system is to minimize the tracking error of the reference trajectory, provided by the guidance system, for the whole descent phase.On the other hand, the control system must guarantee an error at the target landing point less than 250 km.Considering that the capsule follows a quasi-vertical path below the target altitude H = 30 km, the requirement can be converted into a maximum admissible position error at the TAEM interface less than 200 km.
In the hottest phase, in order to mitigate the thermal stress on the capsule, the aerobrake surface is kept fixed equal to the reference value S 5 = 0.770 m 2 .Consequently, the optimal control law can only be calculated until the altitude H = 100 km.
To achieve the above-defined goal, taking into account the control signals' bounds, a Model Predictive Control (MPC) strategy is proposed.Following the MPC approach, the control law is computed by minimizing a given cost function, which weights the difference between the predicted state trajectory and the reference one.According to the Receding Horizon paradigm, only the first move of the discrete sequence is applied to the controlled system.
The predicted state trajectory is computed by using a discrete-time state space model by discretizing Equation (9) with a suitable sample time T c .
Let k be the actual discrete time step, with x(k) = R T ECI (k), V T ECI (k) T the state vector composed by the satellite position and velocity in ECI frame, and let u(k) = S(k) be the control input equal to the aero-brake surface.Consider a time window of amplitude T p composed by N p steps.By using the mathematical model, it is possible to predict the future behaviour of the system.At each tome step k, the predicted state x(k + i|k) (i = 1, . . ., N p ) depends on the control sequence U(k) = [u(k|k), . .., u k + N p − 1 k ] T , which must be optimized, and by the initial state x(k), which is measured or, at least, estimated.In order to reduce the computation burden, only N c (with N c ≤ N p ) control moves are optimized by assuming that the last N p − N c control moves are constant u(k The cost function to be optimized is: where R ECI (k + i|k) and V ECI (k + i|k) are the predicted position and velocity vector, respectively, R ECI (k + i) and V ECI (k + i) are the reference position and velocity vectors, S(k + i) is the reference surface of the aero-brake, and S(k + i|k) represents the predicted value of the control variable.The objective function is then composed of three parts, weighted by the scalar values p(i), q(i), and w(i): the first two terms weight the error between the reference position and speed vectors and the predicted ones; the last term attempts to minimize the difference between the control variable and the reference one.
The optimal control law U * (k) is calculated by solving the following problem: subject to x(k where x(k) is the current measured state at discrete time k.Equation (19) takes into account the discrete-time implicit state space model defined by a suitable discretization of (9).Inequality constraints described by (18) allow taking into account the bounds of control signal and can be rewritten as follows: Given the optimal control sequence, according to the Receding Horizon technique, only the first move u * (k|k) is applied.The reason is that the remaining sequence u * (k + i|k) (with i = 1, . . ., N c − 1), though based on a prediction of the system behaviour, is an open-loop control law, affected by disturbances and uncertainties.
At the next time step k + 1, on the basis of a new current state measurement, new predictions and optimization can be performed.Hence, the optimization process can be iterated.
In order to minimize the computational burden and ensure the existence of a solution in a reasonable time, due to the nonlinear nature of the plant model, it is assumed that the control horizon is made up by a single time step (N c = 1) so that, at each step, a single control move is optimized.In this case, the MPC problem can be numerically solved with the golden section method, which guarantees the convergence of the solution to a local minimum.The maximum number of iterations needed to achieve the solution can be evaluated on the basis of the requested resolution of the optimal control variable, which is related to the quantization of the deployable umbrella control command α [58].
In the proposed guidance and control scheme, two MPC-based trajectory tracking algorithms are considered, based on the same model, but using different weights for the objective function definition, in order to better suit the specific re-entry phases.
In particular, in the former phase, the controller objective is to minimize the error over a finite limited prediction horizon, in order to follow the planned path.During this long phase, the sampling frequency is low (<10 −3 Hz) to minimize the computational burden.Under the altitude of 150 km, during the atmospheric approach, the sampling rate of the MPC-based controller must be increased to be effective (>10 −2 Hz), and the objective of the optimization is the error in terms of position and velocity at TAEM.A shorter prediction horizon would make the control signal too sensitive to model uncertainties and atmospheric perturbations, resulting in a strong deviation from the trajectory and leading to an increase in the target error.

Path Planning
In this section, the results of the path planning phase are shown.A number N h = 5 of altitude ranges was chosen, to reduce the computational burden, providing a good sub-optimal solution.
The orbital parameters used in the optimization are resumed in Table 1.As a result, the guidance algorithm provides five constant values of the optimal aero-brake surface S i for the entire re-entry path, as described in Table 2.In Figure 2, the trajectory provided by the guidance system is shown in terms of altitude variation during time.The target is located at an altitude of 30 km, and its geographic position is above the Great Victoria Desert, Australia.
The main stage is the de-orbiting, whose length is 26 days, 19 h, 2 min, and 30 s.The atmospheric approach, together with the atmospheric entry phase, represents just a small piece of the entire path: the former takes 1 h, 12 min, and 19 s to go from 150 km to the altitude of 100 km; the latter takes just 6 min and 11 s to cover the last 70 km of the remaining height.
It is worth noting that the de-orbiting phase is characterized by oscillations (Figure 2), which are due to the perturbative phenomena caused by the asymmetry of the Earth's gravity potential: indeed, the satellite, while descending towards the Earth's surface, follows an orbit whose semimajor axis gradually reduces, and this results in a slightly elliptical spiral trajectory in which the position vector is no longer constant.As a matter of fact, the period of the single oscillation is comparable to the nominal orbital period (P = 5431 s).
The last portion of the trajectory, below 150 km, appears to be quasi-vertical, but a closer look (Figure 2) shows the typical ballistic re-entry path.At these heights, limiting thermal loads on the structure is of paramount importance.In Figure 3, the heating flux on the capsule for the whole descent phase is shown.As can be seen, the effect of the atmosphere is detectable from an altitude of 150 km, but becomes remarkable only below 100 km.This height marks the beginning of the atmospheric entry, in which the thermal flux on the capsule suddenly increases, reaching its peak value of about 584 kW/m 2 at an altitude of 72 km and then decreases.

Sensitivity Analysis
To provide a deeper understanding of the proposed strategy's behaviour, a campaign of numerical simulations was carried out to evaluate the performance of the proposed control system, varying its main operating parameters.
In these simulations, no model uncertainties or initial condition perturbations were taken into account, except for the sensor noise and the quantization of control signals.More specifically, the GPS error is modelled as the sum of two contributions: a sinusoidally varying error with a frequency equal to the average orbital frequency and millimetre level amplitude for position error and 10 −2 millimetre per second amplitude for velocity error and a Gaussian white noise with a standard deviation of 1.5/ √ 2 m for both horizontal and vertical error and standard deviation of 0.03/ √ 3 m/s for speed error.Those values were defined on the basis of the CubeSat kit GPSRM-1 GPS Receiver, which uses a OEM719 Multi-Frequency GNSS Receive (http://novatel.com/products/receivers/gnssgps-receiver-boards/oem719)(accessed on 27 November 2022).
Since the controlled re-entry is composed of two different phases (de-orbiting and atmospheric approach), two NMPC-based schemes were implemented: 1.
The former must guarantee that the satellite follows the reference trajectory and reaches the transition point (H = 150 km) with a low tracking error; 2.
The latter has the goal of reaching the TAEM interface with a very small position error, in order to land near the desired point.
Both controllers are based on the same model, but with different configuration parameters, i.e., the weights of the objective function, the control frequency, and the length of the prediction horizon.In the atmospheric approach, due to its short time duration compared to the deorbiting phase, the controller must be faster, requiring a smaller prediction horizon and a greater control frequency.
Another important constraint to be taken into account is the duty cycle of the actuator moving the aero-brake surface.Usually, this kind of servo cannot be operating all the time, but needs to be stopped to avoid overheating.The characteristics of the deployable aero-brake system are summarized in Table 3.As shown, the low actuation rate and the limited duty cycle suggest avoiding a too-high control frequency, which may be useless or even dangerous if not constrained, in particular, by considering the time to extend the actuator t ext , given by the ratio between the extension and the actuation rate: t ext = 0.03/0.003= 10s A 10% duty-cycle results in a minimum controller sample time of 100 s.During deorbit, a slow control frequency becomes a safety constraint to cope with the presence of perturbations and uncertainties, which could lead to strong variation in the aero-brake surface.On the other hand, in the atmospheric approach, the different objective makes the control signal smoother, enabling a higher control frequency to be used.To highlight the effects of de-orbiting MPC parameters, Figure 4 shows the position error before the atmospheric approach at the altitude of 150 km, obtained by changing the length of the prediction horizon and the controller sample time.As can be noticed, some combinations of the parameter values are not allowable: by reducing the length of the prediction horizon, the sample time must be decreased as well; otherwise, the controller is unable to bring the capsule towards the target.Moreover, looking at a prediction horizon smaller than 14,400 s (4 h), the results are quite similar and the increased computational burden is not compensated by an increased level of performance.The effect of sample time becomes more significant as the prediction horizon length increases; in this case, as expected, a higher controller frequency leads to better results.Obviously, the computational burden increases linearly with frequency; therefore, a trade-off between performance and sample time must be found with respect to the specific on-board computational power.Starting from the de-orbiting results, a sensitivity analysis on the atmospheric approach controller was conducted, as shown in Figure 5.It is worth noticing that, with a sample time of 300 s, sometimes, the capsule is unable to arrive at the destination.As expected, the MPC frequency must be higher in the last phase to guide the capsule towards the target.However, below a sample time threshold of 300 s, the results are very similar, depending on the specific values selected for T c and T p , giving the opportunity to obtain a lower computational burden and avoid higher actuator duty cycles.
In the de-orbiting phase, the weights p(i), q(i), and w(i) of the objective function can be modulated to obtain different results.Figure 6 shows the results in terms of position and velocity error obtained by considering: To reduce the dimensionality of the figure, on the abscissa, the ratio q 0 /p 0 between two weights is highlighted.The sensitivity analysis suggests that the weight on the control signal is almost useless, without a significant change in the final position error.The best ratio q 0 /p 0 is equal to unity, advising that the position and the velocity error must be balanced.

Robustness Analysis
To assess the effectiveness of the proposed MPC-based trajectory tracking technique, a campaign of numerical simulations was carried out in the presence of model uncertainties, initial state perturbations, and measurement noise.
The simulation scenarios are summarized in Table 4.A maximum initial error of −500 m on the semi-major axis combined with an error on the true anomaly of 0.1 deg was taken as the initial state perturbation.The initial velocity error is related to the position error, because of the satellite's kinetic energy: a reduction in the semi-major axis leads to a rise in the kinetic energy of the satellite, thus resulting in an initial velocity error.The model uncertainties were defined in terms of variations in the ballistic coefficient and air density values.In particular, the following was assumed: The number of time steps for the prediction horizon of the de-orbiting phase was set equal to 5 h, i.e., T p = 18,000 s.The control system sampling time, for this stage, was set equal to T c,d = 3600 s, whereas it was set equal to T c,r = 10 s for the atmospheric approach.
Table 5 shows the results in terms of position and velocity errors in the ECEF reference frame.The initial state perturbation barely affects the final result.The most important sources of errors come from the uncertainties on drag and air density, but the maximum error at TAEM is always under the prescribed objective.The worst situation is the case with a drag coefficient under the nominal value, and an interesting result is shown in sim7 and sim8, where the density uncertainty compensates the decrease of the drag coefficient, improving the final result.Figure 7 shows the capsule final position (marked with a blue cross) compared with the target location (marked with a red cross).It can be seen that, even in the worst case (highest position error), the maximum allowable error of 200 km at TAEM is always satisfied.Figures 8 and 9 show the control signals, i.e., the time variation of the aero-brake surface, for the de-orbiting and atmospheric approach, respectively.It must be noted that sim0 denotes the nominal case, i.e., absence of disturbances.
Despite the control law for the first stage of the mission (Figure 8) seeming to be very noisy, it is important to clarify that the control system has a frequency of 2.78 × 10 −4 Hz, thus leading to a variation in the control law every 3600 s for almost 27 days.The atmospheric approach (Figure 9) is characterized by a much more regular control law, though the optimal surface looks quite different from the reference one in the presence of uncertainties.In Figure 10, the actual trajectories obtained in some simulation scenarios are shown, to highlight the difference between the reference trajectory provided by the path planner and the tracked trajectory in the absence of uncertainties and noise.As can be seen in the zoomed-in detail, the trajectories are almost overlapped in all the test cases shown in the figure .Figure 11 shows the same trajectories, focusing on the atmospheric approach.The starting times are not equal, because they depend on the de-orbiting phase, achieving the altitude of 150 km (corresponding to the separation of the capsule from the main body) in different time instants.However, each trajectory reaches the target at TAEM with an acceptable tracking error.

Conclusions
In this paper, the design of a guidance and control algorithm for the re-entry of an orbiting capsule is described.The capsule is equipped with a deployable aero-brake system in order to increase the efficiency and re-usability of this kind of vehicle.The algorithm implements both off-board and on-board path planners to design the re-entry reference trajectory before the mission and locally update it in case of deviations from the prescribed orbit.
The trajectory optimization model is designed to minimize position and velocity error at TAEM, by considering several constraints: the discrete control signal of the aero-brake surface, the maximum admissible value of the heat flux, and the nonlinear dynamic model.
Once the re-entry begins, NMPC is able to modulate the aero-brake surface deployment to minimize position and speed errors.Such an NMPC uses different parameters depending on the specific phase of the re-entry mission.During the long de-orbiting phase, a greater sample time is selected to reduce the high computational burden that would result due to the long prediction horizon.In the atmospheric approach, the sample time can be decreased, because the prediction horizon is bounded by the proximity with the target, resulting in a lower computational burden.
A sensitivity analysis on the parameters of the NMPC was performed, in order to help designers in their proper selection for balancing performance and computational burden.
A numerical robustness analysis was carried out by simulating the re-entry mission in the presence of model uncertainties, initial state perturbations, and measurement noise.The position error at the target is always compliant with the maximum allowable error at TAEM in all the simulations presented.

Figure 1 .
Figure 1.Guidance and control algorithm architecture.

Figure 2 .
Figure 2. Reference trajectory provided by the path planner.The bottom axis displays the time scale in days, while the upper axis shows it in seconds.The main figure depicts the complete reference trajectory, from the de-orbit point up to the target altitude (H = 30 km); the two minor figures present the details of the de-orbiting path (on the left) and atmospheric entry trajectory (on the right).

Figure 3 .
Figure 3. Variation of the heat flux on the capsule during re-entry.The maximum value, reached at the altitude H = 71.94km, is 583.91 kW/m 2 .

Figure 4 .
Figure 4. Sensitivity analysis on prediction horizon T p and controller sample time T c .Position error at the altitude of 150 km, before atmospheric approach.

Figure 5 .
Figure 5. Sensitivity analysis on atmospheric entry MPC sample time.Position error at the altitude of 30 km (TAEM).

Figure 6 .
Figure 6.Sensitivity analysis on the objective function weights.Position error at the altitude of 30 km (TAEM).

•A
±10% uncertainty on the nominal value of the aerodynamic drag coefficient C D ; • Air density uncertainty defined by a time-dependent perturbative model (NASA/MSFC Global Reference Atmospheric Model) in which the perturbation depends on four random variables.

Figure 7 .
Figure 7. Actual and target position comparison.The target point (red cross) is located in the Great Victoria Desert, Australia.The blue cross denotes the capsule actual position for the simulation scenarios in Table5.
Figure 7. Actual and target position comparison.The target point (red cross) is located in the Great Victoria Desert, Australia.The blue cross denotes the capsule actual position for the simulation scenarios in Table5.

Figure 8 .
Figure 8. Optimized control law for the de-orbiting phase.The bottom axis displays the time scale in days, while the upper axis shows it in seconds.

Figure 9 .
Figure 9. Optimized control law for the atmospheric approach phase.The bottom axis displays the time scale in hours, while the upper axis shows it in seconds.

Figure 10 .
Figure 10.Actual and reference trajectory comparison.The bottom axis displays the time scale in days, while the upper axis shows it in seconds.The zoomed-in detail of the de-orbiting path shows the overlapping of the different trajectories.

Figure 11 .
Figure 11.Actual and reference trajectory comparison for the atmospheric approach and the atmospheric entry phases.The bottom axis displays the time scale in hours, while the upper axis shows it in seconds.

Table 2 .
Aero-brake surface deployment schedule along the re-entry path.

Institutional Review Board Statement:
Not applicable.Not applicable.Not applicable.The authors declare no conflict of interest.Nomenclature H 1 , H 2 , H 3 , H 4 , H 5 altitude boundaries of the four main re-entry phases a ECI satellite position in the Earth-Centred Inertial (ECI) reference frame V ECI satellite velocity vector in the Earth-Centred Inertial (ECI) reference frame R ECEF satellite position in the Earth-Centred Earth Fixed (ECEF) reference frame V ECEF satellite velocity vector in the Earth-Centred Earth Fixed (ECEF) reference frame µ standard gravitational parameter a J 2 , a J 3 , a J 4 perturbation accelerations related to zonal harmonics J 2 , J 3 , and J 4