Integrated Optimization of Ascent Trajectory and Srm Design of Multistage Launch Vehicles
Abstract
This paper presents a methodology for the concurrent firststage preliminary design and ascent trajectory optimization, with application to a Vegaderived Light Launch Vehicle. The reuse as first stage of an existing upperstage (Zefiro 40) requires a propellant grain geometry redesign, in order to account for the mutated operating conditions. An optimization code based on the parallel running of several Differential Evolution algorithms is used to find the optimal internal pressure law during Z40 operation, together with the optimal thrust direction and other relevant flight parameters of the entire ascent trajectory. Payload injected into a target orbit is maximized, while respecting multiple design constraints, either involving the alone solid rocket motor or dependent on the actual flight trajectory. Numerical results for SSO injection are presented.
19424
1 Introduction
The rapid evolution and innovation that the space industry has undergone in the last few years opened the possibility to reduce satellite dimensions and mass, especially for Low Earth Orbit (LEO) applications. As a main consequence of this new trend, several companies and space agencies have focused their attention on the development of a new family of launchers, here referred as Light Launch Vehicles (LLVs), devoted to the orbital injection of lightweight payloads (i.e., with a total mass up to a few hundreds kilograms). The existence of these new simple and small launchers and the possibility of boarding multiple payloads within a single fairing, have permitted, and are going to permit even more in the future, to break down the cost of space launches.
The European launcher Vega responds to the commercial market requirements for a newgeneration lightweight launch vehicle capable of orbiting small to mediumsize satellites. The Vega family is now going to be enlarged with the development of (at least) two new vehicles: Vega C and Vega E. The objective of such rockets is to increase the load capacity of Vega and improve the overall performance, at more competitive costs. At the same time, the aerospace community is thinking at “scale versions” of existing launch vehicles, devoted to smaller payload transportation into orbit. One notable example is the Vega Cderived LLV here referred as Vega Light, which attempts at exploiting, with minimal changes, already existing and tested Vega motors, in order to reduce development cost and increase reliability. Vega Light is here considered to be a three solid stages launch vehicle composed of Zefiro solid rocket motors (SRMs) of the class Z40, Z9 and ZX (X indicates a few tons of propellant), used as first, second and third stage, respectively. The attitude and orbital control system (AOCS) of the last stage has the duty of compensating SRM performance scattering and carrying out the final orbital injection. Z9 and ZX have been originally designed to operate as upper stages, thus they can be employed directly without any modification. Z40, instead, is now being developed as a secondstage motor for Vega C: for this reason, a redesign of Z40 as a first stage (Z40FS) is mandatory to account for operation in mutated environmental conditions. This can be accomplished modifying as less as possible the present groundtested SRM through an adjustment of the original nozzle and a modification of the propellant grain geometry.
The initial geometrical configuration of the grain plays the main role in determining the performance of a SRM. Once the propellant type and nozzle geometry have been assigned, it is the only available means to achieve a suitable evolution of the burning surface, from which, in turn, the desired pressure and thrust history follow. In the case of complex 3D propellant grains, such as finocyl grains, the design process involves (i) the mathematical modeling of the geometry and (ii) the following evaluation of various independent parameters that define the complex geometry. Changes in the value of each of these parameters bring significant effects on the thrust law, due to finocyl grain particular geometry. In order to ensure that the best possible design (according to some reference performance index) against all achievable configurations is being acquired, the need arises for an intelligent parametric optimization that can control all the design variables.
The problem of optimizing the initial geometry of a threedimensional finocyl grain has been faced by many investigations. Past works dealt with the problem by using a variety of optimization algorithms (sequential quadratic programming[NisarConf2008], simulated annealing[Kamran2011], hybrid optimization techniques[Nisar2008, Adami2017] or hyperheuristic approach[Kamran2012]), but always assuming, as performance measure, a quantity referred to the SRM alone (such as average thrust[NisarConf2008, Kamran2011], specific impulse[Adami2017], propellant[Nisar2008] or motor[Kamran2012] gross mass) and considering constraints related only to SRM geometry or operation (burning time, fixed length or outer diameter of the grain, minimum thrust level, and so on). An attempt to the integrated design and optimization of a complete SRMbased aerospace system, such as an airlaunched satellite vehicle[Rafique2010], a groundlaunched interceptor[Zeeshan2010] or a fourstage launcher[Zafar2013], which contemporarily considers solid propulsion, vehicle masses configuration and flight mechanics, was made, but by using simple models for both grain geometry and trajectory. A concurrent motor robust design and launch vehicle trajectory optimization was carried out by using coupled direct/indirect[Casalino2015] or metaheuristic/indirect[Casalino2014] optimization procedures, but, up to now, applied to a hybrid rocket engine upper stage, with 2D cylindrical geometries for the solid fuel. The simultaneous optimization of the thrust profile and ascent trajectory was already applied to Vega launcher[Civek2017] and to a Vegaderived LLV;[Mancini2018] while the employed dynamical model and thrust law parameterization resemble what is presented here, substantial differences exist in the graindesign approach, flight strategy and optimization technique, because of the different aim of the present work.
The main goal of this paper is, indeed, to present a new methodology for the concurrent first stage propellant grain preliminary design and complete ascent trajectory optimization of a multistage launch vehicle, which has the clear advantage to account for the close and mutual relationships that exist between trajectory and firststage thrust law. In the present analysis, Z40 grain redesign will be accomplished by determining the optimal law for internal pressure during first stage operation, together with the optimal thrust direction along the whole trajectory. Relevant flight parameters (kickoff angle, coasting timelengths, etc.) are also optimized. The aim is to maximize the payload injected into a nominal target orbit, while respecting several path constraints, dependent on the first stage operation (e.g., maximum operating chamber pressure) and the actual trajectory (e.g., maximum dynamicalpressure), and terminal constraints (e.g., assigned orbital radius and inclination). An optimization code based on the parallel running of several improved Differential Evolution algorithms is used to solve the constrained multidisciplinary optimization problem which arises from the formulation employed. Numerical results are presented for a mediumaltitude Sunsynchronous orbit mission.
2 Trajectory Analysis
2.1 Dynamical Model
A threedimensional pointmass model (3DoF model) is considered for the launch vehicle. It has, indeed, the clear advantage to be time efficient with respect to a 6DoF model, while allowing to capture the most prominent features of the mission. Under this assumption, at any time the rocket state is defined only by position , velocity and mass . The equations of motion, in an inertial reference frame, are:
(1)  
(2)  
(3) 
A spherical model is assumed for the Earth; the gravitational acceleration vector is simply given by:
(4) 
where is the Earth gravitational constant.
The expression for the drag force is, instead:
(5) 
where denotes the vehicle crosssectional area and is the drag coefficient, considered here to be just a function of the Mach number .
A simplified atmospheric model is considered; air density , pressure and temperature are evaluated as a function of the altitude , according to U.S. Standard Atmosphere 1976 model.^{1}^{1}1https://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19770009539.pdf Moreover, it is assumed that the atmosphere rigidly rotates with the Earth, with spin rate . Thus, the vehicle velocity relative to the atmosphere is calculated as: . The lift force is here neglected, as it is usually done during the preliminary 3DoF trajectory analysis.
Concerning the thrust force generated by SRMs, it can be expressed as [Wiesel1997]:
(6) 
where indicates the combustion products exhaust velocity, with a mass flow rate equal to , the nozzle exit area, and the ejection pressure. denotes, instead, the thrust in vacuum.
The differential equations and all the involved constants and variables are made dimensionless; in this way, relevant quantities are in a small range around unity, reducing the errors related to the finite digit arithmetic of the computer.
2.2 Path Constraints
In order to ensure the structural integrity of both vehicle and payload during the entire ascent trajectory, a variety of constraints must be set during the simulation phase. Such constraints are, in general, nonlinear functions of the values that the state variables assume, and their main effect is the reduction of the space of the viable trajectories.
In particular, the following path constraints must be respected along the atmospheric flight:
Dynamic pressure  (7)  
Bending load  (8)  
Heat flux after fairing jettisoning  (9)  
Axial acceleration  (10) 
By assuming, for simplicity, that the thrust is always aligned with rocket symmetry axis, the angle of incidence can be evaluated as:
(11) 
Actually, the launcher and the payload are subjected, during the flight, to both static and dynamic loads. For this reason, the socalled quasistatic load (QSL), i.e., the most severe combination of dynamic and static axial accelerations encountered during the mission, might be considered as axial load.[Vega2012] The dynamic loads acting of the vehicle are typically higher during the flight of the first stage, because of the greater thrust level and the interaction with the atmosphere. For this reason, two different values and for the maximum reachable level of axial acceleration have been considered during the firststage flight and the rest of the ascent, respectively. The maximum admissible values assumed for the variables in Eqs. from (7) to (10) are reported in Table 1.
54000  78000  4  5  900 
2.3 Terminal Constraints
In order to ensure the payload injection into a circular orbit of given radius and inclination , the following terminal constraints, acting on final periapsis radius , apoapsis radius and inclination , must be considered:
(12) 
In the present analysis, an injection accuracy of on radius, and on inclination are enforced; the right ascension of the ascending node, instead, can be easily adjusted to the desired value by properly choosing the launch time.
2.4 Flight Strategy
In the case of an aerospace vehicle, devising a flight strategy or guidance program means to schedule the instantaneous direction of the generated thrust with respect to a reference frame along the entire trajectory. As guidance variables defining , the thrust flight path angle and azimuth have here been adopted. Such variables will be referred either to the topocentric frame (topocentric flight path angle and azimuth ) or to the RadialTransverseNormal (RTN, sometimes defined “orbital”) frame (orbital flight path angle and azimuth ), according to convenience. Both reference frames are centered in the instantaneous position of the rocket; they are shown, together with the abovementioned angles, in Figure 1.
The adopted reference guidance program is roughly the same used by the majority of the launchers currently on the market, and, in particular, by launcher Vega.[Vega2012] It is used, as well, with light modifications, in many ascent trajectory optimization problems.[Markl2001, Pagano2010, Brusch1976, Cremaschi2012]. This program consists of a series of specific guidance laws, each to be performed in a determined leg of the trajectory; so, it can also be seen as a way to split the ascent trajectory into phases or arcs. The following guidance laws make up the selected guidance scheme for Vega Light LV.
2.4.1 LiftOff.
The natural way to liftoff is, of course, to fly a vertical trajectory; the topocentric flight path angle is fixed to while, because of the vertical direction of the thrust, the topocentric azimuth is not defined. A vertical liftoff is mandatory since, at the main engine ignition, the generated thrust is the only force able to fight gravity and detach the vehicle from the ground. However, since the liftoff phase has high gravity and misalignment losses, it is flown for the minimum time necessary to reach a safe height above the ramp. Then, the liftoff timelength is fixed to few seconds.
2.4.2 PitchOver and Recovery.
Right after liftoff, the vehicle starts to rotate by deflecting the engine nozzle. As a consequence, the pitch angle (equal to thrust flight path angle in a 3DoF model) decreases (pitchover phase). A simple linear pitchover law is used:
(13) 
The kick angle and the pitchover timelength are optimized. The range in which the optimizer looks for them is reported in Table 2. A constraint is set on the maximum admissible value for the pitchover rate, because of the limited maneuverability of the engine nozzle:
(14) 
with . It will be considered as an additional path constraint.
Concerning the azimuth, a constant value is used in . The optimizer will look for it in a range centered in , i.e., the optimal topocentric azimuth in the case of nonrotating Earth:
(15) 
with the latitude of the launch base, supposed to be in the Azores islands.
A recovery phase is here considered at the end of the pitchover maneuver, in order to progressively align thrust and relative velocity before the gravity turn starts. The thrust direction is fixed in , while the best recovery timelength is determined by the optimization process.
2.4.3 Gravity Turn and Coasting.
The Zero Lift Gravity Turn (ZLGT) or, simply, gravity turn, is a clever maneuver, widely exploited by launchers, during which the thrust vector , the roll axis and the relative velocity are kept aligned to preclude the generation of any side force, and the yaw component of gravity is exploited to rotate slowly the velocity vector downward.[Wiesel1997] The ZLGT guidance law is:
(16) 
The first stage ZLGT (ZLGT 1st stage) is performed up to time , where denotes the burning time of the ith stage. Because of the first stage grain redesign, the time is part of the optimization variables. At this point, the engine shutdown occurs, followed by the first stage separation from the rest of the rocket. A short coasting phase (coasting 12) with fixed duration is considered to allow the separation before the second stage ignition. Along this arc, the vehicle performs a ballistic flight, with null thrust. A ZLGT trajectory is flown by the second stage too (ZLGT 2nd stage), for a total time equal to its burning time (now fixed). When the second stage ends its propellant, it separates from the rest of the rocket. At this time, also the fairing jettisoning occurs. The following ballistic coasting phase (coasting 23) takes place close to the atmosphere limit () and it is characterized by a high relative velocity (). Therefore, it can be exploited to make the vehicle gain altitude, at the expense of a kinetic energy decrease. The coasting 23 duration is determined by the optimization algorithm.
2.4.4 Third Stage Flight.
The third stage ignition marks the end of the atmospheric flight of the vehicle. A guidance law for the flight path angle which can be used, with good results, during the extraatmospheric flight, is the socalled BiLinear Tangent law (BLT).[Markl2001, Brusch1979] Strictly speaking, the BLT is the optimal guidance law under the simplifying hypothesis of flat Earth, plane motion and no aerodynamic forces:
(17) 
A good compromise is to express the value of the three constants in terms of the initial and final flight path angles and a measure of the curvature , with and .[Markl2001] Therefore, the expression of the guidance law for the flight path angle is, in RTN frame:
(18) 
where: is the nondimensional time in the chosen interval. Figure 2 shows the trend of the BLT law for different values of , by supposing and . For , the linear tangent law will result. The optimized variables are reported in Table 2.
Concerning the azimuth, a constant value is assumed in the RTN frame, and determined through the optimization process. The timeofflight of the third stage is imposed equal to its (fixed) burning time: .
2.4.5 Orbital Injection.
The third stage of Vega Light remains connected to the payload until its final orbital injection. A third ballistic coasting (coasting 34) lets the vehicle gain altitude before the final burn is performed. Its duration could vary from few seconds to several minutes, according to the altitude of the target orbit. The final orbital injection of the payload is carried out by means of a single shoot of the AOCS of ZX motor. A simple constant thrust direction in RTN frame, determined by the optimization process, is selected for this final arc. The burning time of the AOCS is evaluated so that its propellant is entirely consumed. A summary of the described guidance laws can be seen in Table 2.
Arc  Phase  Time  Optimization  Variables 

interval  variables  boundaries  
Liftoff  none  none  
Pitchover  
Recovery  
ZLGT 1st stage  
Coasting 12  none  none  
ZLGT 2nd stage  none  none  
Coasting 23  
3rd stage  
flight  
Coasting 34  
Orbital  
injection 
3 First Stage SRM Grain Design
3.1 Solid Propulsion Model
A simplified 1D stationary isentropic model for the expansion process in the nozzle is adopted to evaluate Z40FS performance, i.e., mass flow rate and vacuum thrust , as a function of equilibrium chamber pressure and temperature , combustion products properties (molecular mass , specific heat ratio ) and nozzle geometry (throat area and exit area ), as it is typically done during a preliminary phase[Cornelisse1979].
Under such hypothesis, the expressions for vacuum thrust and mass flow rate are:
(19) 
(20) 
with the universal gas constant and . and denote the vacuum thrust coefficient and the characteristic velocity, respectively. The quality factors and have been introduced to account for losses in performance of the motor due to its nonideal behaviour. Pressure ratio between the exit and entry section of the nozzle is a function of propellant type and nozzle expansion ratio through the relation:
(21) 
The vacuum specific impulse is, instead, by definition: , with the reference sealevel gravitational acceleration. A quasisteady zerodimensional internal ballistic model is here adopted for predicting the physical quantities behaviour within the combustion chamber: gas properties (in particular, the internal pressure ) are considered uniform; in addition, the chamber temperature is considered constant in time.
Nozzle throat erosion has been taken into account. A constant throat radius regression rate has been assumed in the present study, as usual during a preliminary performance estimation.[Kamran2012] This allows for deriving analytic expressions for the consumed propellant mass (see next section). Thus, the following relation can be exploited in order to determine the throat area value at a given time during motor operation:
(22) 
with and a dimensionless regression rate and time, respectively, the total throat radius erosion, considered to be independent on the selected chamber pressure law, and the initial throat area.
3.2 Pressure Law Parameterization
The thrust trend overtime for a first stage is mainly dictated by the several path constraints encountered along the atmospheric flight; complying with this constraints demands a particular shape for the thrust law , that, in turn, can be achieved by an adequate burning surface evolution . Figure 3 shows the main features of a typical first stage thrust profile determined by system requirements.
Five different regions can be clearly distinguished:

The maximum thrust value is attained soon, and limited by the maximum pressure admitted within the combustion chamber, because of structural requirements, often called Maximum Expected Operating Pressure (MEOP).

Subsequently, the thrust decreases, in order to comply with the maximum dynamic pressure constraint.

Once the altitude reached by the launcher becomes significant, the dynamic pressure value starts to drop; thus, an increase of the thrust is possible and beneficial.

Near the end of the first stage flight, the vehicle mass is considerably reduced while the thrust is still quite high. A new decrease of the thrust is, therefore, the only way to bound the maximum value of the axial load .

The last part of the thrust law is the so called tailoff phase; the rapid drop of pressure within the chamber occurs as a consequence of the sudden decrease of the grain burning surface. It may be advantageous to separate the engine and its associated dry mass before the propellant has been fully exhausted.
The optimal shape of the first stage thrust law is, here, attained through a parameterization of the chamber pressure law, which greatly resembles that of thrust. Acting on pressure has two clear advantages: (i) it allows, with respect to a direct thrust shaping, to account for the nozzle throat erosion and the effects of the expansion process; (ii) it maintains the computation cost low, with respect to a direct grain geometry optimization, and reduces the impact that necessary approximations for the SRM internal ballistic may have on the final result. The pressure profile will be modeled by patching together different curves, each represented by a simple time function and described by a small number of parameters. The motor burning time is considered as a free variable too. A parametric optimization process is exploited with the aim of determining the design variable values; the obtained profile is then rescaled along pressure axis in order to comply with the assigned value of the total first stage propellant mass, being the Z40 loading capacity substantially fixed by the external case and central mandrel. In this way, only improvements due to the grain shape, rather than to propellant additional mass, are taken into account.
3.2.1 PPLLLT Pressure Model.
First of all, a dimensionless chamber pressure and a dimensionless time are going to be used.
A pressure model that gives rise to a thrust law close to the typical one for a first stage (see Figure 3) is described by a ParabolicParabolicLinearLinearLinearTailoff trend (here referred as PPLLLT model). The pressure time law is split into six consecutive legs: the first two are parabolic, simulating the 3D burning of the starshaped section of the finocyl grain, the next three are linear, each with its own slope, and the last one is, again, parabolic, approximating the tailoff phase, as shown in Figure 4.
The dimensionless pressure trend is, therefore:
(23) 
being by their definition.
Pressure law during tailoff cannot be imposed arbitrarily, since it must reflect a precise behaviour which occurs while the last kilograms of propellant burn. So, a quite accurate model for tailoff may be necessary to consider a realistic wait period to attain a relatively low thrust before stage separation occurs. Such wait period, in general, could have a significant effect on SRM performance.[Panicker1998] In particular, the initial steep descent of the chamber pressure during tailoff is due to a corresponding reduction of the burning surface, being the propellant next to finish. The final part of tailoff, instead, is often related to liner and thermal protections ablation phenomena[Schiariti2015]; the effect of such residual thrust will be neglected in the present work. Experimental analysis conducted on a firstguess finocyl grain for Z40FS show that a parabolic law interpolates actual tailoff pressure with sufficient accuracy for the aim of the present study. In order not to make the SRM behaviour during tailoff phase too much dependent on the obtained chamber pressure time law, it has been assumed that: (i) the same propellant mass is burnt during such a phase; (ii) the tailoff starts at the same fraction of the total burning time.
In order to guarantee that, during the search for optimum, the condition is always respected, each dimensionless time has been defined in terms of the previous time and a fraction of the remaining time interval (until tailoff). On such basis, dimensionless times in Eq. (23) are evaluated as:
(24) 
The dimensionless time has been evaluated so that it guarantees the derivative continuity at the corresponding time.
The total consumed propellant mass , exploiting Eq. (20) and (22), is enforced as:
(25) 
where:
(26) 
Analytic expressions for , as a function of the pressure law design variables and can be easily obtained. At this point, the maximum value for the operating pressure can be determined:
(27) 
It should be remarked that the constraint on the MEOP must be respected:
(28) 
Because of the assumed hypotheses, the tailoff parabolic law can be determined by imposing the following conditions:
(29) 
where is equal to:
(30) 
In particular, for the sake of simplicity, it has been assumed that, because of the fast dropoff of mass flow rate during tailoff, the nozzle throat erosion is negligible in such phase (i.e., ). Consequently, the following expressions for the the tailoff parabolic law coefficients result:
(31) 
Not all the values for are, however, admissible. Indeed, so that the parabolic trend correctly approximates the real pressure trend during tailoff, the following constraints, acting on , must be set:
(32) 
The former guarantees a convex tailoff parabola, the latter places the first zero of the parabola in .
The design variables referred to the first stage chamber pressure with their respective boundaries, are listed in Table 3. Particular attention must be paid in the selection of the boundaries of each variable. In fact, such boundaries must be: (i) sufficiently wide, so that the search space extension is not limited too much, eventually ruling out some potential optimal solutions; (ii) sufficiently close, to avoid that the algorithm finds putative optimal solutions that, in practise, cannot result in a realistic initial geometry of the grain.
Model  Time  Optimization  Variables 

law  variables  boundaries  
PPLLLT  parabolic  
parabolic  
linear  
linear  
linear  
tailoff  
4 Optimization
4.1 Optimization Problem
Solving an optimization problem means to find a set of design variables which maximize an objective function .
The feasible solution space, i.e., the space within which it is admitted to look for the optimal design variables, is often limited by:
(i) the presence of an eligibility interval for each variable of the solution vector :
(33) 
(ii) the existence of inequality constraints of the form:
(34) 
and, in this case, one is dealing with a constrained optimization problem.
In the launcher trajectory optimization problem that is going to be faced, the objective is to maximize the transported payload: . The problem can be formulated as follows:
(35) 
where:
(36) 
In particular:
(37)  
(38)  
(39) 
Lower and upper bounds for flight/guidance variables and motor variables has been defined in Table 2 and Table 3, respectively.
Concerning the constraints, these are:
Path Constraints:
dynamic pressure  (40)  
bending load  (41)  
heat flux at fairing jettisoning  (42)  
axial acceleration  (43)  
pitchover rate  (44) 
Terminal Constraints:
assigned orbital radius  (45)  
assigned orbital inclination  (46) 
Pressure Law Constraints:
MEOP  (47)  
convex tailoff parabola  (48)  
parabola first zero at  (49) 
4.2 Global Optimization Algorithm
An optimization algorithm based on Differential Evolution (DE), that proved to be effective in other space trajectory optimization problems[Federici2018], has been selected in the present application. DE is a populationbased stochastic metaheuristic algorithm, firstly introduced by R. Storn and K. Price in 1997,[Storn1997] featuring simple and efficient heuristics for global optimization problems defined over a continuous space. Being inspired by evolution of species, it exploits the operations of Crossover, Mutation and Selection to generate new candidate solutions. Because of its good performance on several benchmarks and realworld optimization problems, many researchers all over the world have directed their efforts to further improve the effectiveness of the original algorithm, by devising many variants, collected in Reference Das2011.
In the present implementation four different mutation strategies, among those originally proposed[Storn1997], in conjunction with a binomialtype crossover, have been adopted. By using the current nomenclature[Storn1997, Das2011], the DE schemes here exploited are:

DE/rand/1/bin

DE/best/1/bin

DE/targettobest/1/bin

DE/best/2/bin
In particular, strategies based on mutation of the best individual (strategies 2 and 4) typically show a faster converge rate toward an (often local) minimum, whereas strategies based on randomly chosen individuals (strategies 1 and 3) better explore the whole search space. In order to avoid the manual tuning of DE control parameters, a selfadaptation scheme[Brest2007] has been implemented: the values of the scale factor and the crossover probability are encoded into the individuals that undergo the optimization procedure.
In order to achieve a good balance between the search space exploration and a faster convergence to a good solution, the proposed algorithm involves the creation of different populations, or “tribes”, each located on an “island” of an archipelago, arranged in a radial configuration, as Figure 5 shows. Each tribe evolves independently from the others and features one specific mutation strategy among the four proposed variants (as reported in Figure 5) until a migration is performed, during which each tribe passes its best agents to the “following one” (according to the direction of migration), in which these agents replace the same number of worst agents. Migration tides alternate at each event. The proposed islandmodel allows an easy parallelization[Izzo2009] on architectures with cores. The termination criterion of the optimization process is simply based on the generation number. A partialrestart mechanism, named “Epidemic”, has been adopted to maintain diversity between individuals in each population and avoid a premature convergence on a local optimum.
4.2.1 Constraint Handling.
Evolutionary algorithms are, in their “standard” formulation, unconstrained optimization methods; this fact makes it difficult to understand how to efficiently incorporate constraints of any kind into the objective function during the search. The DE approach proposed by T. Takahama and S. Sakai[Takahama2010], obtained by applying the constrained method to differential evolution, is here adopted. The idea is to decide the selection process through a set of simple rules:

between feasible individuals, the one with the maximum value of the objective function wins;

between a feasible and an infeasible individual, the feasible one is selected;

between two infeasible individuals, the one with the lowest value of the maximum constraint violation:
(50) survives in the next generation.
Without further refinements, this simple approach pushes the evolution process to focus first on the constraint satisfaction and, then, to randomly sample points in the identified feasible region.[Coello1999] Therefore, a constraint satisfaction tolerance for is introduced and decreased during generations according to the rule:
(51) 
with the initial and final value of the tolerance, respectively. In this way, “bad” moves (i.e. toward infeasible solutions) are more likely to be allowed at the beginning, when the tolerance is high and the entire search space must be explored in order to identify promising regions; as the generation number increases, such moves became forbidden, and the search concentrates in the feasible part of the identified “optimal” region. The following values for the above parameters have been used in the present application: , , , .
5 Results
The reported solution has been obtained through an 8island optimization engine, with agents per tribe and a maximum number of generations equal to . In order to increase the confidence on the attained result, the optimization is repeated several times, and the best found solution is taken as putative optimum. A circular sunsynchronous orbit with an altitude of km above Earth surface and an inclination of has been considered as target mission.
The results of the optimization process, in term of optimal values for the design variables, are reported in Table 4. Table 5, instead, lists: the initial and final inertial (ECI) velocity of the vehicle and , the velocity increment provided by the propulsion system , and the cumulative amount of velocity losses (gravitational , aerodynamic and misalignment ) along the obtained ascent trajectory.
[°]  
[°]  
[°]  
[°]  
[°]  
[°]  
[°] 
0.37  7.61  9.46  1.39  0.14  0.68 
Figure 6 shows the altitude profile along the trajectory. A direct ascent trajectory is flown until the third stage ignition, which occurs just outside Earth atmosphere. At this point, the velocity gained during ZX operation enters the vehicle at the perigee of an Hohmannlike (because of the finite thrust) ellipse towards the target orbit. A coasting 34 of half an orbit period lets the rocket reach the apoapsis, where the final orbital circularization is carried out by means of the ZX AOCS. This is the typical mission profile for a mediumaltitude final orbit; a complete direct ascent would, indeed, entail high gravitational losses, because of the great value of the flight angle between inertial velocity and the horizon. It should be noticed that, actually, the considered AOCS is capable of a unique shoot and ZX must cover this deficiency by realizing the first of the Homhannlike transfer; on the contrary, during a typical mission of the Vega launcher, the Hohmann maneuver is completely carried out by AVUM, able to perform multiple ignitions.
Figure 7 shows the topocentric flight path angle during pitchover and recovery (a), 3rd stage flight (b) and final orbital injection (c). The linear pitch law and the constant attitude used during pitchover and recovery, respectively, are clearly visible in image (a); a perfect alignment between thrust and relative velocity occurs just before entering the first gravity turn maneuver, as desired. The bilinear tangent law performed during ZX flight is reported in image (b), whereas the constant orbital flight path angle assumed during AOCS operation is shown in image (c). From Figure 7 one can appreciate that the topocentric flight path angle referred to launcher relative (and, then, inertial) velocity reaches a null value both at the end of the 3rd stage flight and at the end of the coasting 34: this is a clear confirmation that a Hohmannlike transfer brings the vehicle to its final orbit, after ZX shutdown.
Thrust and pressure laws for the first stage are depicted in Figure 8. The depicted quantities, as well as the time axis, have been rescaled so that they belong to interval . The great similarity with the typical first stage thrust (see Figure 3) immediately catches the eye.
Through the initial parabolic trend the thrust reaches its maximum value very early along the ascent (i.e., during pitchover); the great acceleration lets the vehicle acquire, very soon, a relatively high velocity, reducing the amount of gravitational losses (it is, indeed, evident in Table 5 how important gravitational losses are among the velocity losses).
The condition of maximum bending load is reached at this point of the flight; as Figure 8 and Figure 9 clearly display, the thrust value is still quite high while touches its limit value; it is, indeed, convenient to limit the magnitude of the lift force by decreasing the angle of incidence rather then by reducing the rocket acceleration, fundamental in these initial phases of the ascent to maintain gravitational losses low, and, consequently, gain in payload mass. The nonnull bending load , together with the high misalignment and aerodynamic losses, pushes the pitchover rate toward its maximum admitted value: in this way, the desired kickangle is attained with the minimum pitchover timelength.
The parabolic region is followed by a linear decreasing law, dictated by the maximum value admitted for dynamic pressure. It is easy to note from Figure 10 the typical increasingdecreasing trend of the dynamic pressure, as a direct result of the exponential reduction in atmospheric density with altitude and the increasing velocity of the launcher. The socalled maximum q condition is encountered during the ZLGT 1st stage phase and approximately occurs at the beginning of the increasing linear region of pressure law; once such condition is left behind, the thrust level is free to grow again. As the propellant mass is approaching its end, the axial acceleration climbs so rapidly that it is necessary to stop the thrust increment in order not to violate the corresponding constraint (see Figure 11). A final parabolic tailoff drives the thrust toward its minimum; Z40 is detached from the rest of the vehicle at the time the pressure reaches an imposed minimum threshold .
It is possible to note that dynamic pressure, bending load and axial acceleration reach their maximum admitted value, depicted in graphs as a dashed straight line: this is a clear evidence of the optimality of the solution found so far and of the good performance of the optimization algorithm in terms of constraint handling. Indeed, it is foreseeable that the best performance in terms of payload mass should be obtained when the launcher works at its maximum structural capacity, thus exploiting at best its propulsive capability along a trajectory that limits losses. The maximum heat flux after fairing jettisoning, instead, is well below its threshold value; this suggests that fairing ejection could be performed during the second stage flight, maybe allowing for a further improvement of the payload.
Even though, in the present analysis, the geographical constraints concerning the falldown points of the first two stages have not been analyzed, it can be noticed that a southward launch from Azores’ base allows for an initial flight over the Atlantic Ocean, where no land is encountered, as shown in Figure 12.
6 Conclusion
This paper proposes an optimization methodology that offers the great advantage, over “classical” ascent trajectory optimization techniques, to account for the close and mutual relationship that exists between SRM thrust law and flight mechanics. The methodology is quite general and can be implemented, with the appropriate changes related to the launcher configuration, every time the determination of the optimal thrust of firststage main motors or strapon boosters is required, in addition to the mandatory optimization of the rocket trajectory, without the need to perform a complex analysis of the 3dimensional SRM internal ballistic and of the real regression rate of the propellant grain surface. As a postprocess step, further studies have to be performed with the aim of determining the real grain geometry that produces the previouslyobtained chamber pressure history. The analysis can be extended also to upper stage SRM grain design, although the strict path constraints which characterize the atmospheric flight of the launcher first stage have a great influence in modelling the shape of its thrust law.
The chamber pressure time law has been modeled using few analytic functions; this allows for an easy enforcement of the total propellant mass to be loaded into the motor, avoiding the use (i) of additional constraints, which would further limit the search space, and (ii) of quadrature formulae for the approximation of mass flow rate integrals, which worsen the computational load and might bring to dissimilarity in the order of tens of kilograms of propellant with respect to the imposed mass, because of their limited accuracy. Despite the exposed optimization methodology is still at its beginning, solutions presented above seem very promising, in terms of both payload mass and constraint enforcement.
The exploitation of an effective multipopulation metaheuristic optimization code allows a wide exploration of the solution space, and an optimal solution is found in a relatively easy and fast manner. Moreover, handling constraints with the proposed technique does not severely affect the search for the optimal solution, by letting the code explore the entire solution space in order to find a promising feasible region where to move. Obviously, much more can be done in order to improve the performance of the proposed methodology and increase both the reliability of the solutions and the accuracy of the employed model. As an example, different steering strategy for both the third stage and the orbital control system can be analyzed, through the use of more realistic guidance laws which resemble the optimal ones.
The Vega Light performance along the obtained ascent trajectory lets one foresee the competitiveness of such a new European light launch vehicle within an increasingly crowded market. The reuse of reliable and proven technologies, the exploitation of already existing facilities and the knowhow of the expert engineers involved in Vega program, will surely make this new launcher one of the strengths of the upcoming Vega family.