Optimization of the Conceptual Design of a Multistage Rocket Launcher

: The design of a vehicle launch comprises many factors, including the optimization of the climb path and the distribution of the mass in stages. The optimization process has been addressed historically from different points of view, using proprietary software solutions to obtain an ideal mass distribution among stages. In this research, we propose software for the separate optimization of the trajectory of a launch rocket, maximizing the payload weight and the global design, while varying the power plant selection. The launch is mathematically modeled considering its propulsive, gravitational, and aerodynamical aspects. The ascent trajectory is optimized by discretizing the trajectory using structural and physical constraints, and the design accounts for the mass and power plant of each stage. The optimization algorithm is checked against various real rockets and other modeling algorithms, obtaining differences of up to 9%.


Introduction and State of the Art
Since the beginning of the space age, in 1957, access to space has become widespread.Although in 1966, only the USA, France, and the USSR had the capacity to send satellites into space using their own launch vehicles, today, the list of countries with such capacity has increased, including those involved in the Europe's space agency, China, Iran, Ukraine, India, Israel, North Korea, Brazil, Japan, and New Zealand [1].Manned missions have also received significant interest from both governments and the private sector: both world powers, such as China or the USA, have shown interest in taking human beings to new frontiers such as the Moon or Mars, and private companies, such as Virgin, SpaceX, or Blue Origin, are offering commercial flights around the Earth.
A rocket or launch vehicle carries a payload beyond Earth's atmosphere, powered by a jet rocket engine.Depending on the mass of the vehicle, they can be classified in the following categories, shown in Table 1 [2]: To increase launch efficiency and reduce costs, the following alternatives have historically been addressed:

•
The use of a multistage rocket, as proposed by, among others, Tsiolkovsky.A SSTO (Single Stage to Orbit) would consume more fuel, but if complete reusability could be achieved, as in the case of an airplane, the cost savings would be significant.
Unfortunately, the efficiency of the engines prevents their use in these scenarios.Multistage rockets increase performance, eliminating part of the mass that is no longer necessary during the ascent to orbit: By dividing the rocket into multiple zones or stages, we can reduce the size of these losses.Each stage will, thus, have a fraction of the total propellant, which can be optimized to maximize performance.Another advantage comes from the type of propellant: the first stages, being subjected to higher atmospheric pressures, have to withstand a greater aerodynamic drag, together with a substantial modification in the propulsive performance of the nozzle.High pressures present in the early stages require a shorter nozzle geometry, while in vacuum or at high altitudes, a more elongated shape is preferable for optimal performance.Therefore, this is one of the most remarkable advantages of the use of several stages since it allows to use different engines, depending on the external conditions, thus improving the overall efficiency of the rocket.Finally, a multistage rocket allows for the use of different propellants, depending on the stage.The first stages tend to use solid or liquid fuel based on hydrocarbons, which can achieve the great thrusts necessary to overcome the gravitational attraction, while being relatively low cost.For the following stages, the required thrust is substantially lower, as sub-orbital speeds are reached.In this way, fuels are chosen that admit a large specific impulse with a low volume.Therefore, in these stages, high efficiencies can be achieved without an excessive increase in cost [3].On the contrary, for multistage rockets, one or more engines will be required at each stage, thus increasing the cost and structural fraction of the whole.The complexity of design and manufacturing also increases, leading to additional reliability issues.

•
Optimization of launch parameters, among others: the mass distribution between the stages, the most suitable engine for each of them, the size and mass of the resulting rocket, the amount of propellant required, and the ascent path.All these factors interact with each other, so the optimal solution will be given by an adequate balance between all the elements.
Optimization of the climb path is of utmost importance in the process.Normally, starting with a vertical ascent from the launch pad, it usually continues with a turn to reach the horizontal speed necessary to start orbiting.The ascent is subject to gravity losses, aerodynamic losses, and steering losses.Controllability is affected by the fact that it is performed only with a vector thrust in most cases, and it undergoes the coupling of rotational dynamics to inertia, natural instability, and the presence of a non-minimum phase behavior [4].In the process of optimizing the climb path, the potential applications include: maximizing the payload that can be sent to orbit, analyzing how design changes affect performance, and determining the launch window [5].The optimization problem is not trivial: on the one hand, it is a non-linear dynamic system; there are discontinuities in state variables, time-dependent forces may exist, and some final conditions may not be explicitly known [6].At the beginning of the space age, the approach to solving this optimization problem was dominated by the indirect methods (using the necessary analytical conditions obtained from the calculus of variations).However, from the end of the 20th century, the growth in the computational capacity of computers has allowed the use of direct methods on increasingly complex problems.The indirect methods use the calculus of variations to determine the optimality conditions of the problem, resulting in a multiplepoint boundary-value problem [7].By the use of the Euler-Lagrange equations, as well as the Pontryagin's Maximum Principle (among others), conditions applying to the adjoint variables are reached.Possible examples of these methods include the indirect transcription and indirect shooting method.
The direct methods state the original optimization problem in an approximative manner, converting such a problem into a non-linear programming problem (NLP).Di-rect collocation, multiple shooting, or the pseudospectral methods are possible examples.Rao [7] and Betts [8] provide comprehensive surveys of the advantages and limitations of direct and indirect methods.
Both methods are still used nowadays, and each one complements the cadences of the other.Focusing on the indirect ones, it can be stated that Gath and Well, who developed an algorithm to provide good initial guesses for subsequent calculations [9], Pontani, with the optimization of an air-launch vehicle [10], and Filatyev and Yanova, with the optimization of branching ascent trajectories in presence of atmospheric disturbances [11], are the most significant, among some others.
Regarding the direct ones, possible examples are Brusch et al., who used the shooting method for the trajectory optimization of the Space Shuttle [12], Jänsch et al. for the Ariane 5 [13], Bollino, with the optimization of the reentry of the Lockheed Martin X-33 [14], or Zotes and Peñas, with the use of a genetic algorithm for the optimization of maneuvers [15].
According to the best of our investigations, only a few software packages provide relevant functionality in this area of optimization of rocket design and its trajectory.FONT-SIZE, developed in the 90s, provides vehicle dimensioning and trajectory determination at the same time [16].INTROS (integrated rocket sizing model), developed for NASA, allows dimensioning of the launch vehicle using MER (mass estimating relationships).POST (program to optimize simulated trajectories), developed in NASA, provides trajectory optimization under a three degrees of freedom approach.These software packages are restricted and cannot be used publicly.ASTOS (aerospace trajectory optimization software) is a project for the ESA that simulates trajectories and optimization for three and six degrees of freedom in many scenarios, including ascent and interplanetary maneuvers; as in the previous cases, it is not public domain.Finally, ASTER form Roscosmos, with an approach based on the indirect optimization of branched trajectories that describe the motion of the rocket, obtains multipoint boundary-value problems to solve.The program has proven to be capable of solving ballistic and control problems, for injection aid, in a flexible way [17].
In this research, we introduce an optimization software for the conceptual design of a launch rocket and its trajectory under various optimization techniques; the key properties are mass minimization for the take-off and payload mass maximization.The design optimization can provide configuration alternatives useful to assess the effect of design choices-especially on the propulsive plant-on the overall performance.The trajectory optimization, on the other hand, is designed to obtain fast results and prove the feasibility of a conceptual design in a preliminary phase, especially for the newer launch vehicles that carry smallsats.The paper is organized as follows: Section 2 provides a description of the mathematical model and algorithms used in this research; Section 3 shows how an implementation of those algorithms applies to different use cases; the conclusions are presented in Section 4, including the subsequent steps that will follow this work.

System Modeling 2.1.1. Cost Model
The Delta-V is used to evaluate the requirement of the mission, which, in turn, allows for calculating the fuel necessary to carry out that mission.This Delta-V needed for a mission can be obtained as: where the ∆V orbit refers to the speed that the rocket must reach in order to orbit the planet in a given orbit, and it is equal to the orbital velocity at the insertion point.Furthermore, the losses or gains related to the rotation of the Earth ∆V r can be expressed as: where V r,h0,φ0 is the linear velocity of the Earth at that latitude and altitude, and A 0 is the launch azimuth.In prograde orbits, the rotation component will have a negative value, implying some help from the Earth to launch.Therefore, these two components can be calculated precisely.
The rest of the components, that is, the aerodynamic loss ∆V d , propulsive losses ∆V p , and associated losses to gravity ∆V g , cannot be calculated until the final trajectory is known.Therefore, due to the inability to calculate the trajectory before knowing the design of the vehicle, they are estimated based on the works of [3,5,[18][19][20]: Lastly, ∆V margin is a safety margin to take into account effects not analyzed or errors in the data.According to [21], it usually uses 2% of the total Delta-V.

Propulsive Model
The propulsive modeling is carried out from the characterization of a rocket engine.The optimization algorithm will have complete freedom to choose the engine, for each phase, from a database of 27 real engines, with its characteristics gathered in Tables 2 and 3  In order to facilitate convergence, the engines are divided among those specific for the early stages of a rocket (characterized by a great thrust) and those typical of the last stages (with less thrust but greater efficiency).Once the algorithm chooses an engine from the data base, its multiplicity must be calculated so that the thrust-to-weight ratio is higher than that indicated by the optimization variable.Therefore, for each stage: T total = T engine n engines (5) where T W is the thrust-to-weight ratio set by the optimization variable and T W e f f ective is the one of the rocket.To calculate the number of motors (n engines ), the ceiling function is used, so the actual thrust-to-weight ratio of the rocket is equal to or higher than that chosen by the optimization variable.
From the number of motors per stage, it is possible to define the diameter of the rocket.This diameter is considered constant and equal to the largest diameter required by a stage.In turn, for each stage, the diameter chosen is the greater of that required by the payload and by the motors, so it is used for the latter results of the packing of circles [25] for the calculation of the minor diameter needed.The diameter chosen is the minimum that the rocket geometrically needs, which does not imply that it is the most desirable.However, in the absence of a common criterion applicable, the smallest possible diameter is chosen.Once the engine is selected, the fuel and oxidant necessary for its operation are defined.This allows the calculation of their masses and volumes, by using the mixing ratio f, as well as the total propellant mass, m p , the densities of the fuel, ρ f uel , and the oxidant, ρ oxidant , of the propulsive plant at each stage.
The propellant tanks are considered in a tandem type of configuration, made up of a cylinder with half an ellipsoid on top of either end.The volume V total is then: where h cilinder is the length of the cylindrical part of the tank, and R is its radius.Moreover, a, b, and c are the semiaxes of the ellipsoid.In our case, a and b are equal to allow for the continuity of the tank.The height h tank is: And the area S is: where e is the eccentricity of the ellipse formed by the cross section.Using Equations ( 11) to ( 14), the height and area of the tanks can be determined, which will be useful for the calculation of the final geometry and the structural mass of the configuration.

Atmospheric and Gravitational Model
For atmospheric modeling, the NRLMSISE-00 (US Naval Research Laboratory Mass Spectrometer and Incoherent Scatter radar for Exosphere) model was chosen.This model was developed by the US Naval Research Laboratory in the year 2000, and it allows for the determination of dozens of atmospheric parameters from the surface up to 1000 km.For this, values for the date, latitude, longitude, and F10.7 flux are needed.Since, in a preliminary design, it is complicated to predict the launch window, the first of July of 2022 is selected as the date and the values of a standard day for the F10.7.
For the variation of gravity, an implementation of zonal harmonic planetary gravity, based on planetary gravitational potential, is used.The zonal harmonic coefficients J2, J3, and J4 are used.

Mass Model
To estimate the structural mass of the different rocket elements, Akin's Mass Estimating Relationships, or MER [26], has been used to heuristically calculate the mass from information on the geometry of the rocket and the propellant tanks, as well as their capacity.

Aerodynamical Model
Aerodynamic characteristics are difficult to predict in the absence of CFD simulations or experimental methods, which are both impractical for an optimization process.
The following simplifications have been considered:

•
The angle of attack is small, typically less than 10 • in atmospheric flight.

•
The flow is uniform around the rocket.

•
The rocket is axisymmetric, which is the case in serial staged vehicles without fins.Thus, only serial rockets are considered on the trajectory analysis.However, if the vehicle has fins, an additional term can be added on (15) to take into account the additional drag, as described in the study of Caporaso [27].

•
The lift is negligible.

•
Induced drag is considered negligible.

•
The variation in drag with the angle of attack is also negligible.
The total drag of the vehicle can then be divided as follows: where C Dfb is the body drag, and C Db is the base drag.C Dbd is the component resulting from the sum of pressure drag and skin friction drag.The first can be expressed according to Caporaso [27]: where l TR is the total length of the rocket, l b the length of the cylinder, l c is the length of the boat tail (not applicable), d b is the maximum diameter, C f is the friction coefficient, S s is the wetted area, and S m the maximum frontal area.The friction coefficient is modeled through empirical laws, such as Schoenherr's, with the compressibility corrections exposed in Hoerner [28].
The base drag can be described, according to Fleeman [29], as: Equation ( 16) is valid only in subsonic regime.Therefore, Caporaso [27] provides the following semi-empirical corrections in the case of a sharp-nosed vehicle: where C DS is the value of the drag coefficient in subsonic regime, and M is the Mach number.
Therefore, in subsonic regime, Equations ( 16) and ( 17) are used to calculate the total drag coefficient.Once Mach 0.9 is reached, Caporaso [27] corrections are used to adapt the coefficient to transonic and subsonic (Equation ( 18)).In regimes greater than Mach 2, Equation ( 18) is still used despite being outside its validity domain, due to the little importance of drag to the heights at which these speeds are reached.
Any standard geometry with the tip not rounded can be considered.This choice affects the total wetted area S S of the rocket.

Algorithm Design 2.2.1. Optimization Algorithm
In the optimization process, the use of optimization algorithms is necessary to solve and speed up the problem.In our research, we chose to use the genetic algorithm to solve the optimization design problem (including the distribution of masses), whilst for the optimization of the trajectory, we selected Patternsearch.This is because both problems are non-smooth and, consequently, the algorithms based on gradients are not a good option.
The genetic algorithm is a metaheuristic algorithm and, therefore, stochastic, and evolutive, which is inspired by the process of evolution that occurs in nature to optimize a certain problem.This algorithm can provide near-optimal solutions, exploring the domain without the need of an initial guess.Due to its nature, a convergence to a minimum is not guaranteed.However, the solution calculated by the genetic algorithm could be refined by the use of a local solver with adequate convergence properties.
Patternsearch is a deterministic algorithm that is not based in gradients.The process is based on analyzing a series of points that are in the neighborhood of the current point, detecting in which one the objective function has a lower value compared to the one in the first point.This algorithm is chosen for the case of trajectory optimization due to its capacity to find optimal, or near optimal, solutions in a reasonable time without providing information about the gradient.According to [30], this algorithm is proven to be convergent under certain conditions.In both cases, MATLAB was used to find the solution.

Optimal Stage Distribution
For each configuration of the design parameters, the optimal distribution of the total mass in each stage is calculated; that is, their size, considering the mass ratio n i , the exhaust velocity V e,i , and the structural ratio i variable at each stage.Under these conditions, a numerical optimization of the problem is needed.Figure 1 displays a summary of the masses described in this section for a parallel configuration: For each configuration of the design parameters, the optimal distribution of the mass in each stage is calculated; that is, their size, considering the mass ratio  , th haust velocity  , , and the structural ratio ϵ variable at each stage.Under these co tions, a numerical optimization of the problem is needed.Figure 1 displays a summa the masses described in this section for a parallel configuration: Considering, for an arbitrary stage ,  , is the payload, which is the part o vehicle that we want to send into space.Normally, there are fixed characteristics s the customer, so from a design-oriented point of view, it is an input. , is the prope mass, consisting of the fuel and oxidant that is consumed to generate thrust.Finally, Considering, for an arbitrary stage i, m pl,i is the payload, which is the part of the vehicle that we want to send into space.Normally, there are fixed characteristics set by the customer, so from a design-oriented point of view, it is an input.m p,i is the propellant mass, consisting of the fuel and oxidant that is consumed to generate thrust.Finally, m e,i is the structural mass, which includes everything that does not belong to any of the previous categories, such as avionics, external structure, or engines.Thus, we may define the mass ratio n and the structural coefficient as follows: The optimization process for a serial configuration can be written as follows: Minimize: With respect to: Subject to: where the first equation of the system (23) is the equation of Tsiolkovsky generalized for a number N of stages.For a parallel configuration, such as the one showed in Figure 1, the vehicle can be divided in the strap-on boosters and in the core vehicle.Thus, a zeroth stage [4] is defined (the parallel boosters and the core first stage are burning simultaneously), and the first stage, in which the boosters have been separated, is when the propellant remaining on the core's first stage combusts.From the second stage onwards, the rocket is considered in the same way as a serial configuration.Therefore, the coefficients and n must be expressed in a different way for stage 0 and 1.For stage 0: For stage 1: For both cases, the subscript b implies the use of propellants, and m p,10 is the propellant burned in the main stage in parallel with the boosters.Moreover, the exhaust speed of the stage 0 is: The optimization scheme now depends on the number of stages.If we consider a three-stage configuration, having the first one boosters, the scheme is now: Minimize: With respect to: where (31) is a system of 11 non-linear equations with 11 unknowns.For the case of two stages, the three equations referring to the third stage are not needed (resulting in an 8 equations system), while for a fourth, we need to add three more equations, resulting in a system of 13 equations.
In our case, the optimization of the number of stages couples with the rest of the optimization modules in order to determine the optimal propulsion plant.The optimization of the mass distribution is performed jointly with the rest of the parameters and not separately.

Design Optimization
The Optimization Process Follows: Minimize: With respect to: n j , I i , TWR i f or j = 1, . . ., N − 1 and i = 1, . . ., N Subject to: where i is the index of the engine used, n Mot is the number of engines in a stage, and TWR the thrust-to-weight ratio.The variable to minimize is the GLOW, since, according to with Koelle [32], it can be related to its cost in the form of CER = aM x , where CER are the cost estimation relationships, x is a constant, and M is the mass at lift-off.The problem is a mixed integer non-linear type.
The flow diagram of the problem is shown in Figure 2.
To initialize the process, the prediction of the structural ratio is required, which will be defined later.The Delta-V is then estimated, which allows us to calculate the mass ratio n by means of Equation (20).Once the type of engines has been selected, we can then determine the number of engines of each stage, along with other parameters, using Equations ( 3)-( 9).This enables the determination of the dimensions of the rocket and propellant tanks using Equations ( 10)-( 13).
Once all the previous steps are completed, the structural mass of the rocket can be estimated more accurately using the MER.At the beginning of the process, arbitrary values for the structural coefficient were chosen, which also impacted the structural mass.Using the MER, we may now check whether this estimation is acceptable (the predicted mass is between 5% and 30% above the MER calculation).If correct, the process is finished, the optimizer is given the GLOW of the chosen configuration, and a different one is tested.If the estimate is not acceptable, it is modified accordingly, and the process starts again with the determination of the mass ratio.
The trajectory is not known, so the Delta-V cannot be recalculated and, therefore, cannot correct the estimates of the gravitational, aerodynamic, and propulsive losses.However, if we impose the same, non-optimal trajectory for all designs, the aerodynamic and propulsive losses due to steering would vary a relatively small amount, since the TWR limits how much different the burnout times will be, specifically, on the first stages, where the losses are higher due to the high value of the flight path angle γ and the thrust T, according to the Equation (35).
Aerospace 2022, 9, x FOR PEER REVIEW 11 of 22 estimation relationships, x is a constant, and M is the mass at lift-off.The problem is a mixed integer non-linear type.
The flow diagram of the problem is shown in Figure 2. To initialize the process, the prediction of the structural ratio ϵ is required, which will be defined later.The Delta-V is then estimated, which allows us to calculate the mass ratio n by means of Equation (20).Once the type of engines has been selected, we can then determine the number of engines of each stage, along with other parameters, using Equations (3)-( 9).This enables the determination of the dimensions of the rocket and propellant tanks using Equations ( 10)-(13).

Trajectory Optimization
This section addresses the task of optimizing the trajectory ascent of a launch vehicle from the ground station to the desired circular orbit.Maneuvers are not considered, so we are free to choose the optimization variables at each discretization point.In addition, a coasting phase, without thrust, is considered on the last stage.This phase will have variable characteristics that are subject to being optimized.The equations that govern motion in three degrees of freedom are the following [32]: In an ECEF (Earth-Centered Rotating Frame) system, θ is the geographic longitude, and ϕ is the geographic latitude.On the other hand, γ is the flight path angle, χ is the heading angle, α is the angle of attack, β is the sideslip angle, ω is the angular velocity of the Earth, r is the radius from the center of the Earth, T is the thrust, and V is the speed.Therefore, we have the following arrangement of variables:

independent monotonically increasing variable: t
With the ligatures between phases: The subindex 0 is for the initial time, and f is for the final moment.The optimization process is, thus: Maximize: With respect to: Subject to: where and M are two integers to choose depending on of the accuracy of the desired discretization and t f is the time when the completion of the ascent happens.t coast,0 and t coast are the initial time of the coast phase, measured since the ignition of the last stage, and the coast duration, respectively.The discretization of the domain follows a linear pattern.
Since the constraints cannot be expressed directly, depending on the optimization variables, they have to be imposed by means of penalties, such as on the Equation (41) or as a set of non-linear constraints to be passed to the solver.
where J is the objective value, η is the penalty, and d is the difference between the actual value and the constraint.The penalty approach was selected since it provides a faster optimization on the selected algorithm.In addition, as the temporal points of verification of the constraints depends on the number of steps used by the integrator, we will have a finite number of relationships.Then, starting with the optimization variables, we must generate a continuous function of class 1 or higher that goes through the aforementioned discrete values and allows the extraction of the value of it, as well as the value of its derivative at an arbitrary point.This interpolation is achieved by means of a piecewise polynomial cubic interpolation "Pchip" (Piecewise Cubic Hermite Interpolating Polynomial).Moreover, in the first 100 m, the rocket performs a vertical ascent since a turn at a lower altitude could make the rocket crash into the launchpad support structure.Therefore, below 100 m, .γ = 0 and .χ = 0. From the system of differential Equation ( 37), we may integrate by means of the MATLAB function "ode45", which uses the Runge-Kutta Dormand-Prince method.In the integration itself, the aerodynamic, atmospheric, and gravitational magnitudes are updated, through their corresponding modeling, as we have discussed previously.In the final step, if we consider a circular orbit, the following conditions must be met: The mechanism to impose the final value state works by conferring a competitive advantage to those solutions that are closer to the desired value.Figure 3  The mechanism to impose the final value state works by conferring a competitiv advantage to those solutions that are closer to the desired value.Figure 3 depicts the bas outline of the algorithm:

Mass Optimization Validation
We proceed to validate the distribution of stages, in parallel, using the Indian GSL Mark I rocket, with data coming from Isakowitz [33].The target orbit is geosynchronou

Mass Optimization Validation
We proceed to validate the distribution of stages, in parallel, using the Indian GSLV Mark I rocket, with data coming from Isakowitz [33].The target orbit is geosynchronous, considering a Hohmann transfer, starting from a circular orbit 200 km high.Under the assumptions from Table 4, GLOW can be calculated as a function of the mass of the thrusters, as seen in Figure 4. Therefore, the optimal value for the mass of the boosters is between 80,000 kg and 120,000 kg.However, on the real rocket, this maximum is not selected for economical or other factors, and instead, larger boosters are preferred.The aforementioned rocket the boosters have a mass of 160,000 kg, obtaining the comparison with respect to the simulation shown in the Table 5.As can be seen in the table, the difference of the calculated GLOW, with resp the real one, is approximately 1000 kg, which means a percentage error of less than The mass ratios are also similar, with a difference of about 9%.The vehicle designed the optimization has the last stage 230 kg larger than the real one, while the rest are sli smaller.

Design Algorithm Application
The design algorithm will be used to generate alternative configurations for ex rockets using serial configuration.A circular orbit at 200 km and a 28.5° inclination i  As can be seen in the table, the difference of the calculated GLOW, with respect to the real one, is approximately 1000 kg, which means a percentage error of less than 0.3%.The mass ratios are also similar, with a difference of about 9%.The vehicle designed from the optimization has the last stage 230 kg larger than the real one, while the rest are slightly smaller.

Design Algorithm Application
The design algorithm will be used to generate alternative configurations for existing rockets using serial configuration.A circular orbit at 200 km and a 28.5 • inclination is considered, and the payload is the maximum available for that orbit in an expendable configuration.Data were obtained from Isakowitz and Spaceflight [32,33].The thrust-toweight ratio during takeoff is between 1.2 and 2.2, so this limit is used in process [21].
The results of the simulation are displayed in Table 6.The algorithm manages to find a configuration that successfully reduces GLOW, which can be related to cost.The 17.9% reduction, in the case of the Falcon, and 32.3%, in the case of the Atlas, it is basically due to the selection of more efficient engines and a distribution of stages with greater importance to the second stages.Variables such as the cost of the engines or propellant are not considered due to the absence of data.
The important reduction of the GLOW, with respect to the real values, is, thus, mainly due to the selection of more efficient motors and propellants.As mentioned, the economic aspects of the launch vehicle are only supposed to be dependent of the GLOW and not of the individual choice of the engine.As can be seen in Table 7, the thrust-to-weight ratio increases in the first stage, accelerating quickly at the beginning, which is a key fact for GLOW reduction.However, if lower g-forces are needed (for satellites, or even with humans, more sensitive to these effects), the value of the corresponding variable of optimization can be further restricted.The structural ratio, which is calculated in an iterative way in the optimization process, finally achieves similar values to the real ones, especially in the Atlas.In the case of the second stage of the Falcon, the value is considerably higher due to the use of four motors in that stage.The diameter, for the same reasons, is almost equal in the Atlas and substantially larger in the case of the Falcon.Lastly, the length is a modified product of the change of the GLOW and diameter.

Trajectory Optimization Application
We proceed, now, to apply the optimization algorithm of the trajectory.The trajectory chosen to be optimized is carried out by the Vega rocket to the Polar Earth Orbit (PEO) at 700 km of altitude.The payload at the end is an optimization variable.The characteristics of the vehicle are shown in Table 8 [34,35].Additionally, it is necessary to establish a series of constraints to prevent unrealistic solutions and to impose boundary conditions on the trajectories.These can be found in Table 9 and are as follows: • The angle of attack and slip angle is limited, for structural reasons, to 6 degrees at altitudes below 100 km and 15 degrees for higher altitudes [36].

•
The dynamic pressure is proportional to the aerodynamic stresses, with the most demanding point being max-Q.For a standard rocket, the limit is 38,300 Pa [5].For the Vega rocket, the limit is 55,000 Pa [5].

•
The aerodynamic load is the product of the dynamic pressure and the angle of attack, and it affects the integrity of the rocket.The considered limit is 230,000 Pa degrees [36].

•
The heat flux .q = 1 2 ρv 3 must not exceed the limit of 1135 W/m 2 [5,37] after the fairing jettisoning.An additional restriction, aimed to address the heating damage on the structure, is also included, with a limit of 40 MW/m 2 [35].The use of these variables to take into account the structural characteristics is a common recurrence in the assessment of trajectories [5,10,35,36,38].For the optimization, a total of 124 total optimization variables (half relative to the angle of attack α and the rest relative to the heading angle β) are considered.We additionally considered three extra variables that correspond to the end time of the ascent path t f , the initial coast time measured since the ignition of the last stage t coast,0 , and the duration of the coast phase t coast .
Figures 5-7 show how the altitude, the velocity, and the trajectory angle converge to their target values at the end point of the path.That is, the height reaches 700 km, the speed reaches close to the orbital speed of 7.50 km/s, and the trajectory angle reaches approximately 0 • , with the payload staying in the target orbit.A comparison of the solutions obtained with the reference data available in the open literature was included [5,35].A digitalization software was used in order to reproduce the results.
their target values at the end point of the path.That is, the height reaches 700 km speed reaches close to the orbital speed of 7.50 km/s, and the trajectory angle reach proximately 0°, with the payload staying in the target orbit.A comparison of the solu obtained with the reference data available in the open literature was included [5, digitalization software was used in order to reproduce the results.The angles of attack and sideslip, as seen in Figures 8 and 9, on the other han under the selected limitation of the Table 9.To obtain more continuous control v The angles of attack and sideslip, as seen in Figures 8 and 9, on the other hand, range under the selected limitation of the Table 9.To obtain more continuous control variables, 246 extra lineal constraints are implemented to have the value of the control variables within a determinate range of the previous value.In this way, big jumps are prevented but not the smaller ones, which are still present.On the coasting phase, there is no thrust, so the angles are irrelevant and will not be considered.
The distribution of the masses with time is, as displayed on Figure 10, not continuous, with low points that correspond to the exhaustion of a stage and its corresponding expulsion.
Table 10 shows the Keplerian elements of the achieved orbit and the target orbit.The differences are small except for the eccentricity that ends in a higher value.An extra maneuver by the spacecraft will have to be considered to circularize the orbit if a circular one is desired.The last three parameters are not specified as objectives and can take an arbitrary value.The angles of attack and sideslip, as seen in Figures 8 and 9 under the selected limitation of the Table 9.To obtain more con 246 extra lineal constraints are implemented to have the valu within a determinate range of the previous value.In this way, but not the smaller ones, which are still present.On the coasting so the angles are irrelevant and will not be considered.phase.
The distribution of the masses with time is, as displayed on Figure 10, not con with low points that correspond to the exhaustion of a stage and its correspondin sion.Table 10 shows the Keplerian elements of the achieved orbit and the target o differences are small except for the eccentricity that ends in a higher value.An e neuver by the spacecraft will have to be considered to circularize the orbit if a circ is desired.The last three parameters are not specified as objectives and can take trary value.Table 11 shows the maximum values that the constraints for the path reach, as well as their limit values.All of them are lower than their respective limits, so the trajectory is structurally and physically viable.The fairing jettisoning should, on the other hand, be done after at least 221 s to be in accordance with the Table 5 constraint.Regarding the available weight for the payload, its value stays at 1431.3 kg, which is a slightly more than 1 kg increase with respect to the real value (a 0.09%).Cos , kun [5] arrives to a value of 1488 kg, whereas Castellini [35] arrives to 1403 kg.

Conclusions and Future Works
This work addresses the optimization of various parts of a launch vehicle design and its mission.The mass distribution between stages is optimized, both in series and in parallel, thus analyzing the effect of the boosters' mass on the GLOW.The overall design of a launch vehicle is also optimized by varying the propulsion plant and mass distribution.The aerodynamic, geometric, propulsive, and mass characteristics are studied and defined from the optimization variables.This algorithm is used in the revision of real rockets, obtaining alternative configurations designed for their usual operation.
Finally, the ascent trajectory of a launch vehicle is optimized by means of direct methods, discretizing the trajectory.Structural and physical constraints are taken into account in order to generate a feasible trajectory.The final conditions are specified in order to achieve an injection into a given circular orbit.Trajectory optimization is applied to the trajectory of the Vega rocket, obtaining close results with respect to other modeling algorithms.
The design of a launch vehicle is a broad area, and therefore, there are multiple lines of development to be pursued.Firstly, the scope of use of the design algorithm can be extended by supporting parallel stages, as well as extending the database of available engines.
As for trajectory optimization, direct methods have been used to carry out the process.The implementation of indirect methods, on the other hand, would allow-through the calculation of variations-a more robust and efficient optimization, increasing the mathematical complexity on the contrary.Additionally, the implementation of partial reusability of the rocket could also be considered by specifying a minimum fuel remaining in the first stage at the moment of ignition of the second stage.
Finally, the logical evolution of such programs is to couple the rocket design with trajectory optimization, so aerodynamic and gravity losses can be calculated dynamically, and the design is optimized for a specific mission.

Figure 1 .
Figure 1.Diagram of a launch vehicle with parallel stages [31].

Figure 1 .
Figure 1.Diagram of a launch vehicle with parallel stages [31].

Figure 4 .
Figure 4. Optimization of the GLOW in function of the boosters' mass.

Figure 4 .
Figure 4. Optimization of the GLOW in function of the boosters' mass.

Figure 5 .
Figure 5. Altitude in the optimized trajectory.The red marks indicate the staging or coasting phase.

Figure 5 .Figure 6 .
Figure 5. Altitude in the optimized trajectory.The red marks indicate the staging or coasting phase.

Figure 7 .
Figure 7. Flight path angle in the optimized trajectory.The red marks indicate th

Figure 6 .Figure 6 .
Figure 6.Velocity in the optimized trajectory.The red marks indicate the staging or coasting phase.

Figure 7 .
Figure 7. Flight path angle in the optimized trajectory.The red marks indicate the staging o phase.

Figure 7 .
Figure 7. Flight path angle in the optimized trajectory.The red marks indicate the staging or coasting phase.

Figure 7 .
Figure 7. Flight path angle in the optimized trajectory.The red marks in phase.

Figure 8 .
Figure 8. Angle of attack in the optimized trajectory.The red marks ind phase.

Figure 8 .Figure 9 .
Figure 8. Angle of attack in the optimized trajectory.The red marks indicate the staging or coasting phase.ce 2022, 9, x FOR PEER REVIEW

Figure 9 .
Figure 9. Sideslip angle in the optimized trajectory.The red marks indicate the staging or coasting phase.

Figure 10 .
Figure 10.Mass evolution in the optimized trajectory.The red marks indicate the staging o phase.

Figure 10 .
Figure 10.Mass evolution in the optimized trajectory.The red marks indicate the staging or coasting phase.

Table 1 .
Classification of launch vehicles according to its mass.

Specific Impulse (sl.) [s] Mixture Ratio Expansion Ratio Height [m] Diameter [m] Mass [kg]
depicts the basic outline of the algorithm:

Table 4 .
Delta-V for the planned mission.

Table 5 .
Comparison between the result of the optimization and the real data for Mb = 160,00

Table 5 .
Comparison between the result of the optimization and the real data for Mb = 160,000 kg.

Table 6 .
Comparison of the optimization variables of the revision, with respect to the real rocket.

Table 7 .
Comparison of relevant quantities of the revised Falcon 9 v1.2 and Atlas V 501, with respect to the real rockets.

Table 10 .
Final orbital elements for the ascend trajectory.

Table 10 .
Final orbital elements for the ascend trajectory.

Table 11 .
Maximum values of the constraints in the trajectory.