A Multidisciplinary Optimization Framework for Ecodesign of Reusable Microsatellite Launchers

Abstract

The commercial space launch sector is currently undergoing a significant shift, with increasing competition and demand for launch services, as well as growing concerns about the environmental impact of rocket launches. To address these challenges, within the New Space Portugal project scope, a multidisciplinary framework for designing and optimizing new launch vehicles is proposed. Creating a more resilient and responsible space industry can be achieved by combining technological innovation and environmental sustainability, as emphasized by the framework. The main scope of the framework was to couple all the disciplines relevant to the space vehicle design in a modular way. Significant emphasis was placed on the infusion of ecodesign principles, including Life Cycle Assessment (LCA) considerations. Optimization techniques were employed to enhance the design and help designers conduct trade-off studies. In general, this multidisciplinary framework aims to provide a comprehensive approach to designing next-generation launch vehicles that meet the demands of a rapidly changing market while also minimizing their environmental impact. A methodology that leverages the strengths of both genetic and gradient-based algorithms is employed for optimizations with the objectives of maximizing the apogee altitude and minimizing the Global Warming Potential (GWP). Despite only being tested at the moment for sounding rockets, the framework has demonstrated promising results. It has illuminated the potential of this approach, leading to the identification of three optimal designs: one for maximizing the apogee, another for minimizing GWP, and a compromise design that strikes a balance between the two objectives. The outcomes yielded a maximum apogee of $6.41$ km, a minimum GWP of $9.06$ kg ${\mathrm{CO}}_{2_{eq}}$, and a balanced compromise design featuring an apogee of $5.75$ km and a GWP of $25.64$ kg ${\mathrm{CO}}_{2_{eq}}$.

1. Introduction

The current article initiates an exploration into optimized design methodologies, following a structured path to redefine the optimization process. Sounding rockets serve as the initial applications in this endeavor. It initiates, in Section 1.1, by outlining the driving motivation for the ongoing work; then, a comprehensive literature review of the core topics is conducted in Section 2. A foundational understanding of the current framework follows in Section 3. Aligned with key objectives, this framework integrates Life Cycle Assessment (LCA), conducts extensive parametric studies, scrutinizes optimization solutions, employs various metrics for result comparison, and proposes refined solutions, reported in Section 4. Given the aspiration for the multidisciplinary optimization framework to advance to the design of orbital vehicles, potential future work is acknowledged and concluding remarks are provided in Section 5.

1.1. Motivation

Terms such as New Space or the democratization of space aptly characterize the global movement towards achieving quicker and more cost-effective access to space. This stands apart from traditional government-driven initiatives that primarily focus on security, politics, or scientific pursuits. The increased accessibility to space has broadened participation to a much larger pool of contributors leading to commercial competition over traditional state-sponsored launch systems, driving the market towards the use of small satellites and dedicated launch services [1]. This has led to the spread of new vehicles and new propellants such as methane and propane [2]. In parallel, concerns about anthropogenic climate change have led to further scrutiny and attention to the environmental impacts of all industries and to the increasing application of tools such as environmental LCA to minimize emissions. While this surge in activity could potentially exacerbate the degradation of the space environment, numerous New Space operators are actively engaged in efforts to minimize their impact [1]. However, New Space represents a significant departure from past practices and the conventional business approach may not be effective in this evolved landscape. Establishing new standards, space policies, and licensing approaches holds substantial potential and responsibility to shape the future of all missions [1,3,4].

In this context, led by GEOSAT, the New Space Portugal initiative [5] is a collaborative consortium comprising 38 entities with the collective goal of conceptualizing, developing, producing, and launching Earth Observation satellites: through the creation of the Atlantic Constellation. The project aims to develop the space sector in Portugal and assume a significant role in the New Space era.

Within the New Space Portugal initiative, universities play a crucial role in the ongoing efforts. Their focus lies on developing new research areas and curricula aimed at preparing future engineers to meet the specific needs of Portuguese companies. Additionally, the initiative involves the preliminary design of two types of vehicles: (i) a conventional reusable launcher for microsatellites, a concept vehicle capable of delivering potential future Portuguese satellites to low Earth orbit (LEO) with partial or complete reusability; and (ii) a reusable vehicle for multipurpose missions, a concept vehicle capable of changing orbit and performing on-orbit servicing with partial or complete reusability. This strategic collaboration between academia and industry contributes to the advancement and alignment of educational programs with the evolving demands of the space sector in Portugal.

2. Literature Review

Given the final objective of the ongoing work, to direct the design of the vehicle towards an innovative and competitive solution, an overview of existing and future microlaunchers should be conducted to identify possible gaps and alternatives.

2.1. Existing and Future Vehicles

This survey was limited to launchers with a payload capability of up to 500 kg into a LEO orbit, usually considered as the threshold to distinguish between small-lift launchers and medium- or heavy-lift launchers [6,7]. However, as the performance of a launcher is related to the desired orbit and comprises a range of masses, vehicles with slightly higher maximum capabilities have been included as well. Moreover, only commercial vehicles have been considered [8].

As of the writing of this article, only nine operational vehicles are dedicated to this range of masses worldwide, reported in Table 1. Notably, none of these vehicles are of European origin.

After the bankruptcy of Virgin Orbit, only three small-lift launchers are now available in the United States of America (USA) for commercial use. The Indian Small Satellite Launch Vehicle (SSLV), although developed by the Indian Space Research Organization, is intended to be a commercial vehicle and has already launched an American payload [9] to reinforce this intention. The remaining operational vehicles are Chinese and, although all of them are for commercial purposes, only two of them have dedicated websites, only one in English (Ceres1 [10]).

Table 1. Existing commercial microlaunchers.

Name	Company	Country	Payload [kg]	Type	Reusable
Pegasus XL* [11]	Northrop Grumman	USA	<460	ALTO	Aircraft
Minotaur 1 * [12]	Northrop Grumman	USA	<584	VT	No
Electron+ [13]	Rocket Lab	USA/New Zealand	<310	VT	No
SSLV [9]	NSIL	India	<500	VT	No
Ceres1 [10]	Galactic Energy	China	<400	VT	No
Hyperbola1 [14]	iSpace	China	<300	VT	No
Jielong 1 [15]	China Rocket Co., Ltd.	China	<150	VT	No
Kuaizhou-1A [16]	ExPace	China	<400	VT	No
OS-M * [17]	OneSpace Technology	China	<553	VT	No

* Although considered operational, no launches occurred in 2023, and no scheduled ones could be identified in 2024.

Table 1. Existing commercial microlaunchers.

Name	Company	Country	Payload [kg]	Type	Reusable
Pegasus XL* [11]	Northrop Grumman	USA	<460	ALTO	Aircraft
Minotaur 1 * [12]	Northrop Grumman	USA	<584	VT	No
Electron+ [13]	Rocket Lab	USA/New Zealand	<310	VT	No
SSLV [9]	NSIL	India	<500	VT	No
Ceres1 [10]	Galactic Energy	China	<400	VT	No
Hyperbola1 [14]	iSpace	China	<300	VT	No
Jielong 1 [15]	China Rocket Co., Ltd.	China	<150	VT	No
Kuaizhou-1A [16]	ExPace	China	<400	VT	No
OS-M * [17]	OneSpace Technology	China	<553	VT	No

* Although considered operational, no launches occurred in 2023, and no scheduled ones could be identified in 2024.

Among the existing launchers, Pegasus XL stands out as the only one that does not utilize vertical take-off (VT). Additionally, as it uses the air-launch-to-orbit (ALTO) method, it is the sole launcher characterized by partial reusability, as the carrier aircraft is a conventional aircraft capable of landing after the launch.

It is clear how reusability is still not widely considered in the design of launch vehicles, even for small masses. This is destined to change, as can be seen from the projects for future microlaunchers, where an increasing number of vehicles strive to include these considerations in their design.

The survey on future projects was focused only on Europe-based vehicles. Although Europe has, at the moment, no operational microlauncher, several vehicles are being developed by various private companies, reported in

Table 2. European projects for microlaunchers.

Name	Company	Country	Payload [kg]	Status	Type	Reusable
Miura1 [18]	PLD Space	Spain	<150	FT	SSTO VTSD	Yes
Bloostar [19]	Zero2infinity	Spain	<140	FT	Rockoon	None
MK3 [20]	Dawn Aerospace	Netherlands New Zealand	<5	Stage 2 FT	HTHL	Yes
NorthStar1 [21]	Nammo	Norway	<10	Stage 1 FT	VT	No info
Skyrora XL [22]	Skyrora	UK/Ukraine	<335	Prototype	VTVL	Stages 1–2
Prime [23]	Orbex	UK	<180	Prototype	VT	Partially *
Black Arrow 2 [24]	Black Arrow Space Technologies	UK	<500	GT	VT	No info
SL1 [25]	HyImpulse	Germany	<500	GT	VT	No
OB1-MK1 [26]	HyPrSpae	France	<200	GT	VT	No
Zephyr [27]	Latitude	France	<100	GT	VT	No
Colibri [28]	B2space	UK	<150	GT	Rockoon	Stage 1
EOS [29]	Sidereus	Italy	<13	GT	SSTO VTHL	Yes

* What does not burn up harmlessly in the atmosphere will be recovered and re-used.

.

As can be seen, half of these designs include partial or complete reusability through vertical or horizontal landing (VL or HL) or splashdown (SD).

Moreover, a much broader range of types of vehicles is present, including single-stage-to-orbit (SSTO), rockoons (combination of rockets and balloons), and horizontal take-off (HT) (airplane-like design).

However, at present, only two of them have been through a flight test (FT), both with limited apogees for safety reasons [30,31]. Some companies are proceeding through flight tests of intermediate stages, while others already have full prototypes prepared for upcoming flight tests. The remaining ones are still undergoing ground testing (GT) of some components.

It is evident that, despite the current limited number of vehicles dedicated to this range of masses, there is a notable increase in development activities, attributed to the growing demand associated with the advent of the New Space era. Furthermore, a noteworthy trend observed in the projects under development is the increasing consideration of reusability in launcher design. This underscores the contemporary integration of sustainability as an integral component of space vehicle design practices.

2.2. Life Cycle Assessment (LCA)

The urgency to address climate change has underscored the need to identify and measure environmental impacts, which, in turn, needs environmental legislation. To quantitatively analyze environmental effects, the one technique most commonly used known as LCA is employed in the formulation, execution, and oversight of global environmental policies. LCA is a well-established approach for evaluating the ecological footprint of a product or process throughout its entire life cycle. The international standards for LCA, described in ISO 14040-14044, are maintained by the International Organization for Standardization (ISO) [32]. LCA also serves to prevent the inadvertent trade-off of one environmental issue for another (burden shifting), which involves transferring impacts from one phase of the life cycle of the product to another, from one geographical area to another, from one generation to the next, or between various categories of environmental effects. It carries substantial weight in decision-making, emphasizing also the potential need to incorporate the Multi-Criteria Decision Analysis (MCDA) technique during the design phase [33].

As emphasized at COP21 (Conference of the Parties) in Paris, 26 of the 50 crucial variables for assessing the climate of Earth are observed through satellite monitoring [34]. As of the writing of this article, COP28 is currently underway, where it is supposed to emphasize the potential role of space-based earth observation data in advancing the objectives of the Paris Agreement. This includes its crucial contribution to the Global Stocktake process, serving as a key tool for monitoring climate-related changes, comprehending the impacts of climate change, and evaluating the effectiveness of mitigation and adaptation measures [35]. Recognizing the orbital space of Earth as an additional ecosystem, deserving meticulous care and thorough regulations underscores the importance of implementing LCA in the Space sector [36,37].

The scarcity of information, particularly regarding altitude-dependent effects of non-${\mathrm{CO}}_2$ emissions and ozone-destroying compounds from launchers and satellites, adds to the complexity [38]. Additionally, environmental impacts, such as emissions into the upper atmosphere, are rarely addressed in traditional LCAs, with the magnitude of their effect being unknown [38,39]. In fact, National Oceanic and Atmospheric Administration (NOAA) scientists investigating the stratosphere have identified measurements revealing that around 10% of aerosol particles in that part of the atmosphere consist of aluminum and other metals generated from the combustion of satellites and rocket stages during reentry. The consequences of this concentration of metallic content on the characteristics of stratospheric aerosol remain uncertain [40,41]. In recent years, the European Space Agency (ESA) has integrated environmental requirements into various projects, including Copernicus, Ariane 6, and Galileo. These initiatives, along with other technological advancements, contribute to the creation of multiple LCA datasets, including the creation of their own. This database emerges as a crucial resource for the space sector to meet environmental standards, also to be shared with all European stakeholders, facilitating their participation in future LCA studies [42,43]. Moreover, the Strathclyde Space Systems Database (SSSD) is another useful tool, featuring over 250 validated space-specific Life Cycle Inventory (LCI) datasets and designed expressly for the streamlined LCA of early space mission concepts [33]. Also, even if not space-oriented, various software options, both open-source [44,45,46] and not [47], are available for conducting LCA. ESA also has a portal dedicated to the user community of its Space Debris Software [48], which acts as a gateway for applying for software licenses and accessing the software.

The ESA LCA working group also has published a handbook with the objective of setting uniform methodological guidelines for the accurate execution of LCA specifically tailored to the space industry [42].

In response to the sustainability challenges in the space sector, where centralized regulatory authority is lacking, the Space Sustainability Rating (SSR) concept emerged from discussions within the World Economic Forum Global Future Council on Space Technologies. The SSR aims to assign a score that assesses the sustainability of a mission, particularly concerning debris management and adherence to international guidelines [49,50,51].

Efforts are underway to develop simulation tools for assessing the direct emissions from rocket launches [52]. With the space sector envisioning increasingly ambitious plans for the future, it becomes crucial to proactively assess the environmental impact of those rocket launches on Earth. This includes work that demonstrates the consequences of space activities, analyzing data, and forecasting that the impacts may become more meaningful with the scaling up of those activities [53]. Other research endeavors focus on estimating the ${\mathrm{CO}}_2$ intensity of the space sector and examining the environmental impact of emissions from space launches [34,54].

Launch vehicles stand as the sole anthropogenic entities emitting directly into every atmospheric layer, and reusability could potentially introduce added burdens. Proper comprehension of these dynamics is essential to make sustainable design decisions in space transportation. Additional research [33,55] showcases the potential optimization of space missions for greater sustainability, addressing adverse environmental, social, and economic impacts early in the design process. Indeed, the studies, adopting a Life Cycle Engineering (LCE) approach, provide valuable insights into the life cycle sustainability impacts of SmallSat missions.

In conclusion, heightened awareness of climate change across all sectors can foster constructive dialogue within the space sector on this subject, facilitating the proper dissemination of methodology and findings. The need for clarity and norms cannot be overstated to increase the reliability of the presented results and avoid the risk of “greenwashing”.

2.3. Multidisciplinary Design Optimization (MDO)

The major challenge of a launch vehicle design process is given by its intrinsic multidisciplinary nature: as each discipline will influence the overall system and the other disciplines, it is impossible to optimize each of them alone.

An improvement to the traditional sequential workflow can be achieved through the Multidisciplinary Design Optimization (MDO) approach. This method employs numerical optimization techniques to identify the best trade-offs between the disciplines, thus aiding the designer in their choices.

The general MDO problem can be stated as follows [56]:

$$\begin{array}{ccc}\hfill \mathrm{minimize}& f(x;u;z)\hfill & \\ \hfill \mathrm{by}\mathrm{varying}& x_i\hfill & \hfill i=1,\dots ,n_x\\ \hfill \mathrm{subject}\mathrm{to}& g_j(x;u)\le 0\hfill & \hfill j=1,\dots ,n_g\\ \hfill & h_k(x;u)=0\hfill & \hfill k=1,\dots ,n_{\mathrm{h}}\\ \hfill & {\underline{x}}_i\le x_i\le {\overline{x}}_i\hfill \\ \hfill \mathrm{while}\mathrm{solving}& r_l(u;x)=0\hfill & \hfill l=1,\dots ,n_u\\ \hfill \mathrm{by}\mathrm{varying}& u_l\hfill \end{array}$$

This can be read as: minimize the objective function $f(x;u;z)$ by varying the design variables (DVs) $x_i$ within their lower ${\underline{\mathrm{x}}}_i$ and upper ${\overline{\mathrm{x}}}_i$ bounds subject to inequality and equality constraints $g_j(x;u)$ and $h_k(x;u)$ while solving the state equations $r_l(u;x)$ by varying the state variables $u_l$.

The design of a launch vehicle can be therefore imposed as the minimization of one or more objectives, like the mass or the emissions, given a set of varying DVs and a set of constraints. Within the framework, different functions can be introduced to solve the state equations relative to the subsystems and the overall system. The output will be a set of DVs that ensure a feasible design that minimizes the imposed objective, finding a set of state variables that solve the state equations for that set of DVs.

The design of a launch vehicle is a multidimensional, constrained problem. The optimization can be single or multi-objective: in the first case, the output will be unique and will be the one minimizing the objective function; if two or more objectives are imposed, there might not be a single optimal design, and the concept of Pareto optimality should be used [56].

The algorithm used to solve the optimization problem can be either heuristic or deterministic. On the one hand, random guesses are given to the optimizer, and their outputs are used to educate the next generation of guesses by sorting, combining, and mutating them, thus proceeding towards the optimum. On the other hand, the deterministic approach uses first or second-order information to identify the minimization direction [56].

Lastly, the architecture of the framework can be either monolithic or distributed. Monolithic architectures solve the optimization problem at the same level of the disciplines, while the distributed ones decompose the problem into sub-optimization problems for each discipline, then combined at the system level by another optimizer [56].

As it is not the scope of this document, the advantages and disadvantages of each choice will not be discussed. Details about the utilized optimizer will be explained and presented later in the document.

The importance of incorporating LCA into the design phase has already been emphasized, which resolves into its inclusion within the MDO framework, because the MDO approach seamlessly integrates into the framework alongside other disciplines, providing a direct and cohesive inclusion. The absence of integrated LCA may introduce biases due to the overlooking of various variables affecting environmental aspects. The integration of LCA within MDO addresses this issue by providing access to all relevant variables, preventing omissions from the MDO model [57,58]. Nevertheless, challenges exist, particularly in handling the complexity of data and variables. The data-intensive nature of LCA increases interdependencies within the model, requiring careful consideration [57].

Primarily, the main contribution to the state-of-the-art of the ongoing work is the development of a multiobjective, multidisciplinary framework for the design of sounding rockets—which will ultimately lead to microsatellite launchers. It incorporates multiple disciplines, most notably LCA, and combines gradient-based and genetic algorithms. With a focus on emphasizing ecodesign principles, this methodology seeks to offer a comprehensive and adaptable approach to optimizing complex aerospace systems.

3. Methods and Tools

The current work is being developed based on a monolithic framework initially designed for sounding rockets, entirely built in MATLAB^®. The original purpose of its creation was the participation in student competitions, such as the European Rocketry Challenge (EuRoC) or Spaceport America^® Cup. The work now involves upgrading the existing modules and introducing new ones to enhance the capability of the framework, enabling it to model orbital vehicles [59,60,61].

The code is able to perform optimization using genetic and gradient-based algorithms.

In the ongoing work, the Non-Dominated Sorting Genetic Algorithm II (NSGA-II) proposed by Deb et al. [62] as an evolution of [63] and documented in [64] is employed as the genetic algorithm. Its design revolves around two key principles: non-dominant sorting to organize solutions based on their priority and crowding distance to identify the optimal solutions within the same front. The optimization process stops when the number of generations reaches the specified maximum $n_{gen}$.

The MATLAB^® built-in algorithm $fmincon$ was instead chosen as the gradient-based optimizer as it is able to deal with constrained nonlinear multivariable functions [65].

3.1. MDO Framework

The framework comprises four integral building blocks: propulsion, mass and sizing, aerodynamics, and LCA. Detailed explanations of each discipline will be provided in the subsequent section.

The framework has undergone a restructuring process from its initial formulation [60]. This transformation aimed to eliminate backward loops of state variables by establishing a sequential flow from one discipline to another. Consequently, this design modification eradicates the necessity for enforcing consistency at the system level or through an iterative solver.

The flowchart of the conceptual methodology framework is illustrated in Figure 1.

To enhance clarity regarding the DVs and inputs within the framework, a schematic representation of the rocket is presented in Figure 2.

The DVs and their associated bounds are sourced from the original formulation by Yamada [60] and are predominantly derived from reference values obtained from existing sounding rockets as enumerated in

Table 3. Original DVs and bounds. $V_{OT}$ is the volume of the oxidizer tank; $C_{inj}$ is the effective area of the injector; $L_{FG}$, $D_{FG}$ are the length and diameter of the fuel grain, respectively; $D_{Th}$, $A_{Th}$ are the diameter and area of the throat; $A_{exit}$ is the exit area; and $B_{fin}$, $b_{fin}$, $S_{fin}$, and ${\theta }_{fin}$ are the root chord, tip chord, span, and sweep angle of the fins, respectively.

Number	Design Variables (DVs)	Lower Bound (LB)	Upper Bound (UB)	Units
1	$V_{OT}$	$4.0\times {10}^{-3}$	$1.4\times {10}^{-2}$	${\mathrm{m}}^3$
2	$C_{inj}$	$1.0\times {10}^{-5}$	$4.0\times {10}^{-5}$	${\mathrm{m}}^2$
3	$L_{FG}$	$2.0\times {10}^{-1}$	$7.0\times {10}^{-1}$	m
4	$D_{FG}$	$3.0\times {10}^{-2}$	$8.0\times {10}^{-2}$	m
5	$D_{Th}$	$2.0\times {10}^{-2}$	$5.0\times {10}^{-2}$	m
6	$A_{exit}/A_{Th}$	$3.0\times {10}^{-0}$	$10\times {10}^{-0}$	-
7	$B_{fin}$	$1.0\times {10}^{-4}$	$3.0\times {10}^{-1}$	m
8	$b_{fin}/B_{fin}$	$1.0\times {10}^{-4}$	$9.0\times {10}^{-1}$	-
9	$S_{fin}$	$1.0\times {10}^{-4}$	$4.0\times {10}^{-1}$	${\mathrm{m}}^2$
10	${\theta }_{fin}$	$0.0\times {10}^{-0}$	$6.0\times {10}^1$	°

.

As will be explained later, the bounds of the DVs play a major role in the behavior of the optimizer. However, releasing them without worrying about the physical implications of such choices might result in unrealistic simulations.

Scalability issues might be related to the loss of validity of taken assumptions within the different functions (e.g., linear regressions or empirical correlations), the need to adapt fixed inputs (e.g., thickness of the walls of the tanks), imposition of constraints related to unconsidered phenomena (e.g., extinguishing of the flame due to excessive liquid drops from the injector), or simply manufacturing issues.

Moving to each discipline, there is a requirement for several fixed inputs. These inputs should be subject to modification in tandem with any changes in the definition of disciplines or DVs.

General inputs relate to the payload (mass and length) and the launch site (ambient temperature and pressure, altitude, and launch rail length and angle). Those are treated as requirements and not changed throughout the whole process.

Within the propulsion block, inputs are associated with the properties of the oxidizer, the fuel, and their interactions. Additionally, initial conditions of the oxidizer tank (OT) are provided, along with efficiency coefficients for the combustion chamber (CC) and the nozzle.

Given the absence of structural analysis within the framework, the sizing block utilizes fixed thicknesses and density to dimension the rocket components. Lengths of the nose cone and nozzle, as well as the number of fins, are also fixed.

A correction factor for the drag coefficient can be introduced to the aerodynamic function for sensitivity studies related to inaccuracies in its calculations.

In the LCA module, SimaPro^® is employed to assess the impact of fuel and oxidizer production, along with combustion products. A dedicated database has been established to cover various oxidizer-to-fuel (OF) ratios and pressures.

Finally, options for propulsion and trajectory simulations, such as time steps, propagators, and tolerances, are customizable.

The state variables are computed in each discipline to be used in the following functions and are summarized in

Table 4. State variables. All symbols are reported at the end of the document.

Variable	Type	Use
$m_{o{x}_{init}}$	Scalar	LCA
$R_{C{C}_{in}}$	Scalar	Sizing, CG
$OF$ ratio	Vector	LCA
$P_{CC}$	Vector	LCA
$m_{out}$	Vector	LCA
$m_{liq}$	Vector	CG
$m_{vap}$	Vector	CG
$d_{liq}$	Vector	CG
$d_{vap}$	Vector	CG
$m_{fuel}$	Vector	CG, LCA
Thrust curve	Table	Trajectory
$D_e$	Scalar	CD, Trajectory
$D_{No{z}_{ex}}$	Scalar	CD, Trajectory
$H_{No{z}_{ex}}$	Scalar	CD
$L_{tube}$	Scalar	CD, Re
$m_{rocket}$	Vector	Trajectory

.

As no variable computed in the aerodynamic function is used in the LCA module, the state variables come exclusively from the propulsion and sizing blocks.

To have a feasible design, several constraints have been imposed, adapted from the original formulation of Yamada [60], and reported in

Table 5. Constraints. All symbols are reported at the end of the document.

Number	Constraint	Purpose
1	$P_{CC}-$ 0.8 ($P_{OT}$ − $\Delta P_{inj}$)	Avoid backflow
2	$G_{max}-$ 650	Combustion stability
3	$fd-$50	Avoid excessively slender bodies
4	5$-fd$	Avoid excessively bulk bodies
5	SM$-5$	Avoid overstabilization issues
6	1−SM	Impose stable configuration
7	$3050-$Apogee	Competition objective
8	$Ac{c}_{max}-$150	Avoid excessive loads

.

Among these constraints, the first and sixth are deemed necessary. The first constraint prevents backflow from the CC to the OT, ensuring a unidirectional flow. The sixth constraint is essential for maintaining a stable rocket configuration.

Constraints two, five, and eight could potentially be relaxed, but careful consideration of their implications is advised. According to Fraters and Cervone [66], higher values of oxidizer mass flux may be feasible, but additional studies are required to precisely assess their impact on engine performance. Regarding the upper limit of the static margin (SM) (constraint five), conventional practice uses a reference value of two to mitigate the excessive influence of lateral winds [67]. However, given the absence of a proper structural design and a more detailed definition of other vehicle components, the decision was made to allow the optimizer to explore configurations within a higher limit, to be fine-tuned in subsequent stages. The eighth constraint is enforced to prevent excessive loads and is anticipated to be replaced with failure modes from the structural analysis module once implemented.

Constraints three and four are imposed to preclude unreasonable configurations, such as extremely slender or bulky vehicles, based on reference values from Niskanen [67].

Lastly, the apogee constraint is linked to the objective of the competition and has been retained to discourage configurations with excessively low apogees.

The initial objectives outlined in Yamada [60] included the take-off mass and the specific impulse. Subsequently, with the incorporation of the LCA module, a third objective was introduced, according to the desired metric. Comprehensive details regarding the various metrics will be expounded upon in the subsequent section.

Additional objectives under consideration encompass the apogee altitude for a given payload mass or, conversely, the payload mass for a specified minimum apogee.

3.2. Disciplines

The conceptual design of the rocket consists of interconnected modules for propulsion, mass and sizing, aerodynamics, and LCA. The vehicle is a single-stage hybrid rocket, as illustrated in Figure 2.

3.2.1. Propulsion Module

The propulsion module comprises four models representing the OT, the injector, the CC, and the nozzle, as illustrated in Figure 3. It simulates the mass and energy fluxes during the combustion of a hybrid motor. The choice of hybrid propulsion systems was based on the numerous advantages they offer when compared to both solid and liquid rocket motors [68,69]. The hybrid system used in this study combines a solid fuel (paraffin (${\mathrm{C}}_n$${\mathrm{H}}_{2n+2}$)) with a liquid oxidizer (nitrous oxide (${\mathrm{N}}_2$O)). Nitrous oxide is a commonly selected oxidizer for small hybrid rocket engines, primarily due to its self-pressurizing property and overall safety.

• Oxidizer Tank (OT)—As mentioned earlier, the OT is self-pressurized, maintaining equilibrium between its liquid and gaseous states. Therefore, the thermal saturation properties are employed in this context [70]. The emptying process is calculated in two distinct stages: the first involves both liquid and gas phases, and the second occurs when the tank depletes its liquid ${\mathrm{N}}_2$O. Given that the assumption of an ideal gas for ${\mathrm{N}}_2$O is not valid near the saturation state, then the Span–Wagner equation of state is applied to simulate the behavior of the real gas [71,72].
• Injector—The injector connects the OT and the CC. In this system, the estimation of mass flow is made with an incompressible single-phase model (SPI) [71].
• Combustion Chamber (CC)—The oxidizer is then injected into the CC, initiating a chemical reaction that results in the complete combustion of the fuel. Steady flow and constant properties are assumed. Hybrid rocket grains typically undergo radial burning, according to the shape of their port. In the current work, a cylindrical grain with a circular port and no aluminum droppage was considered. There is a dependence between the fuel regression rate and the oxidizer mass flow, considering also empirical/experimental coefficients adjusted individually for each fuel-oxidizer pair. The regression rate is used to calculate the fuel mass flow. The pressure in the chamber is calculated using relations for the upstream section of a converging–diverging nozzle. In order to calculate the remaining mixture properties, the Chemical Equilibrium with Applications (CEA) software [73] is used given the varying chamber pressure and OF ratio.
• Nozzle—The nozzle is assumed to be ideal isentropic [74], without shock waves, and with the steady axial flow of an ideal gas in a homogeneous fluid. The flow through the supersonic nozzle is also assumed to be choked. These assumptions help to calculate downstream properties. Additionally, a correction factor is applied to consider losses not taken into account with the ideal model when determining the thrust. Moreover, once the OT runs out of liquid oxidizer, overexpansion at the exit of the nozzle is considered.

3.2.2. Mass and Sizing Module

The Mass and Sizing module delineates the geometries and masses of all rocket components. Additionally, it encompasses stability analyses, not including any structural considerations, which will be addressed in subsequent developments.

Utilizing geometric relationships, the module determines the lengths of the OT and CC. These lengths are then combined with those of the nose cone and nozzle, enabling the estimation of the overall rocket length through linear regressions based on existing sounding rockets. The computed overall length is then used to approximate the dimensions of the avionics and recovery modules using similar regressions. Lastly, the tube is extended to cover the nozzle if its exit area is smaller than that of the tube, or until it makes contact with it.

Similar regression models are employed to estimate the masses of these components. Geometric rules are applied to compute the masses of the cone and tube skins, fins, and CC and OT vessels, computed assuming a cylindrical shape.

To determine the center of gravity (CG) height for cylinders, cones, and cone frustums, geometrical relations are used to find the distances from the cone tip for each element. The fins are divided into simpler geometries, and the same relations are applied. Note that the fins are not controlled.

Subsequently, the total CG is computed as a vector over time, considering the oxidizer phases and propellant depletion. Barrowman equation [75] is utilized to compute the center of pressure (CP) and the SM of the vehicle, using the tube diameter as a reference.

If the calculated SM falls outside the predefined bounds in Table 5, a ballast mass is introduced inside the nose cone or around the portion of the nozzle enclosed by the tube.

In instances where compliance with both constraints is unattainable due to limited available volume on both sides, the loop is terminated without adding mass. The optimizer proceeds normally, and the constraints will be violated.

Upon successful adjustment, the updated mass, CG, and SM values are transferred to the Aerodynamics function.

3.2.3. Aerodynamics

As the stability calculations have already been carried on, the Aerodynamics module assumes responsibility for computing aerodynamic coefficients and trajectory propagation.

Three propagators are at disposal: two fixed-step Runge–Kutta algorithms of 4th and 8th order [76,77] and an embedded 5th-order method featuring variable step sizes [78].

A simplified one degree of freedom (DoF) vertical model is utilized to mitigate the computational cost of the process.

During each computational step, atmospheric properties (pressure, Mach number $Ma$, and Reynolds number $Re$) are determined based on altitude. These properties are then utilized to calculate thrust using the thrust curve and the drag coefficient, where the latter depends on the $Ma$ and $Re$. The drag coefficient is updated at each step, while the lift coefficient remains null.

The computation of the drag coefficient draws on the principles outlined in the open-source code OpenRocket [67]. This approach decomposes the drag coefficient into three components: skin friction, pressure, and base coefficients.

Skin friction drag arises from the viscous flow of air around the rocket. The skin friction coefficient is computed with consideration to the $Re$ and skin roughness (defaulted to 60 $\mathsf{\mu }$m) and adjusted for compressibility effects. Subsequently, it is used to compute the drag coefficient, scaling it appropriately by the reference areas [67].

Pressure drag results from air being forced around the rocket. In the subsonic regime, it can be correlated with the nose joint angle. However, a more detailed model is required for supersonic conditions due to the presence of shockwaves. To simplify computations, a semi-empirical model is employed [67].

Base drag emerges from a low-pressure area at the base of the rocket and is computed through an empirical function of the $Ma$ [67].

The official documentation notes a subsonic regime error in the drag coefficient compared to real data below 10%, with no verification possible in supersonic conditions, where a dedicated model should be implemented [67].

Propagation concludes upon the impact of the rocket with the ground, as no recovery or landing options have been considered at this stage.

3.2.4. LCA Module

To systematically conduct an LCA, a four-step procedure is followed, encompassing goal and scope definition, inventory analysis, impact assessment, and interpretation. The first step of the LCA is crucial to provide a roadmap for the ensuing phases. This methodology was adopted to understand the environmental impacts of the designs resulting from the adoption of NSGA-II and subsequent post-processing outcomes. Thus, the selected functional unit, representing the impact resulting from each design option, serves as a foundational metric in this exploration. Notably, in this initial phase of the ongoing work, the system boundaries include only the use phase of each design, providing a focused perspective. For the inventory analysis, data are collected on the activities identified in the system boundary selected, namely fuel production, oxidizer production, and emissions associated with the combustion process. The combustion products data, essential for a comprehensive analysis, is obtained using CEA software [73]. This software analyzes combustion and rocketry problems for various OF ratios and CC pressures. Additionally, the Ecoinvent 3.8 database is also integrated to model diverse processes, ensuring a comprehensive inventory. Moreover, for the environmental impact assessment, this study resorts to SimaPro^® software version 9.5.0.2 by adopting the ReCiPe 2016 Midpoint (Hierarchist) method [47]. This approach, compared, for example, with an endpoint approach, provides a more transparent and detailed evaluation of different impact categories despite introducing a potential complexity during the analysis.

As already mentioned, since space-oriented databases and environmental impact assessment methods are not available for the present study, an example from a launch vehicle life cycle is presented in

Table A1. Example of midpoint and flow indicators that can be used in the LCA of a launcher, adapted from ESA LCA Working Group [42].

Impact Indicator	Unit	Description
Global Warming Potential (GWP)	kg ${\mathrm{CO}}_{2_{eq}}$	the extent of radiative forcing attributed to greenhouse gas (GHG) emissions
Ozone Depletion Potential	kg CFC$-{11}_{eq}$	depletion of stratospheric ozone resulting from the emission of ozone-depleting substances (ODS)
Human Toxicity Potential	CTUh	impact of emitting toxic substances on human health
Photochemical Ozone Formation Potential	kg NMVOC	creation of tropospheric ozone resulting from the interaction of primary pollutants, such as volatile organic compounds (VOCs), with solar radiation in the atmosphere
Fossil Resource Depletion Potential	GJ Fossil	the reduction in the availability of non-living natural resources, such as crude oil
Gross Water Consumption Potential	${\mathrm{m}}^3$	quantity of water extracted from lakes, rivers, oceans, and groundwater for industrial processes
Mass Disposed in the O cean	kg	the total mass of stages that have been disposed of in the ocean
${\mathrm{Al}}_2$${\mathrm{O}}_3$ Emissions in Air	kg	release of alumina emissions into the air during a launch event
Mass Left in Space	kg	total mass of space hardware that remains in orbit at the conclusion of the mission

, adapted from ESA LCA Working Group [42], showcasing possible environmental midpoint indicators and flow indicators.

Up to this point, the framework has exclusively used the Global Warming Potential (GWP) metric (see Table A1), focusing specifically on the fuel. However, it is imperative to explore and analyze additional metrics to emphasize the comprehensive nature of LCA. It is also crucial to extend this analysis to cover all other components, such as structures, not just limited to fuel.

The function within this discipline computes kg of ${\mathrm{CO}}_{2_{eq}}$ for each time step and integrates it throughout the entire propulsion.

4. Results

Before the following investigations, a sensitivity analysis was conducted regarding the time steps of the propulsion and aerodynamics propagators, leading to an average simulation time on the order of magnitude of seconds for a 13th Gen Intel^® Core™ i9 13900KF 3.00 GHz processor.

To gain confidence in the current problem, various investigations were conducted to scrutinize the hyperparameters of the NSGA-II algorithm, namely $n_{pop}$, $n_{gen}$, and ${\%}_{mut}$, as well as the influence of constraints and the role of the bounds of the DVs.

To mitigate the risk of over-constraining the problem, during the sensitivity analysis focusing on the NSGA-II algorithm, Constraints 3, 4, 5, and 8 were considered non-mandatory for this phase.

Preliminary findings revealed that the optimizer tended to converge towards solutions located at the upper or lower bounds (UB or LB) of the DVs, limiting exploration towards these extremes. Consequently, a decision was made to expand the bounds, aiming to understand the behavior of the optimizer without the constraints imposed by these limits. It is important to note that the expansion of bounds was purely numerical, executed without considering the physical implications, and undertaken bilaterally solely to assess their influence on the final results.

The sensitivity analysis began by varying one parameter at a time while keeping the remaining parameters at their reference values (baseline). A simulation with simultaneously increased $n_{pop}$ and $n_{gen}$ was performed, as well as a simulation with the baseline values and original bounds. Moreover, three additional simulations were performed to clarify the role of randomness along the process.

Following the assessment of hyperparameter influence, the role of constraints was investigated by comparing the unconstrained version with the constrained one. This analysis also considered the relations between original and enlarged bounds and their correlation with the selected objective.

It is crucial to acknowledge that expanding the bounds and eliminating constraints could potentially result in designs that are unfeasible or unrealistic. However, it is imperative to emphasize that these investigations were undertaken solely to comprehend their impact on the optimization process. The primary objective was not to generate practically implementable designs, but rather to gain insights into the behavior of the optimizer.

Subsequently, reverting to the original bounds and constraints, a final optimization was executed with hyperparameters deemed suitable for the problem. The optimization results underwent post-processing through a gradient-based optimizer, culminating in the presentation of the final Pareto front and three proposed designs.

4.1. NSGA-II Sensitivity Analysis

The code initiates by creating an initial set of solutions, forming a population. Following an evaluation process, half of the population undergoes recombination to produce offspring, and these offspring are subject to random mutation based on the specified mutation probability. Each cycle of this procedure is termed a generation.

In the initial generation, the entire population undergoes evaluation, while in each subsequent generation, only the newly generated half is assessed. The total number of function calls $FC$ can thus be computed as follows:

$$FC=\frac{n_{pop}}{2}(n_{gen}+1).$$

(1)

The Latin Hypercube Sampling (LHS) method [56] has been employed to generate the initial population of designs. The parameters for the conducted simulations are summarized in

Table 6. Sensitivity analysis simulations.

Simulation	1	2	3	4	5	6 *	7	8	9
$n_{pop}$	200	200	200	800	800	200	800	500	500
$n_{gen}$	200	200	800	200	800	200	800	500	500
${\%}_{mut}$	30	10	30	30	30	30	10	30	10
$FC$	2.01 $\times {10}^4$	2.01 $\times {10}^4$	8.01 $\times {10}^4$	8.04 $\times {10}^4$	1.60 $\times {10}^5$	2.01 $\times {10}^4$	1.60 $\times {10}^5$	6.28 $\times {10}^4$	6.28 $\times {10}^4$

* Same as Simulation 1 with original bounds (Table 3).

.

To conduct the initial parametric study, only one parameter is altered at a time, while the remaining parameters are maintained at their reference value, i.e., the baseline.

Three objectives have been considered: wet mass (mass of the rocket with all the propellant), specific impulse $I_{sp}$, and GWP measured in kg of emissions of ${\mathrm{CO}}_{2_{eq}}$.

Furthermore, simulations were executed using a population size $n_{pop}$ of 100. Nevertheless, the results obtained suggested that, given the nature of the specific problem and the number of DVs, utilizing such a modest population size would result in an inefficient optimization process, yielding a Pareto front constrained to two or three points. Consequently, it would not provide sufficient information to draw meaningful conclusions or make accurate comparisons.

The initial simulations unveiled a linear relationship between mass and emissions, demonstrating a trend where minimizing one objective correlated with the reduction of the other. This distinct relationship is illustrated in detail in Figure 4, depicting Simulations 1 and 6.

Hence, while all simulations encompassed three objectives, the reported figures will focus solely on the projections concerning $I_{sp}$ and GWP.

Many considerations can be derived from the charts in Figure 5:

• A lower mutation probability leads to a diminished exploration of the design space, while excessive mutation may hinder the exploitation of promising results by overly altering them. These results do not provide a conclusive determination regarding whether a 30% mutation probability represents an optimal compromise.
• Augmenting the population size results in broader coverage of the design space, enhancing exploration. However, without a corresponding increase in the number of generations, exploitation may be compromised.
• Increasing the number of generations enhances exploitation, driving all designs towards the Pareto front and improving its extremes.

Fine-tuning of these parameters is essential to identify the optimal compromise between exploration and exploitation, aligning with the specific problem under study. For improved comparison, the Pareto fronts of Simulations 1–6 have been presented in Figure 6.

The Pareto fronts from Simulations 1 to 5 exhibit similar results and trends. As anticipated, capturing the precise influence of the mutation probability proves challenging. Simulation 3 showcases designs with a favorable compromise between the two objectives, benefiting from increased exploration of the design space while maintaining similar extremes. In contrast, Simulation 4 highlights a design with significantly higher $I_{sp}$ due to improved exploitation of promising designs. Simulation 5 illustrates how simultaneous increases in $n_{pop}$ and $n_{gen}$ can synergistically contribute to achieving a superior optimal front.

Notably, optimization with the original bounds yields a highly restricted Pareto front, yet overall superior results compared to other simulations. This can be attributed to the smaller design space being more accessible for exploration and the original bounds being optimized for the initial framework definition. Any expansion in the design space would necessitate a corresponding increase in population and generations to adequately cover it, as well as an adjustment of the fixed inputs.

It is noteworthy that among the eight configurations comprising the optimal front of Simulation 6, only one featured all DVs within their specified bounds. The remaining configurations exhibited two to five variables reaching the boundary limits.

Finally, for a more in-depth understanding of the role of the mutation probability, Simulations 1 & 2 were compared to 5 & 7, and 8 & 9, to investigate whether its effect became more pronounced with increases in population number and generations.

Analysis of the optimization process, depicted in Figure 7, reveals an intricate relationship between mutation rates, population size, and the number of generations. Notably, for a population of 200 individuals and a maximum of 200 generations, the comparison between mutation rates of 30% and 10% presents an initial ambiguity in identifying the superior rate.

Observing the trend across varying population sizes and generations, it can be said that a pattern emerges, indicating the influence of mutation rates on optimization outcomes. With an increase in both population size and generations, the advantage of a higher mutation rate becomes progressively evident. This observation aligns with the ability of the genetic algorithm to explore a broader solution space and mitigate premature convergence to suboptimal solutions.

Thus, considering the potential for greater exploration and the avoidance of premature convergence, the mutation rate of 30% was selected as the preferred parameter for subsequent optimizations. This choice aims to leverage the inherent capability of the algorithm to diversify and enhance the search for more optimal solutions within the given problem space.

4.2. Unconstrained Problem

To enhance the performances of the optimizer, the objectives were reduced from three to two. Given the linear relationship between mass and emissions, only emissions were considered. Additionally, $I_{sp}$ was replaced by apogee altitude, aligning with the primary objective of designing an orbital vehicle.

To gain deeper insights into the problem, all constraints except the apogee one were removed, and the bounds were further expanded to attempt to obtain a well-defined Pareto front. The number of generations was set equal to the number of designs in a population to prevent under-exploitation, and it was increased for several cases reported in Figure 8. The mutation rate was fixed at 30%.

Additionally, in this version, the ballast loop was enhanced by incorporating an additional mass at the bottom, to set the maximum SM to 5.

Clearly, with an increase in the number of populations and generations, the optimality front becomes progressively more pushed towards the origin, thus leading to better intermediate results.

In a final analysis, the original bounds were reinstated, and an optimization with $n_{pop}$ = $n_{gen}$ = 200 was conducted, then compared to the one with expanded bounds. In this scenario, the obtained fronts exhibit similar trends and values. The front with larger bounds demonstrates superior extremes and compromise designs overall (refer to Figure 9), albeit with fewer points due to the increased challenge of exploring a larger design space.

Given that the primary changes from the previous version were the alteration in objectives and the removal of constraints, a similar comparison was undertaken using two hybrid versions: an unconstrained one with $I_{sp}$ as one of the objectives and a constrained one retaining the apogee altitude.

Upon comparing the unconstrained and constrained versions, it becomes evident that the imposition of constraints introduces a noticeable offset in the identified fronts (Figure 9). Furthermore, this effect is more pronounced in the version with larger bounds, once again indicating poorer performance when utilizing the larger bounds. It is hypothesized that, as the constraints and functions have been specifically defined for sounding rockets, the application of larger bounds increases the likelihood of failures throughout the optimization process, ultimately resulting in inferior performance. This phenomenon is not observed in the unconstrained version, where the optimizer has the freedom to explore the entire design space, facilitating the discovery of superior solutions.

When employing the unconstrained version with $I_{sp}$ as one of the objectives, the version with original bounds once again demonstrated superior performance, as depicted in Figure 10. It is noteworthy that, in this case, the change in $I_{sp}$ along the fronts was relatively small. However, for the original bounds, consistently lower emissions were observed to achieve the same specific impulse.

This behavior presents challenges in understanding, given the intricate computation of $I_{sp}$ that involves various aspects of the design. Further investigation is required to elucidate this phenomenon.

4.3. Optimization and Post-Processing

After gaining insights into the problem, the following decisions were made for the final optimization:

• The constraints in Table 5 were reinstated to ensure a feasible design.
• The original bounds in Table 3 were retained, considering the expectation of better performance when these constraints are applied and acknowledging the absence of investigation into potential physical implications associated with their increase.
• A population and generation size of $n_{pop}$ = $n_{gen}$ = 1000 was chosen to enhance exploration of the design space, exploit promising designs, and obtain a well-defined Pareto front.
• A mutation probability of ${\%}_{mut}$ = 30 was selected, deemed suitable for populations of this order of magnitude.

To validate the simulations, the time steps of the propulsion and aerodynamics propagators were decreased by an order of magnitude, and the points on the obtained Pareto front were reanalyzed. All solutions remained feasible, with maximum relative differences for the apogee altitude and GWP at 0.6224% and 0.0693%, respectively.

To maximize the exploitation of potential results, as well as mitigating the influence of randomness related to genetic algorithms, the results of the NSGA-II algorithm underwent post-processing using the gradient-based optimizer $fmincon$, from MATLAB^®, through single-objective optimizations that used the points on the front as starting points with the scope of obtaining a better Pareto front. Given that the designs are already within the feasible region, the sequential quadratic programming algorithm ($sqp$) was employed, with $StepTolerance$ set to $1\times {10}^{-8}$ and $ConstraintTolerance$ set to $1\times {10}^{-3}$. The selection of the specified values for $StepTolerance$ and $ConstraintTolerance$ resulted from a sensitivity analysis. This analysis demonstrated that there were no changes in the objective function with a $StepTolerance$ smaller than $1\times {10}^{-8}$. Additionally, for 28 out of 30 designs, no changes were observed due to the $ConstraintTolerance$ since the starting points were already within the feasible region. The remaining two designs experienced some iterations in the unfeasible region but with constraint violations below $1\times {10}^{-6}$. A change in this factor was therefore considered unnecessary.

To accommodate the two objectives in a single-objective process, both were normalized, and the objective function was defined as the weighted average of the two.

For normalization, the extremes of the front were extracted, and the respective values were used as minima and maxima using the following relation:

$$Ob{j}^{\prime }=\frac{Obj-min\left(Obj\right)}{max\left(Obj\right)-min\left(Obj\right)}$$

(2)

assigning 0 to the minimum and 1 to the maximum.

However, if the objective function was defined as the sum of the two normalized objectives, the optimizer would tend to optimize both to the same extent, which is not conducive to creating a Pareto front. To address this, defining N as the number of designs on the front and sorting them in descending order according to one of the objectives, the following relation was applied:

$$f_n=\frac{N-n}{N-1}Ob{j}_1^{\prime }+\frac{n-1}{N-1}Ob{j}_2^{\prime }$$

(3)

resulting in $f_1=Ob{j}_1^{\prime }$, $f_N=Ob{j}_2^{\prime }$, and weighted objectives in between.

The optimal compromise was identified as the design with the highest value of $\sqrt{Ob{j}_1^{\prime 2}+Ob{j}_2^{\prime 2}}$. As the best compromise from the multi-objective optimization was not located at the position $(N+1)/2$, its post-processing through the $fmincon$ was not achieved with equal weights for the two objectives. The optimization was repeated with different weightings to assess their effects on the final result. Specifically, the optimization was rerun with equal weights, and the weights were also slightly adjusted to better define that portion of the Pareto front. The following weight pairs were utilized: 0.35–0.65, 0.45–0.55, 0.5–0.5, 0.55–0.45, and 0.65–0.35. The results of these optimizations were incorporated into the post-processed Pareto front. Additionally, optimizations with weight pairs 1–0 and 0-1 were also conducted from that same starting point.

Lastly, it was observed that the Pareto front obtained in this iteration closely resembled, in terms of trend and values, the one obtained in the previous version, achieved using $n_{pop}$ = $n_{gen}$ = 200. A final optimization with these parameters was consequently executed, and the extremes and optimal compromise were once again optimized through the $fmincon$.

All the comparisons mentioned have been reported in Figure 11.

As evident, the application of $fmincon$ on the results of the NSGA-II optimization results in a notable enhancement of the Pareto front. This improvement can be attributed to the distinct functioning principles of the two solvers: while the NSGA-II excels in exploring the design space, its exploitation of promising designs is constrained to recombination and random mutations. In contrast, a gradient-based solver can discern the direction of minimization, thereby refining the results.

The improvement is more conspicuous in the middle part of the front, with limited impact at the extremes. This is presumed to be linked to the utilized bounds for the DVs, and further clarification will be provided later.

Regarding the optimization of the compromise design, the obtained results were very similar to each other, and though differences are reported in the figure, they are not visually discernible. Furthermore, the best compromise was achieved through the post-processing of another design. Optimization to maximize apogee from this design yielded unsatisfactory results, whereas minimizing emissions found a design with the lowest emissions and a superior apogee altitude compared to the others in that region.

While this might seem peculiar initially, it can be explained by the propensity of gradient-based algorithms to become entrenched in specific regions of the design space. This results in significant dependency on the starting point, and therefore in the possibility of obtaining a superior final design by altering the initial set of DVs.

Finally, optimizations based on the $n_{pop}$ = $n_{gen}$ = 200 simulation extremes yielded suboptimal results for apogee and compromise, while emissions minimization mirrored the outcomes of the other optimal results.

In conclusion, it is posited that the optimizer excels in optimizing for emissions due to its reliance solely on the propulsion module in the developed framework. However, achieving optimal apogee involves the interplay of all disciplines, making it more intricate to correlate with the DVs.

Please note that, although the randomness involved in the optimization process affects its outcome, as noted in Figure 6, Figure 7 and Figure 8 performed at the beginning of the study, it did not affect the overall trends of the results. Moreover, as can be seen in Figure 11, where a gradient-based optimizer was used to improve the Pareto front resulting from the final NSGA-II optimization, the randomness is attenuated for such a large number of individuals and generations.

For the extreme configurations and the best compromise, various data are presented in

Table 7. Characteristics of the optimal designs.

Design (NSGA-II)	Max. Apogee	Compromise	Min. GWP
Apogee altitude (km)	6.3337	5.3333	3.2901
Wet mass (kg)	32.0284	24.2428	13.5379
Dry mass (kg)	23.6885	18.5558	11.1655
Diameter (m)	0.1111	0.0964	0.0732
Length (m)	4.0535	3.8368	3.1990
$I_{sp}$ (s)	181.6208	182.8082	183.3461
GWP (kg ${\mathbf{CO}}_{2_{eq}}$)	31.8320	21.7254	9.0900
Ozone formation (HH) (kg ${\mathrm{NO}}_x$eq.)	0.0532	0.0363	0.0152
Terrestrial acidification (kg ${\mathrm{SO}}_2$eq.)	0.2930	0.1901	0.0692
Ozone formation (TE) (kg ${\mathrm{NO}}_x$eq.)	0.0543	0.0371	0.0155
Design (fmincon)	Max. Apogee	Compromise	Min. GWP
Apogee altitude (km)	6.4118	5.7479	3.0500
Wet mass (kg)	31.4131	27.3219	14.5904
Dry mass (kg)	23.1603	20.6237	12.1428
Diameter (m)	0.1091	0.1008	0.0766
Length (m)	4.1568	4.0429	3.0602
$I_{sp}$ (s)	183.0729	182.0266	178.7986
GWP (kg ${\mathbf{CO}}_{2_{eq}}$)	31.7850	25.6419	9.0563
Ozone formation (HH) (kg ${\mathrm{NO}}_x$eq.)	0.0529	0.0428	0.0154
Terrestrial acidification (kg ${\mathrm{SO}}_2$eq.)	0.3263	0.2456	0.0417
Ozone formation (TE) (kg ${\mathrm{NO}}_x$eq.)	0.0540	0.0437	0.0157

, comparing the results of NSGA-II and the outcomes of the post-processing.

The configurations of these points are reported in Figure 12. Indeed, it is notable that all the identified configurations exhibit a high degree of slenderness, characterized by fineness ratios around 40. While these values adhere to the imposed constraints, it is crucial to acknowledge that the absence of a comprehensive structural analysis renders it challenging to assess the risk of bending or buckling. If structural considerations were incorporated, it may potentially result in excessively heavy structures for these designs, consequently diminishing overall performance.

An additional analysis was conducted along the Pareto front to observe the behavior of the DVs within their boundaries. The order of the indices along the Pareto front corresponds to minimizing emissions to maximizing apogee: DVs 6 ($A_{exit}/A_{Th}$), 7($B_{fin}$), 8 ($b_{fin}/B_{fin}$), and 10 (${\theta }_{fin}$), which exhibit unrestricted movement within their bounds, will not be presented.

Furthermore, DVs 2 ($C_{inj}$), 3 ($L_{FG}$), and 5 ($D_{Th}$) consistently remained at their LBs along the Pareto front, irrespective of the objectives, whether minimizing the Global Warming Potential (GWP) or maximizing the apogee.

While it may seem intuitive that increasing the length of the fuel grain would lead to a higher apogee, its interrelations with other DVs are more intricate than they appear. The increase in fuel grain length would elevate the fuel mass flow, consequently reducing the oxidizer-to-fuel ratio from optimal values. To counterbalance this effect, an increase in oxidizer mass flow would be required, achievable through a larger pressure drop or a bigger injector port. However, the DVs are interconnected; a larger injector port leads to a smaller pressure drop due to an increase in CC pressure, necessitating a larger throat diameter. Furthermore, an increase in oxidizer mass flow results in a greater mass flux through the grain port, already optimized to the maximum value allowed by the constraint for apogee maximization. To decrease it, a larger initial port is needed, once again altering the fuel mass flow and creating a loop that cannot be linearly handled.

The natural question arises: why does a configuration with low values for all these design variables (DVs 2 ($C_{inj}$), 3 ($L_{FG}$), and 5 ($D_{Th}$) on the lower bound and DV 4 ($D_{FG}$) limited to a maximum of 0.4) outperform a configuration where all these variables are increased proportionally to each other? It is believed that, as the code computes the external diameter of the rocket based on the total consumption of the fuel grain given the OT volume, the bounds of the latter would lead to an even thinner configuration if a longer grain is used.

This hypothesis has been substantiated through a gradient-based optimization initiated from the design that maximizes the apogee, wherein an increase in the upper bound of the OT volume led to higher values of DVs 2-5 and an elevated apogee altitude. However, as anticipated, any modification in the bounds of the DVs must be meticulously analyzed to understand the physical implications of such choices.

From Figure 13, it can be seen that the volume of the oxidizer tank $V_{OT}$ changes almost linearly with the index. This correlation aligns with the logic that minimizing the volume of the oxidizer tank directly reduces the overall oxidizer mass used during operation. This reduction in oxidizer mass directly impacts emissions, especially because the emissions profile considered is solely related to the combustion process.

DV 9 ($S_{fin}$) exhibits an increasing pattern, with peaks when DV 7 ($B_{fin}$) shows its lowest points (see Figure 14a). This suggests potential trade-offs between aerodynamic and stability considerations, where changes in the root chord might affect the aerodynamic characteristics of the fins, potentially leading to higher aspect ratios with smaller root chords. As these variables related to the fins are interconnected, further analysis to try to correlate and optimize them might be conducted.

Lastly, from Figure 14b, it can be observed that both DVs 4 ($D_{FG}$) and 9 ($S_{fin}$) appear to be confined within a normalized range of approximately 0.1 to 0.4. This suggests that the boundaries set for these variables might be too expansive for the specific problem at hand. Refining these bounds could potentially lead to a more focused and effective exploration of the design space.

From Figure 15 and

Table 8. Comparison with respect to the bounds for the design on the Pareto fronts of NSGA-II and $fmincon$.

Design Variable (DV)	LB (NSGA-II)	UB (NSGA-II)	LB (fmincon)	UB (fmincon)
1	7%	13%	3%	0
2	100%	0	3%	0
3	100%	0	3%	0
4	0	0	0	0
5	100%	0	0	0
6	0	3%	0	0
7	0	40%	0	0
8	30%	23%	0	0
9	0	0	0	0
10	0	0	0	3%

, it is evident that after applying the $fmincon$ on the three results, most of the DVs move away from the bounds (see Table 8). This observation suggests a potential limitation of the NSGA-II. For instance, if the crossover operation predominantly occurs between solutions situated near the bounds due to Pareto dominance or diversity preservation mechanisms, it might perpetuate values at the boundaries across generations, restricting exploration away from those bounds. Additionally, it is noteworthy to mention that the first three DVs remain on the LB in the first solution (minimizing GWP) of the $fmincon$ as well.

To emphasize the comprehensiveness of LCA, it is crucial to consider multiple impact categories. This holistic approach ensures a more comprehensive understanding of environmental impacts, preventing the resolution of one environmental issue from inadvertently giving rise to others. For instance, although they all seem to have the same trend (see Figure 16), the results obtained in the environmental impact assessment presented in Table 7 show that while considering an extreme configuration of minimum GWP, despite a decrease of GWP with the outcome design of the post-processing, the impacts associated with the category Ozone Formation increase compared to the results of NSGA-II.

5. Conclusions

The development of a multidisciplinary framework for the conceptual design of microlaunchers has posed challenges due to the intricate nature of the required building blocks, their interconnectedness, and the scarcity of available data and literature on technologies, methods, and particularly sustainability aspects, so far largely neglected in these types of study.

Therefore, the focus of the work was redirected towards the optimization framework and its performance. This involved maintaining simplified discipline blocks to mitigate the computational cost of optimizations to glean insights into the overall optimization process.

5.1. Limitations

The current framework exhibits several limitations, including:

• Lack of structural considerations: The absence of structural analysis may lead to a potential underestimation of structural mass and volume, consequently resulting in overestimates of vehicle performance. The rationale behind this decision is primarily due to the extensive computational demands and the complexity involved in implementing the Finite Element Method (FEM).
• Oversimplified vertical one degree of freedom (DoF) trajectory module: The trajectory module lacks features such as rotations, wind models, more precise atmospheric models (time and location-dependent), perturbations from external bodies, and considerations for the launch site environment. This decision is justified by the absence of a gravity turn or control techniques to achieve an orbit around the planet. Introducing additional degrees of freedom (DoFs) at this stage would only contribute unnecessary complexities to the current modeling approach.
• Limited LCA analysis: The LCA module is confined to the propulsion module, excluding considerations related to other life cycle phases, such as the production and transformation of the materials used to produce the vehicles, different manufacturing technologies implemented, the transportation of components and materials, or even suitable end-of-life strategies. This limitation arises from the absence of space-oriented databases and environmental impact assessment methods currently accessible for the purpose of this study.
• Missing reliability analysis: As uncertainty is part of any design process, to achieve a realistic design, this should be accounted for in the optimization process. However, an uncertainty-based analysis is a very time-consuming process, and it has therefore not been considered in the current version of the framework [79,80]. This poses a limitation to the resulting Pareto fronts that will be assessed in future versions of the framework.

5.2. Future Improvements

Future improvements under consideration involve:

• Exploration of different fuels and oxidizer combinations, as well as the incorporation of various technologies such as external pressurization.
• Enhancement of propulsion module models: implementation of more precise models for the injector or the nozzle, and the incorporation of topology optimization for fuel grain geometry.
• Upgrade of correlations and reference values: given the objective of designing orbital vehicles, correlations and reference values based on sounding rockets should be upgraded.
• Implementation of a control module: introduction of a control module to perform maneuvers for reaching orbit and managing the eventual reentry phase.
• Consideration of aerodynamic accuracy: evaluation and potential refinement of aerodynamic considerations, particularly in supersonic and hypersonic regimes.
• Definition of different flight phases: clear definition of various flight phases in the mission profile, accompanied by the implementation of aerothermodynamics and thermal considerations for the reentry phase.

5.3. Achievements

Nevertheless, the optimization framework has proven itself to be a valuable tool for the conceptual design of launchers, yielding numerous useful results that illuminate the potential of MDO in the initial phases of space vehicle design.

The proposed methodology seeks to leverage the advantages of two categories of algorithms: genetic algorithms and gradient-based algorithms.

Genetic algorithms excel in mapping and exploring the design space, particularly for multi-objective processes. However, careful attention is required in defining the problem. A thorough investigation using the NSGA-II algorithm explored the impact of the number of individuals in a generation, the number of generations, and the mutation probability, as well as the relations between the definition of constraints and objective functions with the bounds of the DVs.

As expected, the observations indicate that an increase in the number of individuals in a population and the number of generations leads to a more favorable Pareto front. Notably, their simultaneous increase appears crucial for enhanced performance. Additionally, it was noted that a higher mutation probability is more suitable when both the number of individuals and generations are sufficiently high.

The investigation revealed that, for the specific problem under consideration, the optimizer’s tendency to become confined to the bounds of the DVs does not necessarily correlate with a performance improvement.

While this behavior was observed in an unconstrained problem with apogee altitude and Global Warming Potential (GWP) as objectives, it vanishes when constraints are reintroduced or when specific objectives, such as $I_{sp}$, are employed.

The impact of constraints can be attributed to their definition based on the original bounds, influencing the behavior of the optimizer. However, the relationship concerning $I_{sp}$ remains unexplained, and further analysis is required in this sense.

Once sufficient knowledge of the problem was acquired, the optimal set of options was employed to perform a final optimization, serving as a starting point for the gradient-based optimizer $fmincon$.

Gradient-based algorithms facilitate better exploitation due to their functioning principle based on identifying an optimization direction. However, this results in a strong dependence on the starting point and the inability to perform multi-objective optimizations.

Utilizing the results of NSGA-II as starting points led to a discernible improvement in the optimal Pareto front, highlighting the limited exploitation of genetic algorithms due to their intrinsic random nature. Nevertheless, using the $fmincon$ from suboptimal starting points resulted in suboptimal final designs, underscoring the limitations of gradient-based algorithms related to the starting point.

The proposed methodology effectively harnesses the strengths of both techniques, culminating in a well-defined Pareto front from which optimal designs for both objectives can be obtained.

Three distinct yet optimal designs were then especially and critically analyzed: one that maximized the apogee, another that minimizes the GWP, and, finally, the compromise between these two contradictory objectives. With this set of constraints, the maximum apogee obtained was $6.41$ km, the minimum GWP was $9.06$ kg ${\mathrm{CO}}_{2_{eq}}$, and a compromise was found with $5.75$ km of apogee and $25.64$ kg ${\mathrm{CO}}_{2_{eq}}$. Comparing the scenario of maximizing the apogee with the compromise scenario, a 10% decrease in the apogee corresponds to a 19% reduction in GWP. When comparing the scenario of maximizing the apogee with minimizing GWP, a 52% reduction in apogee is met with a significant 71% decrease in GWP.

A more profound understanding of the behavior of the design variables within their specified boundaries was obtained through an extensive trend analysis. This analysis unveiled the interrelations among these variables and elucidated how the utilized inputs and bounds restrict the exploration of potential optimal designs.

As environmental sustainability stands as a pivotal goal in this endeavor, the evaluation extended beyond the GWP to encompass additional metrics. This exploration emphasizes the imperative inclusion of a comprehensive LCA right from the preliminary stages of space vehicle design. Such an approach ensures a holistic consideration of sustainability factors, reinforcing the quest for an environmentally conscious design.