An ant-colony based optimization algorithm for interplanetary trajectories

Proyecto Fin de Carrera Ingeniería de Telecomunicación Formato de Publicación de la Escuela Técnica Superior de Ingeniería Autor: F. Javier Payán Somet Tutor: Juan José Murillo Fuentes Dep. Teoría de la Señal y Comunicaciones Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2013 Trabajo Fin de Grado Ingeniería Aeroespacial An Ant-Colony Based Optimization Algorithm for Interplanetary Trajectories Autor: Juan Antonio Villajos Molina Tutor: Rafael Vázquez Valenzuela, Julio César Sánchez Merino Dpto. de Ingeniería Aeroespacial y Mecánica de Fluidos Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2020 Trabajo Fin de Grado Ingeniería Aeroespacial An Ant-Colony Based Optimization Algorithm for Interplanetary Trajectories Autor: Juan Antonio Villajos Molina Tutor: Rafael Vázquez Valenzuela, Julio César Sánchez Merino Profesor Titular, Contrato Predoctoral Dpto. de Ingeniería Aeroespacial y Mecánica de Fluidos Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2020 Trabajo Fin de Grado: An Ant-Colony Based Optimization Algorithm for Interplanetary Trajectories Autor: Juan Antonio Villajos Molina Tutor: Rafael Vázquez Valenzuela, Julio César Sánchez Merino El tribunal nombrado para juzgar el trabajo arriba indicado, compuesto por los siguientes profesores: Presidente: Vocal/es: Secretario: acuerdan otorgarle la calificación de: El Secretario del Tribunal Fecha: Contents Nomenclature III 1 Introduction 1 1.1 Problem Statement 1 1.2 Scope of the Project 2 1.3 Structure 2 2 Basic Analysis Tools 3 2.1 Julian Days 3 2.2 Frame of reference 3 2.3 Orbital elements determination from position and velocity vectors 4 2.4 Position and velocity vectors determination from orbital elements 6 2.5 Planetary Ephemeris 7 2.6 Lambert’s Problem 8 3 Interplanetary Trajectories 13 3.1 Sphere of influence 13 3.2 Method of patched conics 14 3.2.1 Planetary departure, Geocentric phase 14 3.2.2 Heliocentric phase 15 3.2.3 Planetary arrival phase 16 3.3 Hohmann Transfer 16 3.4 Gravity Assist Maneuver (GAM) 17 3.4.1 Powered maneuver: additional periapse impulse 19 3.5 Deep Space Maneuver (DSM) 20 4 Optimization Procedure 21 4.1 General procedure 21 4.2 Optimization tool: fmincon 21 4.3 Direct Mission 23 4.3.1 Synodic Period 23 4.3.2 Pork Chop Plots 24 4.4 Multiple DSMs 25 4.5 GAM 26 4.5.1 Multiple GAMs (MGA) 26 4.6 DSMs in the GAM mission 27 5 Initial Guess Algorithms 29 5.1 Algorithm 1 29 5.1.1 DSMs 29 5.1.2 GAM 29 I II Contents 5.1.3 DSMs + GAM 30 5.2 Algorithm 2 30 5.3 Algorithm 3 31 6 Simulation Scenarios and Optimization Results 33 6.1 Mission to Mars 33 6.1.1 Direct Mission 33 6.1.2 Deep Space Maneuvers 34 6.2 Mission to Saturn 37 6.2.1 Direct Mission 37 6.2.2 Deep Space Maneuvers 37 6.2.3 Gravity Assist Maneuvers 40 Jupiter flyby 40 Mars flyby 41 Mars flyby + Jupiter flyby 41 Mars flyby + Earth flyby 43 Mars flyby + Earth flyby + Jupiter flyby 44 Venus flyby + Mars flyby + Jupiter flyby 45 Venus flyby + Earth flyby + Jupiter flyby 47 Final comparison 49 6.2.4 DSM + GAM 50 Jupiter flyby 50 Mars flyby 52 6.3 Mission to Jupiter 54 6.3.1 Direct Mission 54 6.3.2 Gravity Assist Maneuvers 55 Mars flyby 55 Venus flyby + Mars flyby 56 Venus flyby + Earth flyby 57 Mars flyby + Earth flyby 58 Venus flyby + Earth flyby + Mars flyby 59 Final comparison 61 7 Conclusions and Future Work 63 7.1 Conclusions 63 7.2 Future Work 64 List of Figures 65 List of Tables 67 Nomenclature The Sun ⊕The Earth ♂Mars YSaturn XJupiter ♀Venus a Semi major axis eEccentricity vector i Inclination θTrue anomaly ωArgument of periapsis ΩRight ascension of the ascending node hSpecific angular momentum µGravitational parameter εSpecific energy nNode line CF H Rotation matrix from Perifocal to Heliocentric frame of reference L Mean longitude M Mean anomaly E Eccentric anomaly RPPlanetary radius LPMean distance from the Sun to planet P RePRadius of the Sphere of Influence of planet P hpark Altitude of the parking orbit V∞, VPHyperbolic excess speed C3Mission characteristic energy rpPeriapse radius δTurn angle of the asymptote of the hyperbola Tsyn Synodic period GAM Gravity Assist Maneuver DSM Deep Space Maneuver MGA Multiple Gravity Assist Maneuvers III 6Chapter 2. Basic Analysis Tools •Case 1 (i = 0◦): θ=cos−1r·e r e ,if r·v≥0.θ=2π−cos−1r·e r e ,if r·v<0.(2.11) ¯ ω=cos−1ex e,if ey≥0.¯ ω=2π−cos−1ex e,if ey<0.(2.12) •Case 2 (e = 0): Ω=cos−1(nx),if ny≥0.Ω=2π−cos−1(nx),if ny<0.(2.13) u=cos−1r·n r n ,if rz≥0.u=2π−cos−1r·n r n ,if rz<0.(2.14) •Case 3 (i = 0◦, e = 0): λT=cos−1rx r,if ry≥0.λT=2π−cos−1rx r,if ry<0.(2.15) 2.4 Position and velocity vectors determination from orbital elements The aim of this section is to calculate the position and velocity vectors at a desired time, given the orbital elements at that time. First, both vectors are obtained on the Perifocal reference system and, after that, they must be transformed to the Heliocentric system. Consequently, in order to express those vectors on the Heliocentric frame, a rotation matrix is needed. Vectors expressed on the Perifocal reference system are denoted using the superscript "P". Similarly, the superscript "H" is used when working with vectors in the Heliocentric frame. rH=CF HTrF,(2.16) vH=CF HTvF.(2.17) Depending on the particular case under analysis, the rotation matrix as well as the position and velocity vectors in the Perifocal frame, may be modified. •General case: rF=p 1+ecos(θ)  cos(θ) sin(θ) 0 ,(2.18) vF=rµ p −sin(θ) e+cos(θ) 0 ,(2.19) CF H="cos(Ω)cos(ω)−sin(Ω)sin(ω)cos(i)sin(Ω)cos(ω)−cos(Ω)sin(ω)cos(i)sin(ω)sin(i) −cos(Ω)sin(ω)−sin(Ω)cos(ω)cos(i)−sin(Ω)sin(ω)+cos(Ω)cos(ω)cos(i)cos(ω)sin(i) sin(Ω)sin(ω)−cos(Ω)sin(i)cos(i)#.(2.20) • Case 1 (i = 0◦ ): The position and velocity vectors in the Perifocal system are the ones shown on equations 2.18 and 2.19. CF H=  cos(¯ ω)sin(¯ ω)0 −sin(¯ ω)cos(¯ ω)0 0 0 1 .(2.21) 2.5 Planetary Ephemeris 7 •Case 2 (e = 0): rF=a  cos(u) sin(u) 0 ,(2.22) vF=rµ a −sin(u) cos(u) 0 ,(2.23) CF H="cos(Ω)cos(u)−sin(Ω)sin(u)cos(i)sin(Ω)cos(u)−cos(Ω)sin(u)cos(i)sin(u)sin(i) −cos(Ω)sin(u)−sin(Ω)cos(u)cos(i)−sin(Ω)sin(u)+cos(Ω)cos(u)cos(i)cos(u)sin(i) sin(Ω)sin(u)−cos(Ω)sin(i)cos(i)#.(2.24) •Case 3 (i= 0◦, e = 0): rF=a  cos(λT) sin(λT) 0 ,(2.25) vF=rµ a −sin(λT) cos(λT) 0 ,(2.26) CF H=  cos(λT)sin(λT)0 −sin(λT)cos(λT)0 0 0 1 .(2.27) 2.5 Planetary Ephemeris For the purpose of designing realistic interplanetary missions, it is necessary to determine the state vector of a planet at any given time. Using the method developed in reference [ 7 ], it is possible to obtain the orbital elements of the planets employing their rates of change per century with respect to the J2000 epoch (1 January 2000, 12:00 UT). The foremost aspect of this procedure, which must be taken into account, is the fact that the valid time interval goes from 1800 AC to 2050 AC. Then, having the orbital elements, the state vector can be inferred as it has been shown. The method consists on a linear model, where the orbital elements are calculated as: (a,e,i,Ω,¯ ω,L) = (a0,e0,i0,Ω0,¯ ω0,L0)+ ∂ ∂t(a,e,i,Ω,¯ ω,L)·T,(2.28) where L=¯ ω+M , is the mean longitude and T0 is the number of centuries past J2000, obtained using the desired Julian day (JD) as: T=JD −2451545 36525 .(2.29) Once the mean anomaly M has been calculated using L and ¯ ω , the eccentric anomaly E is determined via Kepler’s Equation (2.30). M=E−esin(E).(2.30) Finally, the true anomaly θis obtained from equation 2.31. tan(θ/2) = r1+e 1−etan(E/2).(2.31) 8Chapter 2. Basic Analysis Tools Figure 2.3 Schematic of the Ephemeris tool. 2.6 Lambert’s Problem The appropriate analysis of an interplanetary mission leads to a comparison of myriad different transfer orbits. A way of tackling this question is solving Lambert’s Problem, which consists on determining a transfer orbit from two position vectors and a desired time of flight between them. Having said that, a tool which solves Lambert’s Problem becomes essential for the aim of this project. Figure 2.4 Example of a solution to Lambert’s Problem. Nowadays, solutions to Lambert’s Problem are used extensively in the orbital targeting problem, specifically in areas such as preliminary orbit determination and rendezvous analysis. Due to the growing complexity of the space missions, it is necessary to determine an efficient and robust Lambert solution algorithm. A large number of approaches have been formulated to solve this problem, each one involving a particular procedure and having a different associated efficiency. Some of them are the Universal Variable method, Battin’s formulation or Sun’s method. Once having investigated and analysed some of the possible solutions using reference [ 19 ], the method proposed in reference [ 2 ] is selected to solve Lambert’s Problem in this project, due to its simplicity despite the fact that an initial guess is needed. 2.6 Lambert’s Problem 9 Firstly, given the position vectors r1 and r2 of the points that must be linked using a transfer orbit ( P1 and P2as illustrated in figure 2.4), r1and r2determine the change in the true anomaly ∆θ, since: cos(∆θ) = r1·r2 r1r2 ,(2.32) where r1=√r1·r1,r2=√r2·r2.(2.33) Another parameter α , which dictates the transfer angle to be less or greater than 180◦ is determined from the third component of the specific angular momentum of the transfer orbit. That is: α=  0 0 1 ·(r1×r2).(2.34) Moreover, it would be possible to select a direct orbit (i ≤90◦ ) or a retrograde one (i > 90◦ ). Nonetheless, only direct orbits are considered. To start the procedure’s explanation, it must be clarified that using the Lagrange coefficients f and g, if the position and velocity of an orbiting body are known at a given instant, it is possible to obtain them at any later time in terms of the initial values. Consequently: r2=fr1+gv1,(2.35) v2=˙ fr1+˙ gv1.(2.36) Solving these equations and considering that f˙g−˙ f g =1 , both velocities are obtained in terms of both positions and the Lagrange coefficients: v1=1 g(r2−fr1)(2.37) v2=1 g(˙ gr2−r1)(2.38) Therefore, Lambert’s Problem is solved once the Lagrange coefficients f, g and ˙ g are determined. The Lagrange coefficients f and g, as well as their time derivatives are listed as functions of the change in the true anomaly: f=1−µr2 h2(1−cos(∆θ)),g=r1r2 hsin(∆θ),(2.39) ˙ f=µ h 1−cos(∆θ) sin(∆θ)µ h2(1−cos(∆θ))−1 r1−1 r2,˙ g=1−µr1 h2(1−cos(∆θ)).(2.40) It is also possible to express these quantities in terms of the universal anomaly χ: f=1−χ2 r1 C(z),g=∆t−a √µχ3S(z),(2.41) ˙ f=√µ r1r2 χ[zS(z)−1],˙ g=1−χ2 r2 C(z),(2.42) where z=βχ2 ( β is the reciprocal of the semimajor axis of the unknown orbit). A relationship between ∆t and ∆θcan be found by equating the two expressions for g: r1r2 hsin(∆θ) = ∆t−1 √µχ3S(z).(2.43) 10 Chapter 2. Basic Analysis Tools Then, to eliminate the angular momentum h, both expressions for f are equated, obtaining: h=sµr1r2(1−cos(∆θ)) χ2C(z)(2.44) Accordingly, equation 2.43 is simplified : √µ∆t=χ3S(z)+AχpC(z),(2.45) where the symbol A refers to: A=sin(∆θ)rr1r2 1−cos(∆θ).(2.46) In order to discover a relationship between z and χ, the function y(z) is defined as: y(z) = r1+r2+AzS(z)−1 pC(z)(2.47) Hence, χ=py(z)/C(z). Substituting it back into equation 2.48 yields: √µ∆t=y(z) C(z)3 2 S(z)+Apy(z),(2.48) After this, it is possible to iteratively apply Newton’s method to solve for z, using the function F(z): F(z) = y(z) C(z)3 2 S(z)+Apy(z)−√µ∆t.(2.49) Its derivative is defined as: ∂F(z) ∂z=       hy(z) C(z)i3 2h1 2zC(z)−3S(z) 2C(z)+3S(z)2 4C(z)i+A 8h3S(z) C(z)py(z)+AqC(z) y(z)i,if z 6=0. √2 40 y(0)3 2+A 8hpy(0)+Aq1 2y(0)i,if z =0. (2.50) Lastly, the functions C and S are presented: C(z) =            1−cos(√z) z,if z >0. cosh(√−z)−1 −z,if z <0. 1 2,if z =0. (2.51) S(z) =            √z−sin(√z) (√z)3,if z >0. sinh(√−z)−(√−z) (√−z)3,if z <0. 1 6,if z =0. (2.52) At this point, an initial guess is required to start iterating with the optimisation parameter of the algorithm ( z ). In order to improve the efficiency of the algorithm, a process to find a favourable initial guess is included. It consists on, starting from a known solution (it is known that z is positive for elliptic orbits and negative for hyperbolic orbits), increasing its value with a non accurate step, that is to say not as small as the optimisation one, until the objective function changes its sign. Hence the optimal solution would be close to that initial guess. Having calculated an initial guess, the iteration procedure begins. The iteration ends when a solution is found, having established a termination tolerance. In addition to that tolerance, a maximum number of 2.6 Lambert’s Problem 11 iterations is set. Therefore, after having used Lambert’s Algorithm, it is necessary to check whether the solution is correct or the algorithm has stopped by reaching the maximum number of iterations fixed, thus the solution could not be used. Summarising, this algorithm calculates the velocities associated with both points, P1 and P2 , in the transfer orbit. Eventually, using the transformations explained in Section 2.3, the orbital elements are determined. Figure 2.5 Schematic of the Lambert tool. 3 Interplanetary Trajectories The overall target of this project is to optimize interplanetary trajectories, hence the procedure which has been followed to calculate the required orbits must be explained. The method of patched conics allows to solve a problem which involves more than two bodies. In order to apply this method, developed in Section 3.2, it is necessary to introduce a meaningful concept, the Sphere of Influence, in Section 3.1. Moreover, since the work is carried out on the 2-D problem, planetary orbits around the Sun are considered to be on the Ecliptic Plane. Therefore, only two vector components for each variable have to be accounted for. To conduct the analysis of planetary missions, general data of every planet is required. This general data is summarised on table 3.1, where RP is the planet’s radius, LP is the mean distance from the Sun to the planet and µPis the planet’s gravitational parameter. Planet RP(km) LP(AU) µPkm3/s2 Mercury 2439.7 0.3871 22032.1 Venus 6051.8 0.7233 324858.8 Earth 6378.14 1 398600.4 Mars 3397 1.5237 42828.3 Jupiter 71492 5.2033 126711995.4 Saturn 60268 9.5808 37939519.7 Uranus 25559 19.2709 5780158.5 Neptune 24764 30.1927 6871307.8 Table 3.1 Planetary general data. 3.1 Sphere of influence In accordance with reference [ 18 ], the Sphere of Influence is the region in which the gravitational influence of a body prevails over the influence of other bodies. Similarly, the Sphere of Influence radius is a reasonable estimate of the distance beyond which the Sun’s gravitational attraction dominates that of a planet. In order to appraise that radius, the Three-Body Problem is considered. It comprises a planet "p" of mass mp , the Sun "s" of mass ms and a space vehicle "v" of mass mv (represented on figure 3.1). To give an illustration, when the vehicle is close to the Sun, the effect of the planet can be treated as a perturbation on the Keplerian motion of the vehicle around the Sun: ¨ Rv=−µs Rv R3 v +µp r r3−R R3 =gs+g0 p,(3.1) 13 14 Chapter 3. Interplanetary Trajectories where r, R and Rv are the norm of the vectors r , R and Rv , respectively. Likewise, if the vehicle is close to the planet, the effect of the Sun can be modeled as a perturbation: ¨ r=−µp r r3+µs Rv R3 v−R R3 =gp+g0 s.(3.2) Figure 3.1 Three-Body Problem representation, pulled from [2]. As a result, the geometric place of the points where g0p gs=g0 s gp separates the influence zones of the Sun and the planet. Nevertheless, the Sun is the dominant celestial body in the Solar System, over 1000 times more massive than Jupiter and over 300000 times more massive than Earth. Therefore, as ms>> mp , that geometric place is similar to a sphere around the planet, which is called Sphere of Influence. Inside the Sphere of Influence, it is possible to consider that the only body which affects the vehicle is the planet. Consequently, the vehicle travels an unperturbed Keplerian path around the planet, simplifying the problem to a Two-Body Problem. On the other hand, outside the Sphere, the only influence is the Sun, thus it is assumed that the spacecraft follows an unperturbed Keplerian orbit around the Sun. The radius of the Sphere is given by: Rep=Lpµp µs2 5 .(3.3) For every planet, it is true that Rp<< Rep<< Lp , hence it is possible to simplify the interplanetary problem. When working on the planet-centric phase, Rep→∞ , and from the heliocentric point of view, Rep→0, so the Sphere of Influence is considered as a point (planet position). 3.2 Method of patched conics The method of patched conics is employed to divide the mission into three parts: the hyperbolic departure trajectory relative to the home planet, the cruise ellipse relative to the Sun and the hyperbolic arrival trajectory relative to the target planet. In order to illustrate the method, a mission from the Earth to Planet 2 is analysed. 3.2.1 Planetary departure, Geocentric phase The vehicle departs from a park orbit at an altitude hpark1 from the Earth’s surface ( rpark1 =R⊕+hpark1 ). At first, it must escape the gravitational pull of the Earth as it is inside its Sphere of Influence, travelling a hyperbolic trajectory whose focus is the Earth, in a Geocentric reference system. ∆Vini is the impulse required at the park orbit to insert the vehicle in a scape trajectory, V⊕ is the Earth’s velocity and V∞1 is the 3.2 Method of patched conics 15 hyperbolic excess speed (when the spacecraft arrives at the limit of the Earth’s Sphere of Influence). ∆Vini =s2µ⊕ rpark1 +(V∞1)2−sµ⊕ rpark1 .(3.4) Figure 3.2 Departure of a spacecraft on a mission from an inner planet to an outer planet, extracted from [2]. Figure 3.3 Departure of a spacecraft on a mission from an outer planet to an inner planet, extracted from [2]. 3.2.2 Heliocentric phase The Heliocentric phase consists on a trajectory, with the Sun as the focus, which starts when the spacecraft leaves the Earth’s Sphere of Influence and ends when it arrives at the Sphere of Influence of Planet 2. After escaping from the Earth’s gravitational pull, in the Heliocentric reference system, the vehicle travels with a velocity VH 0=V⊕+V∞1 if it goes to a superior planet and VH 0=V⊕−V∞1 if it travels to an inferior planet. Furthermore, it arrives at the Sphere of Influence of Planet 2 having a velocity VH 1 . The semimajor axis of the heliocentric conic is a, and the symbol refers to the Sun. V∞1=± s2µ L⊕−µ a−sµ L⊕!=±VH 0−V⊕.(3.5) Figure 3.4 Heliocentric phase representation. 22 Chapter 4. Optimization Procedure Figure 4.1 Diagram of the general procedure. tool is to find the minimum value of a function, called Objective Function, given the constraints that must be satisfied. Tu put it another way, it is defined as a non-linear programming solver. The motive of using a local optimal searcher is to try to obtain equally feasible and optimal solutions as using a Metaheuristic-based algorithm such as a genetic algorithm but reducing the associated computation time. The syntax of the function fmincon is simple, using xas the vector of the variables of the problem: minxf(x)such that                        c(x)≤0, ceq(x) = 0, A·x≤b, Aeq ·x=beq, lb ≤x≤ub, (4.1) where c(x) is the non-linear constraints vector, ceq(x) is the equality constraints vector, A and b are a real matrix and a real vector, respectively, which refer to the linear inequality constraints, Aeq and beq are a real matrix and a real vector, respectively, which refer to the linear equality constraints, lb and ub are vectors where the lower and upper bounds, respectively, are established for each variable. The required inputs include the Objective Function f, and an initial guess x 0 indeed. Furthermore, it is possible to set a wide range of optimization options like: • Algorithm: The optimization algorithm used. The Interior Point Algorithm is selected since it can handle large and sparse problems, as well as small dense ones. It is a large-scale algorithm. •Tolerances: For instance, it may be defined a tolerance for a constraint violation or for the step. • Maximum number of evaluations or iterations: It is determined for the function evaluation and to define the maximum number of iterations that are allowed. • Display: It provides an output in each iteration which shows information about the optimization procedure. Eventually, the solution found by fmincon might not be a local optimal point, since the solver could have stopped searching due to a maximum number of iterations achieved. Therefore, the quality of the "optimal" solution obtained must be checked. The output exitflag makes it possible, returning a different number depending on the nature of the solution. For instance, if it has finished by achieving the maximum number of iterations it returns a 0 and if it has not found a feasible point it returns a -2. 4.3 Direct Mission 23 4.3 Direct Mission First, the direct mission involving only two impulses is analysed for each scenario. The mission begins from a parking orbit at an altitude hpark1 over the departure planet’s surface. The second impulse is performed at an altitude hpark2 =9RP over the destiny planet, where RP is the radius of that planet. This decision was made in order to define the objective of each mission, which is setting the spacecraft at a circular orbit around the arrival planet. Even though That orbit’s radius would be large enough to have a minor influence, it entails an arrival velocity limitation. The algorithm implemented to solve the direct mission uses only the departure date (Tini) and the arrival date (Tend) as variables: 1. Using the Ephemeris of both planets, their position and velocity at the desired time ( Tini for the departure planet and Tend for the destiny planet) are obtained. 2. Using Lambert tool, given the position of each planet, the departure and the arrival heliocentric velocities are calculated. 3. Using the heliocentric velocities and the planet’s ones, the excess velocities and the parameter C3 related to the departure and the arrival impulses are determined. Then, both impulses are calculated: ∆Vini =s2µPini rpark1 +C3ini −sµPini rpark1  ,(4.2) ∆Vend =s2µPend rpark2 +C3end −sµPend rpark2  ,(4.3) where µPini and µPend are the gravitational parameters of the departure planet and the arrival planet, respectively. C3ini =V2 ∞ini refers to the departure planet and C3end =V2 ∞end , refers to the arrival planet. 4. Determination of the total impulse: ∆Vtot =∆Vini +∆Vend.(4.4) Figure 4.2 Schematic of the method developed for the direct mission. 4.3.1 Synodic Period The synodic period ( Tsyn ) is the amount of time necessary for the repetition of a relative configuration between two planets. In other words, how much time it takes to wait until a relative position of two planets recurs. It is calculated using the heliocentric period of the involved planets ( TP1 and TP2 ), which is defined as the time required to complete a revolution around the Sun. Tsyn =1 |1/TP1−1/TP2|.(4.5) The synodic period becomes useful in order to consider the departure time window, hence the range of the departure time is set including a range of time similar to this period. 24 Chapter 4. Optimization Procedure 4.3.2 Pork Chop Plots A common strategy used for the design of interplanetary missions is to represent the pork-chop plot. It consists on a diagram where contour lines of the mission parameter C3 are represented as a function of the departure and arrival dates. The parameter C3 is a measure of the specific energy, calculated from the excess velocity, V∞ , as C3=V2 ∞ . Thus, it is an appropriate way of comparing trajectories through the required energy as well as through its feasibility. Generally, for the 3-D problem, if the orbits are elliptical and quasi-planar, there are two optimal trajectories, as it is shown in figure 4.3. The first one is known as a trajectory of Type I, faster than Hohmann trajectory ( ∆θ<180 ), and the other is known as a trajectory of Type II, slower than Hohmann’s ( ∆θ>180 ). Even so, it should be noted that Hohmann trajectories only exist in the simplified case of circular planetary orbits. Nevertheless, since the 2-D problem involving elliptical and co-planar orbits is considered, both optimal possibilities may not exist. There might be more than one local minimum solutions, yet it is possible to find the best one, as it is explained in references [ 9 ] and [ 10 ]. In fact, the results calculated, which are shown in chapter 6, prove that one clear optimal solution could be identified. Figure 4.3 Pork chop plot for an Earth to Mars mission (3D problem), pulled from [3]. 4.4 Multiple DSMs 25 4.4 Multiple DSMs Secondly, the mission with multiple Deep Space Maneuvers in addition to the departure and the arrival impulses is analysed. The required variables for this function are the departure date ( Tini ) and the arrival date ( Tend ), as well as the date ( TDSMi ) and the position ( rxi and ryi ) where each Deep Space Maneuver is performed. As a result, the number of variables is increased as the number of DSM is extended. The algorithm used to tackle this problem is simple: 1. Using the dates, the time of flight between each impulse is obtained. In addition, using the Ephemeris tool, the position of the planets is calculated. 2. Lambert tool is employed to obtain the heliocentric velocities at the beginning and the end of each segment, as it is shown in figure 4.4. 3. Using the heliocentric departure and arrival velocities, the excess speeds and the parameters C3 are evaluated. Subsequently, the initial and final impulses are calculated using equations 4.2 and 4.3. 4. Calculation of the DSM impulses: looking at figure 4.4, having the heliocentric velocity before and after the i-DSM, it is possible to obtain the impulse required in that DSM ( ∆VDSMi ) using equation 3.18. 5. Determination of the total impulse: ∆Vtot =∆Vini +∆Vend + i=n ∑ i=1 ∆VDSMi.(4.6) Figure 4.4 Schematic of the multiple DSM problem for the case of n=2. Figure 4.5 Diagram of the multiple DSM mission method. As regards the optimisation parameters, the only special constraint, in addition to the fact that the DSM’s dates must take place between the departure and the arrival dates, is associated with their position. The DSM position range is considered between the inferior planet orbit and 1.2 times the superior planet orbit radius. 26 Chapter 4. Optimization Procedure 4.5 GAM The next stage analysed is a mission which involves a Gravity Assist Maneuver (GAM). A powered maneuver is considered by adding an extra impulse at the periapsis as seen by the flyby planet. The method requires as variables the departure date ( Tini ), the arrival date ( Tend ), the flyby date ( TGAM ) and the periapse radius of the hyperbola (rp) where the extra-impulse is performed. The algorithm developed is based on the ideas that are shown on reference [ 6 ]. The match point is the Gravity Assist Maneuver, thus the trajectory is propagated forward in time from the start point, which is the departure planet, until that match point. Then, the trajectory is propagated backwards in time from the end point, which is the destiny planet, to the match point. There, at the match point, the constraints are applied. The procedure consists on: 1. Forward propagation: using planetary Ephemeris and Lambert tool, the departure heliocentric velocity Vini and the heliocentric velocity when the spacecraft arrives at the flyby planet V+ GAM are determined. 2. Backwards propagation: using planetary Ephemeris and Lambert tool, the arrival heliocentric velocity at the destiny planet Vend and the heliocentric velocity when the spacecraft leaves the flyby planet V− GAM are calculated. 3. Using Vini , the excess velocity and its associated parameter C3 are obtained. Accordingly, the initial impulse is calculated using equation 4.2. 4. Using Vend , the excess velocity and its associated parameter C3 are determined. Then, the final impulse is calculated using equation 4.3. 5. Applying the method described in Section 3.4, the required extra-impulse at periapsis ∆Vp is obtained via V+ GAM and V− GAM. 6. Determination of the total impulse: ∆Vtot =∆Vini +∆Vend +∆Vp.(4.7) Figure 4.6 Diagram of the method implemented to solve the Gravity Assist Maneuver problem. With respect to the optimization parameters, there is a minimum periapse radius which depends on the planet under consideration. 4.5.1 Multiple GAMs (MGA) Including more than one Gravity Assist Maneuver is also studied. Consequently, the number of variables depends on the number of maneuvers, adding the periapse radius and the date of each additional GAM to the variables already defined. In that case, all the propagation is done forward in time. Hence using Lambert tool and the Ephemeris of each pair of planets, the velocities are calculated, and then, the extra-impulses at periapsis are determined for each flyby. 4.6 DSMs in the GAM mission 27 4.6 DSMs in the GAM mission Finally, a problem which involves both DSMs and GAM is considered is solved. It consists on a mission to a destiny planet via another planet flyby, including the possibility of adding n DSMs in the first transfer and m DSMs in the second one. The variables of this problem are the departure date (Tini), the arrival date (Tend), the flyby date ( TGAM ), the periapse radius of the hyperbola associated with the flyby, as well as the date and the position of every DSM included. 1. Multiple DSM from the departure planet to the flyby planet. The following parameters are obtained: •Initial impulse : ∆Vini. •Every DSM impulse of the first phase: ∆VDSM11,...,∆VDSM1n. •Heliocentric velocity when the spacecraft arrives at the flyby planet V− GAM. 2. Multiple DSM from the flyby planet to the destiny planet. The following parameters are determined: •Heliocentric velocity when the spacecraft leaves the flyby planet V+ GAM. •Every DSM impulse of the second phase: ∆VDSM21,...,∆VDSM2m. •Final impulse : ∆Vend. 3. Gravity Assist Maneuver: calculation of the required extra-impulse at periapsis ∆Vp , via V+ GAM and V− GAM. 4. Determination of the total impulse: ∆Vtot =∆Vini +∆Vend + i=n ∑ i=1 ∆VDSM1i+ j=m ∑ j=1 ∆VDSM2j+∆Vp.(4.8) Figure 4.7 Diagram of the procedure developed to solve a problem which involves a Gravity Assist Maneuver including multiple Deep Space Maneuvers in both phases. 5 Initial Guess Algorithms Assuming that the optimization process is carried out using the Matlab tool fmincon, the selection of the initial guess is paramount. As a matter of fact, fmincon tries to find a local minimum of the considered function, that search extremely depends on the initial point used. In this project, different ways of selecting the initial guess are evaluated in order to compare them. Algorithm 1 requires a previous investigation to find a suitable initial guess, yet it is associated with a short computation time. Algorithm 2 does not need a preceding analysis but explores a wide range of possible initial guesses, hence its computational burden increases. Lastly, Algorithm 3 includes an Ant Colony Optimization procedure allowing to focus on the best solutions previously found. Therefore, the last algorithm can improve the solutions obtained by Algorithm 2. 5.1 Algorithm 1 This first method consists on inferring the initial solution of the considered problem using a previous simpler problem. Hence the optimal solution of a basic problem associated with the one under analysis must be evaluated. Then, that solution is used as the initial guess to solve the complex problem. To illustrate this method, the procedure followed to determine the initial guess for each variable of every problem is explained. 5.1.1 DSMs In this case, the variables under consideration are Tini , Tend , as well as TDSMi , rxi and ryi , with i=1,...,n , for each DSM. The initial guess for each variable is determined using this method: •T[0] ini and T[0] end: Using the obtained optimal departure and arrival dates for the direct problem. •T[0] DSMi : The time interval T[0] ini −T[0] end is divided into n+1 sections of the same length, tgap . The initial guess for the date of the i-DSM is T[0] DSMi=T[0] ini +i·tgap. •r[0] xi and r[0] yi : Dividing the distance between both planets involved equidistantly, into n+1 sections of length dgap . The initial guess for the components of the position rx and ry of the i-DSM is r[0] x|yi= LPini ±i·dgap , where LPini is the mean distance from the Sun to the initial planet. The sign depends on whether the vehicle is going to a superior or to an inferior planet. After having set these initial values, if the Objective Function is not defined at the evaluated point or its value is very high, the sign and the value of each variable r[0] xiand r[0] yiare randomly modified. 5.1.2 GAM •T[0] ini: Using the optimal departure date for the direct transfer to the first flyby planet. •T[0] GAM : If it is the first GAM, it is calculated as the optimal arrival date for the direct transfer from the departure planet. In case of multiple flybys, the date of each GAM is obtained dividing the problem into only 1 GAM problems and using their optimal dates, employing previously solved problems or adding the respective Hohmann’s times of flight between each planet. For instance, when studying a problem Earth-Mars-Earth-Saturn, the simpler problems Earth-Mars-Earth and Mars-Earth-Saturn are analysed 29 30 Chapter 5. Initial Guess Algorithms to obtain an initial guess for the complete problem. However, if the mission Earth-Mars-Saturn have already been studied, it might be used the optimal solution for that problem to set the initial guess for the date of the Mars flyby and then the time of flight of the Hohmann transfer Mars-Earth may be added to calculate an initial guess for the Earth flyby. •T[0] end : Adding the time of flight of the Hohmann transfer between the last flyby planet and the destiny planet to the date of the last flyby. •r[0] p: It is set as 1.1 times the minimum periapse radius of each flyby planet. 5.1.3 DSMs + GAM The method consists on a combination of both the DSMs problem and the GAM problem. Each initial guess is obtained using the previously mentioned procedure. 5.2 Algorithm 2 The foremost problem of Algorithm 1 is the fact that the optimization solver finds only the closest local minimum to the initial guess. Nevertheless, as the problem is highly non-linear, thus probably non-convex, there may be many different local minimum points in the search space. Accordingly, if a local minimum is found, other possible solutions are left out of consideration. One way to overcome this problem in order to augment the chances of computing the global minimum, is to calculate multiple local minimums providing the local optimization solver with different initial guesses. After this, they are compared, and the best solution is selected. Moreover, the initial guesses must be diverse so that the found local optimums would be different. The main strength of this algorithm is the fact that it explores myriad varied possibilities. Once having determined the search space for these two variables, Tini and Tend , each range is divided equidistantly (a same number of divisions for each variable is not compulsory). The next stage is combining them in order to create different pairs of initial guesses, and then obtaining the rest of the variables which are required for each pair. Obviously, this method’s computation time is higher than the one associated with Algorithm 1, yet it makes possible to discover new optimal points. Figure 5.1 illustrates an example of the procedure implemented by Algorithm 2, taking only two initial dates and three final dates (denoted as 2-3). As it is shown, if any of the obtained initial guesses are not feasible, in this case due to the fact that the arrival date is previous to the departure one, they are eliminated. The method is focused on calculating several initial guesses via the departure time Tini and the arrival time Tend . The rest of the variables of the problem for the case of multiple DSM are obtained following the same procedure developed in Algorithm 1. Exceptionally, when using this algorithm on Gravity Assist Maneuvers problems, the initial guesses for the dates of the flybys are calculated following this method, depending on the number of GAMs: • 1 GAM case: for each pair of Tini and Tend , the time interval between them is divided into 4 gaps, obtaining 3 intermediate dates. These 3 dates are used as initial guesses for the date of the GAM, hence each pair of Tini and Tend produces 3 different initial guesses. • 2 GAM case: for each pair of Tini and Tend , the time interval between them is divided into 5 gaps, obtaining 4 intermediate dates. The first 2 dates are employed as initial guesses for the date of the first GAM. Then, the last 2 dates are used as initial guesses for the date of the second GAM. After this, they are combined, thus for each pair of Tini and Tend, 4 different initial guesses are determined. • 3 GAM case: for each pair of Tini and Tend , the time interval between them is divided into 7 gaps obtaining 6 intermediate dates. The first 2 dates are used as initial guesses for the date of the first GAM, dates 3 and 4 are employed as initial guesses for the date of the second GAM, and the last 2 dates are used as initial guesses for the date of the third GAM. A couple of different initial guesses is obtained for each pair of Tini and Tend . The first involves the intermediate dates 1, 3 and 5 and the second, dates 2, 4 and 6. 5.3 Algorithm 3 31 Figure 5.1 Illustration of the method developed in Algorithm 2. 5.3 Algorithm 3 A challenging problem which arises when applying Algorithm 2 is the lack of feedback. To put it another way, after having obtained the results via Algorithm 2, it would make sense to exploit the best solutions in order to enhance them. An alternative approach to tackle this problem is to include a Metaheuristic-based Ant Colony Optimization (ACO) adapted to continuous optimization. As explained in reference [ 4 ], the inspiring source of ACO is the foraging behaviour of real ants when searching for food. First, they aimlessly explore the area surrounding their nest. When an ant discovers a food source, it evaluates it and carries a portion back to the nest leaving a pheromone trail on the ground during the return trip. The amount of pheromone deposited depends on the quality and the quantity of the food source, which has been previously evaluated. Therefore, it has been proved that indirect communication among ants enables them to discover the shortest path between their nest and food sources. This capability of real ant colonies has inspired the definition of algorithms that can solve hard combinatorial optimization problems, such as the one presented in this project. The overall idea of ACO is the incremental construction of solutions based on the biased probabilistic choice of solutions components, determined by the pheromone. At each construction step, ants make a probabilistic choice of the solution component. Moreover, ants sample a probability density function (PDF) for each variable which defines the probability of selecting a certain value for that variable. A Gaussian kernel PDF is used due to the fact that it can identify more than one promising solution, as it may have more than one maximum. A Gaussian kernel function is represented in figure 5.2. For the variable i of the problem whose vector of variables is x, it is denoted as Gi(x). Gi(x) = k ∑ l=1 ωlgi l(x) = k ∑ l=1 ωl 1 σi l√2πe−(x−µi l)2 2(σi l)2 .(5.1) There must be defined as many Gaussian kernel PDFs as the number of problem variables (i = 1 ,...,n) and the design parameter k represents the number of ants used to look for a solution (necessarily k>n ). Each of the PDF is parameterized with three vectors: •µi is the vector of means. In each iteration, the solutions which were obtained and ordered in the previous iteration become the means of the actual one. µi=µi 1,..., µi l=si−1 1,...,si−1 l.(5.2) 38 Chapter 6. Simulation Scenarios and Optimization Results Figure 6.4 Pork chop plot: Direct mission to Saturn. ∆Vtotal (km/s) TOF (days) Computation time (s) Algorithm 1 11.948413 2172 33.886305 Algorithm 2 (3-10) 11.965802 6105 261.231823 Algorithm 3 (k=50, q=0.1) 11.948408 2172 638.306790 Table 6.7 Optimal solution for the mission to Saturn via 2 DSMs. In view of the results, it appears that Algorithms 2 and 3 can discover a better solution in terms of the cost than Algorithm 1. Particularly, the optimal solutions using Algorithm 3 are presented on table 6.8. It has been found that the optimal solution is not to include DSMs, since the impulses are below tolerance (reaching mm/s accuracy would be unfeasible for a DSM thus these impulses are considered as 0), and it involves a complete revolution around the Sun (figure 6.6). That solution could not be obtained in Section 6.2.1 since the Lambert’s Problem tool does not have the possibility of calculating multiple revolution orbits. Departure Date Arrival Date ∆Vini (km/s) ∆Vend (km/s) ∆VDSM1(km/s) n=1 25/11/2025 10/01/2045 7.272395 4.434734 4.520 ·10−8 n=2 01/06/2025 13/05/2031 7.292580 4.655827 8.193 ·10−7 DSM1Date ∆VDSM2(km/s) DSM2Date 22/01/2035 / / 07/05/2027 3.702 ·10−710/05/2029 Table 6.8 Optimal solution for the mission to Saturn via 1 and 2 DSMs via Algorithm 3. In brief, the cost of the direct transfer has not been reduced enough. Particularly, the initial impulse is actually large. In conclusion, these impulses are expensive to be performed, hence another alternative method 6.2 Mission to Saturn 39 Figure 6.5 Optimal solution for the direct mission to Saturn. Figure 6.6 Optimal solution for the mission to Saturn using Algorithm 3. must be found to make this mission possible. 40 Chapter 6. Simulation Scenarios and Optimization Results 6.2.3 Gravity Assist Maneuvers The next stage, looking forward to reducing the mission cost, is using a Gravity Assist Maneuver. This may increase the time of flight needed, but hopefully it could diminish the total impulse required. First, diverse cases are separately analysed, and then they are all compared. Generally, the range for the initial time lies between 01/01/2024 and 01/01/2026, since it is the required wider range comparing Jupiter’s and Mars’ flybys. Jupiter flyby The first case under analysis is a mission to Saturn via swing by Jupiter. To apply Algorithm 1, the direct transfer Earth-Jupiter must be evaluated, to obtain an initial flight time estimation. The Earth-Jupiter synodic period is Tsyn =399 days. Consulting reference [ 3 ], it is decided to search an optimal departure date around September 2025. Using fmincon, an optimal solution for the direct mission to Jupiter is detected for the C3 parameter (figure 6.18) as the interest lies in the flyby over Jupiter. The optimal point is almost a Hohmann transfer with a flight time of 1056 days, leaving the Earth on 25 September 2025 and arriving at Jupiter on 16 august 2028. Using these dates, as well as adding 10 years for the trajectory between Jupiter and Saturn, which is the time of flight of a Hohmann transfer between them, the initial guess is determined: Tini =01/09/2025, TGAM =01/08/2028, Tend =01/08/2038. Figure 6.7 Optimal solution for the mission to Saturn via GAM around Jupiter using algorithm 3. To apply Algorithms 2 and 3, the range of possible arrival dates must be set. Having calculated both times of flight associated with the Hohmann transfer for Earth-Jupiter (6 years) and for Jupiter-Saturn (10 years), the range of the arrival date is set as (01/01/2035, 01/01/2045). Then, the optimal solution using 3 divisions for Tini and 5 for Tend is calculated using Algorithm 2. The result is almost the same using each algorithm as it is represented on table 6.9. Although the computation time is extremely larger using Algorithm 3, a similar solution is obtained with every tool. However, Algorithm 1 requires a previous analysis of the direct mission individually, which is not necessary to apply Algorithm 2 and Algorithm 3. In brief, including a Gravity Assist Maneuver reduces the total impulse required although it increases the time of flight of the mission. The optimal fuel consumption solution found using Algorithm 3, which reduces the cost of the direct mission by 13 % and increased the time of flight by 6.2 Mission to Saturn 41 ∆Vtotal (km/s) TOF (days) Computation time (s) Algorithm 1 10.418022 4956 2.554828 Algorithm 2 (3-5) 10.417692 4976 59.048049 Algorithm 3 (k=150, q=0.3) 10.417683 4980 404.416760 Table 6.9 Comparison of the optimal solution for the mission to Saturn via GAM around Jupiter obtained using each algorithm. 2900 days, is presented on table 6.10 and figure 6.7, where the symbol Xrefers to the planet Jupiter. Departure Date Arrival Date GAM Date ∆Vini (km/s) ∆Vend (km/s) ∆VGAM (km/s) 06/11/2025 27/06/2039 07/11/2030 6.603505 3.814178 3.730 ·10−8 rp/RXδ(º) 10.5 97.717 Table 6.10 Optimal solution for the mission to Saturn via GAM around Jupiter using Algorithm 3. Mars flyby Secondly, a Gravity Assist Maneuver around Mars is considered. The optimal solution for the Earth-Mars direct transfer is already known but in a different time interval. Consequently, it is necessary to search an optimal point within the range of dates established. It involves a departure date around August of 2025, an arrival date at Mars around July of 2027 and a time of flight of about 700 days. Calculating the time of flight of the Hohmann transfer for Mars-Saturn, 2390 days (6.5 years) and adding it to the date of the Mars flyby, the initial guess for the arrival date at Saturn is obtained (January 2033). Having calculated both times of flight associated with the Hohmann transfer for Earth-Mars (260 days) and for Mars-Saturn (6.5 years), the range of the arrival date is established as (01/01/2028, 01/01/2040). ∆Vtotal (km/s) TOF (days) Computation time (s) Algorithm 1 10.059794 2703 1.878261 Algorithm 2 (3-5) 9.197015 2954 74.539424 Algorithm 3 (k=150, q=0.3) 9.080567 2846 426.528380 Table 6.11 Comparison of the optimal solution for the mission to Saturn via GAM around Mars obtained using each algorithm. In this case, each algorithm improves the result obtained by the previous one. Algorithm 2 clearly outperforms Algorithm 1, allowing to achieve a total impulse 9 % lower. Algorithm 3 performs better than Algorithm 2, slightly reducing the total impulse on an extra 1,2 % although the computation time is 6 times longer. The optimal solution found using Algorithm 3 is presented on table 6.12 and figure 6.8 Departure Date Arrival Date GAM Date ∆Vini (km/s) ∆Vend (km/s) ∆VGAM (km/s) 07/09/2025 26/03/2033 20/09/2027 4.708733 4.371782 5.2330 ·10−5 hmars (km) δ(º) 290.1054 16.477 Table 6.12 Optimal solution for the mission to Saturn via GAM around Mars obtained using Algorithm 3. The results obtained show better results with a Mars flyby than with a Jupiter flyby, since not only the mission’s cost has been reduced but also the time of flight required. The initial impulse has been diminished by 35%due to the fact that Mars is closer to the Earth than Jupiter. Mars flyby + Jupiter flyby At this point, it is time to include multiple Gravity Assist Maneuvers in the mission. Having studied both the individual fly by Jupiter and by Mars, the first idea which comes to our mind is a mission involving a couple of 42 Chapter 6. Simulation Scenarios and Optimization Results Figure 6.8 Optimal solution for the mission to Saturn via GAM around Mars using Algorithm 3. flybys, the first one around Mars and the second one around Jupiter. The initial guesses for the departure date and the date of the Mars flyby are taken from the optimal solution of the previous problem, 02/11/2025 and 12/01/2028 respectively. Having the required time of flight for a Hohmann transfer between Mars and Jupiter, approximately 3 years, the initial guess for the date of the fly by Jupiter is set, 01/01/2031. Finally, as the time of the Hohmann transfer from Jupiter to Saturn is 10 years, the initial guess for the arrival date is 01/01/41. Then, it must be established the range of possible arrival dates to apply Algorithms 2 and 3. Since the approximate total flight time adding all the Hohmann’s times of flight is 14 years, the range set is (01/01/2035, 01/01/2050). ∆Vtotal (km/s) TOF (days) Computation time (s) Algorithm 1 14.550293 5538 1.488648 Algorithm 2 (3-8) 7.74480682 5699 270.368512 Algorithm 3 (k=96, q=0.3) 7.744707 5700 423.423840 Table 6.13 Comparison of the optimal solution for the mission to Saturn via 2 GAMs (around Mars and Jupiter) obtained using each algorithm. Algorithms 2 and 3 definitely perform better than Algorithm 1 in this case, achieving a solution with half of the cost associated with the solution found by Algorithm 1, which does not reduce the cost of the direct mission. Departure Date Arrival Date GAM♂Date GAMXDate ∆Vini (km/s) ∆Vend (km/s) 16/02/2025 25/09/2040 11/07/2026 12/12/2031 4.048778 3.686066 ∆VGAM♂(km/s) ∆VGAMX(km/s) hp♂(km) rpX/RXδ♂(º)δX(º) 3.605 ·10−49.5016 ·10−3301.7373 10.1 11.806 105.917 Table 6.14 Optimal solution for the mission to Saturn via 2 GAMs (around Mars and Jupiter) obtained using Algorithm 3. 6.2 Mission to Saturn 43 Figure 6.9 Optimal solution for the mission to Saturn via 2 GAMs (around Mars and Jupiter) using Algorithm 3. The optimal solution found by Algorithm 3 improves the ones which were calculated using just one GAM in terms of the cost, even though it requires a higher time of flight. Moreover, the total impulse entails a 35 % reduction of the cost of the direct mission. It seems that a combination of two flybys, one around a planet close to the departure one and another one around a planet "close" to the destiny involves a reduction of both the initial and the final impulse, respectively. Mars flyby + Earth flyby This time, a mission to Saturn which involves a Gravity Assist Maneuver around Mars and then another one around the Earth is analysed. The problem is divided to obtain an initial guess for Algorithm 1. First, a mission with a fly by Mars which then comes back to the Earth is analysed. Then, a mission from Mars to Saturn via an Earth flyby is studied, using as initial guess the optimal arrival date at Mars of the previous mission, then adding the time of flight of the Hohmann trajectories between Mars and the Earth, and the one associated with an Earth-Saturn trajectory. The results obtained for both missions are shown on table 6.15. Departure Date Arrival Date GAM Date TOF (days) Computation time (s) Earth-Mars-Earth 15/03/2025 29/08/2029 22/08/2026 1627 4.135338 Mars-Earth-Saturn 10/08/2026 16/06/2033 02/07/2027 2502 1.502232 Table 6.15 Optimal solution for previous missions required. After having examined these missions, it is possible to set an approximate initial guess for the complete problem using Algorithm 1. The initial date for the initial guess is 15/03/2025, the date of the Mars flyby, 15/08/2026, the date of the second flyby, 01/01/2028, and the arrival date at Saturn, 01/01/2034. For the final date, as the approximate total time of flight adding the time of every Hohmann transfer is 7.5 years, the range set is (01/01/2030, 01/01/2040). As it is shown on table 6.16, once again Algorithm 2 outperforms Algorithm 1 and it is surpassed by Algorithm 3. The optimal solution achieved by Algorithm 3 reduces the cost of the direct mission by 21 % , yet the time of flight is increased by only 810 days. All the details of that optimal solution are shown on table 6.17 and figure 6.12. 44 Chapter 6. Simulation Scenarios and Optimization Results Figure 6.10 Optimal solution for the mission Earth-Mars-Earth using algorithm 1. Figure 6.11 Optimal solution for the mission Mars-Earth-Saturn using algorithm 1. ∆Vtotal (km/s) TOF (days) Computation time (s) Algorithm 1 11.221134 3139 3.433734 Algorithm 2 (3-8) 9.700413 4735 305.911223 Algorithm 3 (k= 96, q=0.6) 9.383823 2897 426.168110 Table 6.16 Comparison of the optimal solution for the mission to Saturn via 2 GAMs (around Mars and around the Earth) obtained using each algorithm. Departure Date Arrival Date GAM♂Date GAM⊕Date ∆Vini (km/s) ∆Vend (km/s) 07/03/2025 10/02/2033 15/08/2026 09/06/2028 4.162791 4.863154 ∆VGAM♂(km/s) ∆VGAM⊕(km/s) hp♂(km) hp⊕(km) δ♂(º)δ⊕(º) 3.103 ·10−50.357846 267.0842 1074.5597 11.358 105.917 Table 6.17 Optimal solution for the mission to Saturn via 2 GAMs (around Mars and around the Earth) calculated using Algorithm 3. Mars flyby + Earth flyby + Jupiter flyby Following, a mission which involves 3 Gravity Assist Maneuvers (Mars, the Earth and Jupiter) is studied. In order to apply Algorithm 1, it is possible to use the optimal solution of the problem with two flybys (Mars and Earth) as initial guess for the initial time (06/03/2024), the time of the Mars flyby (31/08/2026) and the date of the Earth flyby (23/07/2031). For the rest of the variables, a mission from the Earth to Saturn via a Jupiter flyby is analysed, using (23/07/2031) as the initial value for the departure time for Algorithm 1. As that mission have already been studied, even though it was analysed for different dates, the required time of flight for that evaluated mission is used to estimate the date of the fly by Jupiter (adding 1800 days to the date of the Earth flyby) and the arrival date (adding 3100 days to the date of the Jupiter flyby). Then, with respect to Algorithms 2 and 3, the range of the possible arrival dates is set, knowing that the total time adding the Hohmann transfer’s time of flight of every trajectory involved is approximately 14 years, as (01/01/2035, 01/01/2050). The best solution in terms of the cost is reached by Algorithm 3, decreasing the total impulse of the direct mission by 22.5 % . That solution is presented on table 6.19 and figure 6.13. The time of flight of this solution has been remarkably increased, as it is more than 18 years longer than the time of the direct transfer. 6.2 Mission to Saturn 45 Figure 6.12 Optimal solution for the mission to Saturn via 2 GAMs (around Mars and around the Earth) calculated using Algorithm 3. ∆Vtotal (km/s) TOF (days) Computation time (s) Algorithm 1 12.190059 7595 1.847411 Algorithm 2 (3-8) 9.471824 8267 213.900810 Algorithm 3 (k= 80, q=0.3) 9.242740 8920 443.750960 Table 6.18 Comparison of the optimal solution for the mission to Saturn via 3 GAMs (around Mars, the Earth and Jupiter) obtained using each algorithm. Departure Date Arrival Date GAM♂Date GAM⊕Date GAMXDate ∆Vini (km/s) 16/03/2024 17/08/2048 22/09/2026 03/10/2033 25/10/2041 5.077950 ∆Vend (km/s) ∆VGAM♂(km/s) ∆VGAM⊕(km/s) ∆VGAMX(km/s) hp♂(km) hp⊕(km) 3.942244 8.403 ·10−41.675 ·10−40.221538 288.4388 1074.6673 rpX/RXδ♂(º)δ⊕(º)δX(º) 10.10 6.105 40.914 80.704 Table 6.19 Optimal solution for the mission to Saturn via 3 GAMs (around Mars, the Earth and Jupiter) calculated using Algorithm 3. Venus flyby + Mars flyby + Jupiter flyby Subsequently, another mission which involves 3 Gravity Assist Maneuvers (Venus, Mars and Jupiter) is studied. To apply Algorithm 1, the time of flight of the Hohmann transfer from the Earth to Venus (146 days) and from Venus to Mars (217 days) are calculated. To determine an initial solution for the departure time, the direct transfer Earth-Venus is analysed. The optimal initial date found is 01/01/2025. Then, adding the respective times of flight, it is calculated the date of the Venus flyby (01/06/2025), the date of the Mars flyby (01/02/2026), the date of the Jupiter flyby (01/02/2029) and the arrival date (01/01/2040). The range of the possible final dates is lying between 01/01/2035 and 01/01/2050, since the total time adding the Hohmann transfer’s time of flight of every trajectory involved is approximately 14 years. 46 Chapter 6. Simulation Scenarios and Optimization Results Figure 6.13 Optimal solution for the mission to Saturn via 3 GAMs (around Mars, the Earth and Jupiter) calculated using Algorithm 3. ∆Vtotal (km/s) TOF (days) Computation time (s) Algorithm 1 17.429050 5477 2.334620 Algorithm 2 (5-10) 8.521050 5240 441.926640 Algorithm 3 (k= 100, q=0.6) 7.986568 4939 555.826420 Table 6.20 Comparison of the optimal solution for the mission to Saturn via 3 GAMs (around Venus, Mars and Jupiter) obtained using each algorithm. In this case, Algorithm 1 does not perform well at all, because the initial guess may be inaccurate. However, using Algorithm 3 an exceptional solution has been found, reducing by 33 % the cost of the direct mission. That solution is presented on table 6.21 and figure 6.14. Departure Date Arrival Date GAM♀Date GAM♂Date GAMXDate ∆Vini (km/s) 22/03/2025 30/09/2038 07/09/2025 05/03/2028 30/04/2030 4.081123 ∆Vend (km/s) ∆VGAM♀(km/s) ∆VGAM♂(km/s) ∆VGAMX(km/s) hp♀(km) hp♂(km) 3.893304 3.490 ·10−41.707 ·10−30.010086 536.9137 667.2086 rpX/RXδ♀(º)δ♂(º)δX(º) 10.0992 40.479 3.736 94.088 Table 6.21 Optimal solution for the mission to Saturn via 3 GAMs (around Venus, Mars and Jupiter) calculated using Algorithm 3. An important factor which leads to the reduction of the total impulse is the decrease of the initial impulse, due to the proximity of Venus to the Earth. By the same token, the GAM around Jupiter allows to reduce the final impulse. 6.2 Mission to Saturn 47 Figure 6.14 Optimal solution for the mission to Saturn via 3 GAMs (around Venus, Mars and Jupiter) calculated using Algorithm 3. Venus flyby + Earth flyby + Jupiter flyby Finally, a mission which involves 3 Gravity Assist Maneuvers (Venus, the Earth and Jupiter) is analysed. To apply Algorithm 1, the optimal solution from the Earth to Venus previously calculated is employed. The optimal initial date found is 01/01/2025. Then, adding the respective times of flight, the following dates are calculated: the date of the Venus flyby (1/06/2025), the date of the Earth flyby (01/01/2026), the date of the Jupiter flyby (01/01/2029) and the arrival date (01/01/2039). After that, the range of the possible arrival dates is set as (01/01/2035, 01/01/2050). ∆Vtotal (km/s) TOF (days) Computation time (s) Algorithm 1 13.322511 5108 3.760382 Algorithm 2 (5-10) 10.108138 6082 431.139110 Algorithm 3 (k= 100, q=0.6) 8.883265 8408 641.068890 Table 6.22 Comparison of the optimal solution for the mission to Saturn via 3 GAMs (around Venus, the Earth and Jupiter) obtained using each algorithm. The results shown on table 6.22 perfectly define the utility of each algorithm. In order to achieve a better solution, it is necessary to spend a higher computation time. In this case, that time is not extremely longer for Algorithm 3 than for Algorithm 2. The solution calculated by Algorithm 3 is presented on table 6.23 and figure 6.15. The maneuvers included in this mission allow to reduce the cost of the direct mission by 25 % although the time of flight is remarkably increased. Again, the impulse at Jupiter’s periapsis reduces the final impulse. 54 Chapter 6. Simulation Scenarios and Optimization Results 6.3 Mission to Jupiter Eventually, the analysis is focused on a mission to Jupiter. Jupiter is, by far, the largest planet in the Solar System and it is surrounded by dozens of moons. Pioneer 10 was the first spacecraft to visit the Gas Giant, when it flew by Jupiter in 1973. Recently, in 2016, NASA’s mission Juno entered an orbit around Jupiter in order to study its formation and evolution, as well as to explore the planet’s deep atmosphere. The mission evaluated in this section departs from a parking orbit at an altitude hpark1=200 km over the Earth and arrives at another parking orbit at an altitude hpark2=9RX over Jupiter. It is important to remind that the symbol X refers to the planet Jupiter. 6.3.1 Direct Mission First, the direct transfer is evaluated. As there are no references to set a range of possible initial dates for this mission, there has been carried out an investigation to find the optimal departure date. Initial dates from 01/01/2021 to 01/01/2026 have been tested and it has been found an optimal solution on the second half of 2025. Since the Earth-Jupiter synodic period is Tsyn =399 days, the range of the departure time Tini is lying between 01/03/2025 and 01/09/2026. To estimate the necessary time of flight to go from the Earth to Jupiter, the time of flight associated with the Hohmann transfer is calculated (about 998 days, approximately 3 years). Then, the range of the arrival date is established between 01/01/2027 and 01/01/2030. Figure 6.18 Pork chop plot: Direct mission to Jupiter. After having built the pork-chop plot associated with this mission (figure 6.18), an optimal solution is detected. It involves a departure date around the month of October of 2025 and an arrival date around August 6.3 Mission to Jupiter 55 of 2028. Then, fmincon is employed to search the optimal point shown in the pork-chop plot, operating with the initial guess ( Tini =01/10/2025 , Tend =01/08/2028 ). The solver successfully finds the optimal point shown in figure 6.4, which is detailed on table 6.35 and figure 6.19: Departure Date Arrival Date ∆Vtotal (km/s) ∆Vini (km/s) ∆Vend(km/s) Direct Mission 25/09/2025 16/08/2028 12.664042 6.406980 6.257062 TOF (days) Computation time (s) 1056 2.486128 Table 6.35 Optimal solution for the direct mission to Jupiter. Figure 6.19 Optimal solution for the direct mission to Jupiter. The optimal solution found is similar to a Hohmann transfer and the associated cost seems to be very expensive, thus it must find an alternative way to reach Jupiter. 6.3.2 Gravity Assist Maneuvers The cost of the direct mission is considerably high hence, in view of the results of the previous section, it was decided to include Gravity Assist Maneuvers in the mission. It was expected that the time of flight would also be increased, yet it would not be a serious problem. The range of the initial date is broadly establish as (01/06/2024, 01/06/2026). Mars flyby Firstly, a swing by Mars is evaluated. In Section 6.2, an optimal solution for a direct Earth-Mars transfer was found in May 2025, where the arrival date at Mars takes place in July 2027. That solution is used as the initial guess for Algorithm 1, calculating the final date adding the time of flight of the Hohmann transfer of Mars-Jupiter (3 years). The lower bound associated with the final date is set as 01/01/2028, and its upper bound is 01/01/2032. Algorithms 2 and 3 find similar solutions, outperforming Algorithm 1. The optimal solution found using Algorithm 3 reduces the cost of the direct mission by 21.5 % , and it is presented on table 6.37 and figure 6.20 56 Chapter 6. Simulation Scenarios and Optimization Results ∆Vtotal (km/s) TOF (days) Computation time (s) Algorithm 1 10.612632 2080 1.969129 Algorithm 2 (3-5) 9.943996 1694 114.962267 Algorithm 3 (k=150, q=0.3) 9.938240 1654 406.142325 Table 6.36 Comparison of the optimal solution for the mission to Jupiter via GAM around Mars obtained using each algorithm. Departure Date Arrival Date GAM Date ∆Vini (km/s) ∆Vend (km/s) ∆VGAM (km/s) 08/01/2025 20/07/2029 02/05/2026 3.896029 6.042191 1.938 ·10−5 hmars (km) δ(º) 285.7283 15.482 Table 6.37 Optimal solution for the mission to Jupiter via GAM around Mars obtained using Algorithm 3. Figure 6.20 Optimal solution for the mission to Jupiter via GAM around Mars using Algorithm 3. Once again, the use of a Gravity Assist Maneuver allows to reduce the total cost of the mission. In this case, the initial impulse is strongly decreased due to the proximity of Mars to the Earth. However, the final impulse is still high. Venus flyby + Mars flyby Secondly, a GAM around Venus is added before the Mars flyby. Having discovered an optimal solution on January 2025 for the direct mission from the Earth to Venus, it is used as the initial guess for the departure date. Then, obtaining the respective times of flight associated with each Hohmann trajectory, the initial guess for the rest of the dates is determined. Since the estimated total time of flight, calculated using Hohmann tranfers, is about 4 years, the range of the arrival dates lies between 01/01/2028 and 01/01/2033. Algorithms 2 and 3 perform better than Algorithm 1, although they need a large computation time, especially Algorithm 3, since 200 ants are needed to find that solution. The solution discovered by Algorithm 3, represented on table 6.39 and figure 6.21, reduces the total cost of the direct mission by 17 % , adding 500 days to the total time of flight. 6.3 Mission to Jupiter 57 ∆Vtotal (km/s) TOF (days) Computation time (s) Algorithm 1 11.992899 1433 6.473917 Algorithm 2 (5-10) 10.581271 2159 522.631001 Algorithm 3 (k=200, q=0.3) 10.455651 1547 2181.518370 Table 6.38 Comparison of the optimal solution for the mission to Jupiter via 2 GAMs (around Venus and Mars) obtained using each algorithm. Departure Date Arrival Date GAM♀Date GAM♂Date ∆Vini (km/s) ∆Vend (km/s) 03/04/2025 29/06/2029 26/09/2025 22/05/2026 4.430988 6.018253 ∆VGAM♀(km/s) ∆VGAM♂(km/s) hp♀(km) hp♂(km) δ♀(º)δ♂(º) 1.048 ·10−35.361 ·10−3363.7370 286.9712 32.630 20.315 Table 6.39 Optimal solution for the mission to Jupiter via 2 GAMs (around Venus and Mars) calculated using Algorithm 3. Figure 6.21 Optimal solution for the mission to Jupiter via 2 GAMs (around Venus and around Mars) calculated using Algorithm 3. Venus flyby + Earth flyby Next, the mission with a first GAM around Venus and a second one around the Earth is analysed. The same initial guess as in the previous case is used for the initial and the Venus flyby dates. Then, the times of flight of the respective direct transfers are added. The departure date’s range is set as (01/01/2028, 01/01/2032). ∆Vtotal (km/s) TOF (days) Computation time (s) Algorithm 1 11.982431 1829 3.447780 Algorithm 2 (3-8) 10.328592 1770 326.222310 Algorithm 3 (k=200, q=0.3) 10.255098 1608 1542.5592 Table 6.40 Comparison of the optimal solution for the mission to Jupiter via 2 GAMs (around Venus and the Earth) obtained using each algorithm. 58 Chapter 6. Simulation Scenarios and Optimization Results This time, the computation time required by Algorithm 3 to improve the solution found by Algorithm 2 is almost 5 times higher, since 200 ants are needed. However, that solution, presented on table 6.41 and figure 6.22, achieves a total cost 19%lower than the direct transfer, considerably reducing the initial impulse. Departure Date Arrival Date GAM♂Date GAM⊕Date ∆Vini (km/s) ∆Vend (km/s) 07/03/2025 01/08/2029 20/08/2025 22/09/2026 3.933149 6.316399 ∆VGAM♂(km/s) ∆VGAM⊕(km/s) hp♂(km) hp⊕(km) δ♂(º)δ⊕(º) 3.178 ·10−32.372 ·10−3636.7420 1017.0090 50.278 33.759 Table 6.41 Optimal solution for the mission to Jupiter via 2 GAMs (around Venus and around the Earth) calculated using Algorithm 3. Figure 6.22 Optimal solution for the mission to Jupiter via 2 GAMs (around Venus and around the Earth) calculated using Algorithm 3. Mars flyby + Earth flyby Subsequently, the mission with a swing by Mars and another one by the Earth is analysed. The optimal solution previously found for the mission with a GAM around Mars is employed to set the initial and the Mars flyby dates. The rest of the dates are calculated adding the times of flight of the Hohmann transfer of each phase. The range of possible arrival dates is established as (01/01/2028,01/01/2035), since the total estimated Hohmann’s time of flight is slightly superior to 4 years. ∆Vtotal (km/s) TOF (days) Computation time (s) Algorithm 1 13.125564 2341 2.849000 Algorithm 2 (3-8) 11.869616 1994 203.541980 Algorithm 3 (k=400, q=0.3) 10.649167 1443 5485.7480 Table 6.42 Comparison of the optimal solution for the mission to Jupiter via 2 GAMs (around Mars and the Earth) obtained using each algorithm. Finding a great optimal solution was hard for Algorithm 3, with a computational time which went through the roof. This time, Algorithm 1 could not obtain a solution which enhances the one associated with the 6.3 Mission to Jupiter 59 direct mission. The solution achieved by Algorithm 3, shown on table 6.43 and figure 6.23, lessens the total cost of the direct mission by 16%, increasing the time of flight by 400 days. Departure Date Arrival Date GAM♂Date GAM⊕Date ∆Vini (km/s) ∆Vend (km/s) 07/10/2024 19/09/2028 30/03/2025 18/11/2025 4.217310 6.428591 ∆VGAM♂(km/s) ∆VGAM⊕(km/s) hp♂(km) hp⊕(km) δ♂(º)δ⊕(º) 1.701 ·10−46.696 ·10−3286.9468 1051.7345 21.612 20.177 Table 6.43 Optimal solution for the mission to Jupiter via 2 GAMs (around Mars and around the Earth) calculated using Algorithm 3. Figure 6.23 Optimal solution for the mission to Jupiter via 2 GAMs (around Mars and around the Earth) calculated using Algorithm 3. Venus flyby + Earth flyby + Mars flyby Lastly, a mission which involves three flybys (Venus, the Earth and Mars) is analysed. The solution for the direct mission to Venus is employed as the initial guess, adding the required times of flight until reaching Jupiter. The arrival date is lying between 01/01/2028 and 01/01/2033. ∆Vtotal (km/s) TOF (days) Computation time (s) Algorithm 1 17.544543 1624 8.036202 Algorithm 2 (3-8) 15.995125 4839 468.122610 Algorithm 3 (k=200, q=0.3) 11.534919 1733 1745.478800 Table 6.44 Comparison of the optimal solution for the mission to Jupiter via 3 GAMs (around Venus, the Earth and Mars) obtained using each algorithm. Only Algorithm 3 is able to reach a solution with a lower total cost than the direct mission to Jupiter. Although that solution reduces the cost by 9% , it requires a high initial impulse, as it is presented on table 6.45 and figure 6.24. 60 Chapter 6. Simulation Scenarios and Optimization Results Departure Date Arrival Date GAM♀Date GAM⊕Date GAM♂Date ∆Vini (km/s) 30/06/2024 30/03/2029 18/02/2025 06/05/2025 10/10/2026 5.374171 ∆Vend (km/s) ∆VGAM♀(km/s) ∆VGAM⊕(km/s) ∆VGAM♂(km/s) hp♀(km) hp⊕(km) 6.157258 1.235 ·10−41.719 ·10−31.648 ·10−3370.7752 1085.9856 hp♂(km) δ♀(º)δ⊕(º)δ♂(º) 300.8328 78.679 81.578 11.447 Table 6.45 Optimal solution for the mission to Jupiter via 3 GAMs (around Venus, the Earth and Mars) calculated using Algorithm 3. Figure 6.24 Optimal solution for the mission to Jupiter via 3 GAMs (around Venus, the Earth and Mars) calculated using Algorithm 3. 6.3 Mission to Jupiter 61 Final comparison In conclusion, the results reached by each algorithm are included in a table, in order to compare them. ∆Vtotal (km/s) TOF (days) ∆Vini (km/s) ∆Vend (km/s) Direct 12.6640 1056 6.4070 6.2571 Mars flyby 10.6126 2080 .2997 6.3126 Venus + Mars 11.9929 1433 4.0185 6.4396 Venus + Earth 11.9824 1829 4.0119 6.6568 Mars + Earth 13.1256 2341 4.4825 6.9572 Venus + Earth + Mars 17.5445 1624 6.4084 6.4303 Table 6.46 Optimal solution for every GAM mission to Jupiter using Algorithm 1. ∆Vtotal (km/s) TOF (days) ∆Vini (km/s) ∆Vend (km/s) Direct 12.6640 1056 6.4070 6.2571 Mars flyby 9.9440 1694 3.9002 6.0436 Venus + Mars 10.5813 2159 4.1637 6.4132 Venus + Earth 10.3276 1770 3.9656 6.3559 Mars + Earth 11.8696 1994 5.3856 6.4828 Venus + Earth + Mars 15.9951 4839 6.9807 8.8102 Table 6.47 Optimal solution for every GAM mission to Jupiter using Algorithm 2. ∆Vtotal (km/s) TOF (days) ∆Vini (km/s) ∆Vend (km/s) Direct 12.6640 1056 6.4070 6.2571 Mars flyby 9.9382 1654 3.8960 6.0422 Venus + Mars 10.4557 1547 4.4310 6.0183 Venus + Earth 10.2551 1608 3.9331 6.3164 Mars + Earth 10.6998 1443 4.2173 6.4286 Venus + Earth + Mars 11.5349 1733 5.3742 6.1573 Table 6.48 Optimal solution for every GAM mission to Jupiter using Algorithm 3. In brief, every algorithm is able to calculate a solution with a lower cost than the direct mission. Nevertheless, Algorithm 3 and Algorithm 2 improve the direct mission for every case under consideration except for the 3 GAMs problem. The major obstacle of the mission to Jupiter is the final impulse, which has not been reduced enough in any of the missions. That point is the main difference between this case and the mission to Saturn, which allowed to reduce the cost of the direct mission by 35 % using Gravity Assist Maneuvers. In light of the results shown on table 6.48, the optimal solution to reach Jupiter seems to be a fly by Mars, which reduces the cost by 21.5 % adding 638 days to the direct mission. Summarising, a mission to Saturn is cheaper than a mission to Jupiter, yet the time of flight required for the first one is extremely higher. 7 Conclusions and Future Work 7.1 Conclusions This project has presented a novel approach to the preliminary design of interplanetary trajectories. The aim of this study was to explore smart ways of using local-optimization based methods in order to perform a global minimum search. This work has demonstrated the efficiency of such approach both in terms of local optimal solution improvement as well as maintaining a reasonable low computational burden (which is usually not the case for Metaheuristic algorithms). Three algorithms of various complexity have been tested. Algorithm 1 yielded significant speed advances when the problem under consideration was not excessively complicated. Nonetheless, it was actually dependent on the previous analysis. Algorithm 2 has performed pretty well, providing some interesting solutions and allowing to discover unexpected ones although it requires a longer computation time. It should be noted that the use of a multi-start search employing local optimization processes which start from initial guesses randomly distributed over the search space, made it possible to identify new optimal solutions. Finally, the most important advantage of Algorithm 3 is that it can enhance the solutions which were found by Algorithm 2, adding a back-up procedure which allows to focus on the best candidates. Definitely, the Ant-Colony optimization procedure allows to achieve better solutions, deeply exploring the search space and concentrating on the most auspicious contenders. Summarising, the main achievement of the proposed algorithms is the fact that they achieved a significantly beneficial solution which appears to be a global minimum, although a local optimization procedure is employed. Moreover, the use of a local optimal searcher led to a noticeable low computation time, which is less than 30 minutes in the worst case. Four mission analysis classes have been investigated: the two impulse direct planet-to-planet transfer problem, the Deep Space Maneuver problem, the Multiple Gravity Assist interplanetary transfer problem and the problem which involves a Gravity Assist and Deep Space Maneuvers. By considering the objective function value reached at the end of the optimization process and the required computation time to get it, the results of the test phase can be summarised as: • Direct transfer: The optimal solution found is always a kind of Hohmann transfer when trajectories with complete revolutions around the Sun are not evaluated. The pork chop plot shows one clear optimal point in terms of fuel consumption, due to the fact that the 2-D problem is considered. The required computation time is very short, and every algorithm can discover the same solution. • DSM problem: It has been proved in the mission to Mars that a DSM can reduce the necessary total cost. Generally, they help to reduce the final impulse and they increase the required total time of flight indeed. Moreover, the investigation of this problem allows to discover new optimal solutions for the direct one. Since this method splits the trajectory, it is possible to find optimal transfers with complete revolutions around the Sun where the DSM are zero which enhance the previously discovered optimal direct trajectory. • MGA problem: The results demonstrate that the proposed methodology can effectively tackle the complexity of multiple Gravity Assist interplanetary transfers while avoiding being trapped in local minima. Furthermore, it was proved that Gravity Assist Maneuvers permit to substantially decrease the 63 70 Bibliography [18] R.A.N. Araujo, O.C. Winter, A.F.B.A. Prado and R. Vieira Martins, Sphere of Influence and Gravitational Capture Radius: A Dynamical Approach, Monthly Notices of the Royal Astronomical Society, vol. 391, 675-684, 2008. [19] Sam Wagner, Bong Wie and Brian Kaplinger, Computational Solutions to Lambert’s Problem on Modern Graphics Processing Units, Journal of Guidance Control and Dynamics, vol. 38, 1305-1310, 2015.