Trajectory Optimization of Launch Vehicles Using Object-oriented Programming

1.Universidade Federal do ABC – Centro de Engenharia, Modelagem e Ciências Sociais Aplicadas – Departamento de Engenharia Aeroespacial – Santo André/ SP – Brazil. 2.Instituto Nacional de Pesquisas Espaciais – Departamento de Engenharia e Tecnologia Espacial – São José dos Campos/SP – Brazil. 3.Ostbayerische Technische Hochschule Regensburg – Department of Mechanical Engineering – Regensburg – Germany. Correspondence author: Fábio Antônio da Silva Mota | Universidade Federal do ABC – Centro de Engenharia, Modelagem e Ciências Sociais Aplicadas – Departamento de Engenharia Aeroespacial | Av. dos Estados, 5001 | CEP: 09210-580 – Santo André/SP – Brazil | E-mail: mota.fabio@ufabc.edu.br Received: Jun. 8, 2017 | Accepted: Nov. 6, 2017 Section Editor: Luiz Martins-Filho ABSTRACT: The aim of this study is to model launch vehicles with focus on 3-DOF trajectory optimization using a modular approach. Despite the large number of operational launch vehicles, they usually consist of basic components and subsystems. In other words, a launch vehicle is an assembly of stages, which in turn is divided into propellant system and engine, and the engine is an assembly of basic components such as pumps, turbines, combustion chamber, and nozzle. To allow future extension and reuse of the codes, a modular structure using object-oriented programming is used. Two formulations of state equations of the trajectory and two optimization methods are described. The launch vehicle performance will be measured by payload mass for a given mission. The simulations of the VLS-1, Ariane 5 and VLS-Alfa were performed and showed good agreement with the literature.

xx/xx 02/15 the multiple-shooting method, and collocation methods (Brown et al. 1969;Teren and Spurlock 1966;Miele and Wang 2003).Presumably because of the possibility of solving complex problems with a minimum effort of mathematical analysis, the direct method is the one chosen by most researchers (Hargraves and Paris 1987;Seywald 1994;Herman and Conway 1996;Balesdent 2011).One of the most popular software used extensively in many publications is called POST (Program to Optimize Simulated Trajectories) (Brauer et al. 1977).POST is a commercial program that has been used successfully to solve a wide variety of atmospheric ascent and reentry problems, as well as exoatmospheric orbital transfer problems.In the framework of this method, the problem is characterized by a set of parameters that define the control law.This problem is a typical Non Linear Programming Problem (NLP) and can be solved using classical Gradient-based methods (deterministic methods) such as Sequential Quadratic Program (SQP) or by heuristic methods.According to Betts (1998), heuristic optimization algorithms are not computationally competitive with gradient methods.Even though presumably due to ease of implementation without a detailed understanding of the system, in the last two decades lots of papers using Particle Swarm Optimization (PSO), genetic algorithms (GA), among others were applied to solve trajectory optimization problems.As for indirect methods, the direct methods can be categorized in direct (multiple) shooting or collocation.In the case where only the control variables are adjusted by a function, the method is called a shooting method.When both the state and control are parameterized, the method is called a collocation method.A well-known software developed by the University of Stuttgart, which addresses the direct collocation method, is the AeroSpace Trajectory Optimization Software (ASTOS).
For either direct or indirect approaches, perhaps the most important benefit gained from a multiple shooting formulation compared to its precursor (single shooting) is enhanced robustness.To take advantage of both methods previously described, a hybrid method can also be considered (von Stryk and Bulirsch 1992;Pontani and Teofilatto 2014;Gath and Calise 2001).The idea behind this approach is to divide the flight trajectory into two distinct phases, namely atmospheric and exoatmospheric phases, applying the direct method in the first phase and indirect method in the second one.Here, exoatmospheric phase means that the vehicle is virtually in vacuum space, i.e., the aerodynamic effects can be ignored.Pontani and Teofilatto (2014) proposed a simple method to evaluate the performance of multistage launch vehicles for given structural data, aerodynamic and propulsive parameters.
The purpose of this work is to develop a tool, which can be easily reused and extended, to model and simulate a launch vehicle with focus on trajectory.Thus, the tool is intended to be capable of modeling and optimizing flight trajectory until orbit injection.Two mathematical models for determination of performance of a launch vehicle will be described and discussed.The launch vehicle performance will be measured by payload mass for a given mission.

TRAJECTORY MODELING
Since the dawn of the space era launch vehicles are responsible to put satellites into orbit.This makes the cost of a satellite strongly related to the performance of the launch vehicle, which in turn depends on the trajectory profile.

ATMOSPHERE MODEL
The atmosphere can be seen as a layer of gases attached to the surface of the Earth by gravitational attraction.This work makes use of the standard atmosphere, which is modeled as adjacent layers of gases in which temperature depends on the altitude.In different layers, the temperature can be modeled as linear function of the altitude.

AERODYNAMICS
During the flight, a launch vehicle needs to cross the atmosphere in which reacts to the vehicle motion by means of aerodynamics forces.The drag force arises due to friction between the body and the fluid (Eq.1): 5 simple method to evaluate the performance of multistage launch vehicles for given structural data, aerodynamic and propulsive parameters.
The purpose of this work is to develop a tool, which can be easily reused and extended, to model and simulate a launch vehicle with focus on trajectory.Thus, the tool is intended to be capable of modeling and optimizing flight trajectory until orbit injection.Two mathematical models for determination of performance of a launch vehicle will be described and discussed.The launch vehicle performance will be measured by payload mass for a given mission.

TRAJECTORY MODELING
Since the dawn of the space era launch vehicles are responsible to put satellites into orbit.This makes the cost of a satellite strongly related to the performance of the launch vehicle, which in turn depends on the trajectory profile.

Atmosphere Model
The atmosphere can be seen as a layer of gases attached to the surface of the Earth by gravitational attraction.This work makes use of the standard atmosphere, which is modeled as adjacent layers of gases in which temperature depends on the altitude.In different layers, the temperature can be modeled as linear function of the altitude.

Aerodynamics
During the flight, a launch vehicle needs to cross the atmosphere in which reacts to the vehicle motion by means of aerodynamics forces.The drag force arises due to friction between the body and the fluid (Eq.1): (1) 2 ) ( 2where: D = drag force (N); ρ = density of the working fluid (kg/m 3 ); S = reference area (m 2 ); C D = drag coefficient (-); and V = air speed (m/s).
The lift is a reaction force to the angle of attack (Eq.2): combustion of the propellants and for liquid rocket engines there is still the sloshing, which is the movement of fluid within the tanks and pipes and rotating equipment such as turbines and pumps.Especially for large launch vehicles, the deflection of the structure should be considered as well (Cornelisse et al. 1979).However, this research is focused on the reference trajectory, thus treatment of the translational motion is sufficient to fulfill this task.Two mathematical models for the trajectory are presented in the following sections.The first one was taken from Schlingloff (2005) and the second one from Tewari (2007).

First Formulation
In this formulation, the spherical celestial (inertial) coordinates and a moving coordinates in the orbit plane were considered.Both reference frames have the origin on the Earth center (Fig. 1).The vector of state variables is conveniently chosen as , where: Ω = right ascension of the ascending node (rad); ω = argument of periapsis (rad); and T = temperature (K).Thus, the system of equation can be given as (Eqs.3-9): (3) combustion of the propellants and for liquid rocket engines there is still the sloshing, which is the movement of fluid within the tanks and pipes and rotating equipment such as turbines and pumps.Especially for large launch vehicles, the deflection of the structure should be considered as well (Cornelisse et al. 1979).However, this research is focused on the reference trajectory, thus treatment of the translational motion is sufficient to fulfill this task.Two mathematical models for the trajectory are presented in the following sections.The first one was taken from Schlingloff (2005) and the second one from Tewari (2007).

First Formulation
In this formulation, the spherical celestial (inertial) coordinates and a moving coordinates in the orbit plane were considered.Both reference frames have the origin on the Earth center (Fig. 1).The vector of state variables is conveniently chosen as y(t) = [r(t) u(t) v(t) Ω(t) ι(t) ω(t)] T , where: Ω = right ascension of the ascending node (rad); ω = argument of periapsis (rad); and T = temperature (K).Thus, the system of equation can be given as (Eqs.3-9): combustion of the propellants and for liquid rocket engines there is still the sloshing, which is the movement of fluid within the tanks and pipes and rotating equipment such as turbines and pumps.Especially for large launch vehicles, the deflection of the structure should be considered as well (Cornelisse et al. 1979).However, this research is focused on the reference trajectory, thus treatment of the translational motion is sufficient to fulfill this task.Two mathematical models for the trajectory are presented in the following sections.The first one was taken from Schlingloff (2005) and the second one from Tewari (2007).

First Formulation
In this formulation, the spherical celestial (inertial) coordinates and a moving coordinates in the orbit plane were considered.Both reference frames have the origin on the Earth center (Fig. 1).The vector of state variables is conveniently chosen as T , where: Ω = right ascension of the ascending node (rad); ω = argument of periapsis (rad); and T = temperature (K).Thus, the system of equation can be given as (Eqs.3-9): combustion of the propellants and for liquid rocket engines there is still the sloshing, which is the movement of fluid within the tanks and pipes and rotating equipment such as turbines and pumps.Especially for large launch vehicles, the deflection of the structure should be considered as well (Cornelisse et al. 1979).However, this research is focused on the reference trajectory, thus treatment of the translational motion is sufficient to fulfill this task.Two mathematical models for the trajectory are presented in the following sections.The first one was taken from Schlingloff (2005) and the second one from Tewari (2007).

First Formulation
In this formulation, the spherical celestial (inertial) coordinates and a moving coordinates in the orbit plane were considered.Both reference frames have the origin on the Earth center (Fig. 1).The vector of state variables is conveniently chosen as T , where: Ω = right ascension of the ascending node (rad); ω = argument of periapsis (rad); and T = temperature (K).Thus, the system of equation can be given as (Eqs.3-9):
There are many methods to estimate the aerodynamic coefficients in the literature, i.e., the well-known tool Missile DATCOM, interpolation of available data from a given vehicle, closed formulas that consider contributions of shock wave in the rocket nose, body friction and base pressure, or even constant value for certain phases of the flight.Except for interpolation from real flight data, the methods to estimate the drag coefficient are quite inaccurate.Fortunately, because the essential acceleration phase begins in the exoatmospheric phase, usually the drag has little influence on launcher performance (Schlingloff 2005).

Equations of the Translational Motion
The modeling of the trajectory of a launch vehicle is usually performed by means of two reference frames (one with origin on the Earth center and the other one moving with the vehicle) and considerations or idealizations according to the requirements of the mission.To model the translational motion, the vehicle can be treated as a particle, ignoring the size and mass distribution.In modeling the rotational motion, the vehicle can be considered a rigid body, reducing the degrees of freedom from infinity (flexible body case) to just six (Tewari 2007).Strictly speaking, a launch vehicle is far from being considered a rigid body.Mass is continuously expelled due to where: L = lift force (N); S = reference area (m 2 ); C L = lift coefficient (-); and V = air speed (m/s).
There are many methods to estimate the aerodynamic coefficients in the literature, i.e., the well-known tool Missile DATCOM, interpolation of available data from a given vehicle, closed formulas that consider contributions of shock wave in the rocket nose, body friction and base pressure, or even constant value for certain phases of the flight.Except for interpolation from real flight data, the methods to estimate the drag coefficient are quite inaccurate.Fortunately, because the essential acceleration phase begins in the exoatmospheric phase, usually the drag has little influence on launcher performance (Schlingloff 2005).

EQUATIONS OF THE TRANSLATIONAL MOTION
The modeling of the trajectory of a launch vehicle is usually performed by means of two reference frames (one with origin on the Earth center and the other one moving with the vehicle) and considerations or idealizations according to the requirements of the mission.To model the translational motion, the vehicle can be treated as a particle, ignoring the size and mass distribution.In modeling the rotational motion, the vehicle can be considered a rigid body, reducing the degrees of freedom from infinity (flexible body case) to just six (Tewari 2007).Strictly speaking, a launch vehicle is far from being considered a rigid body.Mass is continuously expelled due to combustion of the propellants and for liquid rocket engines there is still the sloshing, which is the movement of fluid within the tanks and pipes and rotating equipment such as turbines and pumps.Especially for large launch vehicles, the deflection of the structure should be considered as well (Cornelisse et al. 1979).However, this research is focused on the reference trajectory, thus treatment of the translational motion is sufficient to fulfill this task.Two mathematical models for the trajectory are presented in the following sections.The first one was taken from Schlingloff (2005) and the second one from Tewari (2007).

First Formulation
In this formulation, the spherical celestial (inertial) coordinates and a moving coordinates in the orbit plane were considered.Both reference frames have the origin on the Earth center (Fig. 1).The vector of state variables is conveniently chosen as y(t) = [r(t) u(t) v(t) Ω(t) ι(t) ω(t)] T , where u = vertical velocity (m/s); v = horizontal velocity (m/s); where: Ω = right ascension of the ascending node (rad) and ω = argument of periapsis (rad).Thus, the system of equation can be given as (Eqs.Equations 3 to 5 of the system of differential equations are the dynamic equations of motion and Eqs.6 to 9 are the kinematic equations.The dynamic equations are derived by application of the Newton's Second Law resolved into components of the moving system.The kinematic equations are deducted into two steps: representation of the rotation velocity of the vehicle in a vector form, and applying the Euler angles (Schlingloff 2005).Equations 3 to 5 of the system of differential equations are the dynamic equations of motion and Eqs.6 to 9 are the kinematic equations.The dynamic equations are derived by application of the Newton's Second Law resolved into components of the moving system.The kinematic equations are deducted into two steps: representation of the rotation velocity of the vehicle in a vector form, and applying the Euler angles (Schlingloff 2005).

Second Formulation
The reference frames adopted in this modeling are the planet-fixed reference (SXYZ) frame and the local horizontal frame (oxyz), both are non-inertial (Fig. 2).Equations 3 to 5 of the system of differential equations are the dynamic equations of motion and Eqs.6 to 9 are the kinematic equations.The dynamic equations are derived by application of the Newton's Second Law resolved into components of the moving system.The kinematic equations are deducted into two steps: representation of the rotation velocity of the vehicle in a vector form, and applying the Euler angles (Schlingloff 2005).

Second Formulation
The reference frames adopted in this modeling are the planet-fixed reference (SXYZ) frame and the local horizontal frame (oxyz), both are non-inertial (Fig. 2).

Second Formulation
The reference frames adopted in this modeling are the planet-fixed reference (SXYZ) frame and the local horizontal frame (oxyz), both are non-inertial (Fig. 2).To derive the dynamic equations Newton's Second Law must be introduced (Eq.

18): (18)
Choosing the wind axes to express the forces on the body and doing the appropriate transformation to perform aI in the wind axes, the remaining equations to model the translational motion are obtained (Eqs.19-21): ( 19) Comparing Eqs. 13 and 14 we finally obtain the kinematic equations of motion (Eqs.15-17): To derive the dynamic equations Newton's Second Law must be introduced (Eq.

18): (18)
Choosing the wind axes to express the forces on the body and doing the appropriate transformation to perform aI in the wind axes, the remaining equations to model the translational motion are obtained (Eqs.19-21): ( 19) Comparing Eqs. 13 and 14 we finally obtain the kinematic equations of motion (Eqs.15-17): To derive the dynamic equations Newton's Second Law must be introduced (Eq.

18): (18)
Choosing the wind axes to express the forces on the body and doing the appropriate transformation to perform aI in the wind axes, the remaining equations to model the translational motion are obtained (Eqs.19-21): Comparing Eqs. 13 and 14 we finally obtain the kinematic equations of motion (Eqs.15-17): To derive the dynamic equations Newton's Second Law must be introduced (Eq.

18): (18)
Choosing the wind axes to express the forces on the body and doing the appropriate transformation to perform aI in the wind axes, the remaining equations to model the translational motion are obtained (Eqs.19-21): where: ωE = Earth rotation (rad/s).
With a convenient rotation matrix, Eq. 11 can be written only in terms of axes of the body as (Eq.12): The relative velocity can also be expressed as (Eqs.13 and 14): Comparing Eqs. 13 and 14 we finally obtain the kinematic equations of motion (Eqs.15-17): To derive the dynamic equations Newton's Second Law must be introduced (Eq.18): Choosing the wind axes to express the forces on the body and doing the appropriate transformation to perform a I in the wind axes, the remaining equations to model the translational motion are obtained (Eqs.19-21): where: ω E = Earth rotation (rad/s).
Equations 15-17 are the kinematical equations of motion and Eqs.19-21 are the dynamic equations.With the integration of the system of differential equations, the vector position and the vector velocity of the vehicle can be determined by the following equations (Eqs.22 and 23): In this section the techniques used to solve the trajectory optimization problem are presented.The first approach was based on Silva (1995).The second one describes a hybrid algorithm that merges the direct and indirect methods.

First Approach -Direct Method
The method applied within the framework of this approach is based on Silva (1995).A polynomial control function is used to model the flight profile.Four parameters, namely the coast time duration tcoast (if it is applied) and three parameters of the polynomial control function, are optimized in order to get the maximum payload mass.A code from Jacob (1972) written in FORTRAN is transcript to C++ language and adapted to solve the problem (Eq.24).The trajectory optimization problem can be formulated as (Eq.25): In this section the techniques used to solve the trajectory optimization problem are presented.The first approach was based on Silva (1995).The second one describes a hybrid algorithm that merges the direct and indirect methods.

First Approach -Direct Method
The method applied within the framework of this approach is based on Silva (1995).A polynomial control function is used to model the flight profile.Four parameters, namely the coast time duration tcoast (if it is applied) and three parameters of the polynomial control function, are optimized in order to get the maximum payload mass.A code from Jacob (1972) written in FORTRAN is transcript to C++ language and adapted to solve the problem (Eq.24).The trajectory optimization problem can be formulated as (Eq.25): (12) In this section the techniques used to solve the trajectory optimization problem are presented.The first approach was based on Silva (1995).The second one describes a hybrid algorithm that merges the direct and indirect methods.

First Approach -Direct Method
The method applied within the framework of this approach is based on Silva (1995).A polynomial control function is used to model the flight profile.Four parameters, namely the coast time duration tcoast (if it is applied) and three parameters of the polynomial control function, are optimized in order to get the maximum payload mass.A code from Jacob (1972) written in FORTRAN is transcript to C++ language and adapted to solve the problem (Eq.24).The trajectory optimization problem can be formulated as (Eq.25): (12) In this section the techniques used to solve the trajectory optimization problem are presented.The first approach was based on Silva (1995).The second one describes a hybrid algorithm that merges the direct and indirect methods.

First Approach -Direct Method
The method applied within the framework of this approach is based on Silva (1995).A polynomial control function is used to model the flight profile.Four parameters, namely the coast time duration tcoast (if it is applied) and three parameters of the polynomial control function, are optimized in order to get the maximum payload mass.A code from Jacob (1972) written in FORTRAN is transcript to C++ language and adapted to solve the problem (Eq.24).The trajectory optimization problem can be formulated as (Eq.25): (12) It's known that if one has an inertial vector position and a velocity vector of a given body in orbit, the orbital elements (or Keplerian elements) can be readily determined.Thus, to get the orbital elements, it is necessary to perform an appropriate matrix rotation to obtain the desired inertial vectors.

OPTIMIZATION
In order to obtain the maximal payload capacity of a given launch vehicle, and consequently, make the access to space less expensive, trajectory optimization techniques have been for decades a subject of intense research.The trajectory optimization can be categorized basically into direct and indirect methods.To take advantage of both methods, a combination of both techniques can also be done, i.e., a hybrid method can also be considered.

METHODOLOGY
In this section the techniques used to solve the trajectory optimization problem are presented.The first approach was based on Silva (1995).The second one describes a hybrid algorithm that merges the direct and indirect methods.

FIRST APPROACH -DIRECT METHOD
The method applied within the framework of this approach is based on Silva (1995).A polynomial control function is used to model the flight profile.Four parameters, namely the coast time duration t coast (if it is applied) and three parameters of the polynomial control function, are optimized in order to get the maximum payload mass.A code from Jacob (1972) written in FORTRAN is transcript to C++ language and adapted to solve the problem (Eq.24).
The trajectory optimization problem can be formulated as (Eq.25): which maximize the pay load mass subjects to the constraints at orbit injection.Thus, for the first formulation of equations of motion (Eqs.26-28): 1 parameters (-).
The trajectory optimization problem can be formulated as (Eq.25):

Second Approach -Hybrid Method
The idea behind the strategy is to divide the flight trajectory into two distinct phases, one while the vehicle ascents the dense atmosphere and the other one when the vehicle is virtually in vacuum space, i.e., where aerodynamic effects can be ignored.In the aforementioned sections it was stated that the sensibility of the indirect shooting or even multiple shooting method depends on the initial guess, however, it is possible to where: Re = equatorial radius of the Earth (km).

Second Approach -Hybrid Method
The idea behind the strategy is to divide the flight trajectory into two distinct phases, one while the vehicle ascents the dense atmosphere and the other one when the vehicle is virtually in vacuum space, i.e., where aerodynamic effects can be ignored.In the aforementioned sections it was stated that the sensibility of the indirect shooting or even multiple shooting method depends on the initial guess, however, it is possible to

Second Approach -Hybrid Method
The idea behind the strategy is to divide the flight trajectory into two distinct phases, one while the vehicle ascents the dense atmosphere and the other one when the vehicle is virtually in vacuum space, i.e., where aerodynamic effects can be ignored.In the aforementioned sections it was stated that the sensibility of the indirect shooting or even multiple shooting method depends on the initial guess, however, it is possible to subject to the constraints at orbit injection.Thus, for the first formulation of equations of motion (Eqs.26-28): (13) ( 14) (28) where: Re = equatorial radius of the Earth (km).

Second Approach -Hybrid Method
The idea behind the strategy is to divide the flight trajectory into two distinct phases, one while the vehicle ascents the dense atmosphere and the other one when the vehicle is virtually in vacuum space, i.e., where aerodynamic effects can be ignored.In the aforementioned sections it was stated that the sensibility of the indirect shooting or even multiple shooting method depends on the initial guess, however, it is possible to subject to the constraints at orbit injection.Thus, for the first formulation of equations of motion (Eqs.26-28): (13) (14) (28) where: Re = equatorial radius of the Earth (km).

Second Approach -Hybrid Method
The idea behind the strategy is to divide the flight trajectory into two distinct phases, one while the vehicle ascents the dense atmosphere and the other one when the vehicle is virtually in vacuum space, i.e., where aerodynamic effects can be ignored.In the aforementioned sections it was stated that the sensibility of the indirect shooting or even multiple shooting method depends on the initial guess, however, it is possible to of motion (Eqs.26-28): (13) (14) (28) where: Re = equatorial radius of the Earth (km).

Second Approach -Hybrid Method
The idea behind the strategy is to divide the flight trajectory into two distinct phases, one while the vehicle ascents the dense atmosphere and the other one when the vehicle is virtually in vacuum space, i.e., where aerodynamic effects can be ignored.In the aforementioned sections it was stated that the sensibility of the indirect shooting or even multiple shooting method depends on the initial guess, however, it is possible to compute the initial co-states variable for optimal thrust arcs in vacuum fairly easy using almost arbitrary initial guess when direct and indirect methods are merged.

Stages Before Coasting -Direct Method
To fulfill this task, the method used in the first method is applied here.In other words, the same control variables will be obtained until the beginning of the coast arc.
Furthermore, this step gives the gross lift-off mass (GLOW) of the vehicle.

Upper Stage Trajectory -Indirect method
In this phase, the vehicle is assumed to be out of the denser layers of the atmosphere, so that the aerodynamics forces can be neglected.Thus, the trajectory of the upper stage is accomplished in the orbit plane.The upper stage flight is divided into two phases: a coast arc and a thrust arc.The equations of motion are presented below (Eqs.32-34) (Schlingloff 1987).
compute the initial co-states variable for optimal thrust arcs in vacuum fairly easy using almost arbitrary initial guess when direct and indirect methods are merged.

Stages Before Coasting -Direct Method
To fulfill this task, the method used in the first method is applied here.In other words, the same control variables will be obtained until the beginning of the coast arc.
Furthermore, this step gives the gross lift-off mass (GLOW) of the vehicle.

Upper Stage Trajectory -Indirect method
In this phase, the vehicle is assumed to be out of the denser layers of the atmosphere, so that the aerodynamics forces can be neglected.Thus, the trajectory of the upper stage is accomplished in the orbit plane.The upper stage flight is divided into two phases: a coast arc and a thrust arc.The equations of motion are presented below (Eqs.32-34) (Schlingloff 1987).
where: R e = equatorial radius of the Earth (km).

SECOND APPROACH -HYBRID METHOD
The idea behind the strategy is to divide the flight trajectory into two distinct phases, one while the vehicle ascents the dense atmosphere and the other one when the vehicle is virtually in vacuum space, i.e., where aerodynamic effects can be ignored.In the aforementioned sections it was stated that the sensibility of the indirect shooting or even multiple shooting method depends on the initial guess, however, it is possible to compute the initial co-states variable for optimal thrust arcs in vacuum fairly easy using almost arbitrary initial guess when direct and indirect methods are merged.

STAGES BEFORE COASTING -DIRECT METHOD
To fulfill this task, the method used in the first method is applied here.In other words, the same control variables will be obtained until the beginning of the coast arc.Furthermore, this step gives the gross lift-off mass (GLOW) of the vehicle.

UPPER STAGE TRAJECTORY -INDIRECT METHOD
In this phase, the vehicle is assumed to be out of the denser layers of the atmosphere, so that the aerodynamics forces can be neglected.Thus, the trajectory of the upper stage is accomplished in the orbit plane.The upper stage flight is divided into two phases: a coast arc and a thrust arc.The equations of motion are presented below (Eqs.32-34) (Schlingloff 1987).

xx/xx 08/15
Applying the Lagrange method, the Hamiltonian H is constructed as (Eq.35): Applying the Lagrange method, the Hamiltonian H is constructed as (Eq.35): (35) The adjoint equations for the co-state (Euler-Lagrange equations) are (Eqs.36- Applying the Lagrange method, the Hamiltonian H is constructed as (Eq.35): (35) The adjoint equations for the co-state (Euler-Lagrange equations) are (Eqs.36- Using the Pontryagin minimum principle, the optimal thrust angle control can be expressed as a function of the co-states (Eq.39): (39) and finally (Eq.40): (40) Thus, to get the optimal trajectory, Eqs.32-34 along with Eqs.36-38 should be integrated.However, as the Lagrange multipliers have no physical meaning and the flight path depends very sensitively on the initial guesses, it is difficult to solve these equations.Schlingloff (1987) developed an analytical method to eliminate the Lagrange Using the Pontryagin minimum principle, the optimal thrust angle control can be expressed as a function of the co-states (Eq.39): (39) and finally (Eq.40): (40) Thus, to get the optimal trajectory, Eqs.32-34 along with Eqs.36-38 should be integrated.However, as the Lagrange multipliers have no physical meaning and the flight path depends very sensitively on the initial guesses, it is difficult to solve these equations.Schlingloff (1987) developed an analytical method to eliminate the Lagrange Using the Pontryagin minimum principle, the optimal thrust angle control can be expressed as a function of the co-states (Eq.39): (39) and finally (Eq.40): (40) Thus, to get the optimal trajectory, Eqs.32-34 along with Eqs.36-38 should be integrated.However, as the Lagrange multipliers have no physical meaning and the flight path depends very sensitively on the initial guesses, it is difficult to solve these equations.Schlingloff (1987) developed an analytical method to eliminate the Lagrange multipliers getting formulas that can be represented by a smaller number of variables.
Using the Pontryagin minimum principle, the optimal thrust angle control can be expressed as a function of the co-states (Eq.39): (39) and finally (Eq.40): (40) Thus, to get the optimal trajectory, Eqs.32-34 along with Eqs.36-38 should be integrated.However, as the Lagrange multipliers have no physical meaning and the flight path depends very sensitively on the initial guesses, it is difficult to solve these equations.Schlingloff (1987) developed an analytical method to eliminate the Lagrange multipliers getting formulas that can be represented by a smaller number of variables.Thus, defining z (Eq.41): Using the Pontryagin minimum principle, the optimal thrust angle control can be expressed as a function of the co-states (Eq.39): (39) and finally (Eq.40): (40) Thus, to get the optimal trajectory, Eqs.32-34 along with Eqs.36-38 should be integrated.However, as the Lagrange multipliers have no physical meaning and the flight path depends very sensitively on the initial guesses, it is difficult to solve these equations.Schlingloff (1987) developed an analytical method to eliminate the Lagrange multipliers getting formulas that can be represented by a smaller number of variables.Thus, defining z (Eq.41): where: z = variable with no physical meaning.Schlingloff (1987)   Using the Pontryagin minimum principle, the optimal thrust angle control can be expressed as a function of the co-states (Eq.39): and finally (Eq.40): Thus, to get the optimal trajectory, Eqs.32-34 along with Eqs.36-38 should be integrated.However, as the Lagrange multipliers have no physical meaning and the flight path depends very sensitively on the initial guesses, it is difficult to solve these equations.Schlingloff (1987) developed an analytical method to eliminate the Lagrange multipliers getting formulas that can be represented by a smaller number of variables.Thus, defining z (Eq.41): where: z = variable with no physical meaning.Schlingloff (1987) got alternative differential equations to represent the control law (Eqs.42 and 43): xx/xx 09/15 Thus, to get the new system of first-order differential equations, the control Eqs.42 and 43 must be joined to the equations of motion Eqs.32-34.To integrate this system, since the initial condition for the state equations are fixed, just the initial guesses to the control equations must be set.To take into account the coast, the optimization problem can be stated as: Find X = [β 0 z 0 t coast ] T which minimize the propellant mass of the upper stage F(X) = m prop .It implies to maximize the payload mass that can be injected into the desired orbit.The constraints at orbit injection are given by Eqs.26-28.To solve this problem we can use heuristic methods such as particle swarm optimization (PSO).

PROGRAM STRUCTURE
As aforementioned stated, a modular approach using object-oriented programming (OOP) is chosen and to allow a better visualization of the codes it is used UML diagrams.UML diagrams are used to visualize the code and the communication between objects enabling a high degree of abstraction.Figure 3 presents a UML diagram for a specific launch vehicle, namely the VLS-alfa launch vehicle.From the diagram we can see some parameters and functions of each component and the relationship between them.In order to make the diagram clearer, some parameters and functions are omitted.

RESULTS
To verify the trajectory program, two launch vehicles are considered, namely: the Brazilian VLS-1 and the European Ariane 5.Both mathematical modeling of the ascent trajectory are used in order to verify the concordance between them.The trajectory optimization will be performed using direct method.

VLS-1 LAUNCH VEHICLE
The under development VLS-1 is the future Brazilian satellite launch vehicle.Its development started in 1984, however due to technical problems, the vehicle could not be qualified up to now.VLS-1 is composed of four solid stages.The first stage is equipped with four solid boosters S43 (Fig. 4).The vehicle is designed to perform a non-powered cost arc between third and upper stages.Key parameters of the vehicle used in the simulation are given in Fig. 4. The mission is to launch a satellite into a reference circular orbit of 500 km of altitude from the Alcântara Launch Center (2°22'39.52"S, 44°23'57.71"W).Figures 5 and 6 show the velocity and altitude profiles with the powered phase described by the red curve, and the non-powered phase by the blue curve.A path constraint (dynamic pressure) of the flight is shown in Fig. 7.The trajectory profile of the VLS-1 was verified with the SKYNAV tool (Schlingloff 1991) and presented a good agreement.Tables 1 and 2 summarize the corner instants between flight phases and Table 3 presents the control parameters obtained by optimization subroutine.As expected, the maximum payload mass found for both formulations are very close.Comparing the results from the Tables 1 to 3 and Figs. 5 to 7, we can verify that indeed both mathematical modeling are equivalent.Although the method presented in this work gives sub-optimal trajectory, for the purpose of a preliminary analysis, this method is sufficiently accurate.

ARIANE 5 LAUNCH VEHICLE
Built under supervision of European Space Agency (ESA), Ariane 5 is a European launch vehicle that is part of the Ariane rocket family.The vehicle is used to deliver payload into low Earth orbit (LEO) and geostationary transfer orbit (GTO).Within the framework of this work, the mission is to launch a satellite from Kourou to a low Earth orbit (LEO) of 200 km of altitude.The key parameters of the vehicle are given in Table 4.The trajectory can be divided into two main phases.To take the vehicle from the ground and minimize the gravitational losses, the first phase is powered by two solid boosters and by the core stage using the propellants LOX/LH2.After 146.84 seconds the booster stages are decoupled from the vehicle and the motion is powered only by the core stage.Differently from the VLS-1, this launch vehicle does not perform a non-powered coast-arc.Figure 8 presents the altitude and relative velocity profiles.
From Fig. 8 we can notice that the behavior of the altitude profile firstly exceeded the desired altitude (LEO with 200 km) reaching a maximum altitude, and then the desired altitude is obtained.The reason for that is presumably because the flight does not perform a non-powered phase, so the vehicle takes longer to get the right inclination in order to injected into orbit.

VLS-ALFA LAUNCH VEHICLE
It is known that the VLS-Alfa will replace the last two solid stages of the former VLS-1 by a single liquid upper stage.Then since the VLS-Alfa is an improvement of the former VLS-1, the VLS-1 will be used as a reference vehicle to perform the simulations.The upper stage of the VLS-Alfa presumably will perform a coast phase, so the L75 is supposed to support restart capability.The mission is to launch a satellite into a reference circular orbit of 500 km of altitude from the Alcântara Launch Center (2°22'39.52"S, 44°23'57.71"W).The parameters of the vehicles are given in Tables 5 and 6.As the mission of the VLS-Alfa is still not totally defined, the propellant mass of the upper stage had to be estimated.Thus, an amount of 6900 kg was estimated based on internal reports.From this value an amount of 1100 kg was taken for the phase after coasting, i.e., for orbit injection.The altitude and relative velocity profiles for both vehicles are presented in Fig. 9.

CONCLUSIONS
This article presented a modeling of a launch vehicle using object-oriented programming.Two models for trajectory optimization were presented.To allow a better visualization of the codes, the Unified Modeling Language (UML) was used.The future Brazilian VLS-Alfa and its predecessor VLS-1 were simulated and compared.Both of them are supposed to support perform coast phase.The performance of the Ariane 5, which does not perform coast phase, was presented and the results showed good agreement with the literature.Since the stages of the launch vehicle is an assembly engine and propellant system, the modular approach is interesting to allow future integration of the propulsion system in order to study the coupling among the disciplines.Accordingly, the modular approach using object-oriented programming can facilitate a multidisciplinary optimization. 2 = thrust angle in flight plane (rad); δ = thrust angle out of flight plane (rad); and F = thrust force (N); ω x = inclination change (rad).

Figure 2 .
Figure 2. Planet-fixed and local horizon frames for atmospheric flight (adapted from Tewari 2007).
the relative velocity v and the local velocity of the local horizontal frame (oxyz) relative to the planet-centered rotating frame (SXYZ) can be expressed as (Eqs.10 and 11): 13 and 14 we finally obtain the kinematic equations of motion (Eqs.15-17): derive the dynamic equations Newton's Second Law must be introduced (Eq.18): (18) Choosing the wind axes to express the forces on the body and doing the appropriate transformation to perform aI in the wind axes, the remaining equations to model the translational motion are obtained (Eqs.19-21): (19) β1, β2, β3 = set of optimization control parameters (-).
constraints at orbit injection.Thus, for the first formulation of equations of motion (Eqs.26-28): Re = equatorial radius of the Earth (km).and for the second one (Eqs.29-31): constraints at orbit injection.Thus, for the first formulation of equations of motion (Eqs.26-28): Re = equatorial radius of the Earth (km).and for the second one (Eqs.29-31): got alternative differential equations to represent the control law (Eqs.42 and 43): get the new system of first-order differential equations, the control Eqs.
get the new system of first-order differential equations, the control Eqs.42 and 43 must be joined to the equations of motion Eqs.32-34.To integrate this for the co-state (Euler-Lagrange equations) are (Eqs.36-38):

Figure 3 .
Figure 3. Schematic of a UML diagram of the VLS-Alfa launch vehicle.

Figure 5 .
Figure 5. Velocity profile of the VLS-1 launch vehicle using (a) first formulation of state equations and (b) the second formulation.

Figure 6 .
Figure 6.Altitude profile of the VLS-1 launch vehicle using (a) first and (b) the second formulation of state equations.

Figure 7 .
Figure 7. Dynamic pressure profile of the VLS-1 launch vehicle using (a) first and (b) the second formulation of state equations.

Figure 8 .
Figure 8.(a) Altitude profile of the Ariane 5 and (b) velocity profile of the Ariane.

Figure 9 .
Figure 9. (a) Relative velocity profile of the former and future Brazilian launch vehicle and (b) altitude profile of the former and future Brazilian launch vehicle.

Table 1 .
Values for state variables at start, inter-stage and end instants.

Table 2 .
Values for state variables at start, inter-stage and end instants.