Open-access Staging Optimization of Satellite Launch Vehicles: Application to a Microsatellite Launcher

ABSTRACT

With modern technology, only a small percentage of the overall mass of a launch vehicle, usually less than 5%, can be used as payload mass. To achieve any acceptable performance, a key aspect that is exhaustively analyzed during vehicle design is the definition of the number and sizes of its stages. The problem of finding the best vehicle configuration for a given mission is called staging optimization. This work presents the formulation and a resolution method for the staging optimization of a launch vehicle. The problem is formulated as an optimal control problem (OCP) and, for its solution, a direct approach, employing a gradient-based algorithm as solver, is proposed. An application is performed for the Brazilian microsatellite launcher VLM-1. This vehicle has a notable characteristic that, in its nominal configuration, the first and second stages are very similar to each other (same solid rocket motor with the same quantity of propellant). Different vehicle configurations are investigated which result in an improvement in its performance. The obtained results show good agreement with those from a commercial optimization tool.

Keywords Staging; Launch vehicles; Optimal control; Nonlinear programming

INTRODUCTION

Space-based services are increasingly becoming more prevalent in our society’s daily activities, naturally accompanied by a growth in the so-called space industry (Perondi 2023; Pessoa Filho 2021). These activities, which include, among others, meteorological forecasts, navigation, telecommunications, and remote sensing, are based on data acquisition and transmission via artificial satellites orbiting the Earth. To fulfill the increasing number of launches of such systems, it is observed today a variety of developments of new launch vehicles, mainly in the category of small and micro launchers (Villas Bôas et al. 2020, 2023).

Launch vehicles are highly complex systems, with elevated development and operation costs. These systems are composed of propulsive units called stages, each of which comprises a quantity of propellant and all the subsystems necessary to transform the propellant into a thrust force via the exhaustion of hot gases produced by burning the propellant through the nozzle. The advantage of multistage launch vehicles resides in the possibility of jettisoning the unnecessary structural mass of stages whose propellant has already been burned. There are, basically, two types of staging architectures used in launch vehicles. In a tandem or serial architecture, each stage, apart from the first, is mounted on top of a previous stage, and its ignition can only occur after the previous stage has been jettisoned. In a parallel architecture, stages are mounted side by side, and both stages can burn simultaneously.

The configuration of a launch vehicle, in terms of the number of stages and the distribution of propellant mass among them, has a major impact on its performance, which can be measured in terms of its payload mass capacity. The optimization of this configuration, also called staging optimization, is a crucial aspect to be considered during the design process of the vehicle. The goal of staging optimization can be stated as finding the launch vehicle configuration with the minimum gross lift-off mass for a given payload mass and a target orbit or, equivalently, finding the vehicle configuration with the maximum payload capacity for a given total propellant mass and a target orbit. To solve this problem, different techniques can be applied, with varying simplifying hypotheses, which results in different types of optimization problems.

One possibility is to use analytical equations to estimate the velocity changes impressed on the payload by different sources: gain in velocity due to stage operation and losses due to aerodynamic drag, gravity, and steering (Adkins 1970; Malina and Summerfield 1947; Srivastava 1966; Tawakley 1967). This approach results in a multivariate function optimization problem, where the optimizable variables describe the vehicle configuration.

The analytical approach necessarily relies on simplifying hypotheses to allow the estimation of the velocity changes due to different effects. In an alternative approach, trajectory simulation might be used to provide these velocity gains/losses. A coupled methodology can be envisioned where the results from the staging optimization are used as input to a trajectory optimization process, which, in turn, provides estimations for all velocity gains/losses mechanisms (Cho et al. 2021; Jamilnia and Naghash 2012; Koch 2019). The solution process is, therefore, comprised of the alternation of a multivariate function optimization, related to the staging optimization, and the solution of an optimal control problem (OCP) for the trajectory optimization.

Finally, in a third approach, staging and trajectory optimization are formulated as a unified problem. In this approach, the parameters that describe the vehicle configuration (propellant and structural masses of the stages, motor burn times, etc.) are incorporated as optimizable parameters into a trajectory optimization problem (Jamilnia and Naghash 2012). The resulting problem belongs to the class of OCP.

In the present work, the staging optimization of a multistage launch vehicle is presented. The unified approach, where the staging optimization is accommodated inside the trajectory optimization problem, is used. Therefore, the results obtained in this work come from the formulation and resolution of an OCP.

Regarding the numerical solution of an OCP, two different approaches can be applied (Betts 1998; Rao 2010). In an indirect method, calculus of variations is used to determine the necessary conditions of optimality, which take the form of a multipoint boundary value problem. A numerical solution for this problem should, then, be obtained. In a direct method, the original problem is transformed via some sort of discretization, also called transcription (Hull 1997). This discretization leads to multivariate function approximations for the state and control variables, as well as for all the functions of the OCP, including the performance measure and the constraints. The resulting problem belongs to the class of nonlinear programming problems (NLPP).

To address the staging optimization problem in this work, the direct approach is employed. Two different transcriptions are used to transform the OCP into a NLPP: Hermite-Simpson collocation (Enright and Conway 1992; Hargraves and Paris 1987) and multiple-shooting (Bock and Plitt 1984). The resulting NLPP is then solved using an extension proposed by the authors of the sequential conjugate gradient-restoration algorithm (SGRA) (da Silveira and da Silva Fernandes 2023). As a case study, the proposed methodology is applied to the staging optimization of the Brazilian microsatellite launcher (veículo lançador de microssatélites) VLM-1.

A Mayer problem of optimal control

A typical problem of optimal control theory in Mayer’s form can be stated as (Betts 1998; Bryson and Ho 1975): consider the system described by the following differential Eq. 1:

x ˙ ( t ) = f ( t , x ( t ) , u ( t ) , p ) (1)

where x(t) denotes the n-vector of state variables, u(t) is a m-vector of control variables, and p is a k-vector of parameters. The problem is to find an admissible control u*(t) which transfers the system from an initial state x*(tI) to a final state x*(tF), where tI and tF are, respectively, the initial and final times of the trajectory, and minimizes the performance index in Eq. 2:

J = J t I , x t I , t F , x t F , p (2)

The initial state is defined by the following boundary conditions (Eq. 3):

b i t I , x t I = 0 (3)

with i = 1, ... , nI, nI < n + 1, and the final state is specified by Eq. 4:

b j t F , x t F = 0 (4)

with j = 1, ... , nF, nF < n + 1. Control and state variables may be subject to the path constraints in Eq. 5:

g L , l g l ( t , x ( t ) , u ( t ) , p ) g U , l (5)

with l = 1, ... , ng, and also subject to simple bounds, such as Eq. 6:

x L , r x r ( t ) x U , r (6)

where r = 1, ... , n and s = 1, ... , m. A control u(t) is admissible if it satisfies the previously constraints. The launch vehicle staging optimization problem will be formulated as a Mayer OCP, as described in the previous paragraphs.

Launch vehicle staging optimization problem

The launch vehicle considered here is composed of three propulsive stages, a fairing, and a payload. The stages are mounted using the serial architecture. The formulation of the staging optimization problem presented next, albeit general in several aspects, is adjusted for these characteristics. This formulation could easily be extended for other types of launchers.

The dynamic model employed here is an extension of the model used by the authors in a previous work to address the problem of a launch vehicle trajectory optimization (da Silveira and da Silva Fernandes 2023). The presentation below closely follows the cited work, with the exception of the mass modeling of the vehicle, which is proposed here to address the staging optimization. The position of the vehicle is described by three parameters defined with respect to the Earth’s center: radial distance R, longitude l, and geocentric latitude λ. Its inertial velocity is described by three cartesian components along the axes of a local frame L, whose origin is attached to the vehicle: the radial component VR along XL, the east component Vl along YL, and the north component Vλ along ZL (Fig. 1 shows the local frame L). According to these definitions, the vehicle equations of motion can be stated as Eq. 7 (da Silveira and da Silva Fernandes 2023; Weigel and Well 2000):

R ˙ = V R l ˙ = V l R cos   λ Ω E λ ˙ = V λ R V ˙ R = F x m + 1 R V λ 2 + V l 2 V ˙ l = F y m + V l R V λ tan   λ V R V ˙ λ = F z m 1 R V l 2 tan λ + V R V λ (7)

where m is the vehicle mass, ΩE is the Earth rotational speed, and Fx, Fy, and Fz are the total force components in L frame.

A body frame attached to the vehicle is defined with its XB axis along the vehicle longitudinal direction, pointing to the vehicle nose, and with YB and ZB axes perpendicular to XB (Fig. 1). The attitude of vehicle is given by the yaw angle Ψ and the pitch angle θ, which describe the angular displacement of frame B with respect to frame L (the vehicle is assumed to be symmetric with respect to its longitudinal axis, therefore the roll angle is not relevant in the present context). The rotation matrix that describes this attitude is (Eq. 8):

Figure 1
Local frame and body frame definitions.
L B L = R 2 π 2 + θ R 1 ( ψ ) (8)

where R1 and R2 are the rotation matrices representing a pure rotation around X and Y axes, respectively (Hughes 2004).

Along the flight, the vehicle is subjected to gravitational, thrust, and aerodynamic forces. The gravitational force is a function of the vehicle radial distance and geocentric latitude and can be written as Eq. 9:

W L = μ E m R 2 1 + 1.5 J 2 R E / R 2 1 3 sin 2   λ 0 3 J 2 R E / R 2 sin   λ cos   λ (9)

where μE is the Earth gravitational parameter, RE is the Earth mean equatorial radius, and J2 is the second-order zonal harmonic. Values for these parameters can be found in the literature (Vallado and Macclain 2013). The superscript in WL indicates that this vector is written in L frame.

The thrust force is easily written in the body frame as Eq. 10:

T B = T corr  0 0 (10)

where Tcorr = TvacpatmAe is the motor corrected thrust, which results from the motor vacuum thrust, Tvac, minus a correction term due to the atmospheric pressure patm acting on the nozzle exit area Ae. The thrust force is written in the local system L as Eq. 11:

T L = L B L T T B (11)

where the superscript T denotes the transpose of the matrix LBL.

Finally, the aerodynamic force is given by Eq. 12:

A B = 1 2 ρ atm  V A 2 C A S ref  0 0 (12)

where it is considered that only the axial component of this force has a meaningful impact on the trajectory. In Eq. 12, ρatm is the atmospheric density, VA is the vehicle speed with respect to the atmosphere (air speed), CA is the axial force coefficient, and Sref is the aerodynamic area of reference. If no wind is considered, the air speed is given by Eq. 13:

V A L = V R V l Ω E R cos   λ V λ (13)
Mass modeling

To complete the dynamic model of the vehicle, it is necessary to define how its mass varies along the trajectory. Since this is a staging optimization problem, it is convenient to model the mass of each stage separately. Here, the consumed propellant of each stage is defined as a state variable. These are denoted by Δmp1, Δmp2, and Δmp3 for the first, second, and third stages, respectively. Omitting the index that denotes the stage, the time derivative for these states is given by Eq. 14:

Δ m ˙ p = T v a c g 0 I s p (14)

where Tvac and Isp are the vacuum thrust and specific impulse of the stage engine, respectively, and g0 is the standard gravity acceleration at sea level. Naturally, when the stage engine is not active, there is no burned propellant being exhausted and the vehicle mass remains constant.

Here, to limit the scope of the analysis, only the first and second stages of the vehicle, along with the payload, are subjected to optimization, while the third stage is considered to have fixed and known characteristics. The optimizable stage size characteristic of the launcher is modeled as follows. The propellant masses of the first and second stages are given by Eq. 15:

m p 1 = k 1 m p , r e f 1 m p 2 = k 2 m p , r e f 2 (15)

where mp,ref1 and mp,ref2 are constant values defined by the user and k1 and k2 are optimizable parameters. If the propellant mass of the stage varies, its structural mass cannot remain the same: it is modeled here as a linear function of the stage propellant mass. The structural masses of the first and second stages are modeled as Eq. 16:

m s 1 = m s 0 1 + r 1 m p 1 m s 2 = m s 0 2 + r 2 m p 2 (16)

where ms01 and ms02 are the portion of the stage mass that is independent of the propellant mass (for example, the thrust vector control [TVC] mass) and r1mp1 and r2mp2 represent the portion of the structural mass that varies linearly with the propellant mass of the stage (for example, the motor case or the tank mass). The parameters ms0 and r are associated with the construction technological level of the stage. Regarding the payload, its mass is modeled as Eq. 17:

m c = k c m c , r e f (17)

where mc,ref is a constant value provided by the user and kc is an optimizable parameter.

With the previous definitions, it is possible to calculate the mass of the vehicle at any instant of the trajectory. At lift-off, the mass is the sum of propellant and structural masses of all stages and the mass of the payload (the fairing is considered to be comprised in the second stage). Denoting the structural and propellant masses of the third stage by ms3 and mp3, respectively, the initial mass of the launcher is given by Eq. 18:

m 0 = m s 0 1 + k 1 r 1 + 1 m p , r e f 1 + m s 0 2 + k 2 r 2 + 1 m p , r e f 2 + m s 3 + m p 3 + k c m c , r e f (18)

where k1, k2, and kc are optimizable parameters of the problem. Along the trajectory, the mass of the vehicle is the sum of the structural and remaining propellant masses of the stages not yet jettisoned, and the payload mass (the remaining propellant of a stage is equal to the total minus the consumed propellant). The discontinuities that appear in the vehicle mass along the trajectory are handled with the introduction of different trajectory phases.

Phases of the trajectory

It is usual that, in complex OCP, the trajectory is divided into phases, allowing the use of a different problem modeling in each phase. The launch vehicle trajectory considered here is composed of seven phases, with a phase change defined whenever one of the following events happens: motor ignition; vehicle structure jettisoning; or control law change. The trajectory phases are:

  • Vertical climb;

  • Pitch-over maneuver – first phase;

  • Pitch-over maneuver – second phase;

  • Burn of first stage motor;

  • Burn of second stage motor;

  • Coast phase; and

  • Burn of third stage motor.

The pitch-over maneuver is divided into two phases due to the necessity of different control laws in each one, as described in Table 1. The initial and final times of each phase are viewed as parameters in the OCP. For phase i, they are denoted by tI(i) and tF(i), respectively. The jettisoning of first stage occurs at the end of phase 4 and the jettisoning of second stage occurs at the end of phase 5, along with the fairing separation. It is interesting to note that the phase sequence previously presented does not change due to the optimization of the sizes of the first and second stages. The phase durations, on the other hand, are directly affected by the amount of propellant in the stage.

Table 1
Phases of the OCP.
Control variables and parameters

The control variables of the OCP are the pitch and yaw attitude angles (as already explained, roll is not relevant here). In each phase, a control law is imposed for each of these angles. The possible forms for the control laws are: constant value, described by a single parameter; linear variation, described by two parameters; or gravity-turn law, which has no associated parameter, since the attitude angle is equal to the flight path angle (this law results in zero angle of attack and is only applicable for pitch angle). The specific law adopted for each phase is shown in Table 1.

The OCP parameters for each phase comprise the phase times, the control law parameters, and the stages and payload sizing parameters. The vector p has, therefore, the following general form (Eq. 19):

p ( i ) = t I ( i ) t F ( i ) θ I ( i ) θ F ( i ) ψ I ( i ) ψ F ( i ) k 1 ( i ) k 2 ( i ) k c ( i ) T (19)

where θI and θF are the parameters used to define the pitch control law, while ΨI and ΨF define the yaw control law. Depending on the control law and on the current vehicle configuration, not all parameters previously shown will be present in the parameter vector of a specific phase.

Constraints definition

Vehicle lift-off takes place from a known location on Earth. In addition, no propellant has been burned before lift-off. To correctly model these characteristics in the OCP, constraints should be imposed on the vehicle initial state. The following initial constraints are considered here (Eq. 20):

R ( 1 ) t I R 0 = 0 l ( 1 ) t I l 0 = 0 λ ( 1 ) t I λ 0 = 0 V R ( 1 ) t I = 0 V l ( 1 ) t I Ω E R 0 cos   λ 0 = 0 V λ ( 1 ) t I = 0 Δ m p 1 ( 1 ) t I = 0 Δ m p 2 ( 1 ) t I = 0 Δ m p 3 ( 1 ) t l = 0 (20)

The target orbit considered in this work is circular and has a predefined altitude and predefined inclination. To assure that these characteristics are fulfilled in the OCP, constraints should be imposed on the vehicle final state. The final constraints considered here are (Eq. 21):

h e q ( 7 ) t F h f = 0 V ( 7 ) t F μ E / R E + h f = 0 V R ( 7 ) t F = 0 i ( 7 ) t F i f = 0 (21)

where heq(t) = R(t) – RE is the vehicle equatorial altitude, V(t) denotes its inertial speed, and i(t) is the inclination of the orbit. The constraints related to the vehicle altitude and velocity are clearly stated as functions of the vehicle states. As for the inclination of the orbit, it can be derived from the states; it is the angle between the vehicle angular momentum vector, given by the vector product between its position and velocity vectors, and the Earth axis of rotation.

Other constraints, necessary to assure desired characteristics of the vehicle configuration and the trajectory, are also imposed (Eq. 22):

θ ( 1 ) t I θ 0 = 0 t F ( 7 ) t I ( 7 ) Δ t 3  stage  = 0 m p 1 Δ m p 1 ( 4 ) t F > 0 m p 2 Δ m p 2 ( 5 ) t F > 0 (22)

The first relation in Eq. 22, with θ0 = 90°, assures that the vehicle takes off vertically. The second relation is imposed to guarantee that the burn time of the third stage motor is respected, since this stage is not part of the staging optimization problem and its burn time, Δt3stage, is known a priory. Finally, the last two inequality constraints are necessary to guarantee that the quantity of burned propellant at the burnout of the first and second stages do not exceed the propellant mass hold by the stage.

The last class of constraints are those necessary to assure that phase times, state variables, control variables, and parameters are continuous between phases. These are called phase connection constraints. Denoting a generic state, control, and parameter by x, u, and p, respectively, the phase connection constraints are (Eq. 23):

t I ( i + 1 ) t F ( i ) = 0 x ( i + 1 ) t I x ( i ) t F = 0 u ( i + 1 ) t I u ( i ) t F = 0 p ( i + 1 ) p ( i ) = 0 . (23)
Performance measure definition

The last element of the OCP to be defined is the performance measure. Since the objective is to maximize the payload mass, the performance measure is defined as Eq. 24:

J = k c ( 7 ) (24)

where the negative sign is necessary since the solver is implemented to minimize the objective function, as shown in the next section. The choice for the parameter kc in phase 7 is arbitrary since this parameter should be continuous between phases. Table 1 resumes all the characteristics of the OCP previously described, except for the phase connection conditions, which are present in all phases. The term “Optm” in Table 1 indicates that the parameter value is optimizable.

A direct approach to resolution of OCP

There are different approaches to solving OCP (Betts 2010; Rao 2010). In the indirect approach, the necessary conditions provided by Pontryagin’s maximum principle (Pontryagin et al. 1962) are applied and a two-point or a multipoint boundary value problem is obtained. Such boundary value problems must be solved by iterative techniques. In the direct approach based on transcription methods, in turn, a discrete approximation of the original problem is defined, and a minimizer is sought based on nonlinear programming. The discretization transforms the infinite-dimensional function space of the original OCP into a finite-dimension parameter space. Variational versions of well-known gradient methods also belong to the class of direct methods (Miele et al. 1970).

The approach used in this work to solve the OCP associated with the staging optimization of a launch vehicle belongs to the class of direct methods based on transcription methods. This approach is described in detail in da Silveira and da Silva Fernandes (2023). To transform the OCP into an NLPP, two different transcription methods have been considered: Hermite-Simpson collocation and multiple-shooting. The NLPP is then solved using a gradient-based algorithm with sequential gradient-restoration phases as proposed by Miele et al. (1969). Transcription methods and the SGRA are briefly described in the next subsections. More details can be found in the original paper, where the method is successfully applied to the problem of a launch vehicle trajectory optimization.

Transcription methods

The method used to transform an OCP into a NLPP is called a transcription method (Hull 1997). It replaces the continuous state and/or control variables with state and/or control parameters at a set of discrete time instants, called nodes. These parameters are then used to generate approximations for the constraints and the performance measure of the OCP in the form of multivariate functions in the NLPP. Numerical integration schemes, either implicit or explicit, are used to guarantee the correct dynamic behavior of the system. The discretization of states and controls has the potential of introducing discontinuities in these variables, and these are prevented by imposing equality constraints between consecutive segments of the discretization mesh.

Two transcription methods are considered here. In the Hermite-Simpson collocation method (Enright and Conway 1992; Hargraves and Paris 1987), states are approximated by cubic polynomials within each segment of the discretization. An implicit scheme, resulting from the application of Simpson’s quadrature associated with a cubic Hermite interpolant, is used for the integration of state variables. The second transcription algorithm is the multiple-shooting method (Bock and Plitt 1984), where an explicit scheme is used to integrate the states within each segment of the discretization, employing piecewise linear polynomials to approximate the control variables in the segment.

SGRA and its extension

To solve the NLPP that results from the transcription of the staging optimization problem, the authors have proposed an extension of the method known as the SGRA (da Silveira and da Silva Fernandes 2023).

The SGRA corresponds to an algorithm to find an optimizer for a multivariate function f(x), subject to the set of equality constraints c(x) = 0. The vector x, called the decision vector, has dimension nx, the function f, called the objective function, is scalar, the vector of constraints c has dimension nc, and nc < nx (Miele et al. 1969). The SGRA is composed of the alternate succession of two distinct phases: the first is a gradient phase and the second is a restoration phase. In the gradient phase of iteration i, starting from a point xi that satisfies the constraints, a step Δxi is calculated such that the value of the objective function at the new point yi = xi + Δxi is lower than at point xi. This is accomplished considering linear approximations for the objective and constraint functions and generally leads to a point yi that does not satisfy the constraints. In the restoration phase, starting from point yi, one or several steps are taken to bring back the constraint satisfaction, leading to a point xi+1, where c(xi+1) = 0 and f(xi+1) < f(xi).

An optimal solution is found when the current estimate xi, considered to satisfy the constraints to a certain accuracy, also satisfies an optimality condition. To state this condition, it is necessary to define the augmented function of the problem (Eq. 25):

F ( x , λ ) = f ( x ) + λ T c ( x ) (25)

where λ is a vector of Lagrange multipliers, calculated at each step of the algorithm. Denoting the vector of partial derivatives of F with respect to the components of x as xF, the optimality condition is stated as Eq. 26:

x T F x i , λ i x F x i , λ i ϵ (26)

where is a small value provided by the user. This relation is a measure of how close the point xi is from the necessary optimality condition for a constrained function optimization problem.

In order to apply the SGRA to more complex problems, the authors have proposed an extension to the algorithm to deal with inequality constraints (da Silveira and da Silva Fernandes 2023). Firstly, the step size calculation during the gradient and restoration phases are modified so that all components xj of a solution candidate xi satisfy the simple bounds xLjxjxUj, where xLj and xUj are the lower and upper boundaries of xj. Then, resorting on this capability of dealing with simple bounds, more general inequality constraints of the type d(x) ≤ 0 are addressed with the use of slack variables, where each inequality is transformed into an equality constraint d(x) + s = 0, with the definition of a new variable s that should satisfy the simple bound s ≥ 0. In the results presented next, the extended algorithm is called ESGRA.

Staging optimization of VLM-1 launch vehicle

This section presents the results of the application of the staging optimization methodology to the Brazilian microsatellite launcher VLM-1 (da Silveira and Carrara 2015). These results are obtained with the aid of a numerical tool implemented using the Matlab programming language. The tool performs the transcription of the OCP and applies the ESGRA to solve the resulting NLPP. Both transcription methods presented earlier are used: Hermite-Simpson collocation and multiple-shooting. To validate the results, the commercial software ASTOS 7 (ASTOS Model Library, version 7.0.2, ASTOS Solutions GmbH 2011), capable of dealing with a multitude of space related optimization problems, is used as reference. ASTOS results are obtained with the direct approach: multiple-shooting is used to generate the NLPP and a sequential quadratic programming algorithm is used to solve the problem (Boggs and Tolle 1995; Gill et al. 1981).

VLM-1 is a three-stage vehicle, each stage being equipped with a solid propellant motor. In its nominal design, the first and second stages are very similar to each other: they use the S50 motor and hold 12 tons of propellant each. Both stages are equipped with TVC. The third stage uses a fixed nozzle S44 motor, which holds 800 kg of propellant. The fact that the first and second stages of VLM-1 are very similar to each other (same motor and same propellant mass) is not optimal in terms of the vehicle payload capacity: to achieve the optimal performance for a fixed quantity of propellant, the size of the stages shall decrease when moving from the first to the last stage (Cornelisse et al. 1979). The objective of the staging optimization presented here is to find the optimal propellant mass distribution between the first and second stages of VLM-1, along with the optimal trajectory, which results in the maximum value for the payload mass that can be inserted into a target orbit.

The target orbit considered here possesses the following characteristics: circular low Earth orbit (LEO) with equatorial altitude and inclination of, respectively, 300 km and 3°. The launch site corresponds to Centro Espacial de Alcântara, northeast of Brazil. Referring to Eqs. 20 and 21, the initial and final trajectory conditions are Eq. 27

R 0 = 6378 k m , l 0 = 44 , 3 , λ 0 = 2 , 3 , V l 0 = Ω E R 0 cos   λ 0 = 464 , 7 m / s , h f = 300 k m , V f = 7725 , 8 m / s , i f = 3 (27)

Two different scenarios regarding the staging optimization of VLM-1 are addressed: the first test case considers the propellant mass optimization of the second stage alone, keeping the original VLM-1 first stage; the second test case regards the joint optimization of the propellant masses of the first and second stages, leading to the optimal propellant distribution between them. The general problem formulation presented earlier is adjusted accordingly for each scenario.

The nominal values considered in this work for mass and propulsive properties of VLM-1 are shown in Table 2, while Fig. 2 presents the aerodynamic coefficient CA, as a function of Mach number, for the flight of the first and second stages. Since the flight of the third stage occurs outside the atmosphere, the aerodynamic influence is neglected. The portion of the structural mass that is independent of the propellant mass is indicated inside parentheses in Table 2, for the second stage, this encompasses a control bay and the fairing. For the staging optimization problem, the mass parameters of Eq. 16 are considered to be the same as for VLM-1. According to the values in Table 2, ms01 = 0 kg, ms02 = 383 kg, r1 = 0.141, and r2 = 0.136. Thrust and specific impulse values of the first and second stages are also the same as for VLM-1. Thrust correction due to the atmospheric pressure is neglected. The vehicle diameter is considered to remain constant for different stage sizes, allowing the simplifying hypothesis that the aerodynamic coefficient does not vary with the stage size. It is, therefore, the same as for VLM-1. The aerodynamic reference area is Sref = 1.67 m2. As for the atmospheric properties patm and ρatm, the 1976 U.S. Standard Atmosphere model is used (United States 1976).

Table 2
Nominal masses and propulsive properties of VLM-1 stages.
Figure 2
Aerodynamic coefficient for the flight of the first and second stages.

To analyze the results of VLM-1 staging optimization, it is necessary to know the payload capacity of the vehicle in its original configuration. This result comes from a trajectory optimization problem and was obtained by the authors in the first test case presented in a previous work (da Silveira and da Silva Fernandes 2023): considering the original 12,040 kg of propellant mass in each of the first and second stages and for the same target orbit and trajectory characteristics as modeled here, VLM-1’s optimum payload capacity is 129.0 kg.

Scenario 1: optimization of the second stage

In the first test case, the propellant mass of the second stage of VLM-1 is considered an optimizable parameter, while the first stage is fixed and maintains its original quantity of propellant. The structural mass of the second stage is a linear function of its propellant mass. The goal of this test case is to confirm that the vehicle’s current configuration, with equal first and second stages, is not optimal in terms of performance.

Table 3 shows the main characteristics of the transcriptions used. Due to the state variables integration scheme of collocation methods, a high number of nodes in the discretization mesh is usually necessary to ensure an acceptable approximation for the states, resulting in a relatively large NLPP.

Table 3
Transcription characteristics of the first scenario.

The payload mass capacity of the vehicle, after staging optimization, is shown in Table 4, along with the optimal propellant mass of the second stage. Reference values are those obtained with ASTOS software. Results from ESGRA regarding the payload capacity are very close to those from ASTOS, for both transcription methods. When compared to the original VLM-1 configuration, the vehicle’s performance increases from 129 kg to 147 kg (a 14% increase) for a decrease of approximately 3,700 kg in the propellant mass of the second stage (a 30% decrease). Figure 3 illustrates the different vehicle configurations and performances. The remaining propellant mass at second stage burnout is close to 0, meaning that the inequality constraint associated with this parameter is active at the optimal solution.

Table 4
Optimal payload and propellant masses for the first scenario.
Figure 3
Propellant distribution and payload mass for different vehicle configurations.

To confirm the validity of the results, it is necessary to verify if the launch vehicle reaches the desired orbit. Table 5 shows some final orbit parameters, comparing them with the requirements of the target orbit, where it is clear that the desired orbit is achieved.

Table 5
Final orbit parameters for the first scenario.

The optimality measure of each solution candidate along the iterations is shown in Fig. 4, while Fig. 5 shows the evolution of the payload mass. The initial guess is indicated by the letter G and the first restoration phase is indicated by the letter R. The number of iterations is similar for both transcription methods: 17 for collocation and 15 for multiple-shooting. A relatively smooth path towards the optimum is observed.

Figure 4
Optimality measure along iterations for the first scenario
Figure 5
Payload mass along iterations for the first scenario.

Time histories for some of the main trajectory parameters are shown in Figs. 6-9. Phase changes are indicated by vertical dashed lines. Initial guess, ASTOS solution and ESGRA solutions using collocation and multiple-shooting are presented. The initial guess corresponds to a trajectory simulation using values provided by the user for all control laws and optimizable parameters. Pitch angle, shown in Fig. 6, varies from 90° at take-off to 0° at orbit insertion. This control is mainly related to the final vehicle altitude and velocity, which determines the orbit eccentricity. After the vertical climb and pitch-over phases, the pitch angle corresponds to the necessary angle to keep the angle of attack equal to 0 during the first and second stage burns. Then, the vehicle is reoriented during the coast phase, and the trajectory ends with the third stage burn. The yaw angle, shown in Fig. 7, is mainly related to the final orbit inclination. Since the target inclination is very close to the launch point latitude, yaw presents a very small variation during flight. In general, the control variables provided by both transcription methods are very similar when compared to the results from ASTOS. The difference observed in the yaw angle obtained with collocation, which appears to be relevant in Fig. 7, is actually of the order of 1°. Altitude (Fig. 8) and velocity (Fig. 9) profiles are also similar to ASTOS solutions. The small difference observed in the velocity during the third stage burn is due to differences in the coast phase duration of each optimal solution.

Figure 6
Pitch angle variation for the first scenario.
Figure 7
Yaw angle variation for the first scenario.
Figure 8
Altitude variation for the first scenario
Figure 9
Inertial velocity variation for the first scenario.
Scenario 2: optimization of the first and second stages

In this test case, both first and second stages are considered to have optimizable propellant masses, with their structural masses modeled as linear functions of the respective propellant mass. With the objective of determining the optimal propellant distribution between the first and second stages of VLM-1, a new constraint on the sum of the propellant of these stages is imposed: the total propellant quantity in the first and second stages of the new configuration should be the same as in the original stages of VLM-1. Therefore, in addition to the constraints described in the problem formulation, the following constraint is imposed (Eq. 28):

m p 1 + m p 2 = m p total  (28)

where mptotal = 24,080 kg corresponds to the sum of the propellant mass in the first and second stages of VLM-1 in its nominal configuration.

Table 6 shows the main characteristics of the transcriptions used. The same discretization meshes from the previous case are used here. The increase in the number of optimizable parameters and equality constraints is due to the new parameter associated with the first stage propellant optimization, and the related phase connection relations, besides the optimizable burn duration of the first stage.

Table 6
Transcription characteristics of the second scenario.

The payload mass capacity of the vehicle, after staging optimization, is shown in Table 7, along with the optimal propellant mass of first and second stages. Results from ESGRA regarding payload mass are very similar to those from ASTOS, used as reference. It is clear that a substantial increase in the payload mass could be achieved with a better propellant distribution between first and second stages of VLM-1 (performance increase of 56%). On the other hand, the optimal propellant quantity in each stage varies among ASTOS and ESGRA results, which could indicate a relatively low sensitivity of the optimal solution with respect to small variations in the propellant masses, or the convergence to different local minima. The largest difference is observed when the propellant mass of second stage obtained with collocation is compared with that from ASTOS (13.4% difference). Figure 3 illustrates the different vehicle configurations and performances. Again, the remaining propellant mass at first and second stage burnouts is close to 0, meaning that the inequality constraints associated with these parameters are active at the optimal solution.

Table 7
Optimal payload and propellant masses for the second scenario.

Regarding the achieved orbit, Table 8 shows some final parameters of the trajectory, comparing them with the requirements of the target orbit. It is clear that the desired orbit is achieved.

Table 8
Final orbit parameters for the second scenario.

Optimality measure and payload mass along the iterations are shown in Figs. 10 and 11, respectively. Considerably more iterations are necessary for the collocation transcription, which also presents a more irregular path towards the optimum.

Figure 10
Optimality measure along iterations for the second scenario.
Figure 11
Payload mass along iterations for the second scenario.

Time histories of the control variables are shown in Figs. 12 and 13. The same observations made for the first scenario are applicable here. The yaw angle obtained with collocation is, again, the result that differs the most from that of ASTOS, with a difference of the order of 4°. Finally, time histories of altitude and velocity are shown in Figs. 14 and 15. It is interesting to note that the velocity gain provided by first and second stages are more evenly distributed now when compared to the first scenario. Again, a good agreement is observed between ESGRA and ASTOS results.

Figure 12
Pitch angle variation for the second scenario.
Figure 13
Yaw angle variation for the second scenario.
Figure 14
Altitude variation for the second scenario.
Figure 15
Inertial velocity variation for the second scenario.

CONCLUSION

This work presents results in the area of staging optimization of launch vehicles. The problem is formulated as an OCP, where staging and trajectory optimization are handled in a unified manner. To search for an optimal solution, the problem is first discretized and transformed into a NLPP. The resolution of the NLPP is then obtained using an extension, proposed by the authors, of the SGRA. This methodology is applied to the staging optimization of the VLM-1 Brazilian microsatellite launch vehicle, in two different scenarios.

The results for the first scenario show that the vehicle’s payload capacity could be improved by around 14% if its second stage were replaced by a smaller one, with approximately 30% less propellant. This has the potential to reduce the payload cost per kilogram to orbit, both due to the increase in payload capacity and to the decrease in the vehicle’s overall cost owing to a lower propellant consumption. In the second scenario, the optimal propellant mass distribution between the first and second stages of VLM-1 is determined. It is clear that, with the current total propellant mass held by these stages, the vehicle’s performance could be considerably improved if this mass were distributed in a better way, resulting in a 56% increase in payload capacity. The results of both scenarios underscore the importance of staging optimization analyses during the conceptual phase of a launch vehicle. For VLM-1, apparently, a choice was made in favor of a smaller development cost (only one motor development for both first and second stages) over a higher payload capacity.

The proposed methodology proved to be adequate for addressing the problem of staging optimization of the VLM-1 vehicle. However, it could be extended to a more general staging optimization problem.

ACKNOWLEDGMENTS

Not applicable

  • Peer Review History: Single Blind Peer Review.
  • FUNDING
    Conselho Nacional de Desenvolvimento Científico e Tecnológico
    Grant No: 305981-2020-0

DATA AVAILABILITY STATEMENT

Data sharing is not applicable

REFERENCES

  • Adkins CN (1970) Optimization of multistage rockets including drag. J Spacecr Rockets 7(6):751-755. https://doi.org/10.2514/3.30032
    » https://doi.org/10.2514/3.30032
  • Betts JT (1998) Survey of numerical methods for trajectory optimization. J Guid Control Dyn 21(2):193-207. https://doi.org/10.2514/2.4231
    » https://doi.org/10.2514/2.4231
  • Betts JT (2010) Practical methods for optimal control and estimation using nonlinear programming. 2nd ed. Philadelphia: Siam.
  • Bock HG, Plitt KJ (1984) A multiple shooting algorithm for direct solution of optimal control problems. IFAC Proc 17(2):1603-1608. https://doi.org/10.1016/S1474-6670(17)61205-9
    » https://doi.org/10.1016/S1474-6670(17)61205-9
  • Boggs PT, Tolle JW (1995) Sequential quadratic programming. Acta Numer 4(1):1-51. https://doi.org/10.1017/S0962492900002518
    » https://doi.org/10.1017/S0962492900002518
  • Bryson AE, Ho YC (1975) Applied optimal control. Washington, DC: Hemisphere.
  • Cho B, Jo BU, Ahn J (2021) Integrated framework for staging and trajectory optimization of a launch vehicle considering range safety operations. Int J Aeronaut Space Sci 22:963-973. https://doi.org/10.1007/s42405-020-00348-6
    » https://doi.org/10.1007/s42405-020-00348-6
  • Cornelisse JW, Schöyer HFR, Wakker KF (1979) Rocket propulsion and spaceflight dynamics. Aerospace Engineering Series. London: Pitman.
  • da Silveira G, Carrara V (2015) A six degrees-of-freedom flight dynamics simulation tool of launch vehicles. J Aerosp Technol Manag 7(2):231-239. https://doi.org/10.5028/jatm.v7i2.433
    » https://doi.org/10.5028/jatm.v7i2.433
  • da Silveira G, da Silva Fernandes S (2023) Trajectory optimization of the Brazilian multistage launch vehicle VLM-1. J Braz Soc Mech Sci Eng 45:564. https://doi.org/10.1007/s40430-023-04490-6
    » https://doi.org/10.1007/s40430-023-04490-6
  • Enright PJ, Conway BA (1992) Discrete approximations to optimal trajectories using direct transcription and nonlinear programming. J Guid Control Dyn 15(4):994-1002. https://doi.org/10.2514/3.20934
    » https://doi.org/10.2514/3.20934
  • Gill PE, Murray W, Wright MH (1981) Practical optimization. London: Academic.
  • Hargraves CR, Paris SW (1987) Direct trajectory optimization using nonlinear programming and collocation. J Guid Control Dyn 10(4):338-342. https://doi.org/10.2514/3.20223
    » https://doi.org/10.2514/3.20223
  • Hughes PC (2004) Spacecraft attitude dynamics and control. Mineola: Dover.
  • Hull DG (1997) Conversion of optimal control problems into parameter optimization problems. J Guid Control Dyn 20(1):57-60. https://doi.org/10.2514/2.4033
    » https://doi.org/10.2514/2.4033
  • Jamilnia R, Naghash A (2012) Simultaneous optimization of staging and trajectory of launch vehicles using two different approaches. Aerosp Sci Technol 23(1):85-92. https://doi.org/10.1016/j.ast.2011.06.013
    » https://doi.org/10.1016/j.ast.2011.06.013
  • Koch AD (2019) Optimal staging of serially staged rockets with velocity losses and fairing separation. Aerosp Sci Technol 88(1):65-72. https://doi.org/10.1016/j.ast.2019.03.019
    » https://doi.org/10.1016/j.ast.2019.03.019
  • Malina FJ, Summerfield M (1947) The problem of escape from the earth by rocket. J Aeronaut Sci. 14(8):471-480. https://doi.org/10.2514/8.1417
    » https://doi.org/10.2514/8.1417
  • Miele A, Huang HY, Heideman JC (1969) Sequential gradient-restoration algorithm for the minimization of constrained functions – ordinary and conjugate gradient versions. J Optim Theory Appl 4:213-243. https://doi.org/10.1007/BF00927947
    » https://doi.org/10.1007/BF00927947
  • Miele A, Pritchard RE, Damoulakis JN (1970) Sequential gradient-restoration algorithm for optimal control problems.J Optim Theory Appl 5:235-282. https://doi.org/10.1007/BF00927913
    » https://doi.org/10.1007/BF00927913
  • Perondi LF (2023) The coming of age of space satellite industry: transitioning from a growth to a maturity life cycle phase.J Aerosp Technol Manag 15:e0323. https://doi.org/10.1590/jatm.v15.1291
    » https://doi.org/10.1590/jatm.v15.1291
  • Pessoa Filho JB (2021) Space age: past, present and possible futures. J Aerosp Technol Manag 13:e3421. https://doi.org/10.1590/jatm.v13.1226
    » https://doi.org/10.1590/jatm.v13.1226
  • Pontryagin LS, Boltyanskii VG, Gamkrelidze RV, Mishechenko EF (1962) Mathematical theory of optimal processes. New York: Interscience.
  • Rao AV (2010) A survey of numerical methods for optimal control. Adv Astronaut Sci 135:497-528.
  • Srivastava T (1966) Optimum staging with varying thrust attitude angle. Def Sci J 16(3):153-164. Available at: https://publications.drdo.gov.in/ojs/index.php/dsj/article/view/7243
    » https://publications.drdo.gov.in/ojs/index.php/dsj/article/view/7243
  • Tawakley V (1967) On the calculation of optimum mass distribution of a multi-stage rocket vehicle. Def Sci J 17(2):53-66. Available at: https://publications.drdo.gov.in/ojs/index.php/dsj/article/view/7373
    » https://publications.drdo.gov.in/ojs/index.php/dsj/article/view/7373
  • United States. Air Force; NOAA; NASA (1976) U.S. standard atmosphere. Washington, D.C.: NASA.
  • Vallado DA, Macclain WD (2013) Fundamentals of astrodynamics and applications. 4th ed. Hawthorne: Microcosm.
  • Villas Bôas DJF, Souza CHM, Silva FM, Moraes A (2020) Proposal of low cost launchers for scientific missions using cubesats. Adv Sp Res 66(1):162-75. https://doi.org/10.1016/j.asr.2019.12.012
    » https://doi.org/10.1016/j.asr.2019.12.012
  • Villas Bôas DJF, Pessoa Filho JB, Moraes AO, Souza CHM (2023) Innovative and low-cost launch systems. Next Generation CubeSats and SmallSats 403-419. https://doi.org/10.1016/B978-0-12-824541-5.00005-4
    » https://doi.org/10.1016/B978-0-12-824541-5.00005-4
  • Weigel N, Well KH (2000) Dual payload ascent trajectory optimization with a splash-down constraint. J Guid Control Dyn 23(1):45-52. https://doi.org/10.2514/2.4485
    » https://doi.org/10.2514/2.4485

Edited by

Publication Dates

  • Publication in this collection
    15 July 2024
  • Date of issue
    2024

History

  • Received
    17 Apr 2024
  • Accepted
    08 May 2024
location_on
Departamento de Ciência e Tecnologia Aeroespacial Instituto de Aeronáutica e Espaço. Praça Marechal do Ar Eduardo Gomes, 50. Vila das Acácias, CEP: 12 228-901, tel (55) 12 99162 5609 - São José dos Campos - SP - Brazil
E-mail: submission.jatm@gmail.com
rss_feed Acompanhe os números deste periódico no seu leitor de RSS
Acessibilidade / Reportar erro