3-Dimensional Continuous-Thrust Trajectory Optimization for Electric Solar Wind Sails using Bezier Curves

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 Grado en Ingeniería Aeroespacial 3-Dimensional Continuous-Thrust Trajectory Optimization for Electric Solar Wind Sails using Bezier Curves Autor: Miguel García Ureña Tutor: Guillermo Pacheco Ramos Dpto. Ingeniería Aeroespacial y Mecánica de Fluidos Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2025 Trabajo Fin de Grado Grado en Ingeniería Aeroespacial 3-Dimensional Continuous-Thrust Trajectory Optimization for Electric Solar Wind Sails using Bezier Curves Autor: Miguel García Ureña Tutor: Guillermo Pacheco Ramos Profesor Ayudante Doctor Dpto. Ingeniería Aeroespacial y Mecánica de Fluidos Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2025 Trabajo Fin de Grado: 3-Dimensional Continuous-Thrust Trajectory Optimization for Electric Solar Wind Sails using Bezier Curves Autor: Miguel García Ureña Tutor: Guillermo Pacheco Ramos 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: Agradecimientos E n los últimos 4 años de mi vida no solo me he formado académicamente sino que también personalmente. No hace mucho ponía por primera vez un pie en esta escuela de ingenieros lleno de ilusión y ganas por dejar mi huella de alguna forma u otra, no tardé en darme cuenta que todo iba a ser muy distinto de ese momento en adelante. Me gustaría dedicar este proyecto a todas esas personas que me han acompañado en mi camino durante esta etapa, especialmente a las increíbles personas que me llevo de esta carrera, desde el principio hasta ahora Blasco, Agustín, Álvaro, Emilio, Tachi, Lucía, Calero, Collado y Fernando, sin ellos nada habría sido lo mismo, ayudándonos en esas noches interminables antes de un examen o simplemente un “aguanta que podemos”, sinceramente no habría podido llegar a donde estoy ahora sin ellos, muy agradecido porque compartan mi “locura”. Quiero hacer especial mención a mi compañera de camino durante la mitad del trayecto, muchas gracias Carmen, por recordarme quién soy y que puedo hacer cuando no veo otra alternativa o simplemente recordarme que soy capaz de lograr lo que me proponga, también por darme el ejemplo de siempre elegir lo que te hace feliz y no lo que se tendría que elegir. También quiero agradecer enormemente a mi familia, sin ellos simplemente nada habría sido como ha sido, gracias por ayudarme en todos los aspectos de mi vida y siempre apoyarme sin importar la decisión que tomara fuera más acertada o menos, gracias a mi hermana Elena, por ser el ejemplo más cercano de que siempre hay que seguir hacia adelante sin importar los contratiempos. Finalmente, quiero expresar mi enorme gratitud a mi tutor, Guillermo, juntos hemos llevado a cabo este proyecto a veces muy complaciente, otras veces un poco frustrante, pero siempre guiándome y ayudándome en los momentos que lo he necesitado, ojalá todo el mundo pusiera las ganas y la pasión como él. Espero que este proyecto cumpla con las expectativas del lector, pues, las mías las ha superado con creces, y que se haya quedado una pequeña parte de mí reflejada en él, siempre con Voluntad, Agradecimiento, Ganas, Ilusión, Amor y Sueños. Miguel García Ureña Sevilla, 2025 I Resumen Título Optimización de Trayectorias Interplanetarias para Velas Eléctricas Solares con Curvas de Bezier Resumen del trabajo En este Trabajo Fin de Grado se presenta un método de optimización de trayectorias de transferencia interplanetaria empleando una vela solar eléctrica (E-sail). Partiendo de un modelo dinámico que incluye las fuerzas del Sol, se formula el problema de control óptimo como una Programación No Lineal (NLP). También, la trayectoria continua se discretiza mediante curvas de Bézier, lo que permite transformar la optimización en un problema de dimensión finita y reducir significativamente el tiempo de cómputo optimizando una seria de coeficientes geométricos del desarrollo polinómico. Los resultados demuestran la eficacia de la aproximación con Bézier en términos de precisión y coste computacional, así como la factibilidad de las soluciones obtenidas para distintas ventanas de lanzamiento. Palabras clave Vela solar eléctrica (E-sail), Programación no lineal (NLP), Curvas de Bézier, Órbita interplanetaria, Optimización de trayectorias, Problema de control óptimo (OCP), Misiones de rendezvouz. Conclusiones El objetivo es demostrar que las curvas de Bezier recrean con precisión las trayectorias de velas solares eléctricas mientras reducen los tiempos de cómputo entre diez y cien veces frente a la transcripción directa; se validan transferencias óptimas a Marte, Júpiter y varios asteroides, evidenciando robustez tanto para órbitas planetarias como para trayectorias muy excéntricas; el análisis de la aceleración característica confirma que los mayores ahorros de tiempo de vuelo se alcanzan en el régimen de bajo empuje, aportando una regla práctica para dimensionar la vela; el barrido de orden polinómico identifica grados 8–12 como equilibrio ideal entre exactitud y coste computacional; las trayectorias suaves obtenidas facilitan la implementación de los controladores LQR y MPC propuestos como trabajo futuro para mitigar la variabilidad del viento solar; la metodología se integra con las efemérides calculadas con SPICE y admite reoptimización en vuelo gracias a su baja dimensionalidad; en conjunto, este proyecto ofrece una herramienta rápida y fiable para el diseño preliminar de misiones para Velas Eléctricas Solares y sienta las bases para futuras extensiones con asistencias gravitatorias y control predictivo. III XNotation RTN Radial-Transverse-Normal YSaturn σSail “clocking” angle. SPICE Spacecraft Planet Instrument C-matrix Events Toolkit. tfFinal time of the transfer. τNormalized time parameter. τcCylindrical reference frame ˜aρ,˜aθ,˜azDimensionless propulsion acceleration components. ˜ ρ,˜ θ,˜zDimensionless position coordinates. ˜vi Corresponding dimensionless velocity components for (˜ ρ,˜ θ,˜z). uΦUnit vector in the φcoordinate uρUnit vector in the ρcoordinate uzUnit vector in the zcoordinate ZUranus UT Time Unit. UV Velocity Unit. vρ,vθ,vzVelocity components in ρ,θ,zdirections. ϕAzimuthal angle in cylindrical coordinates. ♀Venus xvector of optimization variables. ucontrol vector (κ,αn,σ). 1 Introduction "I have loved the stars so passionately that I do not fear the night." Miguel F rom the beginning of time, the sky has been the object of enormous admiration, from the first civilisations to the present day, from the first bright points observed by mankind to the most exotic events in the universe. All this has made mankind ask more and more questions, making the curiosity to know and master the laws of the universe and its constituent parts, grow. All this has led us to take a look at our nearest neighbors. However, there has always been an inevitable barrier, time. 1.1 Motivation Interstellar travel is a subject that has always fascinated readers, the possibility of travelling to new worlds and becoming explorers of the modern age; however, as good as it sounds, the idea of interstellar travel fades away when the mathematics of how to achieve it is analysed. If we look at our "home", the Milky Way Galaxy, which has a diameter of 105,700 light-years (about 10 ·1017 kilometres). Considering that the greatest distance ever reached by a probe is about 24,000 million kilometers, reached by the Voyager 1, which has been on an escape trajectory from the solar system for 30 years. As can be seen, we still have a long way to go to understand where we live. If we focus on a less ambitious aspect, the closest planetary system to the Earth, "Alpha Centauri", at a distance of 4.25 light-years, (light takes 4.25 light-years to travel that distance at a speed of 1, 079, 252, 849.6 km/h). To make it even clearer, we have that the human-made object that has reached the highest speed is the Parker probe, thanks to a gravity-assisted maneuver using the Sun, which was 692.017 km/h. At this speed, it is estimated that the time of arrival at "Alpha Centauri" would be 8000 years, which suggests that it is not feasible today. As has just been demonstrated, the travel to other planetary systems is highly improbable, so we must wait some hundreds of years till we become explorers again. However, we can focus on what is near us, the Solar System. [12] Since the beginning of the space race, there have been numerous interplanetary missions to visit the planets of the Solar System and understand them better. For example, the "Venera" capsules aimed at Venus, either to photograph the planet and to identify the main components of the atmosphere in case habitability was possible (the conclusions were negative), the "Mariner" probes (1964-1971), which were mainly used by NASA for the arrival and descent to Mars. Jupiter has 1 2Chapter 1. Introduction been the most visited planet with about 7 flybys with "Pioneer 10" (1973) and "Pioneer 11" (1974), "Voyager 1" (1979) and "Voyager 2" (1979), "New Horizons" (2007) and "Juno" (2016-present).[ 21 ] The exterior planets, Neptune and Uranus, have only been visited once by the Voyager probes on an escape trajectory from the solar system, taking advantage of their passage by these planets to perform gravity assistance and gain more speed. All these missions have used traditional propulsion systems, which take advantage of the action-reaction effect to generate thrust by using fuel, but even so, in order to reduce flight times, they have had to use gravity assist maneuvers, which provide extra thrust by taking advantage of the speed at which the planets orbit around the Sun. Although these maneuvers save fuel and time, the main problem with them is that the planets have to be in a specific position with respect to the Earth, so there is a waiting time, which defines the launch windows. This time increases enormously if more than one planet is involved, which is normal, as the energy requirements for deep space missions are very high. For example, the planetary configuration followed by ‘’Voyager 2‘’, passing Jupiter, Saturn, Uranus and Neptune, is repeated every 174 years, which is a very long time. [19] Figure 1.1 Vogayer 1 Probe Propulsion System, showing RTG [15]. A possible solution to this problem is that instead of generating a large increase in speed over a ‘short’ period of time, a continuous propulsion with a small thrust, but which over a longer period of time can reach high speeds. This is why the Low-Thrust Continuous Propulsion System was developed. With these systems it would be possible to reduce launch waiting times or even combine them with gravitational maneuvers, but they have a major limitation: the fuel they use runs out. This is where Low-Thrust Propellantless Continuous Propulsion Systems come in. [14] The main idea behind this is to use the sun as the primary source of energy production, specially, using the solar wind. The different technologies that have been studied are Electric, Magnetic and Photonic sails. To demonstrate the feasibility and performance of E-sail, some missions have been developed, such as ESTCube-1 developed by the University of Tartu, which tested this technology in a LEO orbit, as a first approximation, since the relative speed of the plasma is 7 km/s which is much slower than the solar wind (300-800 km/s) but there is a higher concentration, so it can be used as a preliminary phase before a mission outside the Earth’s magnetosphere to fully demonstrate the performance of this system. Another mission is LightSail 2, which is able to change its orbit only using solar energy , and the Ikaros mission, with a similar objective.[3] 1.2 State of the Art 3 In this project focus will be on Solar Electric Sails, as they are a viable solution to all the problems described above. 1.2 State of the Art The E-sail concept was proposed by Pekka Janhunen in 2004 [ 8 ]. He was a Finnish space physicist, astrobiologist and inventor. Before considering E-sails, Pekka Janhunen was analyzing the magnetic solar wind sail, analyzed by Dana Andrews and Robert Zubrin in 1988. The European Space Agency (ESA) had given a project to Janhunen’s team to analyze whether the magnetic sail was something that the Agency should investigate further. Janhunen wasn’t a space technologist but a leading space plasma physicist instead, which was enough to analyze the solar wind sail concept.[24] Figure 1.2 Electric Solar Wind Sail concept [18]. After that, the ESA study concluded that the magnetic sail works in theory, but its implementation would need a kilometres-long(1-20 km), superconducting cable, which would require either some lightweight and reliable sun-shield design or a near-room-temperature superconducting material. However, deep immersion into solar wind plasma physics and interplanetary propulsion provided Janhunen with an idea. It was based on replacing the heavy electromagnets with thin wires carrying an electric charge. The principle of operation is simple, the propulsion in E-sails is based on the conservation of linear momentum caused by the repulsion of the solar wind plasma, which is intercepted by electrically charged wires. The electric field generated by the wires affects the trajectory of the protons in the solar wind, reducing the component of their momentum in this direction. This loss of momentum by the protons is transferred to the wires, which experience a Coulomb force (the charged wire undergoes a force due to the plasma). This force is then transferred to the probe as thrust. To maintain the high voltage bias on each tether, requires emitting collected electrons back into the deep space media via an electron gun on the spacecraft.[11] The operation differs from other previously considered methods, such as photonic and magnetic sails, which are also based on the conservation of linear momentum but use different mechanisms. Photonic sails rely on the momentum transfer from photons in the solar radiation, while magnetic sails use magnetic fields to deflect the trajectory of incoming protons. The main advantage proposed by E-sails lies in their thrust behavior, which decreases as 1/r , where ris the distance to the Sun. In contrast, the thrust of photonic sails decreases as 1/r2. 4Chapter 1. Introduction This means that for a mission planing for distances far away from the Sun, E-sails might be more appropriate spacecraft since they have a smooth gradient whereas the other decreased quickly. The primary reason for this difference is that the dynamic pressure of the solar wind decreases with 1/r2 , but the effective area of the sail—proportional to the envelope surrounding the wires—increases with r , following the same trend as the Debye length of the plasma. (The Debye length is the distance over which the electric field of a charge is shielded,so the particle can’t go through the wires and produce thrust by bouncing.) As a result, the product of these factors leads to a thrust dependence of 1/r. As previously demonstrated, this system proposes crucial elements for utilizing the solar wind, whose influence extends throughout the entire solar system. It can be used in unconventional missions that would be impossible with traditional propulsion systems, such as polar orbits for solar observation or non-Keplerian orbits, which require continuous thrust corrections to be maintained around the Lagrange points. If we focus on the evolution of the E-sails, the first design was based on a square shape design whose characteristics was 10-100 m2 surface, near 2010 and 2018, but last tendencies, in contrast, lean towards a circular shape instead, with this geometry it will be more easy to deploy and maintain the structure stability by inducing a centrifugal force so the wires are maintained stretchedthis is not possible in a square-shape geometryand also, increasing the effective area of the sail to produce more thrust(1000-10000 m2) and decrease the density of the sail 1−2.5g/m2. Although optimizing E-sails for higher performance is conceivable, their reliance on the variable solar wind can jeopardize critical mission phases. To eliminate this dependency, a ground-based laser can continuously illuminate the sail, enabling sustained thrust regardless of solar wind conditions. By beaming energy from Earth to the spacecraft, the effective sail area could be increased to between 10 000 and 1 000 000 m2while reducing sail areal density to as little as 0.1 g/m2dramatically improving thrust generation and mission reliability [15]. Figure 1.3 Evolution of the E-sails throughout the time [10]. 1.3 Aim of the project 5 1.3 Aim of the project As previously mentioned, Electric sails (E-sails) still have significant potential for advancement before they can reliably reach distant points in space. The main goal of this project, however, is not to investigate the physical mechanisms behind the thrust generated by E-sails, nor their operational behavior. Instead, the emphasis is placed on enhancing the mission planning stage. Specifically, this project aims to streamline the preliminary design phase of space missions by employing sophisticated trajectory optimization techniques. Incorporating Bezier curves represents a promising method, which not only improves computational efficiency but also yields precise trajectory solutions and reduces overall computational time required for optimal transfer paths. Ultimately, this approach will significantly reduce design costs associated with space missions, enabling more efficient and cost effective mission planning. • NLP problem planning and solution: To begin, we establish a comprehensive thrust model for the E-sail, detailing the expected thrust levels and the key control parameters governing spacecraft maneuvering. Building upon this foundation, we formulate the underlying nonlinear programming (NLP) problem to optimize the transfer trajectory. We then solve the NLP via direct transcription methods, discretizing the control and state variables to transform the continuous-time optimization into a finite-dimensional nonlinear program. Using an optimization solver, we obtain time optimal trajectories that satisfy all mission constraints. Finally, we present the optimized transfer solutions for missions to Mars, Jupiter and several asteroids, highlighting their feasibility and performance metrics. • OCP solution using Bezier curves: Several refinements to the baseline problem are introduced in this section, followed by an analysis that shows the Bezier curve formulation can reduce the overall computational time by an order of magnitude. Throughout the project the spacecraft’s attitude is treated as prescribed. Although active research is under way on spin plane control using voltage modulation and other means that couple directly into the sail’s dynamics we model the attitude subsystem as a black box, assuming that the desired orientation can be reached instantaneously at any point along the trajectory. Furthermore, the present analysis proceeds under the following assumptions: – Only the gravitational attraction of the Sun is considered. We have neglected the other planet’s forces. – The launch date can be changed with freedom, always being as close as possible to the real launch windows. 1.4 Structure of the project The layout of this project is sketched in the paragraphs that follow. Each paragraph offers an expanded synopsis of a single chapter, clarifying its specific aims, the methods applied, and the way its results dovetail with the overarching investigation. In this manner, the reader can quickly grasp how the study develops from foundational theory through methodological innovation to the final assessment and appreciate how the individual pieces combine to form a coherent whole. Chapter 1: Introduction. In this opening chapter lays the foundation for the project by first explaining the underlying motivation that drives the research. It then surveys the current state of the art, highlighting key advances and open challenges in the E-sails field. Building on this context, the chapter clearly states the project objectives, which serve as guiding milestones for the ensuing work. 6Chapter 1. Introduction Finally, an overview of the document’s structure is provided, outlining how each subsequent chapter contributes to fulfilling these objectives and leading the reader logically from problem formulation to final conclusions. Chapter 2: Background in Orbital Mechanics. This chapter derives the spacecraft’s fundamental equations of motion and formulates a continuous thrust model for electric sails (E-sails), noting key practical limitations. It concludes by establishing the physical context required for the trajectory analysis that follow. Chapter 3: Optimal Mission Planning. This chapter casts the mission as an optimal control problem, translates it into a nonlinear programming (NLP) format through direct transcription, and presents time and propellant optimal transfer trajectories to Mars, Jupiter, and a representative set of asteroids. Chapter 4: Bezier Curve Based Planning. After outlining the mathematical definition of Bezier curves, this chapter reformulates the transfer problem within that parametric framework and benchmarks the resulting trajectories against the direct method, comparing computational speed and solution feasibility. Chapter 5: Conclusions and Future Work. The closing chapter consolidates the project’s key results, evaluates the overall efficiency of the proposed methods, and suggests directions for extending the approach to future mission scenarios. All of our procedures are based on the problem formulation described in [ 18 ], specially the Bezier problem formulation. 2 Background in orbital mechanics I n this chapter a look at the basics of orbital mechanics is taken, developing a continuous thrust model that shapes well with the motion of the E-sails throughout interplanetary space. The aim of this chapter is to propose a rapid shape based method where the concept of Bezier curve is used to efficiently design the three dimensional interplanetary trajectory of a spacecraft propelled by E-sails. 2.1 Equation of motion 2.1.1 Deduction of the Equation of Motion For this deduction, the Law of Universal Gravitation discovered by Sir Isaac Newton (1642-1727) is used. This law describes the gravitational effect that a massive body experiences when it is near another massive body. When more planets are involved, no analytical solution exists for these equations, so numerical results are required; however, for the problem being solved, the two body model is sufficient. The Law of Universal Gravitation states that the gravitational force exerted by two massive bodies is directly proportional to the product of their masses and inversely proportional to the square of the distance between their centres of mass. The equation is presented in vector form, so the attractive vector force always lies along the line that connects their centres of mass. These are: F12 =Gm1m2 r2 R2−R1 |R2−R1|(2.1) F21 =Gm1m2 r2 R1−R2 |R1−R2|(2.2) 7 8Chapter 2. Background in orbital mechanics Figure 2.1 Description of the gravitational Force between two bodies. where G is the Universal Gravitational Constant, G = 6.67430 ×10−11 m3kg−1s−2 , Fi j is the force in the body i exerted by the body j , as it can be seen in Figure 2.1, and R1 , R2 are the position vectors of each body in an inertial frame of reference, as it can be seen in Figure 2.1. In what follows, the vector rdenotes the direction from mass 1 to mass 2: r=R2−R1,r=|r|(2.3) After that, Newton’s second law is introduced to derive the differential equation of motion. F=m·a(2.4) Using (2.1) and (2.2) we have: m1·¨ R1=F12 =G·m1·m2 r2·r r(2.5) m2·¨ R2=F21 =−G·m1·m2 r2·r r(2.6) Dividing by m1and m2yields the equation in terms of ¨ R1and ¨ R2. ¨ R1=G·m2 r2·r r(2.7) ¨ R2=−G·m1 r2·r r(2.8) Recalling the definition of rin (2.3) ¨ R2−¨ R1=¨ r=−G·m1+m2 r2·r r(2.9) This equation has an analytical solution in the case where only two bodies are involved; in that case, the center of gravity of the system always moves with constant speed, so for convenience, the definition of our inertial reference system is made at this point. To simplify the problem, certain hypotheses can be made. When one body has a significantly greater mass than the other, the centre of mass of the system lies very close to the massive body, so this point can be assumed to coincide with the centre of the massive body. Under this simplification, R1is constant and the equation of motion for the second mass is: 2.1 Equation of motion 9 ¨ r=−µ r2·r r(2.10) where µ=Gm1 (standard gravitational parameter). Because m1≫m2 , it follows that m1+m2≈ m1 . For example, in the Sun-Earth system the gravitational parameter has the value µ⊙≈1.327 × 1020 m3s−2. 2.1.2 The 3D-problem For an interplanetary mission to another planet, it is possible to assume some hypotheses to simplify the problem. For instance, if the planets are orbiting in the ecliptic plane (the plane where Earth orbits around the Sun), a two-dimensional, co-planar model can be used that fits well with reality. Also, using a circular model for the orbits is a good approximation because the eccentricity of the planets is not big. However, if a mission is to be developed for other bodies that are orbiting around the Sun (like asteroids), these simplifications may not be as good as they can be in a first approach. A comparison between the eccentricity e and the inclination i for some planets and some asteroids in the Solar System is presented: Table 2.1 Comparison between eccentricities and inclinations for planets in the solar system respect to the Ecliptic plane. Planet '♀♂X Y Z [ i(º) 7.01 3.39 1.85 1.31 2.49 0.77 1.77 e 0.2056 0.0068 0.0934 0.0484 0.0541 0.0472 0.0086 Table 2.2 Comparison between eccentricities and inclinations for asteroids in the Solar System . Asteroid Didymos Apophis Bennu 3200 Phaethon Dionysus i(º) 3.41 3.33 6.035 22.2 13.9 e 0.38 0.19107 0.2037 0.89 0.542 Thus, it is evident that the approximation is not suitable for asteroids if the previously mentioned simplifications are applied, as it can be seen in Table 2.2, so inclinations and eccentricities must be included in the calculations. This project focuses on Mars and Jupiter, as well as on Dionysus and Didymos; however, the study can be extended to other bodies without difficulty. To analyse the three dimensional problem, it is necessary to define a launch date that falls within the planet’s synodic period. With this information, the orbital elements of both Earth and the target planet (or asteroid) can be predicted and the optimal transfer trajectory calculated. Since the orbital elements vary over time, they are obtained from ephemerids, which allows these results to be incorporated seamlessly into the project. 16 Chapter 2. Background in orbital mechanics 2.2.2 Thrust Constraint This section considers the practical limitations of the E-sail, focusing on the sail pitch angle αn , the characteristic acceleration, the velocities attainable by the spacecraft, and the maximum payload that the probe can carry on an interplanetary mission. First, the sail pitch angle is analysed, since it is one of the most critical parameters taken into account in the optimisation procedure (as is κ , discussed later). Current simulation results and studies place this angle in the range αnmax ∈[60◦,70◦] . Using the propulsive-acceleration vector described in (2.18), a mathematical expression can be obtained that relates αnmax to the three components of this vector at every point of the trajectory, thereby determining the feasible range of acceleration components for the mission. A typical Radial–Transverse–Normal (RTN) reference frame is introduced, where ˆ iR is the radial unit vector, ˆ iT the transverse unit vector, and ˆ iN the normal unit vector (with the origin at the spacecraft’s centre of mass). This frame differs from the body axes because the RTN frame is tied to the spacecraft’s orbit rotating slowly as position and velocity change whereas the body frame is rigidly fixed to the spinning vehicle and turns at its spin rate; in essence, RTN specifies the desired thrust direction in orbital space, while the body frame dictates how the craft must orient itself to produce that thrust. Figure 2.6 Definition of the Radial-Transverse-Normal reference frame [26]. The propulsive acceleration vector can be written as a function of the sail pitch angle αn , by projecting in the three axes that we have defined before. This is: aR=acr⊕ 2r(1+cos2αn)(2.31) aT=acr⊕ 2rsin αncos αncos σ(2.32) aN=acr⊕ 2rsin αncos αnsin σ(2.33) Here, r⊕ denotes the Sun Earth distance (1 AU) and r the heliocentric distance of the spacecraft. The equation therefore yields the components of the acceleration vector across the range of αn , but it is more conveniently expressed in the following dimensionless form: ˜aR=r⊕ 2r(1+cos2αn)(2.34) ˜aT=r⊕ 2rsin αncos αncos σ(2.35) 2.2 Continuous Thrust Model 17 ˜aN=r⊕ 2rsin αncos αnsin σ(2.36) In the first case, the unconstrained scenario is analysed, where αn spans 0◦ – 90◦ . This variation enables the description of the “force bubble”, which represents the feasible values of the accelerationvector components at a specific point along the trajectory: Figure 2.7 Representation of the acceleration constraint in terms of αn∈[−π 2,π 2]and σ∈[0,2π]. In this case r=r⊕ is adopted, but the procedure can be generalised straightforwardly. Taking into account the previously mentioned constraints on αn, the resulting force bubble is: Figure 2.8 Representation of the acceleration constraint in terms of αn∈[0,65]and σ∈[0,2π]. As stated above, the restrictions in αn affect the components of the acceleration vector, providing an easy geometrical way to verify whether the control-parameter results obtained in the optimisation 18 Chapter 2. Background in orbital mechanics Table 2.4 Spacecraft mass budget for a characteristic acceleration ac=1 mm/s2. Payload (kg) 100 200 500 1000 Number of tethers N 12 16 24 34 Tethers Length L (km) 4.02 5.77 9.27 12.9 problem are consistent with the propulsion model and discarding those that cannot satisfy these restrictions.[22] In addition, the limitations on the thrust that the E-sail can produce at a given distance from the Sun are analysed; a saturation voltage V exists, approximately defined by the condition that eV surpasses the kinetic energy of the solar-wind protons (about 1keV ), so, for a given heliocentric distance and the average number of eV that pass through the E-sail, the maximum thrust attainable at that distance can be determined.[9] The characteristic acceleration is a fundamental parameter for the E-sail; thus, an analysis is made in order to get an expression in function of the E-sail system parameters. Defining the parameter f as: f=fVV0−f0(2.37) where f0 =24.16 nNm−1 and f v =24.16nNm−1kV−1 . The spacecraft’s characteristic acceleration is therefore: ac=f N L m(2.38) In this equation, the characteristic acceleration is expressed in terms of the design parameters: N (number of tethers), L(tether length), V0(nominal voltage), and m(payload mass). This expression provides a mathematical relation between the characteristic acceleration and the payload budget, enabling the maximum payload for a feasible and realistic E-sail to be estimated for each value of ac: In this table a clear overview of the mass budget for an E-sail is obtained. It can be observed that increasing the payload requires more tethers and greater tether length. Although current E-sails are primarily based on CubeSats, where mass budget is not a significant concern, this factor will become increasingly relevant for future interplanetary and interstellar missions. Such missions will carry payloads of roughly 30–1000 kg and will demand characteristic accelerations of up to about 3ms−2. Therefore, it is important to account for this aspect in the present project.[9] 3 Optimal mission planning I n this chapter, the formulation and resolution of the optimal control problem are addressed, incorporating all the constraints that must be satisfied along the trajectory. The goal is to determine the optimal control parameters as functions of time for a continuous-propulsion model in an interplanetary trajectory, minimising the transfer time. To achieve this objective, Bezier functions are introduced as a shape-based algorithm designed for minimal computational time. Various interplanetary transfers are considered, including Earth–Mars, Earth–Jupiter, and asteroid trajectories around the Sun, such as Dionysus and Didymos. The problem is formulated as a nonlinear-programming (NLP) control problem, expressing the state variables in terms of Bézier functions that depend on the optimisation parameter tf and provide a polynomial approximation of the trajectory points. An optimisation routine is subsequently applied to refine these coefficients, yielding the minimum transfer-orbit time with the lowest computational cost. 3.1 Optimal Control Problem Formulation This section addresses the optimal control problem between two bodies, detailing the problem formulation and its solution by means of a Runge–Kutta procedure. 3.1.1 Problem hypotheses In the resolution of the OCP that has been developed in this section, the hypotheses and simplification which were decribe in the section Section 2.1.2 and 2.2.3, are taken into account. 3.1.2 The Optimal Control Problem Optimal control theory is a mathematical optimization method for finding a control law to minimize certain objective function while simultaneously being subject to a set of constraints. Given a set of nfirst-order differential equations describing generic dynamics: ˙ x=f(x,u,t)(3.1) the m control functions u (t), t ∈ [ ti,tf ] must be determined such that the following performance index J=ϕ(x(tf),tf)+Ztf ti L(x,u,t)dt (3.2) is minimized and q final boundary conditions ψ(x(tf),u(tf),tf) = 0(3.3) 19 20 Chapter 3. Optimal mission planning are satisfied. The solution to this problem is derived by the calculus of variations and its complete investigation is beyond the purposes of this report. Thus, only the derivation of the Euler-Lagrange equations will be briefly recalled. Introducing two kinds of Lagrange multipliers the q-dimensional constant vector ν for the final boundary constraints and the n-dimensional variable vector of adjoint or costate variable for the dynamics the argument performance index is defined as J=ϕ(x(tf),tf)+νTψ(x(tf),u(tf),tf)+Ztf ti [L(x,u,t)+λT((x,u,t)−˙ x)]dt (3.4) It is important to observe that the dynamics (3.1) are included in the performance index (3.2) in the same fashion as a constraint, so the optimal solution must both minimise the objective function and satisfy the dynamics. This viewpoint offers an alternative to that of dynamical system theory. Although an Euler–Lagrange approach could feasibly solve the optimisation problem, the present work adopts a direct transcription method that converts the dynamics into a nonlinear programming (NLP) problem, which is then solved with a fourth order Runge–Kutta scheme. When the cost function J consists solely of the first term ϕ , the formulation is said to be in Mayer form, as is the case here, and the problem is subject to the state equations: ˙ x(t) = f(x(t),u(t),t)(3.5) and also to a series of path constraints h(x(t),u(t),t)≤0(3.6) and boundary conditions g(x(t0),u(t0),x(tf),u(tf)) = 0(3.7) In this formulation the state vector is x= [ρ,θ,z,vρ,vθ,vz]T and the control vector is u= [κ,αn,σ]T. For a minimum time trajectory the cost function reduces to the final time itself: J(t) = tf(3.8) and the dynamics equations are: ˙ ρ=vρ(3.9) ˙ θ=vθ ρ(3.10) ˙z=vz(3.11) ˙vρ=v2 θ ρ−µ⊙ρ (ρ2+z2)3/2+κacr⊕ 2(ρ2+z2)(ρcos2αn+zsin αncos αncos σ+ρ)(3.12) ˙vθ=−vρvθ ρ+κacr⊕ 2(ρ2+z2)cos αnsin αnsin σpρ2+z2(3.13) ˙vz=−µ⊙z (ρ2+z2)3/2+κacr⊕ 2(ρ2+z2)(zcos2αn−ρsin αncos αncos σ+z)(3.14) Three different path constraints are defined: the values of the control parameters, the dynamics equations and the boundary restriction in the departure and in the arrival. These are defined as follows: 3.1 Optimal Control Problem Formulation 21 0≤κ≤1(3.15) −π 2≤αn≤π 2(3.16) 0≤σ≤2π(3.17) ρ(t0) = ρ0 θ(t0) = θ0 z(t0) = z0 vρ(t0) = vρ0 vθ(t0) = vθ0 vz(t0) = vz0 And the boundary constraint in the arrival: ρ(tf) = ρf z(tf) = zf vρ(tf) = vρf vθ(tf) = vθf vz(tf) = vzf At arrival the value of θf is not used because it is unknown. To obtain the values essential for the NLP problem, a high precision ephemeris calculator, SPICE, developed by NASA, is employed, as discussed in (2.1.1). Thus, a launch date for the E-sail must be defined because, when the elliptical and three dimensional problem is considered, rotational symmetry is lost and arbitrary parameters cannot be assigned; the position of Earth at the launch date must therefore be determined. To simplify the problem, two additional observations must be made: Astronomical Units (AU) are adopted to handle the equations conveniently, and the solar gravitational parameter is set to µ⊙=1AUUV2. Hence: 1AU =149.598 ·106km 1UV =29.7847 km/s 1UT =58.1324 days 1UA =5.9301 mm/s2 With Units of Time (UT) and Units of Acceleration (UA) Chapter 2.2.3 defines the “Force Bubble,” the geometrical zone within which the modulus of the acceleration vector must lie; on this basis, a characteristic acceleration that satisfies the governing equations can be derived. For the present project the value ac=1mms−2 is adopted, which corresponds to ac=0.16863UA. In addition, without loss of generality, the final time is assumed to satisfy tf>0 , meaning that the problem considers only the interval after the launch date and excludes any solutions that precede it. With all of these considerations in place, the final OCP formulation is stated as follows: 22 Chapter 3. Optimal mission planning min κ,αn,σtf subject to: ˙ ρ=vρ ˙ θ=vθ ρ ˙z=vz ˙vρ=v2 θ ρ−µ⊙ρ (ρ2+z2)3/2+κacr⊕ 2(ρ2+z2)(ρcos2αn+zsin αncos αncos σ+ρ) ˙vθ=−vρvθ ρ+κacr⊕ 2(ρ2+z2)cos αnsin αnsin σpρ2+z2 ˙vz=−µ⊙z (ρ2+z2)3/2+κacr⊕ 2(ρ2+z2)(zcos2αn−ρsin αncos αncos σ+z) ρ(t0) = ρ0 θ(t0) = θ0 z(t0) = z0 vρ(t0) = vρ0 vθ(t0) = vθ0 vz(t0) = vz0 ρ(tf) = ρf z(tf) = zf vρ(tf) = vρf vθ(tf) = vθf vz(tf) = vzf 0≤κ≤1 −π 2≤αn≤π 2 0≤σ≤2π t∈[0,tf],tf>0 Where ac=0.16863 UA,µ⊙=1AU UV2,r⊕=1AU 3.1 Optimal Control Problem Formulation 23 3.1.3 The Nonlinear Programming Problem Essentially, any numerical method for solving the trajectory optimisation problem incorporates an iterative procedure with a finite set of unknowns. The next section shows how an optimal-control problem can be transformed into a nonlinear programming (NLP) problem; unlike an optimal control formulation, an NLP problem involves no dynamics. Suppose that the n dimensional variable x must be chosen to solve min xF(x) subject to the m equality constraints c(x) = 0(3.18) where m ≤n. The Lagrangian of this problem is L(x,λ) = F(x)−λTc(x),(3.19) which is a scalar function of the n variables x and the m Lagrange multipliers λ . The necessary conditions for a point (x∗,λ∗)to be a constrained optimum require solving the following system ∇xL(x,λ) = g(x)−GT(x)λ=0,(3.20) ∇λL(x,λ) = −c(x) = 0(3.21) where g= ∇xF and G are the gradients of the objective function F(x) and the Jacobian of the equality constraint vector c(x), respectively. The last system of equations can be solved via Newton’s method to find the (n+m) variables ( x∗,λ∗ ) or similar numeric methods. Given a generic initial guess ( x,λ ), its corrections ( ∆x,∆λ ) to construct the new solution ( x + ∆x , λ+∆λ ) are given by solving the linear system HL-GT G0∆x ∆λ=−g −c(3.22) also referred as Karush-Kuhn-Tucker system; Karush-Kuhn-Tucker conditions the term HL is the Hessian of (3.6) in x, namely HL=∇2 xF− m ∑ i=1 λi∇2 xci(3.23) It is important to observe that an equivalent way to define the search direction ∆x is to minimize the quadratic form 1 2∆xTHL∆x+gT∆x(3.24) subject to the linear constraint G∆x=−c(3.25) This is the reason why this problem is also referred to as a quadratic programming (QP) problem. The NLP problem formulated above can be generalized to the case that occurs when inequality constraints are imposed; the mconstraints are of the form c(x)≥0(3.26) Constrains that are strictly satisfied, i.e. ci (x)> 0, are called inactive; the remaining active set of constraints are on their bounds, i.e. ci (x)=0. If the active set of constraints is known, the inactive constraints are ignored and the problem is simply solved using the method for an equality constrained problem discussed above. In summary, the general NLP problem requires finding the n vectors to 24 Chapter 3. Optimal mission planning solve min xF(x)(3.27) subject to the m constraints cL≤c(x)≤cU(3.28) and bounds xL≤x≤xU(3.29) In this formulation, equality constraints can be imposed by setting cj,L=cj,U[25] 3.1.4 Direct Transcription Optimal trajectory design is a continuous optimal control problem that can be solved with the Euler Lagrange equations; this procedure is termed the indirect method. An alternative philosophy translates the continuous optimal control problem into a nonlinear programming (NLP) problem, solving for a finite set of variables a procedure known as direct transcription, or the direct method. In a direct approach the solution of the optimal control problem is strictly connected to the numerical integration of the differential equations: the set of equations governing E-sail motion is transcribed into a finite set of equality constraints, and if the NLP solution satisfies these constraints, the original optimal control problem is solved within the numerical accuracy of the scheme employed (here, a Runge–Kutta method). A finite set of variables is first required; therefore, the problem is discretised. Defining the set of unknowns as x= [x1,x2,x3,x4,...]T , the generic NLP problem is formulated as an optimisation in which a scalar objective function F(x)is minimised. For getting these results, a set of inequality constraints and bounds are defined, represented as: subject to (cL≤c(x)≤cU xL≤x≤xU Equality path constraints and variable bounds are incorporated by imposing upper and lower limits within the optimisation problem, thus ensuring that the system dynamics satisfy the prescribed equations. Working in a discrete time domain requires the definition of a time mesh on which the dynamical equations are evaluated as path constraints, namely: t=[t0,t1,t2,...,tN−1,tf]T, where tk+1=tk+h, The parameter h , designated the step size of the time mesh, is essential because an inappropriate choice can preclude any viable solution; a value scaled to the problem, h=tf (N−1) , is therefore adopted, where N is the number of nodes. An uniform time mesh is selected at this stage[ 9 ], since no apriori information about the solution exists; this assumption is deemed reasonable for an initial approximation and can be revisited to redefine h after preliminary results have been analysed[ 25 ]. Once a finite set of time instants is defined, the unknown coefficients also become finite: each node carries six state variables plus three control parameters, giving a total of N(6+3)+1 optimisation variables, the final term representing the last component of the optimisation vector. nT= (6+3)·N+1(3.30) The “1” corresponds to the time variable, which constitutes the objective to be minimised. In the present formulation a mesh is adopted that adapts to each case; using at most 300 points yields a satisfactory compromise between precision and computational speed when solving the OCP. Once the mesh is established and the problem variables are defined, the dynamics equations are incorporated: a fourth order Runge–Kutta scheme is employed for this purpose, as discussed next. 3.1 Optimal Control Problem Formulation 25 The RK4 method is a numerical technique for solving ordinary differential equations of the form x′(t) = f(t,x(t)),x(t0) = x0. It approximates the solution x(t) by computing intermediate slopes and then combining them to update the value of x. The algorithm is described as follows: Given a step size h, the method computes: k1=f(tn,xn), k2=ftn+h 2,xn+h 2k1, k3=ftn+h 2,xn+h 2k2, k4=f(tn+h,xn+hk3). Then, the next value is given by: xn+1=xn+h 6(k1+2k2+2k3+k4). And this expression is the non-linear constraint that includes the dynamics of the problem, and must be verified by the optimal control state vector. A block diagram is made for understanding the dependence of the procedure with k. The following diagram illustrates the evaluation points in one step of the RK4 method: (tn,xn)k1=f(tn,xn) k2=ftn+h 2,xn+h 2k1 k3=ftn+h 2,xn+h 2k2 k4=f(tn+h,xn+hk3)xn+1 k1 h 2 h 2 h 2 h 2 Figure 3.1 Single-step scheme of the fourth order Runge Kutta method . It should also be noted that the values of the control parameters in the vector uk remain constant throughout each interval, and, in addition to the previously stated equations, inequality path constraints are imposed to enforce the minimum and maximum bounds on κ , αn , and σ within every interval; these constraints are expressed as: 0≤κk≤1 −π 2≤αnk ≤π 2 0≤σk≤2π With these elements in place, the full problem can now be expressed in an NLP formulation: 32 Chapter 3. Optimal mission planning sail propulsion. Unlike impulsive maneuvers, which rely on short, high energy burns, solar sails generate a constant but very low acceleration over long periods of time. In this case, with a characteristic acceleration of ac=1mm/s2 , the spacecraft gradually spirals out from Earth’s orbit to reach Mars. This results in longer transfer durations but significantly reduces the need for propellant, enabling highly efficient missions. From a physical standpoint, the time is valid and expected. The trajectory makes optimal use of the continuous solar radiation pressure, and the longer duration allows the sail to achieve the required heliocentric distance and phase angle alignment with Mars using only the gentle, cumulative push from sunlight. When focusing on the control parameters, a bang–bang behavior is observed, as anticipated; this control strategy applies maximal corrections at each point to ensure efficient spacecraft arrival at the target, in agreement with the behavior seen in the other parameters. 3.2 Simulations and results 33 3.2.2 Earth To Jupiter trajectory The transference between Earth and Jupiter is analyzed. In addition, this case is slightly more complex than before since we have to cope with bigger distances and need to be more accurate. Consequently, the simulation time has increased notably, which makes sense with the scale of the problem. However, we found a good and realistic solution. •Optimal Time ≈60.58705 UT ≈9.643 years. •Launch Date : 2031 JANUARY 01 00 : 00 : 00 UT . It is considered to be a good optimal date because it gives a reasonable time for the transference and several simulations has been done and the results show that the optimal date is close to this month. •Number of points : N=275. In this case we needed more precision than before for this reason the number of nodes have increased. •Computation time:13074.17 s≈3.632 h. Figure 3.7 Earth to Jupiter Transference with a characteristics acceleration of ac=1 mm/s2. 34 Chapter 3. Optimal mission planning Figure 3.8 Control parameters throughout the trajectory from Earth to Jupiter with ac=1 mm/s2 . Figure 3.9 Evolution of the states variables throughout the trajectory from Earth to Jupiter with ac=1 mm/s2. 3.2 Simulations and results 35 Figure 3.10 Evolution of the states velocities throughout the trajectory from Earth to Jupiter with ac=1 mm/s2. Figure 3.11 Trajectory in 2-D between Earth and Jupiter compared with the initial guess. 36 Chapter 3. Optimal mission planning In this case, it is obvious that due to the position of Jupiter (farest than Mars) it will be necessary to obtain more energy to reach the orbit, this is why we notice a change in the shape of the orbit and we see more turns around the Sun before reach to Jupiter with a rendez-vous. This behaviour, can be also seen in the position and velocity variables, in which when the spacecraft is close to the Sun the velocities reach a maximum value. Focusing on the comparison between optimal transfer and initial guess it can be said that in the first stages of the orbit either match each other (so we can assume that the initial guess is good) but in the final ones the difference between guess and final orbit is high notably, showing the feasibility of the method. 3.2 Simulations and results 37 3.2.3 Earth To Dionysus trajectory The transfer between Earth and Dionysus is analyzed using a mesh of 250 points to achieve higher precision on this more eccentric and inclined orbit. Dionysus, an Apollo type near Earth asteroid (NEA) composed primarily of silicates, offers an opportunity to gain new insights into the formation of Earth and the solar system. The results obtained are: •Optimal Time ≈21.95278 UT, this is almost 1277 days. •Launch Date : 2030 JANUARY 01 00 : 00 : 00 UT . This date was chosen after analyzing several options and determining that the optimal transfer date lay close to it. The asteroid’s synodical period is approximately 1.38years, which defines the launch window, so, with an appropriate date (such as the one selected), it is possible to reach Dionysus. •Number of points: N=250. •Computation time:12150.39 s≈3.37 h. Figure 3.12 Earth to Dionysus Transference with a characteristics acceleration of ac=1 mm/s2. 38 Chapter 3. Optimal mission planning Figure 3.13 Control parameters throughout the trajectory from Earth to Dionysus with ac= 1 mm/s2. Figure 3.14 Evolution of the states variables throughout the trajectory from Earth to Dionysus with ac=1 mm/s2. 3.2 Simulations and results 39 Figure 3.15 Evolution of the states velocities throughout the trajectory from Earth to Dionysus with ac=1 mm/s2. Figure 3.16 Trajectory in 2-D between Earth and Dionysus compared with the initial guess. 40 Chapter 3. Optimal mission planning In this section, the comparison between the final optimal trajectory and the initial guess has been added, enabling a clearer understanding of the results. Focusing on the graphical outputs, the trajectory and velocity profiles for the Dionysus mission exhibit physically coherent behavior, reflecting the gradual, continuous changes characteristic of low thrust propulsion. Rather than the abrupt jumps typical of impulsive burns, the spacecraft’s path transitions smoothly from Earth departure through transfer adjustments to final velocity matching with the asteroid. This behavior marked by modest continuous acceleration, oscillatory corrections, and progressive reduction in relative velocity aligns with optimized low thrust trajectories that respect orbital mechanics and efficient fuel management. Moreover, the highest velocities occur at the perihelion points, consistent with increased solar pressure and gravitational influence closer to the Sun. 3.2 Simulations and results 41 3.2.4 Earth To Didymos trajectory The Earth to Didymos transfer is analysed using a mesh of 275 points, offering a suitable balance between precision and computational speed. The inclusion of this small body is justified by NASA’s DART (Double Asteroid Redirection Test) mission to the Didymos-Dimorphos system, which demonstrated asteroid deflection techniques by impacting Dimorphos to alter its orbit around Didymos. The results obtained for this transfer are as follows: •Optimal Time ≈15.09409 UT, this is almost 877 days. •Launch Date : 2031 JANUARY 01 00 : 00 : 00 UT . This date is different from the previous since it is one of the possibles optimal date for the mission. •Number of points : N=275. For this case, the number of mesh points was increased to enhance solution accuracy and enable a sensitivity analysis of the transfer time; consequently, the computational time increased notably. •Computation time:20279.74s≈5.63 h. Figure 3.17 Comparison between Earth to Didymos Transference with a characteristics acceleration of ac=1 mm/s2. 48 Chapter 4. Optimal Mission Planning with Bezier Curves By increasing the number of control points, a quadratic Bezier curve can be made, which mathematical form is: B(τ) = (1−τ)[(1−τ)P0+τP1]+τ[ (1−τ)P1+τP2]where 0 ≤τ≤1(4.2) Which can be interpreted as the linear interpolant of the corresponding points on the linear Bezier curves from P0to P1and from P1to P2, this is: Figure 4.2 Bezier curve of order 2. As can be deduced, increasing the number of points improves the approximation of the original curve. When the number of control points is raised to eight, the resulting fit appears as follows: Figure 4.3 Bezier curve of order 8. The Bezier curve shapes the control points in each case; the main purpose is to determine the optimal control points that satisfy the constraints along the spacecraft’s trajectory, from which the optimal transfer path is computed. 4.1 Definition of Bezier Curves 49 4.1.1 Polynomial form The Bezier curves have been defined as a linear combination of both the parameter τ and the control points; but they can be generalized as follows: B(τ) = n ∑ i=0n i(1−τ)n−iτiPi = (1−τ)nP0+n 1(1−τ)n−1τP1+n n−1(1−τ)τn−1Pn−1+τnPn,0⩽τ⩽1 (4.3) Where n iare the binomial coefficients. When the expression is adapted to the present problem, the state variables can be expressed in terms of Bernstein polynomials and the control points.[5] For the formulation the definition of the spacecraft dimensionless cylindrical components and time are presented: [˜ ρ,θ,˜z]and τ≜t tf ,with τ∈[0,1].(4.4) Using this shape-based approach for the optimal transfer trajectory design; the components of the spacecraft (dimensionless) position vector ˜ ρ,θ,˜z are expanded in the domain of τ using Bezier functions of order n ∈N≥3, given by i(τ) = n ∑ j=0 Bj(τ)Pi,jwith i = [ ˜ ρ,θ,˜z](4.5) where Pi,j are the unknown geometric coefficients (i.e, the so-called control points), and Bj(τ) are the Bernstein basis polynomials of degree n, defined as [1]: Bj(τ) = n!τj(1−τ)n−j j!(n−j)!with j ∈[1,2,...,n](4.6) In general, a different value of n can be chosen for each spacecraft state (that is nρ,nθ,nz . However, for the sake of simplicity, in a first approach, it is considered the same number for each state variable. Taking into account equations (3.36) and (3.37), the first and the second τ− derivatives of the ith− coordinate approximation can be written, in a compact form, as i′(τ) = ˜vi= n ∑ j=0 B′j(τ)Pi,jwith i = [ ˜ ρ,θ,˜z](4.7) and i′′(τ) = ˜vi= n ∑ j=0 B′′ j(τ)Pi,jwith i = [ ˜ ρ,θ,˜z](4.8) where the derivative of the Bernstein polynomial is: B′j(τ) =            −n(1−τ)n−1if j=0 n!τj−1(1−τ)n−j (j−1)!(n−j)!−n!τj(1−τ)n−j−1 j!(n−j−1)!if j∈[1,n−1] nτn−1if j=n (4.9) and 50 Chapter 4. Optimal Mission Planning with Bezier Curves B′′ j(τ) =                                  n(n−1)(1−τ)n−2if j=0 n(n−1)(n−2)τ(1−τ)n−3−2n(n−1)(1−τ)n−2if j=1 n!τj−2(1−τ)n−j (j−2)!(n−j)!−2n!τj−1(1−τ)n−j−1 (j−1)!(n−j−1)! +n!τj(1−τ)n−j−2 j!(n−j−2)!if j∈[2,n−2] n(n−1)(n−2)τn−3(1−τ)−2n(n−1)τn−2if j=n−1 n(n−1)τn−2if j=n (4.10) In particular, the boundary values of Bj and B′j are obtained by substituting τ=0 (initial time instant), or τ=1(final time) in equations (3.40) and (3.41), we get: Bj(0) =    1if j=0 0if j∈[1,n] (4.11) Bj(1) =    0if j∈[0,n−1] 1if j=n (4.12) B′j(0) =            −nif j=0 nif j=1 0if j∈[2,n] (4.13) B′j(1) =            0if j∈[0,n−2] −nif j=n−1 nif j=n (4.14) Using equations (3.36) and (3.38), the boundary value of the generic dimensionless ith− coordinate approximation and its first τ−derivative are given by i(0) = Pi,0,i(1) = Pi,n,˜vi(0) = nPi,1−Pi,0, ˜vi(1) = nPi,n−Pi,n−1(4.15) with i= [˜ ρ,θ,˜z] . The twelve geometric coefficients [Pi,0,Pi,1,Pi,n,Pi,n−1] can be obtained, as a function of the given components of the spacecraft state vector at the initial and the final time tf , by combining P˜ ρ,0=ρ0 r⊕ ,Pθ,0=θ0,P˜z,0=z0 r⊕ (4.16) 4.2 Problem formulation including Bezier Curves 51 P˜ ρ,1=vρ0tf r⊕ +ρ0 r⊕ ,Pθ,1=vθ0tf nr⊕ +θ0,P˜z,1=vz0tf r⊕ +z0 r⊕ (4.17) P˜ ρ,n−1=ρf r⊕−vρftf nr⊕ ,Pθ,n−1=θf−vθftf nr⊕ ,P˜z,n−1=zf r⊕−vz f tf nr⊕ ,(4.18) P˜ ρ,n=ρf r⊕ ,Pθ,n=θf,P˜z,n=zf r⊕ .(4.19) As described above, the Bezier curves method has been defined for the present problem and can now be formulated to obtain the optimal transfer; since the arrival and departure conditions are fixed by the Bezier coefficients, only the free coefficients along the trajectory require optimisation, resulting in 3(n−3)+1unknown parameters.[6] 4.2 Problem formulation including Bezier Curves For a given mission scenario, the problem is to find the minimum time transfer trajectory such that the constraints of the acceleration vector are all satisfied. To that end, the components of the propulsive acceleration vector {ax0,ay0,az0} are expressed as functions of the Bezier curve-based approximation given by equations (4.4) and (4.3) using the following procedure. For a given flight time tf and a set of N unknown geometric coefficients {Pi,2,...,Pi,n−2} , with n>3and i∈{˜ ρ,˜ θ,˜z} the dimensionless components of the propulsive acceleration {˜aρ,˜aθ,˜az} can be obtained from the equations of motion by moving the gravitational and inertial terms from the right-hand side to the left-hand side. The components of the propulsive acceleration {aρ,aθ,az} are then calculated from equations (2.23), whereas the components {ax0,ay0,az0} may be written, according to equations (2.14), (2.22) and (4.26), as   ax0 ay0 az0  =  cosφ0−sinφ 0 1 0 sinφ0 cosφ   aρ aθ az  (4.20)   ax0 ay0 az0  =         aρ˜z−az˜ ρ p˜ ρ2+˜z2 aθ aρ˜ ρ+az˜z p˜ ρ2+˜z2         ,(58) where φis the orbital inclination angle. In order to obtain the components of the dimensionless acceleration vector, it can be written the dynamical equations in a dimensionless form using r⊕ and the tf that define the parameter τ . Doing this, we obtain the following expressions: τ≜t tf ,τ∈[0,1](4.21) ˜ ρ′=˜v˜ ρ≜f˜ ρ(4.22) 52 Chapter 4. Optimal Mission Planning with Bezier Curves θ′=˜vθ ˜ ρ≜fθ(4.23) ˜z′=˜v˜z≜f˜z(4.24) ˜v′˜ ρ=˜v2 θ ˜ ρ−˜ µ⊙˜ ρ (˜ ρ2+˜z2)3/2+κ˜ac 2˜ ρcos2αn+˜zsinαncosαncosσ+˜ ρ≜f˜v˜ ρ(4.25) ˜v′ θ=−˜v˜ ρ˜vθ ˜ ρ+κ˜ac 2cosαnsinαnsinσp˜ ρ2+˜z2≜f˜vθ(4.26) ˜v′˜z=−˜ µ⊙˜z (˜ ρ2+˜z2)3/2+κ˜ac 2˜zcos2αn−˜ ρsinαncosαncosσ+˜z≜f˜v˜z(4.27) Where the variables are defined by: ˜ ρ≜ρ r⊕ ,˜z≜z r⊕ ,˜ µ≜µ⊙t2 f r3 ⊕ ,˜ac≜act2 f r⊕ , ˜v˜ ρ≜vρtf r⊕ ,˜vθ≜vθtf r⊕ ,˜v˜z≜vztf r⊕ . (4.28) With this formulation, the dimensionless components of the spacecraft’s acceleration vector are computed directly from the Bezier coefficients; by enforcing the acceleration constraints at the acceleration level since acceleration is expressed in the same polynomial basis the need to integrate the dynamics to recover velocities and positions is obviated, and the acceleration vector constraints are imposed continuously along the entire trajectory. Focusing on this expressions we can easily identify the components that we are looking for: ˜a˜ ρ=κ˜ac 2(˜ ρcos2αn+˜zsinαncosαncosσ+˜ ρ)(4.29) ˜a˜ θ=κ˜ac 2cosαnsinαnsinσp˜ ρ2+˜z2(4.30) ˜a˜z=k˜ac 2˜zcos2αn−˜ ρsinαncosαncosσ+˜z(4.31) By substituting these expressions into equations (4.24)–(4.26) and transferring the gravitational term to the left hand side, the dimensionless acceleration components can be written explicitly as functions of the Bezier coefficients ˜v′ θ+˜v˜ ρ˜vθ ˜ ρ=˜a˜ ρ(4.32) ˜v′ θ+˜v˜ ρ˜vθ ˜ ρ=˜a˜ θ(4.33) ˜v′˜z+˜ µ⊙˜z (˜ ρ2+˜z2)3/2=˜a˜z(4.34) Since the dimensional components have been expressed in terms of the dimensionless ones, the analytical relations between them follow by using the normalized time parameter τ=t/tf , yielding the equivalences: 4.2 Problem formulation including Bezier Curves 53 d2 dt2=1 t2 f d2 dτ2(4.35) aτ(τ) = t2 f r⊕ a(t)⇐⇒ a(t) = r⊕ t2 f aτ(τ)(4.36) So, it can be obtained the relationship between dimensional and dimensionless components: aρ=r⊕ t2 f ˜aρ(4.37) aθ=r⊕ t2 f ˜aθ(4.38) az=r⊕ t2 f ˜az(4.39) With these results, the dimensional acceleration vector components are obtained in terms of the Bezier coefficients noting that the dimensionless components are computed via the Bezier expansion and the cartesian frame components [ax0,ay0,az0]can then be calculated as follows:   ax0 ay0 az0  =  cosφ0−sinφ 0 1 0 sinφ0 cosφ   aρ aθ az  =         aρ˜z−az˜ ρ p˜ ρ2+˜z2 aθ aρ˜ ρ+az˜z p˜ ρ2+˜z2         .(4.40) With these expressions, the Cartesian components of the acceleration are defined by the Bézier coefficients, enabling us to formulate the acceleration constraint that the spacecraft must satisfy continuously along its trajectory. 4.2.1 Acceleration Constraints In Chapter 2, the constraint was formulated directly from the system dynamics. Now that the differential equations have been removed from the constraint set, those constraints must be rephrased in terms of the acceleration vector, which itself is a function of the Bezier coefficients. In so doing, the trajectory restrictions are imposed directly on the coefficients being optimized. To that end, the transverse unit vector is introduced: ˆ t≜ˆ r׈ n׈ r sinαn with αn=0(4.41) The propulsive acceleration vector lies on the plane spanned by ˆr and ˆ t since it may be written in terms of radial and transverse components: a=acr⊕ r˜arˆ r+˜atˆ t,(4.42) where ˜ar=κcos2αn+1 2,˜at=κsinαncosαn 2.(4.43) 54 Chapter 4. Optimal Mission Planning with Bezier Curves Making a comparison between equations (4.37) and (2.20) reveals that: ar≜˜aracr⊕ r≡az0(4.44) at≜˜atacr⊕ r=qa2 x0+a2 y0(4.45) And combining equations (4.37) it may be verified that: ˜a2 t+˜ar−3 4κ2=κ 42(4.46) As κ varies over the interval [0,1] , the circles defined by equation (4.40) enclose a region S in which the dimensionless thrust acceleration must lie. Because the sail pitch angle αn∈[−π 2,π 2] (so that the normal ˆn always points anti-Sun), the tangent from the origin to any of these circles satisfies ˜ar=2√2 ˜at, regardless of κ . Therefore, the lower boundary of S is precisely this line of slope 2√2 , and one may characterize S=(˜at,˜ar)˜ar⩾2√2 ˜at,(˜ar,˜at)lies on or below the circle for some κ∈[0,1]. Figure 4.4 Dimensionless Propulsive acceleration components and admissible region for the propulsive acceleration. With the aid of the geometry in the figure, the region can be divided into two constraints in terms of ˜ar; thus, the dimensionless constraints to be satisfied by the acceleration vector are: S=˜ar,˜at: ˜at⩽       ˜ar 2√2,˜ar∈[0,2 3], r1 16 −˜ar−3 42,˜ar∈[2 3,1]. (4.47) 4.2 Problem formulation including Bezier Curves 55 Finally, for a given radius r and characteristic acceleration ac , the propulsive acceleration components {ax0,ay0,az0}must satisfy: qa2 x0+a2 y0⩽az0 2√2,az0 ac(r⊕/r)∈[0,2 3],(4.48) qa2 x0+a2 y0⩽ac r⊕ rs1 16 −az0 ac(r⊕/r)−3 42,az0 ac(r⊕/r)∈[2 3,1].(4.49) Finally, the restrictions on the acceleration vector components that must be satisfied by the spacecraft at each point have been obtained; to enforce these constraints, m=150 points are chosen along the trajectory.[6] 4.2.2 Fixed Coefficient for the transference In the optimization, we employ MATLAB’s fmincon to determine the free Bezier control points namely the 3(n−3) unknown parameters while the twelve boundary coefficients remain fixed by the departure and arrival conditions. By combining both fixed and optimized coefficients, we reconstruct the complete interplanetary trajectory. The values of the fixed coefficients were defined previously in equations (4.16)-(4.19) These coefficients are entirely dictated by the departure and arrival boundary conditions. To streamline their use in our optimization procedure, we assemble them into a single vector with the following structure: P=Pρ,0Pρ,1··· Pρ,n−1Pρ,nPθ,0Pθ,1··· Pθ,n−1Pθ,nPz,0Pz,1··· Pz,n−1Pz,nT(4.50) By structuring the fixed boundary coefficients this way, we avoid recomputing the departure and arrival conditions at each iteration. Instead, we impose them once as constraints and let the optimizer solve only for the remaining free parameters, significantly improving computational efficiency. 56 Chapter 4. Optimal Mission Planning with Bezier Curves 4.2.3 Initial Guess The basic idea to initialize the unknown coefficients in the NLP solver is to provided an approximation of the spacecraft coordinates [ ρ,θ,z ] at the m discretization points. The unknown coefficients can them be calculated by fitting the functions of Bezier curves throughout this set of discrete points. Note that the total flight time tf is a variable that must be optimized. However, a guess value of the flight time tfapp may be obtained as the ratio of the norm of the vectorial difference between the angular momentum of target and parking orbits to the E-sail torque induced by thrust. This expression is: tfapp =qµ⊙p0+pf−2µ⊙√p0pfcos∆i κapp acr⊕cosαn1sinαn1 (4.51) where p0/pf are the semilatus rectum of the parking orbit and the target orbit respectively, also, ∆i is the variation in orbital inclination, κapp ≜1 3 is the estimated average thrust coefficients, and αn1 ≜ 54.7 deg is the sail pitch angle that maximizes the thrust cone angle. To get a better impression about this approximation, in the case of Mars, the following results are given: p0=1AU pf=1.5104 AU ac=0.1686 UA tfapp =8.7427 UT (4.52) Comparing this first approximation with the optimal result obtained for the Earth–Mars direct trajectory, topt =8.191067UT , demonstrates that it constitutes a good starting point for the optimisation process. First, after an initial estimate of the transfer time has been computed, a preliminary guess for the free Bezier coefficients is generated by fitting a third-order Bezier curve to the boundary conditions thus determining all coefficients and evaluating this curve at the m discretization points; these values serve as initial estimates for the free coefficients, which are then supplied to the optimizer for iterative refinement. It is necessary to calculate the fixed coefficients: Pρ,0=ρi,Pρ,1=ρi+tfapp vρi 3,Pρ,2=ρf−tfapp vρf 3,Pρ,3=ρf, Pθ,0=θi,Pθ,1=θi+tfapp vθi 3,Pθ,2=θf−tfapp vθf 3,Pθ,3=θf, Pz,0=zi,Pz,1=zi+tfapp vzi 3,Pz,2=zf−tfapp vzi 3,Pz,3=zf. (4.53) Note that the value of θf needs to be updated in each τk step, due to this, value changes through the trajectory. This is: θfk+1=θfk+2nπn=1,2... (4.54) This is specially essential with the initial guess generation, since if this detail is not incorporated, a wrong results could be obtained. With the third order Bezier curve, the approximation of ρapp,θapp,zapp can be written as: ρapp(τ) = (1−τ)3Pρ,0+3τ(1−τ)2Pρ,1+3τ2(1−τ)Pρ,2+τ3Pρ,3, θapp(τ) = (1−τ)3Pθ,0+3τ(1−τ)2Pθ,1+3τ2(1−τ)Pθ,2+τ3Pθ,3, zapp(τ) = (1−τ)3Pz,0+3τ(1−τ)2Pz,1+3τ2(1−τ)Pz,2+τ3Pz,3. (4.55) 4.2 Problem formulation including Bezier Curves 57 Consequently, the discrete approximation data values of ρapp,θapp and zapp can be obtained by evaluating equations at the m discretization points, for this work m=150. Then an initial guess for the unknown parameters in bezier functions can be obtained from equation (4.4). Besides, a flow diagram is presenter in order to get a better comprehension about how to implement these last sections together. Initial Guess Process Input data x0= (ρi,˙ ρi,θi,˙ θi,zi,˙zi), xf= (ρf,˙ ρf,θf,˙ θf,zf,˙zf), ac,r⊕,κapp,n 1. Estimate tfapp tfapp =√µ⊕(p0+pf)−2µ⊕√p0pfcos∆i κapp acr⊕cosαn1sinαn1 2. Fixed coefficients Cubic Bézier ( n=3 ). For each i∈{ρ,θ,z}: Pi,0=i0 Pi,1=i0+tfapp ˙ i0 3 Pi,2=if−tfapp ˙ if 3 Pi,3=if 3. Evaluate approximation Compute {ρapp(τk),θapp(τk),zapp(τk)} at k=1,...,m discretization nodes. 4. Initial Solve Solve for the 3(n−3) free coefficients Pi,2...n−2by imposing n ∑ j=0 Bj(τk)Pi,j= iapp(τk) (k=1,...,m). 5. Assemble vector {tf,app,Pi,0...,n} ready for fmincon. End of Guessing Process The figure below gathers the complete set of initial guesses that will be used as starting points for each celestial body considered in this work. Every entry in this table has been produced with the procedure detailed earlier: the cubic Bezier framework is fixed, boundary coefficients are assigned, and the remaining degrees of freedom are solved by enforcing the approximation conditions, presenting the guesses side-by-side serves a dual purpose. First, it highlights the internal consistency of the method, showing how a single algorithm yields coherent results across a wide spectrum of initial conditions. Second, it offers a qualitative benchmark before the vectors are handed to the fmincon optimizer for further refinement. All symbols, units, and sign conventions match those introduced in the preceding sections, ensuring traceability from the theoretical development to the graphical validation. 64 Chapter 4. Optimal Mission Planning with Bezier Curves 4.3.2 Earth to Jupiter Trajectory with Bezier curves The solution for Earth to Jupiter is analyzed. It is easy to appreciate that the computation time has decreased and, although, the tolerance for constraints has been increased to address the acceleration vector constraints the results are highly realistic. The step and constraint tolerances that are being used have a value of TOL 10−3. •Optimal Time ≈55.02569 UT which is ≈8.76795 years. •Launch Date : 2031 JANUARY 01 00 : 00 : 00 UT . Close to the launch window for Jupiter in 2031 •Bezier curve order: n=8. •Computation time:74.8118s Figure 4.14 Earth to Jupiter Transference using Bezier Curves. 4.3 Simulations and results using Bezier curves 65 Figure 4.15 Control parameters throughout the trajectory from Earth to Jupiter with a Bezier curve of order 8. Figure 4.16 Evolution of the states variables throughout the trajectory from Earth to Jupiter with a Bezier curve of order 8. 66 Chapter 4. Optimal Mission Planning with Bezier Curves Figure 4.17 Evolution of the states velocities throughout the trajectory from Earth to Jupiter with a Bezier curve of order 8. Figure 4.18 Trajectory in 2-D between Earth and Jupiter compared with the initial guess. 4.3 Simulations and results using Bezier curves 67 To obtain a physically feasible solution and for the convergence of the problem, the tolerance was increased. Although this lower tolerance may make the results less reliable, and they can be unexpected from the solution that was though, they still produced a physically plausible trajectory whose final transfer time matches that of the NLP solution. The future work must focus on enhancing this. Note that for this solution, the final optimal curve has a broaden form than in the NPL problem, but a continuous evolution of either, state variables and velocities shows that this trajectory is physically possible. Moreover, the key factor for introducing this method (as it has been mentioned before) is the elimination of rough changes in the acceleration vector, obtaining a softer shape for the thrust parameter κ, which avoid numerical problem, specially if derivation is necessary. Moving on the bidimensional representation, the initial guess curve differs a little bit more than in the previous one, this might be a signal for improving the accuracy of the initial guess method, when further distances are involved. It is likely that thus it has been necessary to decrease the tolerance constraints. The loops that has been added are slightly superior (n=3) this due to the position of Jupiter, about 5.2 AU, so the spacecraft will need more to reach the objective.[4.2.3] 68 Chapter 4. Optimal Mission Planning with Bezier Curves 4.3.3 Earth to Dionysus Trajectory with Bezier curves In this section the transference between Earth and Dyonisus is presented, in order to get these results the constraint tolerance has been increased (TOL 10−4 ), since meeting the acceleration vector constraint has been more difficult this time, however the results that are obtained are physically possibles. In addition, a characteristic acceleration of ac=1 mm/s2 is used. As it has been saw before, the reduction on the computational time is quite significant, so we can prove that the method works. The results are: •Optimal Time ≈20.06444 UT which is 1167.75 days. •Launch Date:2031 JANUARY 01 00 : 00 : 00 UT. •Bezier curve order: n=8. •Computation time:202.4325s. Figure 4.19 Earth to Dionysus Transference using Bezier Curves. 4.3 Simulations and results using Bezier curves 69 Figure 4.20 Control parameters throughout the trajectory from Earth to Dionysus with a Bezier curve of order 8. Figure 4.21 Evolution of the states variables throughout the trajectory from Earth to Dionysus with a Bezier curve of order 8. 70 Chapter 4. Optimal Mission Planning with Bezier Curves Figure 4.22 Evolution of the states velocities throughout the trajectory from Earth to Dionysus with a Bezier curve of order 8. Figure 4.23 Trajectory in 2-D between Earth and Dionysus compared with the initial guess. 4.3 Simulations and results using Bezier curves 71 The velocity profiles remain free of spikes, ensuring that the solar-sail endures only gentle, continuous thrusting; this is reflected in the control-parameter plots, where the “pitch” angle αn switches only once from a high-thrust regime near departure to a shallow, coast-optimized attitude, while the clocking angle σ and curvature κ describe single, well shaped lobes rather than oscillatory bang-bang behavior. In the 2-D projection a better estimation of the initial guess might be done in or to put forward a more efficient optimization process, moreover, a couple of manual loops have been added for getting a good initial guess.[4.2.3] By comparison, the classical direct-transcription approach (N = 250 RK4 nodes) reaches Dyonisus in about 1 277 days nearly 110 days longer while demanding over 12150s(≈3.4h) of CPU time. Thus, the Bezier method not only trims flight time by 9 % but also slashes computation time from hours to just ≈200s , all without compromising the physical realism of distances, velocities, or control demands. 72 Chapter 4. Optimal Mission Planning with Bezier Curves 4.3.4 Earth to Didymos Trajectory with Bezier curves Finally the rendezvous trajectory between Earth and Didymos is presented. This time a slightly inferior optimal time has been obtained but despite that, the result is completely valid in a early stage design mission. However in this case, a solution with a high value of tolerance (TOL 10−6 ) has been found, as an improvement from the previous analysis. The results are: •Optimal Time ≈11.65088 UT which is 678.03 days. •Launch Date:2031 JANUARY 01 00 : 00 : 00 UT. •Bezier curve order: n=8. •Computation time:259.07384s Figure 4.24 Earth to Didymos Transference using Bezier Curves. 4.3 Simulations and results using Bezier curves 73 Figure 4.25 Control parameters throughout the trajectory from Earth to Didymos with a Bezier curve of order 8. Figure 4.26 Evolution of the states variables throughout the trajectory from Earth to Didymos with a Bezier curve of order 8. 80 Chapter 5. Conclusions and future work versus 877 days and 5.6 h for the conventional solver which shows the high effectivity of the polynomial method. Moreover,a dedicated study on Dionysus showed that raising the Bezier degree from n=10 to n=40 shortens tf by 17 % but increases runtime nine fold; diminishing returns appear beyond n≈20 , which captures 97 % of the attainable saving for only 55 % of the maximum CPU cost. This identifies n∈[10,20]as a sweet spot for preliminary E-sail design. Table 4.2 demonstrated that propulsive upgrades have their greatest leverage in the low thrust regime: raising ac from 0.5 to 0.6 mm/s2 shaved almost 400 days off the Earth to Mars cruise, whereas the same 0.1 mm/s2 increment above 1 mm/s2 yielded a mere 13 day improvement. With this approach also a optimal range of the achas been analyzed. The smoothness and explicit accelerations continuity of Bezier trajectories make them ideal reference profiles for the LQR and Model Predictive Control schemes that can be added in future projects. Because the polynomial solution delivers analytic state derivatives, the guidance laws can be tuned without numerical differentiation, reducing onboard computational load and improving robustness to solar wind fluctuations. All trajectories were cross checked against NASA HORIZONS ephemerids; graphical overlays confirmed that the optimised paths remain within a few hundredths of an AU of the published orbital planes. Benchmark comparisons (Appendix A) further showed excellent agreement with flight times reported in recent E-sail studies. In summary, this thesis establishes a reproducible, rapid convergence pipeline for E-sail mission early design due to: • Reduces optimisation runtimes from hours to seconds without sacrificing too much physical accuracy. •Provides smooth, controller-friendly state and control histories. •Quantifies key design trades involving characteristic acceleration and polynomial order. •Demonstrates versatility across planetary and high-eccentricity asteroid targets. Taken together, these achievements mark a significant step towards integrating advanced shape based techniques into routine interplanetary mission planning, paving the way for cost effective exploration of the Solar System with propellantless continuous thrust. All of these studies have heightened our awareness of the numerous optimal trajectories within the Solar System and, moreover, have underscored how crucial it is to account for mission time. 5.2 Future Work This thesis has deliberately relied on a set of simplifying assumptions to render the E-sail optimisation problem tractable. At the same time, Optimal Control Theory is evolving rapidly and is being extended to ever more sophisticated aerospace systems. Consequently, a wide spectrum of enhancements can be envisioned for the methodology developed here. The results presented offer a clear baseline of what can be achieved with current shape based and direct transcription techniques for Electric Solar Wind Sail mission design; the prospective refinements outlined below should therefore be viewed either as means to sharpen the present solutions or as adaptations that enable the framework to tackle new targets and operational scenarios. The most promising avenues for improvement are summarised in the following list. • Improve ConvergenceThe future work should be to guarantee convergence in every tested scenario without decreasing the tolerance with sufficient accuracy. Achieving this goal entails outlining concrete research directions for instance, re examining the initial iterate to improve 5.2 Future Work 81 its quality, experimenting with alternative optimisers, such as CasaDi, and assessing their impact on both solution precision and computation time. These steps will provide a robust path toward reliably meeting the desired convergence criteria in future studies. • Adaptive mesh refinement. Replace the current uniform grid with an adaptive time–mesh that allocates extra nodes only where the dynamics are stiff or switching occurs, preserving accuracy while reducing the NLP size and CPU time. • Advanced model predictive guidance. Integrate receding and shrinking horizon MPC variants complete with horizon adaptation and constraint tightening to forecast solar-wind fluctuations and generate control actions that keep the E-sail within its thrust “force bubble” even under uncertainty. • Concurrent launch window optimisation. Treat the departure date as an optimisation variable, allowing the solver to pick both launch epoch and thrust history that jointly minimise total time of flight, rather than fixing the launch date. • Meta heuristic initial guess generation. Use population based algorithms genetic, differential evolution or particle swarm to deliver diverse, high quality starting points for the NLP, enlarging the basin of convergence and revealing alternative trajectory families. • Gravity assist trajectories. Extend the dynamics to include third body perturbations and design continuous thrust transfers that exploit planetary flybys, aiming for shorter cruise times and lower characteristic acceleration requirements. • Robust control against solar wind variability. Replace deterministic thrust bounds with a stochastic model of dynamic pressure and synthesise robust LQR/MPC controllers that meet terminal constraints with a chosen confidence level. • Adaptive Bezier segmentation. Let the Bezier degree and the number of patches vary along the path, capturing sharp manoeuvres without global over parameterisation and keeping the decision-vector compact. • Indirect method benchmarking. Implement a Pontryagin based shooting solver to cross validate direct-transcription and Bezier results, helping detect possible sub-optimal local minima and quantify optimality gaps. • Onboard reoptimisation capability. Exploit the low-dimensional Bezier parameter set to enable real time, mid course retargeting whenever ephemerids are updated or sail degradation estimates change. Appendix A Code Validation In this section, the results obtained with the MATLAB code developed herein are compared with the solutions commonly reported in the literature; this comparison validates the code for the present scenario. A.1 Graphical Validation Firstly, the plot produced by the developed code for the planetary and asteroid orbits is compared with the actual orbit obtained from the NASA Horizons System. A.1.1 Earth to Mars Figure A.1 Comparison between orbital planes for Earth and Mars with NASA Horizons System . Figure A.2 Comparison between orbital planes for Earth and Mars in the NLP problem . 83 84 Chapter A. Code Validation For this case, it is known that Earth and Mars orbit in the ecliptic plane, which has an inclination of about i=23.5◦ . In the first image (A.1), both planets are seen in the ecliptic plane as a reference, so it can be confirmed that they orbit within the same plane. In the second image, it can be observed that, unless they are rotated, both bodies orbit together, thereby verifying the graphical solution. A.1.2 Earth to Jupiter Figure A.3 Comparison between orbital planes for Earth and Jupiter with NASA Horizons System. Figure A.4 Comparison between orbital planes for Earth and Jupiter in the NLP problem. Both orbits closely resemble that depicted in the top image; consequently, the scenario under consideration accurately represents the formulated problem. As in the previous case, both planetary orbits lie within the ecliptic plane, so no significant inconvenience is expected for a spacecraft to perform a rendezvous with Jupiter. A.1 Graphical Validation 85 A.1.3 Earth to Dionysus Figure A.5 Comparison between orbital planes for Earth and Dionysus with NASA Horizons System. Figure A.6 Comparison between orbital planes for Earth and Mars in the NLP problem. As in the previous case, the orbits appear to be quite similar, and the high eccentricity of Dionysus is reproduced with accuracy. The relative inclination between Earth and Dionysus is likewise consistent in both cases. Additionally, a point at which the orbits approach each other very closely can be observed, and this feature is present in the representation. Therefore, it can be concluded that the model provides a good approximation of the real orbital configuration. 86 Chapter A. Code Validation A.1.4 Earth to Didymos Figure A.7 Comparison between orbital planes for Earth and Didymos with NASA Horizons System. Figure A.8 Comparison between orbital planes for Earth and Didymos in the NLP problem. Both figures show the orbital planes of Earth and Didymos around the Sun. The model (Figure A.8) closely matches the NASA Horizons reference (Figure A.7), accurately capturing the inclination between the orbits and the elliptical shape of Didymos’ trajectory. The spatial proportions are consistent. Despite stylistic differences, the physical representation in our model remains faithful to the real orbital mechanics. This confirms that trajectory data and orbital elements are well implemented and realistic in this new case. A.2 Numerical Validation 87 A.2 Numerical Validation To test the validity of the numerical results, the transfer orbit is computed for the cases described in Chapter 3 and compared with other three dimensional models, providing an overall validation of the obtained results. A characteristic acceleration of 1mm/s2is employed. A.2.1 Earth to Mars In this section, two models are compared to demonstrate the accuracy of the developed code. Although slight differences exist between the models, these are explained herein. Another model, namely [ 18 ], is employed for comparison. The corresponding values are presented in the next table. For consistency, a departure date of 1 February 2029 is selected, acknowledging that this date is not optimal. Table A.1 Flight Times comparative between different models for a Earth to Mars trajectory with ac=1 mm/s2with a departure date on 1st February 2029. Study This Project 3 Dimensional model from [18] Flight Time [Days] 650 637 Both results are quite similar, differing by only 13 days; therefore, the code can be considered valid for the analysis, and confidence may be placed in the resolution method, as it yields realistic results. A.2.2 Earth to Dionysus As in the previous verification, the model is compared for the Earth to Dionysus transfer; a launch date of 20 March 2024 is selected for the analysis, and the following results are obtained. Table A.2 Flight Times comparative between different models for a Earth to Dionysus trajectory with ac=1 mm/s2with a departure date on 20th March 2024. Study This Project 3 Dimensional model from [18] Flight Time [Days] 1362 1078 In this case, results differ from those obtained in the previous example because a reduced number of points was employed. Comparison with the analysis presented in 3.6, where a final time of approximately tf=1200days was obtained, shows a closer agreement with the reference model. Nevertheless, the results are physically plausible, so confidence can be placed in the present analysis. List of Figures 1.1 Vogayer 1 Probe Propulsion System, showing RTG [15] 2 1.2 Electric Solar Wind Sail concept [18] 3 1.3 Evolution of the E-sails throughout the time [10] 4 2.1 Description of the gravitational Force between two bodies 8 2.2 Description of the cylindrical reference frame τc11 2.3 Definition of the position and normal vectors of the E-sail for the propulsive acceleration vector 13 2.4 Definition of the reference frame τ014 2.5 Representation of the ratio a/ac for different values of the control parameter κ14 2.6 Definition of the Radial-Transverse-Normal reference frame [26] 16 2.7 Representation of the acceleration constraint in terms of αn∈[−π 2,π 2]and σ∈[0,2π]17 2.8 Representation of the acceleration constraint in terms of αn∈[0,65]and σ∈[0,2π]17 3.1 Single-step scheme of the fourth order Runge Kutta method 25 3.2 Earth to Mars Transference with a characteristics acceleration of ac=1 mm/s229 3.3 Control parameters throughout the trajectory from Earth to Mars with ac=1 mm/s230 3.4 Evolution of the states variables throughout the trajectory from Earth to Mars with ac=1 mm/s230 3.5 Evolution of the states velocities throughout the trajectory from Earth to Mars with ac=1 mm/s231 3.6 2-D transference in comparison to the initial guess 31 3.7 Earth to Jupiter Transference with a characteristics acceleration of ac=1 mm/s233 3.8 Control parameters throughout the trajectory from Earth to Jupiter with ac=1 mm/s234 3.9 Evolution of the states variables throughout the trajectory from Earth to Jupiter with ac=1 mm/s234 3.10 Evolution of the states velocities throughout the trajectory from Earth to Jupiter with ac=1 mm/s235 3.11 Trajectory in 2-D between Earth and Jupiter compared with the initial guess 35 3.12 Earth to Dionysus Transference with a characteristics acceleration of ac=1 mm/s237 3.13 Control parameters throughout the trajectory from Earth to Dionysus with ac=1 mm/s238 3.14 Evolution of the states variables throughout the trajectory from Earth to Dionysus with ac=1 mm/s238 3.15 Evolution of the states velocities throughout the trajectory from Earth to Dionysus with ac=1 mm/s239 3.16 Trajectory in 2-D between Earth and Dionysus compared with the initial guess 39 3.17 Comparison between Earth to Didymos Transference with a characteristics acceleration of ac=1 mm/s241 3.18 Control parameters throughout the trajectory from Earth to Didymos with ac=1 mm/s242 89 96 Bibliography [14] Colin R. Mcinnes and Matthew P. Cartmell, 7 - orbital mechanics of propellantless propulsion systems, Modern Astrodynamics (Pini Gurfil, ed.), Elsevier Astrodynamics Series, vol. 1, Butterworth-Heinemann, 2006, pp. 189–235. [15] Giovanni Mengali, Alessandro A. Quarta, and Pekka Janhunen, Electric sail performance analysis, Journal of Spacecraft and Rockets 45 (2008), no. 1, 122–129. [16] Sini Merikallio and Pekka Janhunen, The electric solar wind sail (e-sail): Propulsion innovation for solar system travel, (2018). [17] A. Miele, M. Ciarcià, and J. Mathwig, Reflections on the hohmann transfer, Journal of Optimization Theory and Applications 123 (2004), 233–253, Includes numerical studies of Earth–Mars and Earth–Jupiter Hohmann transfers. [18] Alessandro A.Quarta Mingying Huo, Giovanni Mengali, Electric sail trhust model from a geometrical perspective, (2018). [19] Michael A. Minovitch, The invention that opened the solar system to exploration, Planetary and Space Science 58 (2010), no. 6, 885–892. [20] NASA Jet Propulsion Laboratory – NAIF, Spice toolkit,https:// naif.jpl.nasa.gov/ naif/ toolkit. html, 2025, Accessed 28 Jun 2025. [21] Alessandro Peloni, Anil V. Rao, and Matteo Ceriotti, Automated trajectory optimizer for solar sailing (atoss), Aerospace Science and Technology 72 (2018), 465–475. [22] Alessandro Quarta, Impact of pitch angle limitation on e-sail interplanetary transfers, Aerospace 11 (2024), 729. [23] Alessandro A. Quarta, Impact of pitch angle limitation on e-sail interplanetary transfers, Aerospace 11 (2024), no. 9. [24] Jinjun Shan and Yuan Ren, Low-thrust trajectory design with constrained particle swarm optimization, Aerospace Science and Technology 36 (2014), 114–124. [25] Francesco Topputo and Chen Zhang, Survey of direct transcription for low-thrust space trajectory optimization with applications, Abstract and Applied Analysis 2014 (2014), 1–15. [26] Fan Zichen, Mingying Huo, Naiming Qi, Ce Zhao, Ze Yu, and Tong Lin, Initial design of low-thrust trajectories based on the bezier curve-based shaping approach, Proceedings of the Institution of Mechanical Engineers, Part G: Journal of Aerospace Engineering 234 (2020), 095441002092004.