Chinese Journal of Aeronautics, (2013),26(6): 1544-1553

JOURNAL OF

AERONAUTICS

Chinese Society of Aeronautics and Astronautics & Beihang University

Chinese Journal of Aeronautics

cja@buaa.edu.cn www.sciencedirect.com

Reentry trajectory optimization for hypersonic vehicle satisfying complex constraints

Jiang Zhao, Rui Zhou *

School of Automation Science and Electrical Engineering, Beihang University, Beijing 100191, China

Received 19 November 2012; revised 1 April 2013; accepted 20 July 2013 Available online 6 November 2013

KEYWORDS

Hypersonic vehicles; Reentry trajectory optimization;

Multi-phase Gauss pseudospectral method (MGPM); Waypoint; No-fly zone

Abstract The reentry trajectory optimization for hypersonic vehicle (HV) is a current problem of great interest. Some complex constraints, such as waypoints for reconnaissance and no-fly zones for threat avoidance, are inevitably involved in a global strike mission. Of the many direct methods, Gauss pseudospectral method (GPM) has been demonstrated as an effective tool to solve the trajectory optimization problem with typical constraints. However, a series of difficulties arises for complex constraints, such as the uncertainty of passage time for waypoints and the inaccuracy of approximate trajectory near no-fly zones. The research herein proposes a multi-phase technique based on the GPM to generate an optimal reentry trajectory for HV satisfying waypoint and no-fly zone constraints. Three kinds of specific breaks are introduced to divide the full trajectory into multiple phases. The continuity conditions are presented to ensure a smooth connection between each pair of phases. Numerical examples for reentry trajectory optimization in free-space flight and with complex constraints are used to demonstrate the proposed technique. Simulation results show the feasible application of multi-phase technique in reentry trajectory optimization with way-point and no-fly zone constraints.

© 2013 Production and hosting by Elsevier Ltd. on behalf of CSAA & BUAA.

1. Introduction

Global strike and global persistent attack1 have been focused

by many countries since the human started the space era.

The hypersonic vehicle (HV), a new type of aerospace vehicle,

is of great importance due to its remarkable ability of reentering

the atmosphere to attack the ground targets fast. Of the various HV technologies to enhance global reach capability, reentry trajectory optimization has always been considered a difficult problem,2 which involves numerous nonlinear motion equations and nonlinear constraints. Generally, two categories of numerical methods are introduced for such problems: indirect methods and direct methods. The indirect methods3 are based on the Pontryagin's minimum principle that results in a Hamiltonian boundary-value problem (HBVP). A high accuracy in the solution is the primary advantage of the indirect methods, while the HBVP is quite complicated to solve.4 The direct methods mainly convert the optimal control problem to the nonlinear programming problem (NLP).5 It is easier to use due to the larger radii of convergence without deriving the first-order necessary conditions.3

* Corresponding author. Tel.: +86 10 82339232. E-mail addresses: jzhao@asee.buaa.edu.cn (J. Zhao), zhr@buaa. edu.cn (R. Zhou).

Peer review under responsibility of Editorial Committee of CJA.

1000-9361 © 2013 Production and hosting by Elsevier Ltd. on behalf of CSAA & BUAA. http://dx.doi.Org/10.1016/j.cja.2013.10.009

Over the past decade, Gauss pseudospectral method (GPM) has proved to be effective to solve the typical constraints in trajectory optimization for HV, such as heating rate, dynamic pressure, aerodynamic load, control boundaries, and terminal conditions. Previous studies of optimization problem using GPM are summarized as follows. Reddien6 firstly proposed that the optimal control problem can be solved by GPM. Later, Benson7'8 expatiated on the integral GPM and differential GPM, explicitly formulating a mapping between the Kar-ush-Kuhn-Tucker (KKT) conditions and the discretized firstorder necessary conditions that can be used to obtain an accurate estimate of the costate. Huntington3 improved the GPM by a revised pseudospectral transcription for problems with path constraints and differential dynamic constraints, and he also presented a new method to compute the control at boundaries. Fahroo and Ross9,10 formulated a generalized costate mapping theorem that clarified the selection of correct pseudospectral method for solving problems efficiently. The practical applications of GPM to trajectory optimization for HV have already focused on ascent trajectory reconstruction11 and reentry trajectory optimization12-16 with typical constraints.

However, since HVs are being designed versatile, complex constraints are inevitably involved in trajectory optimization problems, such as the waypoints for reconnaissance and no-fly zones for threat avoidance or geopolitical restrictions. To solve these complex constraints using GPM, some difficulties arise as follows: (1) the passage time of waypoints is uncertain, which brings trouble to the choice of discrete nodes for the state and control; (2) the discrete nodes are sparse in the middle of general interval, so that the approximate trajectory is less accurate near the no-fly zones; (3) the increasing order of the interpolating polynomial probably causes large approximate error in the state and control.

The overall objective of the research is to generate an optimal reentry trajectory for HV satisfying waypoint and no-fly zone constraints. The multi-phase Gauss pseudospectral method (MGPM) is proposed to overcome the three difficulties above. Three kinds of specific points are introduced to divide the full trajectory. The main results on analysis of MGPM are presented to validate the feasibility of the multiphase technique.

2. Problem formulation

where r is the radial distance. h and / are the longitude and latitude. V is the Earth-relative velocity. y and W are the flight-path angle and heading angle. r is the bank angle. m is the mass of the vehicle. D and L are the aerodynamic drag force and lift force. X and g are the Earth angular velocity and gravitational acceleration. D and L are given as

D = - q V Cd Sref

1 2 L = 2 qV2CLSref

where Sref is the reference area. q = q0e—b is the atmospheric density with q0 being the atmospheric density at sea level, b the density coefficient constant, and h the altitude from the sea level. Describe the Earth radius as Re, then h is given as h = r — Re. The coefficients of drag and lift are expressed by angle of attack a and Mach number Ma as

Cl = clo + CL1a + CL2eCL3Ma CD = cdo + CD1a2 + CmeCD3'Ma

where the coefficients CLi and CDi (i = 0,1,2,3) are based on the common aero vehicle (CAV) data.18

2.2. Multiple constraints 2.2.1. Heating rate

Qd indicates the heating rate at a specified point on the surface of HV, generally the stagnation point on the nose. It is determined by the atmospheric density and Earth-relative velocity. To satisfy the serious thermal protection constraints, the heating rate must be in the limit as

,0.5 r/3.15

qd = kqp v3 qd - qd

where Kq is the heating rate normalization constant. Note that, all the subscripts "max/min" denote the maximum/minimum allowable values throughout the paper.

2.2.2. Dynamic pressure

The dynamic pressure, Pd, is a typical path constraint, limiting the hinge moment of actuator in the reasonable range. The peak value of dynamic pressure must be less than the maximum value as

This section describes the reentry dynamic model and multiple constraints for HV. The formulation of trajectory optimization problem is also presented.

2.1. Reentry dynamics

Using a spherical rotating Earth model, the point-mass dynamics of reentry vehicle is given as17

' r = Vsin y d =(Vcos y sin W)/(r cos/) / = ( Vcos y cos W)/r

V =—D/m — g sin y + X2r(sin y cos / — cos y sin / sin W) cos / y = (L cos r)/(mV) — (g cos y)/V +(Vcos y)/r

+2X cos / cos W + X2r cos/(cos y cos / + sin W sin / sin y)/V W = (L sin r)/(mVcos y) — (Vcos y cos W tan/)/r +2X(tan y cos / sin W — sin/) — X2r sin / cos / cos W/( Vcos y)

Pd — qV2 Pd

Pdmax 6 0

2.2.3. Aerodynamic load

To protect the mechanical system, the third "hard" constraint during the reentry flight is the aerodynamic load, nL. The instantaneous value and the maximum value are given as

( p l2+d2

I nL -~mor (6)

mg0 nLmax 6 0

where g0 is the gravity coefficient at Earth surface. 2.2.4. Control boundaries

To maintain the stability of vehicle, the two control of HV, the angle of attack and the bank angle, should satisfy the following constraints as

max 7 max;

a 6 ^max

(7) x(s) =

f — t0

f(x(s), u(s), s; to, tf)

2.2.5. Terminal conditions

Generally, the terminal conditions correspond to different flight missions. Let subscript "f" denote the terminal state, then, the terminal altitude, position, velocity, flight-path angle, and heading angle are given as

(x(so), to, x(sf), tf) = 0

hf = hf, Of = Of, /f = /f Vf = Vf, Cf = cî, Wf = Wî

2.2.6. Waypoints

To meet the requirement of reconnaissance mission, flight calibration, and payload delivery, HV must fly directly over a series of waypoints. Here, the position of the ith waypoint is presented by its longitude and latitude (Qh /). Thus, the way-point constraint is given as

p(o(ti),/m = m)- Oif+m)- /)2=o (i=1, 2,..., p)

where P is the total number of waypoints and (O(t), /(t)) the current position of vehicle.

C(x(s), u(s), s; to, tf) 6 0

Generally, the objective function is defined by the user corresponding to the specified mission. For HV, the objective function can be selected as the minimum total heat, maximum downrange, maximum crossrange, minimum arriving time, etc.

3. Methodology

In this section, we first introduce a numerical algorithm based on GPM for reentry trajectory optimization in free-space flight. Later, we propose a multi-phase technique to generate an optimal reentry trajectory satisfying waypoint and no-fly zone constraints. We refer to three kinds of phase breaks between each pair of segmented trajectories: waypoints, quasi-contact points, and turn points.

3.1. Traditional GPM

2.2.7. No-fly zones

In actual flight, additional geopolitical sensitive regions and threats should be imposed. The HV must not violate the boundary of these regions. In this paper, no-fly zones are specified as cylinder zones with infinite altitude. The cross section of the jth no-fly zone is described by the center (Oj, /j) and the radius RZj. The no-fly zone constraints are given as

S(O(tj),/(tj)) = (O(tj)-Oj)2 + (/(tj)— j p R^ (j = 1,2,...,S) (10)

where S is the total number of no-fly zones.

2.3. Reentry trajectory optimization

Subject to the dynamic model, the purpose of reentry trajectory optimization problem is to find the control, a and r, such that the objective function is a minimum (or a maximum), meanwhile satisfying the multiple constraints (Eqs. (4)-(10)).

The reentry trajectory optimization problem can be described as the general optimal control problem, and here, we consider it in the Bolza form.7 First, we map the state in dynamic equations (Eqs. (1)-(3)) to the general interval s 2 [s0, Tf] = [-1,1] by the following affine transformation:

tf + to

tf — to tf — to

where t0 and tf are the initial time and final time, respectively.

The continuous Bolza problem is to determine the control u(s) 2 Rm and the state x(s) 2 Rn that minimizes the Bolza objective function

tf — t0

J = U(x(so), to, x(sf), tf) -

X / g(x(s), u(s), s; to, tf) ds

subject to the dynamic constraints, boundary constraints, and path constraints, which are stated formally as7

In free-space flight, the typical constraints (Eqs. (4)-(8)) are considered. The basic principle of GPM is described as follows. The state and control of the dynamic equations are dis-cretized at Legendre-Gauss (LG) nodes, and then approximated by Lagrange interpolating polynomials. The derivatives of each state can be approximated by differentiating the global interpolating polynomials. The terminal states are presented by the initial states and a Gauss quadrature. The integral parts of the objective function are also approximated by the Gauss quadrature.

The description herein is based on the researches of Ben-son8 and Huntington.3 To solve the continuous Bolza problem, we first collocate the state and control in the dynamic equation (Eq. (13)) at N LG nodes sk (k =1,2,...,N) that are the roots of the Nth degree Legendre polynomial. With the two boundary nodes, s0 and sf, there are N +2 discretized nodes in total. Thus, the state, x(s), is formed with a basis of N +1 Lagrange interpolating polynomials L as

x(s) « X(s) = L (s)X(s. )

L,(s) = H

j=oj-isi sJ

(i = o, 1, ... , N)

Similarly, the control, u(s), is formed with a basis of N Lagrange interpolating polynomials L as

s(s) « U(s) ^Lf(s) U(s,)

Lf(s) = II

j=1j-isi sj

(i = 1, 2, ... , N)

Then, the derivatives of each state at the LG node are described in the form of a differential approximation matrix D as

x(Tk)~ X(sk) = J2L(sk)X(s,) = X(s,) (k = 1, 2,... , N)

j=0 j=0

where the elements of differentiation matrix, D 2 RNx(N+1), are determined by8

(1 + sk)i> N(sk) + PN(sk)

, i—k

D)d = ¿.(St) = 7 (Sk - s,)[(1 + s,)pn(s,) + pn(s,)] ki ' k (1 + sk)P-N(s,) + 2Pn(s,) ,2[(1 + s,)Pn(s,) + pn(s,)] ,

where PN(s) is the Nth degree Legendre polynomial.

Thus, the dynamic equation (Eq. (13)) at the collocation nodes is transcribed into algebraic constraints as

N t - t

YJDk,X(s,)- t^-t0f(X(sk), U(sk), Sk; to, tf)= 0 (k = 1, 2,..., N)

Since the state at the final time is ignored by the state approximation equation (Eq. (16)), an additional constraint with the final state is required as t - t N

X(sf)- X(so)- Xkf(X(sk), U(sk), Sk; to, tf) = 0 (21)

where xk is the Gauss weights.

Finally, the objective function (Eq. (12)) is approximated by a Gauss quadrature as J = U(X(so), to, X(sf), tf)

tf - t0 N

"5Z™kg(X(Tk), U(sk), sk; to, tf)

with the boundary conditions and path constraints in the form of discretization as

/(X(so), to, X(sf), tf)=0 (23)

C(X(sk), U( sk), sk; to, tf) 6 0 (k = 1, 2, ... , N) (24)

The solution to the continuous Bolza problem (Eqs. (11)—(15)) is determined by the solution to the NLP with the objective function (Eq. (22)), dynamic constraints (Eqs. (20)), additional constraints ((21)), boundary constraints (Eq. (23)), and path constraints (Eq. (24)). In addition, the solution must satisfy a set of KKT conditions that are described by Huntington3in detail.

3.2. MGPM

Considering the waypoint and no-fly zone constraints (Eqs. (9) and (1o)), the basic principle of MGPM is described as

follows. First, we divide the reentry trajectory into multiple phases by three kinds of specific points: waypoints, quasi-contact points, and turn points. Then, the state and control are discretized at LG nodes in each phase, separately. Similarly as the GPM, the dynamic equations are transcribed into algebraic constraints, and the terminal states are presented by the initial states and a Gauss quadrature. Finally, the continuity conditions on the time, state, and control are introduced to keep the connection smooth enough between each pair of phases. The objective function is determined by the combination of all the phases to ensure a global sub-optimal trajectory.

To demonstrate how to deal with the waypoint constraints using MGPM, we present an example of segmented trajectories and distribution of LG nodes. As shown in Fig. 1, fifteen LG nodes are selected using GPM that cannot be collocated at the waypoints as a result of the uncertainty of passage time for each waypoint. In contrast, the MGPM divides the full trajectory into four phases by two waypoints and a quasi-contact point. Three LG nodes are collocated in each phase. Thus, each waypoint turns to be the last LG node in the previous phase as well as the first LG node in the next phase. The waypoint constraints can be treated as the terminal conditions in these phases.

Using MGPM to eliminate the inaccuracy of the approximate trajectory near no-fly zones is described as follows. It is known that LG nodes in the general interval are sparse in the middle and dense on the two sides. This kind of distribution greatly satisfies the terminal conditions, while the approximate trajectory close to the no-fly zone is probably less accurate, as shown in Fig. 2. Using MGPM, the reentry trajectory is divided into phases by the quasi-contact points where the trajectory is likely to contact with no-fly zones. This leads to an increasing number of LG nodes around the quasi-contact point. In this way, the inaccuracy of trajectory near no-fly zones can be eliminated.

In general, the quasi-contact points are estimated by the pre-generation of reentry trajectory without divisions. As shown in Fig. 2, the boundary of no-fly zone is divided into two arcs by the flight trajectory penetrating through it. The trajectory is expected to go just along the boundary of the no-fly zone theoretically, such that the redundant flight distance is reduced. Herein, the quasi-contact point is supposed to be located at the boundary circle, obviously the small piece of arc. The selection of quasi-contact point is arbitrary in principle, since the dense distribution of LG nodes around the no-fly zone can ensure an accurate trajectory. The division by the quasi-contact point at the small arc can also lead to a global sub-optimal trajectory. However, the midpoint of the small arc is a better choice which benefits to generate a smoother trajectory and alleviate the demand for maneuverability. There-

Reentry point X, Waypoint 1 Quasi-contact point Waypoint 2 „ -----<• No-fly^\ zone ) \ Target

GPM (iV—15) 1* * * * * * * * * * * * * * * i

MGPM (iV=3) * * * t* * ** * * # ( * * * t

Fig. 1 Example of segmented trajectories and distribution of LG nodes.

GPM MGPM — Quasi-contact point

( No-fly \ ^^s. ( No-fly \

1 zone 1 I zone I

Fig. 2 Example of trajectories near the no-fly zone.

GPM ^ Max value MGPM Turn point /Max value

/ v^^—-

Fig. 3 Example of hard constraints near peak value.

fore, the midpoint of the small arc is selected as the quasi-contact point in the numerical examples.

The turn points, where the peak values of hard constraints stand, are also introduced to divide the full trajectory. During the reentry flight, an abrupt change in the control often leads to the peak value of heating rate, dynamic pressure, or aerodynamic load. However, LG nodes are not exactly distributed at the peak values using GPM, as shown in Fig. 3. To satisfy these hard constraints strictly, more LG nodes should be collocated around the turn points. It can be easily conducted by MGPM.

Unlike the quasi-contact points, the selection of turn points is fixed. The main goal of division by turn points is to ensure that the peak values of the hard constraints are less than the allowable maximum. Therefore, the peak value of each hard constraint is selected as the turn point. Note that, the turn point is not necessary when the peak value of some hard constraints is within the limits.

In addition, the MGPM can eliminate the Runge phenom-enon.3 Generally, the reentry trajectory has quite a long range, and the optimization problem requires a great number of the interpolating nodes to ensure the accuracy. However, the increasing order of interpolating polynomial probably results in large approximate errors in the state and control, which is known as the Runge phenomenon. Using MGPM, multiple phases are used instead, each of which needs less LG nodes than those using GPM, as shown in Fig. 1.

The detailed description of MGPM is presented as follows. Here, the total number of phases is described by M, and each phase has N LG nodes. Let P, S, and T denote the number of waypoints, no-fly zones, and turn points, respectively. Thus, M is determined by

M = P + S + T (25 )

Then, based on Eqs. (20)-(24), the objective function of the mth phase is given as

j(m) = U(m) (X<m) (s0m)) , t0m), (s(m)) , t(m))

t(m) _ t(m) N

+ tf 0 ^w(m)g(m) (Xm) (skm)) , Um) (skm)) , skm); t0m), t(m)) 2 k=1

(m = 1,2,...,M) (26)

Note that each phase has the same number of LG nodes in the objective function (Eq. (26)). If different numbers of LG nodes are collocated in the M phases, N<m) is used instead of N. The overall objective function, determined by the combination of all the phases, is described as

J = ^J<m) (m = 1, 2, ... , M) (27)

with the boundary conditions and path constraints as

/<m)(X<m)(?0m)),t0m),Xm)(s(m)),t(m))=0 (m = 1,2,...,M) (28)

C(m) (Xm) (skm)), Um) (skm)), skm); t0m), t(m)) 6 0 (k = 1, 2, ... , N; m = 1, 2,... , M) (29)

Similarly, the dynamic equation at the collocation points is transcribed into algebraic constraints in the form of multiple phases as

N t(m) — t(m)

YDm X(m) (s(m))—tf t0 f(m) (X<m> (skm)), Um) (skm)), skm); i0m), 4m) )= 0 i=0 2

(k = 1,2,...,N; m = 1,2,...,M) (30)

The additional constraint with final state in each phase is given as

Xm)(?[m)) — Xm)(?0m))

— tf 2 t0 Y1 akm)f <m)(X<m)(skm)), U<m)(skm)),skm); t0m), t(m)) 2 k=1

= 0 (m = 1, 2, ... , M) (31)

To ensure a smooth connection between each pair of phases, continuity conditions must be satisfied seriously. In detail, the time, state, and control at the last node in the previous phase must be the same with those at the first node in the next phase. The discretized forms of the continuity conditions are given as

t(m) — t0m+1)= 0 (m = 1, 2, ... , M — 1) (32)

X*m)(s<m)) — X<m+1)(?0m+1)) = 0 (m = 1, 2, ... , M — 1) (33)

U<m)(s(m)) — U*m+1)(?0m+1)) = 0 (m = 1, 2, ... , M — 1) (34)

Thus, the complete description of MGPM consists of the objective function (Eqs. (26) and (27)), boundary constraints (Eq. (28)), path constraints (Eq. (29)), dynamic constraints (Eq. (30)), additional constraints (Eq. (31)) and continuity conditions (Eqs. (32)-(34)).

3.3. Initial guess

Considering the reentry dynamics and trajectory constraints, two features of such optimization problems are as follows. First, the NLP converted from the continuous Bolza form includes large numbers of variables and constraints. Furthermore, the nonlinear

Table 1 Initial and terminal conditions of maximum crossrange trajectory for HV.

Condition h (km) V (m/s) h (°) / (°) c (°) W (°) a (°) r (°)

Initial 80 6500 0 0 -3 90 26 -45

Terminal 24 1000 - - -3 0 18 0

functions in the objective and constraints, as well as the solutions, are smooth. In general, the sequential quadratic programming (SQP) demonstrates an effective method for such a large-scale non-

linear optimization problem. Herein, the SQP algorithm that converges rapidly with great reliability is used to solve the resulting NLP.

(j) Dynamic pressure (k) Aerodynamic load

Fig. 4 Numerical results of reentry trajectory optimization using GPM.

Table 2 Initial and terminal conditions of minimum time trajectory for HV.

Conditions h (km) V (m/s) h (°) / (°) c (°) W (°) a (°) r (°)

Initial 80 6500 0 0 -3 90 26 -45

Terminal 24 1000 50 10 - - 18 0

A key step to enhance the performance of optimization problem is the choice of initial values for the NLP. An improper initial guess may lead the solutions of NLP to converge locally or even diverge. Traditional strategies for initial guess are summarized as follows.

(1) Running the NLP rapidly based on small number of discrete nodes to generate a reentry trajectory with lower accuracy. Then, interpolating the previous trajectory to obtain the initial values of the NLP based on desired number of discrete nodes.

(2) Integrating the equations of motion based on simplified reentry dynamics and estimation of control for initial values of the NLP. The method is less accurate than the first one.

(3) Estimating the initial values based on the exiting data and engineering experience. The method is generally hard to operate and lacks universality.

The first strategy is commonly adopted, and here, we use it for initial guess in the following examples.

4. Numerical examples

In this section, we present the simulation results of reentry trajectory optimization using GPM and MGPM. The aerodynamic data and characteristics parameters are based on the CAV-H19 data. The control boundaries and hard constraint limits remain fixed throughout all simulations: amax = 3oo, amin = 50, rmax = 890, rmin = -89o, Qdmax = 8.o X 1o5W/ m2, Pdmax = 5.o x 1o4 Pa, nLmax = 2.5. The types of the initial states and objective functions for the two algorithms are specified in each simulation separately. The optimization results are found through MATLAB 7.14.

4.1. GPM (N = 80)

In the simulation of a free-space flight, the initial and terminal conditions of maximum crossrange trajectory are described in Table 1. The objective function in the optimization is given as

J = max /f (35)

Eighty LG nodes are used to discretize the state and control. Fig. 4 shows the simulation results using GPM, including the state, control and path constraints for HV. The circle markers in Fig. 4 represent the approximate values at LG nodes. The total flight time is about 2191 s.

From the numerical results, it is observed that the 3D trajectory is smooth throughout the flight (Fig. 4(a)). In the ground trajectory by longitude and latitude, we can find a maximum crossrange of 39.63370 (about 4411.9 km) (Fig. 4(b)).The terminal conditions, such as the final altitude, velocity, flight-path angle, and heading angle, are satisfied accurately (Fig. 4(c)-(f)). Both the angle of attack and the bank angle are between the given control boundaries, respectively (Fig. 4(g)-(h)). The heating rate, dynamic pressure, and aerodynamic load are less than the maximum allowable values strictly (Fig. 4(i)-(k)). The solution demonstrates that GPM performs well to solve reentry trajectory optimization in free-space flight.

4.2. MGPM (M = 4, N = 10)

In the simulation of a reentry flight with two waypoints and a no-fly zone, the initial and terminal conditions of minimum time trajectory are described in Table 2. Table 3 presents the parameters of waypoint and no-fly zone constraints. The objective function in the optimization is given by

J = min K] t(m)

The full trajectory is divided into four phases, and ten nodes are collocated in each phase. Fig. 5 shows the simulation results using MGPM. The approximate values at LG nodes in the four phases are represented by different shape intervals. The total flight time is about 1623 s, and the flight in each phase lasts about 3o2 s, 377 s, 3o9 s, and 635 s, respectively.

From the 3D trajectory, we can clearly find that the four phases are smoothly connected, and the LG nodes around each waypoint and no-fly zone are denser than those in freespace flight (Fig. 5(a)). From the ground trajectory, it is observed that the two waypoints are directly passed through, and the trajectory near quasi-contact point goes just along the boundary of the no-fly zone (Fig. 5(b)). The numerical results also show that the terminal conditions are satisfied throughout the flight (Fig. 5(c)-(h)). Note that the typical hard constraints are all within the limits, so it is not necessary to divide the trajectory by turn points (Fig. 5(i)-(k)). The solution demonstrates that MGPM performs well to solve reentry trajectory optimization with complex constraints.

4.3. Discussion

Table 4 contains the comparisons on the effectiveness between GPM and MGPM for generating the reentry trajectory with

Table 3 Parameters of waypoint and no-fly zone constraints.

Constraint Parameter

Waypoint 1 Waypoint 2 No-fly zone Position: hj Position: h2 Center: h3 = = 15°, / = 3°. = 40°, /2 = 5° : 40°, /3 = 5°; Radius: RZ3 = 445.3 km

(k) Aerodynamic load

Fig. 5 Numerical results of reentry trajectory optimization using MGPM.

complex constraints. The differences of the two solutions are GPM. The no-fly zone constraints can be satisfied by both obvious in this problem. MGPM is capable of solving the way- MGPM and GPM. However, the generation of approximate point constraints with high accuracy, while not feasible by trajectory with higher accuracy by GPM requires large number

Table 4 Comparisons of the effectiveness between GPM and MGPM.

Method LG node Waypoint constraint No-fly zone constraint CPU time (s)

GPM N =80 Not feasible Higher accuracy 200-300

N =40 Not feasible Lower accuracy 80-150

MGPM M = 4, N =10 Higher accuracy Higher accuracy 50-90

of LG nodes (more than eighty LG nodes) such that the computation time commonly lasts for 3-5 min according to different objective functions. In contrast, fewer LG nodes are required in total for MGPM that practically increases the density of LG nodes around the quasi-contact point. Through extra numerical tests using GPM and MGPM with the same number of LG nodes, we have found that the solutions by MGPM take nearly half computation time compared with that by GPM.

Further discussion focuses on the advantages of the proposed method compared with the existing multiple-phase methods.20'21 The adaptive pseudospectral method proposed in Ref. 20 iteratively adjusts the segment interval and the degree of the polynomial required in each segment to obtain a solution to a specified accuracy. It is found to be successful on a range of problems. However, the mesh can only increase in size20 and the method terminates using more segments, which probably causes inefficiency for reentry trajectory optimization with complex constraints. In contrast, MGPM adjusts the phases according to the number of waypoints and no-fly zones such that only desired segments are utilized without redundant segments. Ref. 21 also proposed a multiple-phase method to solve non-sequential optimal control problems. It based on the Legendre pseudospectral method (LPM) that has a lower convergence rate than GPM.22 The method was demonstrated on an example of entry trajectory optimization with only traditional constraints. Thereby, MGPM is an extension of it aimed at waypoint and no-fly zone constraints for reentry trajectory optimization problems.

5. Conclusions

This paper presents a multi-phase technique based on GPM to generate an optimal reentry trajectory for HV satisfying the waypoint and no-fly zone constraints.

(1) Using MGPM, the waypoints are selected as the phase breaks of reentry trajectory, without considering the uncertainty of passage time for each waypoint. The way-point constraint is transcribed into the terminal condition in each phase, and then, can be easily solved as the typical constraint.

(2) The quasi-contact points are introduced to divide the full trajectory. The MGPM increases the number of LG nodes around each quasi-contact point, which eliminates the inaccuracy of trajectory near the no-fly zones.

(3) With fewer LG nodes in each phase, the MGPM decreases the order of interpolating polynomial, which leads to small approximate errors in the state and control.

Numerical results graphed here demonstrate the feasible application of MGPM to generate an optimal reentry trajectory satisfying waypoint and no-fly zone constraints. These results are also expected to contribute to the study of guidance

approaches for HV, especially when the expensive flying experiments are not available.

Acknowledgements

The authors are grateful to thank the anonymous reviewers for their critical and constructive review of the manuscript. This study was co-supported by Aviation Science Foundation of China (No. 2011ZC13001 and 2013ZA18001), National Natural Science Foundation of China (Nos: 60975073, 61273349, 61175109 and 61203223) and Innovation Foundation of BUAA for PhD Graduates.

References

1. USAF. Air force handbook. Washington D.C.: Department of the Air Force; 2005.

2. Vinh NX. Optimal trajectories in atmospheric flight. New York: Elsevier; 1981.

3. Huntington GT. Advancement and analysis of a Gauss pseudospectral transcription for optimal control problems dissertation. Cambridge, MA: Massachusetts Institute of Technology; 2007.

4. Ross IM, Fahroo F. A perspective on methods for trajectory optimization; 2002. Report No.: AIAA-2002-4727.

5. Betts JT, White EB. Practical methods for optimal control using nonlinear programming. Philadelphia: Society for Industrial and Applied Mechanics; 2001.

6. Reddien GW. Collocation at Gauss points as a discretization in optimal control. SIAM J Control Optim 1979;17(2):298-306.

7. Benson DA, Huntington GT, Thorvaldsen TP, Rao AV. Direct trajectory optimization and costate estimation via an orthogonal collocation method. J Guid Control Dyn 2006;29(6):1435-40.

8. Benson DA. A Gauss pseudospectral transcription for optimal control dissertation. Cambridge, MA: Massachusetts Institute of Technology; 2005.

9. Fahroo F, Ross IM. On discrete-time optimality conditions for pseudospectral methods; 2006. Report No.: AIAA-2006-6304.

10. Fahroo F, Ross IM. Advances in pseudospectral methods for optimal control; 2008. Report No.: AIAA-2008-7309.

11. Zhang D, Lu XF, Liu L, Wang YJ. An online ascent phase trajectory reconstruction algorithm using Gauss pseudospectral method. In: International conference on, modeling, identification and control; 2012. p. 1184-9.

12. Yong E, Tang G, Chen L. Rapid trajectory optimization for hypersonic reentry vehicle by Gauss pseudospectral method. J Astron 2008;29(6):1767-72 [Chinese.

13. Yao YW, Li HB. The generation of three-dimensional constrained entry trajectories for hypersonic vehicle based on the Gauss pseudospectral method. Aerosp Control 2012;30(2):33-8 [Chinese].

14. Sun Y, Zhang M. Optimal reentry range trajectory of hypersonic vehicle by Gauss pseudospectral method. In: International conference on intelligent control and information processing; 2011. p. 545-9.

15. Liu P, Zhao JS, Gu LX. Rapid approach for three-dimensional trajectory optimization of hypersonic reentry vehicles. Flight Dyn 2012;30(3):263-8 [Chinese].

16. Zhou WY, Yang D, Li SL. Solution of reentry trajectory with maximum cross range by using Gauss pseudospectral method. Syst Eng Electron 2010;32(5):1038-42 [Chinese].

17. Josselyn S, Ross IM. Rapid verification method for the trajectory optimization of reentry vehicles. J Guid Control Dyn 2002;26(3):505-8.

18. Phillips TH. A common aero vehicle (CAV) model, description, and employment guide. Technical Report, Schafer Corporation for AFRL and AF-SPC; 2003.

19. Jorris TR. Common aero vehicle autonomous reentry trajectory optimization satisfying waypoints and no-fly zone constraints. Wright-Patterson: Air Force Institute of Technology; 2007.

20. Darby CL, Hager WW, Rao AV. An hp-adaptive pseudospectral method for solving optimal control problems. Optim Control Appl Methods 2011;32(4):476-502.

21. Rao AV. Extension of a pseudospectral Legendre method to nonsequential multiple-phase optimal control problems; 2003. Report No.: AIAA-2003-5634.

22. Huntington GT, Benson D, Rao AV. A comparison of accuracy and computational efficiency of three pseudospectral methods; 2007. Report No.: AIAA-2007-6405.

Zhao Jiang received the B.S. degree in automation and M.S. degrees in control theory and control engineering from Northwestern Polytechnical University, Xi'an, China, in 2009 and 2012, respectively. He is currently working toward the Ph.D. degree in guidance, navigation, and control from Beihang University, Beijing, China. His research interests focus on reentry trajectory optimization and guidance of hypersonic vehicle, UAV path planning, and multiple missiles cooperative guidance.

Zhou Rui received the Ph.D. degree from Harbin Institute of Technology, Harbin, China, in 1997. He is currently a professor and Ph.D. supervisor in guidance, navigation, and control in School of Automation Science and Electrical Engineering, Beihang University, Beijing, China. His research interest focuses on advanced flight control and guidance, multiple UAVs cooperative control and intelligent decision, and multiple missiles cooperative guidance.