Optimal planning and guidance for Solar System exploration using Electric Solar Wind Sails

Acta Astronautica 217 (2024) 116–129 Available online 1 February 2024 0094-5765/© 2024 The Authors. Published by Elsevier Ltd on behalf of IAA. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). Contents lists available at ScienceDirect Acta Astronautica journal homepage: www.elsevier.com/locate/actaastro Research paper Optimal planning and guidance for Solar System exploration using Electric Solar Wind Sails Javier Urrios, Guillermo Pacheco-Ramos∗, Rafael Vazquez Department of Aerospace Engineering, Universidad de Sevilla, Seville, Spain ARTICLE INFO Keywords: Electric Solar Wind Sails Rendezvous missions Optimal planning Guidance strategies Non-Linear Programming Model Predictive Control ABSTRACT Electric Solar Wind Sails (E-Sails) are a new type of spacecraft propellantless propulsion system that gathers its energy from solar wind protons and is potentially useful for interplanetary missions. Although optimal interplanetary trajectories have been the subject of thorough research, the substantial variability of the solar wind necessitates the adoption of active guidance strategies, an area that has received significantly less scholarly attention. This paper proposes guidance algorithms for E-Sails based on Model Predictive Control (MPC), a modern control methodology based on online re-planning of the trajectory. To this end, first, properties of E-sail time-optimal orbits are studied applying Pontryagin’s Minimum Principle, and then timeoptimal orbits for missions to Mars and Jupiter are computed via direct transcription methods. Next, solar wind perturbations are modeled, posing a challenging saturation problem due to their high variability. Guidance strategies based on Shrinking Horizon and Receding Horizon Model Predictive Control (RHMPC) are developed, analyzed and compared using Monte Carlo simulations, successfully implementing MPC to E-sail guidance. Lastly, the RHMPC strategy is successfully tested with accurate historical solar wind data from the WSA-Enlil model. 1. Introduction Since the dawn of the Space Age, the exploration of the Solar System has been one of the main drivers of space technology. However, going beyond Earth’s closest neighbors is exceedingly challenging due to the enormous energetic costs involved in orbital transfers. This limitation motivated the study of more efficient propulsion methods, such as propellantless systems. The Electric Solar Wind Sail (E-sail) falls into the broader category of propellantless propulsion systems, and in particular, belongs to those designed to extract from the Sun the energy needed for orbital transfers. The solar sails, also known as photon sails, conceived by K. Tsiolkovsky [1] and F. Tsander [2] among others, use sunlight photons as the source of energy. They work similarly to a sailing boat, hence their name. The idea is to deploy a large reflective surface around the spacecraft to reflect photons from the solar wind. By Maxwell’s Theory of Electromagnetism, a resultant force and moment (not necessarily perpendicular to the surface) are exerted on the spacecraft. Therefore, both orbit and attitude can be modified. This type of solar sails is the most common to date; LightSail 2 [3] and IKAROS [4] are examples of successful missions carried out using this technology. Looking ahead, the ACS3 [5] or Solar Cruiser [6] missions are representative of the advancements that can be achieved in the coming years. ∗Corresponding author. E-mail addresses: [email protected] (J. Urrios), [email protected] (G. Pacheco-Ramos), [email protected] (R. Vazquez). Later, in 1988, D. Andrews and R. Zubrin, proposed the concept of magnetic solar sail [7]. It works by creating a magnetic field around the sail that deflects the solar plasma, and therefore extracts its momentum. However, technical challenges in regards to the magnetic field source are still unsolved, and no real mission has used this type of sail yet. Lastly, the Electric Solar Wind Sail concept, which is the one considered in this article, was initially conceived by P. Janhunen in 2004 [8]. Inspired by the development of photonic solar wind sails, he studied the feasibility of a similar model that produced thrust by modifying the momentum of protons using an electric field generated by a mesh of wires, thus creating repulsive Coulomb forces. From there, many investigations were conducted around its design [9–11], deployment [12], applications [13], thrust modeling [14,15] and performance [16,17], among other topics. In terms of the orbital mechanics problem, techniques for computing optimal trajectories [18–21], applications to potentially viable missions [22–25] or basic algorithms for trajectory correction [26,27] are examples of research carried out. Investigations around attitude dynamics have also been conducted [28,29], some of the latest related to the development of control laws with realistically shaped tethers [30–33], the study of the effects in dynamics and control of important variables such as spin rate using FEMs [34,35], the deduction of propulsive models taking into account more complex https://doi.org/10.1016/j.actaastro.2024.01.030 Received 7 August 2023; Received in revised form 16 November 2023; Accepted 24 January 2024 Acta Astronautica 217 (2024) 116–129 117 J. Urrios et al. phenomena like coning angle dynamics [36], or the study of tethers’ oscillations using dynamic multibody models [37]. In regards to the suitability of this propellantless system, the ESTCube-2 CubeSat is planned to be launched in the second half of 2023, and will serve to demonstrate the E-sail technology. While the optimal trajectory problem has been extensively studied, the guidance problem due to solar wind uncertainty is a much more challenging problem [38,39]. Although this problem has been previously tackled by repetitively solving the trajectory indirect problem in [26], this work aims to solve this challenge using conventional control strategies based on Model Predictive Control (MPC). This concept, which originated in the late seventies, is a family of control strategies that use a system model to predict future values of state variables up to a future time horizon in terms of the past and present state and control values, and then minimizes a cost function in order to obtain the optimal future control values. Each sampling time, this process is repeated. MPC is therefore very suited for slow systems that require a robust control due to the presence of significant disturbances. In fact, it has already been considered for spacecraft missions such as rendezvous [40] or attitude control for spacecraft rendezvous [41] or E-sail attitude maneuver [32,33]. Due to the high solar wind variability, MPC is well-suited for the development of guidance strategies for the E-sail, which is crucial to its applicability, since thrust produced is directly related to solar wind. This fact is of great importance, because until now, solar wind fluctuations interfered with interplanetary spacecrafts as variations in radiation pressure or similar physical phenomena, which were relatively small in terms of magnitude. However, variations in solar wind can potentially double or halve the thrust produced by the E-sail. This dramatic variability has to be taken into account in the guidance strategy. One of the simplest ways to deal with uncertainty, and the one chosen in this work, is the consideration of worst-case scenarios or safety margins taking into account the distribution of uncertainty. In MPC, similar ideas have been used in the past, by setting the constraints to the worst possible case, see e.g. [42] or using chance constraints from which safety margins are derived [43,44]. In addition, the scope of robust design methodologies is broad and diverse, offering a multitude of avenues for enhancing guidance robustness in the face of uncertainty. Indeed, the most recent literature contains a variety of such methodologies, such as tube-based MPC [45], direct adaptive control [46], robust learning methods [47], reinforcement learning [48], stochastic optimal control [49], belief Markov decision processes and belief optimal control [50], or robust epistemic uncertainty optimization [50]. These approaches are particularly pertinent when facing indeterminate models of uncertainty, wherein the identification and characterization of uncertainty become essential [42,45,51]. Each of these methodologies can contribute in different ways to the robustness of guidance systems, expanding the toolbox available for dealing with uncertainty in optimal trajectory design. The main contributions of this work can be summarized as follows. Firstly, the performance of the guided trajectory problem formulation is evaluated under solar wind uncertainty, employing both shrinking and receding strategies for missions to Mars and Jupiter. Beyond feedback’s own ability to reject disturbances, uncertainties are dealt with by using worst-case scenarios and safety margins (taking into account the distribution of solar wind variability), which inevitably leads to sub-optimal solutions. The proposed strategy is applicable to various planetary or interplanetary missions, including destinations within the inner planets. However, Mars and Jupiter are chosen based on the superior performance attributed to the E-sail when operating farther from the Sun compared to other propellantless systems relying on solar wind, owing to its lower decay [16]. Furthermore, a study of the influence of the safety margin for nominal acceleration on trajectory recalculations and its impact on the success rate through Monte Carlo analysis is presented. For the receding horizon strategy, various objective functions are proposed and compared. Finally, the viability Fig. 1. Conceptual E-sail thrust model illustration. Fig. 2. Variation with respect to 𝛼𝑛of 𝛾(in blue) and 𝛼(in orange). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) of receding horizon MPC is demonstrated by considering time-varying solar wind variations derived from accurate models. In this work, a two dimensional problem is considered for simplicity and circular orbits are assumed for planet’s orbits [16,39]. Sun-based canonical units are used in the formulation. The cited hypothesis allows to keep the focus on the guidance problem and the solar wind variability by eliminating the influence of the planets’ orbits but do not diminish the value of the main contribution of this work. The consideration of a more detailed 3D problem with more realistic orbits has already been assessed for E-sail and its feasibility demonstrated [52]. This paper is organized as follows. Section 2introduces the models used for orbital motion, E-sail thrust and solar wind uncertainty. Next, Section 3presents the optimal trajectory planning under nominal solar wind for interplanetary missions to Mars and Jupiter. Subsequently, in Section 4, the formulation of the Shrinking Horizon and Receding Horizon MPC guidance strategies is discussed, taking into account the effects of solar wind uncertainties. Additionally, the results of the Monte Carlo analysis carried out for both approaches are provided. In Section 5, a successful mission to Mars is computed using actual historical solar wind data and employing the Receding Horizon MPC strategy. Finally, the conclusions drawn from the obtained results are presented in Section 6and the future works are identified. 2. Models for E-Sail trajectories In this section a model for the E-sail thrust is introduced and the equations of motion are presented. Next, solar wind uncertainty is introduced and modeled as a probability distribution. Additionally, model limitations and assumptions are also discussed. For convenience, canonical units are used, with the acronym AU meaning Astronomical Unit, and UV Units of Velocity (that is, AU divided by the canonical Unit of Time UT). In this sense, 1AU ≈ 1.496 ⋅ 1011 m, 1UT ≈ 58.1324 days and 1UV ≈ 29.785 km/s. 2.1. E-Sail thrust model The E-sail thrust model used for this work is the one presented in [14], which is developed from a geometrical standpoint, and Acta Astronautica 217 (2024) 116–129 118 J. Urrios et al. Fig. 3. Gamma distribution of 𝑝⊕(red line) compared with real data gathered from NASA (bars). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) proposes simple yet precise expressions of radial and tangential accelerations produced by the E-sail. It should be noted the thrust model is reliable in the solar wind near 1 AU and for typical E-sail parameters. Outside of those conditions, the accuracy of the model may be diminished. The radial and tangential accelerations are given by 𝑎𝑟=𝜏𝑎𝑐 2(𝑟⊕ 𝑟)(1 + cos2𝛼𝑛), 𝑎𝜃=𝜏𝑎𝑐 2(𝑟⊕ 𝑟)cos 𝛼𝑛sin 𝛼𝑛. (1) Note that the acceleration produced by the E-sail is inversely proportional to the distance from the Sun, 𝑟. In this model, 𝜏∈[0,1], which is a control variable, represents the throttle level of the E-sail, with 0 meaning that no thrust is produced and 1 representing maximum thrust. The value 𝑎𝑐is the spacecraft characteristic acceleration (i.e. the maximum value of acceleration measured at 1AU from the Sun, 𝑟⊕). The sail pitch angle, which is denoted 𝛼𝑛, has an influence not only in the acceleration direction, but also in the modulus. A concept visualization of this thrust model is presented in Fig. 1. In this study, it is assumed that the throttle level and the sail pitch angle can be instantaneously controlled at will, i.e., they are the control variables. In reality, controlling the sail pitch angle presents an attitude problem that is also being addressed (see [30,34]), although it occurs on a much shorter timescale. Defining 𝑎as the acceleration modulus computed from the components described in Eq. (1), the dimensionless propulsive acceleration modulus, considering 𝜏= 1, is defined as 𝛾=𝑎∕(𝑎𝑐(𝑟⊕∕𝑟))and is represented Fig. 2 as a function of the sail angle 𝛼𝑛. The thrust cone angle 𝛼, is also depicted. There is a maximum value of 𝛼𝑛, which recent studies approximate to be in the range of 60◦to 70◦[11], since going beyond that limit implies dealing with non-desirable instabilities. However, looking at the evolution of 𝛾, it can be seen that, assuming an optimal transfer, there is no apparent reason for going beyond the sail pitch angle that yields maximum thrust cone angle 𝛼𝑚𝑎𝑥 (marked in Fig. 2), since the non-dimensional acceleration 𝛾is strictly decreasing with 𝛼𝑛. The results provided in Section 3.2 corroborate that, for optimal mission planning, the sail angle value remains within the feasible range. 2.2. Equations of motion First, the orbital equation of motion including the E-sail acceleration,  𝒓= − 𝜇 𝑟2 𝒓 𝑟+𝒂+𝒂𝑝(2) where 𝜇= 1.327124 ⋅1011 km3/s2is the Standard Gravitational Parameter of the Sun, 𝒓denotes the position vector in the Heliocentric Inertial reference frame (HCI), 𝒂accounts for the E-sail acceleration and 𝒂𝑝 represents the acceleration caused by disturbances, such as third-body perturbations. The equations of motion are obtained by conveniently expressing the equation of motion as a set of first order differential equations expressed in polar coordinates, yielding 𝑟 =𝑣𝑟,  𝜃=𝑣𝜃 𝑟, 𝑣𝑟=𝑣2 𝜃 𝑟−𝜇 𝑟2+𝜏𝑎𝑐 2(𝑟⊕ 𝑟)(1 + cos2𝛼𝑛), 𝑣𝜃= − 𝑣𝑟𝑣𝜃 𝑟+𝜏𝑎𝑐 2(𝑟⊕ 𝑟)cos 𝛼𝑛sin 𝛼𝑛. (3) Note that perturbations are omitted, which will be taken care of by the closed loop guidance strategy. 2.3. Solar wind uncertainty model According to [53], the force per unit length exerted by each one of the tethers that make up the E-sail follows the expression 𝑓𝑙= 0.18 max(0, 𝑉0−𝑉1)√𝜖0𝑝⊕,(4) where 𝑉0denotes the tether voltage, 𝑉1the solar wind electric potential, 𝜖0is the vacuum permittivity and 𝑝⊕denotes the dynamic pressure of the Sun at 1 AU distance from it, all expressed in the International System of Units. Since typically 𝑉0≫ 𝑉1, it can be assumed that the force per unit length generated by a tether is proportional to the square root of the dynamic pressure. Moreover, the total force generated by a tether can be computed knowing the force per unit length and applying a correction factor 𝑘1(𝑓𝑡=𝑘1𝑓𝑙𝑙) which depends on the geometry. Similarly, the total force is obtained by multiplying the force of an individual tether by the number of tethers 𝑛𝑡, and applying another geometric correction 𝑘2(𝐹=𝑘2𝑛𝑡𝑓𝑡). Therefore, the total force is proportional to the force per unit length, and thus acceleration is proportional to the squared root of the dynamic pressure (i.e. 𝑎∝𝑝1∕2 ⊕). Thus, the modeling of dynamic pressure uncertainties is crucial, since it directly impacts the performance of the E-sail. Predicting the dynamic pressure of the Sun is remarkably difficult. The intrinsic chaotic and turbulent processes happening in its interior make the fluctuations of solar wind rather unpredictable [54]. Fig. 3 represents a histogram of the dynamic pressure obtained from data gathered from the NASA OMNIWeb Data Explorer, as proposed in [38]. The mean dynamic pressure at 1 AU can be estimated to be around 𝑝⊕≈ 2 nPa, with an standard deviation of approximately 𝜎𝑝⊕≈ 1.56 nPa. Following [38], the distribution function of the solar wind dynamic pressure at 1AU is approximated by a gamma distribution. A comparative of both the real distribution and gamma approximation can be observed again in Fig. 3. The probability density function of 𝑝⊕is given by 𝑓(𝑝⊕) = 𝛽−𝛼 𝛤(𝛼)𝑝𝛼−1 ⊕exp(−𝑝⊕∕𝛽),(5) where the coefficients 𝛼= 1.6437,𝛽= 1.2168, are considered to adjust the mean and standard deviation. As shown in Fig. 3, the gamma probability density function (PDF) tends to overestimate the likelihood of having values of 𝑝⊕smaller (or slightly greater) than the mean value 𝑝⊕. Simultaneously, it underestimates the peak near 𝑝⊕and the values in the tail. This leads to a certain level of conservatism in the considered solar pressure distribution compared to the actual one. It is worth noting that the guidance strategies proposed in this study could potentially be applied to other models of solar pressure PDFs, and better results could be achieved if a more precise uncertainty model is considered [55]. Finally, since the real data for solar wind dynamic pressure available corresponds to 1 AU, it should be acknowledged that the representativeness of the uncertainty model used is anticipated to diminish with increasing distance to the Sun. Acta Astronautica 217 (2024) 116–129 119 J. Urrios et al. Fig. 4. Saturation probability with respect to control margin. 3. Optimal mission planning under nominal solar wind The optimal orbit control problem is tackled in the current Section. Time-optimal orbits for rendezvous missions to Mars and Jupiter are obtained assuming nominal solar dynamic pressure. However, the solar wind uncertainty described in Section 2.3, is considered during mission planning by applying a control margin of 40%, equivalent to a saturation probability of 20%, to the nominal characteristic acceleration, 𝑎𝑐= 0.6 m∕s2. This results in the characteristic acceleration 𝑎𝑐= 0.36,m∕s2used. The relation between saturation probability and control margin is shown in Fig. 4. Since the time-dependency of the planets’ position is not explicitly considered, this problem can be seen as an orbit-to-orbit transfer. However, due to the rotational symmetry of the problem (which comes from considering circular orbits), the orbit-to- orbit transfer becomes a planet-to-planet transfer simply by choosing the appropriate launch date. First, the Optimal Control Problem (OCP) is formulated in Section 3.1. Then, it is solved by both an indirect (Section 3.2) and a direct (Section 3.3) approach. On the one hand, indirect methods apply optimality conditions first, which makes it possible to obtain local laws for the control variables converting the problem into a Two- Point Boundary Value Problem. On the other hand, direct transcription discretizes the problem in time, converting the set of optimization variables into a finite set of parameters, that is, a Non-Linear Programming (NLP) problem. Additionally, direct methods tend to be more stable and less dependent on initial conditions compared to indirect approaches. Therefore, the optimal orbits are first computed solving the NLP problem. Finally, the results obtained with both methodologies are compared in Section 3.4, in order to validate the solution and implementation of the direct method. This is particularly relevant as the same method is employed for the guidance phase described in Section 4. 3.1. Optimal control problem formulation The classical formulation for an OCP involves a cost function of the type 𝐽[𝒙(⋅),𝒖(⋅), 𝑡0, 𝑡𝑓]=𝜑[𝒙(𝑡0), 𝑡0,𝒙(𝑡𝑓), 𝑡𝑓]+∫𝑡𝑓 𝑡0 [𝒙(𝑡),𝒖(𝑡), 𝑡]d𝑡(6) subject to the following state equations, path constraints and boundary conditions: 𝒇[𝒙(𝑡),𝒖(𝑡), 𝑡]= 𝒙,(7) 𝒉[𝒙(𝑡),𝒖(𝑡), 𝑡]≤0,(8) 𝒈[𝒙(𝑡0), 𝑡0,𝒙(𝑡𝑓), 𝑡𝑓]= 0.(9) The bold notation is used to denote a vector. For this particular case, 𝒙= {𝑟, 𝜃, 𝑣𝑟, 𝑣𝜃}and 𝒖= {𝜏, 𝛼𝑛}. In this manner, the minimum time OCP can be expressed as min 𝜏,𝛼𝑛 𝑡𝑓 s.t. 𝑟 =𝑣𝑟  𝜃=𝑣𝜃 𝑟 𝑣𝑟=𝑣2 𝜃 𝑟−𝜇 𝑟2+𝜏𝑎𝑐 2(𝑟⊕ 𝑟)(1+cos2𝛼𝑛) 𝑣𝜃= − 𝑣𝑟𝑣𝜃 𝑟+𝜏𝑎𝑐 2(𝑟⊕ 𝑟)cos 𝛼𝑛sin 𝛼𝑛 ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟ state equations 0≤𝜏≤1 −𝛼𝑚𝑎𝑥 𝑛≤𝛼𝑛≤𝛼𝑚𝑎𝑥 𝑛 ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟ path constraints 𝑟(0) = 𝑟⊕ 𝜃(0) = 0 𝑣𝑟(0) = 0 𝑣𝜃(0) = 𝑣⊕ 𝑟(𝑡𝑓) = 𝑟𝑝 𝑣𝑟(𝑡𝑓) = 0 𝑣𝜃(𝑡𝑓) = 𝑣𝑝, ⏟⏞⏞⏞⏞⏟⏞⏞⏞⏞⏟ boundary conditions (10) where, without loss of generality, it has been imposed that 𝑡0= 0. The terms 𝑟𝑝and 𝑣𝑝, with 𝑝=♂,X, represent the distance from the Sun and the velocity of Mars and Jupiter respectively. Using astronomical units: 𝜇= 1 AU UV2,𝑟⊕= 1 AU, 𝑟♂= 1.5237 AU, 𝑟X= 5.2043 AU, 𝑣♂= 0.8101 UV, 𝑣X= 0.4384 UV. Lastly, the value of 𝛼𝑛is bounded by the previously introduced value 𝛼𝑚𝑎𝑥 𝑛= 70◦. The characteristic acceleration chosen for optimal planning is denoted as 𝑎𝑐. Note that the path constraints only affect the control variables. These constraints can be simply represented as 𝒖∈, where = {𝒖∈R2∶ 0 ≤𝜏≤ 1,−𝛼𝑚𝑎𝑥 𝑛≤𝛼𝑛≤𝛼𝑚𝑎𝑥 𝑛}is the admissible control set. Notice that in (10), no restriction is imposed on 𝜃𝑓. That is because the launch date becomes ‘‘free to choose’’ under the assumed orbits of the planets (coplanar and circular) due to rotational symmetry. Once the problem is solved, the launch date can be adjusted to effectively meet the corresponding planet when the spacecraft reaches its orbit. As a final remark, it should be acknowledged that due to the non-convexity of the problem constraints, the global nature of the optimal solutions found cannot be ensured, and the obtained solutions may be regarded as local optima. 3.2. Indirect approach The difficulty of solving OCPs lies in the fact that input variables are functionals rather than discrete quantities, which makes it rather different compared to the standard parametric optimization problem. However, local laws that control parameters need to follow for the solution to be optimal can be found using indirect methods, and thus they transform the OCP into a Two Point Boundary Value Problem (TPBVP), which can then be solved numerically, although reaching convergence is not trivial. Indirect approaches have been previously used to characterize the performance of the E-sail [16,19,22,23,26]. To proceed, the state variables are expanded to the costate space, obtaining the set 𝝀= {𝜆𝑟, 𝜆𝜃, 𝜆𝑣𝑟, 𝜆𝑣𝜃}, which satisfies the Euler– Lagrange equations  𝜆𝑖= − 𝜕𝐻 𝜕𝑖 , with 𝑖=𝑟, 𝜃, 𝑣𝑟, 𝑣𝜃. The variable 𝐻is the Hamiltonian of the system, defined as 𝐻=+𝝀⊺𝒇, where is the integrand in (6). The next step is to obtain the optimal control laws 𝒖=𝒖(𝒙,𝝀, 𝑡), which can be done by applying Pontryagin’s Minimum Principle, that establishes a necessary condition (although not sufficient) for 𝒖to be optimal. Pontryagin’s Minimum Principle states that for a cost function 𝐽to be minimized, the control parameters must minimize the Hamiltonian at every instant of time, i.e., 𝒖⋆=argmin𝒖∈𝐻∀𝑡∈ [0, 𝑡𝑓], where ⋆denotes optimality. For this particular system, the Hamiltonian is given by 𝐻=𝜆𝑟𝑣𝑟+𝜆𝜃 𝑣𝜃 𝑟+𝜆𝑣𝑟(𝑣2 𝜃 𝑟−𝜇 𝑟2)−𝜆𝑣𝜃 𝑣𝑟𝑣𝜃 𝑟 +𝜏𝑎𝑐 2(𝑟⊕ 𝑟)[𝜆𝑣𝑟(1 + cos2𝛼𝑛)+𝜆𝑣𝜃cos 𝛼𝑛sin 𝛼𝑛]. (11) Acta Astronautica 217 (2024) 116–129 120 J. Urrios et al. Fig. 5. Evolution of variables (left) and trajectory (right) in a transfer orbit to Jupiter using a direct approach. If the optimal control is not located on the boundary of the admissible set, that is, 𝒖⋆∉𝜕, the minimum value of a Hamiltonian must be a critical point. Therefore, 𝜕𝐻 𝜕𝛼𝑛 is computed, which yields −𝜆𝑣𝑟sin(2𝛼𝑛) + 𝜆𝑣𝜃cos(2𝛼𝑛)=0,(12) obtaining two solutions, 𝛼𝑛1=1 2arctan 2 (𝜆𝑣𝜃, 𝜆𝑣𝑟)and 𝛼𝑛2=1 2arctan 2 (−𝜆𝑣𝜃,−𝜆𝑣𝑟), where arctan 2 is the two-argument inverse tangent function. Computing 𝜕2𝐻 𝜕𝛼2 𝑛 it is easy to see that the solution that corresponds to a minimum is 𝛼𝑛2. However, that point could be outside of the admissible region, so the complete expression is 𝛼⋆ 𝑛={𝛼𝑛2if |𝛼𝑛2|≤𝛼𝑚𝑎𝑥 𝑛, 𝛼𝑚𝑎𝑥 𝑛⋅sgn (𝛼𝑛2)if |𝛼𝑛2|> 𝛼𝑚𝑎𝑥 𝑛.(13) Regarding 𝜏, a linear dependency is seen in 𝐻. Thus, to minimize the Hamiltonian, the best possible strategy for 𝜏is to have minimum value when it is multiplied by a positive expression, and maximum in the opposite case. 𝜏⋆={0 if 𝜆𝑣𝑟(1 + cos2𝛼⋆ 𝑛)+𝜆𝑣𝜃cos 𝛼⋆ 𝑛sin 𝛼⋆ 𝑛>0, 1 if 𝜆𝑣𝑟(1 + cos2𝛼⋆ 𝑛)+𝜆𝑣𝜃cos 𝛼⋆ 𝑛sin 𝛼⋆ 𝑛≤0.(14) It can be computed that 𝜏⋆switches to 0 whenever 𝛼⋆ 𝑛>arccos (1∕√3) ≈ 54.7◦, effectively turning the E-sail off. This value coincides with the value of 𝛼𝑛that yields the maximum thrust cone angle 𝛼𝑚𝑎𝑥 = arctan (1∕2√2)so, as conjectured in Section 2.1, there is no apparent reason for the sail pitch angle to go beyond 54.7◦. Therefore, since 𝛼𝑚𝑎𝑥 𝑛>arccos (1∕√3), saturation never takes place in the sail pitch angle and the optimal control law simply becomes 𝛼⋆ 𝑛=1 2arctan2 (−𝜆𝑣𝜃,−𝜆𝑣𝑟), 𝜏⋆=⎧ ⎪ ⎨ ⎪ ⎩ 0if |𝛼𝑛|>arccos (1∕√3), 1if |𝛼𝑛|≤arccos (1∕√3). (15) What is known as bang–bang control appears in 𝜏due to its linear dependency with 𝐻, meaning that for a solution to be at least locally optimal, the throttle level must be at all times deactivated or at its maximum. Similar results are obtained in [16], although a slightly different thrust model is used. The full TPBVP is formulated below: 𝑟 =𝑣𝑟,  𝜃=𝑣𝜃 𝑟, 𝑣𝑟=𝑣2 𝜃 𝑟−𝜇 𝑟2+𝜏⋆𝑎𝑐 2(𝑟⊕ 𝑟)(1 + cos2𝛼⋆ 𝑛), 𝑣𝜃= − 𝑣𝑟𝑣𝜃 𝑟+𝜏⋆𝑎𝑐 2(𝑟⊕ 𝑟)cos 𝛼⋆ 𝑛sin 𝛼⋆ 𝑛,  𝜆𝑟=𝜆𝜃𝑣𝜃 𝑟2+𝜆𝑣𝑟(𝑣2 𝜃 𝑟2−2𝜇 𝑟3)−𝜆𝑣𝜃 𝑣𝑟𝑣𝜃 𝑟2+𝜏⋆𝑎𝑐 2𝑟(𝑟⊕ 𝑟) ×[𝜆𝑣𝑟(1 + cos2𝛼⋆ 𝑛)+𝜆𝑣𝜃cos 𝛼⋆ 𝑛sin 𝛼⋆ 𝑛],  𝜆𝜃= 0,  𝜆𝑣𝑟= − 𝜆𝑟+𝜆𝑣𝜃𝑣𝜃 𝑟,  𝜆𝑣𝜃= − 𝜆𝜃 𝑟−2𝜆𝑣𝑟𝑣𝜃 𝑟+𝜆𝑣𝜃𝑣𝑟 𝑟, (16) where 𝛼⋆ 𝑛and 𝜏⋆are now optimal and given by (15). The system is subject to the following boundary conditions 𝑟(0)=𝑟⊕ 𝜃(0)= 0 𝑣𝑟(0)= 0 𝑣𝜃(0)=𝑣⊕ 𝑟(𝑡𝑓)=𝑟𝑝 𝜆𝜃(𝑡𝑓)= 0 𝑣𝑟(𝑡𝑓)= 0 𝑣𝜃(𝑡𝑓)=𝑣𝑝 𝐻(𝑡𝑓)= −1. (17) This TPBVP is solved by implementing a shooting algorithm in the full state-costate system that also keeps track of the precise switching times of 𝜏for numerical stability, similar to what is done in [56]. Several initial guesses are tried until a feasible solution is obtained. 3.3. Direct approach The aim of this section is to convert this OCP into a Non-Linear Programming (NLP) problem, i.e. an optimization problem with a cost function subject to a set of algebraic equalities or inequalities. Thus, time is divided into 𝑁intervals 𝒕=[𝑡0, 𝑡1,…, 𝑡𝑁]⊺. This solving method is known for being considerably more robust than indirect approaches, which highly depend on the initial guess. It has already been used for orbit optimization in low thrust spacecraft [57]. In this discretization, state variables are defined for each instant of time and control variables are assumed constant during each time Acta Astronautica 217 (2024) 116–129 121 J. Urrios et al. Fig. 6. Evolution of variables (left) and trajectory (right) in a transfer orbit to Mars. Comparison between the direct and indirect solution. interval, resulting in 4(𝑁+ 1)state variables and 2𝑁control variables. Taking into account that 𝑡𝑓is also a decision variable, and dividing the transfer in 𝑁= 500 intervals, the total number of NLP variables is 3005. Note that the value of N is relatively high because it considers not only the intervals for discretizing control but also determines the number of intervals used for discretizing the dynamic equations, as the direct transcription approach is employed. Hence, a resolution higher than the one proposed in [58] is required. Representing the values of the state variables at the instant 𝑡𝑘 and of the control parameters in the interval [𝑡𝑘, 𝑡𝑘+1]as 𝒙𝑘and 𝒖𝑘 respectively, dynamic constraints of the type 𝒙𝑘+1 =𝒇𝒅(𝑡𝑘,𝒙𝑘,𝒖𝑘)can be constructed by applying some integration scheme. A fourth-order Runge–Kutta is used in this work (details are skipped for simplicity). The boundary conditions can be added by setting 𝒙0and 𝒙𝑁to the appropriate values. The path constraints 0≤𝜏𝑘≤1and −𝛼𝑚𝑎𝑥 𝑛≤𝛼𝑛𝑘≤ 𝛼𝑚𝑎𝑥 𝑛are set for each instant of time. With all these considerations, there are 4𝑁+ 7 equality and 4𝑁inequality constraints. To solve the NLP problem, the algorithm Interior Point OPTimizer (IPOPT) [59] is used by means of the CasADi Toolbox for MATLAB [60]. As the problem constraints are not convex, one cannot expect IPOPT to arrive to the global optimal solution, but rather a local one. 3.4. Optimal planning results This section presents and compares the optimal trajectories obtained with both direct and indirect methods discussed in Sections 3.2 and 3.3. As already mentioned, a locally optimal solution is expected instead of a global, because of the non-convexity of the problem. However, indirect and direct methods can be compared to see if the solutions differ by much, giving insight into the quality of the solution, since both methods reach (local) optimality by means of different paths. It is useful to remember that, as demonstrated using Pontryagin’s Minimum Principle, a bang–bang control in 𝜏is expected. The results obtained for missions to Mars and Jupiter are shown. A typical value of 𝑎𝑐= 0.36 mm/s2has been used for the characteristic acceleration. It should be noted that this value is lower than the maximum acceleration considered in other works [26,61]. This discrepancy can be attributed to a combination of factors. On one hand, anticipating that a higher-than-usual payload fraction might be necessary for the interplanetary missions under consideration, a conservative value for the nominal characteristic acceleration is selected, being 𝑎𝑐= 0.6 m∕s2. On the other hand, taking into account the PDF distribution presented in Section 2.3, a control margin of 40%, equivalent to a probability of saturation of 20%, has been applied to the nominal characteristic acceleration 𝑎𝑐selected. Consequently, the acceleration 𝑎𝑐= 0.36 m∕s2 considered for mission planning is obtained. In Fig. 5, the resulting orbit to Jupiter using direct methods is shown, along with the evolution of state variables, control inputs, semi-major axis 𝑎and eccentricity 𝑒. Bang–bang structure in 𝜏can be observed. The intervals in which thrust is off are called coasting arcs. Moreover, a tendency to return to distances of ∼ 1AU can also be seen, and this can be easily be explained by the fact that the closer the spacecraft is from the Sun, the higher the acceleration produced by the E-sail. The total duration of the mission is around 15 years. A comparison between the direct and indirect approaches for the Mars transfer orbit is included in Fig. 6. As can be seen, the solutions obtained by both methods coincide, reaching Mars orbit in 1007 days. This does not guarantee that the solution is global, but helps validate both planning methods. More specifically, it helps justify the use of the direct method, since it can produce good results while being more reliable. Just for validation purposes, the transfer to Mars is also computed for 𝑎𝑐= 1.00 mm/s2, resulting in 502 days, which is in line with typical transfer times reported in [26,61]. Due to its higher stability, and given that the solutions obtained are similar between both methods, the direct approach is used for the guidance problem. Moreover, this direct transcription scheme makes the model predictive controller easier to implement. The values obtained for final time and polar angle are denoted as 𝑡𝑛𝑜𝑚 𝑓,𝑗 and 𝜃𝑛𝑜𝑚 𝑓,𝑗 (for 𝑗=♂,X), and are used in Section 4. 4. Guidance strategies considering solar wind uncertainties In this section the Sun variability is introduced along with its implications on saturation phenomena and control margins. Two approaches based on Model Predictive Control (MPC) are implemented, involving shrinking and receding time horizons. The selection of the propulsion strategy for control saturation periods and the impact of the cost function definition are studied. Based on the results obtained, the most promising strategy and options are identified. 4.1. The effect of solar wind perturbations Hereinafter, the nominal characteristic acceleration 𝑎𝑐is defined as the acceleration of the spacecraft at 1 AU with nominal wind at maximum tether voltage, whereas the characteristic acceleration 𝑎𝑐 considers a more realistic dynamic pressure in which solar fluctuation is Acta Astronautica 217 (2024) 116–129 122 J. Urrios et al. considered according to the model in Section 2.3. Similarly, 𝑝⊕denotes the dynamic pressure and 𝑝⊕its mean. Recalling (4), it is seen that force per unit length (and therefore the characteristic acceleration) is proportional to the tether grid voltage and to the square root of the dynamic pressure. Also, as seen in Section 2.3, the dynamic pressure randomly changes following a gamma distribution. Assuming that tethers can change their voltage instantaneously, and that the voltage change immediately translates into a change in the Coulomb’s force generated by the E-sail (such a process’ timescale is in the order of minutes, while the duration of interplanetary mission is in the order of years), one can compensate a change in 𝑝⊕instantaneously with an appropriate modification of 𝑉0 in order to maintain the needed spacecraft’s acceleration. However, a problem arises when the dynamic pressure becomes so low that, not even with maximum voltage, the necessary acceleration can be obtained. Therefore, solar wind uncertainty may result in a saturation problem, preventing the optimal trajectory tracking and manifesting the necessity of a guidance strategy. Since the characteristic acceleration of the spacecraft becomes variable, it is no longer useful to use the throttle level as a control input. Instead, an effective acceleration 𝜅is defined, which can be expressed in terms of the throttle level, the nominal characteristic acceleration, and the dynamic pressure ratio as 𝜅=𝜏𝑎𝑐=𝜏 𝑎𝑐√𝑝⊕ 𝑝⊕ .(18) If the same characteristic acceleration used in the optimal mission planning is utilized for guidance, saturation occurs each time 𝑝⊕< 𝑝⊕. Therefore, it is convenient to define a control margin to reduce the probability of saturation. In fact, such a margin has already been considered in Section 3.4. Indeed, for the guidance strategy, a typical value of the characteristic acceleration of 𝑎𝑐= 0.6mm/s2is used [11], meaning that a control margin of 40% was left when computing the optimal planning, in which the characteristic acceleration was set to 𝑎𝑐= 0.36 mm/s2. For reference, according to Fig. 4 which represents the saturation probability for a given control margin, it can be observed that not leaving control margin means that the probability of saturation is about 60%, while considering a 40% control margin (that is, 𝑎𝑐= 0.6𝑎𝑐) reduces the probability to around 20%. More control margin could be considered, but as a trade-off mission times get progressively larger. 4.2. Shrinking Horizon MPC strategy Since the dynamics of the system are known, one could predict the evolution of the spacecraft’s orbit into the future given the values of the control variables, and assuming the dynamic pressure is known during that prediction horizon (it is indeed assumed constant). Then, some error measurement can be used to quantify how far the spacecraft has steered away from the reference trajectory. Of course, the dynamic pressure varies in the actual trajectory as opposed to the constant behavior assumed when predicting into the future, but this method becomes more exact as the trajectory recalculation rate increases. Using such orbit prediction, it is possible to formulate an OCP at a given mission time. The way in which the problem is posed here is by adding final constraints on the arrival time 𝑡𝑛𝑜𝑚 𝑓,𝑝 and phase angle 𝜃𝑛𝑜𝑚 𝑓,𝑝 obtained from the nominal planned orbit. Additionally, the modulus of the final velocity is fixed, while the cost function tries to minimize the square of final radial velocity. Even if such relaxation results in the velocity not exactly being that of the reached body, this small error can be handled by other specific propulsion systems to guarantee the rendez-vous. The aim is to obtain the values for the control variables that minimizes such cost function, and the feedback loop is closed by recomputing the optimal controls with a certain frequency. The term ‘‘Shrinking Horizon’’ Model Predictive Control (SHMPC) comes from the fact that the arrival time is fixed, and the integration of the trajectory into the future is always performed up to that time instant, which implies that the remaining time keeps shrinking, as illustrated in Fig. 10(a). Under the previous conditions and considering 𝑁𝑟,𝑝 recalculations, the successive OCPs can be formulated as min 𝜅,𝛼𝑛 𝑣2 𝑟(𝑡𝑓) s.t. 𝑟 =𝑣𝑟  𝜃=𝑣𝜃 𝑟 𝑣𝑟=𝑣2 𝜃 𝑟−𝜇 𝑟2+𝜅 2(𝑟⊕ 𝑟)(1+cos2𝛼𝑛) 𝑣𝜃= − 𝑣𝑟𝑣𝜃 𝑟+𝜅 2(𝑟⊕ 𝑟)cos 𝛼𝑛sin 𝛼𝑛 ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟ state equations 0≤𝜅≤𝜅𝑚𝑎𝑥 −𝛼𝑚𝑎𝑥 𝑛≤𝛼𝑛≤𝛼𝑚𝑎𝑥 𝑛 ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟ path constraints 𝑟(𝑡𝑖) = 𝑟𝑖 𝜃(𝑡𝑖) =  𝜃𝑖 𝑣𝑟(𝑡𝑖) = 𝑣𝑟𝑖 𝑣𝜃(𝑡𝑖) = 𝑣𝜃𝑖 𝑟(𝑡𝑓) = 𝑟𝑝 𝑣2(𝑡𝑓) = 𝑣2 𝑝 𝑡𝑓=𝑡𝑛𝑜𝑚 𝑓,𝑝 𝜃(𝑡𝑓) = 𝜃𝑛𝑜𝑚 𝑓,𝑝 , ⏟⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏟ boundary conditions (19) where for 𝑖= [1,…, 𝑁𝑟,𝑝], the solutions 𝑟𝑖(𝑡),𝜃𝑖(𝑡),𝑣𝑟𝑖(𝑡)and 𝑣𝜃𝑖(𝑡)define the time optimal partial trajectories for 𝑡∈[𝑡𝑖, 𝑡𝑓]. Additionally, the final time 𝑡𝑓and the final polar angle 𝜃(𝑡𝑓)are now fixed and equal to 𝑡𝑛𝑜𝑚 𝑓,𝑝 and 𝜃𝑛𝑜𝑚 𝑓,𝑝 respectively, values obtained for the optimal planning presented in Section 3. Similarly, the final radius 𝑟(𝑡𝑓)is also fixed to 𝑟𝑝, thus ensuring the spacecraft reaches the planet. Finally, the value of 𝑡𝑖refers to the time at which the orbit is recalculated, an instant at which the state variables 𝑟𝑖, 𝜃𝑖,𝑣𝑟𝑖and 𝑣𝜃𝑖are assumed to be known thanks to the navigation system. It is notable that they generally will differ from the optimal solution associated to the previous recalculation 𝑟𝑖−1(𝑡𝑖),𝜃𝑖−1(𝑡𝑖),𝑣𝑟𝑖−1 (𝑡𝑖)and 𝑣𝜃𝑖−1 (𝑡𝑖). In addition 𝑝=♂,Xdepends on whether the spacecraft aims to reach Mars or Jupiter, respectively. The value of the maximum effective acceleration is 𝜅𝑚𝑎𝑥 =𝑎𝑐= 𝑎𝑐√𝑝⊕∕𝑝⊕, which changes over time. However, it must be known in order to solve the OCP. Thus, some assumption has to be made. In this work, where the trajectory is recomputed a total of 𝑁𝑟,♂= 20 times for Mars and 𝑁𝑟,X= 50 times for Jupiter (so recalculation is performed every few months, hereinafter the instants at which recalculation is performed are noted 𝑡𝑖for 𝑖= 1,…, 𝑁𝑟,♂or𝑁𝑟,X) dynamic pressure is assumed constant and equal to the value measured at 𝑡𝑖for the next recalculation interval. Thus, 𝜅𝑚𝑎𝑥,𝑖 is indeed equal to 𝑎𝑐at the current time interval, but it may be convenient to limit it when predicting the trajectory in future intervals, to encourage the spacecraft to use as much acceleration as possible at the current time, since the problem becomes stiffer when approaching the target. Therefore, four distinct maximum values for effective acceleration in future time intervals (denoted as 𝜅𝑚𝑎𝑥,+) are considered: (1) 𝜅𝑚𝑎𝑥,+= 0.6𝑎𝑐= 0.36 mm/s2(that is, a 40% reduction from the real nominal characteristic acceleration, i.e., the value used in optimal planning), (2) 𝜅𝑚𝑎𝑥,+= 0.8𝑎𝑐= 0.48 mm/s2(a 20% reduction, i.e., the mean between the real value and the used in orbit planning), (3) 𝜅𝑚𝑎𝑥,+=𝑎𝑐= 0.6mm/s2(no reduction), and (4) a linearly increasing value of 𝜅𝑚𝑎𝑥,+between 20% reduction and its full nominal value depending on the proximity to the end of the mission, that, is, 𝜅𝑚𝑎𝑥,+= [0.8 + 0.2 (𝑖− 1)∕(𝑁𝑟,𝑝 − 2)] 𝑎𝑐, where 𝑖is the current recalculation interval and 𝑁𝑟,𝑝 the total number of recalculation intervals depending on whether the mission is to Mars or Jupiter. Note that 𝑖can be 𝑁𝑟,𝑝 −1 at maximum, since at the last recalculation interval there are no future intervals. Additionally, there is another issue that needs to be tackled. If the dynamic pressure at a given instant is low enough, or the spacecraft is sufficiently far from the reference trajectory, the solver may not find a solution. The impossibility of solving the OCP (even when relaxing the final conditions further) suggests the implementation of an alternative Acta Astronautica 217 (2024) 116–129 123 J. Urrios et al. Table 1 Monte Carlo results for SHMPC strategy for different 𝜅𝑚𝑎𝑥,+and 𝛼𝑛. The value 𝑒denotes the mean arrival error distance to the planet and 𝜎𝑒 its standard deviation in units of fail distance [UF]. 𝜅𝑚𝑎𝑥,+𝛼𝑛such that Mars Jupiter S [%] F [%] 𝑒[UF♂]𝜎𝑒[UF♂] S [%] F [%] 𝑒[UFX]𝜎𝑒[UFX] 40% reduction 𝛼𝑚𝑎𝑥 49 21 1.53 4.04 88 7 0.28 1.04 𝑎𝑚𝑎𝑥 𝑟22 55 3.21 4.18 57 19 0.57 1.32 𝑎𝑚𝑎𝑥 𝜃43 33 2.04 3.97 76 14 0.39 1.13 20% reduction 𝛼𝑚𝑎𝑥 52 9 0.31 0.66 90 8 0.28 0.92 𝑎𝑚𝑎𝑥 𝑟32 14 0.45 1.07 73 13 0.40 1.22 𝑎𝑚𝑎𝑥 𝜃56 13 0.29 0.82 91 18 0.29 1.06 No reduction 𝛼𝑚𝑎𝑥 18 27 0.86 1.36 68 28 1.09 1.92 𝑎𝑚𝑎𝑥 𝑟14 34 1.41 2.20 57 32 1.22 2.06 𝑎𝑚𝑎𝑥 𝜃18 26 0.62 0.87 65 28 0.88 1.72 Linear increase 𝛼𝑚𝑎𝑥 46 14 0.37 0.73 85 8 0.30 1.03 𝑎𝑚𝑎𝑥 𝑟28 14 1.00 2.79 65 10 0.83 1.35 𝑎𝑚𝑎𝑥 𝜃41 13 0.48 1.02 80 11 0.55 0.98 sub-optimal control law strategy during that interval. Such a strategy must minimize the deviation of the spacecraft trajectory with respect to the reference, to maximize the likelihood that the solver finds a solution in the next recalculation. It can be argued that setting the tethers voltage to its maximum (that is, applying full thrust) is convenient, since the infeasibility of the problem comes ultimately from a low value of 𝑝⊕. However, the angle in which the E-sail should position itself is not obvious. Therefore, three sub-optimal control laws are considered in this work: (1) full thrust (𝑎𝑚𝑎𝑥) with 𝛼𝑛= arccos(1∕√3), which provides maximum thrust cone angle, (2) full thrust with 𝛼𝑛= 0◦, maximizing radial acceleration (𝑎𝑚𝑎𝑥 𝑟), and (3) full thrust with 𝛼𝑛= 45◦, which maximizes tangential acceleration (𝑎𝑚𝑎𝑥 𝜃). Because applying full thrust is convenient in these intervals, they are hereinafter denoted full thrust arcs. 4.2.1. Results for SHMPC In order to compare the guidance strategies and quantify their performance given the solar wind perturbations, the following success and failure criteria are defined. A mission is considered successful (S) if the spacecraft is able to reach the planet’s sphere of influence, while it is failed (F) depending on the miss distance, which is the distance at which the satellite remains from the planet, or 𝑒=‖𝒓(𝑡𝑓) − 𝒓𝑝‖. In the case of a Jupiter transfer, a fail is a miss of ten times (or more) the radius of Jupiter’s sphere of influence. For a Mars transfer, since its sphere of influence is considerably smaller, failed missions are the ones with a miss distance of fifty times (or more) the radius of Mars’ sphere of influence. These distances are considered high enough that no alternative propulsion system available can correct the error, hence the mission failure. A mission which ends with a miss distance in between the limits to be considered successful and failed is named partially successful, since it may be possible to reach the target, but employing considerable efforts from other propulsion systems. Since acceleration saturation phenomena depends on 𝑝⊕as shown in Section 4.1 and thus are stochastic, no deterministic metric of the performance of the strategies can be given. Instead, in order to statistically evaluate the different strategies presented in Section 4.2, a Monte Carlo analysis is conducted, simulating 100 missions for each of them, giving a total of 4⋅3⋅2⋅100 = 2400 simulations (four different 𝜅𝑚𝑎𝑥,+values and three full thrust arc strategies for missions to both Mars and Jupiter). A summary of the results obtained is presented in Table 1, where the mean miss distance 𝑒(the average of the error distances 𝑒of each simulation) and its standard deviation 𝜎𝑒are included to give an estimation on the precision and dispersion of the different strategies, and are expressed as multiples of the failing distance. Looking at the results, it becomes apparent that similar patterns emerge between missions to Mars and Jupiter regarding the most effective strategies. The key finding is that optimal results are achieved when a 20% reduction is Fig. 7. Histogram of the miss distance on missions to Mars with 20% 𝜅𝑚𝑎𝑥,+reduction and 𝛼𝑛||𝛼=𝛼𝑚𝑎𝑥 . considered for 𝜅𝑚𝑎𝑥,+, and during full thrust arcs, thrust is applied in the direction that maximizes the sail cone angle (𝛼𝑛= arccos(1∕√3)). The number of 100 simulations for each scenario was found to be a middle ground between computational effort and representativeness of the data. To check the representativeness, the outputs were recomputed for several random subsets of 50 of the 100 simulations, and the values obtained were sufficiently similar, justifying the chosen number of simulations. As observed in Table 1, many missions to Mars fail. A histogram of the miss distance 𝑒of the simulations with the obtained optimal configuration is included in Fig. 7, while the spacecraft final position in each of the 100 missions is shown in Fig. 9, and an example of a successful mission using this guidance strategy can be seen in Fig. 8. For the obtained optimal configuration, the final radial velocity mean and standard deviation were found to be 𝑣𝑟(𝑡𝑓) = 0.37 km/s, 𝜎𝑣𝑟(𝑡𝑓)= 1.05 km/s. Although with a high variance, this mean value of the radial velocity is small enough for another more specific propulsion system to handle the rendez-vous. 4.3. Receding Horizon MPC strategy In contrast to SHMPC approach presented in Section 4.2, in which arrival time is imposed, this section explores an alternative formulation that enables reaching the target planet in a greater time than initially planned, (𝑡𝑓> 𝑡𝑛𝑜𝑚 𝑓,𝑝 ) in exchange for robustness. To address this, a new optimization problem can be formulated with a cost function that includes not only the arrival position error, but also the position error at future points in time. Note that the reference orbit of the spacecraft is not only known for 𝑡<𝑡𝑛𝑜𝑚 𝑓, but also for future values, because, given that a rendez-vous is desired, the spacecraft’s reference orbit coincides with the planet’s for 𝑡 > 𝑡𝑛𝑜𝑚 𝑓. Therefore, these errors can be predicted by integrating the system into the future and comparing the Acta Astronautica 217 (2024) 116–129 124 J. Urrios et al. Fig. 8. Evolution of variables (left) and trajectory (right) in a successful mission to Mars in the presence of uncertainty using SHMPC with 20% 𝜅𝑚𝑎𝑥,+reduction and 𝛼𝑛||𝛼=𝛼𝑚𝑎𝑥 . Fig. 9. Final position error of the spacecraft in each of the Mars missions simulated for the 20% 𝜅𝑚𝑎𝑥,+reduction and 𝛼𝑛||𝛼=𝛼𝑚𝑎𝑥 SHMPC strategy. obtained and the reference orbits. That means that now there are no arrival constraints, as they are incorporated in the cost function. The idea therefore is to check for the error at the nominal arrival time; if it is larger than an established value, the program continues and the spacecraft keeps getting closer to the target, until a new check is made, then the process repeats. Once it gets sufficiently close (inside the sphere of influence), the mission ends. This approach, called Receding Horizon Model Predictive Control (RHMPC), is schematically shown in Fig. 10(b) and allows for the predicting horizon to recede rather than shrink. When the trajectory is recomputed for the first time, the only instant at which the error is measured is the nominal arrival time 𝑡𝑛𝑜𝑚 𝑓. When recomputing at an arbitrary time 𝑖=[1,…, 𝑁𝑟,𝑝], the current time is 𝑡𝑖, and the error is measured from 𝑡𝑛𝑜𝑚 𝑓until 𝑡𝑛𝑜𝑚 𝑓+𝑡𝑖. That means that the controller is always computing the orbit up to a time of 𝑡𝑛𝑜𝑚 𝑓into the future, although the cost function accounts for position errors only for 𝑡>𝑡𝑛𝑜𝑚 𝑓. The cost function formulated aims to minimize errors of the state variables. The following three position error functions are proposed: 𝜀1(𝑡) = 𝑟2(𝑡) + 𝑟2 𝑝− 2 𝑟(𝑡)𝑟𝑝cos (𝜃(𝑡) − 𝜃𝑝(𝑡)),(20) 𝜀2(𝑡) = (𝑟(𝑡) − 𝑟𝑝 𝑟𝑝)2 + (𝜃(𝑡) − 𝜃𝑝(𝑡))2,(21) 𝜀3(𝑡) = (𝑟(𝑡) − 𝑟𝑝 𝑟𝑝)2 +𝑐1(𝜃(𝑡) − 𝜃𝑝(𝑡))2.(22) It can be observed that the first error function is the exact error distance, but it is more computationally expensive. The second is similar enough for small error values, and the third is a modification of the second in which radial and angular errors are differently weighted. For the third, a weight coefficient 𝑐1is fine-tuned following a brief study to assess its impact on the success rate. A value of 𝑐1= 5 is demonstrated to yield the optimal results. Moreover, time-dependent weights are considered, to give priority to minimizing errors closer to the current time. The weight distribution 𝑤1(𝑡)assumes a constant distribution, and 𝑤2(𝑡)and 𝑤3(𝑡)decreasing linear and quadratic distributions respectively. Therefore, given the definition of the cost function 𝐽𝑘𝑙 =∫𝑡𝑛𝑜𝑚 𝑓+𝑡𝑖 𝑡𝑛𝑜𝑚 𝑓 𝑤𝑘(𝑡)𝜀𝑙(𝑡)d𝑡, (23) a total of nine different cost functions are considered since 𝑘= 1,2,3 and 𝑙= 1,2,3. With all these considerations the RHMPC strategy may be formulated for the successive 𝑁𝑟,𝑝 OCPs as min 𝜅,𝛼𝑛 𝐽𝑘𝑙 𝑡∈[𝑡𝑖, 𝑡𝑛𝑜𝑚 𝑓,𝑝 +𝑡𝑖] s.t. 𝑟 =𝑣𝑟  𝜃=𝑣𝜃 𝑟 𝑣𝑟=𝑣2 𝜃 𝑟−𝜇 𝑟2+𝜅 2(𝑟⊕ 𝑟)(1+cos2𝛼𝑛) 𝑣𝜃= − 𝑣𝑟𝑣𝜃 𝑟+𝜅 2(𝑟⊕ 𝑟)cos 𝛼𝑛sin 𝛼𝑛 ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟ state equations 0≤𝜅≤𝜅𝑚𝑎𝑥 −𝛼𝑚𝑎𝑥 𝑛≤𝛼𝑛≤𝛼𝑚𝑎𝑥 𝑛 ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟ path constraints 𝑟(𝑡𝑖) = 𝑟𝑖 𝜃(𝑡𝑖) =  𝜃𝑖 𝑣𝑟(𝑡𝑖) = 𝑣𝑟𝑖 𝑣𝜃(𝑡𝑖) = 𝑣𝜃𝑖 𝑡𝑓−𝑡𝑖=𝑡𝑛𝑜𝑚 𝑓,𝑝 , ⏟⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏟ boundary conditions (24) where the value of 𝑡𝑖refers to the time at which the orbit is recalculated, an instant at which the state variables 𝑟𝑖, 𝜃𝑖,𝑣𝑟𝑖and 𝑣𝜃𝑖are assumed to be known thanks to the navigation system. Recalling the results from Section 4.2.1, a 20% reduction in 𝜅𝑚𝑎𝑥 is applied when recomputing the orbit to leave some control margin and minimize saturation-related problems. 4.3.1. Results for RHMPC One hundred simulations are conducted for each of the nine cost functions, and for both Mars and Jupiter, resulting in a total of 1800