Hybrid Rocket Engine Design Optimization at Politecnico di Torino: A Review

: Optimization of Hybrid Rocket Engines at Politecnico di Torino began in the 1990s. A comprehensive review of the related research activities carried out in the last three decades is here presented. After a brief introduction that retraces driving motivations and the most signiﬁcant steps of the research path, the more relevant aspects of analysis, modeling and achieved results are illustrated. First, criteria for the propulsion system preliminary design choices (namely the propellant combination, the feed system and the grain design) are summarized and the engine modeling is presented. Then, the authors describe the in-house tools that have been developed and used for coupled trajectory and propulsion system design optimization. Both deterministic and robust-based approaches are presented. The applications that the authors analyzed over the years, starting from simpler hybrid powered sounding rocket to more complex multi-stage launchers, are then presented. Finally, authors’ conclusive remarks on the work done and their future perspective in the context of the optimization of hybrid rocket propulsion systems are reported.


Introduction
The history of trajectory and propulsion system optimization at the Politecnico di Torino (POLITO) begins thanks to the efforts of two brilliant Professors, Giuseppe Bussi and Guido Colasurdo. Two of the authors (L. Casalino and D. Pastrone) had the privilege to have them as mentors, colleagues and friends. With an in-depth study, Giuseppe Bussi reworked the Theory of Optimal Control to obtain a synthetic formulation suitable for an engineering approach. His colleague and friend Guido Colasurdo, inspired by this work, started the development of an in-house indirect approach method for trajectory optimization. This procedure and its improved versions become the backbone of a great numbers of optimization analysis carried out at POLITO, being able to deal with both impulsive and finite thrust, ascent trajectories, orbit tranfers, rendez-vous, unpowered/powered/aeroassisted flybys, interplanetary missions and debris removal.
The capabilities of this indirect optimization procedure was tested during the participation to different editions of Global Trajectory Optimization Competition (GTOC). This competition was initiated by ESA in 2005, with the purpose of comparing and developing methods for space trajectory optimization. The POLITO team was the winner of the GTOC2, (https://sophia.estec.esa.int/gtoc_portal/?page_id=15, accessed on 24 May 2021) organized by NASA Jet Propulsion Laboratory (JPL) in 2006, and GTOC6 [1] (in team with the Università di Roma Sapienza) organized by JPL in 2012, and obtained several additional good placements. The outstanding results obtained proved the validity of the procedure and the excellent level of know-how and skills reached at POLITO in such field. In parallel, due to the strong interdependence of propulsion system performance and specific application, the same approach has been used to carry out coupled optimization of

Propellant Combination
The selection of propellants combination strongly affects design, performance and operation of the engine. As outlined by [12], practical known propellants do not have all the desired properties and present some drawbacks. Therefore, propellant combination is usually chosen as a compromise of different factors such as performance, availability, handling/storability, hazards, physical properties of unreacted propellants and combustion properties. Among the others, density of stored propellants, characteristic velocity c * and regression correlation parameters are the most relevant properties for HREs.
For the oxidizer, a few options are viable for HRE, and the following oxidizers have been considered: Liquid OXygen (LOX), Hydrogen Peroxide (HP), and Nitrous OXide (NOX). Oxygen is stored in liquid phase at cryogenic temperature and can be injected in the combustion chamber directly as a liquid or can be gasified before ignition in order to enhance its stability. The vaporization process of LOX adds costs and weight, reducing at the same time the system safety. Among selected oxidizers, LOX gives the highest characteristic velocities, even if maximum values of c * are obtained at a quite low mixture ratio α, making the low regression rate a pressing issue. Unlike LOX, which is completely non-toxic, HP (in water solutions) and NOX are slightly toxic and may decompose in their tank. The use of HP usually requires a catalytic bed, which is an expensive component (due to the rare metal employed in the decomposition reaction) and introduces complexity to the feeding system. NOX does not require any additional subsystem and is directly injected in the combustion chamber due to its slower decomposition rate. Uncontrolled decomposition may still lead to tank over-pressure; nonetheless, NOX is assessed as fairly safe. Despite its low performance, NOX is currently considered to be appealing for different applications thank to its low cost and its self-pressurizing behavior when stored at ambient temperature.
For the fuel selection, the regression rate value represents an important driver: the fuel mass flow is proportional to the product of burning surface by regression rate and low regression rates lead to excessively long fuel grains or multi-port fuel grains (as will be discussed in the following). A mitigation strategy consists in the use of propellant combinations which are characterized by a maximum value of c * at high values of the mixture ratio, thus reducing the relative contribution of fuel flow to the overall propellant flow exhausted by the engine. Unfortunately, this behavior is typical of oxidizers which contain also non-oxidizing species, which have a detrimental effect on propellant combination performance. Most straightforward solutions are based on the improvement of the regression rate value. Paraffin-based fuels form a melting layer on the grain surface which may be unstable, under proper conditions, and a new mass transfer mechanism, called entrainment, enhances regression rate: liquid droplets, coming from the melting layer, enter the boundary layer without enhancing the negative effect of surface blowing on the convective heat flux [13,14]. Similar effects may be obtained by freezing propellants (solid cryogenic hybrids), and using fuel additives, unconventional injectors or geometries. Please note that all of these approaches have some drawbacks, such as higher cost, safety, and environmental issues, and costs vs. benefits have to be evaluated for any given mission/application.
In their works, the authors considered and compared a few different propellant combinations [15][16][17]. An interesting and comprehensive comparison of propellants for hybrids applications can be found in [18]. Six propellant combinations were taken into account in order to find the optimal design of a hybrid-powered sounding rocket. The cost function to be optimized was the time spent above 100 km (i.e., in micro-gravity conditions) for given initial and payload masses. The regression rate correlation coefficients a and exponent n are reported in Table 1. LOX/paraffin-based wax resulted in being the best propellant combination due to its large regression rate, allowing for the use of simpler grain designs and mixture ratios really close to the optimal values which maximize the characteristic velocity c * .

Feed System
As outlined in Section 1, a time-dependent link exists between thrust level and specific impulse: the fuel mass flowṁ F is a function of oxidizer mass flowṁ O , burning area A b and port area A p . Thus, the choice of the feeding system and control is even more crucial than in LREs, since it has a stronger impact on performance. Some researchers proposed approaches to control mixture ratio at cost of increased engine architecture complexity [25]. At the moment, researchers at POLITO analyzed the potential of the more classical one-lever control architecture. The simplest feeding system available for HREs is the pressure-feed blow down type, which is particularly suited for applications where a low ∆v gain (i.e., a low total impulse) is required. In this case, the propellant flow is not controlled but a proper blow down ratio is to be selected to obtain acceptable thrust level variations while keeping the mixture ratio in a proper range. Due to its high vapor pressure (50.4 bar at 293 K, 58.9 bar at 300 K), NOX has an appealing self-pressurizing behavior, reducing tank pressure decrease and improving performance. In fact, when NOX is stored at ambient temperature as a saturated liquid, the pressure decrease due to liquid outflow during engine operation induces evaporation of a portion of liquid NOX. The benefit is limited by the fact that phase change lowers the temperatures of the liquid, thus lowering the actual vapor pressure. In this case, in order to evaluate HREs performance, a proper model of NOX behavior is required to evaluate tank pressure history. The simplest homogeneous model may be not sufficiently accurate. The authors [26] developed a two-phase lumped model based on the spinodal line which is able to give an approximated but still quite accurate evaluation of tank pressure at a low computational cost (see Figure 1). Further details may be found in Appendix A. For liquid propellants without self-pressurizing capabilities, performance may be enhanced by actively restoring the tank pressure among blow down phases, by means of rapid recharges from auxiliary pressurizing tanks. A single-shot repressurization was considered in Ref. [28], as depicted in Figure 2. This technique was used in the liquid mono-propellant blow down system of Magellan spacecraft and guarantees a lower thrust level variation at cost of a slightly increased system complexity.
Better performance and control capabilities can be obtained with continuously regulated feed system which, however, introduces additional components and weight.
The following examples show that the optimal solution privileges the reduction of the aforementioned additional weight by means of a partially regulated feed system, in which a constant tank pressure phase is followed by a blow down phase.
Mass flow rate µ (kg/s) All the above-mentioned gas pressure feed systems are not well suited for high ∆v applications, when high pressure levels are needed in the tank to keep high thrust levels and regression rates. In this case, the oxidizer and auxiliary tanks are may be excessively heavy and the use of turbo pump feed systems can improve rocket performance in a substantial way: propellant tank mass is reduced since both tank pressure and pressurizing gas requirements are low. However, in HREs there is only one liquid propellant, thus it is quite hard to produce an high-energy working fluid to feed a turbine to drive the pump. Possible solutions, such as decomposing mono propellants (e.g., HP) or auxiliary liquid propellants, introduce additional complexity to the feeding system (i.e., higher costs). For this reasons, the use of batteries to power an electrically driven turbo pump is really promising. Details about the models developed to describe such feed systems are reported in Appendix B.
In Ref. [29] a HRE powered upper-stage of a Vega-like launcher was optimized comparing blow down, partially regulated and electric pump feed systems. The chosen propellant combination was HP/PE and, due to the low regression rate, a multi-port grain was adopted (see the following Section 2.3 and Appendix C). The initial mass of the third stage was given (14,522 kg) and the mission target was the insertion of the payload into a 700-km circular polar orbit.
In Figure 3, trends of thrust and mixture ratio are shown. The thrust range variation is reduced when passing from blow down to turbo pump. On the contrary, a larger variation of mixture ratio occurs and specific impulse penalties are expected. Nevertheless, dry mass reductions prevails: the regulated case has a 200 kg gain (more than 10%) with respect to the blow down case and further 200 kg are gained by the turbo pump case. Similar results are found for different applications and propellants. For instance, in Ref. [30], the same problem was analyzed with a robust-based approach, but with LOX and paraffin-based wax as propellants and with a single port grain (thanks to the higher regression rate). In this case the gain of the turbo pump solutions with respect to a regulated case was 250-400 kg depending on the considered technological level of batteries and electric motor, as reported in Table 2. Thrust and mixture ratio histories for blow down, regulated and pump-fed optimized solutions. Table 2. Mass budget comparison. The first column reports the feed system employed: regulated (R) or turbo-pump (TP). (A) and (B) refer to the electrical properties used in the mass evaluation (see Table A1).

Grain Geometry
In HREs the initial grain geometry, i.e., the initial port area (A p ) i and initial burning area (A b ) i , determines the initial propellant flow rate, which in turn gives the initial thrust. During engine burn these areas change as a function of time. Thus, the evaluation of HRE operation and performance requires the knowledge of the grain geometry evolution as a function of time. Namely, relations are more easily written as a function of burnt distance y(t).
The simplest grain design adopts a single-port geometry, with axial symmetry and a circular port. A uniform regression rate is assumed along the grain and the port area A p and burning area A b can be computed by means of Equations (1) and (2), where (R c ) i is the initial cylindrical port radius and L b is the length of the fuel grain. Single-port fuel grains require a large value of L b in order to ensure consistent values of A b and thrust because of the typical low regression rate of hybrid fuels. For this reason, multi-port grains are usually employed. In Ref. [31], the authors compared two types of multi-port grain, namely triangular port and quadrangular port design. Details about the multi-port grain design can be found in Appendix C. Six and eight ports were considered for both triangular and quadrangular geometry, whereas HP/PE was the propellant combination. Both a blow down feed system and a partially regulated feed system were taken into account for the liquid propellant. The case study was the maximization of the payload inserted into a 700-km circular polar orbit by the hybrid powered upper stage of a Vega-like launch system. In this specific application, the difference in the optimal design between triangular and quadrangular port design are relevant when N = 8 was considered. Nevertheless, the quadrangular wagon-wheel design always showed better performance than the triangular grain geometry. Quadrangular ports required a smaller initial thrust and a larger tank pressure, allowing for a lighter nozzle and a shorter engine, but also an heavier propellant tank. Quadrangular ports showed a slightly longer burning time, due to the lower thrust level during the whole engine burn. Moreover, they had a lower sliver, thanks to the presence of the central additional port, which resulted in a mass saving from around 10 kg to 30 kg (i.e., around 23% to 44% of the sliver mass), depending on the feeding system and the choice of N. Concerning the mixture ratio, the quadrangular ports exhibited a smaller shifting during engine operation, allowing for operations close to the value that provide the maximum specific impulse (about 6.5). In this way, a further mass saving due to propellant mass reduction, was obtained, which resulted in a payload improvement ranging from about 20 kg (N = 6, blow down) to more than 40 kg (N = 8, regulated). Optimal designs and performance are compared in Tables 3 and 4. Table 3. Optimal design and performance for different port shapes (PS): triangular ports (TP) and quadrangular ports (QP). The second column reports the feed system option (FO): blow down (BD) and regulated (R). The design variables are F i , α i , (p t ) i and E. R g , w, R t and L stand for grain outer radius, web thickness, throat radius and engine length, respectively. F/m is the maximum longitudinal acceleration imparted to the payload µ.  Table 4. Mass budget for different port shapes (PS): triangular ports (TP) and quadrangular ports (QP). The second column reports the feed system option (FO): blow down (BD) and regulated (R). The average mixture ratio and specific impulse are reported as α avg and (I SP ) avg , respectively. m res is the residual propellant mass, m g is the pressurizing gas mass, m t is the tank mass, m gt is the pressurizing gas tank mass, m nz is the nozzle mass, m cc is the combustion chamber mass and m hc is the engine casing mass.

Engine Modeling
Engine design is defined by a limited number of design parameters, and a proper selection of the design variables may help the design process. A convenient set comprises the initial thrust F i , the initial mixture ratio α i , the initial throat-to-port area J and nozzle expansion ratio E. From these variables the initial propellant mass flow (ṁ P ) i can be computed: where the initial values of the characteristic velocity c * i and thrust coefficient (C F ) i are calculated as described in Appendix D. The hydraulic resistance Z can then be computed by means of Equation (A35). The initial throat area (A th ) i is given by Equation (4), once the initial chamber pressure (p c ) i is assigned as will be discussed later.
The initial throat-to-port area ratio J gives the initial port area (A p ) i = (A th ) i /J and the mass conservation principle provides the initial burning area as described by Equation (5).
Once the required (A b ) i is know, one can determine grain geometry directly, when a single circular port is used, or iteratively, when more complex geometries are employed.
After designing the engine, the relations recalled in the Appendixes can be used to define the HRE behavior during operations. When throat erosion effects can be neglected, the throat area remains constant during operations, i.e., A th = (A th ) i = constant. Appendix D shows how throat erosion can be taken into account, when required. At each instant t, known the feeding pressure and the engine geometry, one can calculate fuel and oxidizer mass flow rates, their ratio α, chamber pressure p c and head-end pressure p 1 by means of the equations described in Appendix D. The thrust level can then be evaluated as F = p c A th C F once the value of C F has been computed at the actual altitude by means of Equation (A32). In the deterministic approach (see next Section) the time integration is carried out till contemporary burn out of solid fuel and depletion of liquid oxidizer occurs: the overall propellant mass (exhausted and sliver) can be evaluated, which allows for an estimation of the structural masses of the HRE. On the contrary, in the robust case, additional design parameters are used to define residual fuel and residual oxidizer.

Optimization
The optimization process can be regarded as the seek for a particular set of design variables or parameters which, at the same time, satisfy imposed constraints and maximize (or minimize) some assigned performance index. In the particular case of propulsion system, the optimization has usually to do with mass. For example the problem can be posed in the form of the maximization of payload mass for given initial mass or the minimization of engine dry mass for given mission performance. Obviously, the specific performance index, complexity of the problem, design issues and other aspects are strongly related to the considered applications and test cases.
HREs typically are characterized by variable thrust magnitude, mixture ratio and overall performance during operation due to evolving grain geometry and non-constant oxidizer flow rate. Moreover, the optimization problem can be even more tricky, for example when upper stages are considered, since the performance of the whole launcher have to be taken into account, and the optimization of the hybrid propulsion system cannot be detached from the ascent trajectory optimization. For these reason, the authors developed and refined over the years a coupled optimization procedure of engine design and trajectory.
Once grain geometry, propellant combination and feed system have been chosen, HRE performance at ignition and during operation are determined by a small set of parameters; the most commonly employed are: initial thrust magnitude F i , initial mixture ratio α i , initial nozzle-expansion ratio E i , initial chamber pressure value (p c ) i , initial feeding pressure value (p f s ) i , initial port to initial throat area J. Some of these latter are free and optimized by some sort of algorithm, which will be described in the followings, whereas others are constrained by safety or operation limits. For instance, the authors typically assign the initial chamber pressure by imposing a certain ratio between initial chamber pressure and initial feeding pressure, in order to avoid the coupling of oxidizer feed system and combustion chamber dynamics. In particular, (p c ) i = 0.4(p f s ) i is sufficient to guarantee p f s /p c ≥ 1.5 for the whole engine burn, despite the drop in p f s which can eventually occur when blow down or partially regulated system are employed. Some other variables may be implicitly given by imposed constraints related to the trajectory, such as maximum acceleration level or minimum regression rate. In this case, the relevant variables are dropped from the engine optimization process and added to the variables of the trajectory optimization, in order to grant the fulfillment of the related constraints. As a result, usually less than five design parameters have to be optimized, allowing for the use of efficient procedures. On the other hand, once engine design is assigned, an indirect method is employed for the trajectory in order to compute the optimization merit function.
The authors mainly employed two approaches to deal with engine design optimization: a direct method and evolutionary algorithms. The direct method involves the computation of the numeric derivatives of the performance index with respect to each design parameter, as prescribed by a forward-finite-difference scheme. Then, an iterative procedure based on Newton's method is used to obtain corrections of the tentative values of the design parameters in order to nullify the partial derivatives of the performance index [32]. The evolutionary algorithms, on the other hand, select a random initial population of individuals (i.e., a group of tentative engine designs) and evaluate their performance index. Then, biological-inspired processes are applied to the starting population which evolves towards improved individuals leading the optimization to the optimal design. Details about several evolutionary algorithms developed by the authors can be found in References [33][34][35].
A direct or evolutionary algorithm drives the choice of the set of the design parameters which gives the engine configuration and performance. For each engine design, an indirect method optimizes the trajectory. The authors take advantage of an indirect method developed at Politecnico di Torino in the early 1990s. The details of such procedure have been presented in detail in Refs. [4,36] and are here only summarized. A point mass rocket is considered. The state equations, written in non-dimensional form, in order to improve the accuracy of the numerical integration, provide the derivative of position r (radius, latitude and longitude), velocity v (radial, eastward and northward components) and rocket mass M, as reported in a vectorial form by Equation (6): where F and D are the non-dimensional thrust and drag vector, respectively. Equation (7) presents the inverse-square gravity field assumed for g, where G is the gravitational constant and M * is the planet mass.
The trajectory is split into phases, which depend on the particular problem considered. Examples of phases are vertical ascent, kick phase, zero-lift gravity turn ascent, coast arcs etc. Some time length are usually assigned (e.g., stage separation, fairing jettisoning, thrust rotation) whereas the others are determined during the trajectory optimization (e.g., engine burns and coast arcs). Thrust magnitude and propellant flow rate as functions of time are provided by the selected engine design. The mission to be performed is defined by a set of boundary conditions, such as initial or payload mass (or both), initial position, initial velocity, flight direction, maximum heat flux and target orbit in terms of altitude, eccentricity and inclination.
At this point, an adjoint variable is associated to each equation and the theory of optimal control [37] is able to provide Euler-Lagrange equations, alongside with algebraic equations, which give the control variables (i.e., the thrust direction) and with boundary conditions for optimality, which implicitly define the engine switching times. The multipoint boundary value problem that arises is solved by means of a numerical procedure based on Newton's method. A tentative value is initially assigned to each problem unknown and progressively modified aiming at the fulfillment of the problem boundary conditions. Such unknowns are usually the time length of some trajectory phases, the initial value of the adjoint variables and/or some design parameters that cannot be assigned "a priori" but only checked "a posteriori" after the optimal trajectory has been computed. Examples of the settings of the indirect trajectory optimization methods can be found in the references cited in Section 4. Alternative methods for trajectory optimization have been proposed by the authors in Ref. [38], employing an approximate control law for the thrust angles, and in Ref. [18], using a purely evolutionary algorithm.

Robustness
HRE design and optimization require a complex multidisciplinary optimization approach due to their one-lever control feature and their unique combustion mechanism. Furthermore, mission goal attainment could be jeopardized by severe deviation of the engine performance from nominal ones due to the presence of even small uncertainties in particularly sensible design aspects. In this context, a robust-based design and optimization approach is mandatory in order to reduce the sensitivity of the engine performance to such uncertainties. The definitions of robust design are various, here we summarize the basic concept of "robustness" as "the capability of the system to grant a fixed level of performance" (i.e., to match mission goals), "minimizing the effect of uncertainties in the design parameters without eliminating their causes" [39]. The computational cost of robustbased procedures is consistent because its growth is exponential depending on the number of uncertain parameters which are present in the specific problem. Thus, a preliminary selection of the uncertain quantities which exhibit the largest negative influences on the system performance should be carried out.
In Ref. [40] the authors proposed a sensitivity analysis procedure, based on Morris method [41]. First, the model parameters have to be split between design and input parameters. Design parameters, contained in vector b, are the quantities that can be freely set by the designer and their relative design space is discretized in a certain number of finite levels. Then a pool of starting design points are selected, e.g., by means of random search. An acceptance criterion should be included, in order to evaluate the sensitivity of the system to uncertainties for feasible (and interesting) solutions. For HRE design and optimization, possible criteria are a minimum insertion altitude, minimum payload mass or other performance index limitations. On the contrary, input parameters are all the other problem parameters, contained in vector X that, at this point of the procedure, are all regarded as uncertain and discretized in three levels (higher than nominal, nominal and lower than nominal values). In the specific case of HRE, examples of input parameters are nozzle throat erosion reference values, regression rate correlation coefficient and exponent, structural materials densities, propellant densities, initial tank pressure, trajectory parameters and efficiency related to characteristic velocity and thrust coefficient. The sensitivity analysis prescribes the computation of incremental ratios, called Elementary Effects (EEs), for each of the selected starting point. Once at a time, a variation of each input parameter is imposed, while all the others are fixed, in order to obtain a so called Morris trajectory inside the input parameters space. The process is repeated many times for each starting point. The direction of such variations and their sequence are both random. An EEs is obtained from each trajectory for each input parameter. Morris indices are then calculated as the mean of the absolute value and the standard deviation of these EEs. The mean value measures the influence of an input parameter on model output (i.e., on system performance), whereas the deviation measures its non-linear effects or interactions with other input parameters. Based on their indices, the input parameters, which have a strong effect on model output can be identified, and regarded as uncertain parameters in the subsequent robust-based optimization. On the other hand, input parameters which show a negligible effect on system performance are regarded as constant and their value is fixed on the nominal one.
In the literature, robust optimization problems are formally expressed as reported by Equation (8), where z p is the noise vector of the design parameters vector p, g j is the j-th inequality constraint, b L and b U are, in the general case, the lower and upper boundary of the design parameters, respectively.
In the present approach, the authors take into account uncertainties by means of three different levels for each uncertain variables, analogously to the discretization employed in the sensitivity analysis. The performance index, which depends on both design and uncertain parameters, is evaluated for each combination of the uncertain parameters prescribed by a suitable design of the experiments techniques keeping, on the contrary, the design parameters fixed. The corresponding optimization merit function Φ can be computed in several ways, e.g., as the average of the performance index for the prescribed uncertainties combinations or considering the worst-case scenario. An examples of design of experiments techniques is reported in Table 5 for six uncertain parameters. One can notice that the rows of the example table is equal to eighteen when the number of uncertain parameters (i.e., the columns) are six. A precise relation between the number of parameters, the number of levels and the number of combinations does not exist and strongly depends on the type of design of experiment technique considered. However, the computational effort of a robust-based optimization is many times the one required by the corresponding deterministic optimization, because no matter what technique is employed, the number of combinations required to obtain one merit function value is far greater than one (which can be regarded as the number of uncertain combinations in the deterministic optimization). Further details about sensitivity analysis and robust-based approach are provided in Appendix E. Table 5. Example of design of experiment techniques: Taguchi's L 18 orthogonal array. Each row reports a combination of uncertain parameters to be evaluated. The 1, 2 and 3 are the discrete levels that represent respectively lower than nominal, nominal and higher than nominal values. Each column corresponds to one of the uncertain parameters and shows its value expressed by the aforementioned levels. 3  3  3  3  3  2  1  1  2  3  3  2  2  2  3  1  1  2  3  3  1  2  2  3  1  2  3  2  3  3  2  3  1  3  1  3  3  1  2  1  2  1  1  3  2  2  1  1  2  1  3  3  2  1  3  2  1  1  3  2  1  2  1  3  2  2  2  3  2  1  3  2  3  1  3  2  1  3  1  3  3  1  2  3  2  1  1  2  3  3  3  2  2 3 1

Sounding Rocket for Micro-Gravity Platform
One of the key features of HREs is their relatively low cost when compared to the achievable performance, making a hybrid powered sounding rocket a viable option to transport small payloads to high altitudes, where the absence of the atmosphere allows for micro-gravity conditions that can be experienced in free flight. In such kind of application, a suitable optimization merit function to be maximized is the time spent above 100 km (later called t µg ), when the initial and payload mass, m i and µ respectively, are given. In this context, the authors considered a 500-kg sounding rocket as a case study [16]. The initial mass includes a 100-kg payload, fixed masses, propulsion system and propellants. The optimal mass split between propellants and propulsion, the optimal grain geometry and the corresponding optimal trajectory are sought in order to maximize t µg . A single-port grain design and blow down feed system are adopted to reduce the complexity and the cost of the propulsion system. Several options are taken into consideration as propellant combination, namely HP/PE, LOX/HTPB and NOX/HTPB. The design variables are three when self pressurizing NOX is employed as oxidizer: the initial thrust level F i , the initial mixture ratio α i and the nozzle expansion area ratio E. Erosion effects are here neglected, thus E = constant. On the other hand, a total of five design parameters are required when HP and LOX are used: the initial tank pressure (p t ) i and the ullage volume V u alongside with the aforementioned F i , α i and E. The use of NOX fixes the initial tank pressure due to the liquid-vapor equilibrium (see Appendix A.1 for the details), whereas the ullage volume is fixed at 3% of the tank volume.
The results are reported in Table 6, showing the higher performance of the HP/PE combination. The LOX solution exhibits a low initial mixture ratio, requiring a larger and heavier grain which increases propulsion system mass at the expense of the propellant mass on board. On the contrary, the use of the NOX combination results in worst performance due to its lower specific impulse, albeit interesting features in terms of overall rocket length and diameter are obtained thanks to the larger mixture ratio. Moreover, looking at Figure 4, one can notice that the optimal initial thrust and thrust history are quite different depending on the propellants employed. Table 6. Optimal design and performance for the considered propellant combinations. m p , R t , D and L are the propellants mass boarded, the nozzle throat radius, the engine overall diameter and length, respectively. (V g ) i and (p t ) i are optimization variables only for HP/PE and LOX/HTPB combinations.  The complex influence of the design parameters on the thrust history (i.e., on the optimal trajectory) is well exemplified switching the initial thrust of HP/PE and NOX/HTPB. In these cases, t µg presents a relevant reduction to 150 s for NOX/HTPB (−15%) and to 267 s for HP/PE (−10%). Concerning the tank pressure histories, the NOX solution is able to maintain an high pressure level longer than HP and LOX, due to its self-pressurization capability, resulting in a more regular thrust level. Such behavior holds until the liquid phase exist inside the oxidizer tank (t < 50 s); later the thrust drops and one can observe a change in tank pressure derivative due to the residual gas phase which is fed into the combustion chamber acting as a propellant.

Hypersonic Accelerator
Hypersonic propulsion and the development of hypersonic-capable vehicles is a growing trend in recent years. The test conditions of such systems, subsystems and devices require Mach number greater than five at an altitude in the range 25-50 km. Thus, one more time, a HRE is a viable option for this kind of missions for the same reasons explained in Section 4.1. Ref. [42] considered the following test case: initial mass m i fixed and equal to 1000 kg, payload mass µ equal to either 100 kg or 200 kg and 30-km target altitude. HP/PE, single circular port and blow down are the considered propellants, grain design and feed system, respectively. The optimization aims at the maximization of the Mach number M f at the target altitude. Despite the high altitude operation considered as target, the rocket spends a relevant amount of time at lower altitude, thus requiring an accurate modeling of its aerodynamic interactions (i.e., drag) in order to obtain significant results. A single stage design has been considered, which employs a single HRE, alongside with two different two stage designs: the first uses a single HRE in the first and second stage (1+1 configuration), whereas the second design employs two HRE in the first stage and a single engine in the second stage (2+1 configuration). To reduce development cost and system complexity, the same HRE is used for each multi-engine configuration. In this context, the design parameters are the initial thrust F i , the initial mixture ratio α i , the initial tank pressure (p t ) i , the initial volume of pressurizing gas (V g ) i and the nozzle expansion area ratio E. The value of (V g ) i is evaluated in order to grant a sufficient regression rate at burnout (i.e., avoid grain cooking), despite the blow down operation. Table 7 reports the optimal design and performance of the HRE. Table 7. Optimal design and performance of the hybrid-powered accelerator. D and L are the engine diameter and length. In the first column (Cfg.) the configurations are reported as: (1) single stage, (1+1) two stage/single engine per stage, (2+1) two stage/two engines in the first stage.

Cfg.
µ To reduce propulsion system mass, small thrust-low tank pressure designs should be preferred; however such optimization strategy fails in the 100-kg single stage case and in the 2+1 case (both payload) because the resulting thrust level tends to become too small and unfeasible solutions with negative coast arcs would be obtained. Thus, the initial thrust is dropped from the set of design parameters for such cases and determined by the trajectory optimization procedure itself, allowing for a feasible solution. The results show that a hybrid powered rocket can effectively accelerate a payload in the range 100-200 kg to the hypersonic corridor. Looking at the last column of Table 7, one can notice that, counter intuitively, the two stage design does not result in Mach number M f higher than the single stage configuration (first two rows), with the exception of the 2+1 configuration when the lowest payload is taken into account (second-to-last row).

Mars Ascent Vehicle
HREs are particularly well suited for missions which require shut down and restart capabilities and are less sensitive to low temperature, with respect to liquid rocket engines. These characteristics make a hybrid powered Mars Ascent Vehicle (MAV) a viable alternative to the heritage options for such application. In particular, the use of LOX/Wax propellant combination is really promising, since the oxidizer could be produced in situ from CO2 in Mars atmosphere. In Ref. [43] the authors analyzed a two-stage vehicle for the return of samples from Mars' surface and for a manned mission. A cluster of engines is used in the first stage, whereas a single HRE powers the second stage. The same engine is used in both stages, in order to limit the development costs of the MAV. Both a simple blow down feed system and a more complex regulated system are taken into account as possible architectures. The number of engine design parameters is four, three of them being in common between the two feed system alternatives, namely initial mixture ratio α i , initial tank pressure (p t ) i and nozzle area ratio E. The fourth parameter is the initial volume of pressurizing gas (V g ) i , when the blow down feed system is employed, or the exhausted oxidizer mass during regulated operation m * O , when the regulated system is considered. The number of engines used in the first stage must also be selected, since different configuration are possible, in particular N 1 = 2, 3, 4 are analyzed. The optimization aims at the maximization of the MAV payload µ at given initial mass, namely 500 kg in the sample return mission and 60 tons in the manned mission. A reference target orbit of 170-km altitude is taken into consideration for both the test cases. The most significant results are reported in Tables 8 and 9 for the sample return and manned mission, respectively. Table 8. Hybrid engine design and performance: sample return mission. The first column (Cfg.) reports the launcher configuration, e.g., 2+1 stands for two engines in the first stage and one engine in the second stage. The second column reports the feed system, blow down (BD) or regulated (R). V He is the volume of pressurizing gas D and L are engine diameter and length. m d is the engine dry mass and m p is the propellant mass.  Table 9. Hybrid engine design and performance: manned mission. The first column (Cfg.) reports the launcher configuration, e.g., 2+1 stands for two engines in the first stage and one engine in the second stage. The second column reports the feed system, blow down (BD) or regulated (R). D and L are engine diameter and length. m d is the engine dry mass and m p is the propellant mass. Concerning the sample return mission, the maximum payloads are obtained with the use of two and three engines, when the simple blow down feed or the regulated system are respectively used. The use of a regulated system always improves the resulting payload mass. On the other hand, in the manned mission, the greatest payload is achieved when four engines are used in the first stage since the larger and sustained thrust allows for a more efficient trajectory and, alongside with the dry mass reduction, increases the payload mass. For the same reasons, the payloads of the regulated optimal solutions are always greater than the corresponding blow down ones, analogously to the sample return mission optimal solutions.

Cfg. FO
These results are comparable, in terms of payload, to those provided by traditional liquid rocket eninges envisaged for a Mars return mission. However, the low cost and simplicity of HREs makes the described launcher options worth of consideration.

Small Launcher
Since the early 2000s, the demand for small launch vehicles is constantly growing due to the spread of the small satellite market. Until recently, operative launcher especially designed for such application are an exception and HREs, due to their features, are a promising propulsion system for a small satellite launcher. For these reasons, in Refs. [44][45][46] the authors investigated the optimal design of a hybrid-powered 5000-kg small launcher. LOX/paraffin-based wax is considered as propellant combination for its high regression rate, which allows for a simpler single port grain design. Concerning the feed system, a blow down type is taken into account. The launcher employs the same HRE in different number in each stage: six, three and one in the first, second and third stage, respectively. An airborne launch is considered at given altitude and speed (e.g., Pegasus Launcher).
The optimization is carried out from both a deterministic and a robust point of view, and the number of design parameters are respectively two ((V g ) i and α i ) and five ((V g ) i , α i , E, total oxidizer mass boarded m O,tot and R g ). Different trajectory setups are taken into consideration, namely given initial thrust (called design A) and free initial thrust/fixed maximum acceleration (called design B). In the latter cases F i is an additional design parameters which is optimized. Only the regression rate is regarded as uncertain in the robust-based optimization. The optimization aims at the maximization of the payload inserted in a given-altitude, circular and polar orbit. The optimal designs for the deterministic optimization are reported in Table 10, whereas the time histories of their key performance indexes are shown in Figure 5. Table 10. Performance and design summary: deterministic optimization (800 km altitude). R f stands for the final grain burning distance, m F is the total fuel mass burned and m d is the engine dry mass.

Des.
F In this application, the optimal ascent trajectory prescribes a large acceleration for a suitably long time, in order to reduce gravitation and aerodynamic losses. Thus, the deterministic design A exhibits a large value of (V g ) i because its initial thrust is fixed, whereas the deterministic design B is free to follow a steep trajectory by means of a larger initial thrust F i despite a smaller volume of pressurizing gas. This behavior allows for an engine dry mass reduction, albeit a greater propellant load is required due to the less efficient trajectory and lower nozzle expansion area ratio E (i.e., lower specific impulse I SP ), with respect to the deterministic design A. Thus, the deterministic design B is capable of inserting a 73.6-kg payload into the target orbit, i.e., a 25-kg increase with respect to the deterministic design A.
In the context of the robust-based optimization, the authors decided to maintain the same design philosophy, first seeking for a fixed initial thrust solution (here called robust design A) and then trying to find out an optimized thrust solution (here called robust design B). The hope was to obtain a performance improvement, analogously at what happened in the deterministic optimization. In the context of the robust-based optimization, the performance index to be maximized was Φ avg , a linear combination of the payload µ inserted in the target orbit and the average altitude violation ∆ avg due to uncertainty with respect to such orbit. A large value is here used for the linear combination coefficient in order to force the altitude violation to zero. A detailed definition of the robust-based merit function can be found in Ref. [46]. The number of design parameters grows up to four when robust-based procedure is employed, namely the initial mixture ratio α i , the grain outer radius R g , the initial pressurizing gas volume (V g ) i and the total oxidizer mass boarded. The optimal robust design A is reported in Table 11 and compared to the optimal deterministic design A. One can observe that ∆ avg is equal to 500 km for the deterministic solution in the face of a 800-km target orbit, i.e., more than half of the considered uncertain parameters combinations lead to a rocket crash (insertion altitude null) when a proper robust-based design is not carried out. On the contrary, the ROB A solution exhibits a null altitude violation, granting robustness in the design, at the cost of a relevant payload reduction (−28% with respect to deterministic design A). To improve this result, a robust design B (ROB B ) is sough, by freeing the initial thrust, which become an additional trajectory unknown, and adding a constraint on the final acceleration, as previously done in the deterministic procedure. This approach results in payload increases, with respect to ROB A optimal solution, of 10 kg and 20 kg, when maximum accelerations of 6 g and 5 g are imposed in the indirect trajectory optimization. However, the average altitude violation of these solution is not passably close to zero (131 km and 50 km when 6 g and 5 g are respectively used) and unacceptable in a robust-based optimization.

Upper Stage
In the context of small and low-cost launchers, SRMs have been used extensively used: one can consider the Pegasus launcher and the European Vega launcher as an example for two different payload class. In Pegasus all the stages are powered by SRMs. However, the lack of shut down/restart capability of SRMs imposes the use of at least one stage for each burn required by the ascent trajectory, whereas the absence of throttling determines relevant scattering in the target orbit insertion parameters (e.g., the altitude). Moreover, SRMs are characterized by a lower specific impulse than their liquid counterparts, which results in penalties in the final part of the ascent. These issues are solved in Vega by employing an additional small liquid-powered stage, which is devoted to the final orbital maneuvers but results in higher costs and worse structural efficiency. HREs represents an alternative and viable solution for an upper stage in this context, because they are able to deliver the high performance of a SRM and the fine control features of the aforementioned additional small LRE, but at a lower cost. Furthermore, the use of typical hybrid propellant combinations avoids safety issues related to the chemicals used in the LRE final stage, such as Nitrogen TetrOxide (NTO) and MonoMethylHydrazine (MMH).
For these reasons, the authors focused on the design and optimization of a hybrid powered upper stage suitable for the replacement of the third and fourth stages of a launcher based on the characteristic of the European Vega Launch Vehicle (VLV). The optimization aims at the maximization of the payload inserted into a 700-km polar orbit with a launch from Kourou. A baseline value for the payload mass using three solid-propellant stage and a liquid-powered upper stage is 1430 kg, consistent with the performance of the Vega [47]. The initial mass of the launcher is given and the lower stages are unmodified.    Back in 2010 the authors approached to this topic for the first time [48]. The trajectory is optimized from lift-off to orbit insertion and the same constraints of the reference launch are taken into account, such as heat flux limitation after fairing jettisoning. The considered propellant combination is HP/PE, which requires a multiple-port grain design in order to deliver the required high thrust levels. Several options for grain geometry and feed system are compared. The engine design parameters are four, namely the initial thrust F i , the initial mixture ratio α i , the initial tank pressure (p t ) i and the nozzle area ratio E (erosion effect are neglected). The initial ullage volume (V g ) i (blow down case only), the oxidizer mass at the start of the blow down and its final value (regulated case only), and the grain outer radius R g are additional unknowns of the ascent trajectory optimization. Figure 6 depicts tank pressure, thrust and mixture ratio histories for a blow down and a partially regulated feed system. These results are obtained assuming eight quadrilateral ports for the fuel grain. The use of the blow down feed system results in a large thrust level variation during engine burn. Thus, a really high initial tank pressure is mandatory, in order to maintain a sufficient thrust until the end of the engine burn, also avoiding an increase in pressurizing gas mass. A performance improvement is obtained employing a partially regulated feed system: the introduction of an auxiliary pressurizing gas tank is balanced, in terms of engine dry mass, by the reduction of the oxidizer tank mass due to the lower pressurization.
However, a stronger mixture ratio shifting penalizes the performance of the HRE during the constant pressure operation, with respect to the blow down optimal solution. Further performance improvement are achievable employing an electrically driven turbo-pump feed system [29]. The effect of grain geometry (number of ports and shapes) on launcher performance has been already described at the end of Section 2.3.
In 2014, the authors considered alternative propellant combinations for the same upper stage application, namely HP90%/HTPB, LOX/HTPB and LOX/paraffin-based wax [33]. Moreover, the authors took into account nozzle throat erosion effects by means of the Bartz method (see Appendix D). The use of paraffin-based wax allows for a single port grain design, thanks to its relatively high regression rate value. The optimal solution presents a payload increase of about 350 kg with respect to the best performing conventional propellant combination, as shown in Table 12, being at the same time only slightly affected by performance deterioration due to throat erosion.
In 2015 the authors started to shift their attention from design trade-off between different HRE configuration, regarded strictly from a deterministic point of view, to uncertainty based approaches in the context of the optimization of the aforementioned upper stage. The first attempt to perform a robust-based optimization is presented in Ref. [11]. Regression rate correlation coefficient a and exponent n are varied around their nominal values for the optimal solutions and off-design launcher performance are evaluated. The results show that the mission is unfeasible when the actual regression rate is greater than the nominal one, because a relevant amount of unburned oxidizer remains as an inert mass after fuel grain burn out. On the other hand, off-nominal launcher performance of the optimal design are satisfying when a regression rate larger than nominal is assumed during the optimization procedure. In this case the average altitude is increased and the standard deviation reduced (from 250 km to 12 km).
This management of the intrinsic uncertainties represents a preliminary approach to robust design, which the authors improved in Ref. [49] embedding uncertain parameters variations inside the optimization procedure, in place of an "a posteriori" check of optimal solutions off-design performance. Thus, a new optimization merit function is introduced, combining linearly the payload mass, which is unaffected by the uncertainties, with the average altitude violation with respect to target altitude, which, on the contrary, is affected by uncertainty because the actual insertion orbit depends on the actual value of the regression rate. The number of design parameters rises due to the robust-based formulation of the optimization problem. In the deterministic approach oxidizer and fuel masses are such that they end at the same moment at burnout and are optimized (i.e., minimized) by the indirect procedure to reach exactly the target orbit. On the other hand, in the robust-based optimization, the actual value of regression rate varies resulting in three possible situations: • the regression rate is higher than nominal; in this case all the fuel burns whereas a certain amount of unburned oxidizer remains in the tank; • the regression rate is lower than nominal; in this case all the oxidizer is exhausted whereas a fuel sliver is present in the combustion chamber; • the regression rate is the nominal one; in this case fuel and oxidizer end at the same time.
The surplus in fuel/oxidizer mass required to reach (at least) target orbit, even in the higher/lower than nominal operation, have to be fixed before trajectory optimization, in order to compute launcher payload and define its design. Thus, the fuel mass m F (given by grain outer radius R g only for circular port grain, R g and web thickness w for multiport grain design), and the total oxidizer mass boarded are no longer trajectory unknowns but become engine design parameters. These surplus in propellant masses result in payload penalties for the robust optimal solutions with respect to deterministic one. Obviously the robust designs are able to reach higher altitude than the target when nominal regression rate occurs. The results prove that robustness in the design can be achieved at the expense of a small payload reduction of 82 kg, i.e., less than 5% of the launcher payload, with respect to the deterministic optimal solution for the same grain geometry, propellants and feed system. Figure 7 reports the thrust and mixture ratio histories for the deterministic and optimal design. The robust design shows a slightly larger thrust magnitude variation and mixture ratio shifting during engine burn, due to the different grain geometry which causes enhanced changes of port and burning area. The obtained robust optimal solution is able to grant the fulfillment of launcher mission goals (i.e., 100% target orbit reaching) despite the presence of regression rate uncertainties. In Ref. [30] this robust-based approach is applied to a slightly different upper stage architecture, which employs electrically driven pump to feed the oxidizer into the combustion chamber. Mission goals and all the other aspects of the problem are unmodified. The results show that this feed system option is viable, even when uncertainties in the design are taken into account. The use of electric systems is able to grant a relevant performance boost with respect to conventional gas pressurized feed system, which in terms of payload provides a 12% increase (+250 kg). In the robust-based optimization new core aspects become more relevant, namely the number of uncertain parameters and the algorithm to be used to lead the engine design optimization. The first determines the number of uncertain parameters combination that have to be considered in the robust-based procedure, i.e., the number of trajectory integration required to compute one value of the optimization merit function. The latter is devoted to the choice of the engine designs to be evaluated during the optimization in order to reach an optimal solution. Thus, both this aspects have a huge impact on the computational effort of the whole robust-based optimization method. For these reasons, the authors in Ref. [35] proposed a comparison of the performance of different optimization algorithm, namely a Genetic Algorithm (GA), a Particle Swarm Optimization algorithm (PSO), a hybrid genetic-Taguchi algorithm and an Iterated Local Search (ILS). The optimization of the hybrid powered upper stage is the test case, one more time. PSO and ILS outperformed the other methods, showing similar performance, despite their completely different nature: the first is a population-based algorithm, whereas the second employs a meta-heuristic approach.
In Ref. [40] the authors described an approach, which takes advantage of a sensitivity analysis performed by means of Morris Method, to select the model parameters to be regarded as uncertain [41]. Details about this technique are reported in Section 3.1. The results show that even small variations of six model parameter (regression rate correlation coefficient a and exponent n, throat erosion synthetic parameter K ero , fuel density ρ F , initial tank pressure (p t ) i and initial throat to port area ratio J) have a relevant impact on system performance and thus have to be considered as uncertain in the upper stage test case. In terms of payload, this uncertainty model determines a 12% payload reduction, with respect to the deterministic optimal solution obtained for the same grain geometry and propellants. This reduction is more than twice the one imposed by the simpler uncertain model, which considered only the regression rate uncertainties, mentioned in the previous paragraph, but it is worth taking the robustness gain into account. For the sake of comparison, Table 13 reports the normalized mass ratios for the reference and the robust optimized solutions. Table 13. Normalized ratios for the reference (VLV) and robust optimized solutions. GP2 and EP2 are the solutions for gas pressurized and electric turbo pump feed system employing simpler 2-parameter robust-based approach. GP6 reports the mass ratios for the 6-parameter robust-based approach.

Conclusions
The present overview of the activities on HRE at Politecnico di Torino has been structured in three main parts. First, the preliminary design options, (i.e., propellant combination, feed system and grain design) have been presented and investigated by means of references to specific works which faced one of these fundamental choices at a time. Then, in the second part of the review, the focus of the writers shifted on the description of the mathematical model, which underlies the description of HREs performance and operation, and the developed optimization approaches, concerning engine and trajectory. In this context, special attention has been devoted to the robustness issue, which is the authors' most investigated topic in their recent works. In the last part of the manuscript, relevant examples of HREs applications, analyzed by the authors during the years, have been collected.
This rundown of examples underlined the complex interactions between problem input parameters and system performance. The optimal designs are the results of a tradeoff process between conflicting requirements. In general, high initial thrust levels and sustained thrust are beneficial for the trajectory, reducing gravitational and aerodynamic losses. On the other hand, such level of performance comes together with an increase in the propulsion system mass, because heavier and/or additional components are required, reducing the available payload mass. Moreover, this scenario is further complicated by both the coupling between thrust and mixture ratio, typical of hybrids, and the sensitivity of the system to the ineluctable presence of uncertainties in the design. For these reasons, the true advantages attainable by the use of HREs can be highlighted solely through the employment of proper tools, which are able to couple engine design and trajectory optimization, embedding uncertainty management inside the procedure. In this context, the procedure developed at the Politecnico di Torino and here presented, proved to be flexible and reliable thanks to a well-reasoned use of the most suitable direct/evolutionary/indirect methods.
Nevertheless, issues are still present in the proposed approaches. The seek of more robust solutions struggles against the huge computational costs due to the growth of the number of uncertain parameters and approximated methods may not be suitable, because the introduced simplifications could override the complex parameters interaction intrinsic in HREs design and operation. However, the use of innovative propellants and feed system, such as paraffin-based wax and electric turbo-pump, proved to be able to overcome hybrids performance reduction and methods limitations due to a rigorous robust-based approach.
In the near future, the authors will point their attention upon the generalization of their robust-based coupled engine/trajectory optimization method to even more diverse HREs applications and the improvements in the models for the engine description (dry masses and thermal control), paying special attention to the upcoming technological development in the field of hybrid propulsion systems. Acknowledgments: Computational resources provided by hpc@polito, which is a project of Academic Computing within the Department of Control and Computer Engineering at the Politecnico di Torino (http://www.hpc.polito.it, accessed on 11 August 2021).

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: In Ref. [28], among other oxidizer control strategies, the author analyzed the positive effect on performance of the repressurization, with respect to the uncontrolled blow down operation. In that test study the propellant combination was 85%HP/PE and the designed HRE was suitable as an in-space propulsion system for small satellites. In that context, the optimization was carried out aiming at the minimization of the propulsion system mass, which is the sum of engine dry mass, propellants masses and pressurizing gas mass. The initial mass of the satellite was given (100 kg) and the required ∆V assigned (160 m/s). In such design the repressurization gas was stored at pressure p a in a separate auxiliary tank and released at a given value of the oxidizer tank pressure p rep , restoring the initial value (p t ) i . Two additional design parameters, namely p a and p rep , were required, and the amount of auxiliary gas was evaluated by means of Equation (A2) under some simplifying assumptions (adiabatic process, pressure equilibrium between auxiliary and oxidizer tank at the end of the repressurization and no pressure rise due to heat exchange after the repressurization): where (V g ) rep stands for the ullage volume when p t = p rep and V a is the auxiliary gas tank volume. The effect of p a on engine geometry, mass and performance is practically negligible, thus it was considered as fixed (p a = 150 bar). On the other hand, p rep affects both engine design and operation and had to be constrained above a certain value to prevent solid fuel cooking. Results showed that repressurization is able to reduce thrust magnitude variation at the cost of a small overall propulsion mass and volume increases, due to the introduction of the auxiliary gas tank, whereas the mixture ratio shifting is comparable with the reference blow down operation.

Appendix A.2. Self-Pressurizing Oxidizer
In Ref. [26] the authors faced the problem of modeling the tank pressure history during the blow down of a self-pressurizing oxidizer. A simple homogeneous model was compared to a novel two-phase lumped model considering NOX. Both models neglected heat transfer between the tank walls and the oxidizer. The homogeneous model [50] assumes that both the liquid and vapor phase are at the saturation temperature of the propellant at the tank total pressure. The overall propellant mass (i.e., the sum of liquid and vapor inside the oxidizer tank) are considered in an homogeneous state, whose thermodynamic properties are mass-averaged. The first law of thermodynamics and the law of conservation of mass can then be applied to the control volume enclosed in the internal surface of the tank wall, obtaining the tank pressure derivative ∂p/∂t reported in Equation (A3).ṗ where p t , ρ t , V t , h t and u t are pressure, density, volume, specific enthalpy and specific internal energy of the tank, respectively, whereasṁ O and h O are oxidizer mass flow and specific enthalpy. In the homogeneous model, the liquid is supposed to be well mixed (without stratification) and exchanges energy with the vapor at an infinite rate. For these reasons, the corresponding pressure time evolution represents an upper limit for the propellant tank pressure p t .
During engine operation, the liquid is actually warmer than the vapor and a temperature gradient exists inside the oxidizer tank. Thus the authors considered a two-region (bulk liquid and bulk vapor) simplified lumped model, in order to overcome the assumptions of Appendix A. 3

. Regulated Feed System
In this case a relevant amount of pressurizing gas is stored in an auxiliary tank and the pressure in the oxidizer tank can be controlled by a continuous addition of pressurizing gas through a regulation valve. Different control strategies can be applied, such as constant tank pressure, constant chamber pressure, etc. The regulated operation may be then followed by a conventional blow down phase; the resulting operation is here called partially regulated. The authors compared the performance of blow down and partially regulated system in their works considering various applications: sounding rocket in Ref. [15], in-space propulsion system in Ref. [28], micro-gravity platform in Ref. [16] and upper stage in Refs. [31,48]. In general HRE performance resulted improved when a regulated system was employed with respect to the simpler blow down operation, despite the presence of an auxiliary tank and pressurizing gas. Nevertheless, specific characteristics of the optimal design are in general case-dependent, and thus will not be mentioned in this review.

Appendix B. Electrical Turbo Pump Feed System
The authors analyzed the possible use of an electric pump feed system for a hybridpowerd upper stage in Ref. [29] and later in Ref. [30]. The mass of the feed system components (i.e., pump, electrical systems, batteries) were evaluated by means of their typical power density values δ. Equation (A11) provided motor and pump mass m ep , where P e,max is the maximum electrical power required by the motor which drives the pump. In general, the electrical power required can be computed by means of Equation (A12) as a function of discharge pressure p d , oxidizer density ρ O and mass flowṁ O . Tank pressure p t has been assumed constant and equal to 1 bar in this context [51]. The amount of the required pressurizing gas to maintain this constant pressure was small and it was neglected.
The conversion effectiveness of electrical energy (stored in the battery packs) into flow head rise was accounted for by means of the overall efficiency η ep . The mass of the batteries was constrained by the most stringent requirement between the maximum electrical power P e,max and the total electrical energy E e,tot needed to drive the pump for the whole engine burn duration t burn . The total electrical energy was computed as reported in Equation (A13), whereas battery pack mass m b was evaluated by means of typical values of power density δ bp and energy density δ be . A safety factor of 1.2 has been assumed as reported in Equation (A14).
The pump was assumed to operate at constant power P e = P e,max , where P e,max is determined by the initial pump discharge pressure (p d ) i which is the optimization parameter related to electric feed systems. Thus, Equation (A13) can be easily integrated giving E e,tot = P e,max t burn .
Concerning the batteries mass, the power-based constrain can be evaluated immediately from the choice of (p d ) i , while energy-based constrain has to be computed "a posteriori", once the actual burning time is known. Thus, a characteristic burn time t * burn = δ be /δ bp can be defined, which represents the simultaneous fulfillment of both constraints. If the actual burn is shorter than t * burn , power-constrained mass is larger than the energy-constrained one and battery pack mass is independent of the ascent trajectory. On the other hand, if the engine burn is longer than the characteristic burn time, energy-constrained mass is larger than power-constrained one, and batteries mass must be checked "a posteriori".
To obtain an actual advantage, in terms of engine dry mass reduction with respect to a gas pressurized feed system, the electric components have to be as light as possible, despite the strict requirements on their energy and power levels. Their small size and light weight make Lithium batteries well suited for electric feed system. Moreover, the typical rocket operation requires high thrust (i.e., high electric power) for a relatively short time, making Lithium batteries even more attractive thanks to their high rate capabilities. In Table A1 the electric properties employed for the evaluation of the masses of the electric component are presented [52][53][54]. The first set of properties, namely set A, has been used by the authors in their older work on this topic, whereas a second set of properties, namely set B, has been added in their more recent work for the sake of comparison, in order to take into account the most recent technological development in the field of electric motor and batteries.

Appendix C. Multi-Port Fuel Grain Design
Two different multi-port grain geometries have been considered by the authors: triangular port grains and wagon-wheel grains. In triangular port grains, the design is fixed by three parameters, namely the number of ports N, the web thickness w and the grain outer radius R g . Equation (A15) gives the initial port area (A p ) i , where x, h and β are defined as reported in Equations (A16)-(A18) making reference to Figure A1. Figure A1. Triangular port geometry.
The burning area A b at a generic time t (which correspond to a burning distance 0 ≤ y(t) ≤ w) can be obtained using Equation (A19), whereas Equation (A20) gives the generic port area A p .
The second multi-port grain geometry employs a wagon-wheel grain, in which N quadrangular ports enclose a central circular port with initial radius R c , as shown in Figure A2. To minimize the sliver, the grain external radius has to be R g = R c + h 1 + 3w, where h 1 is the radial extension of the quadrangular ports.
It is convenient to define three dummy variables, namely x 1 , x 2 and h 2 as reported in Equations (A21)-(A23).
To have the same regression rate at engine ignition in all ports (both circular and quadrangular ones), the initial area of each quadrangular port is assumed to be equal to the initial central circular port area. In addition, the regression rate is considered to remain the same in all ports for the whole engine burn duration, despite different port area and burning surface history among circular and quadrangular ports. Thus, the overall initial port area (A p ) i can be computed by means of Equations (A24) and (A25): where (A c ) i is the initial port area of the central circular port and (A q ) i is the initial port area of one of the N quadrangular ports. In the case of wagon-wheel grain, the port area A q at a generic time t, which correspond to a burning distance 0 ≤ y(t) ≤ w, can be obtained using Equation (A26), whereas Equation (A27) gives the generic port area A c of the central circular port: The burning surface A b,q of each quadrangular port can be obtained using Equation (A28), whereas Equation (A29) gives the generic burning surface A b,c of the central circular port.
In the end, the overall port area A p and burning surface A b are expressed by means of Equation (A30) and (A31).

Appendix D. Engine Modeling
The set of equation used to model the HRE are here recalled. Characteristic velocity c * and the specific heat ratio γ for the given propellant combination are evaluated by means of polynomial fittings of the results provided by CEA [55] as function of α at reference values of chamber pressure.
It is worth noting that the actual chamber pressure can span over a wide range during engine operation. Nevertheless, the related performance errors is very small, making a constant pressure assumption acceptable. A c * -efficiency η c * is also included [12]. A frozen expansion is assumed in the evaluation of performance and gas properties, i.e., the gas combustion composition remains the combustion chamber one during expansion in the nozzle. The thrust coefficient C F is then evaluated by means of a one-dimensional isentropic expansion in the nozzle, with constant γ, as reported in Equation (A32), where E is the nozzle-expansion area ratio, p c is the chamber pressure, p e is nozzle exit pressure and p atm is the ambient pressure. A correction factor η F is usually introduced to take into account real effects.
The value of E can be either given or an additional optimization parameter, whereas the ambient pressure p atm is provided by an atmosphere model as altitude varies. The aforementioned set of assumptions gives conservative performance evaluation.
The propellants combinations selection also determines the regression rateẏ, by means of the classical correlationẏ where the coefficient a and n are empirical values (e.g., see Table 1) and G O =ṁ O /A p is the oxidizer mass flux. The fuel mass flowṁ F can be evaluated according to Equation (A34) where ρ F is the fuel density and A b the burning surface.
Equation (A35) gives the oxidizer flow rateṁ F , where p f s is the pressure level due to the selected feeding system (e.g., tank pressure p t or discharge pressure p d when gas pressurized or turbopump systems are respectively employed), p 1 is the head-end pressure and Z is the hydraulic resistance in the oxidizer flow path from the tank/the pump discharge to the combustion chamber.
Please note that the value of Z is evaluated on the basis of engine design parameters and it is considered to be constant during engine operation. Another option, not considered in the authors' works, could be the use of a cavitating venturi, which is able to deliver a mass flow rate independent of downstream pressure changes, when operated at cavitation conditions. The aforementioned head-end pressure p 1 can be related to chamber/nozzle stagnation pressure p c by means of the approximate relation reported in Equation (A36), where A th is the throat area [56].
Overall propellant mass flow isṁ and the mixture ratio α can be computed by means of Equation (A38) to evaluate c * and γ.
The chamber/nozzle stagnation pressure p c is calculated assuming an isentropic expansion.
Throat erosion effects can be neglected, assuming a constant value for the throat area during operation (i.e., A th = (A th ) i = constant), or taken into account. In the latter case, one can employ Bartz method [57], in order to model the influences of chamber pressure p c and throat radius R th on the rate of throat erosionṡ, as reported in Equation (A40). The reference parameters for pressure p c,re f , erosion rateṡ re f and throat radius R th,re f can be found in the literature [58].ṡ =ṡ re f p c p c,re f The integration of Equation (A40) during engine burn allows for the computation of the time evolution of throat area A th and nozzle expansion area ratio E, which are required in Equations (A32) and (A39). Throat erosion makes E(t) ≤ E i because A th (t) ≥ (A th ) i , lowering HRE performance. In Ref. [18], the authors compared HRE performance evaluated considering and neglecting throat erosion effects. The considered application was a microgravity platform and the optimization aimed at the maximization of the time spent above 100 km, i.e., in micro-gravity conditions. The impact of throat erosion on performance resulted in being strongly dependent on the propellant combination employed. Self pressurizing propellant (see Sections 2.1 and 2.2) showed a major reduction in performance, whereas combinations which employed HP or LOX experienced a minor performance decrease (less than 4%). It is worth noting that for the LOX/Wax combination, which is particularly fascinating for its high regression rate, the performance reduction due to throat erosion effects was minor than 1% thanks to the lower pressure and larger throat area required.
Finally, thrust F can be evaluated as

Appendix E. Robust-Based Method
Appendix E.1. Sensitivity Analysis and Screening The proper selection of the uncertain parameters is crucial in robust-based models, because their number greatly affects the computational effort required by the resulting model. In this context, Sensitivity Analysis (SA) and parameter screening are useful tools in the building process of a robust-based model. In Ref. [40] the authors considered the design and optimization of a hybrid powered upper stage and employed Morris method for SA [41]. The present Appendix summarizes the most important steps of that procedure.
First of all, model parameters are split between design and input parameters. Design parameters, collected in vector b, can be freely set by the designer, and their design space may be discretized in levels. A random search is employed in order to select one hundred starting design points and an 100 6 grid is used for the design parameters. An altitude-based acceptance criterion (h ≥ 500 km) is taken into account, granting that the pool of starting design points for the SA is actually related to the optimization aim, which is the insertion of the payload into a 700-km circular polar orbit. In this way, design points, which lead to lower insertion altitude or rocket crashes, are not taken into account in the SA.
On the contrary, input parameters, collected in vector X and reported in the second column of Table A2, are regarded as uncertain and discretized by means of three levels (higher than nominal, nominal and lower than nominal). Then, the evaluation of incremental ratios, called elementary effects (EEs), for each starting point is prescribed by Morris method.
where ∆ is the relative input parameter variation and e j a vector of the canonical base. In the considered application, two objectives are relevant (i.e., payload mass µ and altitude h).
For this reason, an -constraint approach is adopted, aiming at finding the Pareto front of solutions: [59] an altitude constraint violation ∆ i = max(0, h * − h i ) is computed and then the merit function of the i-th run is defined by means of Equation (A43).
where k = 20 kg/km is selected to force the average altitude at . Only the case = h * = 700 km is here considered. At this point, Morris indices can be obtained as the mean of the absolute value µ * j and the standard deviation σ j of the elementary effects. µ * j estimates the influence of the j-th input parameter on the system output, whereas σ j estimates its non-linear effects or interactions. Model input vector X and its nominal/off nominal values are summarized in Table A2 for the considered test case, whereas Table A3 reports the settings of the sensitivity analysis and the relative computational cost in terms of objective function evaluations [13,58,60]. Table A3. SA settings and computational cost.

Parameter Value
Number of steps for design parameters 100 Number of input parameters 19 Number of starting points 100 Number of repetitions 10 Number of function evaluations 100 × (19 + 1) × 10 = 20,000 Figure A3 depicts SA mean outputs for each input parameters (i.e., averaged among all starting points and repetitions ) µ * j and σ j . Table A4 summarizes the same results for the sake of completeness. Input parameters may be split into three groups depending on their behavior:

1.
parameters with negligible or no effect on model output (i.e., both µ * j and σ j are negligible or zero): s cc , ||v i || ||v i,re f || , ||r i || ||r ⊕ || , η C F and η c * (marked in red in Figure A3). Hence, these input parameters may be regarded as model constants in the robust-based design and optimization and their values will be set equal to nominal ones; 2.
parameters with linear effect on model output (i.e., µ * j and σ j are both significant and µ * j σ j ): ρ Al , ρ abl , ρ O , T g , s lin , V ull (marked in black in Figure A3). The effect of these input parameters is expected to be the same everywhere in the design space. Thus, they may be considered as model constant during the robust-based design and optimization, in order to minimize the number of uncertain parameters. However, their actual effect on the optimal solutions found out should be checked "a posteriori"; 3.
parameters with non-linear effect on model output (i.e., both µ * j and σ j are significant): s re f , p c,re f , R th,re f , a, n, ρ F , p t,i and J (given in light blue and green in Figure A3). These input parameters exhibit strong, not uniform and unpredictable effects on system output among the design space. Hence, such sources of uncertainty cannot be neglected and have to be taken into account in the robust-based design and optimization in each objective function evaluation. One can notice that the parameters marked in light blue (ṡ re f , p c,re f and R th,re f ) exhibit smaller (but still significant and not linear) effects with respect to the parameters marked in green (a, n, ρ F , p t,i and J). Since these three parameters are involved only in the nozzle throat erosion model (see Equation (A40)), the authors opt for the definition of a synthetic erosion uncertain input parameters K ero , as reported in Equation (A44), in order to maintain the number of uncertain parameters as low as possible.
K ero =ṡ re f · R 0.2 th,re f /p 0.8 c,re f (A44) In this way, the rate of throat erosionṡ, can be rewritten as: Thus, in the end, uncertain parameters vector p collects six parameters and can be written as: p = [K ero , a, n, ρ F , p t,i , J] (A46) In conclusion, the proposed SA resulted able to reduce the number of uncertain parameters from nineteen down to six, reducing by far the number of experimental runs which will be required in the subsequent robust-based design and optimization procedure.

. Robust Design Model
In the robust-based design and optimization, uncertainties are taken into account by the authors assigning three different levels for each uncertain variables, analogously to what did in the aforementioned SA. The performance index (i.e., altitude of the attained orbit h i and payload µ i ) are evaluated for each i-th combination of uncertain parameters, as prescribed by a suitable design of the experiments (DOE) techniques. In this context, only average performance is considered and therefore the robust-based optimization merit function can be computed as: where P i = ∑ j p j p tot , p j = 0.25 if j-th level is equal to 1 or 3 (i.e., an off-nominal value is considered for the j-th uncertain parameters) and p j = 0.5 if j-th level is equal to 2 (i.e., the nominal value is considered for the j-th uncertain parameters) and p tot = ∑ ij p j for j = 1, . . . , 6 (thus ∑ i P i = 1). Only the case h * = 700 km (that is the most demanding in terms of robustness for the considered test case) is here taken into account. Then, one can evaluate the average payload mass µ avg and average violation ∆ avg as: Figure A3. SA outputs: σ j vs. µ * j (log-log plot).
Once the number of uncertain parameters has been reduced by the SA, a suitable DOE technique may be used to guide efficiently the choice of the "experimental" runs to be performed, in order to take into account uncertainties in an affordable and reasonably accurate way. DOE techniques express uncertain parameters combination by means of arrays of discrete levels, namely 1, 2 and 3, corresponding to lower than nominal, nominal and greater than nominal values for each uncertain parameters (see Table A2). The most common DOE technique is the full factorial design which needs to execute "experimental" runs with all possible combinations of parameters levels. In the considered application, the number of uncertain parameters is equal to six, whereas three levels are employed for the discretization. Thus, a full factorial design would require 3 6 = 729 experimental runs. Such an high number of merit function evaluations is not computationally affordable to compute a single values of Φ avg during the HRE optimization procedure. For this reason, the authors propose the use of three different fractional factorial design (FFD), which consider only a subset of the full factorial design, requiring far less runs [61]. The first design, here called F 12 and reported in Table A5, employs a really simple array, in which each experimental run takes into account only one off-nominal value for a single parameter. The second FFD technique exploits Taguchi's L 18 orthogonal array [39] which is shown in Table A6. The latter technique takes advantage of Box-Behnken's array [62], which is shown in Table A7.