Ascent Trajectories of Multistage Launch Vehicles: Numerical Optimization with Second-Order Conditions Verification

Abstract

Multistage launch vehicles are employed to place spacecraft and satellites in their operational orbits. Trajectory optimization of their ascending path is aimed at defining the maximum payload mass at orbit injection, for specified structural, propulsive, and aerodynamic data. This work describes and applies a method for optimizing the ascending path of the upper stage of a specified launch vehicle through satisfaction of the necessary conditions for optimality. The method at hand utilizes a recently introduced heuristic technique, that is, the particle swarm algorithm, to find the optimal ascent trajectory. This methodology is very intuitive and relatively easy to program. The second-order conditions, that is, the Clebsch-Legendre inequality and the conjugate point condition, are proven to hold, and their fulfillment enforces optimality of the solution. Availability of an optimal solution to the second order is an essential premise for the possible development of an efficient neighboring optimal guidance.

1. Introduction

Multistage rockets are employed to place spacecraft and satellites in their operational orbits. The optimization of their ascending trajectory leads to determining the maximum payload mass that can be inserted in the desired orbit. This goal is achieved by finding the optimal control time history and the optimal thrust and coast durations.

The numerical solution of aerospace trajectory optimization problems has been pursued with different approaches in the past. Indirect methods, such as the gradient-restoration algorithm [1, 2] and the shooting method [3] or direct techniques, such as direct collocation [4, 5], direct transcription [6, 7], and differential inclusion [8, 9], are examples of such techniques. However, only a relatively small number of publications are concerned with trajectory optimization of multistage launch vehicles [1, 2, 10, 11]. A recently published paper [12] describes a simple method for performance evaluation through generation of a near optimal trajectory for a multistage launch vehicle.

This research considers the optimal exoatmospheric trajectory of the upper stage of the latter rocket, whose characteristics are specified. The trajectory arc that precedes orbital injection is composed of two phases: (1) coast (Keplerian) arc and (2) thrust phase. More specifically, for the upper stage the existence and duration of a coast arc (with no propulsion) and the optimal thrust direction are being investigated through the first-order necessary conditions for optimality, that is, the Euler-Lagrange equations and the Pontryagin minimum principle, in conjunction with a heuristic technique, that is, the particle swarm algorithm. This technique is intended to numerically determine the state, control, and adjoint variables (conjugate to the dynamics equations). The availability of the latter quantities allows the direct evaluation of the second-order conditions for optimality, that is, the Clebsch-Legendre condition and the nonexistence of conjugate points. The numerical process for verifying these conditions is based on the sweep method and is illustrated in full detail in this work.

The present paper is organized as follows. Section 2 deals with fundamental concepts related to the multistage rocket characteristics and dynamics. Section 3 is concerned with optimal thrust programming and illustrates the first- and second-order conditions for optimality. Section 4 contains an overview of the numerical technique used to generate the optimal ascending path and presents the numerical solution to the trajectory optimization problem, together with the second-order conditions verification. Section 5 concludes the paper, summarizing its main results and accomplishments.

2. Characteristics and Dynamics of the Multistage Launch Vehicle

This work addresses the problem of optimizing the ascending path of the upper stage of a multistage rocket with specified structural, propulsive, and aerodynamic characteristics, through the determination of the optimal control time history, leading to injecting the payload in the desired operational orbit. This is a circular orbit of radius 𝑅_𝑑 and specified spatial orientation. The launch vehicle is modeled as a point mass, in the context of a three-degree-of-freedom problem.

2.1. Rocket Characteristics

The three-stage rocket that is being considered is the MultiRole Air Launch Missile (Muralm) [13]. It has specified structural, propulsive, and aerodynamic characteristics, which are the same as those described in [12] and are (partially) reported in the following.

For the sake of simplicity, the mass distribution of the launch vehicle can be described in terms of masses of subrockets: subrocket 1 is the entire rocket, including all the three stages, subrocket 2 is the launch vehicle after the first stage separation, and subrocket 3 is the launch vehicle after separation of the first two stages and therefore is represented by the third stage only. Let 𝑚₀^(𝑖) denote the initial mass of subrocket 𝑖. This mass 𝑚₀^(𝑖) is composed of a structural mass 𝑚_𝑆^(𝑖), a propellant mass 𝑚_𝑃^(𝑖), and a payload mass 𝑚_𝑈^(𝑖) (𝑚₀^(𝑖)=𝑚_𝑆^(𝑖)+𝑚_𝑃^(𝑖)+𝑚_𝑈^(𝑖)). For the first two subrockets, 𝑚_𝑈^(𝑖) (𝑖=1,2) coincides with the initial mass of the subsequent subrocket (i.e., 𝑚_𝑈^(𝑖)≡𝑚₀^(𝑖+1)). With regard to the third subrocket, 𝑚_𝑈⁽³⁾ is the (actual) payload mass, which is to be maximized. The mass distribution for the Muralm three-stage rocket is the following: 𝑚₀⁽¹⁾=3608kg𝑚_𝑆⁽¹⁾=306kg𝑚_𝑃⁽¹⁾𝑚=2480kg_𝑈(1)=822kg𝑚₀⁽²⁾=822kg𝑚_𝑆⁽²⁾𝑚=46kg_𝑃(2)=370kg𝑚_𝑈⁽²⁾=406kg≡𝑚₀⁽³⁾. (1) The initial mass of the third stage is specified and this implies that minimizing the third stage propellant expenditure implies maximizing the payload mass that can be injected in the desired orbit.

The propulsive characteristics of the launch vehicle can be described in terms of thrust magnitude 𝑇 and specific impulse 𝐼_SP. Both these quantities are assumed as time-independent for all of the stages. For the upper stage, the following propulsive data are assumed: 𝑇⁽³⁾=7500N,𝐼⁽³⁾_SP=320sec, (2) where the superscript refers to stage 3. For each subrocket 𝑖, it is then straightforward to derive the thrust acceleration as a function of the time 𝑡: 𝑇^(𝑖)𝑚^(𝑖)=𝑛₀^(𝑖)1−𝑛₀^(𝑖)𝑡−𝑡^(𝑖)_𝑖𝑔𝑛/𝑔₀𝐼^(𝑖)_SP, (3) where 𝑛₀^(𝑖) is the thrust acceleration at ignition of stage 𝑖 (which occurs at the time 𝑡^(𝑖)_𝑖𝑔𝑛) and 𝑔₀ is the gravitational acceleration at sea level.

The aerodynamics of the Muralm rocket has been modeled through the Missile DATCOM database for the first two stages, whereas it is irrelevant for the upper stage, because the flight is approximated as entirely exoatmospheric, because of the altitude reached at the second stage separation.

2.2. Equations of Motion

Orbital motion is usually described in an Earth-centered inertial frame (ECI), with origin coincident with the Earth center, and axes associated with the unit vectors (𝑐₁̂,𝑐₂̂,𝑐₃̂); 𝑐₁̂ corresponds to the vernal axis, whereas 𝑐₃̂ is aligned with the Earth rotation axis and directed toward North (cf. Figure 1). The orbital plane contains the unit vectors 𝑁, associated with the ascending node crossing (corresponding to the point at which the space vehicle crosses the equatorial plane from South to North), 𝑒̂, aligned with the direction of perigee, and 𝑟̂, directed toward the instantaneous position of the spacecraft. Two angles suffice to define the orbit plane orientation in the ECI frame, that is, the right ascension of the ascending node (RAAN) Ω(0≤Ω<2𝜋) and the orbit inclination 𝑖(0≤𝑖≤𝜋). The argument of perigee 𝜔 defines the position of perigee on the orbit plane, whereas the instantaneous true anomaly 𝑓(𝑡) identifies the position of the space vehicle.

Figure 1 

ECI-frame and geometry of orbits.

The entire trajectory of the upper stage entirely belongs to the desired orbital plane, which has been selected during the ascending trajectory of the first two stages (cf. [12]). This circumstance implies that the orbit inclination 𝑖 and RAAN Ω do not vary during the last stage propulsion. In contrast, the remaining orbit elements, that is, the orbit semimajor axis 𝑎, eccentricity 𝑒, and argument of perigee 𝜔, change up to reaching the desired final values (denoted with the subscript “𝑑”); that is, 𝑎_𝑑=𝑅_𝑑 and 𝑒_𝑑=0 (𝜔 is not defined for circular orbits).

At separation of the second stage, the previously mentioned orbit elements are known and denoted with 𝑎₂, 𝑒₂, and 𝜔₂. The third stage trajectory is assumed to be composed of two phases: a coast arc and a thrust phase. It is relatively straightforward to prove that during the coast arc the true anomaly variation, denoted with Δ𝑓, suffices to describe the rocket dynamics. In fact, if 𝑡₂ and 𝑡₃ represent, respectively, the time of the second stage separation and the ignition time of the upper stage, then 𝑓₃≜𝑓(𝑡₃)=𝑓(𝑡₂)+Δ𝑓=𝑓₂+Δ𝑓. The orbital elements 𝑎₂, 𝑒₂, and 𝜔₂ do not vary along the coast (Keplerian) arc. As a result, the radius 𝑟, inertial velocity 𝑣_𝐼, and inertial flight path angle 𝛾_𝐼 at 𝑡₃ are simply given by 𝑟𝑡₃−𝑎₂1−𝑒²₂1+𝑒₂cos𝑓₃𝛾=0,_𝐼𝑡₃𝑒−arctan²sin𝑓₃1+𝑒₂cos𝑓₃𝑣=0,_𝐼𝑡₃−𝜇_𝐸𝑎₂1−𝑒²₂1+𝑒²₂+2𝑒₂cos𝑓₃=0, (4)(5)(6) where 𝜇_𝐸(=398600.4km³/sec²) denotes the Earth gravitational parameter. In the powered phase, the third stage motion can be conveniently described through the use of the equations of motion for 𝑟, 𝑣_𝐼, and 𝛾_𝐼: 𝑟̇=𝑣_𝐼sin𝛾_𝐼,𝑣̇_𝐼=𝑇⁽³⁾𝑚⁽³⁾cos𝛼_𝐼𝑣_𝐼−𝜇_𝐸𝑟²sin𝛾_𝐼,𝛾̇_𝐼=𝑇⁽³⁾𝑚⁽³⁾sin𝛼_𝐼𝑣_𝐼+𝑣_𝐼𝑟−𝜇_𝐸𝑟²𝑣_𝐼cos𝛾_𝐼, (7) where the subscript “3” will be omitted henceforth and can be written in compact form as state equations: 𝐱̇=𝐟(𝐱,𝐮,𝐚,𝑡). (8) The thrust direction is specified through the thrust angle 𝛼_𝐼(≡𝐮), which is the only time-varying control variable and is taken clockwise from the direction of the inertial velocity vector 𝐯_𝐼. The operational orbit is assumed to be a circular orbit of radius 𝑅_𝑑. This means that the final conditions at orbit injection (which occurs at 𝑡_𝑓) are 𝑟𝑡_𝑓−𝑅_𝑑=0𝑣_𝐼𝑡_𝑓−𝜇_𝐸𝑅_𝑑=0𝛾_𝐼𝑡_𝑓=0. (9) The left-hand sides of (4) through (6) and (9) are collected in the vector 𝝍(𝐱₀,𝐱_𝑓,𝐚,𝑡₀,𝑡_𝑓)(=𝟎)(𝑡₀≡𝑡₃).

3. Optimal Thrust Programming

The trajectory of the upper stage is obtained by minimizing fuel consumption, which is equivalent to minimizing the thrust duration. Hence, the objective function for the third stage trajectory optimization is 𝐽=𝑡_𝑓−𝑡₃. (10) Hence, the optimization problem for the upper stage is the following: determine the optimal control law 𝛼_𝐼(𝑡) and the optimal true anomaly 𝑓₃ such that 𝐽 is minimized. The ignition time 𝑡₃ is then calculated through the Kepler’s law, 𝑡₃=𝑡₂+𝑎³₂𝜇_𝐸×𝐸𝑡₃−𝐸𝑡₂−𝑒₂sin𝐸𝑡₃+𝑒₂sin𝐸𝑡₂, (11) where 𝐸(𝑡_𝑘) is the eccentric anomaly associated with 𝑓(𝑡_𝑘)(𝑘=2,3).

Letting 𝐱≜[𝑥₁𝑥₂𝑥₃]^𝑇=[𝑟𝑣_𝐼𝛾_𝐼]^𝑇, to obtain the necessary conditions for an optimal solution a Hamiltonian 𝐻 and a function of the boundary conditions Φ are introduced as 𝐻≜𝝀^𝑇𝐟=𝜆₁𝑥₂sin𝑥₃+𝜆₂𝑇⁽³⁾𝑚⁽³⁾cos𝛼_𝐼𝑥₂−𝜇_𝐸𝑥²₁sin𝑥₃+𝜆₃𝑇⁽³⁾𝑚⁽³⁾sin𝛼_𝐼𝑥₂+𝑥₂𝑥₁−𝜇_𝐸𝑥²₁𝑥₂cos𝑥₃,Φ≜𝐽+𝝊^𝑇𝝍=𝑡_𝑓−𝑡₃+𝜐₁𝑥₁₀−𝑎₂1−𝑒²₂1+𝑒₂cos𝑓₃+𝜐₂𝑥₂₀−𝜇_𝐸𝑎₂1−𝑒²₂1+𝑒²₂+2𝑒₂cos𝑓₃+𝜐₃𝑥₃₀𝑒−arctan²sin𝑓₃1+𝑒₂cos𝑓₃+𝜐₄𝑥_1𝑓−𝑅_𝑑+𝜐₅𝑥_2𝑓−𝜇_𝐸𝑅_𝑑+𝜐₆𝑥_3𝑓, (12) where 𝑥_𝑘0=𝑥_𝑘(𝑡₃) and 𝑥_𝑘𝑓=𝑥_𝑘(𝑡_𝑓)(𝑘=1,2,3); 𝝀(≜[𝜆₁𝜆₂𝜆₃]^𝑇) and 𝝊(≜[𝜐₁𝜐₂𝜐₃𝜐₄𝜐₅𝜐₆]^𝑇) represent, respectively, the adjoint variable conjugate to the dynamics equations (7) and to the boundary conditions (4)–(6) and (9). With these definitions, the objective function assumes its extended form 𝐽^′=Φ𝐱₀,𝐱_𝑓,𝑓₃,𝑡₀,𝑡_𝑓+^𝑡𝑓_𝑡0𝐻(𝐱,𝐮,𝑡)−𝝀^𝑇𝐱̇𝑑𝑡. (13)

3.1. First-Order Conditions

The first-order necessary conditions for (local) optimality [14] include the adjoint equations: 𝝀̇=−𝜕𝐻𝜕𝐱^𝑇, (14) which yield the following scalar relations for the costate 𝝀(≜[𝜆₁𝜆₂𝜆₃]^𝑇): 𝜆̇₁=𝑥₂𝜆₃cos𝑥₃1𝑥²₁−2𝜇_𝐸𝜆₂sin𝑥₃+2𝜇_𝐸𝜆₃cos𝑥₃𝑥₂1𝑥³₁,𝜆̇₂=−𝜆₁sin𝑥₃−𝜆₃cos𝑥₃1𝑥₁+𝜇_𝐸𝑥²₁𝑥²₂−𝑇⁽³⁾𝑚⁽³⁾1𝑥²₂sin𝛼_𝐼,𝜆̇₃=−𝑥₂𝜆₁cos𝑥₃+𝜇_𝐸𝜆₂cos𝑥₃𝑥²₁+𝜆₃sin𝑥₃𝑥₂𝑥₁−𝜇_𝐸𝑥²₁𝑥₂, (15)(16)(17) in conjunction with the respective boundary conditions: 𝝀₀=−𝜕Φ𝜕𝐱₀^𝑇⟹𝜆_𝑘0=−𝜐_𝑘()𝝀𝑘=1,2,3_𝑓=𝜕Φ𝜕𝐱_𝑓^𝑇⟹𝜆_𝑘𝑓=𝜐_𝑘+3(𝑘=1,2,3), (18)(19) where the subscripts “0” and “𝑓 ” refer to 𝑡₃ and 𝑡_𝑓, respectively. In the presence of initial conditions depending on a parameter (𝐚≡𝑓₃) a pair of additional necessary conditions must hold [15]: Φ_𝐚=0⟹𝜕Φ𝜕𝑓₃=0,Φ_𝐚𝐚𝜕≥0⟹²Φ𝜕𝑓²₃≥0. (20) The first equation yields a relation that expresses 𝜆₃₀ as a function of 𝜆₁₀, 𝜆₂₀, and 𝑓₃, 𝜆₃₀=𝜆₂₀sin𝑓₃1+𝑒²₂+2𝑒₂cos𝑓₃𝑒₂+cos𝑓₃−𝜆₁₀𝑎₂1−𝑒²₂𝑒₂+cos𝑓₃1+𝑒²₂+2𝑒₂cos𝑓₃1+𝑒₂cos𝑓₃². (21) In addition, the optimal control 𝛼^∗_𝐼 can be expressed as a function of the costate through the Pontryagin minimum principle: 𝛼^∗_𝐼=argmin_𝛼𝐼𝐻⟹sin𝛼^∗_𝐼𝜆=−³𝑥₂𝜆₃𝑥₂²+𝜆²₂^−1/2,cos𝛼^∗_𝐼=−𝜆₂𝜆₃𝑥₂²+𝜆²₂^−1/2, (22) which implies the fulfillment of the stationarity condition of 𝐻 with respect to 𝐮; that is, 𝐻^𝑇_𝐮=𝟎. (23) Lastly, as the final time is unspecified, the following transversality condition must hold: 𝐻𝑡_𝑓+𝜕Φ𝜕𝑡_𝑓𝑛=0⟹⁰⁽³⁾1−𝑛₀⁽³⁾𝑡_𝑓−𝑡₃/𝑔₀𝐼⁽³⁾_SP×𝜆_3𝑓𝑥_2𝑓²+𝜆_2𝑓²−1=0. (24) The necessary conditions for optimality allow translating the optimal control problem into a two-point boundary-value problem involving (15)–(24), with unknowns represented by the initial values of 𝝀, 𝑓₃, and 𝑡_𝑓.

However, the transversality condition can be proven to be ignorable. First, homogeneity of the costate equations (15)–(17) can be easily recognized. This circumstance implies that if an optimization algorithm is capable of finding some initial value of 𝝀 such that 𝜆₁(0)=𝑘_𝜆𝜆^∗₁(0), 𝜆₂(0)=𝑘_𝜆𝜆^∗₂(0), and 𝜆₃(0)=𝑘_𝜆𝜆^∗₃(0)(𝑘_𝜆>0) (where the superscript “∗” denotes the actual optimal value of a variable), then the same proportionality holds between 𝝀 and 𝝀^∗ at any 𝑡, due to homogeneity of (15)–(17). Second, if the control 𝛼_𝐼 is written as a function of 𝝀(=𝑘_𝜆𝝀^∗) through (22), then one can recognize that 𝛼_𝐼 coincides with 𝛼^∗_𝐼: sin𝛼_𝐼𝜆=−³𝑥₂𝜆₃𝑥₂²+𝜆²₂^−1/2𝑘=−^𝜆𝜆^∗₃𝑥₂𝑘_𝜆𝜆^∗₃𝑥₂²+𝑘_𝜆𝜆^∗₂^2−1/2≡sin𝛼^∗_𝐼,cos𝛼_𝐼=−𝜆₂𝜆₃𝑥₂²+𝜆²₂^−1/2=−𝑘_𝜆𝜆^∗₂𝑘_𝜆𝜆^∗₃𝑥₂²+𝑘_𝜆𝜆^∗₂^2−1/2≡cos𝛼^∗_𝐼. (25) This circumstance implies that if 𝝀 is proportional to 𝝀^∗, then the final conditions are fulfilled at the minimum final time 𝑡^∗_𝑓 even when the control is expressed as a function of the nonoptimal 𝝀. In contrast, using 𝝀(=𝑘_𝜆𝝀^∗) the transversality condition is violated, because the value of 𝐻(𝑡^∗_𝑓), due to (24), turns out to be 𝐻𝑡^∗_𝑓𝑛=−⁰⁽³⁾1−𝑛₀⁽³⁾𝑡^∗_𝑓−𝑡^∗₃/𝑔₀𝐼⁽³⁾_SP×𝜆_3𝑓𝑥_2𝑓²+𝜆_2𝑓²=−𝑘_𝜆≠−1. (26) Hence, if the proportionality condition 𝝀=𝑘_𝜆𝝀^∗ holds, then optimal control 𝛼^∗_𝐼 can be determined without considering the transversality condition, which in fact is ignorable in this context.

3.2. Second-Order Conditions

The derivation of the second-order optimality conditions, which have local character as the first-order conditions, involves the definition of an admissible comparison path, located in the neighborhood of the (local) optimal path, associated with the optimal state 𝐱^∗, costate 𝝀^∗, and control 𝐮^∗. The nonexistence of alternative neighboring optimal solutions is to be proven in order to guarantee optimality of the solution of interest.

3.2.1. Weierstrass and Legendre-Clebsch Conditions

An admissible comparison path can be defined by introducing a comparison control 𝐮 such that 𝐮−𝐮^∗=𝟎if𝑡₀≤𝑡≤𝑡_𝑝𝐮−𝐮^∗=const.if𝑡_𝑝≤𝑡≤𝑡_𝑞𝐮−𝐮^∗suchthatΨ_𝐱𝑓𝛿𝐱_𝑓=𝟎if𝑡_𝑞≤𝑡≤𝑡_𝑓, (27) where 𝑡₀<𝑡_𝑝<𝑡_𝑞<𝑡_𝑓 and Δ𝑡(≜𝑡_𝑞−𝑡_𝑝) is small. Under the assumption that Ψ depends only on 𝐱_𝑓, the condition Ψ_𝐱𝑓𝛿𝐱_𝑓=𝟎 ensures fulfillment of the boundary conditions to first order (𝛿𝐱_𝑓 is the variation of the state at the final time). The comparison control is not located in the neighborhood of the optimal solution and for this reason the variation (𝐮−𝐮^∗) belongs to the class of strong variations. In contrast, the resulting state variation 𝛿𝐱(=𝐱−𝐱^∗) is reasonably small, and the comparison state 𝐱 is sufficiently close to the optimal state 𝐱^∗. Let 𝐽^′_∗ denote the local minimum value of the objective (associated with 𝐮^∗) and 𝐽^′ the value corresponding to 𝐮. After expanding Δ𝐽^′(=𝐽^′−𝐽^′_∗) in Taylor series, one obtains (cf. [15]) Δ𝐽^′=𝐻−𝐻_∗𝑝Δ𝑡, (28) where 𝐻 and 𝐻_∗ are associated with 𝐮 and 𝐮^∗. The quantity Δ𝐽^′ must be positive for any choice of 𝐮. Since 𝑡_𝑝 is arbitrary and Δ𝑡>0, the condition 𝐻−𝐻_∗>0 (29) must hold at any point along the minimizing path. Equation (29) represents the Weierstrass condition for a strong local minimum. For weak variations, that is, when 𝐮=𝐮^∗+𝛿𝐮 (with 𝛿𝐮 sufficiently small), (29) can be expanded in Taylor series about the optimal solution: 𝐻−𝐻_∗≃𝐻_𝑢∗1𝛿𝐮+2𝛿𝐮^𝑇𝐻_𝑢𝑢∗𝛿𝐮. (30) Along the optimal path (𝐻_𝑢)_∗=𝟎, this property and the fact that 𝛿𝐮 is arbitrary imply that the Hessian matrix (𝐻_𝑢𝑢)_∗ must be positive definite; that is, 𝐻_𝑢𝑢∗>0. (31) Equation (31) is the Clebsch-Legendre sufficient condition for a minimum; in the necessary (weak) form the sign “≥” replaces the inequality sign (i.e., the Hessian must be positive semidefinite).

3.2.2. Conjugate Point Condition

In general, a comparison path fulfills the state equations and boundary conditions to first order. This means that the state, costate, control, parameter, and final time displacements (𝛿𝐱,𝛿𝐮,𝑑𝐚,𝑑𝑡_𝑓) satisfy the linear equations (deriving from (8), (14), and (23)): 𝛿𝐱̇=𝐟_𝑥𝛿𝐱+𝐟_𝑢𝛿𝐮+𝐟_𝑎𝛿̇𝑑𝐚𝝀=−𝐻_𝑥𝑥𝛿𝐱−𝐻_𝑥𝑢𝛿𝐮−𝐻_𝑥𝜆𝛿𝝀−𝐻_𝑥𝑎𝐻𝑑𝐚_𝑢𝑥𝛿𝐱+𝐻_𝑢𝑢𝛿𝐮+𝐻_𝑢𝑎𝑑𝐚+𝐻_𝑢𝜆𝛿𝝀=𝟎 (32)(33)(34) accompanied by the boundary conditions (deriving from the boundary conditions (4)–(6), (9), (18), (19), and (24)): 𝝍_𝐱𝑓𝛿𝐱_𝑓+𝝍_𝐱0𝛿𝐱₀+𝝍_𝐚𝑑𝐚+𝝍^′𝑑𝑡_𝑓=𝟎𝛿𝝀₀=−Φ_𝐱0𝐱0𝛿𝐱₀−Φ_𝐱0𝐚𝑑𝐚−𝝍^𝑇_𝐱0𝑑𝝊𝛿𝝀_𝑓=Φ_{𝐱𝑓𝐱𝑓}𝛿𝐱_𝑓+Φ_𝐱𝑓𝐚𝑑𝐚+𝝍^𝑇_𝐱𝑓𝑑𝝊+Φ^′𝑇_𝐱𝑓𝑑𝑡_𝑓Λ_𝐱𝑓𝛿𝐱_𝑓+Λ_𝐱0𝛿𝐱₀+Λ_𝐚𝑑𝐚+𝝍^′𝑇𝑑𝝊+Λ^′𝑑𝑡_𝑓Φ=0_𝐚𝐱0𝛿𝐱₀+Φ_𝐚𝐱𝑓𝛿𝐱_𝑓+Φ_𝐚𝐚𝑑𝐚+Φ_𝐚𝑡𝑓𝑑𝑡_𝑓+𝝍_𝐚𝑑𝝊=𝟎, (35)(36)(37)(38)(39) where Λ≜𝐻_𝑓+Φ_𝑡𝑓. After solving (34) for 𝐮, 𝛿𝐮=−𝐻⁻¹_𝑢𝑢𝐻_𝑢𝑥𝛿𝐱+𝐻_𝑢𝑎𝑑𝐚+𝐻_𝑢𝜆𝛿𝝀, (40) and inserting the resulting expression in (32)-(33), one obtains 𝛿̇𝛿̇𝐱=𝐴𝛿𝐱−𝐵𝛿𝝀+𝐷𝑑𝐚𝝀=−𝐶𝛿𝐱−𝐴^𝑇𝛿𝝀−𝐸𝑑𝐚, (41)(42) where the matrices 𝐴, 𝐵, 𝐶, 𝐷, and 𝐸 are defined in terms of the quantities appearing in (32)–(34). The final conditions (35) and (37) justify the definition of the sweep variables, introduced through the following relations: 𝛿𝝀=𝑆𝛿𝐱+𝑅𝑑𝝊+𝑚𝑑𝐚𝟎=𝑅^𝑇𝛿𝐱+𝑄𝑑𝝊+𝑛𝑑𝐚𝟎=𝑚^𝑇𝛿𝐱+𝑛^𝑇𝑑𝝊+𝛼𝑑𝐚. (43)(44)(45) After inserting (43) through (45) into (32)-(33), one obtains 𝑆̇=−𝐶−𝐴^𝑇̇𝑆−𝑆𝐴+𝑆𝐵𝑆𝑅=𝑆𝐵−𝐴^𝑇𝑅𝑚̇=𝑆𝐵−𝐴^𝑇̇𝑚−𝑆𝐷−𝐸𝑄=𝑅^𝑇̇𝐵𝑅𝑛=𝑅^𝑇()̇𝐵𝑚−𝐷𝛼=𝑚^𝑇(𝐵𝑚−𝐷)−𝐷^𝑇𝑚−𝐹 (46) with final conditions that must be consistent with (35) through (39); that is, 𝑆_𝑓=Φ_{𝐱𝑓𝐱𝑓}𝑄_𝑓𝑅=0_𝑓=𝜓^𝑇_𝐱𝑓𝑛_𝑓=𝜓_𝐚𝑚_𝑓=Φ_𝐱𝑓𝐚𝛼_𝑓=Φ_𝐚𝐚. (47) Equations (44) and (45) can be solved simultaneously at 𝑡₀ to yield 𝑑𝝊𝑑𝐚=−𝑉₀⁻¹𝑈^𝑇₀𝛿𝐱₀, (48) where 𝑛𝑈=𝑅𝑚,𝑉=𝑄𝑛_𝑇𝛼. (49)(50) If (48) is used at 𝑡₀, then 𝛿𝝀₀=𝑆₀−𝑈₀𝑉₀⁻¹𝑈^𝑇₀𝛿𝐱₀. (51) Letting 𝑆≜𝑆−𝑈𝑉⁻¹𝑈^𝑇, it is relatively straightforward to demonstrate that the first of (46) holds also for 𝑆 (with 𝑆 in place of 𝑆), with boundary condition 𝑆→∞ as 𝑡→𝑡_𝑓. It is worth remarking that all the matrices introduced in this section are evaluated along the optimal path. From (51), it is then straightforward to conclude that 𝛿𝝀₀→𝟎 as 𝛿𝐱₀→0, unless 𝑆 tends to infinity at an internal time 𝑡(𝑡₀≤𝑡<𝑡_𝑓), which is referred to as conjugate point. If 𝛿𝝀₀→𝟎 and 𝛿𝐱₀→0, then also 𝛿𝐮→0, due to (40)–(42) and (48). This means that no neighboring optimal path exists, unless 𝑆→∞ at a certain time 𝑡(𝑡₀≤𝑡<𝑡_𝑓). In other words, the nonexistence of conjugate points is 𝑆<∞∀𝑡∈𝑡₀,𝑡_𝑓. (52) It is worth remarking that the numerical backward integration for 𝑆 can start by integrating (46) up to a time 𝑡_𝑠𝑤 sufficiently close to 𝑡_𝑓; then one can switch to 𝑆 for the remaining time interval [𝑡₀,𝑡_𝑠𝑤].

4. Optimal Ascending Trajectory

This section is focused on the numerical solution to the trajectory optimization problem defined in the previous sections and regarding the upper stage of the Muralm launch vehicle. To this end, the first-order conditions for optimality are used, in conjunction of a simple implementation of swarming algorithm.

4.1. Outline of the Particle Swarm Optimization Method

The particle swarm optimization (PSO) method is a heuristic technique aimed at finding the optimal values of a set of unknown parameters, for a generic dynamical system. In this research, the optimization problem involves a continuous time-dependent control variable and is translated into a parameter optimization problem through the necessary conditions for optimality, which allow expressing the control variable as a function of the adjoint variables conjugate to the dynamics equations. This means that unknown parameters for the problem at hand are 𝜆₁₀,𝜆₂₀,𝜆₃₀,𝑓₃,𝑡_𝑓. (53)

In general, unconstrained parameter optimization problems can be stated as follows: determine the optimal values of the 𝑛 unknown parameters {𝜒₁,…,𝜒_𝑛} such that the objective function 𝐽 is minimized. The PSO technique is a population-based method, where the population is represented by a swarm of 𝑁 particles. Each particle 𝑖(𝑖=1,…,𝑁) is associated with a position vector 𝝌(𝑖) and with a velocity vector 𝐰(𝑖). The position vector includes the values of the 𝑛 unknown parameters of the problem, 𝝌(𝑖)≜[𝜒₁(𝑖),…,𝜒_𝑛(𝑖)]^𝑇, whereas the velocity vector, whose components are denoted with 𝑤_𝑘(𝑖)(𝑘=1,…,𝑛), determines the position update. Each particle represents a possible solution to the problem and corresponds to a specific value of the objective function. The initial population is randomly generated by introducing 𝑁 particles, whose positions and velocities are (stochastically) uniformly distributed in the respective search spaces. The expressions for position and velocity update determine the swarm evolution toward the location of the globally optimal position, which corresponds to the globally optimal solution to the problem of interest. The algorithm terminates when the maximum number of iterations 𝑁_IT is reached or an alternative convergence criterion is met. The position vector of the best particle is expected to contain the optimal values of the unknown parameters, which correspond to the global minimum of 𝐽.

The central idea underlying the method is contained in the formula for velocity updating. This formula includes three terms with stochastic weights: the first term is the so-called inertial component and for each particle is proportional to its velocity in the preceding iteration; the second component is termed the cognitive component, directed toward the personal best position, that is, the best position experienced by the particle; and finally the third term is the social component, directed toward the global best position, that is, the best position yet located by any particle in the swarm. Further details are reported in [16, 17].

Constrained optimization problems involve equalities and/or inequalities, regarding (directly or indirectly) the unknown parameters. Equality constraints narrow considerably the search space where feasible solutions can be located. This is due to the fact that (nonredundant) equality constraints actually reduce the degree of freedom of the problem according to their number. In fact, 𝑚 equality constraints reduce the degree of freedom by 𝑚. Therefore, in the presence of 𝑛 unknown parameters, at most 𝑚=𝑛 equality constraints are admissible (𝑚≤𝑛): 𝑙_𝑟(𝝌)=0(𝑟=1,…,𝑚). The most popular approach for dealing with these constraints consists in penalizing them by adjoining additional terms to the objective function: 𝐽̃=𝐽+^𝑚_𝑟=1𝜂_𝑟𝑙_𝑟(𝝌). (54) This approach is employed also in this research, with 𝜂_𝑟=1(𝑟=1,2,3,4) and equality constraints given by (9) and (21). As the necessary conditions for optimality are used to express the control variable, the PSO assumes the expression (54) (with 𝐽≡0) as the objective. Satisfaction of the first-order conditions is indicative of optimality of the solution.

4.2. Numerical Results

The ascending trajectories are determined by employing canonical units: the distance unit (DU) is the Earth radius (𝑅_𝐸), whereas the time unit (TU) is such that 𝜇_𝐸=1DU³/TU². Thus, 1DU=6378.136km and 1TU=806.8sec. The desired final radius is 𝑅_𝑑=𝑅_𝐸+280km; the true anomaly 𝑓₂ (at 𝑡₀) is specified.

The swarming optimizer is employed with the following settings: 𝑁=50 (number of particles) and 𝑁_IT=500 (maximum number of iterations). The algorithm terminates prematurely its execution if the objective 𝐽̃ (cf. (35)) decreases below the value 10⁻⁶. The optimal values of the unknown parameters are sought in the following ranges: −1≤𝜆_𝑘0≤1(𝑘=1,2,3), 𝑓₂≤𝑓₃≤𝑓₂+𝜋, and 0.01TU≤𝑡_𝑓≤0.2TU. It is worth noticing that constraint reduction allows arbitrarily defining the search space for the initial values of the Lagrange multipliers. This means that they can be sought in the interval −1≤𝜆_𝑘0≤1 by the PSO algorithm.

Only 126 PSO iterations were needed to generate the optimal solution to the specified accuracy, corresponding to the following in-plane orbit elements and mass at burnout: 𝑎_𝑓=𝑅_𝐸𝑒+280.004km_𝑓=7.260⋅10⁻⁷𝑚_𝑓=140.4kg. (55) Figures 2, 3, and 4 portray the state components, that is, altitude (directly related to the radius 𝑟), inertial velocity, and flight path angle. Figure 5 illustrates the costate variables, whereas Figure 6 depicts the optimal control time history. It is apparent that the thrust direction is nearly aligned with the velocity for the entire time of powered flight of the upper stage.

Figure 2 

Optimal ascending trajectory: altitude time history.

Figure 3 

Optimal ascending trajectory: velocity time history.

Figure 4 

Optimal ascending trajectory: flight path angle time history.

Figure 5 

Optimal ascending trajectory: adjoint variables time history.

Figure 6 

Optimal ascending trajectory: optimal thrust direction time history.

Two second-order conditions are described in Section 3.2, that is, the Clebsch-Legendre condition on the Hessian 𝐻_𝑢𝑢 and the nonexistence of conjugate points (52). These conditions are checked for the optimal ascending path. First, due to (22), 𝐻_𝑢𝑢 assumes the form 𝐻_𝑢𝑢=𝑇⁽³⁾𝑚⁽³⁾𝜆₃𝑥₂²+𝜆²₂, (56) which is nonnegative. From inspecting Figure 7, it is apparent that 𝐻_𝑢𝑢 is positive for the entire time of flight. Moreover, backward integration of (46) yields the time history of 𝑆. Figure 8 reports the absolute value of the determinant of 𝑆 (in logarithmic scale) and proves that 𝑆 is finite at any time in the interval [𝑡₀,𝑡_𝑓[. This suffices to demonstrate the two previously mentioned second-order conditions, which enforce optimality to second order.

Figure 7 

Optimal ascending trajectory: Hessian 𝐻𝑢𝑢.

Figure 8 

Absolute value of the determinant of 𝑆.

5. Concluding Remarks

Trajectory optimization of multistage launch vehicles is a challenging problem, treated with different approaches in the past. In this research, the final objective of minimizing fuel consumption is achieved by finding the optimal control time history and the optimal thrust and coast durations. Specifically, for the upper stage the existence and duration of a coast arc (with no propulsion) and the optimal thrust direction are determined through the first-order necessary conditions for optimality, that is, the Euler-Lagrange equations and the Pontryagin minimum principle, in conjunction with a heuristic method, that is, the particle swarm algorithm. This technique is very intuitive and easy to implement and leads to numerically determining the state, control, and adjoint variables (conjugate to the dynamics equations). A further advantage is represented by the fact that the swarming approach does not need any starting guess to yield an accurate solution. In addition, the availability of the adjoint variables allows the direct evaluation of the second-order conditions for optimality, that is, the Clebsch-Legendre condition and the nonexistence of conjugate points. The numerical process for verifying these conditions is based on the sweep method and is illustrated in full detail in this work. The satisfaction of the second-order conditions enforces optimality of the ascending path found in this study and represents an essential premise for the possible development of an efficient neighboring optimal guidance.