Single-stage to orbit ascent trajectory optimisation with reliable evolutionary initial guess

In this paper, the ascent trajectory optimization of a lifting body Single-Stage To Orbit (SSTO) reusable launch vehicle is investigated. The work is carried out using a Direct Multiple Shooting method to solve the Optimal Control problem. The crucial initialisation of the optimisation process is performed by using a combination of two evolutionary algorithms, namely a Multi-Objective Parzen-based Estimation of Distribution (MOPED) algorithm and a Multi-Population Adaptive Inflationary Differential Evolution Algorithm (MP-AIDEA). Multi-Objective Parzen-based Estimation of Distribution (MOPED) belongs to the class of Estimation of Distribution Algorithms (EDAs) and it is used in the first phase of the initial guess research to explore the search space, then Multi-Population Adaptive Inflationary Differential Evolution Algorithm (MP-AIDEA) is used to refine the obtained results, and better fulfill the imposed constraints. The initial guesses obtained with this evolutionary framework were tested on different multiple shooting configurations. The importance of the continuity properties of the employed mathematical models was also quantitatively addressed.


Introduction
In recent years the aerospace industry has seen an increasing interest towards reusable Space Access Vehicles (SAVs), thanks also to the efforts of private initiatives such as Space X and Blue Origin, justifying new research on Space Access Vehicles (SAVs) that can be meaningful both for the industry and the academic community.
The optimal control study is an important phase in the trajectory design of an ascent vehicle since it allows the engineers to fully understand what is the optimal trajectory to perform a certain type of mission considering the constraints imposed by the structure and other components of the vehicle. Different optimal trajectories can be formulated for different objectives, for example a reduction of the propellant consumption can be desired or an increase in the reached altitude.
Optimal control methods can be divided as direct and indirect and a more detailed description of such categories is given in Becerra (2012) while more on indirect methods can be found in Longuski et al. (2010). The approach considered in this work belongs to the category of direct methods with a Multiple-Shooting (MS) transcription and the numerical optimization method chosen in this work is the Sequential Quadratic Programming (SQP) ( Nocedal and Wright 2006) method implemented in the Matlab fmincon library.
Different numerical optimization methods can be found in the literature and as clearly summarized by Conway (2012), they all have different pros and cons. Derivative based Non-Linear Programming (NLP) solvers working with transcription methods are the most adapt to deal with problems with strict constraints and a large number of variables, as the one considered in this work. Although, by being derivative based, they heavily rely on the provided initial guess which must be close enough to the global optimum for the solver to find it. On the contrary, Evolutionary Algorithms (EAs) are very good at finding the global optimum but their computational load increase dramatically with the number of optimization variables and they cannot deal efficiently with strictly constrained problems.
The aim of the proposed work is twofold: on one hand it provides a comprehensive framework which can be employed for other optimal control studies with special insight on the importance of the continuity properties of the employed mathematical models. Such framework relies on an initial guess generation strategy based on the combination of 2 different Evolutionary Algorithms (EAs). On the other hand, this work will further investigate the optimal control study of the ascent trajectory of a Single-Stage To Orbit (SSTO) reusable space access vehicle with a spaceplane configuration, performing a vertical take off from the ground level. Such a vehicle was chosen to highlight the potential that resides in its design. In fact, a SSTO vehicle is more suitable to be efficiently reused, hence a greater reduction of the costs can be achieved, in contrast to reusing a stage of a launch vehicle. This study will also further investigate a potentially good SSTO vehicle design that was discarded in the past due to the limits in the technology of that time.
The initial guess generation or the solution of the optimal control problem using an EA is not a new topic in the literature. Different flavours of Evolutionary 1 3 Single-stage to orbit ascent trajectory optimisation with… Algorithms (EAs) were used, for example Genetic Algorithm (GA) was used by Yokoyama (2002) to select the initial solutions for the Sequential Quadratic Programming (SQP) solver and by Al-Garni and Kassem (2005) to solve the complete trajectory optimization problem. Giri and Ghose (2010) used a Differential Evolution (DE) algorithm to solve the optimal control problem comparing their approach to a Direct Collocation Transcription Method. While in Maddock and Minisci (2016) an Adaptive Inflationary Differential Evolution algorithm is used to generate an initial guess which is then fed to a Non-Linear Programming (NLP) transcribed in a Multiple-Shooting (MS) manner and solved with a gradient based solver. Aside from evolutionary initial guess generation, different approaches can be found in the literature such as the one presented in Toso and Maddock (2017) where a variant of the multi-start method using randomly generate values with Latin Hypercube Sampling and Mersenne Twister is presented; or as the one in Gath et al. (2002) were an indirect method on a simplified version of the considered problem is used to provide an initial guess for the complete problem.
In this work, a combination of two Evolutionary Algorithms (EAs) and a gradient based solver was used. The approach used in this work consists of the use of two different Evolutionary Algorithms (EAs); the first one is Multi-Objective Parzen-based Estimation of Distribution (MOPED), an algorithm which belongs to the class of Estimation of Distribution Algorithms (EDAs) and it was employed to explore the search space by solving a Non-Linear Programming (NLP) problem transcribed in a Single Shooting manner in order to have a small number of optimization variables. Then (MP-AIDEA) was used to refine the obtained results to better exploit the found minima. The Single Shooting transcription solves the same problem of the Multiple Shooting one, but when Single Shooting was used only inequality constraints were considered and some extra constraints were added to guide the evolutionary algorithms towards more meaningful solutions.
Unlike the majority of the other publications on the topic, the proposed work aims also to statistically demonstrate and analyse the robustness of the employed Evolutionary Algorithms (EAs) based initialisation approach. Moreover, the work aims to practically and quantitatively investigate the influence of the setting of the Multiple-Shooting (MS) approach on the achievable performance of the vehicle.
The paper is organized as follows: in Sect. 2 the problem is presented and the chosen mission profile is described (given the general nature of the journal, all details on models are given to facilitate the reading); in Sect. 3 a brief description of optimal control and the chosen transcription method is provided. In Sect. 4 the employed strategy for the initial guess generation is presented; and the results about the initialisation process and the Multiple-Shooting (MS) setting are then discussed in Sect. 5.

Problem description
The considered vehicle, depicted in Fig. 1, is a lifting-body concept vehicle capable of performing vertical takeoff and horizontal landing powered by an aerospike engine. Its mathematical models were used and described in DAngelo et al. (2000).

3
Although this launcher is a quite old concept, such a vehicle could still be of interest. In fact, by being reusable, it would allow for a cost reduction in comparison to classical staged rockets and to staged rockets with a reusable first stage. Moreover, the aerospike propulsion system would allow to have only one engine without the need for additional boosters as was done for the Space Shuttle. Besides, although not being analyzed in this work, the reentry trajectory would benefit from the lifting body design and the control surfaces, allowing for a greater controllability.
Equation (1) describe the dynamical motion of hypersonic vehicle modeled as a rigid body, flying over a spherical and rotating Earth and no wind is considered in the atmosphere. These equations are slightly different from their standard formulation as presented in Vinh (1993). In fact they are usually obtained considering the origin of the heading angle either in the north or in the east. In this work, to be in line with what was done in DAngelo et al. (2000), the equations of motion were derived considering a heading angle = 0 deg going westward and positive from west to east in clockwise direction. This set of equations is completed by the mass variation equation. In this work, the contribution of the bank angle is considered negligible, thus = 0 deg is assumed. This was done to simplify the considered optimization problem by reducing the number of control variables. 1 3 Single-stage to orbit ascent trajectory optimisation with… In order to avoid numerical problems which arise when the flight path angle = 90 deg , a perfectly vertical flight is considered for the first 2 seconds of the flight. Hence the simplified set of equations shown in Eq. (2) is used for this phase of the flight. Then at 2 seconds was set to 89 deg and the full set of equations of motion in Eq. (1) was used.
The use of these two different sets of equation of motion introduces a discontinuity which is a problem for the Non-Linear Programming (NLP) gradient based solver. Hence, both the initial guess generation and the Multiple-Shooting optimization have been started at t = 2s , with the initial values of the states variables as reported in Table 1. The initial values for t = 2s in Table 1 were obtained by propagating Eq.
(2) with the ODE45 solver implemented in Matlab, starting from the initial values at t = 0s.
About the initial values in Table 1, the longitude and latitude were chosen as the one of the Kourou launch site and the initial heading angle is considered as an optimization variable in order to give the optimizer the possibility to find the optimal initial heading angle. Therefore, the initial value of the heading angle is considered unknown as reported in Table 1.
The time dependant quantities used in Eqs. 1 and 2 and in the rest of this paper are listed in Table 2, while the constant parameters are listed in Table 3.
The control variables of the considered configuration are the angle of attack and the linear throttle . About the employed propulsion model, the linear throttle is used to evaluate the axial thrust T x and the thrust angle listed in Table 2 is evaluated as = tan −1 (T z ∕T x ) + . Regarding the other employed models, the aerodynamic model is defined by a lookup table composed by data obtained experimentally. The original database provides c l , c d and c m as a function of Mach, and the body-flap deflection f . In this work, the body-flap deflection was not used as a Single-stage to orbit ascent trajectory optimisation with… control parameter hence the database was reduced removing those data dependent on f . The atmospheric model used is the U.S. Standard Atmosphere 1962 model. These models are the same used in DAngelo et al. (2000), hence please refer to it for more information.

Mission profile
The mission profile for the considered vehicle, schematised in Fig. 2, is composed by a propelled ascent phase until a circular parking orbit is reached and then a Hohmann transfer is performed to bring the SAV from the parking orbit to the target circular orbit at 400 km altitude with an inclination of 51.6 deg. The Hohmann transfer was used to reduce the computational cost of the optimization problem. In fact, without using the Hohmann transfer, the equations of motion should be propagated from ground up to an altitude of 400 km. On the other hand, using the Hohmann transfer the equations of motion can be propagated only until the initial parking orbit is reached. It must be pointed out that the radius of the parking orbit is considered as an optimization variable and thus it is defined during the optimization process. The Hohmann transfer is calculated according to Eq. (3), where h target is the height of the target orbit and GMe is the gravitational parameter.
The final velocity and final heading angle were considered as constraints in order to have a successful insertion in the desired orbit. In this context, also the contribution of the rotational speed of the Earth was considered in the final orbital velocity. Eq.
(4) describe the velocity and heading angle needed to perform a correct orbit insertion, while Eq. (5) show how the velocity components along the 3 wind axes are evaluated for the evaluation of the final absolute velocity, which was then compared to the needed orbital velocity. The final absolute velocity was evaluated as Regarding the other states quantities that appear in Eqs. (4) and (5), they were evaluated using Eq. (1). (3) Single-stage to orbit ascent trajectory optimisation with… 3 Optimal control

Transcription method
As introduced in Sect. 1, the trajectory optimization problem treated in this work is tackled with a direct optimal control method, namely a direct Multiple Shooting method. To be able to use a direct optimal control method, the original trajectory optimization problem must be transcribed into a nonlinear programming problem, whose generic constrained formulation can be defined as: where J is the objective function to be minimized, x are the states and u the controls. L and L are the lower bounds of the path constraints for states and controls while U and U are the upper bounds. 0 and 0 are the initial conditions for the states and controls while f and f are the final conditions.

Multiple-shooting
The direct multiple-shooting approach was introduced by Bock and Plitt (1984) and it is summarized in Fig. 3. In this approach, the total time interval [t 0 , t f ] , where t 0 is the initial time and t f is the final time, is divided into N subintervals, called legs, and then a single-shooting method is used over each leg [t i , t i+1 ] for i = 0, ..., n − 1 , where n is the number of legs. Such single shooting is performed propagating the equation of motion starting from the values of the states at the beginning of each leg, (t i ) , and using as control inputs the values obtained from a polynomial interpolation done over the control points, provided by the optimization procedure in the considered leg. In this work the interpolation is performed using a Piecewise Cubic Hermite Interpolating Polynomial (PCHIP).
In order achieve the continuity of the states and control variables, the additional equality conditions are applied at the interface of each leg: The continuity conditions of Eq. (6) are then set as equality constraints for the Non-Linear Programming (NLP) solver, which is aimed at reducing the value of the dif- refer to the first points of the next leg and x(t + i ) and u(t + i ) refer the last point of the considered leg. The multiple shooting approach can be considered as a compromise between the single shooting method, which has a limited number of variables but also an excessive sensitivity with respect to control variables at the start of the shooting, and the collocation methods, which do not suffer the excess of sensitivity of some variables, but have to deal with a much higher number of variables ( Rao 2009, Becerra 2012, Böhme and Frank 2017and Kelly 2017. For this work, the continuous problem has been formulated as in Table 4. According to Table 4 and the adopted Multiple Shooting transcription, the set of variables to be optimized are listed in Eq. (7). Where 2s is the initial heading angle; i = [V, , , , , h, m] i with i = 2, ..., n l is the array of the states variables and n l is the total number of legs. jk = [ j , j ] k is the array composed by the control variables inside each leg, where j = 1, ..., n cp , k = 1, ..., n l and n cp is the number of control points inside each leg; t m with m = 1, ..., n l corresponds to the size of each leg, in seconds. Regarding the states variables i , these are composed only by the initial values of the states in the legs, starting from the second one to the second-last. Hence, the bounds on the states variables listed in Table 4 are applied only to the initial values of the states in the aforementioned legs. dim( MS ) in Eq. (7) denotes the size of the array composed by the optimization variables for the Multiple-Shooting optimization. Single-stage to orbit ascent trajectory optimisation with… 4 Adopted strategy for initial guess generation As stated in Sect. 1, by being gradient-based the Non-Linear Programming (NLP) solver requires an initial guess close enough to the global optimum in order to fully exploit it and produce good results. Unfortunately, it is quite challenging to produce a well-posed initial guess for an ascent trajectory optimization since the considered systems are defined by a set of non-linear equation of motion influenced by many other different parameters such the aerodynamic coefficients and the thrust models. Hence, even the random creation of an initial guess would require many iterations to provide a proper solution.
(7) MS = [ 2s , i , jk , t m ] with dim( MS ) = 1 + n states (n l − 1) + n controls n cp n l + n l Because of that, the strategy adopted in this work to find a suitable initial guess make use of Evolutionary Algorithms (EAs), which do not require an initial guess themselves and can solve a simplified version of the considered problem, with less optimization variables and more loose constraints. In particular, the chosen strategy is the same as the one adopted in ( [ N cl,l , N cl,u ]. An MP-AIDEA process is then performed, beginning with the sub-population of each cluster. Each process is stopped only when the population contracts to below a predefined threshold. For this work, the code of the Multi Population version of the Adaptive IDEA (MP-AIDEA) (Di Carlo et al. 2019), integrated with features to deal with constraints has been used.

Results
In this section, the obtained optimal solutions are presented. The importance of the continuity in the aerodynamic models is analysed and the robustness of the employed approach to generate a proper initial guess is assessed. Multi

On the continuity of the used experimental databases
Before proceeding to analyze the initial guesses and the final optimal solutions, it is important to discuss the importance of the continuity properties of the experimental databases used in the models involved in the optimization process. It is common that for aerospace applications, experimental databases are created to describe, for example, the aerodynamic properties of the studied vehicle. The data in such databases are noisy, and most of the times generalised through piecewise linear functions, introducing discontinuities that are an issue for the derivative based Non-Linear Programming (NLP) solvers used to treat direct optimal control problems (Betts 2010). In this work the noise on the data has been neglected and generalisation has been improved by using cubic interpolation. A not-a-knot cubic spline interpolation was performed on the original aerodynamic coefficients data, in order to increase the number of points in the database and improve the smoothness of the models. This approach is defined as Cubic Spline Generalisation (CSG). In particular, the number of points in the original aerodynamic database was increased from 221 (17 × 13) to 10000 (100 × 100) using a not-a-knot cubic spline interpolation. Then, when the database was called, a linear interpolation was performed on these data to extract the desired coefficients. This approach was chosen because it was observed that performing a cubic interpolation on the whole database every time the aerodynamic model was called, increased greatly the computational time. The C2 continuity cannot be guaranteed if the desired coefficients are not extracted by performing directly a cubic interpolation when needed. Nonetheless, it was observed that a significant increase in the number of points in the database improves the smoothness of the models. It is not the aim of this work to thoroughly investigate different methods to guarantee the continuity of the used experimental databases, but this could be considered for future works. Other methods can be found in the literature, as variants of the cubic interpolation for example as the Minimum Curvature Spline presented in Betts (2010), or as surrogate models e.g. the one proposed by Ricciardi et al. (2019).
In the following sections, the results obtained with and without CSG are presented. The analysis of the results obtained with and without the CSG regards the aerodynamic database, giving lift, drag and moment coefficients as function of the angle of attack, , and Mach number, M. In both cases, with and without CSG, when the aerodynamic model was called, a linear interpolation was performed to extract the desired coefficients. Hence the CSG and the nominal approaches differ only in the number of points in the dataset.

Initial guess generation
For the initialisation part, the MOPED optimization was set as: size of the population, N = 100 , the fitness parameter f = 0.5 , the sampling proportion = 1 , and the maximum number of generation, N g = 100 . MP-AIDEA was set as: number of populations, N p = 1 , number of individuals per population, N = 20 , contraction factor, l = 0.25 , and no-adaptivity with the weighting factor, F = 0.6 and the crossover probability, CR = 0.9.
The single shooting transcription used for the initialisation process was set with 25 control points and 25 inequality points, which are the points used to evaluate the inequality constraints. Both the control points and the inequality points are evenly distributed inside the leg. (Note that this is the same also for each leg of the multiple shooting approach). The set of optimization variables is defined as in Eq. (8). These are: the initial value of the heading angle 2s ; the control points over the whole trajectory j ; and the final time t f . As in Sect. 3, j = [ j , j ] is the array composed by the control variables, where j = 1, ..., n cp and n cp is the number of control points. dim( IG ) in Eq. (8) denotes the size of the array composed by the optimization variables for the initial guess generation.
The search space for the optimization variables in Eq. (8) is summarized in Table 5 along with the applied inequality constraints. As briefly stated in Sect. 4, a set of relaxed constraints were considered during the initial guess research. They are relaxed in the sense that no equality constraints were used, but the same inequality constraints used in the MS transcription were used plus few extra inequality constraints formulated to guide the initial guess research towards more meaningful solutions.
(8) IG = [ 2s , j , t f ] with dim( IG ) = 2 + n controls n cp 1 3 Single-stage to orbit ascent trajectory optimisation with… Such constraints were set in a way that allows the EA solvers to find an initial guess close to satisfying the orbital requirements to perform the Hohmann transfer, which was then used to obtain the value of the objective function.
Since the initial guess generation process involves two stochastic algorithms, multiple initial guesses were generated in order to assess the robustness of the approach. Figure 4 shows the comparison between the initial guesses obtained by  . 4 Comparison of the initial guesses produced with and without the CSG approach. The initial guesses are represented in terms of Objective Function and Feasibility values. Since a feasibility of 0 cannot be depicted on a logarithmic scale, a feasibility of 10 −16 was assigned to those initial guesses with a feasibility of 0 Multi-Objective Parzen-based Estimation of Distribution (MOPED)+MP-AIDEA with and without using the CSG approach, in terms of objective function and feasibility. The objective function, to be minimised, is the opposite of the ratio between the final mass and the initial wet mass, as also shown in Table 4. The feasibility represents the maximum value of the constraints violation. Hence, the lower the feasibility the more the solution respects the considered constraints and bounds.
To obtain the 44 initial guesses (22 with CSG, 22 without CSG), depicted in Fig. 4, MOPED+MP-IDEA was run three times using the CSG approach and four times without using it. This difference is due to the stochastic nature of the used algorithms.
As can be seen in Fig. 4, despite in the case without CSG the initial guesses possess a lower fitness value and more initial guesses with a feasibility value of 0 ( 10 −16 in Fig. 4) were produced, they are more scattered in comparison to the initial guesses obtained with the CSG approach. This means that the CSG approach slightly modifies the database to have marginally worst performances, but this is compensated by having also a more reliable initialisation. The initial guesses produced using the CSG approach are listed in Table 6 and depicted in Fig. 5. Here, IG 6 is the best initial guess found since it is the one with the best objective function value among those satisfying the constraints. Therefore, IG 6 was chosen for the subsequent Multiple-Shooting optimization. Single-stage to orbit ascent trajectory optimisation with… To further prove the reliability of the initial guesses generated with MOPED+MP-AIDEA using the CSG approach, in Fig. 6 the optimal solutions obtained after performing a multiple shooting optimization with 2 different randomly chosen multiple shooting configurations starting from the produced initial guesses were compared. The chosen multiple shooting configurations are: one with 7 legs, 5 control points per leg and 30 inequality points per leg, which will be called configuration 7-5-30 from this point onward; and the other with 12 legs, 4 control points per leg and  Table 6 for the details of the compared initial guesses Fig. 6 Optimal solutions obtained with the 7-5-30 and 12-4-30 configurations starting from the initial guesses with and without the CSG approach. Each point is an optimal solution generated by a different initial guess 30 inequality points per leg, which will be addressed as 12-4-30 configuration. As in the previous comparison, also here the optimal solutions found using the CSG approach were grouped in a smaller cluster if compared with the optimal solutions found without using the CSG approach. Moreover, after performing the multiple shooting optimization, the CSG approach produced better optimal solutions with a lower (better) objective function value and lower (better) first order optimality and these observations were confirmed on both the tested configurations.
From these data it is clear how the initial guesses obtained with the CSG approach led to similar final solutions, which can be further improved by optimizing the Multiple-Shooting (MS) configuration.
After this analysis it can be concluded that the CSG approach led to more robust initial guesses, which are listed in Table 6 and depicted in Fig. 5 and better final optimal solutions. Among these initial guesses, the initial guess IG 6 obtained using the CSG approach was selected for the Multiple Shooting optimization. IG 6 is the one highlighted in Table 6 and is the solution represented by the thickest line in Fig. 5. This initial guess was chosen since it is the one with the best objective function value among the initial guesses that completely satisfy the constraints.
For the plots of the state variables and the constrained quantities profiles related to the initial guesses found with the CSG approach see "Appendix A".

Multiple-shooting results
Once the best initial guess was available, an investigation into the best multiple shooting configuration was conducted and the results are listed in Table 7 in "Appendix B". The best Multiple-Shooting (MS) configuration is the one with 8 legs, 5 control points per leg and 30 inequality points per leg, where the inequality points are those points used to evaluate the inequality constraints. The optimal profiles of the states variables, control variables and the constrained quantities are depicted in Figs. 7, 8, 9 and 10 (the different colors correspond to the different legs of the Multiple-Shooting (MS)). As it can be observed in Fig. 7b, the throttle profile shows an atypical drop, compared with the optimal throttle profile for a classic rocket with solid propellant. But, in contrast to the aforementioned configuration, the studied vehicle possess an advantage which is the possibility to control the throttle level during the flight. Such feature, allows the reduction of thrust in the first part of the ascent phase in order to keep the dynamic pressure within its upper limit. This throttle control profile was found to be similar to that of the Space Shuttle (Edge andPowers 1976, Lu et al. 2003).
The best optimal solution found have an objective function value of J = −0.122976 which corresponds to a payload of 2.3% ( ∼ 10,348 kg) of the starting mass, for the mission described in Sect. 2.

Conclusions
In this work a reliable approach to solve an open loop optimal control problem for an SSTO ascent trajectory is provided. Such an approach consists of an EA based initialisation, followed by a Multiple-Shooting (MS) treated via a derivative based  NLP solver. The work investigates the effect of the smoothness of data and models, as well as the influence of the different Multiple-Shooting (MS) settings on the effectiveness and reliability of the suggested approach. In particular, the work points out that the smoothness in the data and the models is of great importance to increase the robustness of the EA based initialisation which as a consequence improves the robustness of the MS solutions found. Moreover, it was pointed out that the setting of the Multiple-Shooting (MS), in terms of number of legs and control points, changes the problem landscape in a way that is not linear and not trivial to predict. As a final result of this work, a physically meaningful optimal trajectory, which respects all the imposed constraints, was found. Such trajectory allows to obtain a mass of the payload of 2.3% ( ∼ 10,348 kg) of the starting mass, which is comparable to the payload mass of an medium launch vehicle like the Soyuz-2 to LEO.
Future works will investigate different approaches for the initialisation and the setting and treatment of the Multiple-Shooting (MS) method. In particular, (1) the MOPED+MP-AIDEA approach will be compared to one where only the multi-population MP-AIDEA will be used, in order to understand which could be the most reliable way, and (2) if the use of a different transcription approaches and optimisation solvers can allow for an actual improvement in the results for systems and solutions considered 20 and more years ago. This latter could justify an interesting reevaluation of old concepts and methodologies that were discarded in the past because the technology and tools were not appropriate. Moreover, more work will be carried out to quantify the effect of other smoothing approaches for modelled data.

Comparison of the initial guesses of states and constraints found using the CSG approach
See Figs. 11, 12 and 13. Single-stage to orbit ascent trajectory optimisation with…  Table 7. Single-stage to orbit ascent trajectory optimisation with… Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/ licenses/by/4.0/.