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

Authors: Mauro Pontani; Giampaolo Cecchetti
Journal: ISRN Operations Research
 (2013)

---

## 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.


---

## Body

## 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. Section2 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 radiusRd 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. Letm0(i) denote the initial mass of subrocket i. This mass m0(i) is composed of a structural mass mS(i), a propellant mass mP(i), and a payload mass mU(i) (m0(i)=mS(i)+mP(i)+mU(i)). For the first two subrockets, mU(i) (i=1,2) coincides with the initial mass of the subsequent subrocket (i.e., mU(i)≡m0(i+1)). With regard to the third subrocket, mU(3) is the (actual) payload mass, which is to be maximized. The mass distribution for the Muralm three-stage rocket is the following:
(1)m0(1)=3608kgmS(1)=306kgmP(1)=2480kgmU(1)=822kgm0(2)=822kgmS(2)=46kgmP(2)=370kgmU(2)=406kg≡m0(3).
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 magnitudeT and specific impulse ISP. Both these quantities are assumed as time-independent for all of the stages. For the upper stage, the following propulsive data are assumed:
(2)T(3)=7500N,ISP(3)=320sec,
where the superscript refers to stage 3. For each subrocket i, it is then straightforward to derive the thrust acceleration as a function of the time t:
(3)T(i)m(i)=n0(i)1-n0(i)(t-tign(i))/(g0ISP(i)),
where n0(i) is the thrust acceleration at ignition of stage i (which occurs at the time tign(i)) and g0 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(c^1,c^2,c^3); c^1 corresponds to the vernal axis, whereas c^3 is aligned with the Earth rotation axis and directed toward North (cf. Figure 1). The orbital plane contains the unit vectors N^, associated with the ascending node crossing (corresponding to the point at which the space vehicle crosses the equatorial plane from South to North), e^, aligned with the direction of perigee, and r^, 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 i(0≤i≤π). The argument of perigee ω defines the position of perigee on the orbit plane, whereas the instantaneous true anomaly f(t) 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 i and RAAN Ω do not vary during the last stage propulsion. In contrast, the remaining orbit elements, that is, the orbit semimajor axis a, eccentricity e, and argument of perigee ω, change up to reaching the desired final values (denoted with the subscript “d”); that is, ad=Rd and ed=0 (ω is not defined for circular orbits).At separation of the second stage, the previously mentioned orbit elements are known and denoted witha2, e2, and ω2. 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 Δf, suffices to describe the rocket dynamics. In fact, if t2 and t3 represent, respectively, the time of the second stage separation and the ignition time of the upper stage, then f3≜f(t3)=f(t2)+Δf=f2+Δf. The orbital elements a2, e2, and ω2 do not vary along the coast (Keplerian) arc. As a result, the radius r, inertial velocity vI, and inertial flight path angle γI at t3 are simply given by
(4)r(t3)-a2(1-e22)1+e2cosf3=0,(5)γI(t3)-arctan[e2sinf31+e2cosf3]=0,(6)vI(t3)-μEa2(1-e22)1+e22+2e2cosf3=0,
where μE(=398600.4km3/sec2) 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 r, vI, and γI:
(7)r˙=vIsinγI,v˙I=T(3)m(3)cosαIvI-μEr2sinγI,γ˙I=T(3)m(3)sinαIvI+(vIr-μEr2vI)cosγI,
where the subscript “3” will be omitted henceforth and can be written in compact form as state equations:
(8)x˙=f(x,u,a,t).
The thrust direction is specified through the thrust angle αI(≡u), which is the only time-varying control variable and is taken clockwise from the direction of the inertial velocity vector vI. The operational orbit is assumed to be a circular orbit of radius Rd. This means that the final conditions at orbit injection (which occurs at tf) are
(9)r(tf)-Rd=0vI(tf)-μERd=0γI(tf)=0.
The left-hand sides of (4) through (6) and (9) are collected in the vector ψ(x0,xf,a,t0,tf)(=0)(t0≡t3).


## 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. Letm0(i) denote the initial mass of subrocket i. This mass m0(i) is composed of a structural mass mS(i), a propellant mass mP(i), and a payload mass mU(i) (m0(i)=mS(i)+mP(i)+mU(i)). For the first two subrockets, mU(i) (i=1,2) coincides with the initial mass of the subsequent subrocket (i.e., mU(i)≡m0(i+1)). With regard to the third subrocket, mU(3) is the (actual) payload mass, which is to be maximized. The mass distribution for the Muralm three-stage rocket is the following:
(1)m0(1)=3608kgmS(1)=306kgmP(1)=2480kgmU(1)=822kgm0(2)=822kgmS(2)=46kgmP(2)=370kgmU(2)=406kg≡m0(3).
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 magnitudeT and specific impulse ISP. Both these quantities are assumed as time-independent for all of the stages. For the upper stage, the following propulsive data are assumed:
(2)T(3)=7500N,ISP(3)=320sec,
where the superscript refers to stage 3. For each subrocket i, it is then straightforward to derive the thrust acceleration as a function of the time t:
(3)T(i)m(i)=n0(i)1-n0(i)(t-tign(i))/(g0ISP(i)),
where n0(i) is the thrust acceleration at ignition of stage i (which occurs at the time tign(i)) and g0 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(c^1,c^2,c^3); c^1 corresponds to the vernal axis, whereas c^3 is aligned with the Earth rotation axis and directed toward North (cf. Figure 1). The orbital plane contains the unit vectors N^, associated with the ascending node crossing (corresponding to the point at which the space vehicle crosses the equatorial plane from South to North), e^, aligned with the direction of perigee, and r^, 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 i(0≤i≤π). The argument of perigee ω defines the position of perigee on the orbit plane, whereas the instantaneous true anomaly f(t) 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 i and RAAN Ω do not vary during the last stage propulsion. In contrast, the remaining orbit elements, that is, the orbit semimajor axis a, eccentricity e, and argument of perigee ω, change up to reaching the desired final values (denoted with the subscript “d”); that is, ad=Rd and ed=0 (ω is not defined for circular orbits).At separation of the second stage, the previously mentioned orbit elements are known and denoted witha2, e2, and ω2. 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 Δf, suffices to describe the rocket dynamics. In fact, if t2 and t3 represent, respectively, the time of the second stage separation and the ignition time of the upper stage, then f3≜f(t3)=f(t2)+Δf=f2+Δf. The orbital elements a2, e2, and ω2 do not vary along the coast (Keplerian) arc. As a result, the radius r, inertial velocity vI, and inertial flight path angle γI at t3 are simply given by
(4)r(t3)-a2(1-e22)1+e2cosf3=0,(5)γI(t3)-arctan[e2sinf31+e2cosf3]=0,(6)vI(t3)-μEa2(1-e22)1+e22+2e2cosf3=0,
where μE(=398600.4km3/sec2) 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 r, vI, and γI:
(7)r˙=vIsinγI,v˙I=T(3)m(3)cosαIvI-μEr2sinγI,γ˙I=T(3)m(3)sinαIvI+(vIr-μEr2vI)cosγI,
where the subscript “3” will be omitted henceforth and can be written in compact form as state equations:
(8)x˙=f(x,u,a,t).
The thrust direction is specified through the thrust angle αI(≡u), which is the only time-varying control variable and is taken clockwise from the direction of the inertial velocity vector vI. The operational orbit is assumed to be a circular orbit of radius Rd. This means that the final conditions at orbit injection (which occurs at tf) are
(9)r(tf)-Rd=0vI(tf)-μERd=0γI(tf)=0.
The left-hand sides of (4) through (6) and (9) are collected in the vector ψ(x0,xf,a,t0,tf)(=0)(t0≡t3).


## 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)J=tf-t3.
Hence, the optimization problem for the upper stage is the following: determine the optimal control law αI(t) and the optimal true anomaly f3 such that J is minimized. The ignition time t3 is then calculated through the Kepler’s law,
(11)t3=t2+a23μE×{E(t3)-E(t2)-e2sin[E(t3)]+e2sin[E(t2)]},
where E(tk) is the eccentric anomaly associated with f(tk)(k=2,3).Lettingx≜[x1x2x3]T=[rvIγI]T, to obtain the necessary conditions for an optimal solution a Hamiltonian H and a function of the boundary conditions Φ are introduced as
(12)H≜λTf=λ1x2sinx3+λ2[T(3)m(3)cosαIx2-μEx12sinx3]+λ3[T(3)m(3)sinαIx2+(x2x1-μEx12x2)cosx3],Φ≜J+υTψ=(tf-t3)+υ1[x10-a2(1-e22)1+e2cosf3]+υ2[x20-μEa2(1-e22)1+e22+2e2cosf3]+υ3[x30-arctane2sinf31+e2cosf3]+υ4[x1f-Rd]+υ5[x2f-μERd]+υ6x3f,
where xk0=xk(t3) and xkf=xk(tf)(k=1,2,3); λ(≜[λ1λ2λ3]T) and υ(≜[υ1υ2υ3υ4υ5υ6]T) 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
(13)J′=Φ(x0,xf,f3,t0,tf)+∫t0tf[H(x,u,t)-λTx˙]dt.


### 3.1. First-Order Conditions

The first-order necessary conditions for (local) optimality [14] include the adjoint equations:
(14)λ˙=-[∂H∂x]T,
which yield the following scalar relations for the costate λ(≜[λ1λ2λ3]T):
(15)λ˙1=(x2λ3cosx3)1x12-(2μEλ2sinx3+2μEλ3cosx3x2)1x13,(16)λ˙2=-λ1sinx3-λ3[cosx3(1x1+μEx12x22)-T(3)m(3)1x22sinαI],(17)λ˙3=-x2λ1cosx3+μEλ2cosx3x12+λ3sinx3(x2x1-μEx12x2),
in conjunction with the respective boundary conditions:
(18)λ0=-[∂Φ∂x0]T⟹λk0=-υk(k=1,2,3)(19)λf=[∂Φ∂xf]T⟹λkf=υk+3(k=1,2,3),
where the subscripts “0” and “f ” refer to t3 and tf, respectively. In the presence of initial conditions depending on a parameter (a≡f3) a pair of additional necessary conditions must hold [15]:
(20)Φa=0⟹∂Φ∂f3=0,Φaa≥0⟹∂2Φ∂f32≥0.
The first equation yields a relation that expresses λ30 as a function of λ10, λ20, and f3,
(21)λ30=λ20sinf31+e22+2e2cosf3e2+cosf3-λ10a2(1-e22)e2+cosf31+e22+2e2cosf3(1+e2cosf3)2.
In addition, the optimal control αI* can be expressed as a function of the costate through the Pontryagin minimum principle:
(22)αI*=argminαIH⟹sinαI*=-λ3x2[(λ3x2)2+λ22]-1/2,cosαI*=-λ2[(λ3x2)2+λ22]-1/2,
which implies the fulfillment of the stationarity condition of H with respect to u; that is,
(23)HuT=0.
Lastly, as the final time is unspecified, the following transversality condition must hold:
(24)H(tf)+∂Φ∂tf=0⟹n0(3)1-n0(3)(tf-t3)/(g0ISP(3))×[λ3fx2f]2+[λ2f]2-1=0.
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 λ, f3, and tf.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 λ1(0)=kλλ1*(0), λ2(0)=kλλ2*(0), and λ3(0)=kλλ3*(0)(kλ>0) (where the superscript “*” denotes the actual optimal value of a variable), then the same proportionality holds between λ and λ* at any t, due to homogeneity of (15)–(17). Second, if the control αI is written as a function of λ(=kλλ*) through (22), then one can recognize that αI coincides with αI*:
(25)sinαI=-λ3x2[(λ3x2)2+λ22]-1/2=-kλλ3*x2[(kλλ3*x2)2+(kλλ2*)2]-1/2≡sinαI*,cosαI=-λ2[(λ3x2)2+λ22]-1/2=-kλλ2*[(kλλ3*x2)2+(kλλ2*)2]-1/2≡cosαI*.
This circumstance implies that if λ is proportional to λ*, then the final conditions are fulfilled at the minimum final time tf* even when the control is expressed as a function of the nonoptimal λ. In contrast, using λ(=kλλ*) the transversality condition is violated, because the value of H(tf*), due to (24), turns out to be
(26)H(tf*)=-n0(3)1-n0(3)(tf*-t3*)/(g0ISP(3))×[λ3fx2f]2+[λ2f]2=-kλ≠-1.
Hence, if the proportionality condition λ=kλλ* holds, then optimal control αI* 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 statex*, costate λ*, and control u*. 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 controlu such that
(27)u-u*=0ift0≤t≤tpu-u*=const.iftp≤t≤tqu-u*suchthatΨxfδxf=0iftq≤t≤tf,
where t0<tp<tq<tf and Δt(≜tq-tp) is small. Under the assumption that Ψ depends only on xf, the condition Ψxfδxf=0 ensures fulfillment of the boundary conditions to first order (δxf 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 (u-u*) belongs to the class of strong variations. In contrast, the resulting state variation δx(=x-x*) is reasonably small, and the comparison state x is sufficiently close to the optimal state x*. Let J*′ denote the local minimum value of the objective (associated with u*) and J′ the value corresponding to u. After expanding ΔJ′(=J′-J*′) in Taylor series, one obtains (cf. [15])
(28)ΔJ′=[H-H*]pΔt,
where H and H* are associated with u and u*. The quantity ΔJ′ must be positive for any choice of u. Since tp is arbitrary and Δt>0, the condition
(29)H-H*>0
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 u=u*+δu (with δu sufficiently small), (29) can be expanded in Taylor series about the optimal solution:
(30)H-H*≃(Hu)*δu+12δuT(Huu)*δu.
Along the optimal path (Hu)*=0, this property and the fact that δu is arbitrary imply that the Hessian matrix (Huu)* must be positive definite; that is,
(31)(Huu)*>0.
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(δx,δu,da,dtf) satisfy the linear equations (deriving from (8), (14), and (23)):
(32)δx˙=fxδx+fuδu+fada(33)δλ˙=-Hxxδx-Hxuδu-Hxλδλ-Hxada(34)Huxδx+Huuδu+Huada+Huλδλ=0
accompanied by the boundary conditions (deriving from the boundary conditions (4)–(6), (9), (18), (19), and (24)):
(35)ψxfδxf+ψx0δx0+ψada+ψ′dtf=0(36)δλ0=-Φx0x0δx0-Φx0ada-ψx0Tdυ(37)δλf=Φxfxfδxf+Φxfada+ψxfTdυ+Φxf′Tdtf(38)Λxfδxf+Λx0δx0+Λada+ψ′Tdυ+Λ′dtf=0(39)Φax0δx0+Φaxfδxf+Φaada+Φatfdtf+ψadυ=0,
where Λ≜Hf+Φtf. After solving (34) for u,
(40)δu=-Huu-1(Huxδx+Huada+Huλδλ),
and inserting the resulting expression in (32)-(33), one obtains
(41)δx˙=Aδx-Bδλ+Dda(42)δλ˙=-Cδx-ATδλ-Eda,
where the matrices A, B, C, D, and E 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)δλ=Sδx+Rdυ+mda(44)0=RTδx+Qdυ+nda(45)0=mTδx+nTdυ+αda.
After inserting (43) through (45) into (32)-(33), one obtains
(46)S˙=-C-ATS-SA+SBSR˙=(SB-AT)Rm˙=(SB-AT)m-SD-EQ˙=RTBRn˙=RT(Bm-D)α˙=mT(Bm-D)-DTm-F
with final conditions that must be consistent with (35) through (39); that is,
(47)Sf=ΦxfxfQf=0Rf=ψxfTnf=ψamf=Φxfaαf=Φaa.
Equations (44) and (45) can be solved simultaneously at t0 to yield
(48)[dυda]=-V0-1U0Tδx0,
where
(49)U=[Rm],(50)V=[QnnTα].
If (48) is used at t0, then
(51)δλ0=(S0-U0V0-1U0T)δx0.
Letting S^≜S-UV-1UT, it is relatively straightforward to demonstrate that the first of (46) holds also for S^ (with S^ in place of S), with boundary condition S^→∞ as t→tf. 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 δλ0→0 as δx0→0, unless S^ tends to infinity at an internal time t-(t0≤t-<tf), which is referred to as conjugate point. If δλ0→0 and δx0→0, then also δu→0, due to (40)–(42) and (48). This means that no neighboring optimal path exists, unless S^→∞ at a certain time t-(t0≤t-<tf). In other words, the nonexistence of conjugate points is
(52)S^<∞∀t∈[t0,tf[.
It is worth remarking that the numerical backward integration for S^ can start by integrating (46) up to a time tsw sufficiently close to tf; then one can switch to S^ for the remaining time interval [t0,tsw].


## 3.1. First-Order Conditions

The first-order necessary conditions for (local) optimality [14] include the adjoint equations:
(14)λ˙=-[∂H∂x]T,
which yield the following scalar relations for the costate λ(≜[λ1λ2λ3]T):
(15)λ˙1=(x2λ3cosx3)1x12-(2μEλ2sinx3+2μEλ3cosx3x2)1x13,(16)λ˙2=-λ1sinx3-λ3[cosx3(1x1+μEx12x22)-T(3)m(3)1x22sinαI],(17)λ˙3=-x2λ1cosx3+μEλ2cosx3x12+λ3sinx3(x2x1-μEx12x2),
in conjunction with the respective boundary conditions:
(18)λ0=-[∂Φ∂x0]T⟹λk0=-υk(k=1,2,3)(19)λf=[∂Φ∂xf]T⟹λkf=υk+3(k=1,2,3),
where the subscripts “0” and “f ” refer to t3 and tf, respectively. In the presence of initial conditions depending on a parameter (a≡f3) a pair of additional necessary conditions must hold [15]:
(20)Φa=0⟹∂Φ∂f3=0,Φaa≥0⟹∂2Φ∂f32≥0.
The first equation yields a relation that expresses λ30 as a function of λ10, λ20, and f3,
(21)λ30=λ20sinf31+e22+2e2cosf3e2+cosf3-λ10a2(1-e22)e2+cosf31+e22+2e2cosf3(1+e2cosf3)2.
In addition, the optimal control αI* can be expressed as a function of the costate through the Pontryagin minimum principle:
(22)αI*=argminαIH⟹sinαI*=-λ3x2[(λ3x2)2+λ22]-1/2,cosαI*=-λ2[(λ3x2)2+λ22]-1/2,
which implies the fulfillment of the stationarity condition of H with respect to u; that is,
(23)HuT=0.
Lastly, as the final time is unspecified, the following transversality condition must hold:
(24)H(tf)+∂Φ∂tf=0⟹n0(3)1-n0(3)(tf-t3)/(g0ISP(3))×[λ3fx2f]2+[λ2f]2-1=0.
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 λ, f3, and tf.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 λ1(0)=kλλ1*(0), λ2(0)=kλλ2*(0), and λ3(0)=kλλ3*(0)(kλ>0) (where the superscript “*” denotes the actual optimal value of a variable), then the same proportionality holds between λ and λ* at any t, due to homogeneity of (15)–(17). Second, if the control αI is written as a function of λ(=kλλ*) through (22), then one can recognize that αI coincides with αI*:
(25)sinαI=-λ3x2[(λ3x2)2+λ22]-1/2=-kλλ3*x2[(kλλ3*x2)2+(kλλ2*)2]-1/2≡sinαI*,cosαI=-λ2[(λ3x2)2+λ22]-1/2=-kλλ2*[(kλλ3*x2)2+(kλλ2*)2]-1/2≡cosαI*.
This circumstance implies that if λ is proportional to λ*, then the final conditions are fulfilled at the minimum final time tf* even when the control is expressed as a function of the nonoptimal λ. In contrast, using λ(=kλλ*) the transversality condition is violated, because the value of H(tf*), due to (24), turns out to be
(26)H(tf*)=-n0(3)1-n0(3)(tf*-t3*)/(g0ISP(3))×[λ3fx2f]2+[λ2f]2=-kλ≠-1.
Hence, if the proportionality condition λ=kλλ* holds, then optimal control αI* 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 statex*, costate λ*, and control u*. 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 controlu such that
(27)u-u*=0ift0≤t≤tpu-u*=const.iftp≤t≤tqu-u*suchthatΨxfδxf=0iftq≤t≤tf,
where t0<tp<tq<tf and Δt(≜tq-tp) is small. Under the assumption that Ψ depends only on xf, the condition Ψxfδxf=0 ensures fulfillment of the boundary conditions to first order (δxf 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 (u-u*) belongs to the class of strong variations. In contrast, the resulting state variation δx(=x-x*) is reasonably small, and the comparison state x is sufficiently close to the optimal state x*. Let J*′ denote the local minimum value of the objective (associated with u*) and J′ the value corresponding to u. After expanding ΔJ′(=J′-J*′) in Taylor series, one obtains (cf. [15])
(28)ΔJ′=[H-H*]pΔt,
where H and H* are associated with u and u*. The quantity ΔJ′ must be positive for any choice of u. Since tp is arbitrary and Δt>0, the condition
(29)H-H*>0
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 u=u*+δu (with δu sufficiently small), (29) can be expanded in Taylor series about the optimal solution:
(30)H-H*≃(Hu)*δu+12δuT(Huu)*δu.
Along the optimal path (Hu)*=0, this property and the fact that δu is arbitrary imply that the Hessian matrix (Huu)* must be positive definite; that is,
(31)(Huu)*>0.
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(δx,δu,da,dtf) satisfy the linear equations (deriving from (8), (14), and (23)):
(32)δx˙=fxδx+fuδu+fada(33)δλ˙=-Hxxδx-Hxuδu-Hxλδλ-Hxada(34)Huxδx+Huuδu+Huada+Huλδλ=0
accompanied by the boundary conditions (deriving from the boundary conditions (4)–(6), (9), (18), (19), and (24)):
(35)ψxfδxf+ψx0δx0+ψada+ψ′dtf=0(36)δλ0=-Φx0x0δx0-Φx0ada-ψx0Tdυ(37)δλf=Φxfxfδxf+Φxfada+ψxfTdυ+Φxf′Tdtf(38)Λxfδxf+Λx0δx0+Λada+ψ′Tdυ+Λ′dtf=0(39)Φax0δx0+Φaxfδxf+Φaada+Φatfdtf+ψadυ=0,
where Λ≜Hf+Φtf. After solving (34) for u,
(40)δu=-Huu-1(Huxδx+Huada+Huλδλ),
and inserting the resulting expression in (32)-(33), one obtains
(41)δx˙=Aδx-Bδλ+Dda(42)δλ˙=-Cδx-ATδλ-Eda,
where the matrices A, B, C, D, and E 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)δλ=Sδx+Rdυ+mda(44)0=RTδx+Qdυ+nda(45)0=mTδx+nTdυ+αda.
After inserting (43) through (45) into (32)-(33), one obtains
(46)S˙=-C-ATS-SA+SBSR˙=(SB-AT)Rm˙=(SB-AT)m-SD-EQ˙=RTBRn˙=RT(Bm-D)α˙=mT(Bm-D)-DTm-F
with final conditions that must be consistent with (35) through (39); that is,
(47)Sf=ΦxfxfQf=0Rf=ψxfTnf=ψamf=Φxfaαf=Φaa.
Equations (44) and (45) can be solved simultaneously at t0 to yield
(48)[dυda]=-V0-1U0Tδx0,
where
(49)U=[Rm],(50)V=[QnnTα].
If (48) is used at t0, then
(51)δλ0=(S0-U0V0-1U0T)δx0.
Letting S^≜S-UV-1UT, it is relatively straightforward to demonstrate that the first of (46) holds also for S^ (with S^ in place of S), with boundary condition S^→∞ as t→tf. 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 δλ0→0 as δx0→0, unless S^ tends to infinity at an internal time t-(t0≤t-<tf), which is referred to as conjugate point. If δλ0→0 and δx0→0, then also δu→0, due to (40)–(42) and (48). This means that no neighboring optimal path exists, unless S^→∞ at a certain time t-(t0≤t-<tf). In other words, the nonexistence of conjugate points is
(52)S^<∞∀t∈[t0,tf[.
It is worth remarking that the numerical backward integration for S^ can start by integrating (46) up to a time tsw sufficiently close to tf; then one can switch to S^ for the remaining time interval [t0,tsw].


## 3.2.1. Weierstrass and Legendre-Clebsch Conditions

An admissible comparison path can be defined by introducing a comparison controlu such that
(27)u-u*=0ift0≤t≤tpu-u*=const.iftp≤t≤tqu-u*suchthatΨxfδxf=0iftq≤t≤tf,
where t0<tp<tq<tf and Δt(≜tq-tp) is small. Under the assumption that Ψ depends only on xf, the condition Ψxfδxf=0 ensures fulfillment of the boundary conditions to first order (δxf 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 (u-u*) belongs to the class of strong variations. In contrast, the resulting state variation δx(=x-x*) is reasonably small, and the comparison state x is sufficiently close to the optimal state x*. Let J*′ denote the local minimum value of the objective (associated with u*) and J′ the value corresponding to u. After expanding ΔJ′(=J′-J*′) in Taylor series, one obtains (cf. [15])
(28)ΔJ′=[H-H*]pΔt,
where H and H* are associated with u and u*. The quantity ΔJ′ must be positive for any choice of u. Since tp is arbitrary and Δt>0, the condition
(29)H-H*>0
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 u=u*+δu (with δu sufficiently small), (29) can be expanded in Taylor series about the optimal solution:
(30)H-H*≃(Hu)*δu+12δuT(Huu)*δu.
Along the optimal path (Hu)*=0, this property and the fact that δu is arbitrary imply that the Hessian matrix (Huu)* must be positive definite; that is,
(31)(Huu)*>0.
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(δx,δu,da,dtf) satisfy the linear equations (deriving from (8), (14), and (23)):
(32)δx˙=fxδx+fuδu+fada(33)δλ˙=-Hxxδx-Hxuδu-Hxλδλ-Hxada(34)Huxδx+Huuδu+Huada+Huλδλ=0
accompanied by the boundary conditions (deriving from the boundary conditions (4)–(6), (9), (18), (19), and (24)):
(35)ψxfδxf+ψx0δx0+ψada+ψ′dtf=0(36)δλ0=-Φx0x0δx0-Φx0ada-ψx0Tdυ(37)δλf=Φxfxfδxf+Φxfada+ψxfTdυ+Φxf′Tdtf(38)Λxfδxf+Λx0δx0+Λada+ψ′Tdυ+Λ′dtf=0(39)Φax0δx0+Φaxfδxf+Φaada+Φatfdtf+ψadυ=0,
where Λ≜Hf+Φtf. After solving (34) for u,
(40)δu=-Huu-1(Huxδx+Huada+Huλδλ),
and inserting the resulting expression in (32)-(33), one obtains
(41)δx˙=Aδx-Bδλ+Dda(42)δλ˙=-Cδx-ATδλ-Eda,
where the matrices A, B, C, D, and E 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)δλ=Sδx+Rdυ+mda(44)0=RTδx+Qdυ+nda(45)0=mTδx+nTdυ+αda.
After inserting (43) through (45) into (32)-(33), one obtains
(46)S˙=-C-ATS-SA+SBSR˙=(SB-AT)Rm˙=(SB-AT)m-SD-EQ˙=RTBRn˙=RT(Bm-D)α˙=mT(Bm-D)-DTm-F
with final conditions that must be consistent with (35) through (39); that is,
(47)Sf=ΦxfxfQf=0Rf=ψxfTnf=ψamf=Φxfaαf=Φaa.
Equations (44) and (45) can be solved simultaneously at t0 to yield
(48)[dυda]=-V0-1U0Tδx0,
where
(49)U=[Rm],(50)V=[QnnTα].
If (48) is used at t0, then
(51)δλ0=(S0-U0V0-1U0T)δx0.
Letting S^≜S-UV-1UT, it is relatively straightforward to demonstrate that the first of (46) holds also for S^ (with S^ in place of S), with boundary condition S^→∞ as t→tf. 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 δλ0→0 as δx0→0, unless S^ tends to infinity at an internal time t-(t0≤t-<tf), which is referred to as conjugate point. If δλ0→0 and δx0→0, then also δu→0, due to (40)–(42) and (48). This means that no neighboring optimal path exists, unless S^→∞ at a certain time t-(t0≤t-<tf). In other words, the nonexistence of conjugate points is
(52)S^<∞∀t∈[t0,tf[.
It is worth remarking that the numerical backward integration for S^ can start by integrating (46) up to a time tsw sufficiently close to tf; then one can switch to S^ for the remaining time interval [t0,tsw].


## 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){λ10,λ20,λ30,f3,tf}.In general, unconstrained parameter optimization problems can be stated as follows: determine the optimal values of then unknown parameters {χ1,…,χn} such that the objective function J is minimized. The PSO technique is a population-based method, where the population is represented by a swarm of N particles. Each particle i(i=1,…,N) is associated with a position vector χ(i) and with a velocity vector w(i). The position vector includes the values of the n unknown parameters of the problem, χ(i)≜[χ1(i),…,χn(i)]T, whereas the velocity vector, whose components are denoted with wk(i)(k=1,…,n), 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 N 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 NIT 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 J.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,m equality constraints reduce the degree of freedom by m. Therefore, in the presence of n unknown parameters, at most m=n equality constraints are admissible (m≤n): lr(χ)=0(r=1,…,m). The most popular approach for dealing with these constraints consists in penalizing them by adjoining additional terms to the objective function:
(54)J~=J+∑r=1mηr|lr(χ)|.
This approach is employed also in this research, with ηr=1(r=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 J≡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(RE), whereas the time unit (TU) is such that μE=1DU3/TU2. Thus, 1DU=6378.136km and 1TU=806.8sec. The desired final radius is Rd=RE+280km; the true anomaly f2 (at t0) is specified.The swarming optimizer is employed with the following settings:N=50 (number of particles) and NIT=500 (maximum number of iterations). The algorithm terminates prematurely its execution if the objective J~ (cf. (35)) decreases below the value 10-6. The optimal values of the unknown parameters are sought in the following ranges: -1≤λk0≤1(k=1,2,3), f2≤f3≤f2+π, and 0.01TU≤tf≤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≤λk0≤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:(55)af=RE+280.004kmef=7.260·10-7mf=140.4kg.
Figures 2, 3, and 4 portray the state components, that is, altitude (directly related to the radius r), 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 Section3.2, that is, the Clebsch-Legendre condition on the Hessian Huu and the nonexistence of conjugate points (52). These conditions are checked for the optimal ascending path. First, due to (22), Huu assumes the form
(56)Huu=T(3)m(3)(λ3x2)2+λ22,
which is nonnegative. From inspecting Figure 7, it is apparent that Huu is positive for the entire time of flight. Moreover, backward integration of (46) yields the time history of S^. Figure 8 reports the absolute value of the determinant of S^ (in logarithmic scale) and proves that S^ is finite at any time in the interval [t0,tf[. This suffices to demonstrate the two previously mentioned second-order conditions, which enforce optimality to second order.Figure 7
Optimal ascending trajectory: HessianHuu.Figure 8
Absolute value of the determinant ofS^.


## 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){λ10,λ20,λ30,f3,tf}.In general, unconstrained parameter optimization problems can be stated as follows: determine the optimal values of then unknown parameters {χ1,…,χn} such that the objective function J is minimized. The PSO technique is a population-based method, where the population is represented by a swarm of N particles. Each particle i(i=1,…,N) is associated with a position vector χ(i) and with a velocity vector w(i). The position vector includes the values of the n unknown parameters of the problem, χ(i)≜[χ1(i),…,χn(i)]T, whereas the velocity vector, whose components are denoted with wk(i)(k=1,…,n), 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 N 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 NIT 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 J.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,m equality constraints reduce the degree of freedom by m. Therefore, in the presence of n unknown parameters, at most m=n equality constraints are admissible (m≤n): lr(χ)=0(r=1,…,m). The most popular approach for dealing with these constraints consists in penalizing them by adjoining additional terms to the objective function:
(54)J~=J+∑r=1mηr|lr(χ)|.
This approach is employed also in this research, with ηr=1(r=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 J≡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(RE), whereas the time unit (TU) is such that μE=1DU3/TU2. Thus, 1DU=6378.136km and 1TU=806.8sec. The desired final radius is Rd=RE+280km; the true anomaly f2 (at t0) is specified.The swarming optimizer is employed with the following settings:N=50 (number of particles) and NIT=500 (maximum number of iterations). The algorithm terminates prematurely its execution if the objective J~ (cf. (35)) decreases below the value 10-6. The optimal values of the unknown parameters are sought in the following ranges: -1≤λk0≤1(k=1,2,3), f2≤f3≤f2+π, and 0.01TU≤tf≤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≤λk0≤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:(55)af=RE+280.004kmef=7.260·10-7mf=140.4kg.
Figures 2, 3, and 4 portray the state components, that is, altitude (directly related to the radius r), 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 Section3.2, that is, the Clebsch-Legendre condition on the Hessian Huu and the nonexistence of conjugate points (52). These conditions are checked for the optimal ascending path. First, due to (22), Huu assumes the form
(56)Huu=T(3)m(3)(λ3x2)2+λ22,
which is nonnegative. From inspecting Figure 7, it is apparent that Huu is positive for the entire time of flight. Moreover, backward integration of (46) yields the time history of S^. Figure 8 reports the absolute value of the determinant of S^ (in logarithmic scale) and proves that S^ is finite at any time in the interval [t0,tf[. This suffices to demonstrate the two previously mentioned second-order conditions, which enforce optimality to second order.Figure 7
Optimal ascending trajectory: HessianHuu.Figure 8
Absolute value of the determinant ofS^.


## 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.



---

*Source: 498765-2013-12-10.xml*