Comparison of Optimization Techniques and Objective Functions Using Gas Generator and Staged Combustion LPRE Cycles

Abstract

This paper compares various optimization techniques and objective functions to obtain optimum rocket engine performances. This research proposes a modular optimization framework that provides an optimum design for Gas Generator (GG) and Staged Combustion (SC) Liquid Propellant Rocket Engines. This process calculates the ideal rocket engine performance by applying seven different optimization techniques: Simulated Annealing (SA), Nelder Mead (NM), Cuckoo Search Algorithm (CSA), Particle Swarm Optimization (PSO), Pigeon-Inspired Optimization (PIO), Genetic Algorithm (GA) and a novel hybrid GA-PSO technique named GA-Swarm. This new technique combines the superior search capability of GA with the efficient constraint matching capability of PSO. This research also compares objective functions to determine the most suitable function for GG and SC cycle rocket engines. Three single objective functions are used to minimize the Gross Lift-Off Weight and to maximize Specific Impulse and the Thrust-to-Weight ratio. A fourth multiobjective function is used to simultaneously maximize both Specific Impulse and Thrust-to-Weight ratio. This framework is validated against a pump-fed rocket, and results are within 1% of the actual rocket engine mass. The results of this research indicate that PSO and GA-Swarm produce optimum results for all objective functions. Finally, the most suitable objective function to use while comparing these two cycles is the Gross Lift-Off Weight.

1. Introduction

A Launch Vehicle (LV) design requires complex decisions to be made from the start of the design phase. Mission designers must decide on inter-related factors, which include the most feasible engine design and propellants to reach the mission’s requirements, environmental impact, manufacturing processes, overall cost and the availability and skillset of human resource required to manufacture, test and deploy a launch vehicle [1]. These factors are all inter-connected, highlighting the intricacies involved in making major decisions at an early stage of the design process.

A typical rocket engine uses solid, liquid or hybrid propellants. Compared to solid or hybrid propellants, liquid propellants have the following advantages: (1) relatively high energy density, which leads to a higher Specific Impulse ($I_{sp}$), resulting in increased engine efficiency; (2) propellant storage tanks require a reasonably low volume; (3) engine throttling (if required) is easily achieved by varying the propellant flow rate in order to control the thrust; (4) allowance for multiple burns as the engine can be repeatedly shut down and restarted [2]. Liquid Propellant Rocket Engines (LPREs) are classified by their propellant feed system: either pressure fed or pump fed. Mission designers prefer pump-fed LPREs due to the higher combustion chamber pressure ($p_c$) produced. Higher $p_c$ results in a relatively smaller engine; therefore, the payload carrying capacity of the overall LV increases. However, pump-fed engines, also known as turbopump engines, come with varying complexities. There are three types of turbopump LPRE cycles: Gas Generator (GG), Staged Combustion (SC) and Expander. GG and SC cycles contain a secondary combustion chamber where tapped-off fuel and oxidizers are ignited. The resulting exhaust gases are used to power the turbine. This turbine powers the pumps, which raise the propellant pressure to the required $p_c$. In SC, the exhaust gases are redirected back to the main chamber for complete combustion. Expander cycles do not contain a secondary combustion chamber. Such engines are only used for cryogenic propellants as their heat of vaporization is greater; therefore, more energy can be extracted and used in mechanical work. Thus, Expander cycle LPREs are only used as upper stage engines in LVs. Hence, for single-stage LVs or first-stage engines, there is often a tradeoff between GG and SC engines: SC engines typically achieve higher $p_c$; therefore, the LV can carry greater payloads. However, this increase in payload comes at the expense of greater design complexity.

The design phase for an LV consists of three stages: conceptual, preliminary and detailed design. The propellant pair and engine cycle are finalized in the conceptual design phase. The mission designer’s main aim is to maximize the LPRE performance while simultaneously minimizing the overall LV weight in order to increase the payload [3]. According to Hammond [4], the Life-Cycle Cost (LCC) of an LV follows the Pareto Principle: 80% of a vehicle’s LCC is finalized at the conceptual design stage.Thus, any changes in the preliminary or detailed design stages do not result in exponential improvements. Therefore, to minimize the risk of major revisions at a later stage in the design cycle, research recommends a thorough analysis of all major sub-systems and components at the conceptual design phase [5]. The use of multiple trade-off studies during the conceptual design stage enables mission designers to quantitatively determine the effect of mission requirements on the engine’s performance and overall LV size. Such studies aid mission designers in understanding the effect of interconnected variables such as propellant mass, LV weight, $p_c$, thrust and manufacturing complexity [6].

Trade-off studies rely on mathematical tools to model the propulsion sub-systems and provide LV mass estimates. These tools are based on Mass Estimate Relationships (MERs), which are empirical models derived from parametric equations between thrust and engine mass [7]. Therefore, the accuracy of the model is proportional to the number of data points used. A review of MERs determined that there is no single set of relationships that is more accurate [8]. Wildvank [9], therefore, recently developed a structural dry mass model for an LV, with the aim of moving away from empirical-based mass estimates.

Optimization techniques initially used on LPRE cycles included non-gradient- and gradient-based methods such as Genetic Algorithm (GA) and Sequential Quadratic Programming (SQP) for both entire cycles as well as specific sub-systems [10,11,12]. Akhtar et al. and Bayley et al. [13,14] used GA-based techniques to optimize a single-stage and a four-stage LV, respectively. Da mota [15] utilized SQP to optimize two LPRE cycles with the objective of maximizing payloads. He concluded that algorithms such as Particle Swarm Optimization (PSO) should also be studied. Zhou et al. [16] modelled the optimal starting and regulating characteristics of an electric pump for LPRE using the GA technique. A fuzzy Multi-Objective GA (MOGA) technique has been applied to a GG cycle [17] as well as on a hybrid rocket [18] with the objectives of maximizing altitude and minimizing Gross Lift-Off Weight (GLOW). Orgeira-Crespo et al. [19] compiled a database of existing LPREs and examined conceptual multistage LVs based on those existing engines. A GA technique was used on these conceptual LVs to optimize the GLOW based on different trajectories. Castellini et al. [20] compared different techniques to optimize the LV ascent trajectory and concluded that GA and PSO techniques were the most feasible techniques to use.

A surrogate model-based two-step method was used by Mack et al. [21] to optimize a radial turbine in an Expander cycle LPRE. First, Response Surface Approximation was used to search the solution space to develop an approximate model. Secondly, GA was used to refine the search area and obtain a global optimum point. An alternative to this method is CFD, which is computationally intensive and time-consuming. Cui et al. [22] also used a surrogate-based model to optimize the trajectory of an air-launched solid rocket. Neural Network (NN)-based models were applied by Dresia et al. [23] to LPRE problems such as estimating the fatigue life of reusable thrust chambers and estimating the maximum temperature in nozzle cooling channels. In both cases, NN models were created as large sets of training data were available: 120,000 data points to estimate fatigue life and 20,000 CFD simulations to calculate the maximum temperature. NN-based controllers are proposed by Yu et al. [24] to ensure steady-state spacecraft formation when in orbit. Therefore, similarly to [25], surrogate and NN models are beneficial when problems are computationally intensive and already have numerous data sets available. For the purposes of this research, the numerical calculations are not computationally intensive; therefore, the benefits of such techniques are minimized and non-gradient- and gradient-based algorithms can be directly applied.

When determining GLOW for an LV, according to the best of our investigations, researchers, apart from [19], used mass models based on MERs to provide estimates for the masses. While minimizing GLOW and maximizing payload have been used as Objective Functions (OFs), to the best of our knowledge, other OFs such as maximizing $I_{sp}$ or Thrust-to-Weight ratio ($T/W$) have not been investigated. In addition, the optimization models mentioned previously require Thrust (F) and $p_c$ as inputs. This results in running the optimizer numerous times at different user defined F and $p_c$ values in order to obtain a Pareto front of optimal points. Furthermore, although newer techniques such as the Cuckoo Search Algorithm (CSA) and Pigeon-Inspired Optimization (PIO) have successfully been used in various aerospace applications [26], to our knowledge, these techniques have not been used in LPREs.

This research has three aims. First, to estimate the mass and dimensions of a pump-fed cycle LPRE using mathematical model and not via statistical or empirical estimates. Second, to compare seven techniques, including a novel hybrid GA-PSO technique (GA-Swarm), and four OFs on GG and SC cycles and determine the most suitable technique and function for each cycle. Thirdly, a comparison is made between the two cycles for a given mission.

The structure of this paper is as follows: a description of the modelling and optimization process used in this research is provided in Section 2; a validation of the mass model and results of the optimization algorithms and OFs are discussed in Section 3; conclusions and future work are provided in Section 4.

2. Methodology

A flow chart describing the process used to estimate the mass and dimensions of a pump fed LPRE is shown in Figure 1. The framework developed in this research has a modular format. Therefore, the propulsion sub-systems are ’assembled’ in various combinations to model different cycles.

The required inputs are propellant pair, the LPRE cycle to be modelled (SC or GG) and the mission. The mission as defined in this process includes the required change in velocity ($\Delta V$) and the required altitude (H). In this research, $p_c$, F and burn time ($\Delta t$) are calculated during the optimization process, unlike previously mentioned models in Section 1 where $p_c$ and F were required inputs.

2.1. LV Model

Once the propellant pair is entered, Chemical Propellant Equilibrium Program (CPropEP), an external software based on NASA calculations [27], is used to calculate propellant thermodynamic properties. Detailed information on CPropEP is given in [28,29]. This software calculates a range of output variables, from which the characteristic velocity ($c^{\ast }$), thrust coefficient ($C_F$), combustion chamber temperature ($T_c$), isentropic expansion factor ($\gamma$), molar mass (M) and the expansion ratio ($ϵ$) are required for further calculations in this research.

As $p_c$ is not fixed and as optimum $p_c$ is required, the above variables would have to be recalculated in every iteration during the optimization process. In order to reduce computing times, an additional code was written to calculate the optimum Mixture Ratio (MR) for the propellant combination at different values of $p_c$. The data were fitted using non-linear curve-fitting. The equations generated are then used to recalculate the thermodynamic variables in each iteration during optimization instead of repeatedly running CPropEP again. Therefore, the advantages of this additional step are as follows:

• Time saved during the optimization process as CPropEP is not repeatedly run for each iteration. Instead, it only runs once to generate equations for each of the required variables.

• The optimum MR for each $p_c$ is calculated to achieve the maximum $c^{\ast }$ at that specific $p_c$. Maximizing $c^{\ast }$ is important as it indicates a higher combustion efficiency and an increase in engine performance.

The thermodynamic properties are then used in the engine model developed for the purpose of this research.

This method also minimizes the optimization processing time as otherwise CPropEP would load and run for each value of $p_c$. MR is calculated via the equation below:

(1)$MR=a+b\times p_c^{0.125}+\frac{c}{p_c}$

The output variables $c^{\ast }$, $C_F$, $T_c$, M and $\gamma$ are calculated in succession using the following equation:

(2)$y=a+b\times p_c^{0.125}-\frac{c}{p_c}$

where y alternately represents $c^{\ast }$, $C_F$, $T_c$, M and $\gamma$.

Equation (3) is used to calculate $ϵ$.

(3)$ϵ=a+b\times p_c^{0.769}$

These equations and the values of the variables a, b and c are unique to the combination of a given propellant pair at a given exit pressure. If the performance of an LPRE cycle needs to be analyzed with a different propellant pair or the performance of the propellant pair needs to be evaluated at a different value of $p_e$, data need to be generated from the start to obtain the required values for each equation.

2.1.1. Propellant Mass and Tank Calculations

A forward calculation is used in the optimization techniques, which requires Thrust F, $p_c$ and burn time ($\Delta t$) to be specified. Therefore, initial values are given; however, the final values are obtained as an optimized output. Using the output values calculated in Section 2.1 and standard equations as defined in [30], specific impulse, $I_{sp}$, propellant mass flow rate, ${\dot{m}}_p$, and overall propellant mass, $M_{prop}$, are calculated using the following equations:

(4)$I_{sp}=c^{\ast }\times \frac{C_F}{g_0}$

(5)${\dot{m}}_p=\frac{F}{I_{sp}\times g_0}$

(6)$M_{prop}={\dot{m}}_p\times \Delta t$

where $g_0$ is the standard value of Earth’s gravity defined as 9.81 m/s${}^2$. The multiplication by $g_0$ is for the purpose of unit conversion only as $I_{sp}$ is defined as total impulse per unit weight of propellant. Using the optimum MR calculated in the previous section, the mass of fuel, $M_f$, is calculated via the following.

(7)$M_f=\frac{M_{prop}}{1+MR}$

As fuel has been predefined, its density is already known; therefore, the volume of fuel and, hence, volume of fuel tank can be calculated as per Equation (9). An additional 10% volume is added to account for ullage volume.

(8)$v_f=\frac{M_f}{{\rho }_f}$

(9)$v_{fTank}=1.1v_f$

Assuming the tank is cylindrical in shape, with a Length:Diameter (L/D) of 7, rearranging the standard equation for volume of a cylinder, the radius and length are calculated in the following equations.

(10)$r_{fTank}=\sqrt[3]{\frac{v_{fTank}}{14\pi }}$

(11)$l_{fTank}=\frac{v_{fTank}}{r_{fTank}^2\pi }$

The selected propellant tank material for both propellants is Stainless Steel 316 [30]. Safety factors of 1.25 and 1.5 are, respectively, applied to the yield strength and ultimate strength to calculate the maximum allowable operating stress, ${\sigma }_{tank}$.

The total operating pressure is calculated as given below where the pressure drops due to the injector head assumed to be 10 bar, and an additional 10% is added due to losses in the fittings.

(12)$p_{tot}=p_c+p_{injDrop}+0.1p_c$

As the tank material and allowable operating stress is known, the equation for Hoop Stress as defined in [30] is rearranged to obtain the tank’s thickness via the following equation.

(13)$t_{fTank}=\frac{p_{tot}\times r_{fTank}}{{\sigma }_{tank}}$

Finally, the propellant fuel tank mass is determined using Equation (14):

(14)$M_{fTank}=2\pi \times r_{fTank}\times l_{fTank}\times t_{fTank}\times {\rho }_{tank}$

where ${\rho }_{tank}$ is already known as the tank material is pre-defined. Equations are then repeated to calculate the oxidizer mass, $M_{ox}$, and oxidizer tank mass, $M_{oxTank}$. Instead of Equation (7), the following equation is used to calculate the oxidizer mass.

(15)$M_{ox}=M_{prop}-M_f$

In this framework, it is assumed both tank radii are the same. Therefore, once the radii of both tanks are calculated, the maximum value is taken as the radii for both, and the tank lengths and thicknesses are recalculated accordingly.

2.1.2. Pressurant Mass and Tank Calculations

The selected pressurant is Helium, He, and the pressurant tank material is Ti 6Al-4V [30]. The initial tank pressure, $p_{He}$, is set at 450 bar. The number of moles of He required to pressurize the fuel and oxidizer tanks is calculated using the Ideal Gas Equation:

(16)$nH{e}_f=\frac{p_{tot}\times v_f}{RT}$

(17)$nH{e}_{ox}=\frac{p_{tot}\times v_{ox}}{RT}$

where R is the molar gas constant, 8.314 J/Kmol, and temperature, T, is assumed at room temperature, 298.16K. The volumes for fuel ($v_f$) and oxidizer, ($v_{Ox}$) were previously calculated using Equation (8).

The ratio of He left in the tank, $rati{o}_{H{e}_{left}}$, is assumed to be 2. The pressure of He left in the tank, $p_{H{e}_{left}}$ is given by the following.

(18)$p_{H{e}_{left}}=p_{tot}\times rati{o}_{H{e}_{left}}$

The minimum number of moles of He, (nHe), therefore, required to pressurize both propellant tanks is calculated as follows.

(19)$nHe=\frac{\frac{(nH{e}_f+nH{e}_{ox})\times p_{He}}{p_{H{e}_{left}}}}{\frac{p_{He}}{p_{H{e}_{left}}}-1}$

The mass of He, $M_{He}$, is obtained via converting the number of moles of He to kg.

(20)$M_{He}=\frac{4nHe}{1000}$

The Ideal Gas Equation is again used to calculate the volume of He, $v_{He}$, in the tank. This is assumed to be the inner tank’s volume. The pressurant tank is assumed to be spherical; therefore, using the standard equation for volume of a sphere, the inner radius is calculated via the following.

(21)$r_{HeTan{k}_i}=\sqrt[3]{\frac{v_{He}}{4\pi }}$

Given the allowable stress, S (assumed 1100 MPa), and joint efficiency, E (0.9), of the material and the equation from [30], the tank’s thickness, $t_{HeTank}$, is determined as follows.

(22)$t_{HeTank}=\frac{p_{He}\times r_{HeTan{k}_i}}{2S\times E-0.2p_{He}}$

The tank’s outer radius, $r_{HeTan{k}_o}$, is simply the following.

(23)$r_{HeTan{k}_o}=r_{HeTan{k}_i}+t_{HeTank}$

The He tank outer, $v_{HeTan{k}_o}$, and inner $v_{HeTan{k}_i}$ volumes are calculated using the volume of a sphere and the respective inner and outer tank radii given in equations. The mass of the pressurant tank, $M_{HeTank}$, is therefore calculated below with an additional 20% mass as the safety margin.

(24)$M_{HeTank}=\frac{1.2(v_{HeTan{k}_o}-v_{HeTan{k}_i})}{{\rho }_{HeTank}}$

2.1.3. Combustion Chamber Mass Calculations

The combustion chamber’s shape is assumed to be cylindrical with a converging-diverging nozzle attached at the lower end. The assumptions, based on those given in [30], used are chamber area ($a_c$) equivalent to three times the throat area ($a_t$), nozzle characteristic length ($L^{\ast }$) of 1.143 m, nozzle converging section is a truncated cone and the chamber and nozzle thickness is 3 mm. A 45°angle is assumed for the converging nozzle section and 15°angle for the diverging section. The nozzle and chamber’s material comprises stainless steel. Therefore, once $a_t$ is calculated as in Equation (25), and as material $\rho$ is known, standard equations for volume of a cylinder and truncated cone are used to calculate the remaining dimensions and, consequently, the nozzle’s mass, $M_{nozzle}$, and combustion chamber’s mass, $M_c$.

(25)$a_t=\frac{c^{\ast }\times {\dot{m}}_p}{p_c}$

2.1.4. Pre-Burner and Gas Generator Dimensions and Masses

The pre-burner and gas generator dimensions and mass ($M_{GG}$ and $M_{PB}$ are calculated in a similar manner to the main thrust chamber given in the previous section). The assumptions about the material, shape, converging angle and nozzle thickness are as given in Section 2.1.3.

The mixture ratio for both GG and PB is fixed at 0.5. CPropEP is used to calculate the subsequent thermodynamic parameters required for mass and dimension calculations.

2.1.5. GG Parameters

For the GG cycle, the combustion chamber pressure of the GG, $p_{cGG}$, is assumed to be the same as the main chamber, $p_c$.

The propellant mass flow rate for GG, ${\dot{m_p}}_{GG}$, is assumed to be 2% of the main chamber flow rate, as given in [30].

2.1.6. SC Parameters

For the SC cycle, the combustion chamber pressure of the PB, $p_{cPB}$, is calculated via the critical pressure ratio, as given in [30]:

(26)$p_{cPB}=\frac{P_{ratiot}\times 1.1\times (p_c+\Delta p_{losses})}{0.57}$

The propellant mass flow rate for PB is assumed to be fuel rich [30]. It is assumed that the entire fuel volume flows through the PB. The total fuel mass flow rate through the PB is assumed at $1.5\times {\dot{m}}_{fPB}$. As the mixture ratio is already known, as is ${\dot{m}}_{fPB}$, the PB oxidizer mass flow rate is then calculated.

2.1.7. Engine Balancing

The design of a turbopump system depends on the power balance between the pumps and turbine. The mission designer must ensure that the power generated by the turbine matches the power required by the pumps.

The power required by each pump is a function of the volumetric flow rate, Q, change in pressure, $\Delta p$, and the pump efficiency, $\eta$, as shown in the following equation.

(27)$P_{pum{p}_{ox,f}}=\frac{Q\times \Delta p}{{\eta }_p}$

The power required by each pump is calculated separately and summed to provide the total power required by the pumps, $P_P$.

The actual power produced by the turbine is calculated via the following:

(28)$P_T={\eta }_t\times \dot{m_T}\times c_p\times TTIn\times (1-{(\frac{1}{PRatioT})}^{\frac{{\gamma }_{GG}-1}{{\gamma }_{GG}}})$

where $et{a}_t$ is the turbine efficiency, $\dot{m_T}$ is the turbine mass flow rate, $c_p$ is the specific heat capacity at constant pressure of the gas entering the turbine, $TTIn$ is the Turbine Inlet Temperature, $PRatioT$ is the turbine pressure ratio ($\frac{TurbineInletPressure}{TurbineDischargePressure}$), and $\gamma$ is the ratio of the specific heats.

The equations given above undergo several iterations until the power produced by the turbine matches the power required by the pumps.

2.1.8. GG Power Balance

GG is an open cycle LPRE. The turbine pressure ratio across the turbine is kept fixed at 20 [30]. In addition, the turbine mass flow rate is assumed equal to the GG mass flow rate. As mentioned in Section 2.1.5, the mixture ratio is fixed as is the GG flow rate. Therefore, the oxidizer and fuel mass flow rates through the GG can be calculated as per Section 2.1.5.

2.1.9. SC Power Balance

For closed cycles, such as SC and Expander, the LPRE tool keeps the turbine pressure ratio as a variable to be optimized. According to [30], this ratio is approximately 1.5 for SC and Expander cycles. The turbine must generate enough power to drive two oxidizer pumps and one fuel pump. The PB mass flow rate and individual propellant mass flow rates through the PB are calculated as mentioned in Section 2.1.6. The volumetric flow rate through the secondary oxidizer pump, $Q_{ox2}$ is taken as half of the total fuel volumetric flow rate, $Q_f$. The turbine mass flow rate is assumed to be equal to the PB mass flow rate; however, this flow rate is greater in magnitude compared to the GG flow rate as the entire fuel is directed to the PB.

2.1.10. System Mass and Dimensions

The overall system mass, $M_{sys}$, is then a sum of the combustion chamber, nozzle, tanks, propellant and pressurant masses, turbopump assembly, Pre-Burner (or Gas Generator) and the payload, as shown in Equation (29). This is also known as Gross Lift-Off Weight or GLOW. Calculations for system mass and dimensions are as per the model previously described in [31].

(29)$\begin{array}{cc}{M}_{sys}=\hfill & M_{HeTank}+M_{He}+M_{ox}+M_{oxTank}+M_f+M_{fTank}\hfill \\ \hfill & +M_{nozzle}+M_c+M_{turbine}+M_{pumps}+M_{GG}+M_{payload}\hfill \end{array}$

For the SC cycle, $M_{PB}$ would replace $M_{GG}$ and the secondary oxidizer pump, $M_{ox2}$, would also be included.

2.2. Flight

This section describes the models used to calculate the forces acting on the LPRE (and hence the amount of Thrust, F, required) as it ascends through the atmosphere to reach the required mission altitude. Changes in gravity, air density and atmospheric pressure and temperature are taken into account. The framework assumes vertical lift-off and flight.

Once engine parameters are calculated, F and Drag (D) are calculated at pre-determined intervals of time to determine the performance of the LPRE from takeoff until reaching the required H.

2.2.1. Atmospheric Model

Atmospheric temperature, ($T_a$), and atmospheric pressure, ($p_a$), vary with altitude, H. This variation influences the air density (${\rho }_a$), which in turn affects the drag on any vehicle passing through the Troposphere, lower Stratosphere and upper Stratosphere. In this research, the US Standard Atmosphere Model, 1976 [32], is used to calculate the pressure, temperature and density with increasing altitudes as per the following equations.

For the Troposphere with an altitude below 11,000 m, $T_a$ and $p_a$ are calculated via the following:

(30)$T_a=15.04-0.00649\times H$

(31)$p_a=101.29\times {(\frac{T_a+273.1}{288.08})}^{5.256}$

When the LPRE is travelling through the lower Stratosphere where the altitude is between 11,000 m and 25,000 m, the following set of equations are used to obtain $p_a$ and $T_a$:

(32)$T_a=-56.46$

(33)$p_a=22.65e^{1.73-(0.000157\times H)}$

The upper Stratosphere lies above 25,000 m. For these higher altitudes, ($T_a$) and ($p_a$) are determined as per the following equations:

(34)$T_a=-131.21+0.00299\times H$

(35)$p_a=2.488\times {(\frac{T+273.1}{216.6})}^{-11.388}$

The three equations are given above for calculating $T_a$ return values in °C. Therefore, an additional step is required to add 273.1 to each value to convert to Kelvin. Similarly, the values of $p_a$ are in kPa and must be converted to Pa in order to ensure all units are uniform throughout the framework.

The air density, ${\rho }_a$, is a function of $p_a$ and $T_a$ and is calculated via the followin.

(36)${\rho }_a=\frac{p_a}{(0.286\times (T_a+273.1)}$

2.2.2. Gravitational Model

Newton’s Law of Universal Gravitation is used to determine the change in acceleration due to gravity, $g_0$, at each change in altitude. The values of the constants G, Gravitational constant, $M_e$, Earth’s mass and $R_e$, Earth’s radius, are given in Table 1. These values are used in the following equation to determine $g_0$ at a particular altitude.

(37)$g_0=\frac{G\times M}{{(R+H)}^2}$

2.2.3. Flight Model

The values calculated in Section 2.2.1 and Section 2.2.2 are used to determine the Thrust, F, and Drag, D, during the LPRE’s flight from takeoff to the required altitude, H. For this framework, a time interval, $\delta t$, of 0.1 s is used when calculating the velocity and height during the flight. The Drag Equation (assuming a fixed drag coefficient ($C_{D0}$) of 0.2) and General Thrust Equation are, respectively, used to calculate the Drag and Thrust, as given in the following equations:

(38)$D=\frac{{\rho }_a\times v^2\times A_s\times C_{D0}}{2}$

(39)$F={\dot{m}}_p\times v_e+a_e\times (p_e-p_a)$

where $v_e$ and $R_g$ are defined as per [30] using the universal gas constant R assumed to be 8.314 J/Kmol. $T_c$, M and $\gamma$ are calculated via Section 2.1 and used as follows.

(40)$R_g=\frac{R}{M}$

(41)$v_e=\sqrt{\left(\frac{2\times \gamma \times R_g\times T_c}{\gamma -1}\right)\times \left(1-{\left(\frac{p_e}{p_c}\right)}^{\left(\frac{\gamma -1}{\gamma }\right)}\right)}$

The exit pressure ($p_e$) is fixed and assumed at 0.7 bar for this research study. Therefore, for each interval of time, $p_a$ is recalculated as per the equations given in Section 2.2.1 and then used in the Thrust equation. The second term in the Thrust Equation is the pressure component, which varies as $p_a$ changes as the LPRE gains altitude. F, therefore, increases when the nozzle exit velocity, $v_e$, is at a maximum point (as propellant mass flow rate, ${\dot{m}}_p$ is fixed).

In the Drag Equation, $A_s$ is assumed as the frontal surface area. This is a fixed value, and as $C_{D0}$ is kept at a fixed value, changes in D are only due to changes in ${\rho }_a$ and velocity, v.

The system mass at each interval of time and D, F, $T_a$, $p_a$, ${\rho }_a$ and gravity, $g_0$, calculated earlier are used to determine the system mass, $M_{sys}$, velocity, v, and subsequently the altitude, H, at each time interval as per the following equations.

(42)$M_{sys}=M_{sys}-{\dot{m}}_p\times \Delta t$

(43)$v=v+\frac{F}{M_{sys}}\times \Delta t-g_0\times \Delta t-\frac{D}{M_{sys}}\times \Delta t$

(44)$h=h+v\times \Delta t$

2.3. Optimization

There is a need to quantitatively compare different optimization techniques to evaluate which methods are most suitable for a particular problem. There may not be a single method that consistently performs better than others for different objective functions. This is why the following techniques are evaluated according to their performance: Simulated Annealing (SA), Nelder Mead (NM), Cuckoo Search Algorithm (CSA), Pigeon-Inspired Optimization (PIO), Particle Swarm Optimization (PSO), Genetic Algorithm (GA) and a novel hybrid GA-PSO technique named GA-Swarm, which is described in further detail in Section 2.3.1.

The importance of selecting a relevant OF for use in multidisciplinary design studies cannot be overemphasized. Da Mota et al. [15] tried to minimize the Gross Lift-Off Weight (GLOW), the structural mass (SM) and the expected total cost to investigate the effect of changing the OF using a single optimization technique. Dresia et al. [33] investigated the effect of minimizing GLOW or Structural Mass (SM) on a GG cycle using GA techniques and Akhtar et al. [13] studied the effect of GLOW using GA on a pressure-fed cycle.

These seven optimization techniques and four Objective Functions (OFs) are used to optimize each LPRE cycle and to determine the most suitable OF and technique. A list of the four OFs, design variables and constraints is provided in Table 2.

2.3.1. Genetic Algorithm-Swarm (GA-Swarm)

GA-Swarm is a novel technique developed in this research by incorporating the ’swarm’ behaviour of PSO within GA. GA-PSO hybrid techniques have been used previously [34,35,36]. However, these techniques either incorporate mutation within PSO to ensure the technique can escape local optima or apply the GA ’fitness’ function to swarm particles in order to divide them and then apply crossover and mutation functions on the ’better’ half of the population.

In GA-Swarm, as shown in Figure 2, PSO is incorporated into the GA algorithm by ensuring that after crossover, the offspring ’swarm’ towards the best solution in order to ensure a more efficient search within the solution space.

The offspring is given a bias toward the best solution if a random number is less than a predefined probability, pSwarm, which is usually equal to 0.2.

(45)$chil{d}_j=(1.0-sSwarm)\ast chil{d}_j+sSwarm\ast gBes{t}_j$

For the strength of bias, sSwarm is taken to be equal to 0.2.

2.3.2. Constraints and Boundaries

There are two constraints for all techniques and OF: $\Delta V$ and H. In order to ensure that constraints are met every time the algorithm is run, a penalty-based method is applied. This weightage factor is provided in Equation (46):

(46)$X=100,000.0\times abs\left(\right(v-\Delta V\left)\right)+10,000.0\times abs\left(\right(h-H\left)\right)$

where there is a penalty for each constraint, as shown in Equation (47):

(47)$ObjectiveFunction=RawCost+X$

where the Objective Function is as given in Table 2: i.e., GLOW, $I_{sp}$, $T/W$ or maximizing both $I_{sp}$ and $T/W$. In every iteration, there is a check whether the velocity and height calculated match the required values of $\Delta V$ and H, respectively. Both constraints are multiplied by a penalty coefficient of 100,000 for $\Delta V$ and 10,000 for H. If the constraints are violated, then the measure of violation will be non-zero. If the constraints are not violated then the measure of violation will be zero as the penalty coefficient multiplication has no effect.

The upper and lower bounds for the three variables: F, $p_c$ and $\Delta t$ are given in Table 3. These variables are the same for both GG and SC cycles.

The number of iterations per run for all techniques was kept at 1000, except for PSO and GA-Swarm where each run had 150 iterations. For PSO, the number of particles in the ‘swarm’ was 50; the number of ‘pigeons’ in PIO was 128; CSA algorithm contained 125 ‘nests’, and both GA and GA-Swarm algorithms had 150 ’chromosomes’.

3. Results

3.1. Model Validation

Shahab-3 is an Iranian-designed single-stage GG rocket. The first documented flight was in 1998; however, this rocket is currently still in use. Table 4 lists the approximate engine performance parameters as sourced from the open literature [37].

The values of the variables used to calculate $c^{\ast }$, $C_F$, $T_c$, $\gamma$, M and $ϵ$ as per Equations (1)–(3) are provided in Table 5. A graphical visualization of this data is included in Figure A1.

Two optimized results are included in Table 6 using PSO as the optimization technique. PSO is an established optimization technique and has provided consistent results for all four OFs used in this research. Hence, PSO is used to optimize the given model. The first result gives the lowest GLOW as per the OF used in the optimization technique. The second result also used GLOW as OF; however, the $p_c$ returned is close to the claimed $p_c$ value of Shahab-3 as per the literature. The reason behind including this second set of results is to highlight that optimization techniques can return a range of values within the solution space. Therefore, according to this particular optimization technique and the OF used, this engine could have a lower GLOW but with a higher $I_{sp}$ for the same propellant pair and cycle. The result with a similar $p_c$ (48 bar) to the actual engine also calculated a similar GLOW compared to the reported engine values.

Differences in results compared to the actual engine are due to a range of factors. The choice of material for tanks and engine, L/D used, turbine speed and actual turbine power produced in the actual engine are possible causes that account for the difference in values.

A visual summary of the optimized parameters is given in Figure 3.

3.2. Gas Generator Cycle: Comparison of Techniques and Objective Functions

Four different objective functions and seven optimization techniques are tested on a RFNA-RP1 GG cycle with mission constraints $\Delta V$ 3800 m/s and an H of 120 km, as shown in Figure 4, with a payload of 1000 kg. A ’box’ is included in each graph to depict the operating range of $p_c$ generally ascribed to GG cycle engines, as compiled by [11]. The results of three OF (minimum GLOW, maximum $I_{sp}$ and maximum $T/W$) are provided in Figure 4, plotting the optimum values of each OF against $p_c$. The results of the fourth OF are displayed separately in Figure 5.

The purpose of testing multiple objective functions and techniques is to determine the most suitable function and technique for the GG LPRE cycle.

3.2.1. Nelder Mead

With reference to Figure 4a, the NM technique returned results that met the mission’s constraints for two out of three OFs. Therefore, the results for maximizing the Thrust-to-Weight ratio have been excluded. A probable reason for this inability to meet the constraints is due to NM’s susceptibility to the initial guess. This technique searches for the local maximum within the vicinity of the initial guess. Therefore, if the initial guess is close to the global maximum, the probability of NM converging on the local maxima is higher than when the initial guess is further away from the global maximum.

Out of these results for NM, the best-performing OF is $I_{sp}$, as the optimum result both clearly lies within the accepted operating range of $p_c$ for a GG cycle, and occurs at a $p_c$ of 71 bar, which is lower than the minima obtained using GLOW OF. For GLOW OF, there are distinct bands of results at 20, 30 and 40 bar $p_c$, with the minima lying at 97 bar. This minima lies within the operating $p_c$ range for a GG cycle; however, it is relatively higher than the $p_c$ obtained using $I_{sp}$ as the OF. Using $I_{sp}$ as an OF provided a relatively larger cluster of results between 11 and 18 bar $p_c$, which is an infeasible operating range for GG cycle. Another smaller band of results between 40 and 45 bar $p_c$ coincided with a similar band of results when using GLOW as the OF. Results of both OFs demonstrate the presence of multiple local minima in the solution space.

3.2.2. Simulated Annealing

Figure 4b shows the results for three OFs using the Simulated Annealing technique. Using GLOW as the OF provides a band of results between 70 and 74 bar $p_c$. The results also follow a fairly linear path within the GG cycle $p_c$ operating range. However, this technique and OF return results between 12 and 13 bar and one result at 479 bar, all outside the accepted $p_c$’s operating range for this LPRE cycle.

Similarly to the NM technique, using $I_{sp}$ as the OF provides a relatively greater number of results below the acceptable $p_c$ range for this cycle as compared to results within the accpetable $p_c$ range. In contrast, the T/W function follows the same trend as GLOW by providing a greater range of results within the acceptable $p_c$ range. The T/W function provides a greater number of results between 43 and 60 bar. These bands demonstrate SA’s strength in avoiding local maxima.

These technique required multiple runs in order to generate a sufficient data set for analysis as the technique did not meet the constraints on every run. Therefore, although results are promising, this technique is not recommended due to the additional computational time involved in generating numerous runs and filtering out erroneous data.

3.2.3. Cuckoo Search Algorithm

Looking at Figure 4c, the majority of CSA results lie between $p_c$ of 29 and 86 bar for all three OF. There is relatively less variation between the minimum and maximum $p_c$ for all three OFs, as compared to NM and SA. The maximum $p_c$ achieved with any of the three OF using CSA is 84 bar, unlike NM and SA where a $p_c$ of greater than 90 bar was noted. Similarly to SA, the three OFs provide an optimum result within a relatively close $p_c$ range: minimum GLOW at 86 bar, maximum T/W at 84 bar and maximum $I_{sp}$ at 77 bar. This narrow range of results highlights CSA’s immunity to the objective function.

The values for GLOW lie on a fairly linear curve with distinct clusters of results: 39–41 bar, 56–63 bar, 70–72 bar and 78–83 bar. $I_{sp}$ also provides a cluster of results in a similar $p_c$ range of 69–72 bar. The curve for T/W highlights two clusters: 44–55 bar and 64–75 bar. Therefore, unlike NM and SA, all three OFs return clusters of results near the global optimum result, highlighting the presence of multiple solutions near the optima within the solution space. Apart from GLOW OF, CSA failed to meet the constraints for the other two OFs on numerous runs; hence, multiple runs were required to generate a sufficient data set of results.

3.2.4. Particle Swarm Optimization

In Figure 4d, the trend for all three OFs using the PSO technique follows that of SA and NM with a wide range of values returned. Although the majority of results for all OFs lie within the GG cycle $p_c$’s operating range, both $I_{sp}$ and T/W OFs return results below and above this range. In particular, the global maxima for both OFs lies above the operating $p_c$ range, with a maximum $I_{sp}$ at 192 bar and maximum T/W at 175 bar, whereas GLOW OF returns a global minima at 112 bar.

The results for GLOW OF are contained in two distinct clusters: 30–33 bar and 57–65 bar. $I_{sp}$ also returns a cluster of results at a similar $p_c$ range of 54–56 bar, underlining the existence of multiple local minima in the solution space. T/W also returns results in three clusters: 23–26 bar, 52–59 bar and 81–89 bar.

It is pertinent to mention that this technique returned results for T/W (111 bar) and $I_{sp}$ (114 bar) functions at a value of $p_c$ near 112 bar (the $p_c$ that resulted in minimum GLOW). Therefore, this indicates the presence of an optimal value for each OF in the search space. CSA also meets the mission constraints for each optimization run.

3.2.5. Pigeon-Inspired Optimization

In Figure 4e, similarly to CSA, the results for all three OFs using the PIO technique vary between a relatively short range: 23 bar to 95 bar. Despite multiple runs, the technique failed to meet both constraints on numerous occasions; therefore, only a limited number of values that met both constraints are included in this research study. However, the majority of these limited values lie within the operating $p_c$ range for a GG cycle. GLOW OF returns a minima at 84.5 bar $p_c$, whereas a maximum $I_{sp}$ is attained at 95.7 bar and maximum T/W is attained at 97 bar.

The data points for each OF are in small clusters, demonstrating multiple local minima. With GLOW OF, there is a cluster near the global minima between 79 and 84 bar. The T/W function returns three groups of values at the following: 30–25 bar, 55–61 bar and 93–97 bar. There is one cluster of values for $I_{sp}$ ranging from 48 to 61 bar, which overlaps with a similar (but smaller) cluster using T/W as the OF.

3.2.6. Genetic Algorithm

The results for all OF using the GA technique, as shown in Figure 4f, are similar to CSA as all results lie within a relatively narrow range of $p_c$ (38–89 bar) compared to other techniques. However, unlike CSA, all results for GA lie within the operating $p_c$ range for this cycle. The minimum GLOW that meets the constraints is at 96 bar $p_c$, the maximum $I_{sp}$ occurs at 89.9 bar and the maximum T/W occurs at a lower $p_c$ value of 75 bar.

For $I_{sp}$ as the OF, the results are in three distinct clusters, 38–48 bar, followed by an equal number in the 58–67 bar and then 80–89 bar range. T/W also returns a dense group of results between 62 and 68 bar, as does GLOW (59–68 bar), reinforcing the presence of multiple local minima. This clustering shows that GA does not become stuck in multiple local minima.

3.2.7. GA-Swarm

The results of the final technique, GA-Swarm, are displayed in Figure 4g. Due to the hybrid nature of this technique, it not only meets the constraints for each run but also provides results with minimal variance and within a relatively narrow range of $p_c$, similarly to GA. The maximum $I_{sp}$ that meets the constraints is at 86 bar, the minimum GLOW is at 87 bar $p_c$, and the maximum T/W occurrs at a lower $p_c$ value of 79 bar. These values are similar to those obtained using the GA technique.

Unlike GA, the values obtained using $I_{sp}$ as OF are not grouped as distinct clusters, however, a dense group of results is observed between 63–66 bar $p_c$ and 71 bar. Using T/W as an OF returns three clusters at: 38–45 bar, 63–68 bar and 71–77 bar. The results for GLOW as OF follow a similar trend as the GA technique, with a relatively large cluster between 50 and 59 bar.

3.2.8. Objective Function: T/W vs. $I_{sp}$

The results of a fourth OF used on GG cycle are shown in Figure 5 where maximizing $I_{sp}$ and T/W are both given equal weightage. Values of the global optima for each technique are given in Table 7. Five techniques return results with a maximum $I_{sp}$ between 262 and 272 s; however, only CSA and SA provide theoretically higher performance values for this Objective Function. F is directly proportional to $I_{sp}$; therefore, a linear trend is expected if $T/W$ is plotted against $I_{sp}$, which is verified by the results as shown in Figure 5. This Objective Function is included in this research to validate the seven techniques. All techniques returned a similar trend line, as shown in Figure 5. If a technique had not returned values other than those shown in Figure 5, the Objective Function would, therefore, have invalidated that particular technique.

3.3. Evaluating Objective Functions and Techniques

Results for all techniques are shown in Table 8, Table 9 and Table 10. When comparing the OF results from each technique, some techniques such as CSA, PIO and GA-Swarm are immune to the influence of one technique over another as the results returned are within a narrow range regardless of the OF. However, PIO (along with NM and SA) and CSA will still be regarded as techniques that fail to meet the mission constraints on numerous occasions. GA-Swarm meets the mission constraints on each run.

Using GLOW as an OF provides relatively few results which lie near the global minimum (as clusters are observed in local minima). However, when comparing the results for each technique, results from T/W provided a maxima at a relatively lower $p_c$ when using GA-Swarm compared to the other OF. Based on experience, the end user prefers an OF with emphasis on low $p_c$. Therefore, using T/W as an OF may be more feasible as compared to using GLOW.

3.4. Staged Combustion Cycle: Comparison of Techniques and Objective Functions

Similarly to Section 3.2, the same objective functions and optimization techniques are tested on an RFNA–RP1 SC cycle. Mission constraints and payload remain the same and the results for three OF are shown in Figure 6. A ‘box’ is included in each graph to depict the operating range of $p_c$ generally ascribed to SC cycle engines [11].

This given mission has returned values (for all techniques and objective functions) that lie below the operating $p_c$ range for SC cycles. However, for the purposes of this research, as the SC cycle is compared to GG cycle for a given mission, the overall trends and results for each technique will be discussed.

3.4.1. Nelder Mead

With reference to Figure 6a, unlike for GG cycle, the NM technique returned results that met the mission constraints for all functions. The minimum GLOW is obtained at a $p_c$ of 122 bar, while with $I_{sp}$ as the OF, the maximum $I_{sp}$ is determined at 99 bar and the maximum T/W occurs at 92 bar.

For $I_{sp}$ OF, there are distinct bands of results from 26 to 33 bar, 47 to 52 bar and 65 to 82 bar. These clusters coincide with similar groups of results when T/W is used as OF, highlighting the presence of local minima in the solution space.The results for GLOW as an OF follow a fairly linear trend without clusters of results. There is also relatively less variance in values as compared to the other two OFs.

3.4.2. Simulated Annealing

Figure 6b shows the results for the three OFs using the Simulated Annealing technique. All three OFs return the optimum result within a relatively close $p_c$ range: minimum GLOW at 94 bar, maximum T/W at 96 bar and maximum $I_{sp}$ at 96 bar.

The results for all three OFs using the SA technique display visible clusters of results, once again highlighting the effectiveness of this technique in evading local minima. However, this technique is unable to meet the mission constraints for all OFs on each run, requiring additional runs in order to obtain a minimum data set of points.

3.4.3. Cuckoo Search Algorithm

The results from CSA technique are shown in Figure 6c. This technique follows the same trend as when run for GG cycle, with results for the three OFs lying within a relatively narrow range of $p_c$, between 46 bar and 92 bar.

Another similarity between the results of both cycles using this technique is the presence of multiple clusters using all three OFs. Once again, the clusters are close to the global optimum’s result, which demonstrates CSA’s superior search ability compared to NM and SA.

3.4.4. Particle Swarm Optimization

The trend for all three OFs using PSO technique is shown in Figure 6d. The graph shows similar results for GG cycle using this technique. Although each run meets exactly the mission’s constraints, however, a wide range of values are returned, which eventually converge on the global optimum solution.

This technique also returns optimum results for each function at a relatively higher $p_c$ compared to other techniques: minimum GLOW at 133 bar, maximum T/W at 111 bar and maximum $I_{sp}$ at 110 bar.

3.4.5. Pigeon-Inspired Optimization

The results for PIO, as shown in Figure 6e, are similar to both CSA and PSO. The similarity with CSA lies in the fact that the results for all three OFs using the PIO technique vary between a relatively short range: 30 bar to 104 bar. Similarly to PSO, there are no distinct clusters of results (as shown in CSA, SA and NM); instead, results show a gradual convergence to the optimum solution for each technique. This is also a contrast to the GG cycle where data were shown in clusters. Moreover, PIO managed to meet the criteria in each run, unlike for the GG cycle.

3.4.6. Genetic Algorithm

For SC cycle, the results of all OF with GA techniques are shown in Figure 6f. For GLOW and $I_{sp}$ functions, the technique failed to meet the constraints on numerous runs. There are only a limited number of data points shown. For T/W, the results show clusters, with a relatively dense cluster between 47 and 70 bar. This cluster coincides with two and one data points using GLOW and $I_{sp}$ functions, respectively, demonstrating that for this particular cycle, T/W is a relatively better function for the GA technique.

3.4.7. GA-Swarm

For the hybrid technique, GA-Swarm results are shown in Figure 6g. There is a marked improvement compared to GA results as this technique met the mission constraints for each run, which is a characteristic of the PSO technique. In addition, unlike PSO, GA-Swarm provided results in a relatively narrow range, demonstrating qualities similar to GA. The maximum $I_{sp}$ that meets the constraints is at 93 bar, and the minimum GLOW and maximum T/W both occur at a lower $p_c$ value of 81 bar and 86 bar, respectively.

Unlike GA, the values obtained using $I_{sp}$ as OF are not grouped as distinct clusters; however, a density of results is observed between 63 and 66 bar $p_c$ and 71 bar. Using T/W as an OF returns three clusters at 38–45 bar, 63–68 bar and 71–77 bar. The results for GLOW as OF follow a similar trend as the GA technique, with a relatively large cluster between 50 and 59 bar.

3.4.8. Objective Function: T/W vs $I_{sp}$

In Figure 7, the results of all seven techniques using the fourth OF on SC cycle are displayed. Table 11 displays values for each variable at max $I_{sp}$ and T/W. For this cycle, NM returns the maximum result for $I_{sp}$ and T/W; however, due to NM’s susceptability to the initial guess, these results are disregarded. SA, CSA and GA-Swarm returned similar values of $I_{sp}$ while SA and GA-Swarm also returned same values of T/W. As mentioned in Section 3.2.8, this Objective Function was specifically included in order to further validate each technique. Similarly to Figure 5, all techniques follow a similar trend with $T/W$ linearly increasing with $I_{sp}$. As the same mission was optimized using both cycles, the magnitude of results for GG and SC using this Objective Function is also similar.

3.5. Evaluating Objective Functions and Techniques

The results for all techniques are shown in Table 12, Table 13 and Table 14. The most promising techniques for SC cycle are CSA, PSO, PIO and GA-Swarm. All four techniques consistently met the mission’s constraints and provided a wide range of results. Apart from PSO, the use of GLOW as an OF returns the lowest value of $p_c$ compared to the other two OFs, with F and $I_{sp}$ values remaining nearly equal. Therefore, if the mission designer requires a low $p_c$ that still provides high engine performance, then GLOW should be used as an OF for the SC cycle.

3.6. GG vs. SC Cycle for Given Mission

The primary aim of comparing GG cycle and SC cycle with the same mission requirements is to aid the mission designer in decision making when faced with the task of selecting a particular engine cycle for a mission.

For GG cycles, all techniques using GLOW OF returned a minima between 52 and 112 bar $p_c$ and an average F of 259 kN. The same techniques on SC, however, resulted in the minima at a relatively higher $p_c$ of 98 to 102 bar with a resulting average increase in F to 272 kN. This higher $p_c$ is still within the operating range for an SC cycle engine.

Therefore, for this mission, a GG cycle would be more beneficial as the same GLOW can be achieved at a relatively lower $p_c$ without a significant loss of F. In addition, the additional design, manufacture and operational complexity of a SC cycle engine do not need to meet the requirements of this mission.

4. Conclusions

This research includes a preliminary design of a conceptual pump-fed LPRE cycle, including mass and dimension estimates of all sub-systems using mathematical models and not empirical models. Unlike previous research, the only inputs required for this model were the mission requirements (H and $\Delta V$) and the propellant pair. The model then calculates the optimum set of engine performance parameters. In addition, seven optimization techniques and four Objective Functions were analyzed in order to provide the maximum performance of the cycle for a given mission, propellant pair and payload. An Objective Function simultaneously maximizing both $I_{sp}$ and $T/W$ was used to validate all techniques. One of the optimization techniques, GA-Swarm, used was a novel technique. This research concludes that the best-performing OF for the GG cycle was T/W as the optimized engine performance values for all variables were better compared to using GLOW as an Objective Function. The best-performing Objective Function for SC was GLOW. However, when comparing two cycles, the benchmark OF to consider for both is GLOW.

Looking at the optimization techniques, Nelder Mead is extremely susceptible to initial guess across both cycles and, therefore, cannot be used for an initial optimization. Simulated Annealing shows variance, provides a wide range of values and also returns results in clusters, highlighting that the results are immune to local minima. However, Simulated Annealing does not match the constraints on each run and, therefore, requires additional computational time to generate a data set of values for analysis. Cuckoo Search Algorithm technique works well for SC cycle only. PSO again shows variance by providing a wide range of results for all cycles and missions; however, this technique meets the mission constraints for each run. Pigeon-Inspired Optimization provided optimum performances for the SC cycle only. Similarly to CSA, for the GG cycle, this technique struggled to meet constraints and provided a greater scatter effect in the results. GA provides results within a wider range of acceptability. However, the constraints were not met on multiple instances. Therefore, although this is a robust technique, the additional computational time negates the use of this technique. GA-Swarm, due to its hybrid nature, combines the search ability of GA by efficiently exploring the solution space with the fast convergence of PSO for all cycles and OFs. In conclusion, GA-Swarm performs better than GA for all cycles and both PSO and GA-Swarm provide a set of optimum results for each objective function for all cycles.

Suggestions for future work would include running the optimizer using different propellant pairs to check if GG is still the more feasible cycle for this mission. This model can be modified to model a pressure-fed LPRE cycle and to incorporate multiple stage LVs. A database can be developed by comparing various LPRE cycles against each other for different missions. In addition, the application of GA-Swarm for a variety of engineering problems is also recommended in order to further estabish its robustness.

Author Contributions

Conceptualization, I.Q. and S.S.K.; methodology, S.S.K. and I.Q.; software, S.S.K. and M.A.A.; validation, S.S.K.; formal analysis, S.S.K.; writing—original draft preparation, S.S.K.; writing—review and editing, M.U.S., R.F.S. and S.R.Q.; supervision, I.Q. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Optimization and modelling flow chart.

Figure 2. Flowchart depicting the GA-Swarm Algorithm.

Figure 3. Optimized Engine Parameters, dimensions and masses for Shahab-3.

Figure 4. Optimization results for the GG Cycle with mission [Forumla omitted. See PDF.] = 3800 m/s and H = 120 km using all 7 techniques and 3 OFs. (a) Nelder Mead (NM). (b) Simulated Annealing (SA). (c) Cuckoo Search Algorithm (CSA). (d) Particle Swarm Optimization (PSO). (e) Pigeon-Inspired Optimization (PIO). (f) Genetic Algorithm (GA). (g) GA-Swarm.

Figure 5. GG Cycle: T/W vs. [Forumla omitted. See PDF.] for all techniques.

Figure 6. Optimization results for SC Cycle with mission [Forumla omitted. See PDF.] = 3800 m/s and H = 120 km using all 7 techniques and 3 OFs. (a) Nelder Mead (NM). (b) Simulated Annealing (SA). (c) Cuckoo Search Algorithm (CSA). (d) Particle Swarm Optimization (PSO). (e) Pigeon-Inspired Optimization (PIO). (f) Genetic Algorithm (GA). (g) GA-Swarm.

Figure 7. SC Cycle: T/W vs [Forumla omitted. See PDF.] for all techniques.

Table 1

Values used in Newton’s Law of Universal Gravitation.

Variable	Value
G	$6.67408\times {10}^{-11}$
$M_e$ (kg)	$5.9722\times {10}^{24}$
$R_e$ (m)	$6.371\times {10}^6$

Table 2

OF and Variable names and descriptions.

Objective Function	Description
GLOW	Gross Lift-Off Weight (total mass of LV including payload)
$I_{sp}$	Specific Impulse
$T/W$	Thrust-to-Weight ratio
TWISP	Simultaneously maximize $I_{sp}$ and $T/W$
Design Variable	Description
F	Thrust
$p_c$	Combustion Chamber Pressure
$\Delta t$	Burn time
Constraints	Description
$\Delta V$	Change in velocity
H	Required altitude

Table 3

Values of Variables used for GG and SC cycles.

Bounds	F (kN)	${\mathit{p}}_{\mathit{c}}$ (bar)	$\mathbf{\Delta }\mathit{t}$ (s)
Upper	300	280	150
Lower	1.0	0.5	1.0

Table 4

Shahab-3 performance parameters and dimensions.

	Shahab-3
Cycle	GG
Fuel	RP1
Oxidizer	RFNA
Range (km)	1350–1500
Payload (kg)	760–1158
GLOW (kg)	15,832–16,250
Diameter (m)	1.32–1.35
Length (m)	15.852–16.0
F (kN)	260–262
$I_{sp}$ (s)	230
$\Delta t$ (s)	110
$p_c$ (bar)	55

Table 5

Fitted data for RFNA and RP1; $p_e$ 0.7 bar.

Variable	a	b	c
MR	2.92	0.62	−0.14
$c^{\ast }$ (m/s)	1493.41	54.97	−40.97
$C_F$	0.77	0.51	−1.07
$T_c$ (K)	2364.71	448.96	−199.64
$\gamma$	1.19	0.03	0.01
M	20.87	2.10	−0.69
$ϵ$	0.37	0.41	0.77

Table 6

Comparison of Shahab-3 performance parameters and dimensions with the research model.

Parameters	Shahab-3	Calculated	Optimized	Optimized with Similar ${\mathit{p}}_{\mathit{c}}$
GLOW (kg)	15,832–16,250	15,650	14,312	15,973
Diameter (m)	1.32–1.35	1.06	1.03	1.07
Length (m)	15.852–16.0	17.45	16.89	17.58
F (kN)	260–262	309	2.91	306
$I_{sp}$ (s)	230	258	266	255
$\Delta t$ (s)	110	110	109	112
$p_c$ (bar)	55	55	79	48
L/D	11.85	16.46	16.39	16.43
T/W	16.12	19.74	20.33	19.15

Table 7

GG Cycle: Values of Variables at max $I_{sp}$ and T/W for each technique.

TWISPOF	GG Cycle
	NM	SA	CSA	PSO	PIO	GA	GA-Swarm
F (kN)	256	251	264	271	257	257	259
$\Delta t$ (s)	101.1	97.3	103.1	95.3	101.5	101.6	102.0
$p_c$ (bar)	105.4	369.7	66.2	1149.5	95.2	94.6	84.2
$I_{sp}$ (s)	272.6	302.2	262.3	332.9	270.3	270.2	267.6
GLOW (tonnes)	11.6	10.1	12.6	9.8	11.8	11.8	12.1
T/W	2.24	2.52	2.13	2.80	2.21	2.21	2.19

Table 8

GG Cycle: Values of Variables at min GLOW for each technique.

GLOWOF	GG Cycle
	NM	SA	CSA	PSO	PIO	GA	GA-Swarm
F (kN)	257	269	259	255	259	257	259
$\Delta t$ (s)	101.4	104.4	101.9	100.9	102.0	101.5	101.8
$p_c$ (bar)	96.9	52.4	86.9	112.6	84.5	96.2	87.8
$I_{sp}$ (s)	270.7	257.1	268.3	274.1	267.7	270.5	268.5
GLOW (tonnes)	11.8	13.3	12.0	11.6	12.1	11.8	12.0
T/W	2.22	2.07	2.19	2.25	2.19	2.21	2.19

Table 9

GG Cycle: Values of Variables at max $I_{sp}$ for each technique.

${\mathit{I}}_{\mathbf{sp}}$OF	GG Cycle
	NM	SA	CSA	PSO	PIO	GA	GA-Swarm
F (kN)	262	257	261	251	257	258	259
$\Delta t$ (s)	102.7	101.5	102.4	99.1	101.5	101.8	101.9
$p_c$ (bar)	71.5	96.0	77.6	192.1	95.7	89.9	86.4
$I_{sp}$ (s)	264.0	270.5	265.8	286.3	270.4	269.0	268.2
GLOW (tonnes)	12.4	11.8	12.2	10.8	11.8	11.9	12.0
T/W	2.14	2.22	2.17	2.37	2.22	2.20	2.19

Table 10

GG Cycle: Values of Variables at max T/W for each technique.

T/WOF	GG Cycle
	NM	SA	CSA	PSO	PIO	GA	GA-Swarm
F (kN)	-	280	259	251	257	261	260
$\Delta t$ (s)	-	106.7	102.1	99.3	101.4	102.6	102.3
$p_c$ (bar)	-	36.5	84.3	175.2	97.4	75.4	79.7
$I_{sp}$ (s)	-	249.1	267.6	284.2	270.8	265.1	266.4
GLOW (tonnes)	-	14.4	12.1	10.9	11.8	12.3	12.2
T/W	-	1.98	2.19	2.35	2.22	2.16	2.17

Table 11

SC Cycle: Values of Variables at max $I_{sp}$ and T/W for each technique.

TWISPOF	SC Cycle
	NM	SA	CSA	PSO	PIO	GA	GA-Swarm
F (kN)	267	273	271	270	269	277	271
$\Delta t$ (s)	99.0	99.7	99.8	99.5	99.4	100.8	99.7
$p_c$ (bar)	110.7	92.1	90.2	97.1	99.9	70.8	91.9
$I_{sp}$ (s)	273.7	269.6	269.1	270.8	271.4	263.7	269.6
GLOW (tonnes)	11.8	12.3	12.3	12.2	12.1	12.8	12.2
T/W	2.29	2.26	2.25	2.26	2.27	2.19	2.26

Table 12

SC Cycle: Values of Variables at min GLOW for each technique.

GLOWOF	SC Cycle
	NM	SA	CSA	PSO	PIO	GA	GA-Swarm
F (kN)	265	270	264	267.3	269	264	274
$\Delta t$ (s)	98.7	99.6	102.6	98.4	99.5	102.7	100.2
$p_c$ (bar)	122.3	94.2	74.2	133.1	97.0	73.7	81.9
$I_{sp}$ (s)	275.9	270.1	264.8	277.9	270.7	264.6	267.0
GLOW (tonnes)	11.7	12.2	12.5	11.5	12.1	12.5	12.5
T/W	2.32	2.59	2.16	2.34	2.27	2.16	2.23

Table 13

SC Cycle: Values of Variables at max $I_{sp}$ for each technique.

${\mathit{I}}_{\mathbf{sp}}$OF	SC Cycle
	NM	SA	CSA	PSO	PIO	GA	GA-Swarm
F (kN)	269	271	271	269	270	270	271
$\Delta t$ (s)	99.4	99.8	99.7	99.0	99.4	99.6	99.7
$p_c$ (bar)	99.7	90.4	92.6	110.7	99.0	95.0	93.0
$I_{sp}$ (s)	271.3	269.2	269.7	273.7	271.2	270.3	269.8
GLOW (tonnes)	12.1	12.2	12.2	11.9	12.1	12.1	12.2
T/W	2.27	2.25	2.26	2.30	2.27	2.26	2.26

Table 14

Values of Variables at max T/W for each technique.

T/WOF	SC Cycle
	NM	SA	CSA	PSO	PIO	GA	GA-Swarm
F (kN)	271	269	272	267	269	269	272
$\Delta t$ (s)	99.6	99.6	99.9	99.0	99.5	99.5	100.0
$p_c$ (bar)	95.5	96.1	88.5	111.4	98.2	97.7	86.3
$I_{sp}$ (s)	270.4	270.5	268.7	273.8	271.0	270.9	268.1
GLOW (tonnes)	12.2	12.1	12.3	11.8	12.1	12.1	12.4
T/W	2.26	2.26	2.25	2.30	2.27	2.27	2.24

Appendix A. CPropEP Graphical Calculations for Shahab-3

The equations for calculating the thermodynamic properties for RFNA-RP1 at the optimum MR are shown in Figure A1.

Figure A1. Curve-Fitting Results for RFNA-RP1 at 0.7 bar [Forumla omitted. See PDF.].

References

1. Hakanen, J.; Allmendinger, R. Multiobjective optimization and decision making in engineering sciences. Optim. Eng.; 2021; 22, pp. 1031-1037. [DOI: https://dx.doi.org/10.1007/s11081-021-09627-x]

2. Sutton, G.; Biblarz, O. Rocket Propulsion Elements; 9th ed. John Wiley & Sons: New York, NY, USA, 2017.

3. Tizon, J.M.; Roman, A. A Mass Model for Liquid Propellant Rocket Engines. Proceedings of the 53rd AIAA/SAE/ASEE Joint Propulsion Conference AIAA Propulsion and Energy Forum; Atanta, GA, USA, 10–12 July 2017; [DOI: https://dx.doi.org/10.2514/6.2017-5010]

4. Hammond, W. Design Methodologies for Space Transportation Systems; American Institute of Aeronautics and Astronautics: Reston, VA, USA, 2001; 2,72

5. Dong, S.; Fang, Y.; Hu, W.; Shen, X.; Chen, G. Research on Disruptive Design Pre-identification in the Preliminary Design Phase of Aero Engine Flow Pass Multidisciplinary Optimization. Chin. J. Aeronaut.; 2022; [DOI: https://dx.doi.org/10.1016/j.cja.2022.01.010]

6. Jain, P.; Patial, S.; Rawat, P. Review and Anatomical Study of Launch Vehicles in use and under-development. Proceedings of the 8th European Conference for Aeronautics and Space Sciences; Madrid, Spain, 1–4 July 2019; [DOI: https://dx.doi.org/10.13009/EUCASS2019-977]

7. Zandbergen, B.T.C. Simple mass and size estimation relationships of pump fed rocket engines for launch vehicle conceptual design. Proceedings of the 6th European Conference for Aeronautics and Space Sciences (EUCASS); Krakow, Poland, 29 June–3 July 2015; Available online: https://www.researchgate.net/publication/279711349_Simple_mass_and_size_estimation_relationships_of_pump_fed_rocket_engines_for_launch_vehicle_conceptual_design (accessed on 16 September 2022).

8. Rohrschneider, R.R. Development of a Mass Estimating Relationship Database for Launch Vehicle Conceptual Design. Master’s Thesis; Georgia Institute of Technology: Atlanta, GA, USA, 26 April 2002.

9. Wildvank, R. Launch Vehicle Structural Mass Prediction Model. Master’s Thesis; TU Delft Aerospace Engineering: Delft, The Netherland, 14 June 2018; Available online: http://resolver.tudelft.nl/uuid:bb943ec1-1a2a-462d-873e-05338b29235f (accessed on 16 September 2022).

10. Cai, G.; Fang, J.; Zheng, Y.; Tong, X.; Chen, J.; Wang, J. Optimization of System Parameters for Liquid Rocket Engines with Gas-Generator Cycles. J. Propuls. Power; 2010; 26, pp. 113-119. [DOI: https://dx.doi.org/10.2514/1.40649]

11. Castellini, F. Multidisciplinary Design Optimization for Expendable Launch Vehicles. Ph.D. Thesis; Politecnico Di Milano: Milan, Italy, 2012.

12. Montazeri, M.; Ebrahimi, R. Multidisciplinary optimization of a pump–fed system in a cryogenic LPE using a systematic approach based on genetic algorithm. Aerosp. Sci. Technol.; 2016; 49, pp. 185-196. [DOI: https://dx.doi.org/10.1016/j.ast.2015.12.003]

13. Saqlain, A.; He, L. Optimization and Sizing for Propulsion System of Liquid Rocket Using Genetic Algorithm. Chin. J. Aeronaut.; 2007; 20, pp. 40-46. [DOI: https://dx.doi.org/10.1016/S1000-9361(07)60005-3]

14. Bayley, D.; Hartfield, R.; Burkhalter, J.; Jenkins, R. Design Optimization of a Space Launch Vehicle Using a Genetic Algorithm. J. Spacecr. Rocket.; 2008; 45, pp. 733-740. [DOI: https://dx.doi.org/10.2514/1.35318]

15. Mota, F.A.d.S. Modeling And Simulation Of Launch Vehicles Using Object-Oriented Programming. Ph.D. Thesis; Instituto Nacional de Pesquisas Espaciais—INPE: Sao Jose dos Campos, Brazil, 2015.

16. Zhou, C.; Yu, N.; Wang, J.; Cai, G. Starting and regulating characteristics of electric pump feed system for LRE under different schemes. Appl. Sci.; 2022; 12, 6441. [DOI: https://dx.doi.org/10.3390/app12136441]

17. Mirshams, M.; Naseh, H.; Taei, H.; Fazeley, H. Liquid propellant engine conceptual design by using a fuzzy-multi-objective genetic algorithm (MOGA) optimization method. Proc. Inst. Mech. Eng. Part G J. Aerosp. Eng.; 2014; 228, pp. 2587-2603. [DOI: https://dx.doi.org/10.1177/0954410014521390]

18. Kosugi, Y.; Oyama, A.; Fujii, K.; Kanazaki, M. Multidisciplinary and multi-objective design exploration methodology for conceptual design of a hybrid rocket. Proceedings of the Infotech@Aerospace 2011; St Louis, MO, USA, 29–31 March 2011.; [DOI: https://dx.doi.org/10.2514/6.2011-1634]

19. Orgeira-Crespo, P.; Rey, G.; Ulloa, C.; Garcia-Luis, U.; Rouco, P.; Aguado-Agelet, F. Optimization of the Conceptual Design of a Multistage Rocket Launcher. Aerospace; 2022; 9, 286. [DOI: https://dx.doi.org/10.3390/aerospace9060286]

20. Castellini, F.; Lavagna, M.R. Comparative analysis of global techniques for performance and design optimization of launchers. J. Spacecr. Rocket.; 2012; 49, pp. 274-285. [DOI: https://dx.doi.org/10.2514/1.51749]

21. Mack, Y.; Haftka, R.; Griffin, L.; Snellgrove, L.; Dorney, D.; Huber, F.; Shyy, W. Radial turbine preliminary aerodynamic design optimization for Expander Cycle Liquid Rocket Engine. Proceedings of the 42nd AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit; Sacramento, CA, USA, 9–12 July 2006; [DOI: https://dx.doi.org/10.2514/6.2006-5046]

22. Cui, S.; Li, J.; Zhang, S.; Bai, X.; Sui, D. Overall Parameters Design of Air-Launched Rockets Using Surrogate Based Optimization Method. Aerospace; 2022; 9, 15. [DOI: https://dx.doi.org/10.3390/aerospace9010015]

23. Dresia, K.; Deeken, J.; Oschwald, M.; Waxenegger-Wilfing, G. Machine Learning Methods for the Design and Operation of Liquid Rocket Engines—Research Activities at the DLR Institute of Space Propulsion. Proceedings of the Space Propulsion 2020+1 Conference (Digital Conference); Estoril, Spain, 17–19 March 2021; [DOI: https://dx.doi.org/10.48550/arXiv.2102.07109]

24. Yu, J.; Li, Z.; Jia, L.; Zhang, Y. Switching Neural Network Control for Underactuated Spacecraft Formation Reconfiguration in Elliptic Orbits. Appl. Sci.; 2022; 2022 12, 5792. [DOI: https://dx.doi.org/10.3390/app12125792]

25. Sohail, M.U.; Hamdani, H.R.; Islam, A.; Parvez, K.; Khan, A.M.; Allauddin, U.; Khurram, M.; Elahi, H. Prediction of Non-Uniform Distorted Flows, Effects on Transonic Compressor Using CFD, Regression Analysis and Artificial Neural Networks. Appl. Sci.; 2021; 2021 11, 3706. [DOI: https://dx.doi.org/10.3390/app11083706]

26. Feng, Q.; Hai, X.; Sun, B.; Ren, Y.; Wang, Z.; Yang, D.; Hu, Y.; Feng, R. Resilience optimization for multi-UAV formation reconfiguration via enhanced pigeon-inspired optimization. Chin. J. Aeronaut.; 2022; 35, pp. 110-123. [DOI: https://dx.doi.org/10.1016/j.cja.2020.10.029]

27. Gordon, S.; Mcbride, B.J. Computer Program for Calculation of Complex Chemical Equilibrium Compositions and Applications: I. Analysis; NASA Reference Publication 1311; NASA: Cleveland, OH, USA, 1944.

28. Lefebvre, A. Welcome to the Rocketworkbench Project. Rocketworkbench Project. Available online: http://rocketworkbench.sourceforge.net/index.phtml (accessed on 4 August 2022).

29. Brown, E.D. An Introduction to PROPEP, A Propellant Evaluation Program for Personal Computers. J. Pyrotech.; 1995; 1, pp. 11-18. Available online: http://www.jpyro.co.uk/wp-content/uploads/j01_11_htfsr.pdf (accessed on 16 September 2022).

30. Huzel, D.; Huang, D. Modern Engineering for Design of Liquid-Propellant Rocket Engines; Knovel: Norwich, NY, USA, 1992.

31. Sadullah Khan, S.; Ahmad, M.A.; Qamar, I. Optimization of Pressure fed Liquid Propellant Rocket Engine using Various techniques and Objectives. J. Space Technol.; 2022; 12, pp. 39-48.

32. NOAA Document S/T 76-1562 U.S. Standard Atmosphere; National Oceanic and Atmospheric Administration, National Aeronautics and Space Administration: Washington, DC, USA, 1976.

33. Dresia, K.; Jentzsch, S.; Waxenegger-Wilfing, G.; Santos, H.R.; Deeken, J.; Oschwald, M.; Mota, F. Multidisciplinary Design Optimization of Reusable Launch Vehicles for Different Propellants and Objectives. J. Spacecr. Rocket.; 2021; 58, pp. 1017-1029. [DOI: https://dx.doi.org/10.2514/1.A34944]

34. Kuo, R.J.; Han, Y.S. A hybrid of genetic algorithm and particle swarm optimization for solving bi-level linear programming problem—A case study on Supply Chain Model. Appl. Math. Model.; 2011; 35, pp. 3905-3917. [DOI: https://dx.doi.org/10.1016/j.apm.2011.02.008]

35. Ali, A.F.; Tawhid, M.A. A hybrid particle swarm optimization and genetic algorithm with population partitioning for large scale optimization problems. Ain Shams Eng. J.; 2017; 8, pp. 191-206. [DOI: https://dx.doi.org/10.1016/j.asej.2016.07.008]

36. Esmin, A.A.A.; Lambert-Torres, G.; Alvarenga, G.B. Hybrid evolutionary algorithm based on PSO and GA mutation. Proceedings of the Sixth International Conference on Hybrid Intelligent Systems (HIS’06); Rio de Janeiro, Brazil, 13–15 December 2006; [DOI: https://dx.doi.org/10.1109/his.2006.264940]

37. Brugge, N. The North-Korean/Iranian Nodong-Shahab Missile Family. Available online: http://www.b14643.de/Spacerockets/Specials/Nodong/index.htm (accessed on 16 September 2022).